基于SWAT模型和SUFI-2算法的碧流河流域徑流模擬

劉 寧，張 霞，祝雪萍，趙志懷，武鵬林

(太原理工大學水利科學與工程學院，山西太原030024)

碧流河流域水資源是供應城市用水的主要水源，徑流則是更直接的水源。研究評價SWAT模型對該流域的徑流模擬情況，結合流域水文要素變化以及利用自動率定法率定參數并較準確的預測徑流、模擬徑流意義重大，將為進一步模擬未來氣候情景下的徑流奠定基礎。

基于物理基礎的SWAT模型在研究流域水循環的影響方面作出了積極貢獻，可以模擬流域內發生的各種物理過程，被應運于很多方面[1-2]。在非點源污染方面，Li等研究SWAT模型在江西省漣水流域進行多尺度的泥沙模擬研究以及氮和磷的輸移過程[3]。在徑流模擬方面，姚蘇紅等[4]以內蒙古閃電河流域為例，應用SWAT模型進行徑流模擬，結果顯示模擬結果受參數影響很大。張余慶[5]等學者以及宋增芳[6]等學者，均利用SWAT結合SUFI-2分別對修水流域、石羊河流域進行徑流模擬，有效控制了參數對模擬結果的影響。曹明亮[7]以碧流河流域為例，利用SWAT模型率定法對多站點進行研究，驗證了模型的適用性；但在參數敏感性分析、模擬結果評價以及不確定分析方面使用SWAT模型自帶的調參功能進行參數校準，該方法不僅耗時長、分析率定的參數較少、對經驗要求高且易受人為因素影響。

故本文將SWAT模型應用于碧流河流域徑流模擬。在分析流域水文要素變化的基礎上，結合自動率定工具SWAT-CUP，利用其中的SUFI-2算法對參數自動率定、敏感性排序，并對模擬結果分析評價，彌補了SWAT模型自帶校準功能的不足，最后多角度分析模擬結果。

1 研究區概況

碧流河屬于北溫帶濕潤氣候區，6月～9月為汛期，尤其集中在7月～8月多年平均降雨量為742.8 mm，多年平均徑流深為294.5 mm，多年多年平均蒸發量為1 404.4 mm，流域面積2 814 km2，壩址以上控制流域面積2 085 km2。1975年興建的碧流河水庫是流域內的主要水庫(總庫容9.34億m3)，主要功能為向大連市供水。流域內以林地和耕地為主，以草甸土、棕壤性土為主要土壤類型。2001年建成了控制流域面積為313 km2的玉石水庫。

流域內有小石硼、桂云花、天益等9個雨量站，有繭場、碧流河、玉石水庫3個水文站。其中，碧流河水文站作為徑流模擬站點(見圖1)。

圖1 碧流河流域示意

2 研究方法

2.1 SWAT分布式水文模型

SWAT模型是流域尺度的分布式動態水文模型[2]，以水量平衡方程為問題研究的驅動力，綜合考慮氣候、下墊面等因素，模擬水文過程。該模型的分布式運算特點如下：先根據DEM及實際水網把流域劃分為多個子流域；而后對子流域內的水文響應單元進行聚類分析，在此基礎上獲得不同時間尺度的模擬結果。

2.2 SUFI-2算法

SWAT-CUP中內嵌了5種算法。李倩楠[8]等對5種不確定性評價方法進行了比較，得出SUFI-2方法運算效率高，運行次數少，卻能得到較好的結果，故本研究選擇此法。

SUFI-2是利用拉丁超立方體抽樣法獲得參數值帶入模型進行模擬[9]。其敏感性分析通過Global以及One-at-a-time敏感性分析實現，其原理可參見SWAT-CUP說明手冊。此外，該算法通過P-factor(95%不確定性區間內的觀測數據的百分數)和R-factor(95PPU上下限的平均距離與標準偏差的比值)評價參數不確定性程度，P-factor的理論值區間為0～100%、R-factor的理論值區間為0～+∞。當P-factor=1，且R-factor=0時，為完全理想的模擬結果；故通常以P-factor→1，且R-factor→0的程度來判斷模擬效果。但P-factor的增大會導致R-factor也增大，因此要找到二者的相對最佳值。

