基于SWAT模型的龍溪河瀘縣境內(nèi)面源污染特征分析

馬亞麗，白祖暉，敖天其

(1.甘肅農(nóng)業(yè)大學(xué)水利水電學(xué)院，蘭州 730000；2.甘肅省水利水電勘測設(shè)計研究院有限責(zé)任公司，蘭州 730000；3.四川大學(xué) 水力學(xué)與山區(qū)河流開發(fā)保護國家重點實驗室，成都 610065)

相關(guān)研究表明，SWAT模型在河灌區(qū)、流域面積僅為6.26 km2的重慶市萬州區(qū)陳家溝小流域、柴河水庫流域、三峽庫區(qū)等多個流域[7-10]的面源污染負(fù)荷模擬研究都可以取得較好效果，因此，本文首次利用SWAT模型在龍溪河瀘縣境內(nèi)的面源污染進(jìn)行研究，探討SWAT模型在該流域的適用性，并建立適合該流域的SWAT模型數(shù)據(jù)庫及確定敏感性參數(shù)值，定量計算該流域的氮磷污染負(fù)荷，分析氮磷污染負(fù)荷特征。其成果將對該流域的面源污染的控制和管理、水源地的保護、城市飲水安全以及流域水資源管理有一定的理論和實踐意義，可以作為龍溪河流域水資源保護、規(guī)劃和面源污染控制、規(guī)劃、決策的依據(jù)。

1 流域概況

龍溪河發(fā)源于四川省登東山，屬于長江左岸一級支流，總干流長110 km，流域面積502 km2，途徑瀘縣立石、云錦、玄灘、奇峰、兆雅、云龍等鎮(zhèn)以及龍馬潭區(qū)長安、石洞、魚塘、羅漢等鎮(zhèn)，在龍馬潭區(qū)匯入長江，其中，龍溪河瀘縣境內(nèi)干流長約69 km，流域面積242 km2，涉及人口14.60萬人。多年平均流量為6.58 m3/s，平均深度為0.4 m，河寬40～60 m，河道地形起伏不大，平均比降為0.17%。研究區(qū)四季分明，具有亞熱帶濕潤氣候的特點。全縣多年平均降水量為1 013.6 mm，降雨量最大的年份可達(dá)1 450.2 mm；多年平均風(fēng)速1.2 m/s，一般以偏北風(fēng)為主；多陰天，日照少；多年平均溫度為17.7 ℃；海拔在316～376 m左右，地形起伏不大。

2 數(shù)據(jù)來源及處理

SWAT模型數(shù)據(jù)庫包括空間數(shù)據(jù)庫和屬性數(shù)據(jù)庫，空間數(shù)據(jù)包括高程地形圖(DEM)、土地利用圖、土壤類型圖、水系圖，屬性數(shù)據(jù)包括土壤、氣象、水文、水質(zhì)、污染源、農(nóng)業(yè)管理方式等，其中空間數(shù)據(jù)必須具有統(tǒng)一的投影坐標(biāo)[11]，本研究采用Albers等積圓錐投影，相應(yīng)投影參數(shù)見表1。由于瀘縣環(huán)保局委托編制《瀨溪河流域(瀘縣境內(nèi))水污染評估項目成果報告》、《瀘州市龍溪河生態(tài)環(huán)境瀘縣段(2013-2015年)保護方案》，由瀘縣環(huán)保局協(xié)調(diào)收集獲得瀘縣相關(guān)資料。各種數(shù)據(jù)的精度、格式和來源具體見表2。

表1 研究區(qū)投影參數(shù)表Tab.1 The projection parameters of the study area

表2 模型所需數(shù)據(jù)分辨率、格式和來源Tab.2 The resolution, format, and source of the data required for the model

2.1 流域內(nèi)子流域劃分

SWAT模型基于DEM圖以及已有的水系圖生成最終水系以及劃分子流域。利用SWAT模型進(jìn)行流域離散化，首先根據(jù)DEM生成的河網(wǎng)分布和出水口的位置，將整個流域劃分為若干個子流域，其中，最小河道集水面積閾值設(shè)置得越小，生成的水系就越詳細(xì)，子流域劃分?jǐn)?shù)目就越多，然后根據(jù)每個子流域內(nèi)的土地利用、土壤和坡度屬性疊加劃分出一個或多個水文響應(yīng)單元。HRU是SWAT模型模擬的基本單位，假定認(rèn)為具有統(tǒng)一的水文效應(yīng)。本研究中設(shè)定最小河道集水面積閾值400 hm2，以龍溪河在瀘縣的出口斷面作為流域斷面劃分17個子流域，具體子流域劃分見圖1。

