基于SWAT模型的古藺河流域面源污染模擬研究

王文章，程 艷，敖天其,2，黎小東

(1.四川大學 水利水電學院，成都 610065；2.四川大學 水力學與山區(qū)河流開發(fā)保護國家重點實驗室，成都 610065)

0 引 言

隨著點源污染得到了有效治理與控制，面源污染已成為國內(nèi)水環(huán)境污染治理研究的熱點[1]。SWAT模型以擴展模塊的形式嵌入與Arcgis軟件中，能夠在不同土壤類型和土地利用的大尺度流域內(nèi)模擬流域產(chǎn)流、產(chǎn)沙以及營養(yǎng)物負荷等的變化情況,已成功應(yīng)用于國內(nèi)多個流域的非點源污染模擬研究[2]。

子流域劃分是分布式水文模型模擬的最基本步驟，也是保證模擬精度的前提。流域的建模過程和模擬結(jié)果都會受到子流域大小、規(guī)模和數(shù)量的影響，不同的子流域劃分的參數(shù)敏感性差異也較大。張雪松等[3](2004)認為模型輸入?yún)?shù)的空間集總程度會隨子流域劃分數(shù)量大小而改變，從而對流域模擬結(jié)果產(chǎn)生一定影響; 喻曉等[4](2017)的研究表面，較低的亞流域劃分水平會造成模型結(jié)果不穩(wěn)定；王國波等(2016)[5]探討了不同子流域數(shù)量與流域關(guān)鍵污染源區(qū)(CSAs)的相關(guān)關(guān)系，表明相對于SWAT模型自動確定的集水閾值，運用合理的方法確定研究流域的最佳集水閾值更能準確的識別CSAs。目前子流域劃分的最佳集水閾值確定方法有水系河網(wǎng)密度法、流域?qū)挾确植纪魄蟆⑺捣中尉S數(shù)法等，其中河網(wǎng)密度法、流域?qū)挾确植纪魄蠓ㄟm用于較大流域，而水系分維值反映流域地形地貌特征，更加適用于小流域的最佳集水閾值推求[6]。考慮到研究區(qū)流域面積，因此本文運用分形理論確定流域最佳集水閾值，基于SWAT模型進行合理的子流域劃分與流域面源污染模擬，以期能夠精確的識別流域重點治理區(qū)域。

1 研究區(qū)概況

古藺縣坐落于四川盆地的南緣，地處東經(jīng)105°30′-106°20′，北緯27°41′-28°20′。古藺河是古藺縣境內(nèi)最大的河流，是赤水河右岸一級支流，流域占古藺縣各鄉(xiāng)鎮(zhèn)總面積的大部分，集水面積為965 km2，干流河長70.75 km，河道平均比降1.84%。根據(jù)古藺縣境內(nèi)國家水文站1959-2007年實測的徑流資料計算，流域多年的平均年徑流深為422 mm，多年的平均流量5.07 m3/s，年徑流量1.6 億m3。徑流在年內(nèi)的分配不均勻，豐水期(5-9月)的徑流占年徑總流的65.8%，枯水期11月-翌年3月約占年徑總流的21.7%。最枯月為12月，僅占年徑流的3.7%。

研究區(qū)包含的8個鄉(xiāng)鎮(zhèn)，包括太平鎮(zhèn)、石屏鄉(xiāng)、魚化鄉(xiāng)、龍山鎮(zhèn)、永樂鎮(zhèn)、護家鄉(xiāng)、古藺鎮(zhèn)、德躍鎮(zhèn)；總?cè)丝跒?44 328，占全縣人口的4.01% ，其中農(nóng)業(yè)人口283 408，城鎮(zhèn)人口60 920，各鄉(xiāng)鎮(zhèn)的平均人口密度為360人/km2，遠遠超過了中國平均人口密度(130 人/km2)。研究區(qū)包含的8個鄉(xiāng)鎮(zhèn)2013年種植小麥7 027 hm2，水稻5 087 hm2，玉米6 840 hm2，飼養(yǎng)牛36 795頭，生豬151 802頭，羊17 221只，雞406 714只，其中包含集中式畜禽養(yǎng)殖與散養(yǎng)畜禽養(yǎng)殖。

