CN102663267B - 一种半湿润区流域面源污染负荷的确定方法 - Google Patents

Abstract

本发明公开了一种半湿润区流域面源污染负荷的确定方法，属于污染控制领域。其步骤为：（A）收集整理水文气象和水环境历史观测资料以及野外实验资料；（B）划分子流域和水文响应单元；（C）改进水文模块；（D）采用MUSLE方程进行土壤侵蚀量的计算；（E）进行面源污染负荷量的计算；（F）将计算得到的水文过程和面源污染负荷量与实测数据进行对比，确定水文模块、MUSLE方程以及水污染负荷计算模块的参数值，完成模型的建立。本方法能够有效地模拟半湿润地区淮河流域典型区域的面源污染负荷估算，给出合理科学的计算结果，提高了模型应用的普及性和有效性，能更好地服务于水环境保护与水污染控制等水环境研究项目。

Description

一种半湿润区流域面源污染负荷的确定方法

技术领域

本发明涉及一种半湿润区流域面源污染负荷的确定方法，具体地说，是一种通过改进原有面源污染负荷模型中的水文模块以精确反映半湿润区水文过程对面源污染计算的影响，从而提高流域面源污染负荷计算精度的数值模拟方法。

背景技术

目前，许多发达国家的研究已经证实，面源污染是导致水环境恶化的主要原因之一。据美国、日本等国报道，即使点源污染得到全面控制以后，江河的水质达标率仅为65％，湖泊的水质达标率为42％，海域水质达标率为78％。研究指出美国的面源污染量占污染总量的2/3，农业面源贡献率占75％左右。按照1987年经济合作与发展组织（OECD）报告，在其成员国中，硝酸盐与杀虫剂为最大的非点源污染，欧洲国家也得到相似的结论，荷兰农业面源提供的总氮、总磷分别占水环境污染的60％和40％-50％。

面源污染主要集中在水土流失严重的地区和主要的农业区。与通过集中排污口排放的点源污染相比，农业面源污染具有污染发生时间的随机性、发生方式的间歇性、机理过程复杂性、排放途径及排放量的不确定性、污染负荷时空变异性及监测、模拟与控制困难性等特点。目前，有关农业面源污染定量评价数据非常缺乏，尤其是淮河流域，对农业面源污染的影响性质、程度、范围缺乏深入、全面、系统研究，水文特征影响下的农业面源污染物到环境水体的整个输移途径及其过程不清、污染贡献不明，农业面源污染控制措施的制定和实施缺乏有效必要的科学依据。

近些年来，人们对农业面源污染识别和治理能力逐渐增强，在面源污染发生机制、影响因素评价、输运特征及入河贡献、治理控制措施方面，国内外均开展了大量研究并取得了重要进展。建立流域时间和空间序列上的水环境模拟模型是流域面源污染负荷量化评价最为有效和直接的方法。

在20世纪70年代后期至90年代，非点源污染机理模型研究兴起，一大批机理复杂、适用性强的模型出现，并形成与计算机和3S技术综合发展的趋势。随着对非点源污染物理化学过程研究的深入和对非点源过程的广泛监测，机理模型逐渐成为非点源模型开发的主要方向，国外的非点源污染模型己由简单的统计分析向机理模型、由平均负荷输出或单场暴雨分析向连续的时间响应分析、由集总式模型向分布式模型发展。

