基于SUFI-2方法的SWAT模型在淮河息縣流域月徑流模擬中的應(yīng)用

程海洲，余麗華

(寧波弘泰水利信息科技有限公司，浙江 寧波 315192)

?

基于SUFI-2方法的SWAT模型在淮河息縣流域月徑流模擬中的應(yīng)用

程海洲，余麗華

(寧波弘泰水利信息科技有限公司，浙江 寧波 315192)

通過對模型輸入?yún)?shù)進行不確定性分析，確定模型關(guān)鍵性參數(shù)和不確定性來源，對于提高模型率定效率和改進模型結(jié)構(gòu)具有重要意義。通過應(yīng)用SWAT分布式水文模型對息縣流域1980—1995年的月徑流過程進行模擬，同時基于SUFI-2算法進行模型參數(shù)敏感性分析、率定驗證以及不確定性分析。結(jié)果表明：率定期與驗證期的相關(guān)性系數(shù)R2和模型效率系數(shù)ENS均高于0.7，月徑流量模擬值與實測流量過程線擬合程度較好；驗證期小于率定期，模型不確定性較??；SUFI-2算法能夠較好地分析SWAT模型參數(shù)的不確定性。

SUFI-2算法；SWAT模型；徑流模擬

近年來，分布式水文模型越來越多地用于為洪水頻發(fā)、土地利用變化和氣候變化以及環(huán)境污染嚴(yán)重地區(qū)提供科學(xué)的決策依據(jù)[1]。20世紀(jì)90年代，在美國農(nóng)業(yè)部(USDA)農(nóng)業(yè)研究中心(ARS)Jeff Arnold博士的主持下，在吸取CREAMS、GLEAMS、EPIC、SWRRB等模型優(yōu)點的基礎(chǔ)上，將SWRRB和ROTO整合為一個新的模型[2-3]，即SWAT。其后，該模型開始廣泛應(yīng)用于徑流模擬[4-7]、水質(zhì)模擬[7-9]、土地覆被變化的影響評估[10-11]、氣候變化的影響評估[11]等科學(xué)研究及實際應(yīng)用中。通過以上資料可知，基于物理過程的分布式水文模型SWAT能夠很好地分析出復(fù)雜大流域上土地管理措施對水、土/沙、營養(yǎng)物及農(nóng)化物等的影響。由于分布式水文模型在模擬水文過程時需使用到大量的參數(shù)，因此，如何對這些參數(shù)進行率定并進行合理的不確定性分析就顯得尤為重要。許多學(xué)者已針對這一問題采取過許多不同的研究方法[12-16]，本文使用由Abbaspour等人開發(fā)的SWAT-CUP中的SUFI-2算法對SWAT模型進行率定及不確定性分析。

1 研究方法

1.1 SWAT模型

SWAT模型是一個具有很強物理機制的分布式流域水文模型，在具有不同的土壤類型、不同的土地利用方式和管理條件下的復(fù)雜流域、資料缺乏地區(qū)，以及干燥、濕潤或介于兩者之間的氣候條件下，對多種不同的水文物理過程，都能取得令人滿意的模擬結(jié)果[17]。

SWAT模型是以研究流域的數(shù)字高程模型為基礎(chǔ)，提取坡度、坡向、坡長、流向、河網(wǎng)和流域邊界等流域特征信息，綜合利用土壤類型、土地植被類型、土地利用等的空間數(shù)據(jù)，將流域離散化為不同的子流域以至水文響應(yīng)單元(HRU)；再以HRU為空間計算單位，對流域的產(chǎn)流、坡面匯流、河道匯流等進行分析計算，并將運算結(jié)果匯集于流域出口。根據(jù)水文循環(huán)原理，SWAT模型的水量平衡基本表達式為：

(1)

式中：SWt為期末的土壤含水率(mm)；SW0為土壤初始含水率(mm)；t為時間(d)；Rday為第i天的降水量(mm)；Qsurf為第i天的地表徑流(mm)；Ea為第i天的蒸發(fā)量(mm)；Wseep為第i天側(cè)向滲流量(mm)；Qgw為第i天的垂向地下水出流量(mm)。

