多情景分析的农业面源污染关键源区识别软件开发及应用

覃苑，胡海棠，淮贺举，李存军*，张巧玲，杨铁利，王佳宇

1.辽宁科技大学

2.北京市农林科学院信息技术研究中心

农业面源污染（agricultural non-point source pollution，ANSP）指由于农业生产中的氮磷营养元素、农药及其他溶解态和固态污染物流失导致的污染。近年来，农业面源污染对水体污染的贡献越来越受到重视。2021年3月生态环境部和农业农村部联合发布了《农业面源污染治理与监督指导实施方案（试行）》，明确了“十四五”至2035年农业面源污染防治的总体要求、工作目标，提出主要任务包括深入推进农业面源污染防治、完善农业面源污染防治政策机制和加强农业面源污染治理监督管理。ANSP的污染负荷估算、关键源区源头识别和治理决策是监督管理的重要内容。

ANSP在空间上分布分散且隐蔽，其形成受到区域的气候、地形、土地利用类型和植被覆盖等多方面因素的影响，导致ANSP在监测、模拟和控制等方面面临较大挑战。GIS技术在空间数据管理、集成、分析方面有优势，已被广泛应用于ANSP的数据管理、空间分析、面源污染模型集成、信息系统开发和可视化等方面[1]。国内外采用年度农业面源污染模型（annualized agricultural nonpoint source pollution model，AnnAGNPS），水土评价模型（soil and water assessment tool，SWAT），输出系数模型（export coefficient model，ECM）等模型和平均浓度法等方法进行ANSP模拟与识别研究。如夏昊[2]结合AnnAGNPS模型，开发了基于GIS的面源信息系统，其拥有空间分析和可视化表达等功能，用于管理ANSP中不同来源的空间和属性数据；乔卫芳等[3]基于研究区域数字高程模型（digital elevation model，DEM）数据，将GIS与SWAT模型结合构建了ANSP基础信息库，分析了研究区ANSP负荷的分布特征，并模拟了区域内主要土地利用类型单位面积ANSP的负荷量。但AnnAGNPS和SWAT等模型所需参数多，计算复杂，在大区域的模拟上略显乏力；生态系统服务和权衡的综合评估（InVEST）模型中的营养物传输率模型（NDR）所需参数少，机理清晰，适用于大区域的ANSP模拟研究[4]。

在进行ANSP治理时，关键源区识别是面源污染管理的重要组成部分。研究表明，少数区域输出的污染物通常是整个流域污染负荷的主要贡献者，对水体环境质量有着至关重要的作用[5]。《2015年中国环境状况公报》显示，海河、巢湖等流域的水质状况存在明显的空间异质性，表现为支流水质在重度污染的情况下，其所属的干流仅为轻度污染[6]。目前主要是针对面源污染负荷空间分布的异质性特征进行ANSP关键源区识别[7]，但是单纯考虑污染负荷易忽略流域面积与产水量对污染负荷贡献的影响，导致关键源区的识别准确度下降，进而影响治理效果。应结合潜在氮、磷径流浓度和入河负荷的空间分布进行ANSP防控的关键源区识别与分析[8]。目前农业面源污染关键源区识别涉及水文分析、农业面源污染模型、栅格矢量运算和统计，对于面源污染非专业人士和基层工作者存在跨平台操作、数据类型多、资料汇总繁琐、计算量大等问题，亟需集成式的软件工具平台。

ANSP包括农田种植业面源污染、畜禽养殖、农村生活污染等，其中农田种植业面源污染是农业面源污染的重要组成部分，污染源相对分散且影响因素多。笔者针对在农田种植业氮、磷化肥施用形成的面源污染关键源区识别过程中情景单一、数据量大、涉及环节多，存在计算复杂、操作繁琐等问题，将GIS与InVEST模型结合，研发农业面源污染关键源区识别与治理模拟的集成计算机软件系统，在多种情景下对关键源区进行识别，并模拟治理后的效果，以期为农田种植业面源污染治理提供有效的辅助手段。

1 ANSP关键源区识别软件系统设计

基于GIS二次开发，集成水文网络拓扑构建和面源污染关键源区识别技术，利用DEM数据提取河流节点、河段与子流域等水文地貌要素，整合InVEST产水量模型和营养物传输率模型，结合基于入河污染负荷、潜在径流污染物浓度和负荷与产水量比值3种情景的面源污染关键源区识别方法，以决策支持为重点，建立农业面污染关键源区识别与治理模拟软件系统。

