基于SWAT模型的涪江流域下墊面對(duì)面源污染負(fù)荷的影響分析

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

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

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

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

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

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

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

選擇流域內(nèi)涪江上游、通河口、梓潼江、老池、大安及香山6個(gè)斷面對(duì)應(yīng)的匯流區(qū)，研究不同下墊面條件下面源污染負(fù)荷的構(gòu)成及其變化機(jī)理。涪江上游斷面為涪江江油市與平武縣的交界斷面，大安斷面為涪江的一級(jí)支流安居河四川—重慶交界斷面，老池?cái)嗝鏋楦⒔闪骶d陽(yáng)市和遂寧市交界斷面，香山斷面為四川省和重慶市交界斷面，梓潼江與通口河斷面為區(qū)內(nèi)典型支流出口斷面，6個(gè)典型匯流區(qū)基本情況如表1所示。

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

1.2 下墊面條件對(duì)面源污染負(fù)荷的影響分析

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

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

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

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

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

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

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

2 結(jié)果與分析

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

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

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

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

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

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

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

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

基于驗(yàn)證后的SWAT模型，將流域內(nèi)的點(diǎn)源排放量全部置零，并將18個(gè)氣象站的均值作為模型輸入，以消除點(diǎn)源對(duì)河道污染通量的貢獻(xiàn)以及流域降雨時(shí)空分布不均勻?qū)γ嬖次廴就康挠绊?。圖3為零點(diǎn)源通量和流域無(wú)差別降雨情況下的面源NH3、TP和IMn負(fù)荷（單位面積通量）比較，豐水期（6—9月）TP負(fù)荷在年總負(fù)荷中占44.4%～73.8%，平水期（4—5月、10月）占22.9%～34.2%，枯水期（11—2月）占5.7%～12.4%；豐水期IMn負(fù)荷占年負(fù)荷的74.4%～92.1%，平水期占18.6%～29.4%，枯水期占6.1%～12.4%，豐水期、平水期和枯水期NH3負(fù)荷分別占年總負(fù)荷的53.3%～71.4%、24.4%～32.4%和15.2%～22.4%。

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

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

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

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

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

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

2.3 下墊面對(duì)面源污染負(fù)荷影響敏感性分析

根據(jù)子流域內(nèi)的農(nóng)業(yè)區(qū)種植作物和化肥使用量的調(diào)查結(jié)果確定化肥中TN和TP施用量，根據(jù)農(nóng)村畜禽養(yǎng)殖數(shù)量以及農(nóng)村位置確定農(nóng)村生活、畜禽養(yǎng)殖污染入田量。按照下墊面實(shí)際情況，基于GIS信息提取土壤滲透系數(shù)、平均坡度、河網(wǎng)密度、土地利用類(lèi)型以及土壤本底信息。面源污染排放量、徑流條件和土壤本底共12個(gè)參數(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，說(shuō)明多因子協(xié)同條件下，徑流條件仍為影響面源污染負(fù)荷的最敏感因子，而污染排放量因子的影響最小。已有研究也表明，區(qū)域尺度下主控流動(dòng)途徑對(duì)污染物遷移和入河過(guò)程的影響程度高于其他因素[23]。

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

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

3 結(jié)論

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

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

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

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