基于磷指数模型的海河流域农田磷流失环境风险评价

吴汉卿，万 炜，单艳军，陈延华，李 强，李存军，胡海棠，张宝贵

基于磷指数模型的海河流域农田磷流失环境风险评价

吴汉卿1，万 炜1，单艳军1，陈延华2，李 强1，李存军3，胡海棠3，张宝贵1※

（1. 中国农业大学土地科学与技术学院，北京 100193；2. 北京市农林科学院植物营养与资源研究所，北京 100097；3. 北京农业信息技术研究中心，北京 100097）

农田面源磷流失是导致水体富营养化的主要原因，识别农田磷流失的关键源区、影响因子是农田面源污染防治的重要环节。该研究以海河流域为研究区，采用磷指数模型，选取土壤有效磷含量、磷肥施用量作为源因子，以土壤侵蚀模数、年径流深、农田和水体间归一化距离指数作为迁移因子，结合GIS技术评估海河流域农田磷流失风险，并利用结构方程模型研究农田磷流失风险指数与各影响因子间关系。结果表明：1）海河流域农田土壤有效磷、磷肥施用量、土壤侵蚀模数、年径流深及归一化距离指数处于中-低、中-高、极低、中-高和高级别风险等级的区域面积占比最高，分别占农田总面积的66.5%、61.1%、99.0%、54.2%和64.8%；2）影响农田磷流失的关键因子为迁移因子，其中关键的迁移因子为年径流深及归一化距离指数；3）源因子与迁移因子间呈极显著负相关（<0.01），土壤性质（包括土壤质地及有机碳含量等）与源因子呈极显著负相关（<0.01），与迁移因子呈极显著正相关（<0.01）；4）海河流域农田磷流失关键源区位于黄河北岸的山东省和河北省东南部的平原农耕区、海河流域西北部的山区地带。该研究结果对流域尺度上农田磷流失研究的方法创新有所裨益，可为海河流域农田面源污染防治提供科学参考。

农田；模型；面源污染；磷流失；磷指数模型；关键源区；结构方程模型；海河流域

0 引 言

磷（P）是作物生长的主要限制性养分，农业生产过程中施用大量磷肥[1]。农田土壤磷素在降水冲刷及土壤侵蚀作用下，以径流、壤中流、吸附于泥沙表面等方式进入受纳水体[2]。据估计，全球每年有3×106~4×106t P2O5从土壤迁移到水体中[3]，其中约19.5 kg/hm2来自中国农田[4]。农业生产引发的磷流失是水体富营养化的最直接原因，已成农田面源污染防控研究的热点问题[5]。

农田面源污染因受土壤类型、地形特征、气候条件和人类活动等诸多因素影响具有高时空异质性，农田面源污染防控难度较大[3,6]。在流域尺度上，识别出高污染风险的关键源区（Critical Source Areas，CSAs）是农田面源污染防控的前提[7-8]。磷指数模型（Phosphorus Index Model，PIM）是识别关键源区及评估污染风险的经典模型之一[9-10]。该模型作为磷素综合管理工具，已被国内外学者广泛应用[10-12]。磷指数模型最初根据源因子与迁移因子特征用以评估田间尺度上的磷污染风险等级，后来Gburek和Sharpley等将农田到河网的水文连通性、距离等因子也考虑在内，对流域尺度上的磷流失潜在风险指数进行了定性化评估[9-10]。Buchanan等开发了基于径流行程时间的磷指数模型，以改进磷指数表征排水对磷流失潜在风险的影响[13]。在中国，张淑荣等最先采取磷指数模型，选取土壤侵蚀模数、地表径流、潜在污染源距水体的距离、土壤磷含量、磷肥施用量和施用方法等因子，识别出了小流域磷流失高风险区[14]。其后，学者们通过将土壤磷吸附指数（Phosphorus Sorption Index，PSI）、磷饱和度（Degree of Phosphorus Saturation，DPS）、溶解态磷、颗粒态磷及地形指数因子等指标纳入评价体系，对磷指数模型进行修正，对流域尺度上磷流失潜在风险进行评价并识别出关键源区[15-16]。综上，修正的磷指数模型通过定性评估流域尺度上农田磷流失潜在风险指数，从而识别关键源区，但对如何识别关键影响因子尚不清楚。

2011年“重点流域‘十二五’农业面源污染防治研究报告”指出，在8个重点流域中，海河流域的总磷总负荷量及排放强度分别位居第一和第二位，总磷总负荷量占8个重点流域总量的43.1%，排放强度是8个重点流域平均值的2.1倍。2018年《中国生态环境状况公报》指出海河流域为中度污染，地表水中劣Ⅴ类占20.0%。因此，海河流域农田磷流失环境风险评价对于磷污染防治具有重要意义。尽管磷指数模型在小流域尺度上农田磷流失关键源区识别方面得到了广泛应用和发展，而海河流域空间范围广，空间分异规律明显，农田面源污染的驱动机制难以较好地把握，因此在海河流域农田面源污染防治中的应用极少，关于农田磷流失关键源区、关键影响因子识别及磷流失潜在风险指数与各影响因子间关系研究亦是鲜见报道。

因此，本研究以海河流域为研究对象，以1 km栅格尺度为最小地域单元（传统方法主要以行政区为最小地域单元进行空间分析），利用磷指数模型并结合GIS技术，通过分析影响海河流域农田磷流失的影响因子及潜在风险指数的空间分布，识别出海河流域农田磷流失的关键源区，在此基础上采用结构方程模型方法定量地探讨海河流域农田磷流失潜在风险与各影响因子间关系，进而识别出关键影响因子。本研究结果希望对流域尺度上农田磷流失研究方法创新有所裨益，并为海河流域农田面源污染防治提供科学参考。

1 研究区概况及数据源

1.1 研究区概况

海河流域地处华北（112°～120°E、35°～43°N），共辖5省（河北省、山西省、山东省、河南省、辽宁省）、2直辖市（北京市、天津市）及1自治区（内蒙古自治区），流域总面积达3.18×105km2。该区域整体地势为西北高、东南低，西北部为燕山-太行山山脉，东南部则为黄河及海河形成的冲积平原（图1a）。境内主要河流为海河与滦河。气候类型以暖温带大陆性季风气候为主，年均温为1.5～14 ℃，多年平均降水量为539 mm。土壤类型主要有华北山地棕壤与褐土，海河平原潮土、黄垆土、盐碱土及风沙土，内蒙古高原栗钙土与绵土。