1.1 系统目标

通过研究ANSP的污染负荷情况与成因，运用GIS技术，研制具有农业面源污染关键源区识别与治理模拟的辅助决策软件系统。对研究区域内的水文网络进行拓扑构建、分析和专题信息可视化表达，快速高效得到研究区内拓扑关系、子流域营养元素分布规律和子流域潜在营养元素径流浓度等信息。在此基础上，根据不同的方法分析面源污染关键源区空间位置，并模拟关键源区治理后的效果，为ANSP治理提供依据。

1.2 功能设计

以灵活的应用集成接口构建GIS软件系统，对研究区内的水网数据、DEM数据和土壤数据等多源异构的资源环境数据进行集成与调度，根据研究区已知的地理数据和资源环境数据，识别面源污染关键源区。系统按照实际功能分为6个模块单元：1)拓扑关系分析单元，用于基于目标区域的DEM数据构建水文网络拓扑关系；2)产水量分析单元，用于利用InVEST模型的产水量子模型，根据水文网络拓扑关系，确定子流域产水量；3)入河负荷估算单元，根据地表营养物负荷、传输率、次表层营养物传输率和各栅格单元的入河负荷，计算子流域营养物入河总负荷；4)潜在径流浓度分析单元，用于根据子流域产水量，确定目标区域内各子流域的潜在氮、磷径流浓度；5)关键源区识别单元，基于氮、磷入河负荷，潜在氮、磷径流浓度，子流域氮、磷负荷与产水量比值3种情景，识别和确定目标区域内的面源污染关键源区；6)关键源区治理模拟单元，用于模拟对识别出的关键源区进行治理后的结果，确定治理效果。

ANSP关键源区识别与治理模拟统一建模语言如图1所示。由图1可知，1种关键源区识别结果对应1个关键源区识别模型，模拟治理结果依赖于关键源区识别结果；1个关键源区识别模型对应若干农田氮、磷入河负荷估算模型，子流域负荷与产水量比值估算模型与潜在氮、磷径流浓度估算模型；农田氮、磷入河负荷估算模型包含1个氮、磷流失负荷分配结果，1组地表营养物传输率和1组次地表营养物传输率；子流域负荷与产水量比值估算模型包含1个产水量模型，依赖于农田氮、磷入河负荷估算模型结果；潜在氮、磷径流浓度估算模型包含1个产水量模型，依赖于水网拓扑关系和农田氮、磷入河负荷估算模型结果。ANSP关键源区识别与治理模拟软件结构如图2所示。

图1 ANSP关键源区识别与治理模拟统一建模语言Fig.1 Unified modelling language for the identification and control simulation of ANSP key source areas

图2 ANSP关键源区识别与治理模拟软件结构Fig.2 Software structure diagram for identification and simulation of the treatment of ANSP key source areas

2 ANSP关键源区识别主要方法与模型

2.1 水文网络拓扑构建

获取研究区的DEM数据，基于DEM数据建立水文网络拓扑关系。利用计算机编程提取河段以及子流域，径流节点是基于几何网络的方法提取的，分为源头节点、汇流节点、流域出口节点[9]。具体流程见图3。

图3 水文地貌要素提取Fig.3 Extraction of hydrological and geomorphic elements

流域水文网络拓扑结构的建立以流域汇流关系建立的基础——河段为纽带，对河网的径流节点、河段、子流域进行矢量要素的空间分析，构建拓扑关系。河段属性表里记录对应河段的起始节点和终止节点，以及对应的子流域编码；子流域属性表里记录与其相对应的河段编码。定义节点的编码为Point-ID，节点类型是NodeType，上下游节点分别为FROM_NODE、TO_NODE；河段编码为 River_ID，长度为River_len，汇入和汇出河段分别为Riv_in、Riv_out；子流域编码为Basin_ID，每个子流域的汇入和汇出流域编码分别为Basin_in、Basin_out。

2.2 流域氮、磷入河负荷计算方法

流域氮、磷入河负荷的计算分为2步，首先计算出流域氮、磷流失负荷，再将其与InVEST模型中的NDR模型进行结合，计算得到流域氮、磷入河负荷。

