考慮大孔隙下滲的HBV模型及其在山洪預(yù)報(bào)中的應(yīng)用

崔逸凡，呂海深,2,3，朱永華,2，蘇建賓

(1.河海大學(xué) 水文水資源學(xué)院，江蘇 南京 210098；2.河海大學(xué) 水文水資源與水利工程科學(xué)國家重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210098；3 河海大學(xué) 全球變化與水循環(huán)國際合作聯(lián)合實(shí)驗(yàn)室，江蘇 南京 210098)

1 研究背景

近年來，由于全球氣候變化導(dǎo)致極端天氣頻發(fā)，突發(fā)性暴雨頻率快速增加，洪水災(zāi)害已成為造成人民生命財(cái)產(chǎn)損失的重大災(zāi)害，嚴(yán)重制約著廣大地區(qū)社會(huì)經(jīng)濟(jì)的發(fā)展[1]。據(jù)統(tǒng)計(jì)，中國是一個(gè)洪水災(zāi)害頻發(fā)的國家，洪水災(zāi)害造成的死亡人數(shù)占自然災(zāi)害死亡人數(shù)的比例呈上升趨勢(shì)，對(duì)人民生命安全造成了嚴(yán)重危害[2]，因此洪水的防治已迫在眉睫。洪水防治主要包括工程措施和非工程措施，其中，洪水預(yù)報(bào)作為洪水防治最重要的非工程措施已上升到國家戰(zhàn)略層面。目前，對(duì)于大江大河的洪水預(yù)報(bào)技術(shù)已經(jīng)比較成熟，而受限于對(duì)山洪形成機(jī)制認(rèn)識(shí)不全面、山區(qū)小流域水文資料匱乏和降水預(yù)報(bào)精度不高等因素，小流域山洪預(yù)報(bào)的精度往往并不能滿足實(shí)際需要。水文模型模擬是進(jìn)行洪水預(yù)報(bào)的重要方法，但是由于小流域山洪具有歷時(shí)短強(qiáng)度大、陡漲陡落等特點(diǎn)，許多模型對(duì)于小流域山洪的模擬精度并不高，容易出現(xiàn)洪峰模擬較差以及洪水過程線模擬不準(zhǔn)確等現(xiàn)象。

HBV模型是由瑞典國家水文氣象局(SMHI)在20世紀(jì)70年代開發(fā)的水文預(yù)報(bào)模型[3]，目前已被廣泛應(yīng)用于全球40多個(gè)國家的數(shù)百個(gè)流域的洪水預(yù)報(bào)和水資源評(píng)估中。近年來，HBV模型在國內(nèi)也已經(jīng)得到廣泛應(yīng)用，如趙彥增等[4]利用HBV模型對(duì)淮河官寨流域進(jìn)行了模擬研究，徑流過程模擬結(jié)果比較理想，可在國內(nèi)推廣應(yīng)用；王有恒等[5]利用HBV模型在白龍江流域進(jìn)行了模擬研究，得出HBV模型能較好地模擬洪水對(duì)降水的響應(yīng)過程。然而HBV模型對(duì)土壤層的處理較為簡單，僅將土壤層分為上、下兩層，其物理過程并不明確。在植被發(fā)育區(qū)，由于土壤脹縮、動(dòng)物活動(dòng)、植被發(fā)育及植被根系死亡等，大孔隙普遍存在。大孔隙雖然在土壤總孔隙中所占比例較小，但由大孔隙引發(fā)的優(yōu)先流對(duì)降雨的快速入滲和土壤水的快速補(bǔ)給產(chǎn)生了重要貢獻(xiàn)[6]，大孔隙的存在還會(huì)改變坡面尺度地表徑流、壤中流和地下徑流的形成過程和比例，提早地下徑流的出流時(shí)間和增加地下徑流在總徑流中的比重[7]，從而對(duì)徑流形成過程產(chǎn)生重要影響。

