熊丁暉, 桂發(fā)二, 張洪賓, 周 泉, 羅 源, 李韓文, 陳 華
(浙江貴仁信息科技股份有限公司, 浙江 杭州 310051)
成屏一級(jí)水庫位于浙江省麗水市遂昌縣甌江支流松蔭溪上游,距遂昌縣城12 km,是松蔭溪梯級(jí)開發(fā)的首級(jí)電站。壩址以上集水面積為185 km2,主流長(zhǎng)為28 km。庫區(qū)多年平均降雨量1 670 mm,多年平均流量5.84 m3/s。庫區(qū)雨量時(shí)空分布不均,在多雨季節(jié),尤其在7-8月份極易造成洪澇災(zāi)害,而在相對(duì)干旱的1-2月卻常處于缺水狀態(tài)。由于庫區(qū)受局地氣候影響較大,水文氣象要素不穩(wěn)定,小范圍內(nèi)地形起伏大,在降雨集中情況下易出現(xiàn)突發(fā)狀況。因此,研究成屏一級(jí)水庫入庫流量的快速、準(zhǔn)確預(yù)報(bào)有著重要的工程應(yīng)用價(jià)值。
水文模型是水庫入庫流量模擬的重要工具之一,國(guó)內(nèi)外已廣泛應(yīng)用水文模型進(jìn)行水庫入庫流量研究。李藝婷等[1]將新安江模型的蓄滿產(chǎn)流模式應(yīng)用于浙江省臺(tái)州市山區(qū)型和平原型兩個(gè)水庫的入庫徑流量模擬之中,取得較好的效果。魏恒志等[2]采用流溪河模型開展了白龜山水庫入庫洪水預(yù)報(bào)研究,結(jié)果表明流溪河模型具有良好的水文模擬性能。Anghileri等[3]將VIC模型應(yīng)用于加利福尼亞奧羅維爾水庫,模擬效果較好。然而,國(guó)內(nèi)關(guān)于HYPE模型[4]在水庫入庫流量模擬方面的研究較少。本文以成屏一級(jí)水庫為研究對(duì)象,應(yīng)用HYPE模型推算長(zhǎng)序列日尺度水庫入庫流量信息,探討HYPE模型在成屏一級(jí)水庫的適用性。
此外,國(guó)內(nèi)外現(xiàn)有研究表明,雖然應(yīng)用水文模型模擬流域日尺度流量過程總體效果良好,但是在洪峰模擬上略有不足,主要體現(xiàn)在模擬值較大幅度地低于或高于洪峰。金君良等[5]建立VIC模型模擬嘉陵江流域日徑流過程,結(jié)果表明日徑流過程的Nash效率系數(shù)在0.7以上,但幾場(chǎng)暴雨過程對(duì)應(yīng)的洪峰較高于模擬值。楊柳等[6]應(yīng)用SWAT模型模擬山美水庫集水區(qū)流域日尺度的徑流過程,其日徑流過程的Nash效率系數(shù)在0.8以上,但是依然有若干洪峰高于實(shí)測(cè)值。Golmohammadi等[7]同時(shí)評(píng)估了MIKE-SHE模型、APEX模型和SWAT模型在格蘭德流域日尺度徑流模擬的效果,結(jié)果表明3個(gè)模型總體效果較好,但是洪峰流量大于或小于模擬值。以上研究表明,水文模型在長(zhǎng)序列日尺度徑流模擬存在洪峰模擬不準(zhǔn)確的問題,本文提出一種方法擬嘗試解決這個(gè)問題。
成屏一級(jí)水庫壩址以上流域地形特征和雨量站分布見圖1,土地利用和土壤類型見圖2。雨量站有桂洋站、石柱站、垵口站、大公坑站和成屏一級(jí)站,流域出口位于成屏一級(jí)站。土地利用類型包括農(nóng)田、森林、草地、水體、城鎮(zhèn)等,其中森林面積占94%。土壤類型包括黏土、壤土和砂壤土,其中壤土由于組成比例不同分為壤土1和壤土2,砂壤土也由于組成比例不同分為砂壤土1和砂壤土2。
使用HYPE模型模擬成屏一級(jí)水庫入庫流量需要的數(shù)據(jù)包括數(shù)字高程圖、土地利用圖、土壤類型圖、氣象數(shù)據(jù)、實(shí)測(cè)入庫流量數(shù)據(jù)等。
(1)數(shù)字高程圖。從SRTM數(shù)據(jù)庫(http://srtm.csi.cgiar.org/)下載分辨率為90 m×90 m的數(shù)字高程圖。
(2)土地利用圖。從USGS數(shù)據(jù)庫(https://www.usgs.gov/)下載分辨率為1 km×1 km 的2010年土地利用圖。
(3)土壤類型圖。從聯(lián)合國(guó)糧農(nóng)組織(FAO)和維也納國(guó)際應(yīng)用系統(tǒng)研究所(IIASA)所構(gòu)建的世界土壤數(shù)據(jù)庫HWSD (http://webarchive.iiasa.ac.at/Research/LUC/External-World-soil-database/HTML/index.html?sb=1)下載分辨率為1km×1km的土壤類型圖。
(4)氣象數(shù)據(jù)。從遂昌縣氣象局收集2013-2015年的降雨數(shù)據(jù)和溫度數(shù)據(jù)。
(5)實(shí)測(cè)入庫數(shù)據(jù)。從遂昌縣氣象局收集2013-2015年的成屏一級(jí)水庫入庫流量數(shù)據(jù)。

圖1 成屏一級(jí)水庫壩址以上流域地形特征和雨量站分布
2.3.1 HYPE模型 HYPE(Hydrological Predictions for the Environment)模型是瑞典國(guó)家氣象水文研究所在HBV模型的基礎(chǔ)上研發(fā)的半分式水文模型,該模型綜合考慮了多種影響流域徑流的因素,將水文運(yùn)動(dòng)過程、氮磷輸移轉(zhuǎn)化過程和模型率定等結(jié)合起來,形成了一套完整的模型預(yù)報(bào)方案。該模型具有易于操作、結(jié)構(gòu)清晰、參數(shù)率定速度快、模擬效果好等優(yōu)點(diǎn),被廣泛應(yīng)用于水量平衡研究、洪水設(shè)計(jì)、氣候變化影響研究以及地下水模擬和水文預(yù)報(bào)等諸多方面[8-10]。
圖3為HYPE模型的結(jié)構(gòu)示意圖。研究區(qū)域被劃分為若干個(gè)子流域,每個(gè)子流域又被劃分為不同的類,主要包括土地類、湖泊類和河流類,土地類由不同的土壤類型和土地利用類型組合而成。圖4為HYPE模型的水文過程原理圖。HYPE模型模擬的主要水文過程包括融雪和積雪、蒸散發(fā)、入滲、大孔隙流、滲流、地表徑流、壤中流、排水管流和地下徑流。除了模擬水文過程之外,HYPE模型還可以模擬徑流及河流水體中營(yíng)養(yǎng)物的濃度,如有機(jī)氮、無機(jī)氮、總氮、溶解態(tài)磷、顆粒態(tài)磷及總磷等。詳細(xì)物理過程參考文獻(xiàn)[4]。
2.3.2 LH-OAT方法 采用拉丁超立方-單次單因素(Latin hypercube-One factor At a Time,LH-OAT)方法進(jìn)行HYPE模型參數(shù)敏感性分析[11]。LH-OAT方法將參數(shù)敏感度表示為一個(gè)無量綱的指數(shù),反映模型輸出結(jié)果隨模型參數(shù)的微小改變而變化的敏感性程度[12]。
具體而言,假設(shè)模型有P個(gè)需要分析的參數(shù),首先將這些參數(shù)的范圍劃分為N層,然后分別在每一層進(jìn)行1次抽樣,得到1個(gè)拉丁超立方抽樣點(diǎn)(包括P個(gè)參數(shù)的集合),總共進(jìn)行N次抽樣。之后對(duì)每個(gè)拉丁超立方抽樣點(diǎn)中的參數(shù)進(jìn)行一定程度的擾動(dòng),每次只改變1個(gè)參數(shù),共有P次擾動(dòng)。因此,N個(gè)拉丁超立方抽樣點(diǎn)最終可產(chǎn)生N×(P+1)個(gè)參數(shù)組,將各參數(shù)的敏感度平均后即可得到其全局敏感度。該方法結(jié)合了拉丁超立方抽樣算法的穩(wěn)定性和單次單因素算法的精確性,大幅度降低了模型運(yùn)行的次數(shù)和時(shí)間,提高了模型的計(jì)算效率。LH-OAT方法步驟如下:(1)在設(shè)定的P個(gè)參數(shù)范圍內(nèi),采用拉丁超立方抽樣并隨機(jī)生成N個(gè)拉丁超立方抽樣點(diǎn);(2)采用單次單因素方法生成N×(P+1)個(gè)參數(shù)組;(3)針對(duì)生成的參數(shù)組,重復(fù)運(yùn)行HYPE模型N×(P+1)次,輸出變量為Nash-Sutcliffe系數(shù);(4)計(jì)算各參數(shù)的全局敏感度,得到各參數(shù)最終的全局敏感度及排序結(jié)果。
全局敏感度計(jì)算方法如公式(1)。

(1)
式中:M為模型輸出變量;N為抽樣次數(shù);γ為縮放因子;ε為擾動(dòng)幅度;GSi為第i個(gè)參數(shù)的全局敏感度值。
2.3.3 DE算法 模型采用微分進(jìn)化蒙特卡洛算法——DE算法[13](Differential Evolution Markov Chain method)進(jìn)行參數(shù)率定。DE算法的基本思想是:對(duì)種群中的每一個(gè)個(gè)體,從當(dāng)前種群中隨機(jī)選擇3個(gè)點(diǎn),以其中1個(gè)點(diǎn)為基礎(chǔ),另外兩個(gè)點(diǎn)作為參照作一次擾動(dòng),所得點(diǎn)與該個(gè)體交叉后進(jìn)行“自然選擇”,保留較優(yōu)者,實(shí)現(xiàn)種群的進(jìn)化。DE算法適用于無約束連續(xù)變量的全局優(yōu)化問題,包括非線性規(guī)劃、線性規(guī)劃、非光滑優(yōu)化等。
DE算法正是貝葉斯推理與微分進(jìn)化兩種算法的結(jié)合,它采用微分進(jìn)化算法和馬爾科夫鏈蒙特卡羅模擬對(duì)后驗(yàn)分布進(jìn)行采樣,獲得參數(shù)的后驗(yàn)邊緣概率密度函數(shù),在此基礎(chǔ)上獲得參數(shù)的數(shù)學(xué)期望等統(tǒng)計(jì)量。
2.3.4 評(píng)價(jià)指標(biāo) 本文采用4個(gè)指標(biāo)[14]評(píng)價(jià)HYPE模型模擬效果,分別是決定系數(shù)(R2)、百分比偏差(percent bias,PBIAS)、納什效率系數(shù)(Nash-Sutcliffe efficiency,NSE)、和均方根誤差與標(biāo)準(zhǔn)差比值(ratio of root mean square error to standard deviation,RSR),其公式如下:
(2)
(3)
(4)
(5)
采用ArcGIS軟件對(duì)成屏一級(jí)水庫以上流域進(jìn)行劃分,網(wǎng)格大小為1 km×1 km,得到80個(gè)子流域和每個(gè)子流域的面積、平均海拔、河長(zhǎng)以及子流域之間的流向,采用泰森多邊形方法得到流域平均降雨量,80個(gè)子流域和泰森多邊形見圖5。根據(jù)土地利用圖和土壤類型圖得到14個(gè)類,這14個(gè)類的土地利用和土壤類型組成見表1,其中土壤共分成3層,第1層深30 cm,第2層深50 cm,第3層深20 cm。

表1 成屏一級(jí)水庫壩址以上流域14個(gè)類信息
注:土地利用中的1為農(nóng)田;2為森林;3為草地;4為水體;5為城鎮(zhèn);土壤類型中的1為黏土;2為砂壤土1;3為壤土1;4為砂壤土2;5為壤土2。

圖2 成屏一級(jí)水庫壩址以上流域土地利用和土壤類型
HYPE模型參數(shù)可分為4類:獨(dú)立參數(shù)、與土地利用類型相關(guān)的參數(shù)、與土壤類型相關(guān)的參數(shù)以及和區(qū)域相關(guān)的參數(shù)。由于HYPE模型涉及到的參數(shù)眾多,首先根據(jù)經(jīng)驗(yàn)進(jìn)行手動(dòng)篩選,然后采用LH-OAT方法進(jìn)行參數(shù)敏感性分析,最終確定的參數(shù)有17個(gè),見表2。

圖3 HYPE模型結(jié)構(gòu)示意圖

圖4 HYPE模型水文過程原理圖

圖5 成屏一級(jí)水庫壩址以上流域80個(gè)
采用LH-OAT方法進(jìn)行參數(shù)敏感性分析時(shí),將HYPE模型的17個(gè)參數(shù)分成10層進(jìn)行LH抽樣,隨機(jī)得到10個(gè)LH抽樣點(diǎn),將各個(gè)參數(shù)值除以相應(yīng)參數(shù)的范圍得到相對(duì)大小,結(jié)果見圖6。
對(duì)抽樣參數(shù)組進(jìn)行OAT分析計(jì)算,計(jì)算次數(shù)為10×(17+1)=180次,使用目標(biāo)函數(shù)NSE計(jì)算每個(gè)參數(shù)的全局敏感度,縮放因子γ取100;擾動(dòng)幅度ε取0.01,結(jié)果見圖7。從圖7中可以看出,極敏感的參數(shù)為cevp,較敏感的參數(shù)為pcluse、cevpcorr、rrcscorr、rrcs1,一般敏感的參數(shù)為lp、preccorr、wcep、pcaddg,較不敏感的參數(shù)為cevpam、tcelevadd、rivvel、wcfc、epotdist、rrcs2、damp和mperc1。

表2 HYPE模型參數(shù)及范圍

圖6 LH抽樣結(jié)果

圖7 HYPE模型各參數(shù)全局敏感度
在確定HYPE模型各參數(shù)敏感度排序的情況下,采用DE算法進(jìn)行參數(shù)率定。率定期選擇2013年5月14日至2014年12月31日,參數(shù)率定結(jié)果見表3。

表3 HYPE模型各參數(shù)率定結(jié)果
注:S表示土壤類型,S1表示黏土;S2表示砂壤土1;S3表示壤土1;S4表示砂壤土2;S5表示壤土2。L表示土地利用,L1表示農(nóng)田;L2表示森林;L3表示草地;L4表示水體;L5表示城鎮(zhèn)。
采用HYPE模型對(duì)成屏一級(jí)水庫入庫流量進(jìn)行模擬,率定期選擇2013年5月14日至2014年12月31日,驗(yàn)證期選擇2015年1月1日至2015年12月31日。率定期和驗(yàn)證期的HYPE模型模擬值和實(shí)測(cè)值對(duì)比分別見圖8和9。
從圖8中可以看出,率定期HYPE模型模擬值和實(shí)測(cè)值吻合程度較高,除了2014年8月20日的一場(chǎng)洪峰模擬值和實(shí)測(cè)值存在較大的差異,模擬值遠(yuǎn)低于實(shí)測(cè)值。從圖9中可以看出,驗(yàn)證期HYPE模型模擬值和實(shí)測(cè)值基本吻合,同樣2015年8月9日的一場(chǎng)洪峰,模擬值遠(yuǎn)低于實(shí)測(cè)值。此外,其他若干場(chǎng)洪峰模擬值也低于實(shí)測(cè)值。
為了進(jìn)一步分析HYPE模型的模擬效果,采用4個(gè)評(píng)價(jià)指標(biāo)進(jìn)行評(píng)估,見表4。從表4中可以看出,率定期和驗(yàn)證期的NSE均大于0.8,R2均大于0.9,RSR均小于0.5,PBIAS的絕對(duì)值均在25%之間,表明不管是在率定期還是在驗(yàn)證期,HYPE模型對(duì)成屏一級(jí)水庫日尺度入庫流量模擬均達(dá)到要求,展現(xiàn)了HYPE模型較強(qiáng)的模擬能力。

表4 率定期和驗(yàn)證期HYPE模型模擬效果的4個(gè)評(píng)價(jià)指標(biāo)
為了進(jìn)一步分析HYPE模擬值與實(shí)測(cè)值之間的誤差,將率定期和驗(yàn)證期的模擬值與實(shí)測(cè)值之差進(jìn)行比較,見圖10。結(jié)果表明在實(shí)測(cè)值較小時(shí),模擬值和實(shí)測(cè)值之差在0附近,表明模擬值很接近實(shí)測(cè)值;在實(shí)測(cè)值較大時(shí),模擬值與實(shí)測(cè)值之差偏離0較大,表明模擬值偏離實(shí)測(cè)值。上述結(jié)果表明HYPE模型能夠較好地模擬實(shí)測(cè)值較小的情況,而在實(shí)測(cè)值較大的情況下模擬效果較差。圖11為率定期和驗(yàn)證期模擬值與實(shí)測(cè)值之差隨時(shí)間的波動(dòng)圖,同樣反映了HYPE模型在洪峰模擬中的不足。究其原因,本文認(rèn)為有兩點(diǎn):(1)輸入數(shù)據(jù)的精確性有待商榷,尤其是降雨量數(shù)據(jù)。降雨量輸入在很大程度上影響模擬結(jié)果,由于采用泰森多邊形法對(duì)各個(gè)降雨站的雨量求平均,能否完全代替實(shí)際流域的雨量分布是存在疑問的[15]。此外,降雨量和實(shí)測(cè)徑流量數(shù)據(jù)也可能存在一定程度的人為誤差。(2)HYPE模型參數(shù)可能并不能反映出暴雨過程中的流域特征。由于在暴雨過程中,流域下墊面等情況發(fā)生了較大的改變,并且水流存在湍流等現(xiàn)象,導(dǎo)致暴雨過程中的流域特征和較小降雨情況下的流域特征存在明顯差異。
HYPE模型參數(shù)是針對(duì)長(zhǎng)序列進(jìn)行率定的,其率定結(jié)果更反映較小降雨情況下的流域特征,而不能刻畫暴雨情況下的流域特征,所以導(dǎo)致在洪峰過程中模擬值比實(shí)測(cè)值較小的現(xiàn)象。也許針對(duì)不同的降雨情況采用不同的參數(shù)進(jìn)行刻畫會(huì)更加符合現(xiàn)實(shí)情況。
為了刻畫暴雨情況下的流域特征,選擇2014年8月1日-31日的一場(chǎng)洪峰進(jìn)行率定,選擇2015年8月1日-31日的一場(chǎng)洪峰進(jìn)行驗(yàn)證,率定的參數(shù)見表5。其他時(shí)段的模擬值保持一致,考慮暴雨情況的率定期和驗(yàn)證期的模擬值和實(shí)測(cè)值對(duì)比圖分別見12和13。結(jié)果表明2014年8月份和2015年8月份的洪峰模擬值更符合實(shí)測(cè)值。
考慮暴雨特征的率定期和驗(yàn)證期HYPE模型模擬效果的4個(gè)評(píng)價(jià)指標(biāo)見表6。對(duì)比表4和6,率定期和驗(yàn)證期的NSE和R2變化不大,而表6的RSR和PBIAS明顯小于表4,表明考慮暴雨特征之后的模擬值和實(shí)測(cè)值的偏差減小。

圖8 率定期HYPE模型模擬值和實(shí)測(cè)值對(duì)比圖

圖9 驗(yàn)證期HYPE模型模擬值和實(shí)測(cè)值對(duì)比圖

圖10 率定期和驗(yàn)證期模擬值和實(shí)測(cè)值之差和實(shí)測(cè)值的比較(縱坐標(biāo)表示模擬值和實(shí)測(cè)值之差)

圖11 率定期和驗(yàn)證期模擬值與實(shí)測(cè)值之差隨時(shí)間的波動(dòng)(縱坐標(biāo)表示模擬值和實(shí)測(cè)值之差)

圖12 考慮暴雨情況的率定期模擬值與實(shí)測(cè)值對(duì)比圖

圖13 考慮暴雨情況的驗(yàn)證期模擬值與實(shí)測(cè)值對(duì)比圖

表5 考慮暴雨情況的HYPE模型各參數(shù)率定結(jié)果
注:S表示土壤類型,S1表示黏土;S2表示砂壤土1;S3表示壤土1;S4表示砂壤土2;S5表示壤土2。L表示土地利用,L1表示農(nóng)田;L2表示森林;L3表示草地;L4表示水體;L5表示城鎮(zhèn)。

表6 考慮暴雨特征的率定期和驗(yàn)證期HYPE模型模擬效果的4個(gè)評(píng)價(jià)指標(biāo)
結(jié)合LH-OAT方法、DE方法和HYPE模型,通過參數(shù)敏感性分析、參數(shù)率定和誤差分析,對(duì)成屏一級(jí)水庫入庫流量模擬效果進(jìn)行了較為全面的研究。得到如下結(jié)論:
(1)HYPE模型中極敏感的參數(shù)為cevp,較敏感的參數(shù)為pcluse、cevpcorr、rrcscorr、rrcs1,一般敏感的參數(shù)為lp、preccorr、wcep、pcaddg,較不敏感的參數(shù)為cevpam、tcelevadd、rivvel、wcfc、epotdist、rrcs2、damp和mperc1。
(2)HYPE模型率在率定期的納什效率系數(shù)(NSE)為0.87,決定系數(shù)(R2)為0.94,和均方根誤差與標(biāo)準(zhǔn)差比值(RSR)為0.42,百分比偏差(PBIAS)為-23.30%,在驗(yàn)證期的納什效率系數(shù)(NSE)為0.83,決定系數(shù)(R2)為0.92,和均方根誤差與標(biāo)準(zhǔn)差比值(RSR)為0.48,百分比偏差(PBIAS)為-20.88%。表明HYPE模型具備較強(qiáng)的模擬能力,能夠用于長(zhǎng)序列日尺度水文模擬。
(3)HYPE模型在洪峰模擬上模擬值低于實(shí)測(cè)值,原因可能是輸入數(shù)據(jù)的精確性不足以及率定的模型參數(shù)并不能較好地反映出暴雨過程中的流域特征。
(4)對(duì)暴雨情況采用新的參數(shù)進(jìn)行流域特征刻畫,模型模擬效果更好。