SWAT模型中，計算地表徑流有2種方法：SCS徑流曲線數(shù)法[18]、Green-Ampt模型法[19]，本文選取SCS曲線法計算地表徑流；計算潛在蒸散發(fā)有3種方法：Hargreaves[20]、Penman-Montieth[21]、Priestly-Taylor[22]，本文因缺少實測太陽輻射數(shù)據(jù)，而Hargreaves方法不需要太陽輻射數(shù)據(jù)，故選取該方法計算潛在蒸散發(fā)。

1.2 SUFI-2算法

SUFI-2[23](Sequential Uncertainty Fitting，ver.2)是1種參數(shù)估計最優(yōu)化方法，能夠考慮模型結(jié)構(gòu)、輸入數(shù)據(jù)、模型參數(shù)等因素對模型計算的不確定性影響。它采用拉丁超立方隨機采樣方法取得參數(shù)值，帶入模型運行摸擬，再計算目標(biāo)函數(shù)值，并用P因子和R因子表示模型的不確定性程度[24]。P因子用95PPU表示，表示模擬的數(shù)據(jù)包括了95%的不確定性，剔除了5%的極壞模擬情況；R因子表示95PPU的上下限的平均距離與標(biāo)準(zhǔn)偏差的取值；P因子的范圍為0～1，R因子的范圍為0～∞，當(dāng)P因子和R因子分別為1和0時表示模擬結(jié)果和實測數(shù)據(jù)一致。

在SUFI-2方法[25]中，敏感性矩陣和參數(shù)協(xié)方差矩陣的計算公式為：

(2)

(3)

參數(shù)bj的95%置信區(qū)間通過中的對角元素計算得到：

(4)

而SUFI-2方法采用的評估模型不確定性程度的P因子和R因子的計算公式為：

(5)

在SUFI-2算法計算中，通常先假設(shè)一個比較大的參數(shù)區(qū)間，使得實測數(shù)據(jù)能夠被包含在95PPU范圍內(nèi)，然后再結(jié)合模型運行結(jié)果逐步縮小區(qū)間范圍，每一次參數(shù)范圍的改變，都重新計算敏感性矩陣和協(xié)方差矩陣，并更新參數(shù)范圍，進行下一步模擬，循環(huán)往復(fù)直至模擬值與實測值基本達到一致。根據(jù)SUFI-2算法，模型結(jié)果由一組模型參數(shù)決定，因此能最大程度地反映模型輸入?yún)?shù)的不確定性。

2 應(yīng)用研究

2.1 研究區(qū)概況

淮河流域地處我國東部，介于長江和黃河兩流域之間，位于東經(jīng)111°55′～121°25′，北緯30°55′～36°36′，面積為27萬km2。息縣位于淮河上游、河南省東南部、信陽市東北部，地形以低平的平原和緩丘為主，呈西北向東南略為傾斜，平均海拔47.00 m?；春訖M貫全境，境內(nèi)流長75.5 km，面積為920.97 km2，最低高程為39.00 m，最高高程為957.00 m，平均高程為141.30 m。氣候?qū)賮啛釒蚺瘻貛н^渡形氣候，年平均氣溫15.2 ℃，年平均降雨量為946 mm，全年無霜期200 d左右。流域內(nèi)主要土壤類型為水稻土和黃褐土，另外還有粗骨土、黃棕壤土和灰潮土等。流域土地利用/土地覆被的空間分布具有明顯的地形和地域的差異。其主要類型以農(nóng)田為主，占流域面積的46.8%；其次是林地和稻米，分別占18.95%和18.79%。息縣素有天下第一縣之稱，是我國重要的農(nóng)業(yè)生產(chǎn)基地。近年來，由于氣候變化及農(nóng)業(yè)耕作方式的改變，使得該地區(qū)的水循環(huán)特征發(fā)生了明顯的改變。

2.2 基礎(chǔ)資料庫建立

SWAT模型所需的數(shù)據(jù)有DEM、土壤類型、土地利用、氣象數(shù)據(jù)、降水?dāng)?shù)據(jù)、蒸發(fā)數(shù)據(jù)、控制站點流量資料等，相關(guān)數(shù)據(jù)的參數(shù)類別及數(shù)據(jù)來源見表1。

表1 SWAT模型輸入數(shù)據(jù)表

