海河流域农田氮磷面源污染的空间分布特征及关键源区识别

2021-04-28张巧玲胡海棠王道芸邱春霞李存军靖亭亭

灌溉排水学报 2021年4期

张巧玲，胡海棠，王道芸，邱春霞，李存军，白翠，靖亭亭

（1.西安科技大学测绘科学与技术学院，西安710054；2.北京农业信息技术研究中心，北京100097；3.国家农业信息化工程技术研究中心，北京100097）

0 引言

【研究意义】海河流域是我国水环境污染最严重的流域之一，其中农业面源污染是海河流域水质污染的重要来源[1]。2015年《中国环境状况公报》显示，海河流域为重度污染，地表水中V类及劣V类占60%。海河流域总面积为3.18×105km2，占全国总面积的3.3%，水资源仅占全国1.3%，却承载全国8%的耕地和10%的人口，贡献了全国粮食总产量的9.4%[2]。作为我国北方粮食主产区之一，海河流域化肥施用量一直居高不下，导致主要河流水质恶化，为该区域的面源污染治理带来了挑战[3]。鉴于此，研究海河流域农田氮磷面源污染的空间分布特征和识别关键源区，对海河流域农业面源污染防控具有重要意义。

【研究进展】面源污染负荷空间分布研究目前多采用风险评价或模型模拟的方法[4]。相比以氮磷指数法为代表的风险评价方法，模型模拟方法能够更为准确地进行面源污染负荷定量化评估。国内外研究主要采用SWAT 模型[5-6]、输出系数模型（ECM）[7-9]、AnnAGNPS[10]和平均浓度法[11]等模型方法。其中SWAT 模型、AnnAGNPS 等机理模型由于对数据要求高，计算复杂，不适合大区域模拟[1]。InVEST 是国际上应用最广泛的生态系统服务评估集成平台，其中的营养物传输率模型（NDR）所需参数比较简单，机理清晰，适宜用于大尺度面源污染模拟研究[12]。

关键源区识别是面源污染管理的重要组成部分[13]。研究表明，少数区域的输出污染物往往贡献整个流域污染负荷的大部分，对水体的环境质量起着决定性作用[14]，这些区域被称为关键源区。目前关键源区识别方法主要是基于面源污染负荷的空间分布异质性特征[15-16]。面源污染负荷产生和输移具有不同的空间分布特征[17-18]。李文超等[19]提出应依据入河负荷的空间分布识别关键源区。关键源区识别在很大程度上受污染源影响，但与水质等因素也是高度相关的[20]。一方面，农田入河负荷高的地方，水质污染不一定严重，后者受景观结构和水资源状况的共同影响。另一方面，水质污染严重的地区，不一定是因为农田污染源，需要将农田对水质的影响剥离出来。【切入点】目前农田面源污染负荷对河流水质的潜在影响和空间分布特征还缺乏研究。依据《中国环境状况公报》，海河流域水质状况存在显著的空间异质性，部分支流河段水质重度污染，而所属干流呈现轻度污染[21]。农田面源污染对水质的影响是否也具有相似的分布规律，还有待进一步研究。在防控资源有限的条件下，为了进一步提高农田面源污染优先防控区域选择的精准程度，需要研究农田面源污染负荷对河流水质影响程度的评估方法。结合面源污染负荷的空间分布特征，选择水质潜在影响严重的河段并追溯污染物贡献率大的源区进行关键源区识别。

【拟解决的关键问题】综上，本研究以海河流域为研究区，利用InVEST 营养物传输率NDR 模型和产水量模型，对海河流域农田氮磷入河负荷和产水量进行空间估算，并基于水文网络拓扑关系估算和分析河流断面氮磷入河通量和潜在氮磷径流质量浓度的空间分布特征，结合GIS 的空间热点分析方法和水文网络方法识别关键源区，以期为海河流域农田面源污染精准防控提供参考依据。

