基于SWAT模型和SUFI-2算法的黃柏河?xùn)|支流域徑流模擬研究

萬 浩，董曉華，彭 濤，劉 冀，喻 丹，薄會娟，陳 亮

(1.三峽大學(xué)水利與環(huán)境學(xué)院，湖北 宜昌 443002；2.水資源安全保障湖北省協(xié)同創(chuàng)新中心，武漢 430072)

0 引 言

黃柏河流域位于湖北省宜昌市西北部，是宜昌市重要的水源地，也是著名的磷礦產(chǎn)區(qū)。近年來隨著經(jīng)濟的不斷發(fā)展和磷礦的大量開采，流域生態(tài)系統(tǒng)的結(jié)構(gòu)和功能遭到嚴重損害，引發(fā)了一系列的生態(tài)與環(huán)境問題，例如流域水環(huán)境惡化、水土流失以及河道斷流等現(xiàn)象，水資源保護形勢日趨嚴峻[1,2]。水作為污染物傳輸和轉(zhuǎn)化的基本載體，是維持生態(tài)平衡的物質(zhì)基礎(chǔ)，因此準確模擬徑流的變化規(guī)律是流域水資源管理和規(guī)劃的重要基礎(chǔ)[3]。流域水文模型是研究水文循環(huán)過程的重要工具?，F(xiàn)有的水文模型可以分為集總式水文模型和分布式水文模型，其中分布式水文模型能夠準確模擬下墊面和降雨的空間差異性，被廣泛用于流域水文過程模擬的相關(guān)研究中[4]。SWAT(Soil and Water Assessment Tool)模型是美國農(nóng)業(yè)部(USDA)農(nóng)業(yè)研究中心(ARS)于20世紀90年代開發(fā)并不斷完善的一個流域尺度的分布式水文模型。SWAT模型具有較強的物理基礎(chǔ)，能夠模擬流域氣候、土壤條件、土地利用以及管理措施變化下的水沙、營養(yǎng)物、殺蟲劑等物質(zhì)的運移[5]。在國內(nèi)外，SWAT模型已得到了廣泛應(yīng)用[6,7]。

SWAT模型參數(shù)眾多，很多參數(shù)難以直接獲取，且同時需要考慮不同子流域參數(shù)的空間差異性，參數(shù)率定是一項繁雜的工作。此外模型是將復(fù)雜的水文現(xiàn)象和過程經(jīng)概化后的近似的數(shù)學(xué)物理模型，其中存在著模型結(jié)構(gòu)、模型輸入和模型參數(shù)這三個方面的不確定性。因此，采用合理的參數(shù)率定方法對于提高模型模擬精度和模型的可靠性格外重要。序貫不確定擬合法(Sequential Uncertainty Fitting ver2，SUFI-2)作為近年來使用較多的不確定性分析方法，可以方便快捷地對SWAT模型進行參數(shù)率定，驗證和不確定性分析，對于模擬精度和運算效率都有很大的提高。

以往對于黃柏河的研究，都主要集中在流域的水質(zhì)評價以及污染成因分析，采用分布式水文模型來研究流域的徑流過程的研究較少。本文構(gòu)建黃柏河?xùn)|支流域SWAT模型，進行流域的月徑流模擬，同時運用SUFI-2算法進行模型的率定、驗證和不確定性分析，以期為黃柏河流域的水資源高效管理規(guī)劃提供科學(xué)依據(jù)。

1 數(shù)據(jù)與方法

1.1 研究區(qū)概況

黃柏河位于宜昌市西北部，是長江葛洲壩庫區(qū)最大的一條支流，干流全長162 km。黃柏河在夷陵區(qū)兩河口以上分為東西兩支，東支發(fā)源于宜昌市夷陵區(qū)黑良山，東支干流長130 km，集雨面積1 165 km2，河道平均坡降0.6%。流域地勢為北高南低，北部海拔1 300～1 900 m，中部海拔1 000～1 200 m，南部海拔200～800 m。黃柏河流域?qū)儆趤啛釒Ъ撅L(fēng)性濕潤氣候，四季分明，雨量豐沛，多年平均氣溫16.9 ℃，多年平均降水量1 138 mm，且降水分布具有明顯的季節(jié)性，降水主要集中在6-9月，暴雨、短時強降水等災(zāi)害性天氣多發(fā)。

