陳俊宏， 張利華*， 馬永明

(1.中國地質(zhì)大學(xué)(武漢) 地理與信息工程學(xué)院， 武漢 430074； 2.昭通學(xué)院地理科學(xué)與旅游學(xué)院， 云南 昭通 657000)

對于流域內(nèi)水質(zhì)污染的研究，采用定量化的水文模型計算和評估，已成為有效量化的手段[1]，各種分布式水文模型如WEPPO模型、AGNPS模型、EPIC模型和SWAT模型等大量出現(xiàn)，在國內(nèi)外的研究中應(yīng)用十分廣泛[2]，但各類模型通常都需要輸入大量的參數(shù)，但目前的水文站點分布尚有不足，如我國就存在水文站點多分布在經(jīng)濟發(fā)達(dá)的中東部地區(qū)的狀況[3]，許多山區(qū)地區(qū)缺少資料或者資料不足，各類模型難以應(yīng)用，治理受到制約，因此無資料地區(qū)的水文研究成為水文預(yù)報研究中的重點和難點.

目前，無資料地區(qū)的研究多應(yīng)用區(qū)域化方法，如Narbondo等[4]以烏拉圭國內(nèi)的13個流域為例研究了區(qū)域化方法，開發(fā)了基于相似性概念的流域分區(qū)方法，建立了降雨-徑流模型參數(shù)與流域物理屬性的函數(shù)關(guān)系，以獲取無資料流域的參數(shù)，從而進(jìn)行模擬.而在區(qū)域化方法中最易使用的是基于水文相似性的參數(shù)移植法[5]，Dumedah等[6]在加拿大安大略省南部的8個流域應(yīng)用數(shù)據(jù)同化進(jìn)行了徑流模擬，根據(jù)得到的結(jié)果，他們認(rèn)為模型參數(shù)值具有高度的共同性，可以將另一個流域的參數(shù)應(yīng)用到給定的流域，丁洋等[7]對水文資料缺失的湖北省松滋市小南海流域應(yīng)用了參數(shù)移植法進(jìn)行研究，將具有相似性的洞庭湖等流域的水文參數(shù)移植到小南海流域，得到了流域的徑流模擬結(jié)果，Cho等[8]將美國喬治亞州西南部沿海平原流域的經(jīng)過了率定優(yōu)化的參數(shù)移植到同一地理區(qū)域的其他五個流域，不再進(jìn)行參數(shù)校準(zhǔn)，得到了較好的模擬結(jié)果，樊琨[9]對甘肅省涇河上游區(qū)、三水河等流域應(yīng)用了參數(shù)移植法，模擬結(jié)果良好.

SWAT模型可應(yīng)用范圍廣，在國內(nèi)外的各個流域均表現(xiàn)出適用性，可用于非點源污染分析[10].有研究應(yīng)用模型參數(shù)移植，實現(xiàn)了缺資料地區(qū)的SWAT模型構(gòu)建[11-13]，SWAT模型在犟河流域的上級流域已經(jīng)表現(xiàn)出了適用性[14-16]，陳昊榮[17]在漢江流域成功應(yīng)用了參數(shù)移植法構(gòu)建了SWAT模型.犟河流域位于南水北調(diào)工程核心區(qū)，在調(diào)水后曾表現(xiàn)出水質(zhì)污染嚴(yán)重狀態(tài)[18]，亟待解決水環(huán)境治理問題，但是由于水文資料缺乏，流域內(nèi)的相關(guān)研究較少，因此基于SWAT模型，應(yīng)用參數(shù)移植法進(jìn)行犟河流域的SWAT模型構(gòu)建和非點源污染模擬，為該流域污染防治提供研究依據(jù).

1 材料和方法

1.1 研究區(qū)概況

犟河是長江流域漢江右岸的二級支流，同時也是堵河的一級支流，其起于十堰市茅箭區(qū)大川鎮(zhèn)，穿越南水北調(diào)核心十堰市，在黃龍鎮(zhèn)東灣村與堵河交匯，再匯入漢江至丹江口水庫，作為十堰市的納污河，長期接受十堰市城市污染排放.犟河干流全長50.2 km，流域范圍為32°29′21″N～32°42′37″N，110°31′38″E～110°43′26″E，流域面積316 km2，流域地勢上南高北低，西南高、東北低，整體海拔139～1 441 m，為亞熱帶季風(fēng)氣候，雨熱同期，季節(jié)變化明顯，生態(tài)基流小，在秋冬季枯水期流量較小.流域地理位置見圖1，其位于上級流域堵河流域的東北角.

圖1 研究區(qū)概況圖 Fig.1 Geographical location

1.2 模型簡介