2.2.1 氮、磷流失负荷模型

流域氮、磷流失负荷计算公式如下：

2.2.2 氮、磷入河负荷模型

通过模拟区域内营养物的来源和转移过程，计算氮、磷入河系数，在考虑氮、磷流失负荷的情况下估算氮、磷入河负荷。首先按照式（2）将氮、磷流失负荷分配到表层和次表层：

最后，汇总计算氮、磷入河负荷，计算公式如下：

2.3 流域产水量计算方法

利用InVEST模型的产水量子模型计算流域产水量，计算公式如下：

为了反映流域多年平均产水量，降低年间的气候差异和降水量差异造成不确定性的影响，潜在的蒸发散失量与降水量采用2001—2015年的平均值。

2.4 ANSP关键源区的识别

由于面源污染易受到地形、植被、气候、施肥和土地利用类型等因素影响[12-15]，不同识别方法通常会呈现出不同程度的空间差异性。本研究采用3种不同情景识别关键源区。

2.4.1 基于入河负荷

该情景由分区统计获得子流域氮、磷入河负荷，根据每个子流域氮、磷入河负荷大小进行关键源区的识别[16-19]，计算公式如下：

2.4.2 基于潜在径流浓度

该情景通过计算河段氮、磷入河通量与径流量的比值，得到河段潜在径流浓度，根据潜在径流浓度的贡献率大小进行关键源区的识别[8]。

2.4.3 基于负荷与产水量比值

将子流域产水量与氮、磷入河负荷同时纳入考虑，计算子流域产水量与氮、磷入河负荷的比值，根据该比值进行关键源区的识别[20-22]，计算公式如下：

2.4.4 关键源区识别

根据式（8）～式（11）计算得到入河负荷、潜在径流浓度和负荷与产水量比值，对子流域或河段对应子流域进行降序排列。按顺序累加子流域入河负荷，每累加1次就通过子流域累计入河负荷与流域总负荷之比（Per）来判断是否满足设定阈值，若大于等于设定阈值则停止累加并将这些区域标记为面源污染关键源区，流程如图4所示。

图4 关键源区识别流程Fig.4 Identification process for key source area

2.5 关键源区治理模拟

针对3种情境识别出的关键源区，设定相同的预期治理幅度，可直观地观察到在不同情景下对关键源区治理的效果。将3种情景下识别出的关键源区，按照与河流源头的接近程度进行排序，根据设定治理幅度依次降低这些区域的氮或磷负荷，再从上游—下游的顺序依次更新流域所有受影响的河段的氮或磷潜在径流浓度。

3 农业面源污染关键源区识别软件应用

3.1 研究区概况及特征

海河流域（112°E～120°E，30°N～43°N）地处我国华北地区，包括北京、天津、河南北部、山东北部、山西东部、内蒙古和辽宁小部分区域与河北省大部分区域。流域地势为西北高、东南低，包含滦河、北三河、大清河、永定河、子牙河、漳卫南运河、黑龙港运东河、徒骇马颊河与冀东沿海八大水系[23]。海河流域是我国三大粮食生产基地之一，流域土地利用类型主要为耕地，占流域总面积的48.9%。ANSP是海河流域主要污染源之一，也是海河水质污染的主要原因之一[24]。

3.2 总氮关键源区识别

研发的农业面源污染关键源区识别软件主界面见图5。以TN为例，使用该软件对海河流域的农业面源污染关键源区进行分析。流域ANSP的面积与TN、TP的关系通常是呈对数上升趋势，往往不到10%的面积就贡献了30%以上的氮、磷污染负荷[25]，据此将软件识别幅度设定为30%，分别基于入河污染负荷、潜在径流污染物浓度、负荷与产水量比值3种情景进行TN面污染关键源区识别，结果如图6、图7所示。

图5 农业面源污染关键源区识别软件用户界面Fig.5 User interface of the key source area identification software for ANSP