黃柏河流域水能資源豐富，多年平均徑流量約8.95 億m3，由上游到下游建有玄廟觀、天福廟、西北口、尚家河水庫，流域供水覆蓋了宜昌市50%的人口，是宜昌市的重要水源地。本研究選擇黃柏河?xùn)|支尚家河水庫以上流域作為研究區(qū)(圖 1)，研究區(qū)面積為932.26 km2。

圖1 研究區(qū)概況Fig.1 Overview of the study area

1.2 數(shù)據(jù)來源及處理

SWAT模型的輸入數(shù)據(jù)主要包括：數(shù)字高程DEM、土地利用數(shù)據(jù)、土壤數(shù)據(jù)以及氣象數(shù)據(jù)。其中DEM數(shù)據(jù)、土地利用和土壤柵格數(shù)據(jù)需要采用相同的投影坐標(biāo)系統(tǒng)。根據(jù)研究區(qū)的地理位置，本文采用WGS-1984地理坐標(biāo)系統(tǒng)和6°分帶的UTM(通用橫軸墨卡托)投影坐標(biāo)，帶號為49N，中央經(jīng)線為111°E。

本研究所用的DEM高程數(shù)據(jù)來源于地理空間數(shù)據(jù)云(http://www.gscloud.cn/)提供的SRTM數(shù)據(jù)，分辨率為90M，研究區(qū)的DEM圖如圖 1所示。土地利用數(shù)據(jù)是基于Landsat8遙感影像，數(shù)據(jù)來源于地理空間數(shù)據(jù)云(http://www.gscloud.cn/)，使用ENVI5.1軟件采用K-Means非監(jiān)督分類[8]的方法進行遙感解譯得到，如圖 2(a)所示。根據(jù)SWAT模型數(shù)據(jù)庫的要求，進行土地利用編碼，各土地利用類型與SWAT模型編碼的對應(yīng)關(guān)系如表 1所示。本研究的土壤數(shù)據(jù)來源于聯(lián)合國糧農(nóng)組織(Food and Agriculture Organization-FAO)和維也納國際應(yīng)用系統(tǒng)研究所(International Institute for Applied Systems Analysis-IIASA)構(gòu)建的世界協(xié)調(diào)土壤數(shù)據(jù)庫HWSD(Harmonized World Soil Database)。本文采用魏懷斌等[9]提出的方法來構(gòu)建研究區(qū)的土壤數(shù)據(jù)庫。土壤類型的采用FAO-90土壤分類系統(tǒng)，研究區(qū)共有5種土壤類型，圖 2(b)顯示的是研究區(qū)的土壤空間分布情況，各土壤類型的面積及對應(yīng)編碼見表1。SWAT需要輸入日尺度的氣溫(日平均、最高和最低)、太陽輻射，風(fēng)速、相對濕度、降水、天氣發(fā)生器等數(shù)據(jù)來驅(qū)動模型。本研究采用美國國家氣象預(yù)報中心(National Centers for Atmospheric Prediction-NCEP)的全球氣候再分析系統(tǒng)(Climate Forecast System Reanalysis-CFSR)提供的全球高精度氣候變量數(shù)據(jù)來構(gòu)建研究區(qū)的天氣發(fā)生器。

圖2 研究區(qū)土地利用圖及土壤類型圖Fig.2 Land use map categories and soil map of the study area

序號土地利用類型SWAT編碼面積/km2比例/%FAO90土壤名稱面積/km2比例/%1荒地BARR1.010.10ATc人為土41.521.772耕地AGRL142.1715.25LVh1高活性淋溶土1163.8217.573水域WATR12.381.33LVh2高活性淋溶土2439.7347.174林地FRST758.9381.41LVh3高活性淋溶土3270.7229.035城鎮(zhèn)URBN17.771.91ALh高活性強酸土16.474.46 合計 932.26100合計932.26100

降雨采用位于流域內(nèi)4座水庫大壩處(圖 1)雨量站的2008-2016年實測日降雨數(shù)據(jù)。因缺乏實測數(shù)據(jù)，其他氣象數(shù)據(jù)則采用構(gòu)建的天氣發(fā)生器模擬。各水庫的日平均出入庫流量來源于黃柏河流域管理局提供的逐日觀測數(shù)據(jù)。

1.3 研究方法

1.3.1 SWAT模型

SWAT是分布式的物理-概念水文模型，采用分布式計算。SWAT模型的原理參見SWAT官方網(wǎng)站的相關(guān)技術(shù)文(http://swat.tamu.edu /documentation)，此處不再贅述。

為了準確模擬研究區(qū)的2008-2016年的月徑流過程，本文基于研究區(qū)的DEM、土地利用、土壤和日尺度的氣象數(shù)據(jù)建立研究區(qū)的SWAT模型。根據(jù)DEM數(shù)據(jù)，并結(jié)合流域?qū)嶋H情況和模型的計算效率，設(shè)定集水面積劃分閾值為2 000 hm2，并手動將4個水庫出口添加為子流域出口，最終將研究區(qū)域劃分為29個子流域，如圖1所示。土地利用類型、土壤類型的劃分閾值為5%，坡度類型劃分閾值為10%，最終生成315個水文響應(yīng)單元。SWAT模型的地表徑流模擬采用SCS曲線法，河道匯流演算采用變動存儲系數(shù)法，潛在蒸散發(fā)采用Penman-Monteith方法。

1.3.2 SWAT-CUP和SUFI-2算法

SWAT-CUP(SWAT Calibration and Uncertainty Programs)是Abbaspour等開發(fā)的一款可實現(xiàn)對SWAT模型進行參數(shù)敏感性分析、率定、驗證和不確定性分析的計算機程序[10]。該程序包含了GLUE、ParaSol、SUFI2、MCMC和PSO五種優(yōu)化算法，其中SUFI-2算法因為運行效率高，模擬效果好等優(yōu)點，應(yīng)用較廣泛[11]。本文采用SUFI-2算法對所構(gòu)建的SWAT模型進行參數(shù)敏感性分析、率定、驗證和不確定性分析。

SUFI-2算法了考慮了水文模型的輸入數(shù)據(jù)、模型結(jié)構(gòu)和模型參數(shù)的不確定性，并用最后率定的參數(shù)分布范圍來反映這些不確定性[12]。算法采用P-factor和R-factor兩個指標(biāo)來表征不確定性，P-factor表示在95PPU區(qū)間內(nèi)包含觀測數(shù)據(jù)的百分比，即95%置信區(qū)間包含的觀測數(shù)據(jù)量。R-factor等于95PPU帶的平均寬度除以觀測值的標(biāo)準偏差，反映95%置信區(qū)間內(nèi)觀測值的密集程度[13]。本研究選取2.5%～97.5%的區(qū)間作為SUFI-2算法的95%的置信區(qū)間(95PPU)，兩個不確定性指標(biāo)(P-factor、R-factor)計算方法參見文獻[13]。理論上，P-factor的范圍是0～1，R-factor的范圍是0～∞。理想的情況是用最小的95PPU帶去包含大部分觀測數(shù)據(jù)(即P-factor等于1，R-factor接近于零)，但是通常一個較大的P-factor也對應(yīng)著較大的R-factor[14]。因此，在率定的過程中，要綜合權(quán)衡兩個指標(biāo)，使P-factor和R-factor均達到可接受的范圍。

SUFI-2算法中提供了兩種參數(shù)敏感性分析方法，One-At-a-Time(OAT)和全局敏感性(global sensitivity)方法。本研究采用全局敏感性分析對參數(shù)進行敏感性分析，敏感性計算公式見文獻[13]。在分析中，采用t-Stat和p-Values來表示參數(shù)的敏感性。其中，t-Stat 表示參數(shù)的敏感性，p-Values表示敏感性的顯著性。即t-Stat的絕對值越大和p-Values值越小，表示參數(shù)越敏感[15]。

2 結(jié)果與分析

2.1 敏感性分析結(jié)果

根據(jù)模型運算原理及各參數(shù)意義，參考Abbaspour等[16]的參數(shù)選擇和總結(jié)，結(jié)合其他研究者的成果，選取與徑流相關(guān)的19個參數(shù)進行全局敏感性分析。參數(shù)的原始范圍是根據(jù)經(jīng)驗和參數(shù)的物理意義確定的，見表2。程序運行500次，得出敏感性分析結(jié)果見表2。

表2 參數(shù)敏感性分析及率定結(jié)果Tab.2 Sensitivity analysis and calibrated values for model parameters

注：v_是指絕對變化，直接替代原參數(shù)值；r_是相對變化，用原始值乘以(1十變化值)。

從表2中可以看出對黃柏河?xùn)|支流域月徑流模擬結(jié)果最為敏感的參數(shù)依次為：水分條件II時的初始SCS徑流曲線數(shù)(CN2.mgt)、土壤蒸發(fā)補償系數(shù)(ESCO.hru)、土壤飽和水力傳導(dǎo)度(SOL_K.sol)等。CN2是采用SCS曲線法計算地表產(chǎn)流的關(guān)鍵參數(shù)，其值直接影響著徑流量的大小。與其他濕潤地區(qū)的研究[17,18]相比，相同的敏感性參數(shù)有CN2、ESCO、SOL-K、GWQMN，說明這些參數(shù)對于濕潤地區(qū)的徑流模擬有重要影響。

