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

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

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

碧流河流域水資源是供應(yīng)城市用水的主要水源，徑流則是更直接的水源。研究評(píng)價(jià)SWAT模型對(duì)該流域的徑流模擬情況，結(jié)合流域水文要素變化以及利用自動(dòng)率定法率定參數(shù)并較準(zhǔn)確的預(yù)測(cè)徑流、模擬徑流意義重大，將為進(jìn)一步模擬未來(lái)氣候情景下的徑流奠定基礎(chǔ)。

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

故本文將SWAT模型應(yīng)用于碧流河流域徑流模擬。在分析流域水文要素變化的基礎(chǔ)上，結(jié)合自動(dòng)率定工具SWAT-CUP，利用其中的SUFI-2算法對(duì)參數(shù)自動(dòng)率定、敏感性排序，并對(duì)模擬結(jié)果分析評(píng)價(jià)，彌補(bǔ)了SWAT模型自帶校準(zhǔn)功能的不足，最后多角度分析模擬結(jié)果。

1 研究區(qū)概況

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

流域內(nèi)有小石硼、桂云花、天益等9個(gè)雨量站，有繭場(chǎng)、碧流河、玉石水庫(kù)3個(gè)水文站。其中，碧流河水文站作為徑流模擬站點(diǎn)(見(jiàn)圖1)。

圖1 碧流河流域示意

2 研究方法

2.1 SWAT分布式水文模型

SWAT模型是流域尺度的分布式動(dòng)態(tài)水文模型[2]，以水量平衡方程為問(wèn)題研究的驅(qū)動(dòng)力，綜合考慮氣候、下墊面等因素，模擬水文過(guò)程。該模型的分布式運(yùn)算特點(diǎn)如下：先根據(jù)DEM及實(shí)際水網(wǎng)把流域劃分為多個(gè)子流域；而后對(duì)子流域內(nèi)的水文響應(yīng)單元進(jìn)行聚類(lèi)分析，在此基礎(chǔ)上獲得不同時(shí)間尺度的模擬結(jié)果。

2.2 SUFI-2算法

SWAT-CUP中內(nèi)嵌了5種算法。李倩楠[8]等對(duì)5種不確定性評(píng)價(jià)方法進(jìn)行了比較，得出SUFI-2方法運(yùn)算效率高，運(yùn)行次數(shù)少，卻能得到較好的結(jié)果，故本研究選擇此法。

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

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

2.3 模型的評(píng)價(jià)標(biāo)準(zhǔn)

前述迭代結(jié)果顯示的5個(gè)指標(biāo)即為模型的評(píng)價(jià)指標(biāo)。其中模型的不確定性評(píng)價(jià)用：P-factor、R-factor；PBIAS、Ens和R23個(gè)指標(biāo)用于評(píng)價(jià)模擬精度，計(jì)算公式如下：

(1)

(2)

(3)

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

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

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

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

4 模型數(shù)據(jù)準(zhǔn)備

SWAT模型需要輸入的數(shù)據(jù)和來(lái)源如表1所示。氣象數(shù)據(jù)包括相對(duì)濕度、溫度、風(fēng)速、太陽(yáng)輻射等。

表1 模型數(shù)據(jù)

4.1 空間數(shù)據(jù)

本文使用的DEM數(shù)據(jù)為SRTM 90 m精度柵格型數(shù)據(jù)，2 000 s 1∶10萬(wàn)的土地利用分布圖，土壤數(shù)據(jù)選擇1∶100萬(wàn)數(shù)字化土壤圖件。將各空間數(shù)據(jù)投影類(lèi)型統(tǒng)一轉(zhuǎn)換為WGS_1984_UTM_ZONE_51N，并對(duì)應(yīng)于SWAT模型的代碼，重新分類(lèi)(見(jiàn)表2、3)。

表2 碧流河水庫(kù)上游土地利用類(lèi)型重分類(lèi)結(jié)果