圖1 古藺河流域Fig.1 Gulin River Basin

2 研究方法

2.1 模型數(shù)據(jù)庫建立

(1)數(shù)字高程DEM：古藺河流域的地形數(shù)據(jù)來源于中國科學院計算機網(wǎng)絡(luò)信息中心國際科學數(shù)據(jù)鏡像網(wǎng)站(http:∥datamirror.csdb.cn)提供的30m分辨率的DEM(Digital Elevation Model)數(shù)據(jù)。在DEM輸入SWAT之前在Arcgis中進行處理，通過填洼處理避免水流逆流或者匯流[7-8]，SWAT中輸入的研究區(qū)DEM見圖2。

圖2 數(shù)字高程Fig.2 DEM

(2)氣象數(shù)據(jù)：古藺河流域有古藺水文站一個國家水文站點，位于古藺鎮(zhèn)，本次研究收集到2000-2013年逐時降雨數(shù)據(jù)，由于分布式水文模型中要使用日降雨數(shù)據(jù)，將其整理為日降雨，并收集到古藺水文站2000-2010年的逐月徑流數(shù)據(jù)；研究區(qū)的太陽輻射、日相對濕度、風速等氣象數(shù)據(jù)使用SWAT內(nèi)置的天氣發(fā)生器模擬獲得。

(3)土壤類型數(shù)據(jù)：研究區(qū)土壤數(shù)據(jù)是從FAO網(wǎng)上下載的HWSD 數(shù)據(jù)，并自行剪切出中國區(qū)的土壤柵格圖，其數(shù)據(jù)值(Value)有效范圍是11 000～11 935，古藺河流域的土壤數(shù)據(jù)就是從中國區(qū)的土壤柵格圖裁剪而來。在下載的HWSD數(shù)據(jù)庫中打開 HWSD_DATA 表，可以查找土壤的部分屬性值，其余屬性值通過《四川土種志》查找與用SPAW軟件計算得到[9]，并新增到SWAT模型的土壤數(shù)據(jù)庫中。

表1 水文氣象站點概況Tab.1 Hydrometeorology station

(4)土地利用數(shù)據(jù)：古藺流域的土地利用/植被覆蓋類型數(shù)據(jù)來源于四川省古藺縣林業(yè)局的相關(guān)土地利用信息，利用ARCGIS軟件制成了古藺縣土地利用圖。原始數(shù)據(jù)中包含了水田、旱地、有林地、草地、村莊、河流水面等各種土地利用分類。在分布式水文模型中，主要將土地利用類型分為9類，其中林地占54.08%(FRST),耕地占35.14%(AGRL),草地占7%(PAST)，因此古藺縣主要以林業(yè)和農(nóng)業(yè)種植為主。由于美國的土地利用類型與我國的屬性相似，因此在本次研究中土地利用屬性直接采用SWAT模型系統(tǒng)默認的土地利用類型屬性數(shù)據(jù)[10]。

圖3 土地利用圖Fig.3 Landuse

圖4 土壤類型圖Fig.4 Soil

(4)污染源數(shù)據(jù)：古藺河流域的污染源中按照不同類型又可劃分為若干種類；其中點源主要包括工業(yè)污染、城鎮(zhèn)居民所產(chǎn)生的垃圾廢水，非點源包括散養(yǎng)畜禽污染、化肥污染、農(nóng)村居民所產(chǎn)生的垃圾廢水，此類統(tǒng)計數(shù)據(jù)來源于古藺縣環(huán)保局、農(nóng)業(yè)局、統(tǒng)計局等。本次研究假設(shè)各統(tǒng)計因子在行政區(qū)內(nèi)等密度分布，在實際計算時，采用面積比例法，獲得古藺河流域的統(tǒng)計數(shù)據(jù)信息。

