基于磷指數模型的海河流域農田磷流失環境風險評價

吳漢卿，萬 煒，單艷軍，陳延華，李 強，李存軍，胡海棠，張寶貴

基于磷指數模型的海河流域農田磷流失環境風險評價

吳漢卿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