SWAT模型是由1 013 個中間變量、701 個方程組成的綜合的、分布式的模型體系[19]，經(jīng)過了多年的開發(fā)和應(yīng)用，已經(jīng)有了比較成熟的集成軟件體系如ArcSWAT、SWAT-CUP方便研究人員進(jìn)行分析.SWAT模型構(gòu)建簡單，輸入高程數(shù)據(jù)、土壤數(shù)據(jù)、土地利用數(shù)據(jù)、氣象數(shù)據(jù)等基礎(chǔ)地理數(shù)據(jù)后，即可初步構(gòu)建流域的SWAT模型，根據(jù)研究方向，再輸入其他管理數(shù)據(jù)，并對模型進(jìn)行率定和驗證，便得到模擬結(jié)果.

1.3 數(shù)據(jù)處理

從中國氣象同化驅(qū)動數(shù)據(jù)集CMADS[20]中裁剪獲得了研究區(qū)域的十堰站氣象數(shù)據(jù)；從地理空間數(shù)據(jù)云(http：//www.gscloud.cn)，下載到了包含犟河流域的DEM數(shù)據(jù)，空間分辨率為30 m，裁剪得到的DEM數(shù)據(jù)見圖2(a)；土地利用數(shù)據(jù)主要通過目視遙感影像解譯獲取，下載獲取了包含研究區(qū)域的Landsat 8遙感影像，并結(jié)合GoogleEarth在線地圖，在ENVI 5.3軟件中進(jìn)行遙感解譯和裁剪等，得到了研究區(qū)域的土地利用類型數(shù)據(jù)(圖2(b))；土壤數(shù)據(jù)參考姜曉峰等[21]論述的SWAT模型土壤數(shù)據(jù)庫構(gòu)建方法，基于HWSD數(shù)據(jù)庫，建立研究區(qū)的土壤數(shù)據(jù)庫(圖2(c))，犟河流域土壤主要有四種類型：石灰性沖積土(FLc)、兩種發(fā)育程度不同的弱發(fā)育淋溶土(LVh1、LVh2)、不飽和始成土(CMd).實測徑流數(shù)據(jù)及實測水質(zhì)數(shù)據(jù)收集自十堰市水文資源局，徑流數(shù)據(jù)包括了竹山站和黃龍灘站兩個水文站2014—2017年的數(shù)據(jù)，水質(zhì)數(shù)據(jù)主要包括氨氮和總磷.其他管理數(shù)據(jù)如農(nóng)業(yè)數(shù)據(jù)收集自十堰市農(nóng)業(yè)農(nóng)村局，包括張灣區(qū)各鄉(xiāng)鎮(zhèn)的農(nóng)業(yè)化肥施用量數(shù)據(jù)；禽畜、水產(chǎn)養(yǎng)殖數(shù)據(jù)、人口生活污水、工業(yè)排放數(shù)據(jù)收集自十堰市統(tǒng)計年鑒及十堰市生態(tài)環(huán)境局.

此外，輸入SWAT模型土壤數(shù)據(jù)庫的部分參數(shù)使用了軟件SPAW(Soil Plant Atmosphere Water)進(jìn)行輔助計算，如土壤濕密度、土壤層有效持水量、飽和水利傳導(dǎo)系數(shù).數(shù)據(jù)庫中USLE土壤可蝕性因子K，利用土壤粒徑分布進(jìn)行計算得到：

Kusle=fcsand×fcl-si×forgc×fhisand，

(1)

其中

fcsand=

fhisand=

式中，ms代表在美國制土壤粒徑分類下，土壤中粒徑為0.05～2.00 mm的砂礫的百分含量，msilt代表粒徑為0.002～0.05 mm的粉砂的百分含量，mc代表粒徑小于0.002 mm的粘粒的百分含量，Corg代表土壤中有機碳的含量.

圖2 基礎(chǔ)地理數(shù)據(jù)Fig.2 Basic geographic data

1.4 模型的率定和驗證

分布式水文模型往往有多種需要結(jié)合實測資料進(jìn)行率定的參數(shù)[22-23]，在缺乏資料的地區(qū)，為了提高模型模擬的精度，目前國內(nèi)外研究者一般采用區(qū)域化方法[24-25]解決參數(shù)問題，參數(shù)移植法[26]是其中的一種.朱阿興等[27]提出地理環(huán)境越相似，地理特征越相近.參數(shù)移植應(yīng)用了相同的思想，Parajka等[28]對奧地利的320個流域進(jìn)行了參數(shù)移植的試驗，他們的結(jié)果顯示，當(dāng)流域的自然屬性如地貌、土壤、土質(zhì)等相似時，除了高寒地區(qū)和平坦流域，區(qū)域化處理之后的模擬的結(jié)果良好.Kokkonen等[29]研究參數(shù)移植時得到結(jié)論，如果流域水文特征相似，那么可以將有實測數(shù)據(jù)的流域的全部率定參數(shù)都移植到相似的流域.

犟河流域是堵河流域的子流域，兩者位于相似氣候區(qū)，有相似的降水、地貌等特征，因此，應(yīng)用參數(shù)移植法，將堵河流域參數(shù)率定結(jié)果用于犟河流域進(jìn)行犟河流域的模型構(gòu)建，利用SWAT-CUP進(jìn)行參數(shù)校準(zhǔn)以降低模型不確定性[30].

基于實測徑流數(shù)據(jù)，先進(jìn)行參數(shù)敏感性分析，再對敏感性參數(shù)進(jìn)行率定.敏感性分析采用GLUE算法，設(shè)定模擬次數(shù)5 000 次；率定采用SUFI-2算法，單次率定設(shè)定模擬次數(shù)500 次，迭代，以線性回歸決定系數(shù)R2和納什效率系數(shù)NSE最大化作為目標(biāo)函數(shù)，參數(shù)選擇參考牛利強[16]及戴楓勇[31].堵河流域的徑流模擬以2014—2016年為率定期，以2017年為驗證期，犟河流域的非點源污染模擬以2014—2015年為率定期，以2016—2017年對氨氮模擬結(jié)果進(jìn)行驗證，以2018年對總磷模擬結(jié)果進(jìn)行驗證.

2 結(jié)果和分析

2.1 徑流模擬

SWAT模型中對模擬結(jié)果的評價指標(biāo)，采用納什效率系數(shù)NSE及線性回歸決定系數(shù)R2[12， 32-33]，一般認(rèn)為在進(jìn)行徑流模擬時，R2>0.6和NSE>0.5時，模擬的結(jié)果可以接受[12].堵河流域有黃龍灘、竹山兩個水文站點，基于SUFI-2算法進(jìn)行迭代模擬，每次設(shè)定模擬次數(shù)為500 次，迭代率定計算了12 次，R2和NSE不再增大，且滿足精度要求，率定結(jié)束，最終得到的徑流模擬率定參數(shù)值見表1.

表1 堵河流域徑流模擬參數(shù)取值范圍及最佳值

將經(jīng)過率定的參數(shù)返回輸入模型，進(jìn)行2017年的徑流模擬.以2017年的竹山站數(shù)據(jù)作為驗證，計算得到的R2=0.575，徑流模擬與實測值相關(guān)性見圖3，徑流模擬結(jié)果見圖4.

圖3 2017年徑流模擬與實測相關(guān)性Fig.3 Correlation between runoff simulation and measurement in 2017

圖4 月徑流量模擬結(jié)果Fig.4 Monthly runoff simulation results

2.2 非點源污染模擬

堵河流域的徑流模擬完成，基于參數(shù)移植法，將堵河流域的徑流模擬參數(shù)輸入到犟河流域的模型之中，并進(jìn)行污染物的參數(shù)率定，率定方法與徑流模擬相似，只在參數(shù)選取上有所不同.以犟河流域東灣橋監(jiān)測斷面的氨氮和總磷水質(zhì)數(shù)據(jù)作為實測值，由于輸入模型的的污染負(fù)荷為總量，單位為kg，需要將濃度值轉(zhuǎn)換為總量，計算公式為：

L=F×c×10-3×t，

(2)

式中，L為污染負(fù)荷總量，單位為kg，F(xiàn)為斷面徑流量，單位為m3·s-1，t為時間，單位為秒，計算得到的數(shù)據(jù)，在率定時，輸入SWAT-CUP，以2016—2017年數(shù)據(jù)進(jìn)行氨氮的模擬驗證，以2018年數(shù)據(jù)進(jìn)行總磷的模擬驗證，得到的率定參數(shù)值見表2，模擬結(jié)果精度見表3，氨氮模擬結(jié)果見圖5，總磷模擬結(jié)果見圖6.

表2 犟河流域污染物模擬參數(shù)取值范圍及最佳值

表3 污染物模擬結(jié)果精度

圖5 氨氮模擬Fig.5 Ammonia nitrogen simulate

圖6 總磷模擬 Fig.6 Total phosphorus simulate

氨氮模擬，驗證期的R2達(dá)到0.76，滿足精度要求，效果較好；而總磷模擬，驗證期的R2為0.43，效果不夠好，但可以看到其年際變化與氨氮模擬結(jié)果相似，到2018年，污染量明顯減少.氨氮在2017年汛期，9、10月出現(xiàn)了高值，可能是由于降水量增大，土壤侵蝕、對地表的沖刷等，更多的污染物匯入到了犟河之中.從總體上來看，犟河流域在2014年、2015年的污染物排放量較大，水質(zhì)較差，在十堰市政府的逐漸重視和投資治理下，2018年后，水體中污染物含量已經(jīng)有了明顯降低.

圖7 犟河流域污染負(fù)荷與徑流量的月周期變化 Fig.7 Monthly variation of pollution load and runoff in the Jiang River Basin

2.3 氮磷污染時間分布特征分析

基于模擬結(jié)果進(jìn)行氨氮和總磷的時間分布特征分析，從圖7可以看出，犟河流域的氨氮量和總磷量與徑流量隨時間的變化趨勢相似，在汛期，徑流量大，氨氮和總磷污染量也相應(yīng)較大，至非汛期，徑流量和非點源污染負(fù)荷變小，其中氨氮的污染負(fù)荷最大值出現(xiàn)在2014年10月，總磷最大值出現(xiàn)在2017年9月，均對應(yīng)于汛期.2017年至2018年全年氨氮污染總量下降了70.63%，下降幅度最大(圖8)，2014年至2015年下降幅度為40.15%，2015年至2016年下降幅度為43.58%，從下降幅度來看，2018年治理效果良好，從氨氮減少的量來看，2014年至2015年，氨氮排放減少了45 738.27 kg，減少的總量最大.2016年至2017年，氨氮污染總量增加了24.73%，根據(jù)查到的《漢江流域氣候公報(降水)》數(shù)據(jù)，2017年，十堰市汛期的降水量較常年降水量偏多60%，其中8月28日—10月18日，全市降水量達(dá)到了527.8 mm，達(dá)到歷史同期的3.19倍，可能是降水造成了徑流量增大，更多的污染物通過徑流，匯入到了犟河之中.總體來看，經(jīng)過多年時間，污染量每年以約40%的速率下降，犟河的污染得到了有效治理.

圖8 犟河流域氨氮負(fù)荷年際變化Fig.8 Interannual variation of ammonia nitrogen load in Jiang River Basin

2.4 氮磷污染空間分布特征分析

利用SWAT模型模擬輸出的氨氮和總磷污染量，以2014年和2018年的數(shù)據(jù)作為對比，計算到每個子流域的單位面積上的污染負(fù)荷量.計算結(jié)果見圖9.對比結(jié)果顯示，無論是氨氮還是總磷的單位面積負(fù)荷，高值都主要集中在1、4、5、7、9、10、13、14、15、17、18、25、29等子流域，根據(jù)谷歌地球在線地圖及遙感土地利用分類結(jié)果，這些子流域正好對應(yīng)于城市建筑群，是居民居住的主要區(qū)域.此外，結(jié)合犟河流域的河網(wǎng)分布(圖1)，在財神溝和大西溝，自上游至下游，污染逐漸加重.犟河干流南岸污染高于北岸.同時，將2014年與2018年的氨氮及總磷單位面積負(fù)荷量在各子流域進(jìn)行對比(圖10)，氨氮的單位面積負(fù)荷在高值區(qū)有明顯的降低，而總磷值略有降低.兩者的模擬結(jié)果在空間分布上明顯相近.

圖9 犟河流域污染單位面積負(fù)荷的空間分布 Fig.9 Spatial distribution of pollution load per unit area in the stubborn River Basin

圖10 犟河流域各子流域污染物單位面積負(fù)荷變化Fig.10 Variation of pollutant unit area load in each sub basin of the Jiang River Basin

3 結(jié)果與討論

以缺少徑流資料的犟河流域為例，基于SWAT模型，運用參數(shù)移植法探索了徑流參數(shù)的移植.

1) 對犟河流域的上級流域——堵河流域進(jìn)行了徑流模擬，模擬期精度達(dá)到要求，驗證期竹山站結(jié)果R2=0.62，并將堵河流域的徑流模擬參數(shù)移植到犟河流域.