2.2 參數(shù)率定和驗證結(jié)果

選取2008年為模型的預(yù)熱期，以降低模型在運行初期初始條件的影響。2009-2013年為模型率定期，2014-2016年為驗證期，進行黃柏河?xùn)|支流域的月徑流模擬。選取國內(nèi)外文獻[19,20]中普遍使用的Nash-Sutcliffe Efficiency納什系數(shù)(NSE)、均方根誤差與標(biāo)準差比值(ratio of root mean square error to standard deviation，RSR)和百分比偏差(Percent Bias，PBLAS)來綜合評價模型的模擬精度，各評價指標(biāo)的計算方法見文獻[20]。一般認為NSE>0.5，RSR<0.7，|PBLAS|<25%時，模型模擬結(jié)果是可以接受的。

SUFI-2算法是通過不斷迭代來獲得最佳參數(shù)范圍，每次迭代都會計算目標(biāo)函數(shù)和參數(shù)之間的相關(guān)性，然后為下一次迭代提出新的參數(shù)范圍。需要注意的是，在迭代過程中，一些建議的參數(shù)范圍超出了參數(shù)的物理意義，需要對這些參數(shù)進行手動調(diào)整，使它們不超過最大/最小絕對值的范圍。本研究共進行了3次迭代，每次迭代運行500次，得出了最終的結(jié)果，參數(shù)率定結(jié)果見表2。