土地利用类型以耕地为主，其中耕地面积占区域总面积的48.9%（图1b）。该区典型的农作物种植方式为冬小麦-夏玉米轮作。该区域是中国三大粮食生产基地之一和京、津等大中城市“菜篮子”主要生产基地，人口密度较高，人为作用强度较大，农田面源污染严重，引发水体富营养化等环境问题[17]。

1.2 数据源及预处理

1）县域主要作物的粮食种植面积及磷肥施用量数据：从研究区各市统计年鉴（2019年）提取出化肥折纯量，再根据2019年中国统计年鉴中氮肥、磷肥及钾肥施用量分别为2065.43、728.88和590.28万t，参考全国平均比例折算出研究区磷肥施用量数据。

2）土地利用、DEM及土壤质地数据：由中国科学院资源环境数据云平台（http://www.resdc.cn/）下载得到。土地利用数据为2015年的栅格数据产品。研究区数字高程模型（DEM）利用ArcGIS水文分析模块提取研究区河网空间分布[18]，并用于提取出研究区内的坡度及坡长因子。土壤质地数据用于提取研究区土壤的砂粒、粉粒和黏粒含量。上述数据的空间分辨率均为1 km。

3）土壤全磷含量和土壤有机质含量数据：空间分辨率为1 km，从国家地球系统科学数据中心：土壤分中心（http://soil.geodata.cn/）下载得到。本研究在前人研究基础上[19-20]，通过将有机质含量除以1.724[19]和土壤全磷含量乘以转化系数（6%）[20]的方式分别得到土壤有机碳及有效磷含量。

4）1980—2015年逐月降水量数据：由中国气象数据网（http://data.cma.cn/）下载得到，先通过ANUSPLIN气象插值软件[18]经薄板样条函数处理得到研究区1 km空间分辨率的研究区逐月降水量数据，并处理得到多年平均逐月降水量。

5）2019年逐月的MODIS-NDVI数据：通过美国国家航空航天局（National Aeronautics and Space Administration，NASA，https://ladsweb.modaps.eosdis.nasa.gov/search/）下载得到。通过MRT软件对其进行镶嵌、投影转化、重采样及裁剪等处理，再通过最大值合成法[21]得到研究区NDVI数据，用以计算作物生育期的植被覆盖度。

2 研究方法

2.1 磷指数模型

磷指数模型是一种综合考虑农田磷流失的主要影响因子，进而评价流域内不同区域发生磷流失潜在风险高低的经典模型[9-11,22]。本研究基于磷指数模型来评估海河流域农田磷流失风险，参照国内外相关研究成果[9-11,22]，将影响磷流失的因子划分为源因子与迁移因子。在本研究中，选取了土壤有效磷含量及磷肥施用量作为源因子，这是因土壤有效磷能反映土壤本身磷含量及磷流失潜在能力，而磷肥施用量作为主要的外源磷投入，对土壤磷素含量具有重要影响[11]；选取土壤侵蚀模数、年径流深及农田与河网间归一化距离指数作为迁移因子，这是由于此三者为农田磷流失至受纳水体过程中最主要的影响因素，分别反映了土壤侵蚀、地表径流及农田至受纳水体的距离[11-12]。

采用ArcGIS中的自然断点法对源因子、迁移因子进行等级划分[12]，分为5级（极低、低、中、高、极高），并对各等级进行风险赋值，因子风险值及权重值确定，均参照前人研究结果进行确定[9-11]，具体数值见表1，进而采用磷指数模型计算出海河流域农田磷流失风险指数来反映研究区农田磷流失风险等级情况，再根据风险等级识别出海河流域农田磷流失关键源区[9-11]。

表1 农田磷流失风险等级划分

注：RV代表风险值。

Note: RV represents Risk Value.

2.2 源因子计算

2.2.1 土壤有效磷含量

表层土壤有效磷含量作为土壤磷流失的物质基础，含量越高则磷流失风险越大[14]。本研究利用ArcGIS生成海河流域1 km栅格土壤有效磷含量空间分布。

2.2.2 磷肥施用量

2015年农业农村部339个国家级基层肥料零售网点农户调查资料显示，海河流域磷肥施用类型主要以化肥为主，施肥方式以基肥形式为主。因此，本研究利用ArcGIS将统计年鉴中的县域磷肥施用量数据按照土地利用类型对研究区内耕地进行空间插值，精确至1 km空间分辨率的耕地栅格中，从而生成海河流域磷肥施用量空间分布图。

2.3 迁移因子计算

2.3.1 土壤侵蚀模数

本研究选用修正通用土壤流失方程（Revised Universal Soil Loss Equation, RUSLE）来估算海河流域土壤侵蚀模数，作为农田土壤磷流失风险的重要迁移因子之一[9-11]。利用如下公式进行计算，并获取海河流域土壤侵蚀模数（）分布

式中为土壤侵蚀模数，t/(km2·a)；为降雨侵蚀力因子，MJ·mm/(hm2·h·a)；为土壤可蚀性因子，t·hm2·h/(MJ·mm·hm2)；为地形因子，无量纲；为植被覆盖与管理因子，无量纲；P为水土保持措施因子，无量纲。

1）降雨侵蚀力因子

降雨侵蚀力因子（）是用来评估降雨潜在引起的土壤侵蚀能力的动力学因子。本研究基于海河流域1980–2015年的平均逐月降雨量数据，采用Wischmeier的公式计算[23]，并获取海河流域分布

式中为多年平均降雨量，mm；p为第月多年平均降雨量，mm。

2）土壤可蚀性因子

土壤可蚀性因子（）定义为单位降雨侵蚀力造成的单位面积侵蚀土壤数量，反映了不同土壤的侵蚀速率，是影响土壤侵蚀的内在因子，与土壤物理性质如质地和有机质含量等紧密有关。本研究利用研究区1 km栅格土壤砂粒、粉粒、黏粒及有机碳数据，采用Williams等建立的EPIC（Erosion Productivity Impact Calculator）模型计算[24]，并获得海河流域分布