在數(shù)據(jù)準(zhǔn)備階段，空間數(shù)據(jù)必須具有相同的投影格式和坐標(biāo)系統(tǒng)，同時需要定義空間數(shù)據(jù)的投影。基于研究面積和區(qū)域，本文采用GCS_WGS_1984坐標(biāo)系，Transverse_Mercator投影，中央子午線為111°。其次，本文應(yīng)用ArcGIS的Archydro Tools模塊對鄞江流域的DEM數(shù)據(jù)(見圖1a)進行分析處理，包括洼地填充、流向分析、流量累積網(wǎng)格計算，并通過設(shè)定集水面積閥值自動提取數(shù)字河網(wǎng)，劃分集水區(qū)域，獲得44個子流域(見圖1b)。根據(jù)不同土地利用(見圖1c)和土壤類型(見圖1d)的組合，在每一個子流域內(nèi)進一步劃分水文響應(yīng)單元HRUs，劃分為309個HRUs。

(a)DEM圖

(b)子流域劃分圖

(c)1∶20萬土壤類型圖

(d)1∶20萬土地利用圖圖1 研究區(qū)域的數(shù)字地圖

3 結(jié)果分析

3.1 參數(shù)敏感性分析

SWAT模型參數(shù)眾多，有的參數(shù)對模型模擬結(jié)果影響較小，有的參數(shù)細(xì)微變化對模型模擬結(jié)果有舉足輕重的作用。因此，在采用SUFI-2方法進行參數(shù)不確定性分析之前，首先需要選出敏感性較強的模型參數(shù)，并對參數(shù)范圍進行校正。本文使用SWAT模型自帶的拉丁超立方隨機抽樣LH-OAT方法，對與徑流有關(guān)的26個參數(shù)進行敏感性分析，表2列出等級前10的參數(shù)，以及這些參數(shù)的物理意義及初始范圍。

表2 參數(shù)敏感性排名表

在模型率定的過程中僅考慮這些敏感的參數(shù)，其余的參數(shù)對徑流的模擬并沒有顯著的影響，故其改變不會對模擬效果產(chǎn)生顯著的改變。

3.2 不確定性分析

根據(jù)選定的參數(shù)范圍，利用SUFI-2方法對研究區(qū)域模型輸入?yún)?shù)進行不確定性分析，模型目標(biāo)函數(shù)采用觀測值和模擬值間的相關(guān)系數(shù)(R2)和Nash-Sutcliffe效率系數(shù)(NSE)，以便綜合評價SWAT模型的模擬效果：

(6)

(7)

本文選取1980年作為預(yù)熱年，1981—1990年為率定期，1991—1995年作為驗證期；先對模型進行人工率定，使得模擬值與觀測值間的差值減小，之后再使用SUFI-2算法對模型進行率定及不確定性分析。其中，人工調(diào)整的參數(shù)主要有基流系數(shù)(Alpha-Bf)、土壤蒸發(fā)補償系數(shù)(Esco)、徑流曲線系數(shù)(CN)、土壤可利用水量(Sol-AWC)。人工調(diào)參雖然很費時，但卻有助于提高自動率定的精度。

3.2.1 率定期不確定性分析

文獻[1]模擬12 a日徑流數(shù)據(jù)時，設(shè)置SUFI-2算法模擬次數(shù)為2 000次，模擬效果較好。本文模擬11 a月徑流數(shù)據(jù)，模擬次數(shù)設(shè)置為1 000次，通過將模擬數(shù)據(jù)與實測數(shù)據(jù)進行對比發(fā)現(xiàn)，模擬結(jié)果具有較高的相關(guān)系數(shù)(R2)和Nash-Sutcliffe系數(shù)(NSE)。同時，根據(jù)Moriasi[26]提出的模型計算評價標(biāo)準(zhǔn)：若R2>0.7，NSE>0.5，就認(rèn)為模擬效果滿意。本文的模擬結(jié)果滿足要求，表3列出了率定期月徑流模擬的評價。圖2列出了率定期月徑流量實測值和模擬值。

表3 率定期月徑流模擬評價因子表

圖2 率定期月徑流量實測值和模擬值圖