淮河流域位于我國南北氣候過渡帶，氣候變化復(fù)雜，降水時(shí)空分布不均勻，淮河上游干流及其南部山區(qū)的河道坡度大，匯流時(shí)間短，極易發(fā)生短歷時(shí)高強(qiáng)度的山洪[8]。同時(shí)，土壤大孔隙和大孔隙流在淮河流域普遍存在[9]，對(duì)徑流形成過程產(chǎn)生重要影響。因此，本文以淮河流域上游的黃泥莊小流域?yàn)檠芯繀^(qū)，引入變動(dòng)滲漏面積方法為HBV模型增添大孔隙模塊，探究改進(jìn)后的HBV模型在該流域山洪預(yù)報(bào)中的適用性。

2 研究區(qū)域及數(shù)據(jù)來源

2.1 研究區(qū)域

選擇淮河上游黃泥莊水文站以上集水區(qū)(下文簡稱黃泥莊流域)為研究區(qū)，該研究區(qū)介于北緯31°06′～31°38′和東經(jīng)115°21′～115°43′之間，流域控制面積805 km2。該流域地處我國南北氣候過渡帶，多年平均氣溫為11～16 ℃，流域多年平均年降水量為1 077 mm[10]，但降水的年內(nèi)分布不均。流域年降水量的50%～80%集中在汛期(6-9月)[10]，且多以暴雨形式出現(xiàn)。汛期降水歷時(shí)短、強(qiáng)度大、暴雨中心籠罩范圍小，因此汛期時(shí)由暴雨引發(fā)的山洪頻發(fā)。流域內(nèi)植被發(fā)育良好，森林覆蓋率高達(dá)65%，存在豐富的地下水和壤中流。流域平均海拔為479 m，但地形復(fù)雜，西北和東南以山地為主，河道坡度大，雙扇形河網(wǎng)使區(qū)域產(chǎn)流快速向流域出口匯集。流域概況如圖1所示。

圖1 黃泥莊流域概況

2.2 數(shù)據(jù)來源

如圖1所示，黃泥莊流域有相對(duì)密集的降水觀測(cè)站網(wǎng)，流域及其周邊的12個(gè)雨量站點(diǎn)的降水觀測(cè)數(shù)據(jù)可以從淮河流域水文年鑒資料中獲取。但是獲取的長歷時(shí)資料主要以日尺度觀測(cè)為主，僅在必要時(shí)進(jìn)行小時(shí)尺度加密觀測(cè)。因此，該項(xiàng)研究僅選擇了2010 -2015年間有小時(shí)尺度觀測(cè)數(shù)據(jù)的特定洪水事件進(jìn)行研究。考慮到12個(gè)雨量站點(diǎn)在黃泥莊流域分布相對(duì)均勻，在對(duì)雨量站的降水?dāng)?shù)據(jù)仔細(xì)檢查后，使用泰森多邊形面積權(quán)重方法獲取小時(shí)尺度流域面平均雨量。另外，由于雨量站缺乏氣溫觀測(cè)且流域范圍內(nèi)沒有國家氣象站，因此使用了高分辨率的歐洲中期天氣預(yù)報(bào)中心再分析產(chǎn)品(ECMWF)的逐小時(shí)氣溫?cái)?shù)據(jù)作為HBV模型的驅(qū)動(dòng)，該產(chǎn)品已經(jīng)被證明在中國區(qū)域有很好的應(yīng)用效果[11-12]。黃泥莊水文站的徑流觀測(cè)資料來自淮河流域水文年鑒洪水要素摘錄表。為了與模型輸出結(jié)果進(jìn)行匹配，研究中使用了線性插值方法對(duì)加密觀測(cè)的洪水摘錄資料進(jìn)行插值處理，進(jìn)而得到洪水事件期間逐小時(shí)的流量觀測(cè)資料。

考慮到下游的防汛承受能力，研究中以500 m3/s為閾值，在2010-2015年間共選擇了6次洪水事件(詳細(xì)信息如表1所示)，且在選定的洪水事件期間，12個(gè)雨量站都有完善的加密降水觀測(cè)資料。從表1中可以看出，2010-2015年期間每年都至少有1次較大洪水，且洪水事件主要集中在每年的6-7月間。其中，最嚴(yán)重的1次發(fā)生在2013年，洪峰流量高達(dá)1 628.49 m3/s。

