基于SWAT模型的巴勒更河流域降雨-徑流關(guān)系

李 凱，王永強(qiáng)，許繼軍，吳志俊,2，許 翔,2

(1.長江科學(xué)院 水資源綜合利用研究所，武漢 430010； 2.河海大學(xué) 水文水資源學(xué)院，南京 210098)

1 研究背景

徑流的主要來源是降雨，降雨的變化通常會(huì)直接導(dǎo)致徑流改變。自19世紀(jì)以來，全球范圍內(nèi)的氣溫普遍升高，加劇陸地和海洋水循環(huán)過程，同時(shí)導(dǎo)致區(qū)域降雨和空間分布格局發(fā)生變化。此外，人類活動(dòng)的不斷增強(qiáng)，一方面導(dǎo)致蒸發(fā)、入滲、產(chǎn)流和匯流等水文過程發(fā)生改變，另一方面導(dǎo)致流域降雨-徑流關(guān)系發(fā)生顯著性的改變[1]。為揭示巴勒更河流域降雨-徑流關(guān)系變化特征，研究其在時(shí)間上的演變機(jī)理，定量分離導(dǎo)致徑流變化的影響因素的貢獻(xiàn)率，對合理進(jìn)行流域水資源管理開發(fā)、掌握流域水循環(huán)機(jī)理等方面有著重要的理論價(jià)值和意義。

自1940年代以來，國內(nèi)外水文學(xué)家對降雨-徑流關(guān)系開展了大量研究[2-3]。通過蓄水容量曲線能夠有效研究降雨與徑流關(guān)系，該方法以嚴(yán)格的物理學(xué)概念為基礎(chǔ)，且與流域水平衡方程相吻合[4]。另有研究指出，降雨-徑流在時(shí)空上具有動(dòng)態(tài)和復(fù)雜的結(jié)構(gòu)，建立多變量廣義自回歸條件異方差模型同樣可以有效表示降雨-徑流的方差-協(xié)方差隨時(shí)間變化的關(guān)系[5]。李致家等[6]通過在全國15個(gè)典型流域建立地面徑流相關(guān)圖，分析區(qū)域降雨-徑流關(guān)系。陳利者等[7]在上述工作的基礎(chǔ)上做了進(jìn)一步研究，在不對地下水進(jìn)行分割情況下，首先構(gòu)建流域降雨-徑流相關(guān)圖，重新分析其降雨-徑流關(guān)系。近年來，降雨、徑流研究熱點(diǎn)主要集中在氣候變化和人類活動(dòng)影響下的降雨-徑流關(guān)系，分析方法主要有2種，即定性方法和定量方法。其中，徐東霞等[8]利用嫩江流域近50 a的資料定性分析降雨-徑流變化因素，得出的結(jié)論是徑流變化主要是受到降雨變化的影響。郭愛軍等[9]采用累積斜率變化率比較方法定量分離氣候變化與人類活動(dòng)對涇河徑流的影響。然而，對一些無徑流資料或缺少序列徑流資料的地區(qū)如何研究降雨-徑流關(guān)系是一個(gè)有待解決的難題[10]。目前，無資料地區(qū)徑流估算主要是利用水文比擬法、集總式水文模型、分布式水文模型開展水文模擬，水文比擬法估算精度并不高[11-12]。而集總式水文模型中部分參數(shù)仍然是經(jīng)驗(yàn)方法，其中大多數(shù)依據(jù)需要通過觀測到的降雨和徑流數(shù)據(jù)反演獲得。并且，集總式水文模型與實(shí)際降雨和徑流過程不匹配，因此通常難以獲得令人滿意的結(jié)果[13]。

SWAT模型是搭建在地理信息系統(tǒng)(GIS)上的半分布式水文模型，可以用來模擬各種不同的水文過程[14]，模型組成部分主要包括氣象、地表地下徑流、蒸散發(fā)等[15]。該模型通過分析數(shù)字高程、土地利用類型、土壤類型以及水文氣象數(shù)據(jù)建立主要的源匯特征，可以模擬流域內(nèi)的徑流變化特征。該模型已經(jīng)被應(yīng)用于大中小流域、干旱、半干旱地區(qū)，且都能夠得到令人滿意的模擬效果[16-19]。

本文的研究區(qū)域是位于青海德令哈縣內(nèi)的巴勒更河流域，首先研究水文模型能否適用于巴勒更河流域，然后利用該模型還原出流域1970—2015年的徑流過程。最后從周期、趨勢、突變等方面研究降雨-徑流關(guān)系，并進(jìn)一步探究氣候因素和人類活動(dòng)對該流域徑流變化的貢獻(xiàn)率。

