海河流域農田氮磷面源污染的空間分布特征及關鍵源區識別

張巧玲，胡海棠，王道蕓，邱春霞，李存軍，白 翠，靖亭亭

（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 黑色圓圈內）與入河負荷熱點區的空間分布格局不一致。因此，面源污染防控的關鍵源區識別應將入河負荷和潛在氮磷徑流質量濃度的空間分布結合進行分析。對于潛在氮磷徑流質量濃度高的河段，即使其集水區即源區位于入河負荷的非熱點區，其治理優先序靠前。例如大清河水系東段支流的氮徑流質量濃度很高，其源區以入河負荷次熱點和不顯著區為主，但面源污染防控時仍應重點治理。在選擇優先防控區域時，應著重考慮對內源污染為主的潛在氮磷徑流質量濃度高的區域進行源頭控制，而對外源污染為主的潛在氮磷徑流質量濃度高的區域可以采用過程攔截和末端治理的方式，并追溯其上游的關鍵源區。

4 結論

1）2015年海河流域農田氮、磷入河負荷分別為2.41、0.56 kg/hm2，總負荷分別為3.493 4萬、0.807 7萬t，潛在氮、磷徑流質量濃度分別為5.97、1.47 mg/L，有55%的河段潛在氮磷徑流質量濃度超過國家地表V類水質TN、TP標準。

2）農田氮磷污染物呈明顯的沿西北部山區向中南部平原方向上升的分布特征。徒駭馬頰河水系、漳衛河水系、子牙河水系和黑龍港運東河水系是農田氮磷入河負荷和總負荷較高的水系。

3）農田面源污染引起的河流斷面潛在氮磷徑流質量濃度高值區主要分布在徒駭馬頰河全線、黑龍港運東河中下游、大清河支流和子牙河上游部分河段等，流域水系上游或支流的潛在氮磷徑流質量濃度普遍高于所屬水系干流。

4）農田氮磷入河負荷關鍵源區主要分布在流域中部和南部、石家莊、秦皇島等周邊地區，需要重點防控。