1 材料与方法

1.1 研究区概况

海河流域地处我国华北地区（112°—120°E、30°—43°N），包括北京、天津、河南北部、山西东部、山东北部、内蒙古和辽宁小部分和河北的绝大部分。流域总的地势是“西北高、东南低”，包括滦河及冀东沿海、北三河、永定河、大清河、子牙河、黑龙港运东河、漳卫南运河、徒骇马颊河八大水系[22]（图1）。整个流域的水系从南到北呈扇形分布，水系分散、支流较多。海河流域属温带大陆性季风气候，属半湿润半干旱地带，年均降水量为540 mm[23]。土地利用类型以耕地为主，耕地面积占流域面积的48.9%。海河流域北部及西部山区为一熟农区，中部和南部平原地区以冬小麦-玉米二熟种植为主，是中国三大粮食生产基地之一。

1.2 数据来源与预处理

本研究海河流域农田氮磷入河负荷和潜在氮磷径流质量浓度采用InVEST的NDR模型和产水量模型估算。NDR模型运行所需数据包括海河流域月降水量、土地利用、数字高程模型（DEM）、坡度、农区分布、坡度分布、子流域分布图、地表氮磷负荷等。产水量模型运行所需数据包括土地利用、月降水量、潜在蒸散量、土壤有效含水率等。具体数据如下：

1）气象数据。在中国气象数据中心下载海河流域相关气象站点2001—2015年降水量、气温、风速、日照时间等数据。利用python计算海河流域所有气象站点多年平均降水量，并利用克里格插值方法得到流域多年平均降水量空间分布图。由Penman-Monteith公式计算潜在蒸散量，利用反距离加权空间插值方法得到多年平均潜在蒸散发空间分布图。

图1 海河流域水系概化Fig.1 The river system of Haihe River Basin

2）DEM和2015年土地利用数据来源于中国科学院地理科学与资源研究所，空间分辨率均为30 m。基于DEM数据，利用ArcGIS水文分析模块提取研究区的子流域及水文网络空间分布。整个海河流域被划分为263个子流域，子流域出口设置河流断面，共273个断面，259个河段。

3）中国1∶100万土壤数据库。来自中国西部环境与生态科学数据中心。土壤有效含水率在SPAW软件中基于土壤质地及粒径组成数据，利用土壤有效含水率经验公式计算。

4）农区数据。海河流域农区分布数据来源于《中国农作制》[24]，在ArcGIS中矢量化得到海河流域农区的空间分布图。

5）统计数据。统计数据包括农作物播种面积、氮磷肥施用量、耕地面积等，来源于相关各省和县市统计年鉴。

按照《全国农田面源污染排放系数手册》，遵循“农区-坡度-种植模式”原则，根据农区分布图和各县农作物播种面积数据，将海河流域农田划分为13种面源污染模式（表1），并进行空间匹配（详见文章[25]），得到海河流域面源污染模式分布图。

表1 海河流域农田面源污染模式Table 1 Farmland non-point source pollution models in the Haihe River Basin

1.3 研究方法

1.3.1 农田氮磷入河负荷估算模型

农田面源污染氮磷入河负荷的估算分2步进行，首先基于农田面源污染模式及借鉴输出系数模型，估算农田氮磷流失负荷；然后基于InVEST 3.8.7的NDR模型，计算氮磷输移率和氮磷入河负荷。

1）氮磷流失负荷模型

借鉴输出系数模型，农田化肥施用氮磷流失负荷估算在栅格单元上进行，计算式为：

式中：Loadi,j为化肥施用氮磷流失负荷（kg）；i代表第i个县市，j代表第j类面源污染模式；Ej为栅格单元对应的面源污染模式j的肥料氮或磷流失系数（%），在《全国农田面源污染排放系数手册》中查找面源污染模式j对应的农田氮磷流失系数；Qi,j为栅格单元对应的县市i的面源污染模式j对应的氮磷肥施用折纯量（kg/hm2），根据调研各种植模式的氮磷肥常规用量和各县市农作物播种面积、氮磷肥总量等进行区域化修正。