2 研究區(qū)域概況

2.1 研究區(qū)域

巴勒更河流域位于青海海西內(nèi)蒙古自治州德令哈盆地(圖1)，位于96.6°E—96.83°E，37.33°N—37.6°N，海拔最高處為4 280 m，最低處為2 920 m。發(fā)源于宗務(wù)隆山巴勒更河，穿過懷頭他拉灌區(qū)，最后注入克魯克湖中，河流全長66.5 km。流域的集水面積882 km2，河道上下游比降1/66，多年年均流量0.75 m3/s，年徑流總量0.237億m3。該地區(qū)屬于高原山地氣候，多年平均氣溫4.2 ℃ ，相對濕度40%、多年平均風(fēng)速2.3 m/s、多年平均降雨量181.8 mm，多年最小降水量120 mm，降雨主要集中在5—8月份，約占全年降水量的60%。懷頭他拉水庫位于巴勒更河流域出口處，修建于1959年，多年平均入庫徑流量0.6 m3/s，水庫庫容890萬m3。

圖1 巴勒更河流域Fig.1 Map of Balegen River basin

2.2 數(shù)據(jù)來源與處理

本次研究需要輸入SWAT模型中的數(shù)據(jù)及其來源如表1所示，其中數(shù)字高程模型(DEM)和土地利用數(shù)據(jù)利用ArcGIS作投影變化、裁剪出所需研究區(qū)域；借助SPAW軟件計(jì)算所需的土壤參數(shù)；氣象數(shù)據(jù)來自德令哈、大柴旦、諾木洪3個(gè)氣象站；太陽輻射利用氣象站的日照數(shù)據(jù)進(jìn)行計(jì)算[20]。

表1 原始數(shù)據(jù)說明Table 1 Description of original data

3 SWAT模型的降雨-徑流模擬

3.1 SWAT模型原理

SWAT模型是Dile等[21]基于ArcGIS構(gòu)建和開發(fā)的半分布式水文模型，模型將流域細(xì)分為子流域，再根據(jù)土壤類型、土地利用類型以及坡度將子流域細(xì)分為水文響應(yīng)單元。模型需要輸入包括土地利用類型(圖2(a))、土壤類型(圖2(b))、坡度、水文氣象等數(shù)據(jù)。

圖2 土地利用類型和土壤類型Fig.2 Land use patterns and soil types

模型的主要原理是水量平衡原理(式(1))，它綜合考慮了流域內(nèi)所有的水文過程[22]。

(1)

式中：QWt、QW0分別表示土壤最終狀態(tài)和初始狀態(tài)土壤含水量(mm)；t表示天數(shù)(d)；Pi表示某天的降雨量(mm)；Q表示地表徑流深(mm)；ETi表示蒸散發(fā)量(mm)；Wi及QRi分別表示進(jìn)入包氣帶的水量(mm)和回歸流量(用深度表示，mm)；i表示某天。

(2)

式中：S為滯留系數(shù)(mm)；Ri為降雨深度(mm);Ia為降雨的初損(mm)。初損Ia與滯留系數(shù)S之間存在Ia=S，取值為0.2。當(dāng)Ri>0.2S時(shí)徑流出現(xiàn)，滯留系數(shù)S與土地利用、土壤類型、坡度有關(guān)，其關(guān)系式為

(3)

式中CN是與土壤類型、土地利用和土地管理?xiàng)l件相對應(yīng)的SCS曲線編號。

3.2 參數(shù)敏感性分析與校準(zhǔn)

為了提高模型校準(zhǔn)的效率，本文使用SUFI-2算法的Global Sensitivity Analysis模塊對流域出口處的23個(gè)參數(shù)進(jìn)行敏感性分析，刪除那些對模型模擬結(jié)果影響較弱的參數(shù)，最終得到12個(gè)有效參數(shù)，參數(shù)名稱見表2。SUFI-2算法是將采用拉丁超立方隨機(jī)采樣方法取得的參數(shù)值[23]，代入模型進(jìn)行模擬，計(jì)算得到目標(biāo)函數(shù)值。使用t-stat和p-value指標(biāo)值來評價(jià)參數(shù)是否足夠敏感，當(dāng)前者的絕對值越遠(yuǎn)離0時(shí)，則說明該參數(shù)在本次率定中對率定結(jié)果影響越顯著；當(dāng)p<0.05時(shí)表示該參數(shù)對結(jié)果的影響極為顯著。敏感性分析結(jié)果見表2。

