柳 強，王 康，羅 彬，陳 鵬，楊 淵，陳玲玲（.四川省生態(tài)環(huán)境監(jiān)測總站，四川成都6009；.武漢大學水資源與水電工程科學國家重點實驗室，湖北武漢 007；.四川省地質工程勘察院，四川成都 6007；.內江市環(huán)境監(jiān)測中心站，四川內江 600）

流域下墊面是大氣與地面或水面的分界面[1]。農田化肥、畜禽養(yǎng)殖和農村生活污染排放物遷移、轉化最終進入河道的過程在很大程度上受下墊面徑流條件影響[2-4]。大量研究表明，由于各種面源污染源空間分布的差異性、污染排放的隨機性以及下墊面徑流過程的時空變異性，不同土地利用條件下徑流過程以及面源污染負荷表現(xiàn)出明顯差異性[5]和不確定性[6-8]。分布式水文和面源污染模型（如SWAT[9-10]、AGNPS[11]、CASC2D[12]和 DWSM[13-14]等）通常根據(jù)下墊面土壤性質、土地利用、地形等方面的差異，將匯流區(qū)劃分為一系列基本水文單元，采用不同的參數(shù)描述對應下墊面中的水文過程和污染物遷移、轉化、入河過程。然而，在流域尺度上難以獲取詳細的面源污染負荷及通量監(jiān)測數(shù)據(jù)，通常僅能夠通過若干控制斷面的水文和水質監(jiān)測數(shù)據(jù)進行參數(shù)率定和模型分析，模型中參數(shù)的敏感性在很大程度上受監(jiān)測點位置、監(jiān)測時間和長度、匯流區(qū)大小等多種因素的影響。了解下墊面對面源污染負荷的影響機制，對于分布式模型計算精度的提升以及面源污染監(jiān)控都具有重要的意義。

筆者以涪江流域為研究對象，基于SWAT模型模擬流域水文和污染負荷，采用SWAT-Cup進行參數(shù)率定和模型校驗，基于Sobol全局性方法[15]分析流域不同下墊面參數(shù)對面源污染負荷的影響，以期為流域面源污染防控決策提供依據(jù)。

1 面源污染負荷監(jiān)測及分析方法

1.1 基于SWAT模型的涪江流域面源污染通量過程模擬

涪 江 全長 670 km，流域（29.30°～33.05°N，103.73°～106.27°E）面積36 400 km2，位于三峽水庫上游，四川省綿陽市、遂寧市和重慶市潼南區(qū)、銅梁區(qū)境內。海拔介于200～5 500 m之間，流域多年平均氣溫17.8～18.1℃，多年平均年降水量800～1 400 mm。流域內水系、水文站、氣象站和水質自動監(jiān)測站位置如圖1（a）所示，流域內共有平武、冶城、江油、涪江橋、梓潼、羅江、三臺、潼川、射洪、南北堰和小河壩11個水文監(jiān)測站，逐日測定河道流量，時間序列為2007年1月—2016年12月。流域內及周邊有九寨溝、平武、安縣、北川、江油、中江、青川、劍閣、蒼溪、梓潼、三臺、鹽亭、西充、射洪、蓬溪、安岳、潼南和武勝18個氣象站，測定每日的日均溫度、最高溫度、最低溫度、日照時數(shù)、日均風速、日均濕度和日降雨量7個參數(shù)，時間序列為2006年1月—2016年12月。流域內香山、老池和大安3個斷面為四川省控水質斷面，設置有自動水質監(jiān)測站，測定氨氮、TP和高錳酸鹽指數(shù)（IMn）3種污染物濃度，每日進行6次測定，3個站的水質監(jiān)測時間分別為2012年1月—2016年12月、2014年1月—2016年12月和2006年1月—2016年12月。

流域30 m精度數(shù)字高程模型（DEM）如圖1（b）所示，土地利用分類如圖1（c）所示。流域內包括農田、農村生活、畜禽養(yǎng)殖和水土流失4類面源污染源，基于農村人口數(shù)量、農田類型及面積、畜禽種類及養(yǎng)殖數(shù)量、土壤類型、土地利用方式和地形數(shù)據(jù)，采用HJ 884—2018《污染源源強核算技術指南準則》核算面源污染物排放量。點源污染包括工業(yè)點源和城鎮(zhèn)污水處理廠2類，均采用四川省環(huán)境監(jiān)測總站環(huán)境統(tǒng)計數(shù)據(jù)。流域內點源分布和單位面積面源污染物排放量如圖1（d）所示，涪江流域土壤TN和TP含量分布分別如圖1（e）和圖1（f）所示。