式中SAN、SIL和CLA分别为0~20 cm土层砂粒、粉粒和黏粒含量，%；和SOM分别为土壤有机碳和有机质含量，%；S为粉粒和黏粒含量总和，%。

3）地形因子

地形因子（）亦称坡长坡度因子，主要反映坡长因子（）与坡度因子（）对土壤侵蚀的影响。定义为其他条件一致情况下，特定坡长的坡地土壤流失量与标准小区坡长的坡地土壤流失量比值；定义为其他条件一致情况下，特定坡度的坡地土壤流失量与标准小区坡度的坡地土壤流失量比值。本研究利用ArcGIS水文分析功能从数字高程模型中提取出坡长与坡度，再根据下列公式依次计算出、与[25]，并获得海河流域地形因子（）分布

式中和分别代表坡长因子和坡度因子，无量纲；是由DEM提取出的水平投影坡长，m；22.13是RUSLE采用的标准小区坡长，m；是由DEM提取出的坡度，(°)；是坡度坡长指数，无量纲。

4）植被覆盖与管理因子

植被覆盖与管理因子（）是指在一定条件下，有作物或植物覆盖及实施管理情况下坡地土壤流失量与耕种过后连续休闲地土壤流失量之比。本研究先利用ENVI软件获取归一化植被指数（NDVI），再计算植被覆盖度（），最后利用估算[26]，获得海河流域分布

式中为植被覆盖度，%；NDVI为归一化植被指数；NDVImax和NDVImin分别为研究区内NDVI的最大值和最小值。

5）水土保持措施因子

水土保持措施因子（P）是指采取特定水土保持措施下的土壤流失量与相应未采取措施顺坡耕作下的土壤流失量之比。本研究采用被广泛应用于大尺度流域P值计算的经验模型[27-28]，获得海河流域P分布

式中是由DEM提取出的坡度，%。

2.3.2 年径流深

年径流深（R）定义为每年将径流量均匀地铺在整个流域面积上所相当的水层深度，是一种度量地表径流的经典方法。本研究利用海河流域近35 a（1980–2015年）的多年平均逐月降雨量数据计算出多年平均降雨量，以下列公式计算R，并获得海河流域R分布

式中R为年径流深，mm/a；为年径流系数。参考《海河流域水利手册》[29]和刘春蓁等研究结果[30]，本研究中取值为0.12。

2.3.3 归一化距离指数

农田与受纳水体的距离是影响磷流失的重要因子，距离越远，在迁移过程中被稀释及截留的可能性越大，即潜在的磷流失风险愈低，反映了农田与河网间的距离及潜在的磷流失风险[31]。本研究在前人[32]的研究基础上，将研究区内农田与河网间距离进行归一化处理，首次提出了归一化距离指数（Normalized Differential Distance Index，NDDI）这一概念，得到农田与河网间距离指数。NDDI指数与传统的主观设定农田与水体间距离相比，具有客观定量评价的优势，使研究结果更加精准与可信。

式中为各像元中农田到河网间距离，max和min分别为农田与河网间距离的最大值和最小值。

2.4 磷流失风险指数计算

本研究在前人提出的磷指数模型评价指标体系基础上[11-12]，根据研究区源因子及迁移因子的分级标准与权重，计算出磷流失风险指数。

式中PI表示磷流失风险指数；S和W分别为源因子评价指标对应的风险值及权重，T和W分别为迁移因子评价指标对应的风险值及权重。PI的计算通过对源因子指标与迁移因子指标相乘，以此反映迁移因子对源因子的限制作用。

2.5 统计分析

结构方程模型（Structural Equation Modeling，SEM）是一种验证性的多元统计方法，能反映潜变量与外显变量之间的关系，同时探讨不同因子对同一过程直接或间接的影响，判别不同因子的相对重要性。因此，本研究利用Amos 21.0构建SEM，分析土壤性质（包括土壤质地及有机碳含量等）与磷流失风险指数、源因子（包括土壤有效磷含量、磷肥施用量等）及迁移因子（包括土壤侵蚀模数、年径流深及耕地与水体间距离因子等）间的关系。在构建SEM之前，采用SPSS 21.0对PI及各因子进行相关分析。使用卡方值/自由度（CHI/DF）、配适度指数（Goodness-of-Fit Index，GFI）、调整后配适度指数（Adjusted Goodness-of-Fit Index，AGFI）、赤池信息量准则（Akaike Information Criteria，AIC）和渐进误差均方根（Root Mean Square Error of Approximation，RMSEA）对SEM的拟合度进行评价。

3 结果与分析

3.1 农田磷流失源因子空间分布

海河流域内耕地土壤有效磷质量分数范围为0～193 mg/kg，均值为45 mg/kg。高值区主要分布在燕山-太行山山脉地区及渤海湾滩涂地区（图2a），占农田总面积的33.5%，前者可能是由于山区土壤有机质含量高，在分解过程中会产生有机酸、腐殖酸等物质，减少土壤对磷的吸附，导致土壤中磷含量高[33]，后者可能是由于渤海湾滩涂地区沉积物中磷含量较高，且黏粒、粉粒含量较高，对磷吸附能力较强[34]；中值区主要分布在海河及黄河冲积平原，且与耕地分布基本吻合，这可能是人为施肥及土壤磷素累积所致；低值区主要分布在燕山-太行山山麓地区，这可能是由于地形起伏，坡度较大，土壤侵蚀严重，导致含量低。整体上，研究区内土壤有效磷含量中等的区域占农田总面积为66.5%。

海河流域内磷肥施用量范围为0～115 kg/hm2，均值为52 kg/hm2。空间分布与土壤有效磷空间分布相反（图 2b）。高值区主要分布在海河与黄河冲击平原区域的农耕区，其中山东省与河北省东南部磷肥施用量显著高于其它省份，这是因为黄河及海河的冲积平原地势平缓，适宜耕作，因此耕作强度大，化肥施用量高；低值区主要分布在燕山-太行山山脉地区，这主要是由于山区耕地较少，且土壤有效磷含量较高，因此磷肥施用量较低。整体上，海河流域内磷肥施用量较高的区域占农田总面积的61.1%。

