海河流域農(nóng)田氮磷面源污染的空間分布特征及關(guān)鍵源區(qū)識(shí)別

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

（1.西安科技大學(xué) 測(cè)繪科學(xué)與技術(shù)學(xué)院，西安710054；2.北京農(nóng)業(yè)信息技術(shù)研究中心，北京100097；3.國家農(nóng)業(yè)信息化工程技術(shù)研究中心，北京100097）

0 引 言

【研究意義】海河流域是我國水環(huán)境污染最嚴(yán)重的流域之一，其中農(nóng)業(yè)面源污染是海河流域水質(zhì)污染的重要來源[1]。2015年《中國環(huán)境狀況公報(bào)》顯示，海河流域?yàn)橹囟任廴?，地表水中V類及劣V類占60%。海河流域總面積為3.18×105km2，占全國總面積的3.3%，水資源僅占全國1.3%，卻承載全國8%的耕地和10%的人口，貢獻(xiàn)了全國糧食總產(chǎn)量的9.4%[2]。作為我國北方糧食主產(chǎn)區(qū)之一，海河流域化肥施用量一直居高不下，導(dǎo)致主要河流水質(zhì)惡化，為該區(qū)域的面源污染治理帶來了挑戰(zhàn)[3]。鑒于此，研究海河流域農(nóng)田氮磷面源污染的空間分布特征和識(shí)別關(guān)鍵源區(qū)，對(duì)海河流域農(nóng)業(yè)面源污染防控具有重要意義。

【研究進(jìn)展】面源污染負(fù)荷空間分布研究目前多采用風(fēng)險(xiǎn)評(píng)價(jià)或模型模擬的方法[4]。相比以氮磷指數(shù)法為代表的風(fēng)險(xiǎn)評(píng)價(jià)方法，模型模擬方法能夠更為準(zhǔn)確地進(jìn)行面源污染負(fù)荷定量化評(píng)估。國內(nèi)外研究主要采用SWAT 模型[5-6]、輸出系數(shù)模型（ECM）[7-9]、AnnAGNPS[10]和平均濃度法[11]等模型方法。其中SWAT 模型、AnnAGNPS 等機(jī)理模型由于對(duì)數(shù)據(jù)要求高，計(jì)算復(fù)雜，不適合大區(qū)域模擬[1]。InVEST 是國際上應(yīng)用最廣泛的生態(tài)系統(tǒng)服務(wù)評(píng)估集成平臺(tái)，其中的營養(yǎng)物傳輸率模型（NDR）所需參數(shù)比較簡單，機(jī)理清晰，適宜用于大尺度面源污染模擬研究[12]。

關(guān)鍵源區(qū)識(shí)別是面源污染管理的重要組成部分[13]。研究表明，少數(shù)區(qū)域的輸出污染物往往貢獻(xiàn)整個(gè)流域污染負(fù)荷的大部分，對(duì)水體的環(huán)境質(zhì)量起著決定性作用[14]，這些區(qū)域被稱為關(guān)鍵源區(qū)。目前關(guān)鍵源區(qū)識(shí)別方法主要是基于面源污染負(fù)荷的空間分布異質(zhì)性特征[15-16]。面源污染負(fù)荷產(chǎn)生和輸移具有不同的空間分布特征[17-18]。李文超等[19]提出應(yīng)依據(jù)入河負(fù)荷的空間分布識(shí)別關(guān)鍵源區(qū)。關(guān)鍵源區(qū)識(shí)別在很大程度上受污染源影響，但與水質(zhì)等因素也是高度相關(guān)的[20]。一方面，農(nóng)田入河負(fù)荷高的地方，水質(zhì)污染不一定嚴(yán)重，后者受景觀結(jié)構(gòu)和水資源狀況的共同影響。另一方面，水質(zhì)污染嚴(yán)重的地區(qū)，不一定是因?yàn)檗r(nóng)田污染源，需要將農(nóng)田對(duì)水質(zhì)的影響剝離出來?！厩腥朦c(diǎn)】目前農(nóng)田面源污染負(fù)荷對(duì)河流水質(zhì)的潛在影響和空間分布特征還缺乏研究。依據(jù)《中國環(huán)境狀況公報(bào)》，海河流域水質(zhì)狀況存在顯著的空間異質(zhì)性，部分支流河段水質(zhì)重度污染，而所屬干流呈現(xiàn)輕度污染[21]。農(nóng)田面源污染對(duì)水質(zhì)的影響是否也具有相似的分布規(guī)律，還有待進(jìn)一步研究。在防控資源有限的條件下，為了進(jìn)一步提高農(nóng)田面源污染優(yōu)先防控區(qū)域選擇的精準(zhǔn)程度，需要研究農(nóng)田面源污染負(fù)荷對(duì)河流水質(zhì)影響程度的評(píng)估方法。結(jié)合面源污染負(fù)荷的空間分布特征，選擇水質(zhì)潛在影響嚴(yán)重的河段并追溯污染物貢獻(xiàn)率大的源區(qū)進(jìn)行關(guān)鍵源區(qū)識(shí)別。

