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

馬亞麗，白祖暉，敖天其

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

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

1 流域概況

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

2 數據來源及處理

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

表1 研究區投影參數表Tab.1 The projection parameters of the study area

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

2.1 流域內子流域劃分

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

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

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

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

2.3 模型參數的寫入及編輯

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

3 模型參數率定與驗證

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

采用龍溪河瀘縣出口監測斷面的2006-2010年實測流量、負荷數據，以2005年的數據作為預熱，2006-2008年數據進行率定，2009-2010年的數據進行驗證。

3.1 參數敏感性分析

采用模型自帶的敏感性分析模塊對徑流和氮、磷負荷進行敏感性分析，選取幾個排序在前的參數作為調整參數，確定敏感性參數的最優值，從而得到很好的模擬結果，并且節省工作量。敏感性參數排序見表3。

3.2 徑流的調試與驗證

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

表3 敏感性參數表Tab.3 The sensitivity parameters table

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

3.3 氮、磷污染負荷的調試與驗證

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

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

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

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

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

4 面源污染的特征分析

4.1 面源污染時間變化規律

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

表6 研究區各時期產量情況表Tab.6 The production situation in each period of the study area

4.2 面源污染空間分布規律

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

圖5 年均TN產量分布Fig.5 The annual yield distribution of TN

圖6 年均TP產量分布Fig.6 The annual yield distribution of TP

4.3 重要污染源識別

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

4.3.1 主要污染源識別

從實際污染排放量與年均模擬值可以看出，面源污染是各子流域的主要污染排放類型，點源污染比重很小，且面源污染排放量大小與年均模擬值大小一一對應，排放量大的子流域相應的模擬值也大，因各子流域面源污染排放量明顯高于點源污染排放量，是各個子流域最主要的污染類型，應該重點防治面源污染。具體見圖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 主要面源污染類型識別

本文研究區考慮的面源污染類型包括農村生活污染、畜禽養殖污染、化肥污染，從TN的面源污染分布看，10、13、14子流域面源污染排放量最大，各面源污染排放比重分別為：農村生活污染50%、化肥污染40%、畜禽養殖污染10%；從TP面源污染分布看，10、13、14子流域面源污染排放量最大，相應的子流域年均模擬產值也最多，各面源污染排放比重分別為：化肥污染55%、農村生活污染23%、畜禽養殖污染22%，可見，化肥污染對面源污染產生影響最大，其次為農村生活污染和畜禽養殖污染。農村生活污染對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 污染關鍵區域識別

依據年均模擬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年均產量分級Fig.11 The classification of TN annual output

圖12 TP年均產量分級Fig.12 The classification of TP annual output

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

5 結 語

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

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

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

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

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

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