表1 黃泥莊流域2010-2015年洪水事件概況

3 HBV模型改進(jìn)

3.1 概述

HBV模型為概念性水文模型，它具有輸入?yún)?shù)少、適用性強(qiáng)、模擬精度高等優(yōu)點(diǎn)，可以適用于各種復(fù)雜氣候條件下的水文模擬，目前已被廣泛應(yīng)用于水文預(yù)報(bào)、氣候變化研究、水資源評(píng)估、地下水模擬和水量平衡研究等方面[13-14]。HBV模型以降水量、氣溫和潛在蒸散發(fā)為主要驅(qū)動(dòng)，輸出出口斷面徑流量。如圖2所示，HBV模型主要由3個(gè)部分組成，分別是：積雪與融雪模塊、土壤濕度計(jì)算模塊和徑流響應(yīng)模塊。考慮到大孔隙對(duì)HBV徑流響應(yīng)過程的影響，研究中使用變動(dòng)滲漏面積方法對(duì)上層水庫進(jìn)行改進(jìn)，詳細(xì)描述如下。

圖2 考慮大孔隙的HBV模型結(jié)構(gòu)

3.2 積雪和融雪模塊

HBV模型采用閾值方法進(jìn)行雨雪分割(此次研究中使用0 ℃為雨雪分割閾值)，且模型假定融雪和積雪速率與溫度成正比。當(dāng)溫度低于閾值溫度時(shí)，降水主要以雪的形式累積為積雪，反之則為降雨。當(dāng)溫度一直保持在閾值溫度以下時(shí)，則不會(huì)產(chǎn)生徑流，但溫度一旦超過閾值溫度，融雪徑流開始產(chǎn)生，其中融雪速率如下式：

Sm=DD(T-Tt)

(1)

式中：Sm為融雪速率，mm/h；DD為度日因子，mm/(℃·h)；T為小時(shí)平均氣溫，℃；Tt為閾值氣溫，℃。

3.3 土壤濕度計(jì)算模塊

流域的降水通常分為兩部分：一部分滲入土壤；另一部分用于形成地表徑流，即為時(shí)段徑流深ΔQ。一般來說，土壤含水量越高，降水對(duì)徑流形成的貢獻(xiàn)就越大，當(dāng)土壤含水量接近最大土壤含水量時(shí)，下滲量減少，降水對(duì)產(chǎn)流量的貢獻(xiàn)增加。時(shí)段徑流深的計(jì)算如下：

(2)

式中：ΔQ為時(shí)段徑流深，mm；SM為實(shí)際土壤含水量，mm；FC為最大土壤含水量，mm；β為形狀系數(shù)。

由于蒸散發(fā)與氣溫息息相關(guān)，模型采用氣溫折算方法獲取小時(shí)尺度潛在蒸散發(fā)。如公式(3)所示，利用小時(shí)平均氣溫與日平均氣溫之間的差值來修訂日平均潛在蒸散發(fā)，進(jìn)而獲取小時(shí)尺度的潛在蒸散發(fā)。考慮到實(shí)際蒸散發(fā)受土壤含水量影響很大且與修正后的潛在蒸散發(fā)呈線性關(guān)系，因此當(dāng)土壤含水量大于或等于土壤凋萎含水量時(shí)，實(shí)際蒸散發(fā)就等于修正后的潛在蒸散發(fā)；當(dāng)土壤含水量小于土壤凋萎含水量時(shí)，受限于土壤水分不足，實(shí)際蒸散發(fā)將使用土壤含水量與凋萎含水量的比例進(jìn)行折算。

PEa=(1+C(T-Tm))PEm

(3)

(4)

式中：PEa為修正后的潛在蒸散發(fā)，mm；Tm為日平均氣溫，℃；PEm為日平均潛在蒸散發(fā)，mm；C為模型參數(shù)；Ea為實(shí)際蒸散發(fā)，mm；PWP為土壤凋萎含水量，mm。

3.4 改進(jìn)后的徑流響應(yīng)模塊