表2 SUFI-2算法敏感性分析結(jié)果Table 2 Sensitivity analysis results of SUFI-2 algorithm

利用SWAT模型的配套軟件對模型進(jìn)行參數(shù)率定，模型將2006—2008年3 a時(shí)間設(shè)置為模型預(yù)熱期，2009—2013年作為模型參數(shù)校準(zhǔn)期，2014—2015年作為模型參數(shù)的驗(yàn)證期。執(zhí)行多組迭代，第一組迭代1 000次，其余每組迭代100次，每組迭代結(jié)束得到新的參數(shù)范圍，使參數(shù)范圍不斷縮小，直至得到最優(yōu)解。最優(yōu)值見表2。

3.3 降雨-徑流模擬方法與效果評價(jià)

通過使用納什效率系數(shù)NSE、R2、PBIAS對模型模擬結(jié)果進(jìn)行評價(jià)，詳細(xì)計(jì)算公式如下：

(4)

(5)

(6)

NSE一般用于驗(yàn)證模擬徑流系列與實(shí)測系列的過程符合度，取值范圍為[-∞，1]，值越靠近于1，表示模型模擬質(zhì)量越好。R2表示模擬值與實(shí)測值的相關(guān)性，系數(shù)值越接近1，表明實(shí)測與模擬的相關(guān)性越高。PBIAS偏差百分比衡量模擬數(shù)據(jù)的平均趨勢大于或小于觀測值，最佳情況下為0，低幅度值表示模擬效果很好，值>0時(shí)表示模型被低估，值<0表示模型被高估[24]。當(dāng)PBIAS<25%認(rèn)為模型模擬結(jié)果是可信賴的； PBIAS<15%時(shí)認(rèn)為模型能得到較好的模擬結(jié)果； PBIAS<10%時(shí)模型能夠得到令人滿意的結(jié)果。

3.4 參數(shù)校準(zhǔn)與模型驗(yàn)證結(jié)果

對模型中12個(gè)主要參數(shù)先利用2009—2013年徑流資料進(jìn)行率定，然后再利用2014—2015年的資料進(jìn)行驗(yàn)證，結(jié)果如圖3所示。在校準(zhǔn)期，模型徑流模擬曲線與實(shí)測的徑流曲線擬合效果較好；但在驗(yàn)證期，模型模擬的徑流曲線與實(shí)測的徑流曲線擬合效果較差。同時(shí)，從表3中模擬評價(jià)結(jié)果可知：校準(zhǔn)期內(nèi)R2和NSE均在0.76以上且PBIAS<3.8%，驗(yàn)證期R2和NSE分別只有0.58、0.51；但總體來說，模型模擬的月徑流過程符合R2和NSE>0.5、相對誤差<15%的要求，表明模型可以用于巴勒更河流域的徑流還原[25]。

圖3 實(shí)測與模擬徑流過程Fig.3 Measured and simulated runoff process

表3 模擬評價(jià)結(jié)果Table 3 Simulation evaluation results

3.5 徑流還原結(jié)果

將經(jīng)過校準(zhǔn)及驗(yàn)證后的12個(gè)參數(shù)返回到模型中，修改每個(gè)時(shí)期對應(yīng)的土地利用類型(表4)，并驅(qū)動(dòng)SWAT模型模擬相應(yīng)時(shí)期的逐月徑流過程，模擬結(jié)果如圖4所示。

表4 不同土地利用類型年份對應(yīng)徑流還原時(shí)期Table 4 Runoff restoration periods corresponding toyears of different land use patterns

圖4 徑流量還原過程Fig.4 Runoff restoration process

4 降雨-徑流關(guān)系及影響因素分析

為揭示巴勒更河流域降雨-徑流關(guān)系隨時(shí)的變化特征，研究其在時(shí)間上的演變機(jī)理，定量分析導(dǎo)致徑流變化的影響因素具有重要意義。本文用Kendall秩次檢驗(yàn)[26]和5點(diǎn)滑動(dòng)平均[27]對流域降雨、徑流趨勢展開分析；利用Mann-Kendall(M-K)突變檢驗(yàn)[28]識(shí)別降雨、徑流的突變位置;用小波分析[29]識(shí)別降雨-徑流的周期性;利用累積距平方法[30]判斷降雨-徑流關(guān)系的突變點(diǎn)。采用累積斜率變化率比較方法[31]定量分離出氣候變化和人類活動(dòng)對徑流量變化的貢獻(xiàn)率。

4.1 降雨、徑流年際變化特征分析