圖1 子流域劃分圖Fig.1 The sub-watershed division diagram

2.2 各子流域中水文響應(yīng)單元(HRUs)分配

根據(jù)土地利用類型、土壤類型和坡度的分類，將每個子流域劃分為一個或多個水文響應(yīng)單元(HRUs)，HRU作為SWAT模型模擬的基本單位，計算每個HRU的徑流量、輸沙量和營養(yǎng)物質(zhì)污染負(fù)荷。為適應(yīng)SWAT模型非點源模擬中土地利用類型的分類，根據(jù)流域內(nèi)原來的分類，土地利用類型重新分類編碼[12]，土壤類型圖與土地利用類型圖的處理方法相似。將重分類的土地利用數(shù)據(jù)、土壤數(shù)據(jù)以及坡度分類進(jìn)行疊加，本研究中設(shè)置土地利用、土壤和坡度的閾值分別為20%、25%、15%，最終將流域劃分為65個水文響應(yīng)單元。

2.3 模型參數(shù)的寫入及編輯

建立各類氣象信息測站位置表和實測數(shù)據(jù)表，寫入各類氣象數(shù)據(jù)，創(chuàng)建模型所需的各種輸入文件，包括：結(jié)構(gòu)文件(.fig)，土壤文件(.soil)、氣象文件(.wgn)，子流域文件(.sub)、水文響應(yīng)單元文件(.hru)、主河道文件(.rte)、地下水文件(.gw)、水利用文件(.wus)、農(nóng)業(yè)管理文件(.mgt)、土壤化學(xué)文件(.chm)、池塘數(shù)據(jù)文件(.pnd)和河流水質(zhì)文件(.swq)，輸入完畢后可以通過模型界面修改編輯。

3 模型參數(shù)率定與驗證

模型參數(shù)的率定，從空間上來說，流域從上到下，從支流到干流；從時間上來說，時段從粗到細(xì)；從校準(zhǔn)的內(nèi)容來說，先徑流后污染負(fù)荷的順序，依次校準(zhǔn)。SWAT模型可以進(jìn)行人工校準(zhǔn)也可以進(jìn)行自動校準(zhǔn)，本研究區(qū)采用兩種方法相結(jié)合進(jìn)行參數(shù)校準(zhǔn)。在確定模型的敏感性參數(shù)之后，利用模型自帶的自動校準(zhǔn)分析方法SCE-UA (Shuffled Complex Evolution Algorithm )來確定這些敏感性參數(shù)的初始取值，將模型自動校準(zhǔn)優(yōu)化后的參數(shù)代入模型模擬，比較分析模擬值與實測值，如果3個校準(zhǔn)指標(biāo)滿足要求，采用自動校準(zhǔn)之后敏感性參數(shù)值；如果校準(zhǔn)指標(biāo)不滿足要求，則需要對敏感性參數(shù)采用人工校準(zhǔn)的方法進(jìn)行相應(yīng)的調(diào)整，直到評價指標(biāo)在允許的范圍之內(nèi)。

采用龍溪河瀘縣出口監(jiān)測斷面的2006-2010年實測流量、負(fù)荷數(shù)據(jù)，以2005年的數(shù)據(jù)作為預(yù)熱，2006-2008年數(shù)據(jù)進(jìn)行率定，2009-2010年的數(shù)據(jù)進(jìn)行驗證。

3.1 參數(shù)敏感性分析

采用模型自帶的敏感性分析模塊對徑流和氮、磷負(fù)荷進(jìn)行敏感性分析，選取幾個排序在前的參數(shù)作為調(diào)整參數(shù)，確定敏感性參數(shù)的最優(yōu)值，從而得到很好的模擬結(jié)果，并且節(jié)省工作量。敏感性參數(shù)排序見表3。

3.2 徑流的調(diào)試與驗證

評價模型結(jié)果是否能夠被接受，不存在統(tǒng)一的標(biāo)準(zhǔn)，一般認(rèn)為，當(dāng)R2>0.6，Ens>0.5，Re≤0.3%時，模型是準(zhǔn)確的[13]。月均徑流率定及驗證結(jié)果見表4，模擬值與實測值月均流量對比如圖2所示。徑流模擬結(jié)果可以被接受，引起模擬值與實際值存在的差異的原因是多方面的，包括流域的空間差異性分布所帶來的參數(shù)難確定性，以及模型結(jié)構(gòu)本身的模擬誤差和數(shù)據(jù)的準(zhǔn)確性等原因。

表3 敏感性參數(shù)表Tab.3 The sensitivity parameters table

表4 徑流模擬率定及驗證結(jié)果Tab.4 The runoff simulation results of calibration and verification