HBV模型將徑流形成過程概化為一個(gè)響應(yīng)函數(shù)，該函數(shù)由一個(gè)上層非線性水庫和一個(gè)下層線性水庫組成，即把徑流成分概化為兩個(gè)響應(yīng)盒子。上層盒子有兩個(gè)徑流出口，其響應(yīng)函數(shù)通過消退系數(shù)K0、K1來體現(xiàn)，下層盒子有一個(gè)徑流出口，其響應(yīng)函數(shù)通過消退系數(shù)K2來體現(xiàn)。當(dāng)上層盒子水位超過閾值L時(shí)，上層盒子迅速形成地表徑流(Q0)，并從第一個(gè)出口流出。上層盒子的第二個(gè)出口和下層盒子出口的徑流響應(yīng)相對(duì)較慢，形成壤中流(Q1)和基流(Q2)。流域出口徑流深(Q)即為各部分徑流深之和。

(5)

(6)

Q=Q0+Q1+Q2

(7)

式中：Su為上層水庫水位，mm；Sl為下層水庫水位，mm；L為閾值水位，mm；K0為地表徑流消退系數(shù)；K1為壤中流消退系數(shù)；K2為基流消退系數(shù)；Q0為地表徑流，mm；Q1為壤中流，mm；Q2為基流，mm；Q為總徑流深，mm。

但是，HBV模型中并未考慮大孔隙對(duì)降雨入滲的促進(jìn)作用以及大孔隙的存在對(duì)徑流組成成分的影響，所以模擬出的洪水過程線可能不準(zhǔn)確，尤其是對(duì)洪峰的模擬精度較差?？紤]到這一點(diǎn)，研究中將研究區(qū)劃分為設(shè)置模擬大孔隙下滲的變動(dòng)滲漏面積和一般透水面積兩部分，且變動(dòng)滲漏面積與土壤含水量有關(guān)。參考肖慶元等[9]的研究結(jié)果，變動(dòng)滲漏面積和一般透水面積劃分如下：

(8)

式中：AA為變動(dòng)滲漏面積所占比例；BB為一般透水面積所占比例；IA為AA的上限；n為經(jīng)驗(yàn)指數(shù)。

在增加大孔隙模塊后，改進(jìn)的HBV模型中，將上層盒子分為兩個(gè)部分：變動(dòng)滲漏部分和一般透水部分。在變動(dòng)滲漏部分，下滲能力很強(qiáng)，一般情況下水可以沿著大孔隙直接進(jìn)入下層盒子而無損失；在一般透水部分，仍然依據(jù)HBV模型的滲透常數(shù)方法進(jìn)行，改進(jìn)后的上層土壤滲流可以表述為：

Qprec=Su·Kprec+AA·ΔQ

(9)

式中：Kperc為滲流消退系數(shù)；Qprec為滲流，mm。

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

根據(jù)《水文情報(bào)預(yù)報(bào)規(guī)范》(SL 250-2000)，采用Nash-Sutcliffe效率系數(shù)即Nash系數(shù)、洪峰流量相對(duì)誤差、洪量相對(duì)誤差和峰現(xiàn)時(shí)間誤差來衡量模型的模擬效果，它們的計(jì)算公式和最優(yōu)值如表2所示。Nash系數(shù)反映的是徑流模擬值與觀測(cè)值之間的吻合度，洪峰流量相對(duì)誤差反映的是對(duì)洪峰的模擬效果，洪量相對(duì)誤差反應(yīng)的是對(duì)洪水總量的模擬效果，峰現(xiàn)時(shí)間誤差反映的是對(duì)洪峰出現(xiàn)時(shí)間的模擬效果。

表2 模型評(píng)價(jià)指標(biāo)

4 模擬結(jié)果和分析

4.1 參數(shù)敏感性分析