2）氮磷入河负荷模型

NDR模型是基于质量守恒方法模拟氮、磷营养物在空间上的迁移过程，估算氮磷输移率和氮磷入河负荷。该模型可以模拟流域营养物的来源和输移过程，计算农田氮磷养分通过地表径流和壤中流迁移进入河流的传输率（即入河系数），从而估算海河流域农田氮磷入河负荷。本文对模型进行修正，将流失负荷模型和NDR模型结合起来。模型基本原理如下：

首先，对农田氮磷流失负荷按照一定比例进行地表和次表层分配。

式中：load_xi为农田氮磷流失负荷；prportionsubsurface为次表层营养物占比参数；loadsurf,i为地表营养物负荷；loadsubsurf,i为次表层营养物负荷。

再次，计算地表营养物传输率：

式中：NDR0,i被下游像素保留的营养物传输率，ICj为地形指数，IC0和k是校准参数。

计算次表层营养物传输率：

式中：NDRsubs,i为地下营养物传输率；effsubs是能够达到次表层径流的营养物最大截留效率；lsubs是次表层的截留长度，即假设土壤保持营养物的最大距离，l是从像素到流的距离。

氮磷入河负荷计算式为：

式中：xexp,i为各栅格单元i的入河负荷；loadsurf,i为地表营养物负荷；NDRsurf,i为地表营养物传输率；loadsub,i为次表层营养物负荷；NDRsub,i为次表层营养物传输率；xexptot为子流域营养物入河总负荷。

其中模型参数取值参考了InVEST用户手册[26]。NDR模型重要参数取值及依据如下：①次表层营养物来源占比，根据相关文献中海河流域地表和次表层径流分配状况，取值为0.5；②氮磷营养物最大截留效率由土地利用类型决定，本文林地、草地、耕地、水域、建设用地、未利用地的氮最大截留效率分别为0.8、0.75、0.5、0.05、0.05、0.01，磷最大截留效率分别为0.7、0.6、0.48、0.05、0.05、0.01，具体取值参考InVEST用户手册；③次表层的截留长度按照InVEST模型用户手册建议设定为1个栅格单元大小，即30 m；④被下游像素保留的营养物传输率和校准参数IC0具体计算参见InVEST用户手册，校准参数k默认值为2。

1.3.2 水文网络拓扑关系构建和氮磷入河通量

流域水文网络构建以河段为纽带，是流域汇流关系建立的基础。利用ArcGIS 的水文分析对流域的河段、径流节点和子流域等要素提取。参考蒋好忱[27]的方法，借助Matlab 编程构建整个流域的河网拓扑关系。滦河水系中上游区域水文网络拓扑关系的过程见图2。

基于氮磷入河负荷栅格数据，采用ArcGIS 的分区统计功能，得到每个子流域的氮磷入河总负荷，并基于水文网络拓扑关系，利用Matlab 查询每个河流断面处汇入的子流域清单，得到断面处的氮磷入河通量。

图2 水文网络（左）及子流域汇流拓扑关系（右）Fig.2 Hydrological network(left)and topological relationship between sub-basin（right）

1.3.3 潜在氮磷径流质量浓度估算

本研究通过构建潜在氮磷径流质量浓度指标来表征农田氮磷流失对河流断面处水质的影响程度。因缺乏海河流域各河流断面的径流量数据，故采用InVEST的产水量模型估算产水量，相当于潜在的最大径流量。潜在氮磷径流质量浓度定义为河流断面处氮磷入河通量与断面上游产水量的比值。

产水量模型基于水量平衡原理，是基于Budyko假设和年平均降水量的模型。将每个栅格单元的产水量定义为该栅格的降水量减去蒸散发量（包括地表蒸发和植物蒸腾）[28]。模型主要算法为：