【擬解決的關(guān)鍵問題】綜上，本研究以海河流域?yàn)檠芯繀^(qū)，利用InVEST 營養(yǎng)物傳輸率NDR 模型和產(chǎn)水量模型，對(duì)海河流域農(nóng)田氮磷入河負(fù)荷和產(chǎn)水量進(jìn)行空間估算，并基于水文網(wǎng)絡(luò)拓?fù)潢P(guān)系估算和分析河流斷面氮磷入河通量和潛在氮磷徑流質(zhì)量濃度的空間分布特征，結(jié)合GIS 的空間熱點(diǎn)分析方法和水文網(wǎng)絡(luò)方法識(shí)別關(guān)鍵源區(qū)，以期為海河流域農(nóng)田面源污染精準(zhǔn)防控提供參考依據(jù)。

1 材料與方法

1.1 研究區(qū)概況

海河流域地處我國華北地區(qū)（112°—120°E、30°—43°N），包括北京、天津、河南北部、山西東部、山東北部、內(nèi)蒙古和遼寧小部分和河北的絕大部分。流域總的地勢(shì)是“西北高、東南低”，包括灤河及冀東沿海、北三河、永定河、大清河、子牙河、黑龍港運(yùn)東河、漳衛(wèi)南運(yùn)河、徒駭馬頰河八大水系[22]（圖1）。整個(gè)流域的水系從南到北呈扇形分布，水系分散、支流較多。海河流域?qū)贉貛Т箨懶约撅L(fēng)氣候，屬半濕潤半干旱地帶，年均降水量為540 mm[23]。土地利用類型以耕地為主，耕地面積占流域面積的48.9%。海河流域北部及西部山區(qū)為一熟農(nóng)區(qū)，中部和南部平原地區(qū)以冬小麥-玉米二熟種植為主，是中國三大糧食生產(chǎn)基地之一。

1.2 數(shù)據(jù)來源與預(yù)處理

本研究海河流域農(nóng)田氮磷入河負(fù)荷和潛在氮磷徑流質(zhì)量濃度采用InVEST的NDR模型和產(chǎn)水量模型估算。NDR模型運(yùn)行所需數(shù)據(jù)包括海河流域月降水量、土地利用、數(shù)字高程模型（DEM）、坡度、農(nóng)區(qū)分布、坡度分布、子流域分布圖、地表氮磷負(fù)荷等。產(chǎn)水量模型運(yùn)行所需數(shù)據(jù)包括土地利用、月降水量、潛在蒸散量、土壤有效含水率等。具體數(shù)據(jù)如下：

1）氣象數(shù)據(jù)。在中國氣象數(shù)據(jù)中心下載海河流域相關(guān)氣象站點(diǎn)2001—2015年降水量、氣溫、風(fēng)速、日照時(shí)間等數(shù)據(jù)。利用python計(jì)算海河流域所有氣象站點(diǎn)多年平均降水量，并利用克里格插值方法得到流域多年平均降水量空間分布圖。由Penman-Monteith公式計(jì)算潛在蒸散量，利用反距離加權(quán)空間插值方法得到多年平均潛在蒸散發(fā)空間分布圖。

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