月徑流率定和驗證結(jié)果見圖3和表3。由圖3可以看出，玄廟觀、天福廟、西北口和尚家河站的SWAT模型模擬的月徑流過程線與實測月徑流過程線較為吻合，對于洪峰均有響應(yīng)。各站的均方根誤差與標(biāo)準差比值RSR均小于0.7，在誤差范圍內(nèi)。納什效率系數(shù)NSE除玄廟觀站為0.64外，其余各站均在0.75以上，西北口和尚家河站都達到了0.9以上，模擬精度較高。在驗證期，各站點的納什效率系數(shù)NSE均在0.85以上，相對于率定期均有所提高，這可能是由于驗證期相對于率定期時間較短，參數(shù)擬合相對較容易。同時發(fā)現(xiàn)，在率定期和驗證期下游站點的模擬精度都高于上游站點，上游站點所控制的流域面積較小，徑流受人類活動的影響較大。最上游的玄廟觀站在率定期和驗證期納什效率系數(shù)NSE均最小，分析原因可能是由于玄廟觀站受流域內(nèi)磷礦開采影響較大，流域水循環(huán)模式改變顯著，一部分地表徑流進入磷礦開采所形成的礦坑，這增加了模型模擬的難度。另外，根據(jù)調(diào)查玄廟觀流域內(nèi)建有小水電站，模型未能考慮小水電調(diào)洪等人類活動的影響，對于模型模擬精度也有影響。

表3 黃柏河?xùn)|支流域月徑流模擬結(jié)果統(tǒng)計Tab.3 Evaluation of monthly runoff simulation results

圖3 黃柏河流域各站月徑流觀測值和模擬值對比及95PPUFig.3 Observed,simulated and 95PPU of monthly runoff in calibration period

對于百分比偏差PBIAS而言，大部分站點在率定期和驗證期的PBIAS值均大于0且小于25%，說明模型在率定期和驗證期月徑流模擬結(jié)果較觀測值偏低。整體來看，3個評價指標(biāo)均在誤差范圍內(nèi)，SWAT模型能夠較好地模擬研究區(qū)的月徑流過程，在流域的月徑流模擬中具有較好的適應(yīng)性。

2.3 不確定性分析結(jié)果

不確定性結(jié)果見表3。在驗證期和率定期各站點的P-factor均大于0.65，R-factor均小于1、表明大多數(shù)觀測值落在95PPU區(qū)間內(nèi)，所構(gòu)建的SWAT模型徑流模擬的不確定性較小。

