基于SWAT模型和SUFI-2算法的黃柏河東支流域徑流模擬研究

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

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

0 引 言

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

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

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

1 數據與方法

1.1 研究區概況

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

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

圖1 研究區概況Fig.1 Overview of the study area

1.2 數據來源及處理

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

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

圖2 研究區土地利用圖及土壤類型圖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城鎮URBN17.771.91ALh高活性強酸土16.474.46 合計 932.26100合計932.26100

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

1.3 研究方法

1.3.1 SWAT模型

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

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

1.3.2 SWAT-CUP和SUFI-2算法

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

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

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

2 結果與分析

2.1 敏感性分析結果

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

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

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

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

2.2 參數率定和驗證結果

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

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

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

表3 黃柏河東支流域月徑流模擬結果統計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%，說明模型在率定期和驗證期月徑流模擬結果較觀測值偏低。整體來看，3個評價指標均在誤差范圍內，SWAT模型能夠較好地模擬研究區的月徑流過程，在流域的月徑流模擬中具有較好的適應性。

2.3 不確定性分析結果

不確定性結果見表3。在驗證期和率定期各站點的P-factor均大于0.65，R-factor均小于1、表明大多數觀測值落在95PPU區間內，所構建的SWAT模型徑流模擬的不確定性較小。

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

3 討 論

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

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

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

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

4 結 語

本文在黃柏河東支流域構建了SWAT模型，運用SWAT-CUP軟件中的SUFI-2算法進行了模型參數的敏感性分析、率定和不確定性分析，對研究區的月徑流進行了模擬。主要結論如下。

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

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

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

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

□