2.2 水系分維數(shù)計算

計盒維數(shù)法是求解分形維數(shù)的計算方法之一[11]，其基本思路是用不同尺寸的盒子覆蓋分形體，并統(tǒng)計對應(yīng)尺寸的盒子數(shù)量。當水系作為分形體時，可以利用Arcgis的空間分析模塊選定某一集水閾值Si，提取研究區(qū)的柵格河網(wǎng)，然后不斷地改變柵格邊長r的大小，統(tǒng)計不同邊長的柵格的數(shù)量N，當r→0 時，可以得到研究區(qū)水系分形維數(shù)Di：

(1)

3 結(jié)果與分析

3.1 最佳集水閾值確定

本文選取300 hm2到4 500 hm2之間共15個集水閾值，計算對應(yīng)于10種柵格邊長下的分維數(shù)。圖5(a)為利用計盒維數(shù)法求得的分維數(shù)與集水閾值關(guān)系曲線，可以看出，隨著集水閾值的不斷增大，水系分維值開始減小，變化速度也逐漸降低。

圖5 水系分維值與集水閾值關(guān)系Fig.5 The relationship of fractal dimension and threshold of drainage area

由于集水閾值越大，河網(wǎng)密度越小，所能反映的流域信息就越少，而集水閾值過小，又會導(dǎo)致流域出現(xiàn)過多的偽河道，使模擬結(jié)果與實際情況相差過大，因此需要保證流域不出現(xiàn)過多偽河道的情況下，河網(wǎng)密度盡量的大。在此基礎(chǔ)上，建立分維數(shù)與集水閾值的二階導(dǎo)函數(shù)關(guān)系(即分維值的變化加速度)如圖5(b)，當集水閾值小于1 500 hm2時，分維值變化加速度極大，1 500 hm2后分維值變化加速度趨于穩(wěn)定，說明該值即是相關(guān)研究中提到的閾值“拐點”[12,13]，因此本文選取集水閾值為1 500 hm2進行合理的子流域劃分。

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

通過模型自帶的敏感性分析模塊對研究區(qū)參數(shù)進行敏感性分析，確定了流域中影響徑流的6個敏感性參數(shù)，以及影響氮磷負荷敏感的6個參數(shù)，具體見表2。本次研究收集到古藺河流域古藺站2000-2010年實測月徑流和古藺河水質(zhì)監(jiān)測斷面2012-2013年的水質(zhì)數(shù)據(jù)。水質(zhì)模擬需要控制斷面的徑流數(shù)據(jù)，因此首先率定和驗證研究區(qū)的徑流過程，調(diào)試模型使其精度滿足模擬精度要求，一般要求校準期以及驗證期徑流模擬值和實測值的Nash效率系數(shù)Nash≥0.5，決定系數(shù)R2≥0.6[6]。

表2 敏感性參數(shù)表Tab.2 Sensitive parameters

(1)徑流模擬：將2000年作為預(yù)熱期，2001-2007年作為率定期，2008-2010年作為驗證期，古藺站徑流模擬精度以及月徑流實測值與模擬值對比見表3。可見SWAT模型在古藺河流域有較好的適用性，率定期Nash效率系數(shù)0.78，決定系數(shù)R2為0.92；驗證期Nash效率系數(shù)為0.71，決定系數(shù)R2為0.88；模擬結(jié)果令人滿意。

表3 研究區(qū)徑流模擬率定及驗證結(jié)果Tab.3 Results of runoff simulation and calibration

(2)總氮、總磷模擬：將古藺河流域出口太平渡斷面2012年實測總氮、總磷數(shù)據(jù)(古藺縣環(huán)保局提供)用于模型率定，2013年的數(shù)據(jù)用于模型驗證。氮、磷污染負荷模擬的精度要求對徑流模擬的相對較低，只要保證氮、磷負荷的模擬值與實測值得Nash≥0.5，決定系數(shù)R2≥0.6即可。通過模型的調(diào)試，古藺河流域水質(zhì)監(jiān)測斷面的氮、磷污染負荷率定與驗證結(jié)果滿足精度要求，曲線擬合較好。TN與TP模擬值與實測值得月平均負荷對比見圖6。