3.2 农田磷流失迁移因子空间分布

3.2.1 土壤侵蚀模数

海河流域降雨侵蚀力因子（）范围为54～278 MJ·mm/(hm2·h·a)。东西向差异显著，表现出东高西低的规律。高值区主要分布于河北省东北部区域，低值区主要分布于研究区西北部（图3a）。海河流域内值较高的区域占农田总面积的51.9%。而土壤可蚀性因子（）范围为0～0.5 t·hm2·h/(MJ·mm·hm2)（图3b）。高值区主要分布于海河与黄河冲积平原及燕山-太行山山麓地区，海河流域内值较高的区域占农田总面积的67.0%。

地形因子（）范围为0～362。高值区主要分布在燕山-太行山山脉地区，低值区主要分布在平原区域（图 3c）。整体上，海河流域内值较低的区域占农田总面积的91.6%，绝大部分区域地势平缓。而植被覆盖与管理因子（）高值区主要分布于东南部及西北部（图3d）。其中东南部区域值空间分布主要与建设用地的空间分布吻合，具有较强地人为扰动；西北部区域主要是农牧交错区，以草地为主，植被覆盖度比平原耕地与林地低。整体上，海河流域内值较低的区域占农田总面积的77.9%。水土保持措施因子与地形相一致，高值区主要分布在燕山-太行山山脉地区（图3e）。整体上，海河流域内P值较低区域占农田总面积的75.0%。

土壤侵蚀模数（）的分布状况与地形因子（）及植被覆盖与管理因子（）有关。和越大，值愈高。高值区主要分布在研究区的西北部（图3f），其中山地区域由于地势起伏大，易发生侵蚀；草地区域植被覆盖度相对较低，亦易发生侵蚀。整体上，海河流域值较低的区域占农田总面积的99.0%。

3.2.2 年径流深和归一化距离指数

海河流域年径流深（R）的空间分布与降雨侵蚀力基本一致（图4a和图3a），R范围为44.2～99.7 mm/a。R高值区主要分布于东部沿海地区及太行山山脉地区，前者由于沿海，降水量大，从东部沿海到内陆呈降低趋势，这与海陆地带性规律相契合；后者主要是受山地丘陵的地形因素影响，在降水后因山地径流的汇水作用致使径流深度较多。整体上，海河流域R值较高的区域占农田总面积的54.2%。而海河流域归一化距离指数（NDDI）较低的区域占农田总面积的64.8%，即绝大多数的耕地与河网间距离较近（图4b）。

3.3 农田磷流失关键源区及因子识别

源因子指数（Source Factor Index，SI）与磷肥施用量相一致，迁移因子指数（Transport Factor Index，TI）与归一化距离指数（NDDI）相一致，磷指数（Phosphorus Index，PI）与TI相一致（图5）。海河流域SI和TI空间分布图（图5a和5b）中农田磷流失高风险及以上区域分别占农田总面积的40.2%和30.9%（表2）。海河流域农田磷流失高风险区域（即关键源区）占农田总面积的20.5%（表2），主要分布于黄河北岸的山东省和河北省东南部的平原农耕区、海河流域西北部的山区地带（图5c）。

PI、SI及TI与各源因子、迁移因子相关分析结果如表3所示，PI与年径流深、降雨侵蚀力因子显著正相关（<0.05），与有机碳极显著正相关（<0.01）；SI与土壤有效磷、砂粒含量及土壤侵蚀模数、水土保持措施因子呈极显著正相关（<0.01），与年径流深、降雨侵蚀力因子、土壤可蚀性因子、粉粒含量、磷肥施用量极显著负相关（<0.01），与黏粒含量显著负相关（<0.05）；TI与磷肥施用量、年径流深、归一化距离指数、降雨侵蚀力因子、土壤可蚀性因子、黏粒含量呈极显著正相关（<0.01），与粉粒含量呈显著正相关（<0.05），与水土保持措施因子、砂粒含量和有机碳含量均呈极显著负相关（<0.01）。

图3 海河流域土壤侵蚀模数及其影响因子空间分布

图5 海河流域农田磷流失风险指数空间分布

表2 海河流域农田磷流失不同风险等级面积占比

表3 农田磷流失风险指数与各影响因子间相关分析

注：*和**分别表示在0.05和0.01水平上显著相关，ns表示显著不相关（>0.05）。下同。

Note: * and ** indicate significant correlation at 0.05 and 0.01 levels, respectively. And ns represents no significant correlation (>0.05). Same as below.

结构方程模型（SEM）配适度指标（CHI/DF=0.304；GFI=1.000；AGFI=0.999；AIC=18.304；RMSEA=0.000）表明模型拟合效果较好（图6）。SEM结果表明，PI与迁移因子显著正相关（<0.05），源因子与迁移因子间呈极显著负相关（<0.01），土壤性质与源因子呈极显著负相关（<0.01），与迁移因子呈极显著正相关（<0.01）。

注：箭头的宽度表示标准化路径系数强度。实线表示正路径系数，虚线表示负路径系数。

4 讨 论

评估农田土壤磷流失潜在风险对于制定有效的磷污染防控措施至关重要[7]。本研究识别出的农田磷流失关键源区主要集中在海河流域东南部区域的平原农耕区及山区河谷地带，与孔佩儒等研究结果相吻合[35]。这是由于：1）东南部区域的平原农耕区主要为海河与黄河冲积平原，该区域地势平坦（图3c），土壤肥沃，耕作强度及化肥投入量大（图2b），加之农田与河网间距离较近（图 4b），在中高强度的降雨（径流）作用下，土壤及肥料中磷素易随地表径流向受纳水体迁移，导致农田磷流失风险高；2）海河流域西北部的山区地带，由于其地势起伏，坡度较大，土壤侵蚀较为严重，因此农田磷流失风险高（图3c和图3f）。

