本发明涉及一种灌区水田非点源污染负荷估算方法。尤其涉及一种基于水位调控的灌区水田非点源污染负荷估算方法，属于农业环境保护技术。
背景技术：
流域农业非点源污染是指在降水、融雪及农田退水等过程中，部分土壤颗粒物、营养元素及有毒害作用的农药和重金属等污染物随径流从农田进入水体，进而形成的水体污染。而径流作为非点源污染形成的重要驱动因子，在非点源污染生成、运移及转化过程中起着举足轻重的作用。水田作为人工营造的湿地系统，其表层的溢流和排水及地下渗透过程贯穿水稻整个生长过程，同时也势必加剧氮磷等营养元素的流失，威胁脆弱的流域水环境。
传统的田间采样技术虽然能够用于分析田间尺度的水田非点源污染负荷，但该技术成本高，且监测周期较长，在采样区选取方面也具有一定的技术难度。此外，在流域尺度下，非点源污染的输出过程具有普遍的时空异质性规律，因此也限制田间采样法的应用和推广。相较而言，模型技术在流域非点源污染时空序列绘制及风险源区识别方面有着突出的优势。在模型率定和模型验证技术的支持下，模型技术能够直观的对流域非点源污染进行量化分析。但目前主流的流域非点源预测模型(如swat，hspf和annagnps)在灌区的应用过程中也存在着明显的缺陷：(1)主流流域非点源预测模型并未内置水田模块，仅通过径流曲线法、green-ampt入渗模型等将水田简化为旱田处理；(2)未考虑水田的淹没状态，在灌溉发生时，除部分灌溉水按下渗量计算，其余均视为地表径流量；(3)不能模拟水田的水位的年季波动特性，无法在非点源污染与水田水位调控之间建立起联系。
目前的非点源污染的预测方法正在不断发展，但考虑水田实际灌溉特征的预测方法却十分匮乏，例如申请号为cn201410619286.1，名称为一种基于gis的分布式水田氮素径流流失负荷估算方法，公开号cn104361523a，仅利用gis系统、土地利用情景和肥料施用情况构建了简单的输出系模型，用于评价水田非点源污染输出负荷，缺乏足够的模型机理支持，也无助与建立水位调控与非点源污染的对应关系。又如申请号为cn201310232100.2，名称为水田非点源污染流出负荷的观测方法及系统，公开号cn103294923a，构建了一种田间尺度的水田非点源污染分析及估算方法，但因成本问题难以在流域尺度展开。再如申请号为cn201210326701.5，名称为一种基于遥感像元的非点源污染计算方法，公开号为cn102867120a，该方法根据土地利用、植被盖度、土壤属性及氮磷平衡规律估算非点源污染释放规律，但该方法难以表征浸没状态下水田的氮磷释放规律。鉴于此种背景，有必要对现有的流域非点源污染预测模型进行修正，使其能够高效、准确的模拟水位变化条件下的水田非点源污染输出负荷。
技术实现要素：
1、目的：本发明的目的是提供一种基于水位调控的灌区水田非点源污染负荷估算方法，进而可以模拟水田水位变化条件下水田的非点源污染输出规律，有助于精准、有效的灌区的流域非点源污染进行分析。
2、技术方案：本发明一种基于水位调控的灌区水田非点源污染负荷估算方法，其特征在于利用fortran语言对swat模型的壶穴模块进行改良，建立起能够表征水田水循环过程和非点源污染之间关联的水田模块。其具体方法步骤如下：
步骤一：壶穴构型改良
利用fortran77语言对swat模型的壶穴模块的源代码进行编译，并据此构建符合水田灌溉和排水特征的水田模型。本专利对壶穴模型的构型进行了必要改良，将其由原swat模型中的圆锥体设定改良为截面积恒定的多边形结构，并基于改进方法增设了田埂面积的补充算法。田埂面积和改进的水田有效储水面积的计算方法如公式1和2所示。
sr＝shd·fr(1)
sa＝shd·(1-fr)(2)
式中：sr表示田埂面积，ha；fr表示田埂面积与水文响应单元面积之比，％；shd表示水田所处的水文响应单元的总面积，ha；sa表示水田有效储水面积，ha。
步骤二：蒸散发算法的改良
利用编译器对壶穴模块的蒸散发算法进行修正，由原模型提出的基于叶面积指数与最大潜在蒸散发量之间关系的概化算法，改进为根据水田淹没或排空状态，分别计算叶片截留蒸发量、水稻蒸散发量、田面水蒸散发量、土壤蒸散发量的优化算法。当叶片截留量大于水田日潜在蒸散发量时，水田模块将先计算叶片截留蒸发量，见公式(3)。当叶片截留量小于水田日潜在蒸散发量时，其差值由水稻蒸散发量弥补，见公式(4)。当水稻蒸散发量及叶片截留量之和仍小于水田日潜在蒸散发量时，水田模块将根据水位条件，分别进行计算：在浸没状态下，蒸发量的差额由田面水蒸散发量补足，见公式(5)；在排空状态下，蒸发量的差额由土壤蒸散发量补足，见公式(6)。
其中，水田日潜在蒸散发量由penman-monteith方程计算。水稻蒸散发量可根据稻田的叶面积指数，分别由公式(7)和公式(8)计算。当水田处于排空状态下，其土壤蒸散发规律与旱田相近，因此仍沿用原swat模型的计算方法，不进行改动。
ea＝ecanecan>pet(3)
ea＝ep+ecanecan<pet<ecan+ep(4)
ea＝ep+ecan+ewpet>ecan+ep(5)
ea＝ep+ecan+espet>ecan+ep(6)
ep＝(pet-ecan)·l/30≤l≤3(7)
ep＝pet-ecanl>3(8)
式中：pet表示水田日潜在蒸散发量，mm；ea表示水田实际蒸发量，mm；ecan表示叶片截留蒸发量，mm；ep表示水稻蒸散发量，mm；ew表示田面水蒸发量，mm；es表示土壤蒸散发量，mm；l表示叶面积指数。
步骤三：下渗量计算方法的改良
利用编译器对壶穴模型的渗透算法进行编译，由原壶穴模型提出的基于土壤最大持水能力和土壤含水量的计算方法，改进为基于水田储水量与水田日稳定下渗量关系的优化算法。当水田处于浸没状态时，若水田储水量大于水田日稳定下渗量时，水田实际下渗量与水田日稳定下渗量在数值上相等，如公式(9)所示；若水田储水量小于水田日稳定下渗量时，则水田储水量即为水田实际下渗量，如公式(10)。当水田处于排空状态时，其水文循环过程与旱田类似，下渗量的计算方法与原swat模型旱田的下渗量计算方式一致，不做修改。
vinf＝ks·sa＝vmaxvs≥vmax(9)
vinf＝vsvs<vmax(10)
式中：vinf表示水田实际下渗量，m3；ks表示土壤稳定入渗速率，cm/day；sa表示土壤有效储水面积，ha；vs表示水田储水量，m3；vmax表示水田日稳定下渗量。
步骤四：灌溉及排水控制模块的构建
基于水稻的实际设置规程，对水位控制单元进行编译。共设置最低三条水位线用于控制水田的排水及灌溉过程：最低水位、最适水位和最高水位。当稻田水位小于最低水位时，水田模块将触发灌溉功能，直至水位达到最适水位，如公式(11)。在降雨事件中，当水田的水位上升至最高水位时，将触发水田的排水功能，多余的田面水将以溢流的形式排入河网，直至水位回落至最高水位以下，如公式(12)。三条水位线均可设定逐日数据，并以txt文本格式保存，供模型读取。
ri＝(h-pd,i)·sapd,i<hmin(11)
oi＝(pd,i-hmax)·sapd,i>hmax(12)
式中：ri表示灌溉用水量，m3；pd,i表示水田的浸没深度，mm；oi表示水田的溢流量，m3；hmin表示最低水位，mm；h表示最适水位，mm；hmax表示最高水位，mm。
步骤五：水田土壤侵蚀量计算方法的调整
为保持新水田模块与原swat模型的参数一致性，仍采用改进通用流失方程作为水田土壤侵蚀量的计算公式；并根据水田的水位条件，对计算方法进行了优化调整。当水田处于浸没状态时，若未发生溢流或排水，则不计算土壤侵蚀过程。若有溢流或排水发生，则仅计算田埂侵蚀量，水田表土由于受到自由液面保护，不计算其侵蚀量。当水田处于排空状态下，水田田埂及田面均按可侵蚀面计算。
步骤六：水田氮、磷流失量的计算方法
氮、磷的流失量按其流失形态分别计算，氮的流失形态包括颗粒态氮和硝态氮，磷的流失形态包括颗粒态有机磷、颗粒态无机磷和溶解态磷。其中颗粒态氮、磷的流失量依照williamsandhann提出的公式(13)计算，溶解态氮磷流失过程的计算方法与原swat模型的方法一致。
sbp＝0.001con·(sed/sa)·ε(13)
式中：sbp表示吸附在颗粒物上的氮、磷含量，kg/ha；con表示表层10mm深度土壤中的氮磷含量，g/t；ε表示氮、磷在颗粒物上的吸附比率；sed表示土壤侵蚀量，t。
步骤七：改进模型的运行
录入模型所需的基础数据，根据模型率定结果对参数进行调整；设置与改进模型相关的参数：包括壶穴面积、田埂占比、水田逐日水位数据及运行模型。其中，模型所需的基础输入数据包括：数字高程数据dem、土地利用数据、土壤属性数据及气象数据。此外，壶穴面积的参数为pot_fr可在模型中直接调整，田埂占比和逐日水位数据需生成txt文件录入。
其中田埂占比以1％～3％为宜。
步骤八：计算水田非点源污染负荷
在模型输出文件中，读取水田的氮磷流失数据，分别计算流域中各水文响应单元的总氮及总磷的非点源污染负荷。
其中，步骤八中的水田总氮非点源污染负荷为地表流失的有机氮、地表流失的无机氮、土壤侧向流中的无机氮及由地下水进入河道中的无机氮的总负荷；水田的总磷非点源污染负荷为由地表流失的颗粒态无机磷、颗粒态有机磷及溶解态磷的总负荷。
3、优点及功效：
本发明所述的一种基于水位调控的灌区水田非点源污染负荷估算方法，以改进swat模型为基础，能够实现流域尺度的水田非点源污染负荷的估算，并有助于反映非点源污染的时空异质性特征，也有助于实现稻区的流域水环境的优化管理。
本发明所述的一种基于水位调控的灌区水田非点源污染负荷估算方法，充分考虑了水稻生长周期内的水位变化过程，有助于建立稻田水位控制与非点源污染之间的对应关系，同时也可以为水田的科学管理提供指导和帮助。
本发明所述的一种基于水位调控的灌区水田非点源污染负荷估算方法，分别考虑了水田在浸没状态下和排空期内的不同水文特征，同时也对各种算法进行了优化改良，有助于得到更为精确的稻田水文循环过程和养分流失过程。
【附图说明】
图1为改良水田模型运行流程的示意图。
图2为改良水田模型河道径流率定结果图。
图3为0cm～15cm原模型与改良模型的水田土壤水率定结果比较图。
图4为15cm～30cm原模型与改良模型的水田土壤水率定结果比较图。
图5为改良水田模型年均氮流失负荷模拟值与前人监测结果比较图。
图6为改良水田模型年均磷流失负荷模拟值与前人监测结果比较图。
图7为改良模型水田月均氮流失负荷模拟值的结果图。
图8为改良模型的水田月均磷流失负荷模拟值的结果图。
【具体实施方式】
见图1～图8，下面结合附图和东北三江平原的阿布胶灌区的实际案例，对本发明做进一步说明。
本发明一种基于水位调控的灌区水田非点源污染负荷的计算方法，该方法具体步骤如下：
步骤一：改进壶穴形状设定
基于图1中所示的改进模型运行框架，利用intelvisualfortran中的fortran77语言对swat模型的壶穴模块的源代码进行编译，建立能够符合水田的水文循环及污染输出特征的改进模型单元，即水田模块。对原模型中壶穴构型进行修正，由原swat模型中的圆锥体设定修正为截面积恒定的多边形结构，保证水田面积为恒定值，且不随水位变化而出现波动。此外,在原模型的基础上，增加了田埂面积的计算方法。其具体计算方法见公式(1)、公式(2)。
步骤二：改进水田蒸散发计算方法
在保证水田面积恒定的基础上，进一步修改壶穴模型源代码，对其蒸发量的计算方法进行调整。根据稻田的水位变化特征，分别编译排空状态和浸没状态下的两种水田蒸发量计算方法。先利用penman-monteith方程计算水田日潜在蒸散发量，此后根据各项蒸发量与水田日潜在蒸散发量的差值，对计算方法不断进行调整。当叶片截留量大于水田日潜在蒸散发量时，水田模块将只计算叶片截留蒸发量，如公式(3)；当叶片截留量小于水田日潜在蒸散发量时，其差值由水稻蒸散发量弥补，如公式(4)，其中水稻蒸散发量的具体计算方法如公式(7)、公式(8)所示。当水稻蒸散发量及叶片截留蒸发量之和仍小于水田日潜在蒸散发量时，水田模块将根据水位条件，分别进行计算：在浸没状态下，蒸发量的差额由田面水蒸发量补足，如公式(5)；在排空状态下，蒸发量的差额由土壤蒸散发量补足，如公式(6)，其中土壤蒸散发的计算方法与原swat模型一致。
步骤三：改进水田下渗量的计算方法
在步骤一和步骤二的基础上，进一步改进壶穴模块的代码，构建基于水位变化的水田下渗量计算方法。如公式(9)所示，当水田处于浸没状态时，若水田储水量大于水田日稳定下渗量，则水田实际下渗量等于水田日稳定下渗量；若水田储水量小于水田日稳定下渗量，水田储水量均视为水田实际下渗量，如公式(10)。当水田处于排空状态时，其水田实际下渗量的计算方法与原模型旱田下渗量的计算方式相一致。
步骤四：灌溉及排水控制单元的构建
基于水田的逐日水位控制需求，结合原swat模型的灌溉模块，构建三条水位控制线。最低水位为灌溉触发水位，低于此水位时，触发模型灌溉功能，如公式(11)；最适水位为水稻的最适宜生长水位，灌溉操作在触及此水位时停止；最高水位为排水触发水位，在降雨过程中，当水田的水位上升至此高度时，模型将触发排水功能，直至水位回落至最高水位线以下，如公式(12)。三条水位线均需设置逐日的水位值，并分别将对应水位输入至rice_paddy_depth.txt文件中。稻田水位在设置过程中应充分参考水稻的品种及各生长阶段的水位需求，本研究所模拟的水稻品种为三江1号，并根据泡田期、排水期、种植期、复青期、分蘖期、抽穗期、灌浆期、成熟期、晒田期、休耕期的水稻生产技术规程，对水稻生长周期内的水位进行了分别设定。
步骤五：水田土壤侵蚀量计算方法的调整
在完成步骤一至步骤四中水文循环模块的改良的基础上，对水田污染负荷的计算方法进行改良。为保持改进水田模块与原swat模型参数的一致性，仍采用改进通用流失方程musle对水田的土壤侵蚀量的进行计算。此处所用方法与原模型一致，只在计算方式上有所调整。当水田处于浸没状态时，田面由于受到液面的保护，不会发生土壤侵蚀，此时，仅考虑田埂部分的侵蚀计算；当水田处于排空状态时，即最高水位设置值为0，其侵蚀量为田面侵蚀量与田埂侵蚀量之和。其中田埂占比在rice_paddy_depth.txt文件中设定，并以水田面积的1％～3％为宜，本实施例中设置为1％。
步骤六：水田氮、磷流失量的计算方法
基于土壤侵蚀量的计算结果，通过设置氮、磷在土壤颗粒上的分配系数，计算水田颗粒态氮磷的流失量，如公式(13)；田面水中溶解态氮、磷的含量通过表层10mm土壤与田面水之间的交换量计算。完成上述所有编译工作后，利用编译器生成可执行的exe文件，并将其命名为swatuser.exe。
步骤七：改进水田模型的运行
运行swat模型，并输入dem数据、河网数据、土地利用数据、土壤属性数据及气象数据，并将土地利用类型为水田的水文响应单元的pot_fr参数值设置为1。完成所有模型基础设置后，将步骤五和步骤七中准备好的rice_paddy_depth.txt文件和swatuser.exe文件复制至swat工程文件下的textinout文件夹下。完成设置工作后，利用swatcup软件对模型参数进行率定，得到符合流域特征的基本参数。如图2所示，本案例中河道径流的率定和验证结果的决定系数和纳什系数均大于0.6，表明模型具有较好的适用性。在水田土壤水模拟结果方面，在0cm～15cm和15cm～30cm深度土壤范围内，模拟值与监测值的均方根误差分别为14.56mm和24.55mm(图3和图4)，显著优于原模型的模拟结果。在模型的可靠性得到验证之后，用swatcup所得的模型参数替换原模型中的参数，并再次运行模型，输出模拟结果。
步骤八：计算水田非点源污染负荷
根据输出结果，分别计算水田的氮、磷年均流失负荷，并结合其他研究中的水田非点源负荷的场地监测值，对模拟结果的合理性进行验证。如图5及图6所示，改进水田模型的氮、磷年均流失负荷的模拟结果均处于置信区间范围内。且从流域水田月均氮、磷流失负荷的模拟结果分析(图7和图8)，改进水田模型的模拟结果也符合案例研究区水田的氮、磷流失规律。其模拟结果显示，氮、磷流失的关键时期为排水和溢流频繁的7月和8月。因此，基于水位调控的灌区水田非点源污染负荷估算方法将有助于准确、快速的对流域水田的非点源污染进行计算。