圖6 總氮、總磷模擬結(jié)果Fig.6 results of SWAT in simulating TN and TP

3.3 總氮、總磷輸出系數(shù)的空間分布

對古藺河各子流域單位面積內(nèi)的非點源污染負荷(即各子流域TN、TP的輸出系數(shù))進行空間分析(見圖7)，由圖中可以看出子流域TN、TP的輸出系數(shù)空間分布特征一致，TN、TP輸出系數(shù)最大的為6、11、14、18、30子流域，值分別為46.73～67.67、5.72～8.29 kg/hm2。根據(jù)研究區(qū)的概況及流域劃分，11、14、30子流域?qū)儆诠盘A鎮(zhèn)的行政區(qū)內(nèi)，子流域18屬于石屏鄉(xiāng)的行政區(qū)內(nèi)，子流域6屬于德躍鎮(zhèn)的行政區(qū)內(nèi)；而古藺鎮(zhèn)是古藺河流域內(nèi)面積最大，農(nóng)業(yè)人口、畜禽養(yǎng)殖最多的鄉(xiāng)鎮(zhèn)，石屏鄉(xiāng)的農(nóng)業(yè)人口、畜禽養(yǎng)殖較多，但面積最小，因此造成行政區(qū)內(nèi)的子流域污染物輸出系數(shù)較大。

圖7 各子流域單位面積TN、TP輸出Fig.7 Unit area output of TN and TP in Subbasin

3.4 污染源貢獻率分析

本次研究中的非點源污染包括農(nóng)村居民生活所產(chǎn)生的污水、散養(yǎng)畜禽養(yǎng)殖以及化肥，分別模擬有相應(yīng)污染源輸入和無輸入時研究區(qū)出口斷面負荷的變化，根據(jù)變化的差值可以估算出各種污染源對研究區(qū)污染負荷的貢獻。從圖8中可以看出，大多數(shù)子流域中各類污染源的排列順序為化肥產(chǎn)生的污染>畜禽養(yǎng)殖產(chǎn)生的污染>居民生活產(chǎn)生的污染；尤其是TP圖中可以看出，化肥產(chǎn)生的污染TP產(chǎn)量遠遠高于畜禽養(yǎng)殖污染和化肥污染，可見化肥污染與畜禽養(yǎng)殖所造成的污染是古藺河流域非點源污染的主要原因。

圖8 各類污染源輸出Fig.8 Output of each pollution source

4 結(jié) 論

(1)基于分形理論方法，計算出古藺河流域最佳集水閾值為1 500 hm2，該閾值下子流域劃分個數(shù)為37個。

(2)在古藺河流域建立SWAT模型對氮磷負荷進行研究，模型率定期與驗證期的徑流模擬精度Nash效率系數(shù)值大于0.70，決定系數(shù)R2大于0.88；總氮、總磷模擬率定期與驗證期Nash效率系數(shù)值大于0.5，決定系數(shù)R2大于0.8，模擬評價參數(shù)均滿足精度要求，表明SWAT模型在古藺河流域農(nóng)村面源污染模擬研究中具有一定的適用性。

(3)古藺河流域2013年總氮輸出量為2 761.72 t，總磷輸出量為338.27 t，古藺鎮(zhèn)、德躍鎮(zhèn)、石屏鄉(xiāng)污染負荷輸出強度占前三，是需要重點治理的鄉(xiāng)鎮(zhèn)。

(4)農(nóng)業(yè)化肥是流域主要污染源，農(nóng)田中化肥施用的比例失調(diào)會大大降低農(nóng)作物對化肥的利用率，同時增大化肥的流失，所以嚴格控制農(nóng)業(yè)用肥量，進行測土配方施肥是研究區(qū)面源污染治理的需要采取重要措施。

□