由图6、图7可知，TN入河负荷较高的区域主要分布在滦河、北三河、璋卫河水系与黑龙港运东河下游，子牙河、大清河水系中部以及徒骇马颊河全域，呈现自西北向东南方向增长的分布特征。潜在TN径流浓度升高的区域主要分布中部—东南部地区的黑龙港运东河、徒骇马颊河与大清河水系。TN负荷与产水量比值高的区域集中分布在中部—南部地区的徒骇马颊河上游、大清河水系及黑龙港运东河。基于入河负荷情景下识别出的关键源区较为分散，除了永定河与子牙河水系外皆有分布。基于潜在径流浓度情景下的关键源区分布较为集中，多为中部—南部分布，其余在东南部也有分布。基于流域负荷产水比值情景下识别出的关键源区分布非常集中，以中部—南部分布为主，识别出的区域与基于潜在径流浓度情景下的区域重合度较高。

图6 3种情景下的TN污染区域分布情况Fig.6 Regional distribution of total nitrogen pollution under three scenarios

图7 3种情景下识别的关键源区Fig.7 Key source areas identified by three scenarios

3.3 治理模拟结果与分析

在识别出关键源区的基础上，对关键源区进行模拟治理。为了更好地反映治理效果，将识别出的关键源区TN、TP入河负荷清零，并在考虑水网的拓扑关系下，对河段TN、TP潜在径流浓度进行重新统计。模拟结果如图8所示。由图8可知，基于入河负荷情景下，模拟治理的效果欠佳，流域中部水系的潜在径流浓度下降不明显；基于潜在径流浓度、负荷与产水量比值2种情景下，治理效果相近，治理区域的潜在径流浓度均有较为明显的下降，但对东南部和南部的部分水系治理效果不明显。

图8 原始径流浓度与3种情景下治理效果对比Fig.8 Comparison chart of original runoff concentration and that of three methods after treatment

计算了3种情景识别出的面污染关键源区模拟治理后TN、TP负荷的下降率，结果见表1。

由表1可知，在研究区的数据与条件下，基于入河负荷情景识别出的面污染关键源区的TN治理效果弱于基于潜在径流浓度和基于负荷产水量比值情景的治理效果。由于海河流域西北地区入河负荷高，但林地和草地面积较大，拦截和稀释了大量的氮、磷污染物，而中南部地区耕地面积占比高，且采取了一年两熟的种植模式，化肥的施用量大且耗水严重，导致地表径流量偏低，所以3种情景下识别出的关键源区有着较显著的空间差异性。由于基于潜在径流浓度情景反映的是河段与其上游输入污染物的汇集和径流过程共同作用的结果，而子流域负荷与产水量比值反映子流域产水量与氮、磷污染负荷比值，二者更能反映出ANSP对水质的影响程度。

表1 3种情景下关键源区模拟治理效果Table 1 Simulation of governance effects in key source areas under three scenarios

对海河流域的治理模拟结果表明，仅针对面源污染入河负荷识别出的氮和磷污染关键源区在治理的结果上并不理想，基于潜在径流浓度和基于负荷与产水量比值情景下的治理效果优于前者，且结果接近。这是因为面源污染受到地形、植被和土地利用类型等因素影响较大，单一情景下识别效果并不一定能准确反映面源污染对水质的影响程度。考虑到不同地区地理环境的复杂性，单一识别方法具有一定的局限性，在实际使用中，ANSP防控关键源区的识别应将入河负荷、潜在径流浓度和负荷产水量比值的空间分布情况结合分析，以期得到更好的治理效果，支持农业面源污染治理科学决策。

4 结论

(1)针对农业面源污染关键源区识别过程中情景单一、数据量大、涉及环节多，以及模型复杂、操作繁琐和治理效果不明确等问题，设计开发了多情景分析的农业面源污染关键源区识别软件。软件具备拓扑关系分析、产水量分析、入河负荷估算、径流浓度分析、关键源区识别和关键源区模拟治理等功能，可模拟3种常见情景下农业面源污染分布的特点，为农业面源污染治理决策提供有效的辅助手段。

(2)以海河流域为例，进行了农业面源污染关键源区的识别和治理模拟应用，关键源区识别结果显示，流域TN负荷呈现自西北向东南方向增长的分布特征，中部—东南部地区的潜在TN径流浓度较高，中部—南部地区的TN负荷与产水量比值明显高于其他地区。模拟治理的结果表明，基于入河负荷情景对氮或磷的面源污染关键源区治理效果较不理想，而基于潜在径流浓度、负荷与产水量比值情景的治理效果较好，且后者略优于前者。