表3 碧流河水庫(kù)上游土壤類(lèi)型

4.2 屬性數(shù)據(jù)

屬性數(shù)據(jù)的處理有：土壤屬性數(shù)據(jù)庫(kù)的建立和天氣發(fā)生器相關(guān)參數(shù)的計(jì)算。

由于土壤分類(lèi)系統(tǒng)記錄的土壤屬性的參數(shù)與SWAT模型輸入要求的土壤屬性標(biāo)準(zhǔn)不同，在建立土壤屬性數(shù)據(jù)庫(kù)前要將國(guó)際制標(biāo)準(zhǔn)轉(zhuǎn)換為美制標(biāo)準(zhǔn)。對(duì)土壤水特性參數(shù)的獲得方式為：借助軟件SPAW以及直接應(yīng)用。

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

5 結(jié)果分析

對(duì)碧流河水庫(kù)流域構(gòu)建SWAT模型，流域共劃分為29個(gè)子流域，碧流河水文站位于第21號(hào)子流域，限于當(dāng)前碧流河水文站缺乏日徑流數(shù)據(jù)，對(duì)該水文站進(jìn)行月模擬后輸出結(jié)果，并進(jìn)一步在SWAT-CUP中按算法要求建立SUFI-2工程、輸入相應(yīng)的要求文件得到模擬評(píng)價(jià)結(jié)果，進(jìn)而評(píng)價(jià)分析。

5.1 參數(shù)率定及敏感性分析

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

表4 參數(shù)敏感性及參數(shù)率定分析結(jié)果

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

按照相應(yīng)的修改方案設(shè)定參數(shù)初始范圍，經(jīng)過(guò)多次迭代(每次迭代次數(shù)設(shè)置為200次)逐步滿足評(píng)價(jià)指標(biāo)，得到參數(shù)敏感性分析以及率定結(jié)果如圖3及表4所示。

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

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

5.2 模型模擬結(jié)果的分析和評(píng)價(jià)

按照碧流河降雨和徑流變化趨勢(shì)分析結(jié)果，將1978年～2011年模擬期的前兩年作為預(yù)熱期，1981年～1995年作為率定期，1996年～2011年作為模型的驗(yàn)證期，建立模型并獲得模擬結(jié)果。利用前述3個(gè)評(píng)價(jià)指標(biāo)對(duì)模擬結(jié)果進(jìn)行評(píng)價(jià)。模型在不同時(shí)期的月值模擬結(jié)果如圖4、5所示，評(píng)價(jià)指標(biāo)見(jiàn)表5。

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

再?gòu)谋塘骱釉聫搅鞣植紒?lái)分析，如圖6給出了兩個(gè)時(shí)段(1981年～1995年、1996年～2011年)多年平均月徑流變化，6月～9月是碧流河豐水期，8月達(dá)到徑流峰值，之后開(kāi)始下降，10月至次年5月為枯水期，豐枯時(shí)段區(qū)分明顯，這與碧流河年內(nèi)降水時(shí)間分布有關(guān)。一年降水峰值主要集中在主汛期(7月～8月附近)，圖中同樣可見(jiàn)模型對(duì)峰值模擬效果較差。

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

圖5 月徑流驗(yàn)證期模擬結(jié)果(1996年～2011年)

時(shí)期P-factorR-factorR2EnsPBIAS/%率定期0.880.610.920.91-3.80驗(yàn)證期0.910.610.910.901.30

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

5.3 模型不確定性分析

水文過(guò)程包含著各種不確定性因素。對(duì)于水文模型，通常是運(yùn)用大量的概化或經(jīng)驗(yàn)公式對(duì)實(shí)際流域的水文循環(huán)過(guò)程的一個(gè)簡(jiǎn)化，各方面的不確定性往往相互作用而影響模型輸出結(jié)果的不確定性。

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

6 結(jié) 語(yǔ)

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