在改進(jìn)的HBV模型中，IA表示變動(dòng)滲漏區(qū)在流域中所占的最大比例；n為經(jīng)驗(yàn)指數(shù)，反映隨著土壤濕度變化，變動(dòng)滲漏區(qū)中實(shí)際大孔隙區(qū)域的變化。由公式(8)可知，當(dāng)IA為定值時(shí)，隨著n的增加，流域中大孔隙實(shí)際所占的面積逐漸減少；反之，當(dāng)n為定值時(shí)，隨著IA的增加，實(shí)際大孔隙面積比例增加。當(dāng)IA=0時(shí)，變動(dòng)滲漏面積所占比例為0，可認(rèn)為土壤中不存在大孔隙；當(dāng)IA=1，n=0時(shí)，變動(dòng)滲漏面積所占比例為1，可認(rèn)為土壤中全是大孔隙。為分析模型改進(jìn)效果，并確定新增大孔隙模塊主要控制參數(shù)的敏感性，不失一般性的選擇20100705次洪水為例，采用控制變量法分析參數(shù)IA和n對(duì)洪峰流量和洪水總量的影響，即使IA在0.1～1之間以0.1的步長均勻變化，n在0～0.8之間以0.2的步長均勻變化，統(tǒng)計(jì)IA和n的變化對(duì)洪峰流量和洪水總量的影響，結(jié)果如圖3和4所示。

由圖3和4可知，洪峰流量和洪水總量對(duì)新增加大孔隙模塊的兩個(gè)參數(shù)都很敏感。隨著IA的增加和n的減小，流域中大孔隙所占面積逐漸增加，洪峰流量和洪水總量均隨之減小且洪峰流量的變率遠(yuǎn)超洪水總量。這主要是由于流域孔隙度的增加促進(jìn)了土壤的快速入滲進(jìn)而產(chǎn)生了更多的壤中流，而壤中流的匯流時(shí)間較長使得模擬洪水過程線坦化，降低了洪峰的模擬值。深入分析可以發(fā)現(xiàn)，IA和n之間有著復(fù)雜的相互關(guān)聯(lián)。當(dāng)IA較小時(shí)，變動(dòng)滲漏區(qū)域在流域中的占比較小，洪峰流量和洪水總量對(duì)n的變化響應(yīng)強(qiáng)度較弱，隨著IA的增加，曲線的分布更加分散，表明此時(shí)主要的控制參數(shù)為經(jīng)驗(yàn)參數(shù)n。總的來說，對(duì)于IA而言，其高值區(qū)間比低值區(qū)間更加敏感，而對(duì)n而言則是低值區(qū)間更加敏感。

4.2 模擬結(jié)果分析

基于改進(jìn)后的HBV模型，對(duì)2010-2015年間選定的6次洪水事件進(jìn)行模擬，其中前4次洪水事件用來校準(zhǔn)模型，采用試錯(cuò)法獲取模型的最優(yōu)參數(shù)[15]，后2次洪水事件用來檢驗(yàn)?zāi)Ｐ徒Y(jié)果。

表3匯總了6次洪水事件模擬結(jié)果的統(tǒng)計(jì)指標(biāo)，總的來說，改進(jìn)后的HBV模型在黃泥莊流域洪水模擬中模擬效果較好。率定期的4次洪水事件Nash系數(shù)分別為0.82、0.92、0.91和0.89；驗(yàn)證期的兩次洪水事件Nash系數(shù)分別為0.92和0.76，能較好地匹配兩個(gè)時(shí)期洪水過程。除明顯低估了20130705次洪水的洪峰流量外(洪峰流量相對(duì)誤差為-28.38 %)，其余洪峰流量和洪量相對(duì)誤差均在20 %以下，能較好地捕獲洪峰流量和洪量，對(duì)災(zāi)害預(yù)警和災(zāi)情評(píng)估有重要意義。此次率定期和驗(yàn)證期的峰現(xiàn)時(shí)間誤差均在2 h以內(nèi)，改進(jìn)后的HBV模型能夠準(zhǔn)確地再現(xiàn)洪峰出現(xiàn)時(shí)間。依據(jù)《水文情報(bào)預(yù)報(bào)規(guī)范》(SL 250-2000)，6次洪水事件的模擬精度均達(dá)乙級(jí)及以上標(biāo)準(zhǔn)，滿足實(shí)際工程應(yīng)用精度，可用于該流域洪水模擬與預(yù)報(bào)。