巴勒更河流域1970—2015年降雨、徑流變化過程及5點(diǎn)滑動(dòng)平均如圖5所示，降雨極大值發(fā)生在2009年，為307.39 mm，徑流極大值發(fā)生在1986年，為0.27×108m3/s，極大值出現(xiàn)年份不統(tǒng)一，且偏離較大。降雨和徑流出現(xiàn)極小值的時(shí)間大致相同，降雨極小值出現(xiàn)在1972年，徑流發(fā)生在1971年,其值分別為66.73 mm和0.11×108m3。降雨、徑流的趨勢變化也不相同，其中降雨具有明顯增加的趨勢，而徑流的變化趨勢并不明顯。由降雨、徑流的滑動(dòng)平均曲線同樣可以明顯看出，降雨、徑流波動(dòng)變化程度差別較大。同樣地，降雨、徑流發(fā)生突變的發(fā)生位置有差別，通過M-K檢驗(yàn)得到降雨的突變位置發(fā)生在1981年、徑流發(fā)生突變位置在1995年、2010年(圖6)。從圖7、圖8中可以明顯看出，降雨主要存在4、9、28 a的周期，其中以28 a的周期性變化最為顯著。在這個(gè)周期段內(nèi)，降雨經(jīng)歷“枯—豐—枯—豐—枯”5次豐枯交替，并且在未來的一段時(shí)間降雨仍然是偏枯的。徑流主要存在4、10、23 a的周期，其中23 a的周期最為明顯，在這個(gè)時(shí)期內(nèi)徑流經(jīng)歷了“枯—豐—枯—豐—枯”5次豐枯交替，同樣地，在未來一段時(shí)間內(nèi)，徑流也將是偏枯，這與降雨的變化是一致的。

圖5 降雨量和徑流量隨時(shí)間變化過程Fig.5 Time-history curves of rainfall and runoff

圖6 降雨量和徑流量的M-K檢驗(yàn)Fig.6 M-K test of rainfall and runoff series

圖7 小波分析系數(shù)實(shí)部等值線Fig.7 Wavelet coefficient real part of rainfall andrunoff series

圖8 降雨-徑流小波方差Fig.8 Wavelet variance of rainfall and runoff series

4.2 降雨-徑流關(guān)系變化分析

利用累積距平方法判斷降雨-徑流關(guān)系發(fā)生顯著性改變的年份，并得到1970—2015年降雨累積距平年際變化(圖9(a))，在1989年前后降雨累積距平分布表現(xiàn)出增大和減小的趨勢，顯然在該時(shí)段降雨發(fā)生突變。從1989—2015年降雨量累積距平的年際變化(圖9(b))可以看出，在2001年前后降雨量累積距平曲線的變化趨勢同樣也是先減小后增大。

計(jì)算得到1970—2015年徑流累積距平曲線，徑流累積曲線在1989年(圖9(c))、2001年(圖9(d))前后都表現(xiàn)出不一樣的趨勢，這與同期降雨量累積距平的變化趨勢是一樣的。通過徑流系數(shù)來進(jìn)一步判斷降雨-徑流關(guān)系變化情況，首先采用5點(diǎn)滑動(dòng)平均法分析徑流系數(shù)變化趨勢(圖10)。從圖10可知徑流系數(shù)隨著時(shí)間變化有明顯的減小趨勢。各個(gè)時(shí)期的徑流系數(shù)均值能很好地表現(xiàn)出徑流系數(shù)減少趨勢。1989—2000年相較于1970—1988年，徑流系數(shù)減少26.7%。推測徑流系數(shù)減小的原因可能是德令哈地區(qū)在1988年撤鎮(zhèn)建市，工業(yè)得到迅猛發(fā)展，城市取用水量明顯增大。2001—2015年相較于1989—2000年，徑流系數(shù)減少10%，這一時(shí)期相較于1989—2000年，徑流減少更為顯著，推測主要原因是在該時(shí)期懷頭他拉撤鄉(xiāng)建鎮(zhèn)，城鎮(zhèn)建立，該地區(qū)的需水量更大，河道內(nèi)取水量增大。

圖9 降雨、徑流累積距平變化特征Fig.9 Changes in the cumulative anomalies ofrainfall and runoff series

圖10 徑流系數(shù)變化Fig.10 Change of runoff coefficient

4.3 徑流變化對氣候和人類活動(dòng)的響應(yīng)

氣候、人類活動(dòng)以及下墊面都是影響徑流形成和變化的主要因素。本文利用累積斜率變化率比較方法計(jì)算影響研究區(qū)域徑流變化因素的貢獻(xiàn)率。此方法將氣候因素簡化為降雨，綜合考慮人類活動(dòng)的影響，并忽略了流域下墊面的影響[32]。