由表3看出，雖然率定期R2為0.85，NSE為0.82，r-因子0.67，說明模擬結(jié)果令人滿意，但P-因子只達到0.28。P-因子是95%預(yù)測不確定性(95PPU)范圍內(nèi)實測數(shù)據(jù)的百分比，0.28說明僅有28%的實測數(shù)據(jù)落在95PPU范圍內(nèi)，說明模型的不確定程度較高。原因可能是輸入數(shù)據(jù)如降雨、溫度等的不準(zhǔn)確，上游水庫放水等。

表4列出了結(jié)合人工率定及用SUFI-2算法率定后的敏感性參數(shù)最適值及范圍，這些參數(shù)值在SWAT模型中已加以修改，以便于進行下一步的模型驗證及對流域進行其他方面的研究。

表4 SWAT徑流敏感參數(shù)及最適值表

3.2.2 驗證期不確定性分析

由表3看出，驗證期1991—1995年R2為0.73，NSE為0.71，R-因子0.65，由于率定期P-因子較低，故驗證期P-因子0.26也較低。圖3列出了驗證期徑流量實測值和模擬值，可看出模擬效果較好，由驗證期率定得的參數(shù)能夠很好地代表流域及其產(chǎn)匯流的特征，故而模型能夠在今后對流域的研究中起到重要作用。

圖3 驗證期月徑流量實測值和模擬值圖

4 結(jié) 論

本文利用淮河息縣流域的DEM、土壤類型、土地利用、水文氣象等資料，應(yīng)用SWAT模型對該流域的1980—1995年16 a的月徑流過程進行模擬，并采用SUFI-2算法對模型參數(shù)進行率定驗證和不確定性分析，主要結(jié)論如下：

(1)應(yīng)用SWAT模型進行的月徑流模擬在率定期和驗證期都能取得比較滿意的模擬結(jié)果，其中，月徑流量率定期NES為0.82，R2為0.85，驗證期NES為0.71，R2為0.73，表明SWAT模型在研究區(qū)域具有較好的適用性。

(2)參數(shù)敏感性分析結(jié)果表明，與徑流有關(guān)的26個參數(shù)中，Alpha-BF、GWQMN、ESCO、CN2、CANMX等對徑流模擬結(jié)果影響最為顯著。

(3)不確定性分析結(jié)果顯示，率定期和驗證期的P-因子分別為0.28和0.26，表明SWAT模型在研究區(qū)域的月徑流模擬中具有較高的不確定性。

(4)采用SUFI-2方法得到的參數(shù)組，可以取得較好的模擬效果，表現(xiàn)在較高的NSE值和值。

(5)率定期及驗證期徑流模擬結(jié)果曲線圖顯示出良好的模擬效果，驗證了SUFI-2方法可以較好的應(yīng)用于SWAT摸型不確定性分析之中。

(6)綜合考慮率定期和驗證期的模擬結(jié)果，可知在人工調(diào)參的基礎(chǔ)上對SWAT模型采用SWATCUP中SUFI-2算法進行率定能達到較高的模擬精度。

[1]Shimeils G.Setegn，Ragahavan Srinivasan，Bijian Dargahi.Hydrological modelling in the Lake Tana Basin，Ethiopia using SWAT model[J].The Open Hydrology Journal，2008(2):49-62.

[2]張銀輝.SWAT模型及其應(yīng)用研究進展[J].地理科學(xué)進展，2005，24(5):121-130.

[3]王中根,劉昌明,黃友波.SWAT模型的原理、結(jié)構(gòu)及應(yīng)用研究[J].地理科學(xué)進展，2003，22(1):79-86.

[4]Arnold J G and Allen P M.Estimating hydrologic budgets for three Illinois watersheds[J].Journal of Hydrology，1996，176(1/2/3/4):57-77.

[5]Arnold J G，Srinivasan R，Muttiah R S，Allen P M.Continental scale simulation of the hydrologic balance[J].Journal of the American Water Resources Association，1999，35(5):1037-1051.

[6]Eckhardt K and Arnold J G.Automatic calibration of a distributed catchment model[J].Journal of Hydrology，2001，251(1/2):103-109.

[7]Vassilios Pisinaras，Christos Petalas，Georgios D.Gikas，et al.Hyrologyical and water quality modeling in a medium-sized basin using the Soil and Water Assessment Tool[J].Desalination，2010，250(1):274-286.