3.3 氮、磷污染負(fù)荷的調(diào)試與驗證

與徑流模擬相比，氮磷負(fù)荷的模擬影響因素更多，不確定性將更大，因此，氮磷負(fù)荷模擬的精度要求低于徑流模擬的要求。氮磷負(fù)荷率定及驗證結(jié)果見表5，模擬值與實測值TN月日均負(fù)荷對比見圖3，TP月日均負(fù)荷對比見圖4。氮磷負(fù)荷模擬值與實測值差異，與徑流模擬的精確程度密不可分，同時污染負(fù)荷的遷移轉(zhuǎn)化過程所引起的不確定性也是造成模擬值與實測值不一致的重要因素。

圖2 月均流量對比圖Fig.2 The comparison of average monthly flow

圖3 TN月日均負(fù)荷對比圖Fig.3 The comparison figure of TN monthly average load

圖4 TP月日均負(fù)荷對比圖Fig.4 The comparison figure of TP monthly average load

表5 氮磷負(fù)荷模擬率定及驗證結(jié)果Tab.5 The calibration and verification results of nitrogen and phosphorus load simulation

4 面源污染的特征分析

4.1 面源污染時間變化規(guī)律

依據(jù)研究區(qū)3個氣象站的降雨數(shù)據(jù)以及瀘縣1960-2010年多年平均降雨資料，將一年分成3個時期進(jìn)行分析，5-9月作為豐水期，占年降水量的73%，12-3月作為枯水期，占年降水量的9%，其他月份作為平水期。從5年的產(chǎn)量均值，可以看出產(chǎn)水量、TN產(chǎn)量、TP產(chǎn)量都呈現(xiàn)豐水期顯著多于平水期，平水期多于枯水期，豐水期產(chǎn)量比重大約為70%～90%，枯水期產(chǎn)量比重在5%以內(nèi)，可見，豐水期是面源污染產(chǎn)生的關(guān)鍵時期，決定了面源污染的嚴(yán)重程度。這與面源污染的產(chǎn)生機理有關(guān)，面源污染是以徑流沖刷為主要的動力來源，降雨量大的月份造成的土壤侵蝕就嚴(yán)重，隨之?dāng)y帶的氮磷負(fù)荷量就大。具體見表6。

表6 研究區(qū)各時期產(chǎn)量情況表Tab.6 The production situation in each period of the study area

4.2 面源污染空間分布規(guī)律

從空間角度分析面源污染的分布情況，子流域5、9、11、13的TN輸出系數(shù)最大，為5 599～5 903 kg/km2，且東部地區(qū)產(chǎn)量比西部地區(qū)大；子流域12、13、14的TP輸出系數(shù)最大，為184～197 kg/km2，且下游地區(qū)比上游地區(qū)大。綜合來看，5、11、13三個子流域氮磷負(fù)荷都較為嚴(yán)重的，4號子流域最輕，TN負(fù)荷東部地區(qū)較西部地區(qū)更嚴(yán)重，TP負(fù)荷下游比上游更嚴(yán)重。這不僅與子流域的下墊面差異性有關(guān)，還與各子流域污染排放情況密切相關(guān)。具體見圖5、圖6所示。

圖5 年均TN產(chǎn)量分布Fig.5 The annual yield distribution of TN

圖6 年均TP產(chǎn)量分布Fig.6 The annual yield distribution of TP

4.3 重要污染源識別

SWAT模型是基于流域水文地理特征的分布式水文模型，能夠模擬出各子流域污染負(fù)荷總產(chǎn)量，但是污染確定污染來源，本文為了確定點源與面源污染哪個是主要污染來源，哪種面源污染類型比重較大，采用將排放量與模擬結(jié)果進(jìn)行對比分析的方法。

4.3.1 主要污染源識別

從實際污染排放量與年均模擬值可以看出，面源污染是各子流域的主要污染排放類型，點源污染比重很小，且面源污染排放量大小與年均模擬值大小一一對應(yīng)，排放量大的子流域相應(yīng)的模擬值也大，因各子流域面源污染排放量明顯高于點源污染排放量，是各個子流域最主要的污染類型，應(yīng)該重點防治面源污染。具體見圖7、圖8。

圖7 TN污染源類型識別Fig.7 the source type identification of TN

圖8 TP污染源類型識別Fig.8 the source type identification of TP

4.3.2 主要面源污染類型識別