選擇流域內涪江上游、通河口、梓潼江、老池、大安及香山6個斷面對應的匯流區(qū)，研究不同下墊面條件下面源污染負荷的構成及其變化機理。涪江上游斷面為涪江江油市與平武縣的交界斷面，大安斷面為涪江的一級支流安居河四川—重慶交界斷面，老池斷面為涪江干流綿陽市和遂寧市交界斷面，香山斷面為四川省和重慶市交界斷面，梓潼江與通口河斷面為區(qū)內典型支流出口斷面，6個典型匯流區(qū)基本情況如表1所示。

基于ArcGIS軟件將流域基本信息進行模型展布，確定SWAT模型的輸入信息后，以2005—2006年作為模型預熱期，消除初始條件對模擬結果的影響，采用2007—2014年的水文流量以及3種污染物的濃度監(jiān)測結果進行模型參數(shù)調整，采用SWAT-Cup進行參數(shù)率定。將2015—2016年作為模型驗證期，采用率定后的參數(shù)進行模型校驗。

1.2 下墊面條件對面源污染負荷的影響分析

下墊面對污染負荷的影響可以表現(xiàn)為3個方面：（1）面源污染物排放量，包括農田化肥TN和TP施用量、畜禽養(yǎng)殖和農村生活污染排放進入田間的污染物4個污染指標；（2）下墊面徑流條件，包括土壤性質、地表坡度、河網(wǎng)密度以及土地利用類型4個參數(shù)；（3）土壤本底，包括土壤氮磷本底含量、有機質含量以及C/N比4個指標。

圖1 涪江流域基本情況Fig.1 Basic information of the Fujiang River Basin

Sobol方法基于方差分解，通過將目標函數(shù)的總方差分解為單個參數(shù)的方差和參數(shù)間相互作用的方法，實現(xiàn)參數(shù)的全局性敏感性分析。對于所選擇的3類12個下墊面指標參數(shù)中的任一參數(shù)，根據(jù)其分布特性（流域內下墊面徑流參數(shù)和土壤本底參數(shù)空間分布基本符合對數(shù)正態(tài)分布，污染物排放量變化基本符合正態(tài)分布），采用Monte Carlo方法[16-17]，基于參數(shù)調查值的均值和方差，隨機產生2個不同的低差異隨機序列：

式（1）～（2）中，m為矩陣行數(shù)，即抽樣重復數(shù)；n為下墊面參數(shù)個數(shù)。MA和MB這2個隨機序列具有不同的統(tǒng)計特征，反映下墊面狀態(tài)的變化，利用矩陣MB中的第i列數(shù)據(jù)（下墊面參數(shù)xi的抽樣值）替換矩陣MA中的第i列數(shù)據(jù)（下墊面參數(shù)xi的抽樣值），形成一個新的聯(lián)合矩陣

對每一個下墊面參數(shù)都生成一個聯(lián)合矩陣，最終形成n+2個矩陣，基于SWAT模型對所有參數(shù)矩陣下的水文及面源污染通量過程進行計算。下墊面參數(shù)xi的一階敏感性指數(shù)總敏感性指數(shù)分別為

式（4）～（5）中，yA(k)、yB(k)和yCi(k)分別為矩陣 MA、MB、MCi第k行下墊面參數(shù)對應的模型輸出；V(yA)為矩陣MA輸出的方差。

由式（4）可以看出，參數(shù)的一階敏感性指數(shù)反映參數(shù)xi變化對流域徑流或污染負荷模擬結果變化方差所產生的貢獻，因而參數(shù)的一階敏感性指數(shù)越大，其敏感性越強。同樣，總敏感性指數(shù)反映該參數(shù)與其他參數(shù)之間交互作用對徑流或污染負荷模擬結果的敏感性，敏感性越強，則表明該參數(shù)與其他參數(shù)之間的交互作用對徑流或負荷的影響程度越大。

2 結果與分析

2.1 基于SWAT模型的涪江流域污染負荷模擬

SWAT模型參數(shù)率定分為2個步驟：首先根據(jù)下墊面性質的實驗結果以及相關性質確定參數(shù)的初始值；其次采用SWT-Cup確定參數(shù)的最終取值，使實測所有監(jiān)測數(shù)據(jù)與模擬數(shù)據(jù)之間的誤差達最低。SWAT-Cup 分析表明，CN2、ALPHA_BF、GW_DELAY、OV_N和CNCOEF等參數(shù)對徑流過程具有顯著影響（P＜0.05），CH_K對于徑流過程和污染物濃度均有明顯影響，USLE_P、CMN、PPERCO和ERORGP等參數(shù)對NH3、IMn和TP濃度的影響最為顯著（P＜0.05），涪江流域SWAT模型關鍵參數(shù)取值如表2所示。