2）DEM和2015年土地利用數(shù)據(jù)來源于中國科學(xué)院地理科學(xué)與資源研究所，空間分辨率均為30 m?；贒EM數(shù)據(jù)，利用ArcGIS水文分析模塊提取研究區(qū)的子流域及水文網(wǎng)絡(luò)空間分布。整個(gè)海河流域被劃分為263個(gè)子流域，子流域出口設(shè)置河流斷面，共273個(gè)斷面，259個(gè)河段。

3）中國1∶100萬土壤數(shù)據(jù)庫。來自中國西部環(huán)境與生態(tài)科學(xué)數(shù)據(jù)中心。土壤有效含水率在SPAW軟件中基于土壤質(zhì)地及粒徑組成數(shù)據(jù)，利用土壤有效含水率經(jīng)驗(yàn)公式計(jì)算。

4）農(nóng)區(qū)數(shù)據(jù)。海河流域農(nóng)區(qū)分布數(shù)據(jù)來源于《中國農(nóng)作制》[24]，在ArcGIS中矢量化得到海河流域農(nóng)區(qū)的空間分布圖。

5）統(tǒng)計(jì)數(shù)據(jù)。統(tǒng)計(jì)數(shù)據(jù)包括農(nóng)作物播種面積、氮磷肥施用量、耕地面積等，來源于相關(guān)各省和縣市統(tǒng)計(jì)年鑒。

按照《全國農(nóng)田面源污染排放系數(shù)手冊(cè)》，遵循“農(nóng)區(qū)-坡度-種植模式”原則，根據(jù)農(nóng)區(qū)分布圖和各縣農(nóng)作物播種面積數(shù)據(jù)，將海河流域農(nóng)田劃分為13種面源污染模式（表1），并進(jìn)行空間匹配（詳見文章[25]），得到海河流域面源污染模式分布圖。

表1 海河流域農(nóng)田面源污染模式Table 1 Farmland non-point source pollution models in the Haihe River Basin

1.3 研究方法

1.3.1 農(nóng)田氮磷入河負(fù)荷估算模型

農(nóng)田面源污染氮磷入河負(fù)荷的估算分2步進(jìn)行，首先基于農(nóng)田面源污染模式及借鑒輸出系數(shù)模型，估算農(nóng)田氮磷流失負(fù)荷；然后基于InVEST 3.8.7的NDR模型，計(jì)算氮磷輸移率和氮磷入河負(fù)荷。

1）氮磷流失負(fù)荷模型

借鑒輸出系數(shù)模型，農(nóng)田化肥施用氮磷流失負(fù)荷估算在柵格單元上進(jìn)行，計(jì)算式為：

式中：Loadi,j為化肥施用氮磷流失負(fù)荷（kg）；i代表第i個(gè)縣市，j代表第j類面源污染模式；Ej為柵格單元對(duì)應(yīng)的面源污染模式j(luò)的肥料氮或磷流失系數(shù)（%），在《全國農(nóng)田面源污染排放系數(shù)手冊(cè)》中查找面源污染模式j(luò)對(duì)應(yīng)的農(nóng)田氮磷流失系數(shù)；Qi,j為柵格單元對(duì)應(yīng)的縣市i的面源污染模式j(luò)對(duì)應(yīng)的氮磷肥施用折純量（kg/hm2），根據(jù)調(diào)研各種植模式的氮磷肥常規(guī)用量和各縣市農(nóng)作物播種面積、氮磷肥總量等進(jìn)行區(qū)域化修正。

2）氮磷入河負(fù)荷模型

NDR模型是基于質(zhì)量守恒方法模擬氮、磷營養(yǎng)物在空間上的遷移過程，估算氮磷輸移率和氮磷入河負(fù)荷。該模型可以模擬流域營養(yǎng)物的來源和輸移過程，計(jì)算農(nóng)田氮磷養(yǎng)分通過地表徑流和壤中流遷移進(jìn)入河流的傳輸率（即入河系數(shù)），從而估算海河流域農(nóng)田氮磷入河負(fù)荷。本文對(duì)模型進(jìn)行修正，將流失負(fù)荷模型和NDR模型結(jié)合起來。模型基本原理如下：

首先，對(duì)農(nóng)田氮磷流失負(fù)荷按照一定比例進(jìn)行地表和次表層分配。