式中：Yxj为年产水量（mm）；Px为年均降水量（mm）；AETxj为土地利用类型j上栅格单元x的年平均实际蒸散量（mm）。

式中：Rxj为Bydydo 干燥指数，无量纲；kxj为植被蒸散系数；ET0为潜在蒸散发，由Penman-Monteith 公式计算[29]。Z为季节因子，根据研究区的降水分布情况确定，降水分布均匀和夏季降水为主的地区，Z取值接近1，若冬季降水为主，则Z值接近10[30]。根据海河流域降水量分布情况Z取值为2。AWCx为栅格单元x的植被可利用水量，由根系深度、土层厚度和土壤有效含水率共同决定。

为了反映海河流域产水量多年平均水平，减少年际间气候差异，尤其是降水量差异对模型不确定性的影响，降水量和潜在蒸散量采用2001—2015年的多年均值数据。

1.3.4 氮磷关键源区识别

本文农田面源污染关键源区识别是在氮磷入河负荷和潜在氮磷径流质量浓度的热点分析基础上结合管理单元贡献率确定关键源区。

氮磷入河负荷的热点分析是利用ArcGIS 中热点分析（hot spots）工具，首先以子流域或县域等为单元进行分区统计得到总负荷或单位面积平均负荷，然后基于Getis-Ord Gi*指数进行冷点和热点类型分区[25]，Getis-Ord Gi*是一种识别高值或低值要素在空间上发生聚类的位置的空间统计方法，原理见式（8）。其中Gi*指数值为3 代表热点区、2 代表次热点区、1代表不显著、-2～0 代表次冷点区、-3 代表冷点区，本文将热点区视为关键源区。

式中：Xj为斑块j的值；为所有斑块的均值；Wij为斑块i和j的空间权重；n为斑块总数。

河段氮磷潜在径流质量浓度的热点分析，采用GIS 的Jenks 自然断点分级法划分不同潜在氮磷径流质量浓度等级。此外，采用贡献率方法定量评估各个子流域区对河流断面水质的贡献。贡献率方法是计算管理单元的污染贡献率，即针对某个河流断面，首先通过水文网络分析方法回溯其上游子流域或行政单元清单，计算断面上游污染物总负荷、各管理单元的总贡献率或单位面积贡献率，然后依据贡献率大小进行关键源区识别。

2 结果与分析

2.1 海河流域农田氮磷入河负荷空间分布格局

海河流域农田氮磷入河负荷见图3，农田氮、磷入河负荷均值分别为2.41 kg/hm2和0.56 kg/hm2，总负荷分别为3.493 4 万、0.807 7 万t。

在空间分布上，氮入河负荷整体呈现出“西北低，东南高”的特征，流域的中部和南部地区的氮磷入河负荷明显高于西北部山区；磷入河负荷空间分布无明显规律性，高值分布在海河流域南部和东北部的秦皇岛和唐山。西北部山区的氮、磷入河负荷均值分别为1.58、0.58 kg/hm2，中部和南部平原地区的氮、磷入河负荷均值分别为2.75、0.54 kg/hm2，即平原地区氮入河负荷显著高于山区，而磷入河负荷稍低于山区。山区种植模式以一年一熟为主，磷肥施用量远低于一年二熟为主的平原地区，但因山区土壤侵蚀强度高，磷流失系数和迁移系数较高，导致其磷输出仍然偏高。

图3 农田氮磷入河负荷Fig.3 The export load of nitrogen and phosphorus from farmland

图4 河段氮磷入河通量的空间分布Fig.4 Spatial distribution of export of nitrogen and phosphorus flux

以子流域为单元统计农田氮磷入河负荷，并基于水文网络拓扑关系，汇总每个河流断面的上游所有子流域入河总负荷，即得到各河流断面的氮磷入河通量（图4）。氮磷入河通量随河流方向自西向东累积，位于流域中东部的黑龙港运东河和子牙河下游河段的入河通量较高。