涪江流域全部監(jiān)測斷面的流量和污染物濃度誤差標準均值如表3所示。模型參數(shù)率定期和模型校驗期涪江上游江油站、中游射洪站、下游小河壩站流量模擬和監(jiān)測結果的比較，以及大安站、香山站、老池站NH3、TP和IMn實測濃度和模擬濃度的比較見圖2。

采用 Nash-Sutcliffe系數(shù)[18]、相對均方根誤差RRMSE、相對偏差 FB和相對總誤差 FGE[19]4 個指標對流量以及NH3、TP和IMn濃度模擬誤差進行評價。SWAT應用于流域水文和污染負荷模擬的研究結果表明，流域尺度在102km2級別時，Nash-Sutcliffe系數(shù)普遍能夠達到 0.58～0.72；在 103km2級別時，能夠達到 0.51～0.55；在 104km2級別時，較 難超過 0.55[5，20-22]。

表3SWAT模擬誤差分析Table 3 Indexes to quantify the accuracy of SWAT simulations

圖2 流量、污染物濃度模擬值和實測值的比較Fig.2 Comparison of simulated and measured flow rate and NPS pollution concentrations

由表3可知，流量以及3種污染物濃度模擬結果的Nash-Sutcliffe系數(shù)普遍超過0.5，具有較好的模擬效果。FB用于描述計算值和模擬值之間的系統(tǒng)誤差，若IMn的模擬值FB/FGE在±0.3/0.5范圍內，則說明模擬值精度能夠被接受[22]。流量、NH3和TP的模擬值FB/FGE在±0.15/0.3范圍內，表明模擬精度良好。

2.2 不同下墊面的面源污染負荷特性分析

基于驗證后的SWAT模型，將流域內的點源排放量全部置零，并將18個氣象站的均值作為模型輸入，以消除點源對河道污染通量的貢獻以及流域降雨時空分布不均勻對面源污染通量的影響。圖3為零點源通量和流域無差別降雨情況下的面源NH3、TP和IMn負荷（單位面積通量）比較，豐水期（6—9月）TP負荷在年總負荷中占44.4%～73.8%，平水期（4—5月、10月）占22.9%～34.2%，枯水期（11—2月）占5.7%～12.4%；豐水期IMn負荷占年負荷的74.4%～92.1%，平水期占18.6%～29.4%，枯水期占6.1%～12.4%，豐水期、平水期和枯水期NH3負荷分別占年總負荷的53.3%～71.4%、24.4%～32.4%和15.2%～22.4%。

圖3 不同匯流區(qū)面源污染負荷比較Fig.3 Comparison of NPS pollution loads of the catchment

6個匯流區(qū)NH3、TP和IMn負荷的比較如表4所示。對比6個匯流區(qū)平均面源污染風險因子（表1）和面源污染負荷（表4）可知，涪江上游和通口河污染物排放量和下墊面徑流條件均明顯低于其他4個匯流區(qū)，大安站NH3、IMn和TP這3種污染物負荷最小，梓潼江斷面NH3負荷明顯超過其他斷面，而香山和通口河斷面的TP和IMn負荷值最大，且反映污染負荷變化幅度的標準差表現(xiàn)出隨污染負荷增加而增大的趨勢。相比其他匯流區(qū)，涪江上游不同水文期NH3變化幅度最大，而IMn變化幅度最小，相同徑流條件下，NH3和IMn負荷變化表現(xiàn)出明顯差異。這與IMn主要來源于農村生活和畜禽養(yǎng)殖，而NH3則來源于多種污染源有關。

表4 匯流區(qū)面源污染負荷Table 4 Non-point source pollution loads of the catchment kg·km-2

在下墊面以山丘區(qū)、自然林草為主的匯流區(qū)（如大安和涪江上游），下墊面徑流條件中的坡度和坡長因子對面源污染負荷的影響最大。IMn負荷主要來源于土壤中有機質以及分散分布在流域各處的畜禽養(yǎng)殖和農村生活污染，下墊面徑流條件中坡度和坡長對其的影響最為大。比較表1和表4可以看出，盡管老池匯流區(qū)的污染物排放量因子和下墊面徑流條件均高于香山匯流區(qū)，然而老池斷面匯流區(qū)內NH3、TP和IMn這3種污染負荷及其變化幅度均在一定程度上小于香山匯流區(qū)，老池匯流區(qū)內農田（特別是水田區(qū)）面積比例明顯增加，由于平原區(qū)農田對徑流過程具有較強的滯留能力，很大程度上減少了徑流入河量，因此降低了匯流區(qū)內的污染負荷，其他匯流區(qū)的情況基本相同。