2) 對犟河流域的非點源污染進(jìn)行了模擬，驗證期非點源污染模擬結(jié)果(氨氮R2=0.76，總磷R2=0.43)，總磷的模擬結(jié)果精度偏低，但其和氨氮的模擬結(jié)果時空分布表現(xiàn)出相似的特征，表明SWAT模型在犟河流域是適用的.

3) 犟河流域經(jīng)過多年治理，污染負(fù)荷量在時空兩層次上均表現(xiàn)出下降的趨勢.2017年至2018年治理效果表現(xiàn)最佳，全年氨氮污染總量下降幅度超過70%，下一步可以進(jìn)一步關(guān)注模擬出現(xiàn)高值的4、5、7、9等子流域，做更精確的污染防治.

目前的參數(shù)移植，相似流域的選取主要是基于水文相似性，如汪銀龍等[1]研究選擇了在空間上相鄰板澗河與毫清河流域，陳磊等[23]研究大寧河流域時，選擇了流域內(nèi)部的子流域進(jìn)行對比，這是兩種根據(jù)空間鄰近性尋找水文相似流域的方法.犟河流域的模擬表明，無資料流域的水文模擬參數(shù)從空間上考慮，可以參考上級流域的參數(shù)，對其他缺乏資料區(qū)域的研究具有參考價值.