除玄廟觀站外，驗證期和率定期相比，各站點的不確定性均有所降低，觀測值落在95PPU區(qū)間內(nèi)的個數(shù)增多。在率定期，天福廟的P-factor最小，徑流模擬的不確定性最大，但是天福廟站的納什效率系數(shù)NSE并不是最小的，這表明模型的不確定性結(jié)果并不完全和模型模擬精度相一致。較大的參數(shù)范圍能夠反映出參數(shù)對模擬結(jié)果的影響，但由此產(chǎn)生較寬的不確定性區(qū)間，降低了模擬的置信水平。在參數(shù)率定的過程中，綜合考慮了模型的模擬精度和模型率定的不確定性，來提高模型模擬的可靠性。

3 討 論

由表3可知，玄廟觀、天福廟兩個水庫壩址處的徑流模擬結(jié)果比另外兩個水庫壩址處的模擬結(jié)果要差。原因可能是因為人類大量的采礦活動改變了流域下墊面條件。為了進一步分析這兩個水庫壩址處徑流模擬結(jié)果較差的原因，圖4給出了流域降雨的空間分布[圖 4(a)]、實測徑流深分布[圖 4(b)]以及磷礦位置。

圖4 降雨徑流空間分布Fig.4 The spatial distribution of annual rainfall, runoff depth in the study area

由圖 4(a)可以看出，磷礦均位于天福廟水庫以上的區(qū)域。流域的年降雨量從上游到下游呈遞減趨勢，河流右岸地區(qū)高于左岸地區(qū)，年降雨量中游最高，下游流域出口地區(qū)降雨量最小。比較圖 4(a)和圖 4(b)可以發(fā)現(xiàn)，天福廟以上區(qū)域的降雨量較大，但是并沒有產(chǎn)生相應(yīng)較大的徑流深。特別是在上游磷礦密集分布區(qū)域，降雨量大而徑流深卻小。由于玄廟觀、天福廟以上區(qū)域磷礦分布密集，我們認為造成徑流深減少的原因是由于磷礦開采活動引起的。磷礦開采過程中形成的地下裂隙、礦洞和礦坑，下滲的降水通過這些地下通道大量流失，從而減少了地表徑流。這與文獻[21]的研究結(jié)果相同，即磷礦開采減少了徑流量。

磷礦開采活動改變了流域下墊面條件，使地表徑流通過礦洞、裂隙等通道進入地下。由于SWAT模型對地下水過程的模擬比較簡單，僅僅分為淺層和深層含水層，因此使得玄廟觀、天福廟兩個水庫壩址處的徑流模擬精度遠低于另外兩個水庫壩址處的徑流模擬精度。如果需要提高天福廟以上區(qū)域的徑流模擬精度，需要改進SWAT模型中的地下水模擬模塊，或?qū)WAT模型與VISUAL MODFLOW等地下水模型耦合，使其能夠模擬磷礦開采對徑流的影響，提高流域徑流模擬的精度。

4 結(jié) 語

本文在黃柏河?xùn)|支流域構(gòu)建了SWAT模型，運用SWAT-CUP軟件中的SUFI-2算法進行了模型參數(shù)的敏感性分析、率定和不確定性分析，對研究區(qū)的月徑流進行了模擬。主要結(jié)論如下。

(1)參數(shù)敏感性分析結(jié)果表明，在所選取的19個與徑流相關(guān)的參數(shù)中，CN2、ESCO、SOL_K最為敏感，對研究區(qū)的月徑流模擬結(jié)果影響最為顯著。

(2)應(yīng)用SWAT模型進行黃柏河流域的月徑流模擬可以取得較滿意的結(jié)果，多數(shù)站點在率定期和驗證期的NSE均大于0.75，流量過程線擬合程度較好，表明SWAT模型在研究區(qū)具有較好的適用性。

(3)采用SUFI-2方法的不確定性分析表明，各個站點的不確定性都較小，大多數(shù)觀測值落在95PPU區(qū)間內(nèi)，驗證期的不確定性小于率定期。其中，天福廟站的徑流模擬精度雖然較高，但是其模擬結(jié)果的不確定性最大，說明模擬精度和不確定性結(jié)果并不一致，需要在模擬精度和模擬結(jié)果的不確定性之間取得平衡。

(4)模擬結(jié)果表明，流域上游徑流模擬精度低于下游，原因是上游區(qū)域大量存在磷礦開采活動，這一人類活動改變了流域的下墊面條件，增加了SWAT模型的模擬難度，預(yù)計通過改進SWAT模型的地下徑流模塊可以提高徑流模擬精度。

□