本研究中，海河流域磷指数（PI）与迁移因子显著正相关（<0.05），而与源因子未达到显著相关水平（图6，>0.05），且PI与TI空间分布高度相关。这说明，海河流域农田磷流失主要影响因子为迁移因子，这与Li等研究结果相一致[11]。从迁移因子角度来看，海河流域土壤侵蚀模数整体上处于极低水平，这与前人研究结果相一致[17]。土壤侵蚀模数值范围较广，这可能是由于自然条件、空间和时间尺度、生态系统相互作用的复杂性及土壤侵蚀地块等因素导致土壤侵蚀模数的高度变异性[36]。而对于中高值区域仍需采取一定的水土保持措施，研究发现保护性耕作和等高种植能增强了土壤对高强度降水的适应力，降低土壤侵蚀模数[37-38]。研究区内年径流深及农田与河网间归一化距离指数高值区所占比例较高，分别达31.0%和64.8%（图4）。且TI和降雨侵蚀力因子、年径流深及归一化距离指数极显著正相关（表 3，<0.01）。这说明径流及农田与河网间距离为主要迁移因子。相关研究表明，降水变化会影响磷投入与输出间的平衡[39]，主要是由于高强度降水会影响年径流深、降雨侵蚀力及磷的相对释放。而海河流域不同区域年降水量变化较大，最大差异可达462.5 mm。本研究中降水量较少的区域（西北部）土壤有效磷含量高于高降水量区域（东南部），这可能是由于低降水区域径流导致的磷损失较低[39]。有研究报道采取一定的水土保持措施及管理实践，有助于减少水土流失及径流作用[40]。从源因子角度来看，SI空间分布结果显示40.2%的农田潜在磷流失达到高风险及以上级别（图5a）。因此，源因子也是影响土壤磷流失的重要因子。SI空间分布与磷肥施用量一致，而处于海河与黄河冲击平原区域的农耕区（山东省及河北省东南部）磷肥施用量普遍过量，因此，在不影响粮食安全的前提下可以通过采取适量减肥、肥料替代及综合养分管理等措施来减少农田磷流失风险[41]。

与此同时，本研究结果表明土壤性质（土壤质地及有机碳含量）与源因子、迁移因子间均呈极显著相关（图6，<0.01）。表3中SI与黏粒含量显著负相关（<0.05），与粉粒含量极显著负相关（<0.01），且SI与土壤有效磷含量、砂粒含量极显著正相关（<0.01）。这可能是由于黏粒或粉粒含量较高的土壤中磷淋洗量少，因此有效磷含量通常比砂粒含量高的土壤中磷含量高[34]。值得注意的是，表3中SI与磷肥施用量呈极显著负相关，然而通过图2b与图5a结果对比分析，二者应呈正相关关系，这一矛盾结果的出现可能是由于：1）SI与土壤有效磷呈极显著正相关（表3，<0.01），而土壤有效磷与磷肥施用量间呈极显著负相关（表3，<0.01）。农业实际生产过程中，土壤有效磷含量较高的地区，磷肥施用量相对来说会较低。反之，土壤有效磷含量较低的地区，磷肥施用量相对来说较高。磷肥施用量一定程度上受制于土壤有效磷含量；2）海河流域范围广，空间异质性较强，诸多地表要素所造成的不确定性因素也会对结果产生一定偏差。综上两点导致相关分析中出现这一结果矛盾现象。TI与土壤可蚀性因子、质地（砂粒、黏粒和粉粒含量）、有机碳和水土保持措施因子显著相关（表 3，<0.05），而这些指标均与土壤性质紧密相关。这说明土壤性质对农田磷流失具有重要影响[42]。大量研究报道，提升土壤质量一定程度上也能有效减少水土流失[42-43]。土壤理化性质会影响土壤中磷形态及分布[44]；较低的耕作频率及土壤质量增强提高了土壤的入渗能力和土壤结构（增加团聚体稳定性），使土壤抗水土流失能力增强[45]。因此，通过一定的保护性耕作等水土保持措施提升土壤质量，缓冲降水（径流）影响及增加土壤粗糙度，减少在一定降雨（径流）强度下农田土壤磷流失对于农田面源污染防治具有重要意义。

农田磷流失的空间变异性通常归因于多个空间变量的综合作用，如土壤磷含量、地理位置、水文条件、降水因素及农业耕作、管理措施等。一些变量间的协同作用会加剧磷流失，而另一些变量间的拮抗作用会减少磷流失。如山区地形和降水因子通常会导致高水土流失和磷流失风险，而茂密植被又可能会减少土壤侵蚀，同时低施肥量也可以通过降低磷盈余来减少磷流失潜力，农田与河网间距离越远越能通过提供更多的磷滞留机会来减少磷流失至受纳水体。平原区域地势起伏较小、坡度平缓，具有较低的侵蚀潜力，然而由于磷投入量较大，在中高强度的降雨（径流）作用下，仍面临高风险级的磷流失[15]，这与本研究中结果相一致。本研究结果发现，源因子与迁移因子间呈极显著负相关关系（图6，<0.01），这可能是由于：1）磷在农田中为作物必需养分，并非污染物；2）影响海河流域农田磷流失的主导因子为迁移因子，而年径流深及归一化距离指数作为迁移因子中的关键因子。地表径流强度越大、农田与河网间距离越近，农田土壤磷流失风险越高，越易使土壤及肥料中的磷迁移到受纳水体造成污染，进而降低SI。此外，在TI与SI双高区域如本研究中的关键源区（黄河北岸山东省区域、河北省东南部的平原农耕区），农田磷流失受源因子及迁移因子的综合作用，因此需要整合该区域关键的源和迁移因子并进行综合考虑，制定出TI与SI双减的最佳管理措施以减少农田磷流失[46]。综上，对于海河流域农田磷流失关键源区中的平原农耕区，可通过采取适量减肥、综合养分管理、保护性耕作及在农田与河网间布设缓冲区等措施来降低农田磷流失风险；对于山区地带，可采取如等高种植、梯田种植等水土保持措施来减少土壤侵蚀，进而降低农田磷流失风险。