[8]Debele B，Srinivasan R，Yves Parlange J.Coupling upland watershed and downstream waterbody hydrodynamic and water quality models(SWAT and CE-QUAL-W2)for better water resources management in complex river basins[J].Environ Model Assess，2008，13(1):135-153.

[9]Abbaspour KC，Yang J，Maximov I，et al.Modelling hydrology and water quality in the pre-alpine/alpine Thur watershed using SWAT[J].Journal of Hydrology，2007，333(2/3/4):413-430.

[10]胥彥玲,王蘇艦,李懷恩.土地覆被變化對流域非點源污染的影響[J].水土保持研究，2010，17(3):250-253.

[11]張圣微,雷玉平,姚琴.等.土地覆被和氣候變化對拉薩河流域徑流量的影響[J].水資源保護，2010，26(2):39-44.

[12]C.H.Green and A.van Griensven.Autocalibration in hydrologic modeling:Using SWAT2005 in small-scale watersheds[J].Environmental Modelling & Software，2008，23(4):422-434.

[13]A.van Griensven，T Meixner，S.Grunwald，et al.A global sensitivity analysis tool for the parameters of multi-variable catchment models[J].Journal of Hydrology，2006，324(1/2/3/4):10-23.

[14]Zhang Li，Zongxue Xu，Quanxi Shao，et al.Parameter estimation and uncertainty analysis of SWAT model in upper reaches of the Heihe river basin[J].Hydrological Processes，2009，23(19):2744-2753.

[15]秦富倉,張麗娟,余新曉,等.SWAT模型自動校準(zhǔn)模塊在云州水庫流域參數(shù)率定研究[J].水土保持研究，2010，17(2):86-89.

[16]朱麗,秦富倉,姚云峰,等.SWAT模型靈敏性分析模塊在中尺度流域的應(yīng)用[J].水土保持研究，2011，18(1):161-165.

[17]Vanliew M W and Garbrecht J.Hydrologic simulation of the Little Washita River Experimental Watershed using SWAT[J].Journal of the American water resource association，2003，39(2):413-426.

[18]USDA.Urban Hydrology for small watersheds[R].Engineering Division，Soil Conversation Service，USDA，Technical Release 55，1986.

[19]Green W H ，Ampt G A.Studies on Soil Physics，1.The Flow of Air and Water through Soils[J].Journal of Agricultural Sciences，1981(4)：1-24.

[20]Hargreaves，G.L.；Hargreaves，G.H.；Riley，J.P.:Agricultural benefits for Senegal River Basin[J].Journal of Irrigation and Drainage Engineering，1985，111(2):113-124.

[21]Penman．H．L．Evaporation：An introductory survey[J]．Netherlands Journal of Agricultural Science，1956(4)：7-29．

[22]Priestley C.H.B.and Taylor R.J.On the assessment of surface heat flux and evaporation using large-scale parameters[J].Month.Wea.Rev，1972，(100)：81-92.

[23]Abbaspour KC，Johnson CA，van Genuchten MT.Estimating uncertain flow and transport parameters using a sequential uncertainty fitting procedure[J].Vadose Zone Journal，2004，3(4):1340-1352.

[24]魏丹,劉智勇,李小冰.SWAT模型及SUFI-2算法在禿尾河上游流域徑流模擬中的應(yīng)用[J].干旱地區(qū)農(nóng)業(yè)研究，2013，30(6):200-206.

[25]Gupta，H.V’,S.Sorooshian,P.O.Yapo.Toward improved calibration of hydrologic models:Multiple and noncommensurable measures of information[J].Water Resources Research，1998，34(4):751-763.

[26]Moriasi DN，Arnold JG，Van Liew MW，Binger RL，Harmel RD,Veith T.Model evaluation guidelines for systematic quantification of accuracy in watershed simulations[J].Transactions of the ASABE，2007，50(3):885-900.

(責(zé)任編輯 姚小槐)

2015-06-06

程海洲(1983-)，男，工程師，碩士，主要從事水文水資源模型及水利信息化研究。

E-mail:chenghzhyun@126.com

TV121

B

1008-701X(2016)01-0061-05

10.13641/j.cnki.33-1162/tv.2016.01.019