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

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

(1.四川大學(xué) 水利水電學(xué)院，成都 610065；2.四川大學(xué) 水力學(xué)與山區(qū)河流開發(fā)保護(hù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，成都 610065)

0 引 言

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

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

1 研究區(qū)概況

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

研究區(qū)包含的8個(gè)鄉(xiāng)鎮(zhèn)，包括太平鎮(zhèn)、石屏鄉(xiāng)、魚化鄉(xiāng)、龍山鎮(zhèn)、永樂鎮(zhèn)、護(hù)家鄉(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，遠(yuǎn)遠(yuǎn)超過了中國(guó)平均人口密度(130 人/km2)。研究區(qū)包含的8個(gè)鄉(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ù)庫(kù)建立

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

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

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

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

表1 水文氣象站點(diǎn)概況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è)種植為主。由于美國(guó)的土地利用類型與我國(guó)的屬性相似，因此在本次研究中土地利用屬性直接采用SWAT模型系統(tǒng)默認(rèn)的土地利用類型屬性數(shù)據(jù)[10]。

圖3 土地利用圖Fig.3 Landuse

圖4 土壤類型圖Fig.4 Soil

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

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

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

(1)

3 結(jié)果與分析

3.1 最佳集水閾值確定

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

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

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

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

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

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

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

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

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

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

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

對(duì)古藺河各子流域單位面積內(nèi)的非點(diǎn)源污染負(fù)荷(即各子流域TN、TP的輸出系數(shù))進(jìn)行空間分析(見圖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 污染源貢獻(xiàn)率分析

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

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

4 結(jié) 論

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

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

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

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

□