本研究风险评估采用经典成熟的磷指数模型，并在参考国内外研究的基础上进行了一定的修正，如改进了磷指数模型中农田到水体间距离因子，首次提出归一化距离因子的概念，并将源因子、迁移因子与磷指数归一化后进行相关分析、建立结构方程模型，研究了土壤性质（如质地及有机碳）、源因子、迁移因子及农田磷流失风险指数间关系，得出最终研究结果。由于磷指数模型中各因子权重的确定存在一定主观性，权重赋值的差异会显著影响PI数值大小，但对于定性化识别关键源区影响不大。研究过程中由于数据受限，在源因子中未考虑施肥时间及方式，这些可能会给研究结果带来一定偏差，因此，对于各因子及磷流失风险评价有待结合实地情况进一步研究。本研究数据均来源于公共平台，易获取，这为流域尺度上直观地识别农田磷流失关键源区和因子，进而采取有针对性的防控措施提供了快速实用的方法。

5 结 论

本研究利用磷指数模型评估了海河流域农田磷流失环境风险，识别出海河流域磷流失的关键源区、因子。研究结果对流域尺度上农田磷流失研究方法创新有所裨益，为海河流域农田面源污染防治提供科学参考。本研究得出结论如下：

1）海河流域农田土壤有效磷、磷肥施用量、土壤侵蚀模数、年径流深及归一化距离指数的空间分布结果表明，大多数区域分别处于中-低、中-高、极低、中-高和高风险级别，分别占农田总面积的66.5%、61.1%、99.0%、54.2%和64.8%；

2）影响海河流域农田磷流失的关键因子为迁移因子，其中关键迁移因子为年径流深及归一化距离指数；

3）源因子与迁移因子间呈极显著负相关（<0.01），土壤性质（包括土壤质地及有机碳含量等）与源因子呈极显著负相关（<0.01），与迁移因子间呈极显著正相关（<0.01）；

4）海河流域农田磷流失关键源区主要位于黄河北岸的山东省和河北省东南部的平原农耕区、海河流域西北部的山区地带，需要重点防控。

[1] Werner F, de la Haye T R, Spielvogel S, et al. Small-scale spatial distribution of phosphorus fractions in soils from silicate parent material with different degree of podzolization[J]. Geoderma, 2017, 302: 52-65.

[2] Granger S J, Harris P, Peukert S, et al. Phosphate stable oxygen isotope variability within a temperate agricultural soil[J]. Geoderma, 2017, 285: 64-75.

[3] Schoumans O F, Chardon W J, Bechmann M E, et al. Mitigation options to reduce phosphorus losses from the agricultural sector and improve surface water quality: A review[J]. Science of the Total Environment, 2014, 468/469: 1255-1266.

[4] Cao D, Cao W, Fang J, et al. Nitrogen and phosphorus losses from agricultural systems in China: A meta-analysis[J]. Marine Pollution Bulletin, 2014, 85(2): 727-732.

[5] Conley D J, Paerl H W, Howarth R W, et al. Controlling eutrophication: Nitrogen and phosphorus[J]. Science, 2009, 323(5917): 1014-1015.

[6] Sun B, Zhang L, Yang L, et al. Agricultural non-point source pollution in China: Causes and mitigation measures[J]. Ambio, 2012, 41(4): 370-379.

[7] Huang L, Wang H, Li Y, et al. Spatial distribution and risk assessment of phosphorus loss potential in urban–suburban soil of Lishui, China[J]. Catena, 2013, 100: 42-49.

[8] Wang F E, Sun Z L, Zheng S Y, et al. An integrated approach to identify critical source areas of agricultural nonpoint-source pollution at the watershed scale[J]. Journal of Environmental Quality, 2018, 47(4): 922-929.

[9] Gburek W J, Sharpley A N, Heathwaite L, et al. Phosphorus management at the watershed scale: A modification of the phosphorus index[J]. Journal of Environment Quality, 2000, 29(1): 130-144.

[10] Sharpley A N, Weld J L, Beegle D B, et al. Development of phosphorus indices for nutrient management planning strategies in the United States[J]. Journal of Soil and Water Conservation, 2003, 58(3): 137-152.

[11] Li Z, Zhang R, Liu C, et al. Phosphorus spatial distribution and pollution risk assessment in agricultural soil around the Danjiangkou reservoir, China[J]. Science of the Total Environment, 2020, 699: 134417.