式中：load_xi為農(nóng)田氮磷流失負(fù)荷；prportionsubsurface為次表層營養(yǎng)物占比參數(shù)；loadsurf,i為地表營養(yǎng)物負(fù)荷；loadsubsurf,i為次表層營養(yǎng)物負(fù)荷。

再次，計(jì)算地表營養(yǎng)物傳輸率：

式中：NDR0,i被下游像素保留的營養(yǎng)物傳輸率，ICj為地形指數(shù)，IC0和k是校準(zhǔn)參數(shù)。

計(jì)算次表層營養(yǎng)物傳輸率：

式中：NDRsubs,i為地下營養(yǎng)物傳輸率；effsubs是能夠達(dá)到次表層徑流的營養(yǎng)物最大截留效率；lsubs是次表層的截留長度，即假設(shè)土壤保持營養(yǎng)物的最大距離，l是從像素到流的距離。

氮磷入河負(fù)荷計(jì)算式為：

式中：xexp,i為各柵格單元i的入河負(fù)荷；loadsurf,i為地表營養(yǎng)物負(fù)荷；NDRsurf,i為地表營養(yǎng)物傳輸率；loadsub,i為次表層營養(yǎng)物負(fù)荷；NDRsub,i為次表層營養(yǎng)物傳輸率；xexptot為子流域營養(yǎng)物入河總負(fù)荷。

其中模型參數(shù)取值參考了InVEST用戶手冊(cè)[26]。NDR模型重要參數(shù)取值及依據(jù)如下：①次表層營養(yǎng)物來源占比，根據(jù)相關(guān)文獻(xiàn)中海河流域地表和次表層徑流分配狀況，取值為0.5；②氮磷營養(yǎng)物最大截留效率由土地利用類型決定，本文林地、草地、耕地、水域、建設(shè)用地、未利用地的氮最大截留效率分別為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用戶手冊(cè)；③次表層的截留長度按照InVEST模型用戶手冊(cè)建議設(shè)定為1個(gè)柵格單元大小，即30 m；④被下游像素保留的營養(yǎng)物傳輸率和校準(zhǔn)參數(shù)IC0具體計(jì)算參見InVEST用戶手冊(cè)，校準(zhǔn)參數(shù)k默認(rèn)值為2。

1.3.2 水文網(wǎng)絡(luò)拓?fù)潢P(guān)系構(gòu)建和氮磷入河通量

流域水文網(wǎng)絡(luò)構(gòu)建以河段為紐帶，是流域匯流關(guān)系建立的基礎(chǔ)。利用ArcGIS 的水文分析對(duì)流域的河段、徑流節(jié)點(diǎn)和子流域等要素提取。參考蔣好忱[27]的方法，借助Matlab 編程構(gòu)建整個(gè)流域的河網(wǎng)拓?fù)潢P(guān)系。灤河水系中上游區(qū)域水文網(wǎng)絡(luò)拓?fù)潢P(guān)系的過程見圖2。

基于氮磷入河負(fù)荷柵格數(shù)據(jù)，采用ArcGIS 的分區(qū)統(tǒng)計(jì)功能，得到每個(gè)子流域的氮磷入河總負(fù)荷，并基于水文網(wǎng)絡(luò)拓?fù)潢P(guān)系，利用Matlab 查詢每個(gè)河流斷面處匯入的子流域清單，得到斷面處的氮磷入河通量。

圖2 水文網(wǎng)絡(luò)（左）及子流域匯流拓?fù)潢P(guān)系（右）Fig.2 Hydrological network(left)and topological relationship between sub-basin（right）

1.3.3 潛在氮磷徑流質(zhì)量濃度估算

本研究通過構(gòu)建潛在氮磷徑流質(zhì)量濃度指標(biāo)來表征農(nóng)田氮磷流失對(duì)河流斷面處水質(zhì)的影響程度。因缺乏海河流域各河流斷面的徑流量數(shù)據(jù)，故采用InVEST的產(chǎn)水量模型估算產(chǎn)水量，相當(dāng)于潛在的最大徑流量。潛在氮磷徑流質(zhì)量濃度定義為河流斷面處氮磷入河通量與斷面上游產(chǎn)水量的比值。

產(chǎn)水量模型基于水量平衡原理，是基于Budyko假設(shè)和年平均降水量的模型。將每個(gè)柵格單元的產(chǎn)水量定義為該柵格的降水量減去蒸散發(fā)量（包括地表蒸發(fā)和植物蒸騰）[28]。模型主要算法為：

式中：Yxj為年產(chǎn)水量（mm）；Px為年均降水量（mm）；AETxj為土地利用類型j上柵格單元x的年平均實(shí)際蒸散量（mm）。

式中：Rxj為Bydydo 干燥指數(shù)，無量綱；kxj為植被蒸散系數(shù)；ET0為潛在蒸散發(fā)，由Penman-Monteith 公式計(jì)算[29]。Z為季節(jié)因子，根據(jù)研究區(qū)的降水分布情況確定，降水分布均勻和夏季降水為主的地區(qū)，Z取值接近1，若冬季降水為主，則Z值接近10[30]。根據(jù)海河流域降水量分布情況Z取值為2。AWCx為柵格單元x的植被可利用水量，由根系深度、土層厚度和土壤有效含水率共同決定。

為了反映海河流域產(chǎn)水量多年平均水平，減少年際間氣候差異，尤其是降水量差異對(duì)模型不確定性的影響，降水量和潛在蒸散量采用2001—2015年的多年均值數(shù)據(jù)。

1.3.4 氮磷關(guān)鍵源區(qū)識(shí)別

本文農(nóng)田面源污染關(guān)鍵源區(qū)識(shí)別是在氮磷入河負(fù)荷和潛在氮磷徑流質(zhì)量濃度的熱點(diǎn)分析基礎(chǔ)上結(jié)合管理單元貢獻(xiàn)率確定關(guān)鍵源區(qū)。

氮磷入河負(fù)荷的熱點(diǎn)分析是利用ArcGIS 中熱點(diǎn)分析（hot spots）工具，首先以子流域或縣域等為單元進(jìn)行分區(qū)統(tǒng)計(jì)得到總負(fù)荷或單位面積平均負(fù)荷，然后基于Getis-Ord Gi*指數(shù)進(jìn)行冷點(diǎn)和熱點(diǎn)類型分區(qū)[25]，Getis-Ord Gi*是一種識(shí)別高值或低值要素在空間上發(fā)生聚類的位置的空間統(tǒng)計(jì)方法，原理見式（8）。其中Gi*指數(shù)值為3 代表熱點(diǎn)區(qū)、2 代表次熱點(diǎn)區(qū)、1代表不顯著、-2～0 代表次冷點(diǎn)區(qū)、-3 代表冷點(diǎn)區(qū)，本文將熱點(diǎn)區(qū)視為關(guān)鍵源區(qū)。

式中：Xj為斑塊j的值；為所有斑塊的均值；Wij為斑塊i和j的空間權(quán)重；n為斑塊總數(shù)。

河段氮磷潛在徑流質(zhì)量濃度的熱點(diǎn)分析，采用GIS 的Jenks 自然斷點(diǎn)分級(jí)法劃分不同潛在氮磷徑流質(zhì)量濃度等級(jí)。此外，采用貢獻(xiàn)率方法定量評(píng)估各個(gè)子流域區(qū)對(duì)河流斷面水質(zhì)的貢獻(xiàn)。貢獻(xiàn)率方法是計(jì)算管理單元的污染貢獻(xiàn)率，即針對(duì)某個(gè)河流斷面，首先通過水文網(wǎng)絡(luò)分析方法回溯其上游子流域或行政單元清單，計(jì)算斷面上游污染物總負(fù)荷、各管理單元的總貢獻(xiàn)率或單位面積貢獻(xiàn)率，然后依據(jù)貢獻(xiàn)率大小進(jìn)行關(guān)鍵源區(qū)識(shí)別。

2 結(jié)果與分析

2.1 海河流域農(nóng)田氮磷入河負(fù)荷空間分布格局

海河流域農(nóng)田氮磷入河負(fù)荷見圖3，農(nóng)田氮、磷入河負(fù)荷均值分別為2.41 kg/hm2和0.56 kg/hm2，總負(fù)荷分別為3.493 4 萬、0.807 7 萬t。

在空間分布上，氮入河負(fù)荷整體呈現(xiàn)出“西北低，東南高”的特征，流域的中部和南部地區(qū)的氮磷入河負(fù)荷明顯高于西北部山區(qū)；磷入河負(fù)荷空間分布無明顯規(guī)律性，高值分布在海河流域南部和東北部的秦皇島和唐山。西北部山區(qū)的氮、磷入河負(fù)荷均值分別為1.58、0.58 kg/hm2，中部和南部平原地區(qū)的氮、磷入河負(fù)荷均值分別為2.75、0.54 kg/hm2，即平原地區(qū)氮入河負(fù)荷顯著高于山區(qū)，而磷入河負(fù)荷稍低于山區(qū)。山區(qū)種植模式以一年一熟為主，磷肥施用量遠(yuǎn)低于一年二熟為主的平原地區(qū)，但因山區(qū)土壤侵蝕強(qiáng)度高，磷流失系數(shù)和遷移系數(shù)較高，導(dǎo)致其磷輸出仍然偏高。

圖3 農(nóng)田氮磷入河負(fù)荷Fig.3 The export load of nitrogen and phosphorus from farmland

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

以子流域?yàn)閱卧y(tǒng)計(jì)農(nóng)田氮磷入河負(fù)荷，并基于水文網(wǎng)絡(luò)拓?fù)潢P(guān)系，匯總每個(gè)河流斷面的上游所有子流域入河總負(fù)荷，即得到各河流斷面的氮磷入河通量（圖4）。氮磷入河通量隨河流方向自西向東累積，位于流域中東部的黑龍港運(yùn)東河和子牙河下游河段的入河通量較高。

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

海河流域八大水系的氮磷入河負(fù)荷和入河通量的空間統(tǒng)計(jì)結(jié)果詳見表2。灤河水系和永定河水系農(nóng)田氮入河負(fù)荷均值分別為1.24 kg/hm2和1.35 kg/hm2，其他六大水系均高于2.48 kg/hm2，其中徒駭馬頰河水系高達(dá)2.92 kg/hm2。磷入河負(fù)荷高值區(qū)分布在漳衛(wèi)河水系和徒駭馬頰水系，分別為1.07 kg/hm2和1.43 kg/hm2。氮入河通量較高的水系包括子牙河水系、黑龍港運(yùn)東河水系、永定河水系、漳衛(wèi)河水系，均超過5 000 t；而磷入河通量較高的區(qū)域主要分布在永定河水系、黑龍港運(yùn)東河水系、漳衛(wèi)河水系、子牙河水系、徒駭馬頰河水系，均超過1 000 t。

綜上，徒駭馬頰河水系、漳衛(wèi)河水系、子牙河水系是海河流域八大水系中農(nóng)田氮磷入河負(fù)荷較高的水系，漳衛(wèi)河水系、子牙河水系、永定河水系、黑龍港運(yùn)東河水系是農(nóng)田氮磷污染物總量較大的水系。

2.2 海河流域潛在氮磷徑流質(zhì)量濃度空間分布格局

海河流域潛在氮徑流質(zhì)量濃度（圖5）范圍為0.07～23.96 mg/L，均值為5.97 mg/L，超過國家地表水V類限值（2 mg/L）3倍；潛在磷徑流質(zhì)量濃度范圍為0.02～5.86 mg/L，均值為1.47 mg/L，超過國家地表水V類限值（0.4 mg/L）3.7倍，約55%的河段超過了國家地表水V類水質(zhì)標(biāo)準(zhǔn)，其分布與2015年《中國環(huán)境狀況公報(bào)》中海河流域水質(zhì)分布狀況較為一致[21]。這表明農(nóng)田面源污染是海河流域重要的污染源，即使不考慮其他污染源的影響，農(nóng)田氮磷流失仍然可能造成近1/2的河段處于嚴(yán)重污染狀態(tài)。

由圖5可知，海河流域潛在氮磷徑流質(zhì)量濃度呈沿山區(qū)向平原方向上升趨勢(shì)，流域中部和南部的污染比較嚴(yán)重，高值區(qū)主要分布在徒駭馬頰河全線、黑龍港運(yùn)東河中下游、大清河支流和子牙河上游部分河段等。流域水系上游或支流的潛在氮磷徑流質(zhì)量濃度普遍高于干流。以國家地表水V類限值為標(biāo)準(zhǔn)，支流潛在氮磷徑流質(zhì)量濃度超標(biāo)河段個(gè)數(shù)分別占支流總數(shù)的62%、71%，而干流超標(biāo)數(shù)分別占干流總數(shù)的43%、47%。這主要是因?yàn)樗瞪嫌位蛑Я魉诩畢^(qū)，都是耕地且占比高、化肥施用強(qiáng)度較高和耗水較多的農(nóng)耕區(qū)。

圖5 潛在氮磷徑流質(zhì)量濃度空間分布Fig.5 Spatial distribution of the potential runoff concentrations of nitrogen and phosphorus

2.3 海河流域農(nóng)田面源污染關(guān)鍵源區(qū)識(shí)別

海河流域氮磷入河負(fù)荷的熱點(diǎn)區(qū)（圖6）主要分布在流域南部、石家莊和秦皇島的周邊地區(qū)。氮、磷熱點(diǎn)區(qū)分別占流域總面積的14.45%、16.02%，貢獻(xiàn)入河總負(fù)荷的23.39%、27.71%。

強(qiáng)調(diào)關(guān)鍵源區(qū)識(shí)別重要性的典型案例是Gburek等[31]得出流域20%區(qū)域貢獻(xiàn)了80%的磷流失量。本研究關(guān)鍵源區(qū)的貢獻(xiàn)率較低，主要原因是海河流域氮磷面源污染的空間異質(zhì)性不顯著，尤其是作為流域農(nóng)田主體的平原地區(qū)，種植模式、施肥量、地形條件、耕地占比都很相似。因此，基于污染源負(fù)荷識(shí)別關(guān)鍵源區(qū)的方法在海河流域提高農(nóng)田面源污染防控效率的作用不是很顯著。

河段潛在氮磷徑流質(zhì)量濃度按照自然斷點(diǎn)法劃分為5個(gè)等級(jí)，從低到高分別為低值區(qū)、中低區(qū)、中等區(qū)、中高區(qū)和高值區(qū)。其中潛在氮、磷徑流質(zhì)量濃度高值區(qū)分別為44、35個(gè)河段（圖6），分別占河段總數(shù)的17%、14%。

圖6 氮磷入河負(fù)荷熱點(diǎn)分析Fig.6 The hotspots analysis of export load of nitrogen and phosphorus

徒駭馬頰河是海河流域農(nóng)田面源污染最嚴(yán)重的水系，以其為例統(tǒng)計(jì)各縣域?qū)λ党隹谔幩|(zhì)的貢獻(xiàn)率（圖7）。該水系中上游各縣的貢獻(xiàn)率較高，在3%～6%之間，上游比中游略高，下游各縣貢獻(xiàn)率范圍為1%～3%，說明徒駭馬頰河水系上游是其主要污染源，需要重點(diǎn)治理。通過貢獻(xiàn)率可以迅速判別對(duì)流域出口貢獻(xiàn)程度最大的單元及貢獻(xiàn)程度，也可作為流域生態(tài)補(bǔ)償時(shí)的參考依據(jù)。

圖7 徒駭馬頰河各個(gè)縣域的氮磷貢獻(xiàn)率Fig.7 The contribution rate of nitrogen and phosphorus in each county of Tuhaimajia River

3 討 論

3.1 模擬結(jié)果的準(zhǔn)確性及對(duì)比

由于本研究針對(duì)的是整個(gè)海河流域農(nóng)田氮磷負(fù)荷的大尺度模擬，且NDR 模型估算的是總氮、總磷指標(biāo)，而我國生態(tài)環(huán)境部數(shù)據(jù)中心可公開獲取的水質(zhì)數(shù)據(jù)不包括總氮、總磷指標(biāo)，而且水體污染物除了農(nóng)田面源污染物，還包括點(diǎn)源污染物和其他面源污染物，因此很難全面系統(tǒng)地評(píng)價(jià)模型模擬結(jié)果的準(zhǔn)確度與可信性。