SUFI-2算法：首先，需定義較大的參數范圍(結合文獻分析以及SUFI-2推薦的參數范圍確定)；然后，按照SUFI-2算法運行要求，設定各參數范圍、模擬次數及實測值等文件，并選擇目標函數；最后，由程序進行多次迭代，依據迭代結果納西效率系數(Ens)、相關系數(R2)、相對偏差(PBIAS)、P-factor和R-factor的值)調整參數范圍，直至獲得理想的迭代結果。

2.3 模型的評價標準

前述迭代結果顯示的5個指標即為模型的評價指標。其中模型的不確定性評價用：P-factor、R-factor；PBIAS、Ens和R23個指標用于評價模擬精度，計算公式如下：

(1)

(2)

(3)

當PBIAS<25%，Ens>0.5，R2>0.6時，模型的模擬結果可信；Ens≥0.65時，認為模擬結果很好[10-11]。P-factor>0.5且R-factor<1.5時，模擬不確定性可接受[12]；當P-factor>0.5且P-factor>0.7時，模擬的不確定性較小[13]。

3 碧流河水文要素變化分析

水文要素主要包括降雨、徑流、溫度、等水文循環相關的時間序列。其中，降雨是徑流的直接來源，而徑流則是降雨最直觀的輸出形式。為了通過歷史趨勢了解水文要素的總體趨勢，參照文獻[14]給出的1958年～2011年的分析結果，對研究區間1978年～2011年的降水及徑流變化趨勢進行了分析(見圖2)。其分析可得出相似的結論：降雨和徑流都為下降趨勢，徑流的下降趨勢更明顯。此外，與多年平均徑流量對比可知，豐水年在1979年、1985年及1995年出現、枯水年在2002年、2011年出現；同時，1995年之前，偏豐水年多于偏枯水年，而1995年之后偏枯水年多于偏豐水年，同樣驗證了徑流的下降趨勢，其下降的原因可能與當下氣候全球變暖導致的降水、氣溫等氣象要素有關，有待下一步研究。

圖2 碧流河水庫流域1958年～2011年降雨徑流系列變化

4 模型數據準備

SWAT模型需要輸入的數據和來源如表1所示。氣象數據包括相對濕度、溫度、風速、太陽輻射等。

表1 模型數據

4.1 空間數據

本文使用的DEM數據為SRTM 90 m精度柵格型數據，2 000 s 1∶10萬的土地利用分布圖，土壤數據選擇1∶100萬數字化土壤圖件。將各空間數據投影類型統一轉換為WGS_1984_UTM_ZONE_51N，并對應于SWAT模型的代碼，重新分類(見表2、3)。

表2 碧流河水庫上游土地利用類型重分類結果

表3 碧流河水庫上游土壤類型

4.2 屬性數據

屬性數據的處理有：土壤屬性數據庫的建立和天氣發生器相關參數的計算。

由于土壤分類系統記錄的土壤屬性的參數與SWAT模型輸入要求的土壤屬性標準不同，在建立土壤屬性數據庫前要將國際制標準轉換為美制標準。對土壤水特性參數的獲得方式為：借助軟件SPAW以及直接應用。

SWAT模型天氣發生器的參數計算由于工作量較大，運用計算程序SwatWeather.exe(北京師范大學建)進行，將氣象數據整理后導入程序，得到模型需要的數據[15]。

5 結果分析

對碧流河水庫流域構建SWAT模型，流域共劃分為29個子流域，碧流河水文站位于第21號子流域，限于當前碧流河水文站缺乏日徑流數據，對該水文站進行月模擬后輸出結果，并進一步在SWAT-CUP中按算法要求建立SUFI-2工程、輸入相應的要求文件得到模擬評價結果，進而評價分析。

5.1 參數率定及敏感性分析

SWAT模型能夠相對準確地模擬徑流的關鍵在于將眾多的SWAT模型中與徑流有關的參數選出，去除對模擬結果影響較小的參數，進行敏感性分析、率定。本文通過前述SWAT-CUP的敏感性分析方法最終選出14個最為敏感的參數校準，部分參數含義具體如下：CN2為平均濕度時的SCS曲線數，值越大，徑流總量越大；SOL_AWC為田間持水量越大，可下滲水量越少，滲透時間越短；GW_DELAY用于計算每天流入地下水層的回灌水量，與水平面的深度和滲漏區及地下水區地質組成的水力特性有關；GW_REVAP值越大，再蒸發量越多，當等于1時，蒸發量等于潛蒸發量；CH_K2：默認為0，值越大，河道輸移損失越大；REVAPMN：只有當淺層地下水層的含水量超過閾值aqshthr，rvp時，再蒸發過程才發生。CH_N2：值越大，匯流流速越小；RCHRG_DP值越大，回灌水流入深層地下水層的比例越大；SOL_K為值越大，由土壤層流入下一土壤層的水量就越多，同時產生的壤中流也越多。