[12] 刘洁，庞树江，何杨洋，等. 基于小流域产流特征的磷流失关键源区识别[J]. 农业工程学报，2017，33(20)：241-249. Liu Jie, Pang Shujiang, He Yangyang, et al. Critical area identification of phosphorus loss based on runoff characteristics in small watershed[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2017, 33(20): 241-249. (in Chinese with English abstract)

[13] Buchanan B P, Archibald J A, Easton Z M, et al. A phosphorus index that combines critical source areas and transport pathways using a travel time approach[J]. Journal of Hydrology, 2013, 486: 123-135.

[14] 张淑荣，陈利顶，傅伯杰，等. 农业区非点源污染潜在危险性评价：以于桥水库流域磷流失为例[J]. 第四纪研究，2003，23(3)：262-269. Zhang Shurong, Chen Liding, Fu Bojie, et al. The risk assessment of nonpoint pollution of phosphorus from agricultural lands: A case study of Yuqiao Reservoir Watershed[J]. Quaternary Sciences, 2003, 23(3): 262-269. (in Chinese with English abstract)

[15] Zhou H, Gao C. Assessing the risk of phosphorus loss and identifying critical source areas in the Chaohu Lake Watershed, China[J]. Environmental Management, 2011, 48(5): 1033-1043.

[16] 李文超，刘申，雷秋良，等. 高原农业流域磷流失风险评价及关键源区识别：以凤羽河流域为例[J]. 农业环境科学学报，2014，33(8)：1591-1600. Li Wenchao, Liu Shen, Lei Qiuliang, et al. Assessing risk of non-point source phosphorus loss and identifying critical source areas in a Chinese highland agricultural watershed[J]. Journal of Agro-Environment Science, 2014, 33(8): 1591-1600. (in Chinese with English abstract)

[17] 李晓松，吴炳方，王浩，等. 区域尺度海河流域水土流失风险评估[J]. 遥感学报，2011，15(2)：372-387. Li Xiaosong, Wu Bingfang, Wang Hao, et al. Regional soil erosion risk assessment in Hai Basin[J]. Journal of Remote Sensing, 2011, 15(2): 372-387. (in Chinese with English abstract)

[18] 万炜，李含微，王佳莹，等. 基于空间平滑法的旱作区粮食产量时空变化与影响因素研究[J]. 农业工程学报，2019，35(16)：284-296. Wan Wei, Li Hanwei, Wang Jiaying, et al. Spatio-temporal changes and influencing factors of grain yield based on spatial smoothing method in dryland farming regions[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2019, 35(16): 284-296. (in Chinese with English abstract)

[19] Pribyl D W. A critical review of the conventional SOC to SOM conversion factor[J]. Geoderma, 2010, 156(3/4): 75-83.

[20] 鲁如坤. 土壤磷素水平和水体环境保护[J]. 磷肥与复肥，2003，18(1)：4-8. Lu Rukun. The phosphorus level of soil and environmental protection of water body[J]. Phosphate & Compound Fertilizer, 2003, 18(1): 4-8. (in Chinese with English abstract)

[21] 万炜，邓静，王佳莹，等. 基于潜力衰减模型的东北-华北平原旱作区耕地生产力评价[J]. 农业工程学报，2020，36(5)：270－280. Wan Wei, Deng Jing, Wang Jiaying, et al. Evaluation of cultivated land productivity based on potential attenuation model in the dryland farming regions of Northeast and North China Plain[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2020, 36(5): 270－280. (in Chinese with English abstract)

[22] Sharpley A N, Mcdowell R W, Peter J A K. Phosphorus loss from land to water: Integrating agricultural and environmental management[J]. Plant and Soil, 2001, 237(2): 287-307.

[23] Wischmeier W H, Smith D D. Predicting Rainfall Erosion Losses. A guide to conservation planning[M]. Washington: United States Department of Agriculture, 1978.

[24] Williams J, Renard K, Dyke P. EPIC: A new method for assessing erosion’s effect on soil productivity[J]. Journal of Soil and Water Conservation, 1983, 38: 381-383.

[25] Phinzi K, Ngetar N S. The assessment of water-borne erosion at catchment level using GIS-based RUSLE and remote sensing: A review[J]. International Soil and Water Conservation Research, 2019, 7(1): 27-46.

[26] 蔡崇法，丁树文，史志华，等. 应用USLE模型与地理信息系统IDRISI预测小流域土壤侵蚀量的研究[J]. 水土保持学报，2000，14(2)：19-24. Cai Chongfa, Ding Shuwen, Shi Zhihua, et al. Study of applying USLE and geographical information system IDRISI to predict soil erosion in small watershed[J]. Journal of Soil and Water Conservation, 2000, 14(2): 19-24. (in Chinese with English abstract)

[27] Fu B J, Zhao W W, Chen L D, et al. Assessment of soil erosion at large watershed scale using RUSLE and GIS: A case study in the Loess Plateau of China[J]. Land Degradation & Development, 2005, 16(1): 73-85.

[28] Lufafa A, Tenywa M M, Isabirye M, et al. Prediction of soil erosion in a Lake Victoria basin catchment using a GIS-based Universal Soil Loss model[J]. Agricultural Systems, 2003, 76(3): 883-894.

[29] 任宪韶，户作亮，曹寅白. 海河流域水利手册[M]. 北京：中国水利水电出版社，2008.

[30] 刘春蓁，刘志雨，谢正辉. 近50年海河流域径流的变化趋势研究[J]. 应用气象学报，2004，15(4)：385-393. Liu Chunzhen, Liu Zhiyu, Xie Zhenghui. Study on trends in runoff for the Haihe River Basin in recent 50 years[J]. Journal of Applied Meteorological Science, 2004, 15(4): 385-393. (in Chinese with English abstract)

[31] 张世文，王朋朋，叶回春，等. 基于数字土壤系统的县域土壤磷素流失风险简化评估[J]. 农业工程学报，2012，28(11)：110-117. Zhang Shiwen, Wang Pengpeng, Ye Huichun, et al. Risk simplified assessment on phosphorus loss risk based on digital soil system at county scale[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2012, 28(11): 110-117. (in Chinese with English abstract)

[32] Bechmann M E, Stålnacke P, Kværnø S H. Testing the Norwegian phosphorus index at the field and subcatchment scale[J]. Agriculture, Ecosystems & Environment, 2007, 120(2/3/4): 117-128.

[33] Chacón N, Dezzeo N, Muǹoz B, et al. Implications of soil organic carbon and the biogeochemistry of iron and aluminum on soil phosphorus distribution in flooded forests of the lower Orinoco River, Venezuela[J]. Biogeochemistry, 2005, 73(3): 555-566.

[34] Momtaz H R, Ghasemi N, Sepehr E, et al. Spatial distribution of soil-available phosphorus and single-point phosphorus buffering index in the Khoy region, Iran[J]. Soil Use and Management, 2017, 33(1): 13-24.

[35] 孔佩儒，陈利顶，孙然好，等. 海河流域面源污染风险格局识别与模拟优化[J]. 生态学报，2018，38(12)：4445-4453. Kong Peiru, Chen Liding, Sun Ranhao, et al. Identification and optimized simulation of non-point source pollution risk patterns in the Haihe River Basin[J]. Acta Ecologica Sinica, 2018, 38(12): 4445-4453. (in Chinese with English abstract)

[36] Boix F C, Martínez M M, Calvo C A, et al. Concise review of interrill erosion studies in SE Spain (Alicante and Murcia): Erosion rates and progress of knowledge from the 1980s[J]. Land Degradation & Development, 2005, 16(6): 517-528.

[37] Liu R M, Zhang P P, Wang X J, et al. Assessment of effects of best management practices on agricultural non-point source pollution in Xiangxi River watershed[J]. Agricultural Water Management, 2013, 117: 9-18.

[38] Beniston J W, Shipitalo M J, Lal R, et al. Carbon and macronutrient losses during accelerated erosion under different tillage and residue management[J]. European Journal of Soil Science, 2015, 66(1): 218-225.

[39] Liu Z, Shao M, Wang Y. Spatial patterns of soil total nitrogen and soil total phosphorus across the entire Loess Plateau region of China[J]. Geoderma, 2013, 197/198: 67-78.

[40] Eekhout J P C, de Vente J. Assessing the effectiveness of Sustainable Land Management for large-scale climate change adaptation[J]. Science of the Total Environment, 2019, 654: 85-93.

[41] Heathwaite L, Sharpley A, Gburek W. A conceptual approach for integrating phosphorus and nitrogen management at watershed scales[J]. Journal of Environment Quality, 2000, 29(1): 158-166.

[42] Martínez-Mena M, Carrillo-López E, Boix-Fayos C, et al. Long-term effectiveness of sustainable land management practices to control runoff, soil erosion, and nutrient loss and the role of rainfall intensity in Mediterranean rainfed agroecosystems[J]. Catena, 2020, 187: 104352.

[43] Ouyang W, Huang H B, Hao F H, et al. Evaluating spatial interaction of soil property with non‐point source pollution at watershed scale: The phosphorus indicator in Northeast China[J]. Science of the Total Environment, 2012, 432: 412-421.

[44] Bechmann M, Stalnacke P, Kvoerno S, et al. Integrated tool for risk assessment in agricultural management of soil erosion and losses of phosphorus and nitrogen[J]. Science of the Total Environment, 2009, 407(2): 749-759.

[45] Haruna S I, Anderson S H, Nkongolo N V, et al. Soil hydraulic properties: Influence of tillage and cover crops[J]. Pedosphere, 2018, 28(3): 430-442.

[46] Rao N S, Easton Z M, Schneiderman E M, et al. Modeling watershed-scale effectiveness of agricultural best management practices to reduce phosphorus loading[J]. Journal of Environmental Management, 2009, 90(3): 1385-1395.

Environmental risk assessment of phosphorus loss from farmland based on phosphorus index model in the Haihe River Basin

Wu Hanqing1, Wan Wei1, Shan Yanjun1, Chen Yanhua2, Li Qiang1, Li Cunjun3, Hu Haitang3, Zhang Baogui1※

(1.,,100193,; 2.,,100097,; 3.,100097,)

The non-point phosphorus (P) loss from farmland is the main cause of eutrophication, and identifying the Critical Source Areas (CSAs) and Critical Influence Factors (CIFs) of the P loss from farmland is a pivotal part of controlling agricultural non-point source pollution. The objective of this study was to assess the environmental risk of P loss from farmland based on P index model in the Haihe River Basin. In this study, we took the Haihe River Basin as the research region, and applied the P index model. Thereinto, the soil available P content and fertilizer-P application rate were selected as the source factors. Synchronously, the soil erosion modulus, annual runoff, the Normalized Differential Distance Index (NDDI) between farmland and river network were chosen as the transport factors. Additionally, this research aims at identifying the CSAs combined with GIS technology, distinguishing the CIFs and exploring the relation between P index and influence factors such as soil properties (mainly including soil texture and organic carbon), source and transport factors via Structural Equation Modeling (SEM). The spatial distributions of soil available P content, fertilizer-P application rate, soil erosion modulus, annual runoff and NDDI in the Haihe River Basin indicated that the most areas were at the risk levels of low-medium, medium-high, very low, medium-high and high in the Haihe River Basin, respectively. Namely, these areas accounted for 66.5%, 61.1%, 99.0%, 54.2%, and 64.8% of the total farmland area in the Haihe River Basin, respectively. Moreover, we revealed that the CIF affecting P loss from farmland in the Haihe River Basin is the transport factors in which the key transport factors are the annual runoff and NDDI. Furthermore, SEM prescribed that there was an extremely significant negative correlation between the source factors (<0.01) and extremely significant positive correlation with transport factors (<0.01), and soil properties had highly significant negative correlations with source factors and transport factors (<0.01). Subsequently, the findings illustrated that the CSAs of P loss from farmland in the Haihe River Basin were mainly located in the plain farming areas of Shandong Province and the southeast of Hebei Province on the north bank of the Yellow River, and the mountain areas in the northwestern of the Haihe River Basin. Due to the plain farming areas have high cultivation intensity and fertilizer-P application rate, flat terrain, and close distance to river network, P of soil and fertilizer could easily migrate into the rivers following the surface runoff with the medium-high intensity. Thus, for the plain farming areas, some practices should be taken to decrease the risk level of P loss from farmland, such as appropriate reduction of the fertilizer-P application rate, comprehensive nutrient management, conservation tillage and deployment of buffer zone between farmland and river network and so on. Because mountain areas have severe soil erosion, therefore, the risk level of P loss from farmland is relatively high. For this, proper soil and water conservation measures should be taken to reduce the soil erosion and risk level of P loss from farmland, such as contour tillage and terrace cropping and so on. In conclusion, the results are beneficial to the method innovation of P loss from farmland at the watershed scale, and can provide a scientific reference for controlling agricultural non-point source pollution in the Haihe River Basin.

farmland; models; non-point source pollution; phosphorus loss; phosphorus index model; critical source area; structural equation modeling; Haihe River Basin

吴汉卿，万炜，单艳军，等. 基于磷指数模型的海河流域农田磷流失环境风险评价[J]. 农业工程学报，2020，36(14)：17-27.doi：10.11975/j.issn.1002-6819.2020.14.003 http://www.tcsae.org

Wu Hanqing, Wan Wei, Shan Yanjun, et al. Environmental risk assessment of phosphorus loss from farmland based on phosphorus index model in the Haihe River Basin[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2020, 36(14): 17-27. (in Chinese with English abstract) doi：10.11975/j.issn.1002-6819.2020.14.003 http://www.tcsae.org

2020-06-22

2020-07-10

国家重点研发计划项目（2016YFD0800906）

吴汉卿，博士生，主要从事农田面源氮磷污染研究。Email：wuhanqing920705@163.com

张宝贵，教授，博士生导师，主要从事农田面源污染和土壤重金属污染对生态系统服务影响研究。Email：zhangbg@cau.edu.cn

10.11975/j.issn.1002-6819.2020.14.003

X52；K903

A

1002-6819(2020)-14-0017-11