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

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

0 引 言

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

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

1 研究區(qū)概況

古藺縣坐落于四川盆地的南緣，地處東經105°30′-106°20′，北緯27°41′-28°20′。古藺河是古藺縣境內最大的河流，是赤水河右岸一級支流，流域占古藺縣各鄉(xiāng)鎮(zhèn)總面積的大部分，集水面積為965 km2，干流河長70.75 km，河道平均比降1.84%。根據古藺縣境內國家水文站1959-2007年實測的徑流資料計算，流域多年的平均年徑流深為422 mm，多年的平均流量5.07 m3/s，年徑流量1.6 億m3。徑流在年內的分配不均勻，豐水期(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)；總人口為344 328，占全縣人口的4.01% ，其中農業(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 模型數據庫建立

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

圖2 數字高程Fig.2 DEM

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

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

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

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

圖3 土地利用圖Fig.3 Landuse

圖4 土壤類型圖Fig.4 Soil

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

2.2 水系分維數計算

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

(1)

3 結果與分析

3.1 最佳集水閾值確定

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

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

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

3.2 模型參數率定與驗證

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

表2 敏感性參數表Tab.2 Sensitive parameters

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

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

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

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

3.3 總氮、總磷輸出系數的空間分布

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

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

3.4 污染源貢獻率分析

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

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

4 結 論

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

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

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

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

□