CN115188434A - 一种流域尺度水体非点源污染分类分区识别方法 - Google Patents

一种流域尺度水体非点源污染分类分区识别方法 Download PDF

Info

Publication number
CN115188434A
CN115188434A CN202210527218.7A CN202210527218A CN115188434A CN 115188434 A CN115188434 A CN 115188434A CN 202210527218 A CN202210527218 A CN 202210527218A CN 115188434 A CN115188434 A CN 115188434A
Authority
CN
China
Prior art keywords
point source
load
data
flow
source pollution
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202210527218.7A
Other languages
English (en)
Inventor
何圣嘉
俞珂
严琰
姜培坤
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Zhejiang A&F University ZAFU
Original Assignee
Zhejiang A&F University ZAFU
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Zhejiang A&F University ZAFU filed Critical Zhejiang A&F University ZAFU
Priority to CN202210527218.7A priority Critical patent/CN115188434A/zh
Publication of CN115188434A publication Critical patent/CN115188434A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • GPHYSICS
    • G16INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16CCOMPUTATIONAL CHEMISTRY; CHEMOINFORMATICS; COMPUTATIONAL MATERIALS SCIENCE
    • G16C20/00Chemoinformatics, i.e. ICT specially adapted for the handling of physicochemical or structural data of chemical particles, elements, compounds or mixtures
    • G16C20/70Machine learning, data mining or chemometrics
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F16/00Information retrieval; Database structures therefor; File system structures therefor
    • G06F16/20Information retrieval; Database structures therefor; File system structures therefor of structured data, e.g. relational data
    • G06F16/24Querying
    • G06F16/245Query processing
    • G06F16/2458Special types of queries, e.g. statistical queries, fuzzy queries or distributed queries
    • G06F16/2462Approximate or statistical queries
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F16/00Information retrieval; Database structures therefor; File system structures therefor
    • G06F16/20Information retrieval; Database structures therefor; File system structures therefor of structured data, e.g. relational data
    • G06F16/29Geographical information databases
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06QINFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q50/00Information and communication technology [ICT] specially adapted for implementation of business processes of specific business sectors, e.g. utilities or tourism
    • G06Q50/10Services
    • G06Q50/26Government or public services
    • GPHYSICS
    • G16INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16CCOMPUTATIONAL CHEMISTRY; CHEMOINFORMATICS; COMPUTATIONAL MATERIALS SCIENCE
    • G16C20/00Chemoinformatics, i.e. ICT specially adapted for the handling of physicochemical or structural data of chemical particles, elements, compounds or mixtures
    • G16C20/20Identification of molecular entities, parts thereof or of chemical compositions
    • GPHYSICS
    • G16INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16CCOMPUTATIONAL CHEMISTRY; CHEMOINFORMATICS; COMPUTATIONAL MATERIALS SCIENCE
    • G16C20/00Chemoinformatics, i.e. ICT specially adapted for the handling of physicochemical or structural data of chemical particles, elements, compounds or mixtures
    • G16C20/90Programming languages; Computing architectures; Database systems; Data warehousing
    • YGENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02ATECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A20/00Water conservation; Efficient water supply; Efficient water use
    • Y02A20/152Water filtration

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Databases & Information Systems (AREA)
  • Computing Systems (AREA)
  • Business, Economics & Management (AREA)
  • Software Systems (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Data Mining & Analysis (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Chemical & Material Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Crystallography & Structural Chemistry (AREA)
  • General Engineering & Computer Science (AREA)
  • Probability & Statistics with Applications (AREA)
  • General Health & Medical Sciences (AREA)
  • Tourism & Hospitality (AREA)
  • Health & Medical Sciences (AREA)
  • Educational Administration (AREA)
  • Computational Linguistics (AREA)
  • Medical Informatics (AREA)
  • Economics (AREA)
  • Fuzzy Systems (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Evolutionary Computation (AREA)
  • Development Economics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Mathematical Physics (AREA)
  • Human Resources & Organizations (AREA)
  • Marketing (AREA)
  • Primary Health Care (AREA)
  • Strategic Management (AREA)
  • General Business, Economics & Management (AREA)
  • Remote Sensing (AREA)
  • Artificial Intelligence (AREA)

Abstract

本发明公开了一种流域尺度水体非点源污染分类分区识别方法,步骤1:包括基础数据获取和非点源污染负荷定量;步骤2:基流分割定量和基流负荷分割定量;步骤3:获得所求地区负荷量数据、土地利用类型、某种土地利用类型的单位面积;步骤4:构建非点源污染源解析模型;步骤5:通过步骤1获得数据,完成非点源污染源解析模型的参数率定;步骤6:源解析模型结果验证及评价:将步骤1中的数据按时间顺序分为两个部分;步骤7:非点源污染源分类识别;步骤8:非点源污染源分区识别。本发明对输出系数模型进一步优化,实现流域尺度非点源污染更加科学、合理和准确的分类分区定量源解析,有利于更好地进行源头监测、责任认定,并对其进行有效削减。

Description

一种流域尺度水体非点源污染分类分区识别方法
技术领域
本发明涉及污水控制治理技术领域,尤其涉及一种流域尺度水体非点源污染分类分区识别方法。
背景技术
随着工业废水、城市污水等点源污染的控制和治理水平的不断提高,由农业生产、水土流失以及大气沉降等引起的非点源污染已经成为水环境的首要污染源。就不同类型非点源污染对水环境质量退化的贡献而言,又以农业非点源污染的贡献最大,因而由农业生产活动所引起的地表水环境非点源污染一直是国内外水体污染研究所关注的一个焦点问题。非点源污染具有成分复杂,类型多样,发生模式间歇性,发生时间不确定,排放和排放模式多样性、难以确定,时空变化迅速等特点。由于这些特点,使得非点源污染的监测和控制更加困难。
基流通常指河道中来源于地下水或其他滞后性水源的那部分径流。长期的农业集约化不仅导致地表水体中氮(N)、磷(P)浓度的增加,同时也造成了诸多地区地下水严重的富营养化倾向,使基流成为农业流域非点源N、P污染物输出的一个重要途径。然而,现有流域非点源污染定量溯源的方法通常将所有的污染负荷归结于地表径流,将严重误导流域非点源污染治理的科学认识和后续的工作策略;加之传统输出系数方法自身的不足(如缺少考虑污染物迁移损失和生物降解等),使得基于传统方法得到的流域非点源污染分类分区识别将存在较大的不确定性。
发明内容
1.要解决的技术问题
本发明的目的是为了解决流域尺度充分考虑基流非点源污染贡献的前提下,通过优化和改进传统输出系数模型,而建立的一种更加科学合理的水体非点源污染分类分区识别方法。
2.技术方案
为了实现上述目的,本发明采用了如下技术方案:
一种流域尺度水体非点源污染分类分区识别方法,包括以下步骤:
步骤1:包括基础数据获取和非点源污染负荷定量;其中基础数据获取包括给定流域出口处的水文(流量)、水质(污染物浓度)数据,整个流域的气象数据以及土地利用数据;非点源污染负荷定量以监测获得的离散水文(流量)、水质(污染物浓度)数据为基础,利用美国地质调查局开发的LOADEST模型计算得到流域出口处逐日的径流负荷量:
Figure BDA0003644841120000021
式中,
Figure BDA0003644841120000022
为利用修正的最大似然估计方法计算得到的瞬时河道总径流负荷量(kg·d-1);a0~aj为污染物通量回归方程的参数;H(a,b,s2,α,k)为无穷级数的似然逼近函数;a和b为解释变量的函数;α、κ为gamma分布的函数;m为自由度;s2为残差。
步骤2:基流及其负荷分割定量;
其中,所在流域的基流分割定量采用Eckhardt于2005年提出的双参数递归滤波(ERDF):
Figure BDA0003644841120000031
式中,qt、Qt分别为t时刻的总径流量和基流量;△t为时间步长(天);α为反映河流退水速度的退水常数(0<α<1);BFImax为最大基流指数。由于BFImax无法直接测定,Eckhardt(2005)给出了三个经验值:0.80、0.50和0.25,分别对应多孔介质含水层地区的常年性河流与季节性河流以及硬质岩介质含水层地区的常年性河流。
在实现基流分割定量的基础上,以与基流负荷消退过程密切相关的水文、气象因子为基础,构建基流负荷消退系数的回归统计模型,见公式3,用于估算逐日的基流负荷消退系数,在确定负荷消退系数之后,将其回代入递归数字滤波方程,见公式2,即可对LOADEST模型模拟得到的日负荷量进行分割,得到基流负荷量,这种基流负荷定量分割方法简称为递归滤波基流负荷分割算法(RFLSA),
τ(t)=γ1×a(t)+γ2×P(t)+γ3×E(t)+C (3)
式中,τ为基流负荷消退参数;a为退水参数;P和E分别为降雨量和蒸发量;t为时间(以天为单位);γ1~γ3为拟合参数;
步骤3:获得所求地区负荷量数据、土地利用类型、某种土地利用类型(水田、林地等)的单位面积,区域气象数据、大气沉降进入河流水体的污染物数据、点源污染产生的污染物;
步骤4:构建非点源污染源解析模型:
Figure BDA0003644841120000041
式中,L为非点源污染物入河量;BL为基流非点源污染负荷量;Xi为第i种土地利用类型某一时段(如1天、1周)单位面积某种污染物的排放量,Si为Xi的修正参数,可以表示为气象因子的一个函数;W(t)为某时段地表径流量与研究时段平均地表径流量的比值;M(t)为大气干沉降进入河流水体的污染物;D(t)为点源污染产生的污染物;n为土地利用类型数量;t为时间,一般为天或周,P为降雨数据、E为蒸发数据、T为平均温度数据、R为相对湿度数据,均为气象数据,η表示除函数中所包含因子以外其他所有因素的综合影响,θ1~θ4为数值优化求解得到的方程拟合参数;
步骤5:通过步骤1获得数据,完成非点源污染源解析模型的参数率定;
步骤6:校准与验证:即将步骤1中的数据按时间顺序分为两个部分,前一个部分用于步骤2中模型参数的确定和校准,后一部分则用于模型模拟结果的验证和评价,选用的评价指标为NSE、均方根误差-实测值标准差比率(RSR)和决定系数(R2),具体的计算方法如下:
Figure BDA0003644841120000043
Figure BDA0003644841120000051
Figure BDA0003644841120000052
式中,L(i,m)和L(i,s)分别表示第i个污染土地负荷量的实测值和模拟值;Lavg表示污染物负荷量的实测值的平均值;n为实测值个数;
所述实测值可以通过上述L=KCQ公式基于水流瞬时流量Q,以及该流量条件下对应的污染物浓度C连续数据计算获得,模拟值为分类分区识别的公式出来的;
L=KCQ公式中,L为负荷;K为单位换算系数;C浓度,Q为河水流量,该公式为水体污染负荷的估算公式,污染负荷量是指在一定时段内,由点污染源和面污染源,进入水体的污染物总量;
步骤7:非点源污染源分类识别,利用ARCGIS软件和相应的土地利用图完成该地区各类土地利用面积的统计汇总;结合已经完成参数优化与校验的输出周尺度的输出系数模型,计算得到所在流域不同地类给定时段内地表径流的非点源污染排放量。
步骤8:非点源污染源分区识别,利用ARCGIS软件的统计分析功能对流域内各种数据进行统计分析,通过全国乡镇分区将研究区划分成所需要等级,利用上述得到的周尺度的输出系数模型,将各行政区土地面积及所需时间段内气象数据代入,分别得到研究区内各行政区块给定时段内地表径流的非点源污染排放量。
优选地,所述步骤1中逐日水文、气象数据由政府相关部门提供;水质数据通过定期监测得到,具体指标的分析测定方法参考国家相关标准,根据2017最新版《土地利用现状分类》(GBT 21010-2017),利用ARCGIS软件和遥感影像,完成给定流域各类土地利用的分类统计。
优选地,所述步骤2中由于BFImax目前无法通过直接测定来获取,BFImax的值主要是在结合流域水文地质特征的基础上,利用数学算法计算模拟得到:多孔介质含水层地区的常年性河流与季节性河流以及硬质岩介质含水层地区的常年性河流的BFImax建议值分别为0.80、0.50和0.25。
优选地,其中ERDF参数的确定方法包括以下步骤:
S1:根据经验公式N=0.83A0.2(N为洪水峰值后地表径流完全停止所需天数;A为流域面积,单位为km2)确定纯基流退水的起始点;
S2:将满足条件y1>y2>…>yk>yk+1>yk+2的流量数据yk和yk+1从日流量序列中筛选出来,得到纯基流退水过程;所选退水过程的最小长度为5天,在降雨发生频率较高的地区,退水分析时需尽量排除降雨的干扰;
S3:利用过原点的线性方程拟合散点图(yk vs yk+1)的上边界点,由此得到的方程斜率即为ERDF的退水常数;
S4:根据此方法计算全过程的平均退水常数和月尺度退水常数;
S5:为了减少不同基流退水过程基流消退规律变异性以及参数BFImax经验性取值造成的基流分割结果的不确定性,以月退水常数为基础,得到日退水常数的一阶傅里叶拟合函数;然后以此函数为基础,构建ERDF的日退水常数的计算方程,见公式8,利用退水分析筛选得到的日基流量和遗传算法实现方程参数的优化求解,计算得到逐日的退水常数以及BFImax最优值,
Figure BDA0003644841120000071
式中,a(t)为日基流退水常数优化方程;Fa(t)为月退水常数傅里叶拟合函数;R(t)为修正函数;E、P分别为蒸发和降雨;t为时间,通常以1天为步长;dtime为校正之后的分数时间(dtime=decimal time-center of decimal time);β0~β3为拟合参数,其优化目标函数如下:
Figure BDA0003644841120000072
式中,Qi为第i天基流量;上标obs、sim和mean分别表示实测值、模拟值和平均值,其中基流模拟值是由ERDF分割得到,基流实测值则是根据退水分析,从日实测径流量序列中筛选得到。
优选地,所述步骤4中为了增加模型精度,使其更适用于研究区各土地利用类型,Xi的取值在通过查阅文献确定范围区间的条件下,通过遗传算法进行拟合。
优选地,所述步骤5中利用基于MATLAB软件的遗传算法进行识别方程的参数优化求解。
优选地,所述步骤7中根据2017最新版《土地利用现状分类》(GBT 21010-2017),利用ARCGIS软件和遥感影像,完成给定流域各类土地利用的分类统计。
3.有益效果
相比于现有技术,本发明的优点在于:
本发明中,通过对输出系数模型作出进一步的改进,有利于从通过模型模拟,实现研究区非点源污染更加科学、合理和准确的分类分区定量源解析,有利于更好地进行源头监测、责任认定,并对其进行有效削减。
附图说明
图1为本发明中实施例3中步骤Ⅰ的示意图;
图2为本发明中实施例3中步骤Ⅱ的示意图;
图3为本发明中实施例3中步骤Ⅲ的示意图;
图4为本发明中实施例3中步骤Ⅲ的示意图;
图5为本发明中实施例3中步骤Ⅳ的示意图;
图6为本发明中实施例4中基流分割结果验证图;
图7为本发明中实施例4中基流TN负荷量结果验证图;
图8为本发明中实施例4中上梧溪流域2020年11月—2021年10月逐日河道总径流及基流量TN负荷量图;
图9为本发明中实施例4中上梧溪地表径流TN各土地利用类型地表径流非点源污染输出系数图;
图10为本发明中实施例4中2020年11月—2021年10月上梧溪流域非点源氮污染来源图。
具体实施方式
下面将结合本发明实施例中的附图,对本发明实施例中的技术方案进行清楚、完整地描述,显然,所描述的实施例仅仅是本发明一部分实施例,而不是全部的实施例。
实施例1:
一种流域尺度水体非点源污染分类分区识别方法,包括以下步骤:
步骤1:包括基础数据获取和非点源污染负荷定量;其中基础数据获取包括给定流域出口处的水文(流量)、水质(污染物浓度)数据,整个流域的气象数据以及土地利用数据,其中,逐日水文、气象数据由政府相关部门提供;水质数据通过定期(如以月为步长)监测得到,具体指标的分析测定方法参考国家相关标准。根据2017最新版《土地利用现状分类》(GBT 21010-2017),利用ARCGIS软件和遥感影像,完成给定流域各类土地利用的分类统计;
非点源污染负荷定量以监测获得的离散水文(流量)、水质(污染物浓度)数据为基础,利用美国地质调查局开发的LOADEST模型计算得到流域出口处逐日的径流负荷量:
Figure BDA0003644841120000091
式中,
Figure BDA0003644841120000092
为利用修正的最大似然估计方法计算得到的瞬时河道总径流负荷量(kg·d-1);a0~aj为污染物通量回归方程的参数;H(a,b,s2,α,k)为无穷级数的似然逼近函数;a和b为解释变量的函数;α、κ为gamma分布的函数;m为自由度;s2为残差。
步骤2:基流及其负荷分割定量;
其中,所在流域的基流分割定量利用Eckhardt于2005年提出的双参数递归滤波(ERDF):
Figure BDA0003644841120000101
式中,qt、Qt分别为t时刻的总径流量和基流量;△t为时间步长(天);α为反映河流退水速度的退水常数(0<α<1);BFImax为最大基流指数。由于BFImax无法直接测定,Eckhardt(2005)给出了三个经验值:0.80、0.50和0.25,分别对应多孔介质含水层地区的常年性河流与季节性河流以及硬质岩介质含水层地区的常年性河流。
其中基流负荷分割定量以与基流负荷消退过程密切相关的水文、气象因子为基础,构建基流负荷消退系数的回归统计模型,见公式3,用于估算逐日的基流负荷消退系数,在确定负荷消退系数之后,将其回代入递归数字滤波方程,见公式2,即可对LOADEST模型模拟得到的日负荷量进行分割,得到基流负荷量,这种基流负荷定量分割方法简称为递归滤波基流负荷分割算法(RFLSA),
τ(t)=γ1×a(t)+γ2×P(t)+γ3×E(t)+C (3)
式中,τ为基流负荷消退参数;a为退水参数;P和E分别为降雨量和蒸发量;t为时间(以天为单位);γ1~γ3为拟合参数;
步骤3:获得所求地区负荷量数据、土地利用类型、某种土地利用类型(水田、林地等)的单位面积,区域气象数据、大气沉降进入河流水体的污染物数据、点源污染产生的污染物;
步骤4:构建非点源污染源解析模型:
Figure BDA0003644841120000111
式中,L为非点源污染物入河量;BL为基流非点源污染负荷量;Xi为第i种土地利用类型某一时段(如1天、1周)单位面积某种污染物的排放量,Si为Xi的修正参数,可以表示为气象因子的一个函数;W(t)为某时段地表径流量与研究时段平均地表径流量的比值;M(t)为大气干沉降进入河流水体的污染物;D(t)为点源污染产生的污染物;n为土地利用类型数量;t为时间,一般为天或周,P为降雨数据、E为蒸发数据、T为平均温度数据、R为相对湿度数据,均为气象数据,η表示除函数中所包含因子以外其他所有因素的综合影响,θ1~θ4为数值优化求解得到的方程拟合参数;
为了增加模型精度,使其更适用于研究区各土地利用类型,Xi的取值在通过查阅文献确定范围区间的条件下,通过遗传算法进行拟合;
步骤5:通过步骤1获得数据,完成非点源污染源解析模型的参数率定;
步骤6:源解析模型结果验证及评价:即将步骤1中的数据按时间顺序分为两个部分,前一个部分用于步骤2中模型参数的确定和校准,后一部分则用于模型模拟结果的验证和评价,选用的评价指标为NSE、均方根误差-实测值标准差比率(RSR)和决定系数(R2),具体的计算方法如下:
Figure BDA0003644841120000121
Figure BDA0003644841120000122
Figure BDA0003644841120000123
式中,L(i,m)和L(i,s)分别表示第i个污染土地负荷量的实测值和模拟值;Lavg表示污染物负荷量的实测值的平均值;n为实测值个数;
所述实测值可以通过上述L=KCQ公式基于水流瞬时流量Q,以及该流量条件下对应的污染物浓度C连续数据计算获得,模拟值为分类分区识别的公式出来的;
L=KCQ公式中,L为负荷;K为单位换算系数;C浓度,Q为河水流量,该公式为水体污染负荷的估算公式,污染负荷量是指在一定时段内,由点污染源和面污染源,进入水体的污染物总量;
步骤7:非点源污染源分类识别,利用ARCGIS软件的统计分析功能对流域内各种数据进行统计分析,得到研究区流域各类土地利用面积;
根据2017最新版《土地利用现状分类》(GBT 21010-2017),利用ARCGIS软件和遥感影像,完成给定流域各类土地利用的分类统计。将研究区域土地利用类型分为水田、草地、林地、水域、旱地、人居地。利用上述得到的改进后的输出系数模型,将各地类土地面积及所需时间段内气象数据代入,分别得到研究区各地类一定时间内地表径流非点源污染排放量;
步骤8:非点源污染源分区识别,利用ARCGIS软件的统计分析功能对流域内各种数据进行统计分析,通过全国乡镇分区将研究区划分成所需要等级,利用上述得到的周尺度的输出系数模型,将各行政区土地面积及所需时间段内气象数据代入,分别得到研究区内各行政区块一定时间内地表径流排放量。
本实施例中,通过分类分区识别可针对性的,更加合理、科学、准确的进行非点源污染的治理。
实施例2
其具有上述实施例的实施内容,其中,对于上述实施例的具体实施方式可参阅上述描述,此处的实施例不作重复详述;而在本申请实施例中,其与上述实施例的区别在于:
本实施例中,其中ERDF参数的确定方法包括以下步骤:
S1:根据经验公式N=0.83A0.2(N为洪水峰值后地表径流完全停止所需天数;A为流域面积,单位为km2)确定纯基流退水的起始点;
S2:将满足条件y1>y2>…>yk>yk+1>yk+2的流量数据yk和yk+1从日流量序列中筛选出来,得到纯基流退水过程;所选退水过程的最小长度为5天,在降雨发生频率较高的地区,退水分析时需尽量排除降雨的干扰;
S3:利用过原点的线性方程拟合散点图(yk vs yk+1)的上边界点,由此得到的方程斜率即为ERDF的退水常数;
S4:根据此方法计算全过程的平均退水常数和月尺度退水常数;
S5:为了减少不同基流退水过程基流消退规律变异性以及参数BFImax经验性取值造成的基流分割结果的不确定性,以月退水常数为基础,得到日退水常数的一阶傅里叶拟合函数;然后以此函数为基础,构建ERDF的日退水常数的计算方程,见公式8,利用退水分析筛选得到的日基流量和遗传算法实现方程参数的优化求解,计算得到逐日的退水常数以及BFImax最优值,
Figure BDA0003644841120000141
式中,a(t)为日基流退水常数优化方程;Fa(t)为月退水常数傅里叶拟合函数;R(t)为修正函数;E、P分别为蒸发和降雨;t为时间,通常以1天为步长;dtime为校正之后的分数时间(dtime=decimal time-center of decimal time);β0~β3为拟合参数,其优化目标函数如下:
Figure BDA0003644841120000142
式中,Qi为第i天基流量;上标obs、sim和mean分别表示实测值、模拟值和平均值,其中基流模拟值是由ERDF分割得到,基流实测值则是根据退水分析,从日实测径流量序列中筛选得到。
实施例3
其具有上述实施例的实施内容,其中,对于上述实施例的具体实施方式可参阅上述描述,此处的实施例不作重复详述;而在本申请实施例中,其与上述实施例的区别在于:
参照图1-5,MATLAB进行识别方程参数数值优化求解的过程如下步骤:
Ⅰ:打开MATLAB程序,在主工具栏找到应用程序一栏,打开Optimization选项;
Ⅱ:在Optimization程序中,首先根据不同的问题类型选择不同的模型,同时输入约束等;
Ⅲ:选择优化程序运行的条件。在最中间的一栏中添加优化程序运行的条件,如优化截止的标准等等;
Ⅳ:开始运算,点击“Start”即可实现优化程序的运行,在图示的框中即可出现运行结果。
本实施例中,纳什系数越接近于1,代表该模型的模拟结果越可靠,将纳什系数最大时对应的运行结果框中各未知数值代入至识别方程中,便得到了相应的改进后的输出系数模型。
实施例4
其具有上述实施例的实施内容,其中,对于上述实施例的具体实施方式可参阅上述描述,此处的实施例不作重复详述;而在本申请实施例中,其与上述实施例的区别在于:
本实施例中,在充分考虑上梧溪流域水系分布、河道地形、污染物排放状况、植被与水土流失等因素的基础上,布设水质监测站点,于2020年11月至2021年10月,对上梧溪流域出口处的水文水质状况进行定期监测。
本实施例中,水质监测频率为每周一次;水文流量监测采用流速仪法,以日为步长,具体方法参考《河流流量测验规范》(GB50179—2015)。水样从河道中间离水面30cm处采集,装于2.5L聚乙烯瓶之中。水样采集完成以后立即将水样运回实验室,在采水后24小时内,利用碱性过硫酸钾消解紫外分光光度法测总氮,以监测获得的离散水文水质数据为基础,利用Load Estimator(LOADEST)模型计算本研究区日尺度上的非点源氮污染负荷量。
本实施例中,通过淳安县气象局得到2020.11-2021.11淳安县气象站逐日气象数据。利用ARCGIS软件的统计分析功能对流域内各种数据进行统计分析,得到上梧溪流域各类土地利用面积。
本实施例中,根据淳安县政府研究报告,因此不考虑农村居民生活的点源污染对本流域TN负荷量的影响。通过查阅相关文献资料,并结合研究区实际情况,千岛湖流域TN的大气干沉降量为2.24kg·(hm·a)-1
本文调用LOADEST模型的模型自动选择项,得到候选模型的模拟结果如表1、2所示。
表1上梧溪流域逐日TN负荷定量回归方程
Figure BDA0003644841120000161
本实施例中,注:公式中L为污染物日负荷量(t·d-1);Q为径流量(mm);dtime为中心化之后的分数日期(dtime=分数日期—中心分数日期);a0、a1、a2、a3、a4、a5和a6为方程系数;SCR为残差序列相关系数,用于源解析模型结果验证及评价残差是否存在序列相关性,其值越小,说明残差之间相互独立,方程中涉及的变量不存在序列相关性;PPCC为概率曲线相关系数,源解析模型结果验证及评价模型的残差是否服从标准正态分布。
表2上梧溪流域逐日TN负荷估算的LOADEST模型参数率定及评价
Figure BDA0003644841120000171
本实施例中,通过用基于拟合傅里叶函数和分数日期反推得到的日退水常数替换单一的年退水常数,使其基流分割结果有较高的准确性和较好的模拟精度(NSE=0.81,RSR=0.44,R2=0.94,图6)。
本实施例中,基于RFLSA的基流TN负荷定量结果不仅能够较好地反映上梧溪流域日尺度上基流TN实测负荷量的变化趋势(R2=0.91),同时还具有较高的模拟精度(NSE=0.89,RSR=0.33,图7)。上梧溪流域2020年11月—2021年10月逐日河道总径流量、基基流量及其相应的TN负荷量变化如图8所示。
本实施例中,将上梧溪流域2020.11-2021.10地表径流TN排放量分为单周与双周,将单周数据用来模型参数的率定,将双周数据用来进行模型的验证与评价。
本实施例中,运用MATLAB软件的遗传算法对上梧溪流域2020年11月—2021年10月每周地表径流TN排放量和降雨(P)、蒸发(E)、平均温度(T)、相对湿度(R)等气象因子做相关性分析,得到关于Si(t)的拟合方程:
Si(t)=0.000097×P(t)-0.000223×E(t)-0.000287×T(t)+0.000145×R(t)+0.001818
本实施例中,通过查阅江浙地区为研究区的参考文献并结合研究区的实际情况,大致确定各土地利用类型输出系数的取值范围,保证输出系数的取值具有相对的精确性。利用改进的输出系数模型,通过MATLAB遗传算法进行参数拟合,得到以周为步长的各土地利用输出系数值,如表3所示:
表3上梧溪地表径流TN各土地利用类型周输出系数
Figure BDA0003644841120000181
本实施例中,通过将双周数据代入进行验证,运用通过MATLAB软件的遗传算法进行参数优化得到的识别方程计算出的地表径流TN负荷量能够较好地反映上梧溪流域周(时间)尺度上地表径流TN实测负荷量的变化趋势(R2=0.80),同时还具有较高的模拟精度(NSE=0.80,RSR=0.44)。
本实施例中,利用改进的输出系数模型计算得到研究区各土地利用类型的非点源污染TN排放量。2020年11月—2021年10月,上梧溪流域地表径流的总氮负荷量为56.06吨,如图10所示,其中,水田、草地、林地、旱地、水域、人居地、大气沉降的非点源氮污染负荷量为5761.23kg、5507.77kg、28668.62kg、750.32kg、5.83kg、26.79kg、15335.04kg,贡献率分别为10.28%、9.83%、51.14%、1.34%、0.01%、0.05%、27.36%。本实施例中,其中林地(含经济林)的地表径流TN负荷贡献率最高,是由于在该流域中林地占总面积的比例高达79.57%
以上所述,仅为本发明较佳的具体实施方式,但本发明的保护范围并不局限于此,任何熟悉本技术领域的技术人员在本发明揭露的技术范围内,根据本发明的技术方案及其发明构思加以等同替换或改变,都应涵盖在本发明的保护范围之内。

Claims (7)

1.一种流域尺度水体非点源污染分类分区识别方法,其特征在于,包括以下步骤:
步骤1:包括基础数据获取和非点源污染负荷定量;其中基础数据获取包括给定流域出口处的水文(流量)、水质(污染物浓度)数据,整个流域的气象数据以及土地利用数据;非点源污染负荷定量以监测获得的离散水文(流量)、水质(污染物浓度)数据为基础,利用美国地质调查局开发的LOADEST模型计算得到流域出口处逐日的径流负荷量:
Figure FDA0003644841110000011
式中,
Figure FDA0003644841110000012
为利用修正的最大似然估计方法计算得到的瞬时河道总径流负荷量(kg·d-1);a0~aj为污染物通量回归方程的参数;H(a,b,s2,α,k)为无穷级数的似然逼近函数;a和b为解释变量的函数;α、κ为gamma分布的函数;m为自由度;s2为残差;
步骤2:基流及其负荷分割定量;
其中,所在流域的基流分割定量采用Eckhardt于2005年提出的双参数递归滤波(ERDF):
Figure FDA0003644841110000013
式中,qt、Qt分别为t时刻的总径流量和基流量;△t为时间步长(天);α为反映河流退水速度的退水常数(0<α<1);BFImax为最大基流指数,由于BFImax无法直接测定,Eckhardt(2005)给出了三个经验值:0.80、0.50和0.25,分别对应多孔介质含水层地区的常年性河流与季节性河流以及硬质岩介质含水层地区的常年性河流;
在实现基流分割定量的基础上,以与基流负荷消退过程密切相关的水文、气象因子为基础,构建基流负荷消退系数的回归统计模型,见公式3,用于估算逐日的基流负荷消退系数,在确定负荷消退系数之后,将其回代入递归数字滤波方程,见公式2,即可对LOADEST模型模拟得到的日负荷量进行分割,得到基流负荷量,这种基流负荷定量分割方法简称为递归滤波基流负荷分割算法(RFLSA),
τ(t)=γ1×a(t)+γ2×P(t)+γ3×E(t)+C (3)
式中,τ为基流负荷消退参数;a为退水参数;P和E分别为降雨量和蒸发量;t为时间(以天为单位);γ1~γ3为拟合参数;
步骤3:获得所求地区负荷量数据、土地利用类型、某种土地利用类型(水田、林地等)的单位面积,区域气象数据、大气沉降进入河流水体的污染物数据、点源污染产生的污染物;
步骤4:构建非点源污染源解析模型:
Figure FDA0003644841110000021
式中,L为非点源污染物入河量;BL为基流非点源污染负荷量;Xi为第i种土地利用类型某一时段(如1天、1周)单位面积某种污染物的排放量,Si为Xi的修正参数,可以表示为气象因子的一个函数;W(t)为某时段地表径流量与研究时段平均地表径流量的比值;M(t)为大气干沉降进入河流水体的污染物;D(t)为点源污染产生的污染物;n为土地利用类型数量;t为时间,一般为天或周,P为降雨数据、E为蒸发数据、T为平均温度数据、R为相对湿度数据,均为气象数据,η表示除函数中所包含因子以外其他所有因素的综合影响,θ1~θ4为数值优化求解得到的方程拟合参数;
步骤5:通过步骤1获得数据,完成非点源污染源解析模型的参数率定;
步骤6:校准与验证:即将步骤1中的数据按时间顺序分为两个部分,前一个部分用于步骤2中模型参数的确定和校准,后一部分则用于模型模拟结果的验证和评价,选用的评价指标为NSE、均方根误差-实测值标准差比率(RSR)和决定系数(R2),具体的计算方法如下:
Figure FDA0003644841110000031
Figure FDA0003644841110000032
Figure FDA0003644841110000033
式中,L(i,m)和L(i,s)分别表示第i个污染土地负荷量的实测值和模拟值;Lavg表示污染物负荷量的实测值的平均值;n为实测值个数;
所述实测值可以通过上述L=KCQ公式基于水流瞬时流量Q,以及该流量条件下对应的污染物浓度C连续数据计算获得,模拟值为分类分区识别的公式出来的;
L=KCQ公式中,L为负荷;K为单位换算系数;C浓度,Q为河水流量,该公式为水体污染负荷的估算公式,污染负荷量是指在一定时段内,由点污染源和面污染源,进入水体的污染物总量;
步骤7:非点源污染源分类识别,利用ARCGIS软件和相应的土地利用图完成该地区各类土地利用面积的统计汇总;结合已经完成参数优化与校验的输出周尺度的输出系数模型,计算得到所在流域不同地类给定时段内地表径流的非点源污染排放量;
步骤8:非点源污染源分区识别,利用ARCGIS软件的统计分析功能对流域内各种数据进行统计分析,通过全国乡镇分区将研究区划分成所需要等级,利用上述得到的周尺度的输出系数模型,将各行政区土地面积及所需时间段内气象数据代入,分别得到研究区内各行政区块给定时段内地表径流的非点源污染排放量。
2.根据权利要求1所述的一种流域尺度水体非点源污染分类分区识别方法,其特征在于,所述步骤1中逐日水文、气象数据由政府相关部门提供;水质数据通过定期(如以月为步长)监测得到,具体指标的分析测定方法参考国家相关标准,根据2017最新版《土地利用现状分类》(GBT 21010-2017),利用ARCGIS软件和遥感影像,完成给定流域各类土地利用的分类统计。
3.根据权利要求1所述的一种流域尺度水体非点源污染分类分区识别方法,其特征在于,所述步骤2中由于BFImax目前无法通过直接测定来获取,BFImax的值主要是在结合流域水文地质特征的基础上,利用数学算法计算模拟得到:多孔介质含水层地区的常年性河流与季节性河流以及硬质岩介质含水层地区的常年性河流的BFImax建议值分别为0.80、0.50和0.25。
4.根据权利要求3所述的一种流域尺度水体非点源污染分类分区识别方法,其特征在于,其中ERDF参数的确定方法包括以下步骤:
S1:根据经验公式N=0.83A0.2(N为洪水峰值后地表径流完全停止所需天数;A为流域面积,单位为km2)确定纯基流退水的起始点;
S2:将满足条件y1>y2>…>yk>yk+1>yk+2的流量数据yk和yk+1从日流量序列中筛选出来,得到纯基流退水过程;所选退水过程的最小长度为5天,在降雨发生频率较高的地区,退水分析时需尽量排除降雨的干扰;
S3:利用过原点的线性方程拟合散点图(yk vs yk+1)的上边界点,由此得到的方程斜率即为ERDF的退水常数;
S4:根据此方法计算全过程的平均退水常数和月尺度退水常数;
S5:为了减少不同基流退水过程基流消退规律变异性以及参数BFImax经验性取值造成的基流分割结果的不确定性,以月退水常数为基础,得到日退水常数的一阶傅里叶拟合函数;然后以此函数为基础,构建ERDF的日退水常数的计算方程,见公式8,利用退水分析筛选得到的日基流量和遗传算法实现方程参数的优化求解,计算得到逐日的退水常数以及BFImax最优值,
Figure FDA0003644841110000061
式中,a(t)为日基流退水常数优化方程;Fa(t)为月退水常数傅里叶拟合函数;R(t)为修正函数;E、P分别为蒸发和降雨;t为时间,通常以1天为步长;dtime为校正之后的分数时间(dtime=decimal time-center of decimal time);β0~β3为拟合参数,其优化目标函数如下:
Figure FDA0003644841110000062
式中,Qi为第i天基流量;上标obs、sim和mean分别表示实测值、模拟值和平均值,其中基流模拟值是由ERDF分割得到,基流实测值则是根据退水分析,从日实测径流量序列中筛选得到。
5.根据权利要求1所述的一种流域尺度水体非点源污染分类分区识别方法,其特征在于,所述步骤4中为了增加模型精度,使其更适用于研究区各土地利用类型,Xi的取值在通过查阅文献确定范围区间的条件下,通过遗传算法进行拟合。
6.根据权利要求1所述的一种流域尺度水体非点源污染分类分区识别方法,其特征在于,所述步骤5中利用基于MATLAB软件的遗传算法进行识别方程的参数优化求解。
7.根据权利要求1所述的一种流域尺度水体非点源污染分类分区识别方法,其特征在于,所述步骤7中根据2017最新版《土地利用现状分类》(GBT 21010-2017),利用ARCGIS软件和遥感影像,完成给定流域各类土地利用的分类统计。
CN202210527218.7A 2022-05-16 2022-05-16 一种流域尺度水体非点源污染分类分区识别方法 Pending CN115188434A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210527218.7A CN115188434A (zh) 2022-05-16 2022-05-16 一种流域尺度水体非点源污染分类分区识别方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210527218.7A CN115188434A (zh) 2022-05-16 2022-05-16 一种流域尺度水体非点源污染分类分区识别方法

Publications (1)

Publication Number Publication Date
CN115188434A true CN115188434A (zh) 2022-10-14

Family

ID=83513120

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210527218.7A Pending CN115188434A (zh) 2022-05-16 2022-05-16 一种流域尺度水体非点源污染分类分区识别方法

Country Status (1)

Country Link
CN (1) CN115188434A (zh)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115712813A (zh) * 2022-11-23 2023-02-24 北京中科三清环境技术有限公司 确定非点源的污染负荷输出系数的方法、装置及电子设备
CN116564431A (zh) * 2023-06-02 2023-08-08 江苏捷利达环保科技有限公司 一种基于大数据处理的污染源在线分析***及方法
CN117674293A (zh) * 2023-12-07 2024-03-08 华能西藏雅鲁藏布江水电开发投资有限公司 一种梯级水电站的长期发电优化调度方法及装置
CN118152890A (zh) * 2024-04-02 2024-06-07 中国电建集团西北勘测设计研究院有限公司 流域相似性的分类方法及装置、存储介质、设备

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115712813A (zh) * 2022-11-23 2023-02-24 北京中科三清环境技术有限公司 确定非点源的污染负荷输出系数的方法、装置及电子设备
CN116564431A (zh) * 2023-06-02 2023-08-08 江苏捷利达环保科技有限公司 一种基于大数据处理的污染源在线分析***及方法
CN116564431B (zh) * 2023-06-02 2024-01-09 江苏捷利达环保科技有限公司 一种基于大数据处理的污染源在线分析***及方法
CN117674293A (zh) * 2023-12-07 2024-03-08 华能西藏雅鲁藏布江水电开发投资有限公司 一种梯级水电站的长期发电优化调度方法及装置
CN118152890A (zh) * 2024-04-02 2024-06-07 中国电建集团西北勘测设计研究院有限公司 流域相似性的分类方法及装置、存储介质、设备

Similar Documents

Publication Publication Date Title
CN115188434A (zh) 一种流域尺度水体非点源污染分类分区识别方法
CN112966926B (zh) 一种基于集成学习的洪水敏感性风险评估方法
CN114168906B (zh) 一种基于云计算的测绘地理信息数据采集***
CN112001610A (zh) 农业面源污染的处理方法及装置
CN109933637B (zh) 一种洪水风险动态展示及分析***
CN114723179A (zh) 一种基于水量水质联合管控的水环境预警溯源***和方法
CN114254802B (zh) 气候变化驱动下植被覆盖时空变化的预测方法
CN111090831A (zh) 一种湖泊面积变化关键驱动因子识别方法
CN113011993A (zh) 基于标准数据的农业污染源入水体负荷量测算方法
CN117974404B (zh) 一种水陆协同的陆域污染源分析方法及***
CN111626103A (zh) 一种基于夜晚灯光遥感影像的城市化梯度度量方法
CN116644889A (zh) 一种水质污染溯源方法及终端
CN112287299A (zh) 河流健康变化定量归因方法、装置及***
CN115409374A (zh) 城镇化差异对流域水生态服务功能的影响预警方法
CN114595631A (zh) 一种基于efdc模型和机器学习算法的水质预测方法
Wu et al. Impact of surface and underground water uses on streamflow in the upper-middle of the Weihe River basin using a modified WetSpa model
CN116050703B (zh) 村镇建设与资源环境承载力协调评价及协调模式分析方法
Pomari et al. A new tool to assess ecosystem health in large subtropical reservoirs: Development and validation of a Planktonic Index of Biotic Integrity
CN115936545A (zh) 一种基于水热互补方程的两参数月尺度水文模型建立方法
CN112001562A (zh) 一种灌溉用水数据校验方法及装置
CN111310124A (zh) 一种城市径流数据处理方法
CN115453664A (zh) 一种适用于无资料地区的降雨径流预报方法
CN115186960A (zh) 一种城镇区域污水有效收集处理能力的核算方法及装置
CN113361114A (zh) 一种基于径流路径的多尺度面源污染物入河系数测算方法
Zhang et al. Load quantification and effect evaluation of urban non-point source pollution in the Licun River based on SWAT model

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination