多情景分析的農(nóng)業(yè)面源污染關(guān)鍵源區(qū)識別軟件開發(fā)及應用

覃苑，胡海棠，淮賀舉，李存軍*，張巧玲，楊鐵利，王佳宇

1.遼寧科技大學

2.北京市農(nóng)林科學院信息技術(shù)研究中心

農(nóng)業(yè)面源污染（agricultural non-point source pollution，ANSP）指由于農(nóng)業(yè)生產(chǎn)中的氮磷營養(yǎng)元素、農(nóng)藥及其他溶解態(tài)和固態(tài)污染物流失導致的污染。近年來，農(nóng)業(yè)面源污染對水體污染的貢獻越來越受到重視。2021年3月生態(tài)環(huán)境部和農(nóng)業(yè)農(nóng)村部聯(lián)合發(fā)布了《農(nóng)業(yè)面源污染治理與監(jiān)督指導實施方案（試行）》，明確了“十四五”至2035年農(nóng)業(yè)面源污染防治的總體要求、工作目標，提出主要任務包括深入推進農(nóng)業(yè)面源污染防治、完善農(nóng)業(yè)面源污染防治政策機制和加強農(nóng)業(yè)面源污染治理監(jiān)督管理。ANSP的污染負荷估算、關(guān)鍵源區(qū)源頭識別和治理決策是監(jiān)督管理的重要內(nèi)容。

ANSP在空間上分布分散且隱蔽，其形成受到區(qū)域的氣候、地形、土地利用類型和植被覆蓋等多方面因素的影響，導致ANSP在監(jiān)測、模擬和控制等方面面臨較大挑戰(zhàn)。GIS技術(shù)在空間數(shù)據(jù)管理、集成、分析方面有優(yōu)勢，已被廣泛應用于ANSP的數(shù)據(jù)管理、空間分析、面源污染模型集成、信息系統(tǒng)開發(fā)和可視化等方面[1]。國內(nèi)外采用年度農(nóng)業(yè)面源污染模型（annualized agricultural nonpoint source pollution model，AnnAGNPS），水土評價模型（soil and water assessment tool，SWAT），輸出系數(shù)模型（export coefficient model，ECM）等模型和平均濃度法等方法進行ANSP模擬與識別研究。如夏昊[2]結(jié)合AnnAGNPS模型，開發(fā)了基于GIS的面源信息系統(tǒng)，其擁有空間分析和可視化表達等功能，用于管理ANSP中不同來源的空間和屬性數(shù)據(jù)；喬衛(wèi)芳等[3]基于研究區(qū)域數(shù)字高程模型（digital elevation model，DEM）數(shù)據(jù)，將GIS與SWAT模型結(jié)合構(gòu)建了ANSP基礎(chǔ)信息庫，分析了研究區(qū)ANSP負荷的分布特征，并模擬了區(qū)域內(nèi)主要土地利用類型單位面積ANSP的負荷量。但AnnAGNPS和SWAT等模型所需參數(shù)多，計算復雜，在大區(qū)域的模擬上略顯乏力；生態(tài)系統(tǒng)服務和權(quán)衡的綜合評估（InVEST）模型中的營養(yǎng)物傳輸率模型（NDR）所需參數(shù)少，機理清晰，適用于大區(qū)域的ANSP模擬研究[4]。

在進行ANSP治理時，關(guān)鍵源區(qū)識別是面源污染管理的重要組成部分。研究表明，少數(shù)區(qū)域輸出的污染物通常是整個流域污染負荷的主要貢獻者，對水體環(huán)境質(zhì)量有著至關(guān)重要的作用[5]?！?015年中國環(huán)境狀況公報》顯示，海河、巢湖等流域的水質(zhì)狀況存在明顯的空間異質(zhì)性，表現(xiàn)為支流水質(zhì)在重度污染的情況下，其所屬的干流僅為輕度污染[6]。目前主要是針對面源污染負荷空間分布的異質(zhì)性特征進行ANSP關(guān)鍵源區(qū)識別[7]，但是單純考慮污染負荷易忽略流域面積與產(chǎn)水量對污染負荷貢獻的影響，導致關(guān)鍵源區(qū)的識別準確度下降，進而影響治理效果。應結(jié)合潛在氮、磷徑流濃度和入河負荷的空間分布進行ANSP防控的關(guān)鍵源區(qū)識別與分析[8]。目前農(nóng)業(yè)面源污染關(guān)鍵源區(qū)識別涉及水文分析、農(nóng)業(yè)面源污染模型、柵格矢量運算和統(tǒng)計，對于面源污染非專業(yè)人士和基層工作者存在跨平臺操作、數(shù)據(jù)類型多、資料匯總繁瑣、計算量大等問題，亟需集成式的軟件工具平臺。

ANSP包括農(nóng)田種植業(yè)面源污染、畜禽養(yǎng)殖、農(nóng)村生活污染等，其中農(nóng)田種植業(yè)面源污染是農(nóng)業(yè)面源污染的重要組成部分，污染源相對分散且影響因素多。筆者針對在農(nóng)田種植業(yè)氮、磷化肥施用形成的面源污染關(guān)鍵源區(qū)識別過程中情景單一、數(shù)據(jù)量大、涉及環(huán)節(jié)多，存在計算復雜、操作繁瑣等問題，將GIS與InVEST模型結(jié)合，研發(fā)農(nóng)業(yè)面源污染關(guān)鍵源區(qū)識別與治理模擬的集成計算機軟件系統(tǒng)，在多種情景下對關(guān)鍵源區(qū)進行識別，并模擬治理后的效果，以期為農(nóng)田種植業(yè)面源污染治理提供有效的輔助手段。

1 ANSP關(guān)鍵源區(qū)識別軟件系統(tǒng)設(shè)計

基于GIS二次開發(fā)，集成水文網(wǎng)絡(luò)拓撲構(gòu)建和面源污染關(guān)鍵源區(qū)識別技術(shù)，利用DEM數(shù)據(jù)提取河流節(jié)點、河段與子流域等水文地貌要素，整合InVEST產(chǎn)水量模型和營養(yǎng)物傳輸率模型，結(jié)合基于入河污染負荷、潛在徑流污染物濃度和負荷與產(chǎn)水量比值3種情景的面源污染關(guān)鍵源區(qū)識別方法，以決策支持為重點，建立農(nóng)業(yè)面污染關(guān)鍵源區(qū)識別與治理模擬軟件系統(tǒng)。

1.1 系統(tǒng)目標

通過研究ANSP的污染負荷情況與成因，運用GIS技術(shù)，研制具有農(nóng)業(yè)面源污染關(guān)鍵源區(qū)識別與治理模擬的輔助決策軟件系統(tǒng)。對研究區(qū)域內(nèi)的水文網(wǎng)絡(luò)進行拓撲構(gòu)建、分析和專題信息可視化表達，快速高效得到研究區(qū)內(nèi)拓撲關(guān)系、子流域營養(yǎng)元素分布規(guī)律和子流域潛在營養(yǎng)元素徑流濃度等信息。在此基礎(chǔ)上，根據(jù)不同的方法分析面源污染關(guān)鍵源區(qū)空間位置，并模擬關(guān)鍵源區(qū)治理后的效果，為ANSP治理提供依據(jù)。

1.2 功能設(shè)計

以靈活的應用集成接口構(gòu)建GIS軟件系統(tǒng)，對研究區(qū)內(nèi)的水網(wǎng)數(shù)據(jù)、DEM數(shù)據(jù)和土壤數(shù)據(jù)等多源異構(gòu)的資源環(huán)境數(shù)據(jù)進行集成與調(diào)度，根據(jù)研究區(qū)已知的地理數(shù)據(jù)和資源環(huán)境數(shù)據(jù)，識別面源污染關(guān)鍵源區(qū)。系統(tǒng)按照實際功能分為6個模塊單元：1)拓撲關(guān)系分析單元，用于基于目標區(qū)域的DEM數(shù)據(jù)構(gòu)建水文網(wǎng)絡(luò)拓撲關(guān)系；2)產(chǎn)水量分析單元，用于利用InVEST模型的產(chǎn)水量子模型，根據(jù)水文網(wǎng)絡(luò)拓撲關(guān)系，確定子流域產(chǎn)水量；3)入河負荷估算單元，根據(jù)地表營養(yǎng)物負荷、傳輸率、次表層營養(yǎng)物傳輸率和各柵格單元的入河負荷，計算子流域營養(yǎng)物入河總負荷；4)潛在徑流濃度分析單元，用于根據(jù)子流域產(chǎn)水量，確定目標區(qū)域內(nèi)各子流域的潛在氮、磷徑流濃度；5)關(guān)鍵源區(qū)識別單元，基于氮、磷入河負荷，潛在氮、磷徑流濃度，子流域氮、磷負荷與產(chǎn)水量比值3種情景，識別和確定目標區(qū)域內(nèi)的面源污染關(guān)鍵源區(qū)；6)關(guān)鍵源區(qū)治理模擬單元，用于模擬對識別出的關(guān)鍵源區(qū)進行治理后的結(jié)果，確定治理效果。

ANSP關(guān)鍵源區(qū)識別與治理模擬統(tǒng)一建模語言如圖1所示。由圖1可知，1種關(guān)鍵源區(qū)識別結(jié)果對應1個關(guān)鍵源區(qū)識別模型，模擬治理結(jié)果依賴于關(guān)鍵源區(qū)識別結(jié)果；1個關(guān)鍵源區(qū)識別模型對應若干農(nóng)田氮、磷入河負荷估算模型，子流域負荷與產(chǎn)水量比值估算模型與潛在氮、磷徑流濃度估算模型；農(nóng)田氮、磷入河負荷估算模型包含1個氮、磷流失負荷分配結(jié)果，1組地表營養(yǎng)物傳輸率和1組次地表營養(yǎng)物傳輸率；子流域負荷與產(chǎn)水量比值估算模型包含1個產(chǎn)水量模型，依賴于農(nóng)田氮、磷入河負荷估算模型結(jié)果；潛在氮、磷徑流濃度估算模型包含1個產(chǎn)水量模型，依賴于水網(wǎng)拓撲關(guān)系和農(nóng)田氮、磷入河負荷估算模型結(jié)果。ANSP關(guān)鍵源區(qū)識別與治理模擬軟件結(jié)構(gòu)如圖2所示。

圖1 ANSP關(guān)鍵源區(qū)識別與治理模擬統(tǒng)一建模語言Fig.1 Unified modelling language for the identification and control simulation of ANSP key source areas

圖2 ANSP關(guān)鍵源區(qū)識別與治理模擬軟件結(jié)構(gòu)Fig.2 Software structure diagram for identification and simulation of the treatment of ANSP key source areas

2 ANSP關(guān)鍵源區(qū)識別主要方法與模型

2.1 水文網(wǎng)絡(luò)拓撲構(gòu)建

獲取研究區(qū)的DEM數(shù)據(jù)，基于DEM數(shù)據(jù)建立水文網(wǎng)絡(luò)拓撲關(guān)系。利用計算機編程提取河段以及子流域，徑流節(jié)點是基于幾何網(wǎng)絡(luò)的方法提取的，分為源頭節(jié)點、匯流節(jié)點、流域出口節(jié)點[9]。具體流程見圖3。

圖3 水文地貌要素提取Fig.3 Extraction of hydrological and geomorphic elements

流域水文網(wǎng)絡(luò)拓撲結(jié)構(gòu)的建立以流域匯流關(guān)系建立的基礎(chǔ)——河段為紐帶，對河網(wǎng)的徑流節(jié)點、河段、子流域進行矢量要素的空間分析，構(gòu)建拓撲關(guān)系。河段屬性表里記錄對應河段的起始節(jié)點和終止節(jié)點，以及對應的子流域編碼；子流域?qū)傩员砝镉涗浥c其相對應的河段編碼。定義節(jié)點的編碼為Point-ID，節(jié)點類型是NodeType，上下游節(jié)點分別為FROM_NODE、TO_NODE；河段編碼為 River_ID，長度為River_len，匯入和匯出河段分別為Riv_in、Riv_out；子流域編碼為Basin_ID，每個子流域的匯入和匯出流域編碼分別為Basin_in、Basin_out。

2.2 流域氮、磷入河負荷計算方法

流域氮、磷入河負荷的計算分為2步，首先計算出流域氮、磷流失負荷，再將其與InVEST模型中的NDR模型進行結(jié)合，計算得到流域氮、磷入河負荷。

2.2.1 氮、磷流失負荷模型

流域氮、磷流失負荷計算公式如下：

2.2.2 氮、磷入河負荷模型

通過模擬區(qū)域內(nèi)營養(yǎng)物的來源和轉(zhuǎn)移過程，計算氮、磷入河系數(shù)，在考慮氮、磷流失負荷的情況下估算氮、磷入河負荷。首先按照式（2）將氮、磷流失負荷分配到表層和次表層：

最后，匯總計算氮、磷入河負荷，計算公式如下：

2.3 流域產(chǎn)水量計算方法

利用InVEST模型的產(chǎn)水量子模型計算流域產(chǎn)水量，計算公式如下：

為了反映流域多年平均產(chǎn)水量，降低年間的氣候差異和降水量差異造成不確定性的影響，潛在的蒸發(fā)散失量與降水量采用2001—2015年的平均值。

2.4 ANSP關(guān)鍵源區(qū)的識別

由于面源污染易受到地形、植被、氣候、施肥和土地利用類型等因素影響[12-15]，不同識別方法通常會呈現(xiàn)出不同程度的空間差異性。本研究采用3種不同情景識別關(guān)鍵源區(qū)。

2.4.1 基于入河負荷

該情景由分區(qū)統(tǒng)計獲得子流域氮、磷入河負荷，根據(jù)每個子流域氮、磷入河負荷大小進行關(guān)鍵源區(qū)的識別[16-19]，計算公式如下：

2.4.2 基于潛在徑流濃度

該情景通過計算河段氮、磷入河通量與徑流量的比值，得到河段潛在徑流濃度，根據(jù)潛在徑流濃度的貢獻率大小進行關(guān)鍵源區(qū)的識別[8]。

2.4.3 基于負荷與產(chǎn)水量比值

將子流域產(chǎn)水量與氮、磷入河負荷同時納入考慮，計算子流域產(chǎn)水量與氮、磷入河負荷的比值，根據(jù)該比值進行關(guān)鍵源區(qū)的識別[20-22]，計算公式如下：

2.4.4 關(guān)鍵源區(qū)識別

根據(jù)式（8）～式（11）計算得到入河負荷、潛在徑流濃度和負荷與產(chǎn)水量比值，對子流域或河段對應子流域進行降序排列。按順序累加子流域入河負荷，每累加1次就通過子流域累計入河負荷與流域總負荷之比（Per）來判斷是否滿足設(shè)定閾值，若大于等于設(shè)定閾值則停止累加并將這些區(qū)域標記為面源污染關(guān)鍵源區(qū)，流程如圖4所示。

圖4 關(guān)鍵源區(qū)識別流程Fig.4 Identification process for key source area

2.5 關(guān)鍵源區(qū)治理模擬

針對3種情境識別出的關(guān)鍵源區(qū)，設(shè)定相同的預期治理幅度，可直觀地觀察到在不同情景下對關(guān)鍵源區(qū)治理的效果。將3種情景下識別出的關(guān)鍵源區(qū)，按照與河流源頭的接近程度進行排序，根據(jù)設(shè)定治理幅度依次降低這些區(qū)域的氮或磷負荷，再從上游—下游的順序依次更新流域所有受影響的河段的氮或磷潛在徑流濃度。

3 農(nóng)業(yè)面源污染關(guān)鍵源區(qū)識別軟件應用

3.1 研究區(qū)概況及特征

海河流域（112°E～120°E，30°N～43°N）地處我國華北地區(qū)，包括北京、天津、河南北部、山東北部、山西東部、內(nèi)蒙古和遼寧小部分區(qū)域與河北省大部分區(qū)域。流域地勢為西北高、東南低，包含灤河、北三河、大清河、永定河、子牙河、漳衛(wèi)南運河、黑龍港運東河、徒駭馬頰河與冀東沿海八大水系[23]。海河流域是我國三大糧食生產(chǎn)基地之一，流域土地利用類型主要為耕地，占流域總面積的48.9%。ANSP是海河流域主要污染源之一，也是海河水質(zhì)污染的主要原因之一[24]。

3.2 總氮關(guān)鍵源區(qū)識別

研發(fā)的農(nóng)業(yè)面源污染關(guān)鍵源區(qū)識別軟件主界面見圖5。以TN為例，使用該軟件對海河流域的農(nóng)業(yè)面源污染關(guān)鍵源區(qū)進行分析。流域ANSP的面積與TN、TP的關(guān)系通常是呈對數(shù)上升趨勢，往往不到10%的面積就貢獻了30%以上的氮、磷污染負荷[25]，據(jù)此將軟件識別幅度設(shè)定為30%，分別基于入河污染負荷、潛在徑流污染物濃度、負荷與產(chǎn)水量比值3種情景進行TN面污染關(guān)鍵源區(qū)識別，結(jié)果如圖6、圖7所示。

圖5 農(nóng)業(yè)面源污染關(guān)鍵源區(qū)識別軟件用戶界面Fig.5 User interface of the key source area identification software for ANSP

由圖6、圖7可知，TN入河負荷較高的區(qū)域主要分布在灤河、北三河、璋衛(wèi)河水系與黑龍港運東河下游，子牙河、大清河水系中部以及徒駭馬頰河全域，呈現(xiàn)自西北向東南方向增長的分布特征。潛在TN徑流濃度升高的區(qū)域主要分布中部—東南部地區(qū)的黑龍港運東河、徒駭馬頰河與大清河水系。TN負荷與產(chǎn)水量比值高的區(qū)域集中分布在中部—南部地區(qū)的徒駭馬頰河上游、大清河水系及黑龍港運東河?；谌牒迂摵汕榫跋伦R別出的關(guān)鍵源區(qū)較為分散，除了永定河與子牙河水系外皆有分布?；跐撛趶搅鳚舛惹榫跋碌年P(guān)鍵源區(qū)分布較為集中，多為中部—南部分布，其余在東南部也有分布?；诹饔蜇摵僧a(chǎn)水比值情景下識別出的關(guān)鍵源區(qū)分布非常集中，以中部—南部分布為主，識別出的區(qū)域與基于潛在徑流濃度情景下的區(qū)域重合度較高。

圖6 3種情景下的TN污染區(qū)域分布情況Fig.6 Regional distribution of total nitrogen pollution under three scenarios

圖7 3種情景下識別的關(guān)鍵源區(qū)Fig.7 Key source areas identified by three scenarios

3.3 治理模擬結(jié)果與分析

在識別出關(guān)鍵源區(qū)的基礎(chǔ)上，對關(guān)鍵源區(qū)進行模擬治理。為了更好地反映治理效果，將識別出的關(guān)鍵源區(qū)TN、TP入河負荷清零，并在考慮水網(wǎng)的拓撲關(guān)系下，對河段TN、TP潛在徑流濃度進行重新統(tǒng)計。模擬結(jié)果如圖8所示。由圖8可知，基于入河負荷情景下，模擬治理的效果欠佳，流域中部水系的潛在徑流濃度下降不明顯；基于潛在徑流濃度、負荷與產(chǎn)水量比值2種情景下，治理效果相近，治理區(qū)域的潛在徑流濃度均有較為明顯的下降，但對東南部和南部的部分水系治理效果不明顯。

圖8 原始徑流濃度與3種情景下治理效果對比Fig.8 Comparison chart of original runoff concentration and that of three methods after treatment

計算了3種情景識別出的面污染關(guān)鍵源區(qū)模擬治理后TN、TP負荷的下降率，結(jié)果見表1。

由表1可知，在研究區(qū)的數(shù)據(jù)與條件下，基于入河負荷情景識別出的面污染關(guān)鍵源區(qū)的TN治理效果弱于基于潛在徑流濃度和基于負荷產(chǎn)水量比值情景的治理效果。由于海河流域西北地區(qū)入河負荷高，但林地和草地面積較大，攔截和稀釋了大量的氮、磷污染物，而中南部地區(qū)耕地面積占比高，且采取了一年兩熟的種植模式，化肥的施用量大且耗水嚴重，導致地表徑流量偏低，所以3種情景下識別出的關(guān)鍵源區(qū)有著較顯著的空間差異性。由于基于潛在徑流濃度情景反映的是河段與其上游輸入污染物的匯集和徑流過程共同作用的結(jié)果，而子流域負荷與產(chǎn)水量比值反映子流域產(chǎn)水量與氮、磷污染負荷比值，二者更能反映出ANSP對水質(zhì)的影響程度。

表1 3種情景下關(guān)鍵源區(qū)模擬治理效果Table 1 Simulation of governance effects in key source areas under three scenarios

對海河流域的治理模擬結(jié)果表明，僅針對面源污染入河負荷識別出的氮和磷污染關(guān)鍵源區(qū)在治理的結(jié)果上并不理想，基于潛在徑流濃度和基于負荷與產(chǎn)水量比值情景下的治理效果優(yōu)于前者，且結(jié)果接近。這是因為面源污染受到地形、植被和土地利用類型等因素影響較大，單一情景下識別效果并不一定能準確反映面源污染對水質(zhì)的影響程度?？紤]到不同地區(qū)地理環(huán)境的復雜性，單一識別方法具有一定的局限性，在實際使用中，ANSP防控關(guān)鍵源區(qū)的識別應將入河負荷、潛在徑流濃度和負荷產(chǎn)水量比值的空間分布情況結(jié)合分析，以期得到更好的治理效果，支持農(nóng)業(yè)面源污染治理科學決策。

4 結(jié)論

(1)針對農(nóng)業(yè)面源污染關(guān)鍵源區(qū)識別過程中情景單一、數(shù)據(jù)量大、涉及環(huán)節(jié)多，以及模型復雜、操作繁瑣和治理效果不明確等問題，設(shè)計開發(fā)了多情景分析的農(nóng)業(yè)面源污染關(guān)鍵源區(qū)識別軟件。軟件具備拓撲關(guān)系分析、產(chǎn)水量分析、入河負荷估算、徑流濃度分析、關(guān)鍵源區(qū)識別和關(guān)鍵源區(qū)模擬治理等功能，可模擬3種常見情景下農(nóng)業(yè)面源污染分布的特點，為農(nóng)業(yè)面源污染治理決策提供有效的輔助手段。

(2)以海河流域為例，進行了農(nóng)業(yè)面源污染關(guān)鍵源區(qū)的識別和治理模擬應用，關(guān)鍵源區(qū)識別結(jié)果顯示，流域TN負荷呈現(xiàn)自西北向東南方向增長的分布特征，中部—東南部地區(qū)的潛在TN徑流濃度較高，中部—南部地區(qū)的TN負荷與產(chǎn)水量比值明顯高于其他地區(qū)。模擬治理的結(jié)果表明，基于入河負荷情景對氮或磷的面源污染關(guān)鍵源區(qū)治理效果較不理想，而基于潛在徑流濃度、負荷與產(chǎn)水量比值情景的治理效果較好，且后者略優(yōu)于前者。