具體步驟為設(shè)累積降雨-年份線性關(guān)系式在拐點(diǎn)前后2個(gè)時(shí)期的斜率分別為SPa和SPb(mm)，累積徑流量-年份線性關(guān)系的斜率在拐點(diǎn)前后2個(gè)時(shí)期的徑流量為SRa和SRb(億m3)，則累積降雨斜率變化率為(SPb-SPa)/SPa；同樣地，累積徑流斜率變化率為(SRb-SRa)/SRa，那么氣候變化導(dǎo)致徑流量變化的貢獻(xiàn)率Cp(%)為

(7)

由此簡化得到人類活動(dòng)對徑流變化的貢獻(xiàn)率CH(%)，即

CH=100-Cp。

(8)

經(jīng)過計(jì)算得到巴勒更河流域的累積降雨、徑流斜率及其變化率如表5、表6所示。第1時(shí)期對應(yīng)1970—1988年、第2時(shí)期對應(yīng)1989—2000年、第3時(shí)期對應(yīng)2001—2015年。第2時(shí)期與第1時(shí)期相比，累積降雨量-年份線性關(guān)系式的斜率減小了-8.72 mm/a，變化率為140.19%(表5)。與同期相比，累積徑流-年份線性關(guān)系式的斜率減小了0.000 4億m3/a，變化率為400%(表6)，利用式(7)、式(8)計(jì)算得到第2時(shí)期比第1時(shí)期的氣候變化對徑流量變化的貢獻(xiàn)率為35%、人類活動(dòng)影響的貢獻(xiàn)率為65%(表7)。

表5 累積降雨斜率及其變化率Table 5 Cumulative rainfall slope and its rate of change

表6 累積徑流斜率及其變化率Table 6 Cumulative runoff slope and its rate of change

表7 氣候變化和人類活動(dòng)對降雨、徑流貢獻(xiàn)率Table 7 Contribution rates of climate change andhuman activities to rainfall and runoff

第3時(shí)期與第2時(shí)期相比，累積降雨-年份關(guān)系式的斜率減小-2.29 mm/a，變化率為91.6%(表5)。同期相比累積徑流-年份關(guān)系式的斜率減少-0.000 5億m3/a，變化率為180%(表6)。計(jì)算得到第3時(shí)期相對于第2時(shí)期氣候變化和人類活動(dòng)對徑流變化的貢獻(xiàn)率分別為51%、49%(表7)。

5 結(jié)論與討論

利用SWAT模型還原巴勒更河流域46 a徑流過程，并分析其降雨-徑流關(guān)系的變化及原因，結(jié)論如下。

(1) 將SWAT模型應(yīng)用于巴勒更河流域，參數(shù)校準(zhǔn)期(2009—2013年)的評價(jià)指標(biāo)R2=0.78、NSE=0.76，說明模型在校準(zhǔn)期模擬效果較好。在驗(yàn)證期(2014 — 2015年)R2=0.58、NSE=0.50，模型結(jié)果可信。進(jìn)而可利用模型還原1970—2008年的徑流數(shù)據(jù)。

(2)從趨勢、突變、周期以及相關(guān)關(guān)系角度分析巴勒更河的降雨序列和經(jīng)過SWAT模型延長以后的徑流序列，降雨具有顯著的上升趨勢，而徑流的上升趨勢不顯著；降雨-徑流的突變點(diǎn)上也存在差異，降雨突變發(fā)生在1981年，而徑流在1995年、2010年都發(fā)生突變；兩者周期性上存在著相同的周期，但主周期不同，降雨主周期為23 a，徑流則為28 a。通過降雨量、徑流量累積距平曲線可以看出，降雨量與徑流量在1989年和2001年降雨-徑流的關(guān)系發(fā)生兩次顯著性變化，變化前后徑流系數(shù)分別減小26.7%、10%。

(3)以1970—1988年為基準(zhǔn)期，1989—2000年降雨變化和人類活動(dòng)的相對貢獻(xiàn)率為35%和65%，在2001—2015年的相對貢獻(xiàn)率為19%和81%。以1989—2000年為基準(zhǔn)期，2001—2005年降雨量和人類活動(dòng)對徑流量的影響貢獻(xiàn)率分別為51%、49%。

(4)降雨-徑流關(guān)系在1989年、2001年發(fā)生了明顯的改變，主要是因?yàn)樵?988年德令哈建市、2001年懷頭他拉撤鄉(xiāng)建鎮(zhèn)都導(dǎo)致人類活動(dòng)增強(qiáng)從而影響徑流變化。