圖3IA-n-洪峰流量相關(guān)圖 圖4IA-n-洪水總量相關(guān)圖

表3 黃泥莊流域洪水模擬結(jié)果

圖5和6分別展示了率定期和驗(yàn)證期洪水過程線的觀測(cè)值和改進(jìn)后HBV模型的模擬結(jié)果。如圖1所示，黃泥莊流域由南北兩個(gè)扇形子流域構(gòu)成，洪水漲落迅速，主要以單峰的尖細(xì)洪水過程為主(圖5和6)。盡管改進(jìn)后的HBV模型能較好地模擬洪水過程線，但依然傾向于低估洪峰。這可能是受限于集總式模型結(jié)構(gòu)特征，改進(jìn)后的HBV模型依然難以刻畫徑流過程對(duì)流域降水空間異質(zhì)性的響應(yīng)。同時(shí)，黃泥莊流域的雨量站點(diǎn)均分布在靠近河道的相對(duì)平緩的谷地，對(duì)坡面地形雨及短歷時(shí)暴雨的捕獲率不足，難以準(zhǔn)確反映洪峰對(duì)應(yīng)期流域降水，這也進(jìn)一步加劇了HBV模型對(duì)洪峰流量的低估。另外，可能是受到前期土壤濕度的影響，黃泥莊流域水文過程線對(duì)小雨并不敏感，這也導(dǎo)致HBV模型高估了洪水過程的起漲點(diǎn)(例如20100705和20130705洪水)，后續(xù)的研究需要針對(duì)這一過程進(jìn)行進(jìn)一步的改進(jìn)。集總式HBV因其具有模型結(jié)構(gòu)相對(duì)簡單、計(jì)算量小、模型參數(shù)和所需驅(qū)動(dòng)數(shù)據(jù)較少等優(yōu)點(diǎn)，已經(jīng)被廣泛應(yīng)用于洪水模擬與預(yù)報(bào)，但對(duì)降雨徑流過程概化不足也進(jìn)一步增加了模型的不確定性。總體而言，添加對(duì)物理過程的精細(xì)化描述能進(jìn)一步提高模型的模擬精度。

圖5 模型率定期觀測(cè)和模擬的洪水過程線

圖6 模型驗(yàn)證期觀測(cè)和模擬的洪水過程線

5 結(jié) 論

無資料小流域的山洪預(yù)報(bào)一直是洪水預(yù)報(bào)中的一個(gè)難點(diǎn)，土壤大孔隙及其引發(fā)的大孔隙流通過改變流域產(chǎn)匯流過程，進(jìn)一步增加了這類地區(qū)洪水模擬與預(yù)報(bào)的難度。本文使用變動(dòng)滲漏面積法改進(jìn)HBV模型，并以淮河流域上游黃泥莊小流域?yàn)槔?，?duì)2010- 2015年的6次洪水事件進(jìn)行模擬，得出了以下結(jié)論：

(1)大孔隙模塊的控制參數(shù)IA和n的變化會(huì)對(duì)洪峰流量和洪水總量產(chǎn)生較大影響。

(2)率定期的Nash系數(shù)在0.82～0.92之間，驗(yàn)證期的Nash系數(shù)在0.76～0.92之間，洪峰流量相對(duì)誤差和洪量相對(duì)誤差基本控制在20 %以內(nèi)，峰現(xiàn)時(shí)間誤差不超過2 h。3次洪水事件的模擬精度達(dá)到洪水預(yù)報(bào)甲級(jí)標(biāo)準(zhǔn)，3次達(dá)到乙級(jí)標(biāo)準(zhǔn)，說明改進(jìn)后的HBV模型能夠很好地模擬黃泥莊流域的洪水過程，滿足該流域的實(shí)際應(yīng)用需求。

(3)本研究可為淮河流域山區(qū)的防洪減災(zāi)、水資源管理調(diào)度等提供參考。同時(shí)，引入變動(dòng)面積方法對(duì)HBV模型的改進(jìn)和對(duì)存在大孔隙的相似小流域山洪預(yù)報(bào)研究有重要的借鑒意義。