本文研究區(qū)考慮的面源污染類型包括農(nóng)村生活污染、畜禽養(yǎng)殖污染、化肥污染，從TN的面源污染分布看，10、13、14子流域面源污染排放量最大，各面源污染排放比重分別為：農(nóng)村生活污染50%、化肥污染40%、畜禽養(yǎng)殖污染10%；從TP面源污染分布看，10、13、14子流域面源污染排放量最大，相應(yīng)的子流域年均模擬產(chǎn)值也最多，各面源污染排放比重分別為：化肥污染55%、農(nóng)村生活污染23%、畜禽養(yǎng)殖污染22%，可見，化肥污染對面源污染產(chǎn)生影響最大，其次為農(nóng)村生活污染和畜禽養(yǎng)殖污染。農(nóng)村生活污染對TN模擬值影響最大，化肥污染對TP模擬值影響最大。具體見圖9、圖10。

圖9 TN面源類型識別圖Fig.9 the non-point source identification of TN

圖10 TP面源類型識別圖Fig.10 the non-point source identification of TP

4.4 污染關(guān)鍵區(qū)域識別

依據(jù)年均模擬TN、TP值分別分為5級，TN：Ⅰ級41～45 kg/hm2，Ⅱ級45～49 kg/hm2，Ⅲ級49～53 kg/hm2，Ⅳ級53～57 kg/hm2，Ⅴ級57～61 kg/hm2；TP：Ⅰ級1.0～1.2 kg/hm2，Ⅱ級1.2～1.4 kg/hm2，Ⅲ級1.4～1.6 kg/hm2，Ⅳ級1.6～1.8 kg/hm2，Ⅴ級1.8～2.0 kg/hm2。具體見圖11、圖12。

圖11 TN年均產(chǎn)量分級Fig.11 The classification of TN annual output

圖12 TP年均產(chǎn)量分級Fig.12 The classification of TP annual output

子流域3、9、11、13、15的TN負(fù)荷屬于Ⅴ級，子流域產(chǎn)量最大；子流域11、12、13、14、17的TP負(fù)荷屬于Ⅴ級，因此，11、13是非點源污染的關(guān)鍵區(qū)，是值得重點治理的區(qū)域。

5 結(jié) 語

本文首次以龍溪河瀘縣境內(nèi)的面源污染為研究對象，以地理信息系統(tǒng)(GIS)和SWAT模型為工具，探討了SWAT模型在研究區(qū)的適用性；定量計算研究區(qū)內(nèi)面源污染負(fù)荷，并分析了其特征，主要結(jié)論如下。

(1)采用5年資料徑流和水質(zhì)監(jiān)測數(shù)據(jù)，選取2000-2005年的數(shù)據(jù)作為預(yù)熱期，2006-2008年數(shù)據(jù)進(jìn)行率定，2009-2010年數(shù)據(jù)進(jìn)行驗證，依據(jù)相對誤差(Re)、相關(guān)系數(shù)(R2)和Nash-Suttcliffe系數(shù)(Ens)這3個指標(biāo)作為模型校準(zhǔn)的評價指標(biāo)，來表征模型模擬結(jié)果的優(yōu)劣，結(jié)果表明，模擬值的精度在允許范圍之內(nèi)，敏感性參數(shù)的取值具有一定的可信性，所建立的SWAT模型在研究區(qū)具有適用性。

(2)年產(chǎn)水量、TN產(chǎn)量、TP產(chǎn)量三者之間存在著一定的相關(guān)關(guān)系，均呈現(xiàn)豐水期顯著多于平水期，平水期多于枯水期，豐水期產(chǎn)量比重大約為70%～90%，枯水期產(chǎn)量比重在5%以內(nèi)，可見豐水期是面源污染產(chǎn)生的關(guān)鍵時期，決定了非點源污染的嚴(yán)重程度。

(3)各子流域的污染排放量與模擬年均總值是基本對應(yīng)的，并且面源污染是各個子流域最主要的污染類型?；饰廴九c農(nóng)村生活污染在面源污染中所占比重更大，影響也更大，對于TN模擬值影響最大的是農(nóng)村生活污染，對TP模擬值影響最大的是化肥污染。

(4)從年均模擬值看，由于11、13是面源污染的關(guān)鍵區(qū)，位于云錦鎮(zhèn)，是值得重點關(guān)注與治理的區(qū)域。

作為分布式水文模型代表的SWAT模型，盡管考慮了時空的異質(zhì)性，將每個子流域劃分為多個水文響應(yīng)單元進(jìn)行模擬，但模擬單元劃分粗略決定輸入數(shù)據(jù)向模擬單元尺度轉(zhuǎn)換時參數(shù)集總程度的高低[14]，同時，輸入數(shù)據(jù)的分辨率及精度、SWAT模型結(jié)構(gòu)本身的不完善、模型參數(shù)難以準(zhǔn)確獲取都造成模擬結(jié)果的不確定性，需要后續(xù)進(jìn)一步研究。