表2 八大水系氮磷入河负荷及通量Table 2 The export of nitrogen and phosphorus and flux of eight major water systems

海河流域八大水系的氮磷入河负荷和入河通量的空间统计结果详见表2。滦河水系和永定河水系农田氮入河负荷均值分别为1.24 kg/hm2和1.35 kg/hm2，其他六大水系均高于2.48 kg/hm2，其中徒骇马颊河水系高达2.92 kg/hm2。磷入河负荷高值区分布在漳卫河水系和徒骇马颊水系，分别为1.07 kg/hm2和1.43 kg/hm2。氮入河通量较高的水系包括子牙河水系、黑龙港运东河水系、永定河水系、漳卫河水系，均超过5 000 t；而磷入河通量较高的区域主要分布在永定河水系、黑龙港运东河水系、漳卫河水系、子牙河水系、徒骇马颊河水系，均超过1 000 t。

综上，徒骇马颊河水系、漳卫河水系、子牙河水系是海河流域八大水系中农田氮磷入河负荷较高的水系，漳卫河水系、子牙河水系、永定河水系、黑龙港运东河水系是农田氮磷污染物总量较大的水系。

2.2 海河流域潜在氮磷径流质量浓度空间分布格局

海河流域潜在氮径流质量浓度（图5）范围为0.07～23.96 mg/L，均值为5.97 mg/L，超过国家地表水V类限值（2 mg/L）3倍；潜在磷径流质量浓度范围为0.02～5.86 mg/L，均值为1.47 mg/L，超过国家地表水V类限值（0.4 mg/L）3.7倍，约55%的河段超过了国家地表水V类水质标准，其分布与2015年《中国环境状况公报》中海河流域水质分布状况较为一致[21]。这表明农田面源污染是海河流域重要的污染源，即使不考虑其他污染源的影响，农田氮磷流失仍然可能造成近1/2的河段处于严重污染状态。

由图5可知，海河流域潜在氮磷径流质量浓度呈沿山区向平原方向上升趋势，流域中部和南部的污染比较严重，高值区主要分布在徒骇马颊河全线、黑龙港运东河中下游、大清河支流和子牙河上游部分河段等。流域水系上游或支流的潜在氮磷径流质量浓度普遍高于干流。以国家地表水V类限值为标准，支流潜在氮磷径流质量浓度超标河段个数分别占支流总数的62%、71%，而干流超标数分别占干流总数的43%、47%。这主要是因为水系上游或支流所在集水区，都是耕地且占比高、化肥施用强度较高和耗水较多的农耕区。

图5 潜在氮磷径流质量浓度空间分布Fig.5 Spatial distribution of the potential runoff concentrations of nitrogen and phosphorus

2.3 海河流域农田面源污染关键源区识别

海河流域氮磷入河负荷的热点区（图6）主要分布在流域南部、石家庄和秦皇岛的周边地区。氮、磷热点区分别占流域总面积的14.45%、16.02%，贡献入河总负荷的23.39%、27.71%。

强调关键源区识别重要性的典型案例是Gburek等[31]得出流域20%区域贡献了80%的磷流失量。本研究关键源区的贡献率较低，主要原因是海河流域氮磷面源污染的空间异质性不显著，尤其是作为流域农田主体的平原地区，种植模式、施肥量、地形条件、耕地占比都很相似。因此，基于污染源负荷识别关键源区的方法在海河流域提高农田面源污染防控效率的作用不是很显著。

河段潜在氮磷径流质量浓度按照自然断点法划分为5个等级，从低到高分别为低值区、中低区、中等区、中高区和高值区。其中潜在氮、磷径流质量浓度高值区分别为44、35个河段（图6），分别占河段总数的17%、14%。