已有研究中，朱梅等[32]利用輸出系數(shù)模型估算海河流域2007年種植業(yè)面源污染負(fù)荷，得出地表徑流總氮、總磷分別為5.20 萬、0.80 萬t；馮愛萍等[33]利用DPeRS 模型估算的2016年海河流域農(nóng)田徑流型總氮和總磷排放量分別為12.7 萬t 和0.7 萬t。與已有研究結(jié)果對(duì)比發(fā)現(xiàn)，本文估算的農(nóng)田氮、磷入河總負(fù)荷與上述結(jié)果在數(shù)量級(jí)上是一致的，具體估算值有所差別，這可能與不同的模型應(yīng)用的數(shù)據(jù)源，研究時(shí)期以及模擬精度等有關(guān)。將潛在氮磷徑流質(zhì)量濃度的空間分布與2015年《中國環(huán)境狀況公報(bào)》海河流域水環(huán)境分布圖進(jìn)行對(duì)比，較為一致，說明本研究采用的面源污染估算方法具有一定的區(qū)域適用性和可行性。

3.2 農(nóng)田面源污染輸移引起的空間差異性

潛在氮磷徑流質(zhì)量濃度（圖5）與氮磷入河負(fù)荷（圖3）呈顯著的空間差異性。主要表現(xiàn)在：①西北部山區(qū)中氮磷入河負(fù)荷偏高的地區(qū)，其潛在氮磷徑流質(zhì)量濃度卻較低，主要是因?yàn)檫@些地區(qū)耕地占比偏低，大量林地和草地對(duì)農(nóng)田流失的氮磷污染物有攔截和稀釋作用。②海河流域中部和南部平原地帶的大清河、黑龍港運(yùn)東河、徒駭馬頰河的河段，潛在氮磷徑流質(zhì)量濃度較高，主要是因?yàn)檫@些地方種植模式以一年二熟為主，化肥施用強(qiáng)度大，且耗水嚴(yán)重導(dǎo)致徑流量偏低，兩者綜合作用下加重了水質(zhì)污染程度。潛在氮磷徑流質(zhì)量濃度是本地污染源和上游輸入源污染物匯集與徑流過程共同作用的結(jié)果。潛在氮磷徑流質(zhì)量濃度比氮磷入河負(fù)荷指標(biāo)更能準(zhǔn)確反映農(nóng)田氮磷面源污染對(duì)水環(huán)境影響的程度。

3.3 關(guān)鍵源區(qū)識(shí)別

關(guān)鍵源區(qū)識(shí)別一般針對(duì)面源污染負(fù)荷進(jìn)行[4,15-16,34]。本研究構(gòu)建了潛在氮磷徑流質(zhì)量濃度指標(biāo)，反映面源污染對(duì)水質(zhì)的影響。潛在氮磷徑流質(zhì)量濃度高值區(qū)（圖6 黑色圓圈內(nèi)）與入河負(fù)荷熱點(diǎn)區(qū)的空間分布格局不一致。因此，面源污染防控的關(guān)鍵源區(qū)識(shí)別應(yīng)將入河負(fù)荷和潛在氮磷徑流質(zhì)量濃度的空間分布結(jié)合進(jìn)行分析。對(duì)于潛在氮磷徑流質(zhì)量濃度高的河段，即使其集水區(qū)即源區(qū)位于入河負(fù)荷的非熱點(diǎn)區(qū)，其治理優(yōu)先序靠前。例如大清河水系東段支流的氮徑流質(zhì)量濃度很高，其源區(qū)以入河負(fù)荷次熱點(diǎn)和不顯著區(qū)為主，但面源污染防控時(shí)仍應(yīng)重點(diǎn)治理。在選擇優(yōu)先防控區(qū)域時(shí)，應(yīng)著重考慮對(duì)內(nèi)源污染為主的潛在氮磷徑流質(zhì)量濃度高的區(qū)域進(jìn)行源頭控制，而對(duì)外源污染為主的潛在氮磷徑流質(zhì)量濃度高的區(qū)域可以采用過程攔截和末端治理的方式，并追溯其上游的關(guān)鍵源區(qū)。

4 結(jié)論

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

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

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

4）農(nóng)田氮磷入河負(fù)荷關(guān)鍵源區(qū)主要分布在流域中部和南部、石家莊、秦皇島等周邊地區(qū)，需要重點(diǎn)防控。