表4 參數敏感性及參數率定分析結果

注：①R、V為兩種修改參數的方法，分別是乘以某值(百分百浮動)、賦值；②t-stat越大表示參數越敏感，p-value值表示t值的顯著程度，越接近0值越顯著。

按照相應的修改方案設定參數初始范圍，經過多次迭代(每次迭代次數設置為200次)逐步滿足評價指標，得到參數敏感性分析以及率定結果如圖3及表4所示。

圖3 參數敏感性分析(P-Value及t-stat)結果

圖3及表4呈現率定后參數的最終范圍、最佳值及敏感性判定原則，排在第一位敏感的參數為：CN2，其P值最大，t檢驗值最小，同樣的判斷原理，接下來依次為：SOL_AWS、GW_REVAP、GW_REVAP、GWQMN等，其對模擬結果影響逐漸減小。

5.2 模型模擬結果的分析和評價

按照碧流河降雨和徑流變化趨勢分析結果，將1978年～2011年模擬期的前兩年作為預熱期，1981年～1995年作為率定期，1996年～2011年作為模型的驗證期，建立模型并獲得模擬結果。利用前述3個評價指標對模擬結果進行評價。模型在不同時期的月值模擬結果如圖4、5所示，評價指標見表5。

模擬結果顯示：①碧流河水文站的模擬值與觀測值變化趨勢一致性良好、流量過程擬合良好，率定期、驗證期R2分別為0.92、0.91；Ens分別為0.91、0.90；PBIAS分別均控制在5%以內，且全部落在較小的不確定區間內(95PPU)。②但對峰值的模擬效果較差，模擬值普遍小于實測值，可能由于未考慮融雪影響而致。③驗證期(1995年～2011年)的模擬效果略低于率定期(1980年～1995年)的模擬效果，可能與SWAT模型對豐、枯水年的模擬精度不同有關，偏枯水年模擬精度低于偏豐水年模擬精度[16]。

再從碧流河月徑流分布來分析，如圖6給出了兩個時段(1981年～1995年、1996年～2011年)多年平均月徑流變化，6月～9月是碧流河豐水期，8月達到徑流峰值，之后開始下降，10月至次年5月為枯水期，豐枯時段區分明顯，這與碧流河年內降水時間分布有關。一年降水峰值主要集中在主汛期(7月～8月附近)，圖中同樣可見模型對峰值模擬效果較差。

圖4 月徑流率定期模擬結果(1981年～1995年)

圖5 月徑流驗證期模擬結果(1996年～2011年)

時期P-factorR-factorR2EnsPBIAS/%率定期0.880.610.920.91-3.80驗證期0.910.610.910.901.30

圖6 觀測和模擬的多年平均月徑流變化

5.3 模型不確定性分析

水文過程包含著各種不確定性因素。對于水文模型，通常是運用大量的概化或經驗公式對實際流域的水文循環過程的一個簡化，各方面的不確定性往往相互作用而影響模型輸出結果的不確定性。

本文不確定性結果顯示，率定期P-factor為0.88、R-factor為0.61，依據判別條件，滿足P-factor>0.7、R-factor<1不確定性較小的條件，模型模擬結果較適用。驗證期與率定期基本一致，而略高的P-factor表明：90%的觀測值落在95%置信水平的不確定性區間內，模擬的不確定性相對更小。但從圖4、5可看出，豐水期模擬的不確定性區間依然較大，而峰值模擬的不確定性更大。

6 結 語

SWAT模型在碧流河流域的徑流模擬中在前人手動率定較少參數的基礎上，延長模擬時期，利用自動率定SUFI-2算法提升了模擬效果，并驗證了相對較小的不確定性，為下一步的氣候模式與SWAT模型耦合筑牢了基礎；同時也發現，該模型對碧流河的峰值模擬效果不佳、不確定性只是綜合考慮，應在后續的研究中深入綜合考慮更多氣候變化與人類活動等因素。另外限于水文資料的限制，僅利用一個碧流河水文站序列進行分析，代表性方面有一定欠缺，需要進一步獲取更多水文站數據深入研究；此外，敏感性分析得到的14個參數其各參數對結果的影響貢獻率多大，有待繼續探討。