图6 氮磷入河负荷热点分析Fig.6 The hotspots analysis of export load of nitrogen and phosphorus

徒骇马颊河是海河流域农田面源污染最严重的水系，以其为例统计各县域对水系出口处水质的贡献率（图7）。该水系中上游各县的贡献率较高，在3%～6%之间，上游比中游略高，下游各县贡献率范围为1%～3%，说明徒骇马颊河水系上游是其主要污染源，需要重点治理。通过贡献率可以迅速判别对流域出口贡献程度最大的单元及贡献程度，也可作为流域生态补偿时的参考依据。

图7 徒骇马颊河各个县域的氮磷贡献率Fig.7 The contribution rate of nitrogen and phosphorus in each county of Tuhaimajia River

3 讨论

3.1 模拟结果的准确性及对比

由于本研究针对的是整个海河流域农田氮磷负荷的大尺度模拟，且NDR 模型估算的是总氮、总磷指标，而我国生态环境部数据中心可公开获取的水质数据不包括总氮、总磷指标，而且水体污染物除了农田面源污染物，还包括点源污染物和其他面源污染物，因此很难全面系统地评价模型模拟结果的准确度与可信性。

已有研究中，朱梅等[32]利用输出系数模型估算海河流域2007年种植业面源污染负荷，得出地表径流总氮、总磷分别为5.20 万、0.80 万t；冯爱萍等[33]利用DPeRS 模型估算的2016年海河流域农田径流型总氮和总磷排放量分别为12.7 万t 和0.7 万t。与已有研究结果对比发现，本文估算的农田氮、磷入河总负荷与上述结果在数量级上是一致的，具体估算值有所差别，这可能与不同的模型应用的数据源，研究时期以及模拟精度等有关。将潜在氮磷径流质量浓度的空间分布与2015年《中国环境状况公报》海河流域水环境分布图进行对比，较为一致，说明本研究采用的面源污染估算方法具有一定的区域适用性和可行性。

3.2 农田面源污染输移引起的空间差异性

潜在氮磷径流质量浓度（图5）与氮磷入河负荷（图3）呈显著的空间差异性。主要表现在：①西北部山区中氮磷入河负荷偏高的地区，其潜在氮磷径流质量浓度却较低，主要是因为这些地区耕地占比偏低，大量林地和草地对农田流失的氮磷污染物有拦截和稀释作用。②海河流域中部和南部平原地带的大清河、黑龙港运东河、徒骇马颊河的河段，潜在氮磷径流质量浓度较高，主要是因为这些地方种植模式以一年二熟为主，化肥施用强度大，且耗水严重导致径流量偏低，两者综合作用下加重了水质污染程度。潜在氮磷径流质量浓度是本地污染源和上游输入源污染物汇集与径流过程共同作用的结果。潜在氮磷径流质量浓度比氮磷入河负荷指标更能准确反映农田氮磷面源污染对水环境影响的程度。

3.3 关键源区识别

关键源区识别一般针对面源污染负荷进行[4,15-16,34]。本研究构建了潜在氮磷径流质量浓度指标，反映面源污染对水质的影响。潜在氮磷径流质量浓度高值区（图6 黑色圆圈内）与入河负荷热点区的空间分布格局不一致。因此，面源污染防控的关键源区识别应将入河负荷和潜在氮磷径流质量浓度的空间分布结合进行分析。对于潜在氮磷径流质量浓度高的河段，即使其集水区即源区位于入河负荷的非热点区，其治理优先序靠前。例如大清河水系东段支流的氮径流质量浓度很高，其源区以入河负荷次热点和不显著区为主，但面源污染防控时仍应重点治理。在选择优先防控区域时，应着重考虑对内源污染为主的潜在氮磷径流质量浓度高的区域进行源头控制，而对外源污染为主的潜在氮磷径流质量浓度高的区域可以采用过程拦截和末端治理的方式，并追溯其上游的关键源区。