在6個匯流區(qū)中，以農田為主的大安匯流區(qū)3種污染物負荷均為最小。此外，在以農田為主的匯流區(qū)，TP以懸移質攜帶入河為主要形式，由于農田對徑流具有一定的調節(jié)能力，因此不同水文期TP負荷變化最?。▓D3）。

涪江上游和通口河下墊面類型和面積占比較為一致，然而由于面源污染排放量因子差別明顯，3種污染物負荷存在較大差異，面源污染排放量因子高的污染負荷變幅均顯著大于面源污染排放量因子低的變幅。相比污染排放量因子，下墊面徑流條件對于通量的影響更為明顯。

2.3 下墊面對面源污染負荷影響敏感性分析

根據(jù)子流域內的農業(yè)區(qū)種植作物和化肥使用量的調查結果確定化肥中TN和TP施用量，根據(jù)農村畜禽養(yǎng)殖數(shù)量以及農村位置確定農村生活、畜禽養(yǎng)殖污染入田量。按照下墊面實際情況，基于GIS信息提取土壤滲透系數(shù)、平均坡度、河網(wǎng)密度、土地利用類型以及土壤本底信息。面源污染排放量、徑流條件和土壤本底共12個參數(shù)的Sobol一階和總敏感性指數(shù)如表5所示。

NH3、TP和IMn的一階敏感性指數(shù)均值分別為0.076、0.054和0.119，其中徑流條件所占比例分別為71.1%、88.8%和55.6%。單因子的敏感程度比較而言，徑流條件的影響最大。面源污染物排放量、徑流條件和土壤本底的總敏感性指數(shù)均值分別為0.048、0.161和0.077，說明多因子協(xié)同條件下，徑流條件仍為影響面源污染負荷的最敏感因子，而污染排放量因子的影響最小。已有研究也表明，區(qū)域尺度下主控流動途徑對污染物遷移和入河過程的影響程度高于其他因素[23]。

NH3、TP和IMn這3種污染物的總敏感性指數(shù)平均值分別為0.113、0.072和0.140，可以看出，TP對于下墊面的敏感性最弱，而IMn對下墊面的敏感性最強。這一結果表明，涪江流域TP來源于農田以及其他下墊面表層土壤，以懸移質攜帶的形式進入河道中，由于其來源和入河途徑較為單一，因此在3種污染物中敏感性最弱，且徑流條件的單因素敏感性和總敏感性在全部因子敏感性中占比均為最大。由于IMn來源高度分散，加上不同來源中徑流入河條件存在差異，因此其敏感性最強，但其徑流條件單因素敏感性和總敏感性在全部因子敏感性中占比均為最小。一些研究成果表明，流域的不同污染物通量過程之間及其與徑流過程之間均存在明顯差異，造成這種差異的主要原因在于下墊面徑流過程與污染負荷之間的非線性響應過程[24-25]。筆者對各種因子的敏感程度進行進一步分析，NH3、TP和IMn污染物分別有5、3和6個因子的總敏感性指數(shù)超過0.1，達到顯著水平[26]，多因子作用是造成污染負荷與徑流過程非線性響應的重要原因。

表5 下墊面對面源污染負荷影響的敏感性分析Table 5 Sensitivity analysis of the underlying surface to the NPS pollution loads

3 結論

（1）基于SWAT模型模擬了涪江流域徑流過程和面源污染負荷，采用SWAT-Cup進行了參數(shù)率定。采用Nash-Sutcliffe系數(shù)、相對均方根誤差、相對偏差和相對總誤差4種指標進行檢驗，結果表明模擬徑流過程和污染負荷與實測值符合良好。

（2）各種下墊面條件下，NH3、TP和IMn負荷均隨著負荷均值的增加而增大，在以山丘區(qū)、自然林草為主的匯流區(qū)下墊面條件下，徑流條件因子中的坡度和坡長對于面源污染負荷的影響最為明顯，由于農田對于徑流過程具有一定的調節(jié)能力，面源污染負荷隨匯流區(qū)內農田面積比例增加而降低。

（3）采用Sobol方法分析下墊面中污染物排放量、徑流條件和土壤本底對面源污染負荷的影響，結果表明3類因子中，徑流條件對面源污染負荷的影響最為敏感，而排放量的影響最弱。對于3種污染物，下墊面對TP負荷的敏感性最弱，對IMn負荷的敏感性最強。

（4）涪江流域面源污染分析表明，徑流條件對于河道污染負荷的影響最大，控制下墊面徑流條件對于降低面源污染負荷具有積極作用。