本世纪以来，模型与GIS、遥感等工具的结合日益紧密，非点源污染模型功能和运行效率更高，模型在非点源污染研究中应用更加广泛。在模型的应用方面，以SWAT模型为例，迄今为止，SWAT模型的有效性已经得到了国外许多研究项目和研究者的证明，模型已经广泛地应用到美国国家项目HUMUS(Hydrological Unit Modeling of United States)、大的区域性项目(如：NOAA`s Coastal Assessment Framework)和许多不同尺度的研究项目中，研究内容涉及流域的水平衡、河流流量预测和非点源污染控制评价等诸多方面。

我国20世纪90年代中后期，逐渐加强面源输移过程及入河贡献的定量化研究，不断引入国外面源污染模型，并在其基础上结合水文模型以及GIS技术，评价面源污染的发生区域、污染物的迁移变化以及受纳水体响应等。但鉴于流域水文特征的复杂性及基础性监测数据的缺乏，已进行的研究或考虑的因素相对较少、或模拟的介质体系相对较为单一。以淮河沙颍河流域为例，该流域位于半湿润区域，降雨量空间分布极其不均匀，水文分区复杂，农业面源污染突出。不管使用何种模型对其面源污染负荷进行评估，都必须考虑半湿润地区的水文过程特征对面源污染输移过程的影响，以精确反映面源污染物在水循环中的变化过程，最终获取准确的模拟计算结果。

SWAT（Soil and water assessment tools）模型被广泛的应用于国内外面源污染重点源区的识别以及负荷的估算等领域。但该模型中的水文模块采用SCS曲线数法来进行降雨径流模拟。SCS曲线方程自20世纪50年代逐渐得到广泛使用，是根据美国小流域降水与径流关系20多年的研究成果总结得到的经验模型。很明显，该模型能很好地表征美国流域水文循环特征，但是否适用于我国的水文模拟，还有待进一步的实践和论证，此外，该模型能否准确的反映半湿润地区的水文过程变化规律，目前尚未有研究给出可信的答案。

发明内容

1.发明要解决的问题

针对现有技术中对于我国半湿润地区的流域面源污染负荷的研究模型尚不能准确反映我国半湿润地区的水文过程变化规律的不足，本发明提供了一种适用于半湿润地区的流域面源污染负荷的确定方法，该方法能够有效地模拟半湿润地区淮河流域典型区域的面源污染负荷估算，针对半湿润地区特有水文特征，给出合理科学的计算结果，提高了模型应用的普及性和有效性，能更好地服务于水环境保护与水污染控制等水环境研究项目。

2.技术方案

发明原理：本发明通过建立能反映半湿润区水文过程变化规律的水文模型，并将其与SWAT模型中的土壤侵蚀模块以及污染负荷模块进行嵌套，以改进SWAT模型对水文过程特征的描述，进而获取更为精确的污染负荷模拟计算结果。

本发明的技术方案如下：

一种半湿润区流域面源污染负荷的确定方法，其步骤为：

（A）收集整理水文气象和水环境历史观测资料以及野外实验资料，具体包括水文气象数据、水环境数据、社会经济数据以及水利工程调度数据；

（B）根据DEM（Digital Elevation Model）数据、土地利用数据以及土壤类型数据划分子流域和水文响应单元；其中，子流域的划分是基于DEM数据利用D8算法设定流向，进而获取坡度、集水面积以及水系拓扑关系，最终完成子流域的划定；水文响应单元的生成是根据各栅格土地利用类型以及土壤类型属性，进而将具有相同属性的栅格视为具有相同水文响应的单元；

（C）改进水文模块，包括：先设定上层张力水容量WUM、下层张力水容量WLM、深层张力水容量WDM、蒸散发折算系数K、张力水蓄水容量曲线指数B、自由水蓄水库容量SM、壤中流出流系数KI、地下径流出流系数KG模型参数的初始值，建立基于栅格的产流模块，将径流量进行水源划分，针对不同的径流成分进行流域汇流计算，获取流量过程时间序列；

a)产流模块的构建，半湿润地区多以蓄满产流机制为主；以栅格作为单元，根据蓄满产流的概念，采用下式进行产流计算：

当PE≤0时，不产流，即产流量R=0；

当PE＞0时产流，若PE+A＜W’_mm；

$R=\mathrm{PE}-\mathrm{WM}+W_0+\mathrm{WM}{[1-\frac{\mathrm{PE}+A}{{W^{\′}}_m}]}^{B+1}---\left(1\right)$

若PE+A≥W’_mm

R＝PE-(WM-W₀)                            （2）

其中PE为净雨量，W’_mm为栅格单元内最大的点蓄水容量，WM为栅格单元平均蓄水容量，W₀栅格单元初始平均蓄水量；A为初始土壤含水量对应的纵坐标，B为土壤蓄水容量曲线指数，R为产流量，W’_m为单点最大蓄水容量。

实际蒸散发量采用三层蒸散发模型进行计算。

b)水源划分是引入自由水蓄水库的概念，假定自由水蓄水库设置两个出口，其出流系数分别为KI和KG，产流量R进入自由水库内，通过两个出流系数和溢流的方式把它分成地面径流RS、壤中流RI和地下径流RG，计算流程如下：

当PE＋AU＜SMMF时，地面径流RS为：

$\mathrm{RS}=\mathrm{FR}*\{\mathrm{PE}-\mathrm{SMF}+S+\mathrm{SMF}{[1-\frac{(\mathrm{PE}+\mathrm{AU})}{\mathrm{SMMF}}]}^{\mathrm{EX}+1}\}---\left(3\right)$

当PE＋AU≥SMMF时

RS＝FR*（PE+S-SMF）                        (4)

壤中流RI、地下径流RG的计算式为：

RI＝[((PE+S)*FR-RS]*KI                     (5)

RG＝[(PE+S)*FR-RS]*KG                        (6)

其中，AU为初始含水量对应的纵坐标值，FR为产流面积，PE为净雨量，S为自由水蓄水深，SMF为自由水平均蓄水容量，SMMF为最大点自由水蓄水容量，KI为壤中流出流系数，KG为地下径流出流系数，FR为产流面积百分比，EX为自由水蓄水容量曲线指数；

c)汇流计算分为坡地汇流计算和河网汇流计算；汇流模块又包含坡面流演算以及壤中流和地下水流演算，坡面流演算采用运动波的方法，壤中流和地下水流演算采用串联的线性水库的方法；河网汇流计算采用马斯京根-康吉方法；

（D）基于改进水文模块计算得到的水文过程结果，先设定土壤侵蚀因子、植被覆盖和管理因子、保持措施因子、地形因子和粗碎屑因子模型参数初始值，然后采用MUSLE方程进行土壤侵蚀量的计算，公式如下：

m_sed＝11.8×(Q_surf·q_peak·A_hru)^0.56·K_USLE·C_USLE·P_USLE·LS_USLE·CFRG  （7）

式中：m_sed为土壤侵蚀量，t；Q_surf为地表径流，mm/h；q_peak为洪峰径流，m³/s；A_hru为水文响应单元（HRU）的面积，hm²；K_USLE为土壤侵蚀因子；C_USLE为植被覆盖和管理因子；P_USLE为保持措施因子；LS_USLE为地形因子；CFRG为粗碎屑因子；

（E）基于水文过程以及土壤侵蚀量的计算结果，首先设定氮渗透系数NPERCO、生物混合效率系数BIOMIX、氨氮氧化为亚硝氮的速率常数BC1、亚硝氮氧化为硝氮的速率常数BC2和有机氮氧化为氨氮的速率常数BC3水污染负荷模块的模型参数的初始值，然后采用水污染负荷模块，包括溶解态氮(硝态氮)污染负荷模块、吸附态氮(有机氮)污染负荷模块、溶解态磷污染负荷模块、吸附态磷(有机磷和矿物质磷)污染负荷模块以及河道中各种形态氮的转化模块进行面源污染负荷量的计算，溶解态氮(硝态氮)污染负荷估算方程如下式：

${\mathrm{\ρ}}_{\mathrm{NO}3,\mathrm{mobile}}=\frac{{\mathrm{\ρ}}_{\mathrm{NO}3\mathrm{ly}}\·\exp \left[\frac{-w_{\mathrm{mobile}}}{(1-{\mathrm{\θ}}_e){\mathrm{SAT}}_{\mathrm{ly}}}\right]}{w_{\mathrm{mobile}}}---\left(8\right)$

式中为自由水中硝态氮浓度(以N计)，kg/mm；为土壤中硝态氮的量（以N计)，kg/hm²；W_mobile为土壤中自由水的量，mm；θ_e为孔隙度；STA_ly为土壤饱和含水量；

吸附态氮(有机氮)污染负荷估算方程如下式：

${\mathrm{\ρ}}_{\mathrm{orgNsurf}}=0.001\×{\mathrm{\ρ}}_{\mathrm{orgN}}\·\frac{m}{A_{\mathrm{hru}}}\·{\mathrm{\ϵ}}_N---\left(9\right)$

式中：ρ_orgNsurf为有机氮流失量(以N计)，kg/hm²；ρ_orgN为有机氮在表层（10mm）土壤中的浓度（以N计），kg/t；m为土壤流失量，t；A_hru水表响应单元的面积，hm²，ε_N为氮富集系数，氮富集系数是随土壤流失的有机氮浓度和土壤表层有机氮浓度的比值；

溶解态磷污染负荷估算方程如下式：

$P_{\mathrm{surf}}=\frac{P_{\mathrm{solution},\mathrm{surf}}\·Q_{\mathrm{surf}}}{{\mathrm{\ρ}}_b\·h_{\mathrm{surf}}\·k_{d,\mathrm{surf}}}---\left(10\right)$

式中：P_surf为通过地表径流流失的溶解态磷（以P计)，kg/hm²；P_{solution，surf}为土壤中(表层10mm)溶解态磷（以P计)，kg/hm²；Q_surf为地表径流量；ρ_b为土壤溶质密度，（以P计)，kg/m³；h_surf为表层土壤深度，mm；k_d，surf为土壤磷分配系数，表层土壤（10mm）中溶解态磷的浓度和地表径流中溶解态磷浓度的比值；

吸附态磷(有机磷和矿物质磷)污染负荷估算方程如下式：

$m_{\mathrm{Psurf}}=0.001\×{\mathrm{\ρ}}_P\·\frac{m}{A_{\mathrm{hru}}}\·{\mathrm{\ϵ}}_P---\left(11\right)$

式中：m_Psurf为有机磷流失量(以P计)，kg/hm²；ρ_P为有机磷氮在表层（10mm）土壤中的浓度（以P计），kg/t；m为土壤流失量，t；A_hru为水表响应单元的面积，hm²，ε_P为磷富集系数；

河道中各种形态氮的转化量估算，计算流程如下：

一天内有机氮的变化可以用下式来描述：

Δρ_orgNstr＝(α₁·ρ_a·ρ_algae-β_N,3·ρ_orgNstr-σ₄ρ_orgNstr)·TT    （12）

式中：Δρ_orgNstr为有机氮浓度的变化量(以N计)，mg/L：α₁为藻类生物量中的氮含量(以N计)，mg/mg；ρ_a为当地藻类的死亡速度，d^-1，；ρ_algae为一天中开始时藻类生物量的含量，mg/L；β_N，3为有初氮转化为氨的速度常数，d^-1；ρ_orgNstr为一天中开始时有机氮的含量(以N计)，mg/L；σ₄为有机氮的沉淀系数，d^-1；TT为在该河段的运动时间，d；

一天内氨的变化可以用下式来描述：

${\mathrm{\Δ\ρ}}_{{\mathrm{NH}}_{4\mathrm{str}}}=({\mathrm{\β}}_{N,3}\·{\mathrm{\ρ}}_{\mathrm{orgNstr}}-{\mathrm{\β}}_{N,1}\·{\mathrm{\ρ}}_{{\mathrm{NH}}_{4\mathrm{str}}}+\frac{{\mathrm{\σ}}_3}{100\×h}-{\mathrm{fr}}_{{\mathrm{NH}}_4}\·{\mathrm{\α}}_1\·{\mathrm{\μ}}_a\·{\mathrm{\ρ}}_{\mathrm{algae}})\·\mathrm{TT}---\left(13\right)$

式中：为氨含量的变化量(以N计)，mg/L；β_N,3为有机氮转化为氨氮的速度常数，d^-1；ρ_orgNstr为一天中开始时有机氮的含量(以N计)，mg/L；β_N,1为氨氮的氧化速度常数，d^-1；为一天开始时氨氮的含量(以N计)，mg/L；σ₃为沉淀物的氨释放速度(以N计)，mg/(m³d)；h为河道中的水深，m；为藻类的氨氮吸收系数；α₁为藻类生物量中氮含量(以N计)，mg/mg；μ_a为藻类的生长速度，d^-1；ρ_algae为一天开始时藻类生物量的含量，mg/L；TT为在该河段的运动时间，d;

一天内亚硝酸盐的变化用下式来描述：

${\mathrm{\Δ\ρ}}_{{\mathrm{NO}}_{2\mathrm{str}}}=({\mathrm{\β}}_{N,1}\·{\mathrm{\ρ}}_{{\mathrm{NH}}_{4\mathrm{str}}}-{\mathrm{\β}}_{N,2}\·{\mathrm{\ρ}}_{{\mathrm{NO}}_{2\mathrm{str}}})\·\mathrm{TT}---\left(14\right)$

式中：为亚硝酸盐的改变量(以N计)，mg/L；β_N,1为氨氮的生物氧化速度常数d^-1；为一天中开始时的氨氮含量(以N计)，mg/L；β_N，2为由亚硝酸盐到硝酸盐的氧化速度常数，d^-1；为一天开始时亚硝酸盐的含量(以N计)，mg/L；TT为在河段中的运动时间，d;

一天内硝酸盐的变化量为：

${\mathrm{\Δ\ρ}}_{{\mathrm{NO}}_{3\mathrm{str}}}=[{\mathrm{\β}}_{N,2}\·{\mathrm{\ρ}}_{{\mathrm{NO}}_{2\mathrm{str}}}-(1-{\mathrm{fr}}_{{\mathrm{NH}}_4})\·{\mathrm{\α}}_1\·{\mathrm{\μ}}_a\·{\mathrm{\ρ}}_{\mathrm{algae}}]\·\mathrm{TT}---\left(15\right)$

式中：为硝酸盐的改变量(以N计)，mg/L；β_N,2为氨氮的生物氧化速度常数d^-1；为一天中开始时的氨氮含量(以N计)，mg/L；为藻类的氨氮吸收系数；α₁为藻类生物量中氮含量(以N计)，mg/mg；μ_a为藻类的生长速度，d^-1；ρ_algae为一天开始时藻类生物量的含量，mg/L；TT为在该河段的运动时间，d;

（F）将计算得到的水文过程和面源污染负荷量与实测数据进行对比，确定水文模块、MUSLE方程以及水污染负荷计算模块的参数值，完成模型的建立，其中，模型参数的确定是分两个层次进行，首先利用水文气象资料对水文模块的参数进行率定，然后利用水质数据对水污染负荷模块进行率定，率定的方法就是先设定上层张力水容量WUM、下层张力水容量WLM、深层张力水容量WDM、蒸散发折算系数K、张力水蓄水容量曲线指数B、自由水蓄水库容量SM、壤中流出流系数、地下径流出流系数、氮渗透系数NPERCO、生物混合效率系数BIOMIX、氨氮氧化为亚硝氮的速率常数BC1、亚硝氮氧化为硝氮的速率常数BC2和有机氮氧化为氨氮的速率常数BC3等这样一组水文模块、MUSLE方程和水污染负荷模块参数的初始值，然后通过（C）-（E）步骤的计算，得到模拟的结果，并与实测数据进行比较，当相对误差不超过20％时，就认为值是合理可靠的，如果大于20%，继续假定另一组参数值，重复上述步骤（C）-（E）的计算过程，直到参数值符合要求为止；参数的验证方法与率定相同，只是选取与率定期不同期的数据；土壤特性数据、土壤中的氮以及有机质含量参数，通过野外试验资料获取。

3.有益效果

相比于现有技术，本发明的优点在于：

（1）本发明基于栅格单元，构建了适用于半湿润地区的水文模拟方法，并将上述方法与土壤侵蚀模块以及面源污染负荷计算模块紧密嵌套，以充分考虑水文过程对面源污染负荷计算的影响，提高模拟结果的可靠性，为面源污染控制提供管理工具；

（2）针对大多数面源污染负荷模型都是基于经验水文模型以反映水文过程对面源污染负荷输移的影响，难以反映流域内的真实情况，本发明构建的半湿润地区面源污染负荷数值模拟方法，通过对面源污染负荷计算模型中的水文模块的改进，以准确表征水文循环过程对面源污染物输移规律的影响，提高了模型的计算结果精度；

（3）实际应用表明，本发明能够有效地模拟半湿润地区淮河流域典型区域的面源污染负荷估算，针对半湿润地区特有水文特征，给出合理科学的计算结果，提高了模型应用的普及性和有效性，能更好地服务于水环境保护与水污染控制等水环境研究项目。

说明书附图

图1为本发明的系统结构示意图；

图2为中牟站2001-2002年实测流量和改进后水文模块模拟月流量过程对比图；

图3为中牟站2004-2005年实测流量和改进后水文模块模拟月流量过程对比图；

图4为中牟站2007年实测流量和改进后水文模块模拟月流量过程对比图；

图5为中牟站2008年实测流量和改进后水文模块模拟月流量过程对比图；

图6为扶沟站2001-2002年实测流量和改进后水文模块模拟月流量过程对比图；

图7为扶沟站2004-2005年实测流量和改进后水文模块模拟月流量过程对比图；

图8为扶沟站2007年实测流量和改进后水文模块模拟月流量过程对比图；

图9为扶沟站2008年实测流量和改进后水文模块模拟月流量过程对比图；

图10为中牟站2008年氨氮实测值与改进后面源污染负荷模型模拟值对比图；

图11为扶沟站2008年氨氮实测值与改进后面源污染负荷模型模拟值对比图；

图12为贾鲁河及颍河各水文响应单元面源负荷产生量；

图13为中牟站2001-2002年实测流量和原模型水文模块模拟月流量过程对比图；

图14为中牟站2004-2005年实测流量和原模型水文模块模拟月流量过程对比图；

图15为扶沟站2001-2002年实测流量和原模型水文模块模拟月流量过程对比图；

图16为扶沟站2004-2005年实测流量和原模型水文模块模拟月流量过程对比图；

图17为中牟站2008年氨氮实测值与原模型模拟值对比图；

图18为扶沟站2008年氨氮实测值与原模型模拟值对比图。

具体实施方式

以下通过实施例对本发明作进一步描述。

实施例

选取沙颍河流域贾鲁河子流域为应用区域，该流域面积为5900km²，控制水文站为扶沟闸站，流域内分布有8个雨量站。结合附图1，具体说明本例实施方式，包括前期数据准备、改进模型的率定和验证以及与原模型结果的比较。

其中，如图1中的基础资料准备步骤所示，流域面源污染负荷数值模拟模型建立的数据准备包括：

a、土地利用类型图

为了适应模型的要求，对土地利用进行重新编码，划分为耕地、林地、园地、草地、水域、农村居民点、城镇用地和裸地8类，土地利用图通过arcgis将Shape格式，转化成网格为100＊100m的栅格格式，得到区域内土地利用类型图，不同类型土壤类型面积比重如表1所示。

表1贾鲁河流域土地利用面积分布情况

序号	土地利用类型	面积	占比例	代码
					1	耕地	4984.03	78.14	AGRL
2	林地	411.45	6.45	FRST
					3	园地	49.80	0.78	ORCD
4	草地	93.31	1.46	PAST
					5	水域	57.04	0.89	WATR
6	城镇用地	255.43	4.00	URHD
					7	农村居民用地	526.17	8.25	URLD
8	裸地	0.86	0.01	BALD

b、土壤类型图

根据收集的土壤图（格式为shape文件），并将其直接转换为栅格的格式，其栅格大小为100＊100m，贾鲁河子流域土壤类型分布情况见表2。

表2贾鲁河子流域土壤类型面积分布情况及模型代码

c、气象数据库

气象数据分为以下两类。气象月平均值：从气象资料界面录入以构建模型所需站点数据，用来模拟生成没有观测资料的气象因子；气象站点的实测数据：按照模型的格式要求建立.dbf或.txt格式的文件进行输入，并建立.dbf格式的链接表将气象资料导入模型。

d、土壤数据库

土壤属性数据主要包括两类：物理性数据库和化学性数据。以上参数通过对采集的研究区域土壤样品进行分析，直接获取。

表3土壤物理化学分析指标一览表

根据表3中测定值并参考河南省土壤地理，确定土壤容重、粒径分级及有机质的含量，其余参数的计算参照《SWAT模型土壤数据库建立方法》中提供的方法利用SPAW模型进行水文组分划分、饱和导水率、土壤侵蚀因子K等物理参数的计算与转换。

土壤化学属性决定着土壤中营养物质的赋存状态和含量，主要包括土壤的硝态氮、有机氮和有机磷的浓度，用来给模型进行污染模拟赋初始值。根据所测定土壤的化学参数，整理得部分土壤化学参数见表4。

表4土壤化学参数表

如图1中的模型构建步骤中所示，其中模型参数率定和验证包括：

a、水文模块参数率定与验证

首先，基于DEM数据划分子流域，按照土地利用类型和土壤类型生成水文响应单元；然后设定模型参数（上层张力水容量WUM、下层张力水容量WLM、深层张力水容量WDM、蒸散发折算系数K、张力水蓄水容量曲线指数B、自由水蓄水库容量SM、壤中流出流系数KI和地下径流出流系数KG）的初始值；将2001~2002年和2004-2005年的水文气象资料作为各响应单元的输入，利用公式（1）和（2）计算产流量，利用三层蒸散发模型计算实际蒸散发量，采用公式（4）-（6）将产流量划分为地表径流、壤中流以及地下径流，采用一维运动波方法将地表径流汇至河道，采用串联的线性水库将壤中流和地下径流汇至河道，河道汇流计算采用马斯京根-康吉方法，最终将流量演算至流域出口断面；最后，将计算得到流域出口断面流量与实测流量进行对比，计算各年的径流相对误差，如果年径流相对误差小于20%，认为设定的参数值为所求，如果大于20%，继续假定另一组模型参数，重复上述步骤的计算，直到参数值符合要求为止。

利用2007-2008年的水文气象资料验证上述率定结果，验证过程中将参数值设定为率定获得结果，计算过程与参数率定过程中的流程相同。本例最终获得的水文模块参数结果见表8。表5~表6以及附图2~图9展示了实测流量与计算流量的比较情况，图表中数据表明参数基本符合要求，能反映研究流域的水文特征。

表5改进后的水文模型中牟站日均流量模拟值与实测值对比情况表

年份	实测值(m3/s)	模拟值(m3/s)	相对误差(%)
				2001	9.26	8.70	-6.05
2002	9.36	10.19	8.87
				2004	13.03	13.03	0.00
2005	20.14	22.51	11.77
				2007	11.81	11.82	0.08
2008	17.44	17.42	-0.11

表6改进后的水文模型扶沟站日均流量模拟值与实测值对比情况表

年份	实测值(m3/s)	模拟值(m3/s)	相对误差(%)
				2001	6.79	7.88	16.01

2002	5.46	7.22	32.25
				2004	20.31	20.05	-1.28
2005	19.61	19.57	-0.20
				2007	15.82	15.27	-3.46
2008	23.91	21.58	-9.76

b、污染负荷模块参数率定与验证

以2008年为基准年，基于率定好的改进的水文模块输出的水量计算结果，设定MUSLE方程和污染负荷模型参数(土壤侵蚀因子、植被覆盖和管理因子、保持措施因子、地形因子和粗碎屑因子模型以及氮渗透系数NPERCO、生物混合效率系数BIOMIX、氨氮氧化为亚硝氮的速率常数BC1、亚硝氮氧化为硝氮的速率常数BC2和有机氮氧化为氨氮的速率常数BC3)初始值，利用公式（7）计算土壤侵蚀量，基于已知的水文过程和土壤侵蚀量利用公式（8）-（15）计算面源污染负荷量，将计算得到的面源污染负荷量与水环境资料得到的污染负荷量进行对比，计算污染负荷量相对误差，如果该误差小于20%，模型参数（上层张力水容量WUM、下层张力水容量WLM、深层张力水容量WDM、蒸散发折算系数K、张力水蓄水容量曲线指数B、自由水蓄水库容量SM、壤中流出流系数KI、地下径流出流系数KG、氮渗透系数NPERCO、生物混合效率系数BIOMIX、氨氮氧化为亚硝氮的速率常数BC1、亚硝氮氧化为硝氮的速率常数BC2和有机氮氧化为氨氮的速率常数BC3）即为所求，如果大于20%，重新设定一组参数值，重复水文过程、土壤侵蚀量以及水污染负荷的计算过程，直到参数符合要求。参数率定结果见表8，表7及附图10-11展示了2008年氨氮模拟值与实测值比较情况，结果表明率定得到的参数基本可靠。

表7改进后的面源污染负荷模型2008年氨氮模拟值与实测值比较  单位：吨

站点	实测值(t/年)	模拟值(t/年)	相对误差(%)
				中牟站	4439.8	4679.08	5.39
扶沟站	4852.3	4552.9	-6.17

c、模型参数值的确定

通过对水文模块和污染负荷模块参数进行率定和验证后，模型参数的最终取值见表8。参数一旦确定，研究流域改进的面源污染负荷模型建立完成。

表8改进的面源污染负荷模拟模型校准参数值

参数	模拟过程	典型范围	模型最终值
				上层张力水容量WUM	径流	5-30	20
下层张力水容量WLM	径流	50-100	70
				深层张力水容量WDM	径流	20-50	30
蒸散发折算系数K	径流	0.5-1.5	1.2
				张力水蓄水容量曲线指数B	径流	0.1-0.5	0.3
自由水蓄水库容量SM	径流	10-50	18
				壤中流出流系数	径流	0.4	0.1-0.7
地下径流出流系数	径流	0.3	0.1-0.7

壤中流消退系数	径流	0.15	0-0.9
				地下径流消退系数	径流	0.98	0.95-0.99
氮渗透系数NPERCO	氮	0.2-0.6	0.40
				BIOMIX	氮	0-1.0	0.40
BC1	氮	0.1-1.0	0.10
				BC2	氮	0.2-2.0	1.10
BC3	氮	0.2-0.4	0.40

如附图1的对比分析步骤所示，利用改进后的、适用于半湿润地区的面源污染负荷数值模拟方法对贾鲁河流域的面源污染负荷进行估算和分析，为了展现该模型对面源污染负荷模拟的优势，将该模型的计算结果与原模型结果进行对比，对比方法是通过分别将原模型和改进模型的模拟结果与实测值进行比较，与实测值越接近，就表明模型效果越好。

结果表明，在现状情景设计水文年型下（见表9），贾鲁河流域氨氮面源入河污染负荷1128.64吨，氨氮入河污染负荷面源贡献率约为12.28％；颍河流域氨氮面源入河污染负荷1695.36吨，氨氮入河污染负荷面源贡献率约为35.40％，颍河氨氮面源入河污染负荷大于贾鲁河。贾鲁河及颍河各水文响应单元面源负荷产生量统计见附图12，其中，图12（a）为可溶性磷总量，图12（b）为有机磷总量，图12（c）为硝氮总量，图12（d）为有机氮总量，图12（e）为颍河流域可溶磷总量，图12（f）为颍河流域有机磷总量，图12（g）为颍河流域硝氮总量，图12（h）为颍河流域有机氮总量。表10-11展示了原模型中SCS曲线模型计算得到的径流量（附图13-16），与改进后水文模块计算得到的结果（表5-6、附图2-3以及附图6-7）进行比较后，发现改进后的模块径流相对误差大幅降低；2001-2001年以及2004-2005年期间最大提高幅度达到53.33％，最小提高幅度为1.51％，中牟站年平均提高幅度为7.2％，扶沟站年平均提高幅度高达25.5％。表12和附图17-18展示了原模型计算得到不同控制站的面源污染负荷量，与改进后面源污染负荷模型计算结果相比，中牟站的输出结果的相对误差由14.57％降低到5.39％，扶沟站输出结果相对误差由-9.16％提高至-6.17％。上述结果充分表明改进后模型相比原模型能更好的描述半湿润地区的水文循环特征，能充分凸显研究区域水文特性对面源污染输移过程的影响，获得更为科学合理的计算结果，为行政部门水环境管理与决策提供更有力的支撑，是十分科学实用的工具。

表92008年贾鲁河流域非点源污染负荷  单位：吨

河流	有机氮	硝酸盐氮	氨氮	有机磷	溶解态磷
						贾鲁河	1003.77	1237.47	1128.64	370.10	2.93
颍河	2522.90	11711.77	1695.36	765.89	34.85

表10SWAT模型中牟站日均流量模拟值与实测值对比情况表

年份	实测值(m3/s)	模拟值(m3/s)	相对误差(%)
				2001	9.26	8.36	-9.72

2002	9.36	10.31	10.10
				2004	13.03	13.23	1.51
2005	14.98	20.10	34.10

表11SWAT模型扶沟站日均流量模拟值与实测值对比情况表

年份	实测值(m3/s)	模拟值(m3/s)	相对误差(%)
				2001	6.79	7.13	4.90
2002	5.46	10.13	85.58
				2004	20.31	15.05	-25.91
2005	19.61	22.07	12.55

表12原模型2008年氨氮模拟值与实测值比较  单位：吨

站点	实测值(t/年)	模拟值(t/年)	相对误差(%)
				中牟站	4439.8	5086.5	14.57
扶沟站	4852.3	4407.84	-9.16

Claims (1)

1.一种半湿润区流域面源污染负荷的确定方法，其步骤为：

(A)收集整理水文气象和水环境历史观测资料以及野外实验资料，具体包括水文气象数据、水环境数据、社会经济数据以及水利工程调度数据；

(B)根据DEM数据、土地利用数据以及土壤类型数据划分子流域和水文响应单元；其中，子流域的划分是基于DEM数据利用D8算法设定流向，进而获取坡度、集水面积以及水系拓扑关系，最终完成子流域的划定；水文响应单元的生成是根据各栅格土地利用类型以及土壤类型属性，进而将具有相同属性的栅格视为具有相同水文响应的单元；

(C)改进水文模块，包括：先设定上层张力水容量WUM、下层张力水容量WLM、深层张力水容量WDM、蒸散发折算系数K、张力水蓄水容量曲线指数B、自由水蓄水库容量SM、壤中流出流系数KI、地下径流出流系数KG的初始值，建立基于栅格的产流模块，将径流量进行水源划分，针对不同的径流成分进行流域汇流计算，获取流量过程时间序列；

a)产流模块的构建，以栅格作为单元，根据蓄满产流的概念，采用下式进行产流计算：

当PE≤0时，不产流，即产流量R＝0；

当PE＞0时产流，若PE+A＜W’_mm；

$R=\mathrm{PE}-\mathrm{WM}+W_0+\mathrm{WM}{[1-\frac{\mathrm{PE}+A}{{W^{\′}}_m}]}^{B+1}---\left(1\right)$

若PE+A≥W’_mm

R＝PE-(WM-W₀)   (2)

其中PE为净雨量，W’_mm为栅格单元内最大的点蓄水容量，WM为栅格单元平均蓄水容量，W₀栅格单元初始平均蓄水量；A为初始土壤含水量对应的纵坐标，B为土壤蓄水容量曲线指数，R为产流量，W’_m为单点最大蓄水容量；

b)水源划分是引入自由水蓄水库的概念，假定自由水蓄水库设置两个出口，其壤中流出流系数、地下径流出流系数分别为KI和KG，产流量R进入自由水库内，通过两个出流系数和溢流的方式把产流量R分成地面径流RS、壤中流RI和地下径流RG，计算流程如下：

当PE+AU<SMMF时，地面径流RS为：

$\mathrm{RS}=\mathrm{FR}*\{\mathrm{PE}-\mathrm{SMF}+S+\mathrm{SMF}{[1-\frac{(\mathrm{PE}+\mathrm{AU})}{\mathrm{SMMF}}]}^{\mathrm{EX}+1}\}---\left(3\right)$

当PE+AU≥SMMF时

RS＝FR*(PE+S-SMF)   (4)

壤中流RI、地下径流RG的计算式为：

RI＝[(PE+S)*FR-RS]*KI   (5)

RG＝[(PE+S)*FR-RS]*KG   (6)

其中，AU为初始含水量对应的纵坐标值，FR为产流面积，PE为净雨量，S为自由水蓄水深，SMF为自由水平均蓄水容量，SMMF为最大点自由水蓄水容量，KI为壤中流出流系数，KG为地下径流出流系数，EX为自由水蓄水容量曲线指数；

c)汇流计算分为坡地汇流计算和河网汇流计算；汇流模块又包含坡面流演算以及壤中流和地下水流演算，坡面流演算采用运动波的方法，壤中流和地下水流演算采用串联的线性水库的方法；河网汇流计算采用马斯京根-康吉方法；

(D)基于改进水文模块计算得到的水文过程结果，先设定土壤侵蚀因子、植被覆盖和管理因子、保持措施因子、地形因子和粗碎屑因子模型参数初始值，然后采用MUSLE方程进行土壤侵蚀量的计算，公式如下：

m_sed＝11.8×(Q_surf·q_peak·A_hru)^0.56·K_USLE·C_USLE·P_USLE·LS_USLE·CFRG   (7)

式中：m_sed为土壤侵蚀量，t；Q_surf为地表径流，mm/h；q_peak为洪峰径流，m³/s；A_hru为水文响应单元的面积，hm²；K_USLE为土壤侵蚀因子；C_USLE为植被覆盖和管理因子；P_USLE为保持措施因子；LS_USLE为地形因子；CFRG为粗碎屑因子；

(E)基于水文过程以及土壤侵蚀量的计算结果，首先设定氮渗透系数NPERCO、生物混合效率系数BIOMIX、氨氮氧化为亚硝氮的速率常数BC1、亚硝氮氧化为硝氮的速率常数BC2和有机氮氧化为氨氮的速率常数BC3的参数初始值，然后采用水污染负荷模块，包括溶解态氮污染负荷模块、吸附态氮污染负荷模块、溶解态磷污染负荷模块、吸附态磷污染负荷模块以及河道中各种形态氮的转化模块，进行面源污染负荷量的计算，溶解态氮污染负荷估算方程如下式：

${\mathrm{\&rho;}}_{\mathrm{NO}3,\mathrm{mobile}}=\frac{{\mathrm{\&rho;}}_{{\mathrm{NO}}_3\mathrm{ly}}\&CenterDot;\exp \left[\frac{-w_{\mathrm{mobile}}}{(1-{\mathrm{\&theta;}}_e){\mathrm{SAT}}_{\mathrm{ly}}}\right]}{w_{\mathrm{mobile}}}---\left(8\right)$

式中为自由水中硝态氮浓度，kg/mm；为土壤中硝态氮的量，kg/hm²；W_mobile为土壤中自由水的量，mm；θ_e为孔隙度；STA_ly为土壤饱和含水量；

吸附态氮污染负荷估算方程如下式：

${\mathrm{\&rho;}}_{\mathrm{orgNsurf}}=0.001\&times;{\mathrm{\&rho;}}_{\mathrm{orgN}}\&CenterDot;\frac{m}{A_{\mathrm{gru}}}\&CenterDot;{\mathrm{\&epsiv;}}_N---\left(9\right)$

式中：ρ_orgNsurf为有机氮流失量，kg/hm²；ρ_orgN为有机氮在表层土壤中的浓度，以N计，kg/t；m为土壤流失量，t；A_hru水表响应单元的面积，hm²，ε_N为氮富集系数，氮富集系数是随土壤流失的有机氮浓度和土壤表层有机氮浓度的比值；

溶解态磷污染负荷估算方程如下式：

$P_{\mathrm{surf}}=\frac{P_{\mathrm{solution},\mathrm{surf}}\&CenterDot;Q_{\mathrm{surf}}}{{\mathrm{\&rho;}}_b\&CenterDot;h_{\mathrm{surf}}\&CenterDot;k_{d,\mathrm{surf}}}---\left(10\right)$

式中：P_surf为通过地表径流流失的溶解态磷，kg/hm²；P_{solution，surf}为土壤中溶解态磷，kg/hm²；Q_surf为地表径流量；ρ_b为土壤溶质密度，kg/m³；h_surf为表层土壤深度，mm；k_d，surf为土壤磷分配系数，表层土壤中溶解态磷的浓度和地表径流中溶解态磷浓度的比值；

吸附态磷污染负荷估算方程如下式：

$m_{\mathrm{Psurf}}=0.001\&times;{\mathrm{\&rho;}}_P\&CenterDot;\frac{m}{A_{\mathrm{hru}}}\&CenterDot;{\mathrm{\&epsiv;}}_P---\left(11\right)$

式中：m_Psurf为有机磷流失量，kg/hm²；ρ_P为有机磷氮在表层土壤中的浓度，kg/t；m为土壤流失量，t；A_hru为水表响应单元的面积，hm²，ε_P为磷富集系数；

河道中各种形态氮的转化量估算，计算流程如下：

一天内有机氮的变化用下式来描述：

Δρ_orgNstr＝(α₁·ρ_a·ρ_algae-β_N,3·ρ_orgNstr-σ₄ρ_orgNstr)·TT   (12)

式中：Δρ_orgNstr为有机氮浓度的变化量，mg/L：α₁为藻类生物量中的氮含量，mg/mg；ρ_a为当地藻类的死亡速度，d^-1；ρ_algae为一天中开始时藻类生物量的含量，mg/L；β_N,3为有初氮转化为氨的速度常数，d^-1；ρ_orgNstr为一天中开始时有机氮的含量，mg/L；σ₄为有机氮的沉淀系数，d^-1；TT为在河段的运动时间，d；

一天内氨的变化用下式来描述：

${\mathrm{\&Delta;\&rho;}}_{{\mathrm{NH}}_{4\mathrm{str}}}=({\mathrm{\&beta;}}_{N,s}\&CenterDot;{\mathrm{\&rho;}}_{\mathrm{orgNstr}}-{\mathrm{\&beta;}}_{N,1}\&CenterDot;{\mathrm{\&rho;}}_{{\mathrm{NH}}_{4\mathrm{str}}}+\frac{{\mathrm{\&sigma;}}_3}{1000\&times;h}-{\mathrm{fr}}_{{\mathrm{NH}}_4}\&CenterDot;{\mathrm{\&alpha;}}_1\&CenterDot;{\mathrm{\&mu;}}_a\&CenterDot;{\mathrm{\&rho;}}_{\mathrm{algae}})\&CenterDot;\mathrm{TT}---\left(13\right)$

式中：为氨含量的变化量，mg/L；β_N,3为有机氮转化为氨氮的速度常数，d^-1；ρ_orgNstr为一天中开始时有机氮的含量，mg/L；β_N,1为氨氮的氧化速度常数，d^-1；为一天开始时氨氮的含量，mg/L；σ₃为沉淀物的氨释放速度，mg/(m³d)；h为河道中的水深，m；为藻类的氨氮吸收系数；α₁为藻类生物量中氮含量，mg/mg；μ_a为藻类的生长速度，d^-1；ρ_algae为一天开始时藻类生物量的含量，mg/L；TT为在河段的运动时间，d；

一天内亚硝酸盐的变化用下式来描述：

${\mathrm{\&Delta;\&rho;}}_{{\mathrm{NO}}_{2\mathrm{str}}}=({\mathrm{\&beta;}}_{N,1}\&CenterDot;{\mathrm{\&rho;}}_{{\mathrm{NH}}_{4\mathrm{str}}}-{\mathrm{\&beta;}}_{N,2}\&CenterDot;{\mathrm{\&rho;}}_{{\mathrm{NO}}_{2\mathrm{str}}})\&CenterDot;\mathrm{TT}---\left(14\right)$

式中：为亚硝酸盐的改变量，mg/L；β_N,1为氨氮的生物氧化速度常数d^-1；为一天中开始时的氨氮含量，mg/L；β_N,2为由亚硝酸盐到硝酸盐的氧化速度常数，d^-1；为一天开始时亚硝酸盐的含量，mg/L；TT为在河段中的运动时间，d；

一天内硝酸盐的变化量为：

$\mathrm{\&Delta;}{\mathrm{\&rho;}}_{{\mathrm{NO}}_{3\mathrm{str}}}=[{\mathrm{\&beta;}}_{N,2}\&CenterDot;{\mathrm{\&rho;}}_{{\mathrm{NO}}_{2\mathrm{str}}}-(1-{\mathrm{fr}}_{{\mathrm{NH}}_4})\&CenterDot;{\mathrm{\&alpha;}}_1\&CenterDot;{\mathrm{\&mu;}}_a\&CenterDot;{\mathrm{\&rho;}}_{\mathrm{algae}}]\&CenterDot;\mathrm{TT}---\left(15\right)$

式中：为硝酸盐的改变量，mg/L；β_N,2为氨氮的生物氧化速度常数d^-1；为一天中开始时的氨氮含量，mg/L；为藻类的氨氮吸收系数；α₁为藻类生物量中氮含量，mg/mg；μ_a为藻类的生长速度，d^-1；ρ_algae为一天开始时藻类生物量的含量，mg/L；TT为在河段的运动时间，d；

(F)将计算得到的水文过程和面源污染负荷量与实测数据进行对比，确定水文模块、MUSLE方程以及水污染负荷计算模块的参数值，完成模型的建立，其中，模型参数的确定是分两个层次进行，首先利用水文气象资料对水文模块的参数进行率定，然后利用水质数据对水污染负荷模块进行率定，率定的方法就是先设定上层张力水容量WUM、下层张力水容量WLM、深层张力水容量WDM、蒸散发折算系数K、张力水蓄水容量曲线指数B、自由水蓄水库容量SM、壤中流出流系数KI、地下径流出流系数KG、氮渗透系数NPERCO、生物混合效率系数BIOMIX、氨氮氧化为亚硝氮的速率常数BC1、亚硝氮氧化为硝氮的速率常数BC2和有机氮氧化为氨氮的速率常数BC3、MUSLE方程和水污染负荷模块参数的初始值，然后通过(C)-(E)步骤的计算得到模拟的结果，与实测数据进行比较，当相对误差不超过20％时，就认为值是合理可靠的，如果大于20％，继续假定另一组参数值，重复上述步骤(C)-(E)的计算过程，直到参数值符合要求为止；参数的验证方法与率定相同，只是选取与率定期不同期的数据；土壤特性数据、土壤中的氮以及有机质含量参数，通过野外试验资料获取。