劉楊合,程 磊,張?jiān)品?,張艷軍,衛(wèi)曉婧
(1.武漢大學(xué) 水資源與水電工程科學(xué)國家重點(diǎn)實(shí)驗(yàn)室, 武漢 430072; 2.武漢大學(xué) 水利水電學(xué)院, 武漢 430072)
山洪是突發(fā)性降雨導(dǎo)致的山區(qū)徑流陡漲陡落的自然現(xiàn)象,受降雨強(qiáng)度和下墊面等因素影響,可能伴有泥石流和滑坡等災(zāi)害的發(fā)生,嚴(yán)重威脅著國民經(jīng)濟(jì)增長和人民生命財(cái)產(chǎn)安全[1]。小流域山區(qū)洪水具有突發(fā)性、局地性和頻次多等特點(diǎn)[2],不同流域山洪致災(zāi)機(jī)理差異性較大,受地形條件和資料缺乏所限,以及人為因素影響,小流域山區(qū)洪水一直存在著模擬精度低、預(yù)報(bào)精度差的問題[3-4]。因此,開展小流域山區(qū)洪水模擬的相關(guān)研究對保障社會經(jīng)濟(jì)增長、維護(hù)人民生命財(cái)產(chǎn)安全等具有至關(guān)重要的作用。
秦巴地區(qū)是我國重要的山洪災(zāi)害防治區(qū),1950—2000年防治區(qū)內(nèi)共有小流域2 411個,山洪溝1 765條,山洪災(zāi)害次數(shù)共8 251次,其山洪災(zāi)害嚴(yán)重程度僅次于江南丘陵和華南丘陵地區(qū)[3-4]。本研究選取的官山小流域,以中小丘陵為主,多暴雨和長歷時降雨,具有典型的秦巴地區(qū)氣候和地形特征。官山小流域內(nèi)山洪頻發(fā),自1975年以來,發(fā)生了12場危害性較大的山洪災(zāi)害。其中1975年山洪沖毀河堤24.2 km,受淹農(nóng)田面積2 426 hm2;2012年8月5日山洪沖毀公路160 km,毀壞橋梁29座,轉(zhuǎn)移人口約7 000人,直接經(jīng)濟(jì)損失達(dá)2.3億元。
作為秦巴地區(qū)的典型流域,官山小流域洪水模擬相關(guān)研究已初步開展,韓培等[5]利用GIS技術(shù)獲取了官山小流域的地形、坡度、植被覆蓋度、土地利用類型等下墊面特征,為洪水模擬提供了基礎(chǔ)數(shù)據(jù)支撐。王以逵等[6]研究了官山河兩河口河段的暴雨洪水頂托效應(yīng)。黃艷等[7]比較分析了改進(jìn)的SCS(Soil Conservation Service)、TOPMODEL和新安江模型在官山河流域洪水模擬中的應(yīng)用效果。然而當(dāng)前官山小流域洪水模擬仍然存在以下2個問題:①傳統(tǒng)水文模型洪水模擬精度較低;②各模型模擬不同場次洪水的能力差異大,水文模型挑選難。
本研究利用官山小流域2009—2015年10場降雨徑流關(guān)系較好的場次洪水,分別采用新安江、TUWMODEL和TOPMODEL等模型進(jìn)行模擬?;诟髂P湍M結(jié)果,采用分流量區(qū)間BMA(Bayesian Model Averaging)方法[8]綜合計(jì)算官山小流域場次洪水過程,并與3種單一模型模擬精度進(jìn)行比較分析,以期提高官山小流域山洪模擬精度,為提高秦巴地區(qū)洪水預(yù)報(bào)精度提供技術(shù)支撐。
官山小流域坐落于漢江中上游、丹江口市西南部,位于110°48′00″E—111°34′59″E,32°13′16″N—32°58′20″N,出口孤山水文站控制的流域面積為332 km2,流域人口密度達(dá)241人/km2。以山地丘陵為主,海拔240~1 606 m,平均高程690 m,河道平均坡降5.7%,多年平均流量7.78 m3/s。流域降雨以局地性暴雨和連陰雨為主,年降水量為1 100 mm。山洪及其造成的泥石流和滑坡災(zāi)害頻發(fā)。研究區(qū)降雨由孤山站、袁家河站、大馬站和西河站的降雨觀測數(shù)據(jù)通過泰森多邊形法計(jì)算得到,流量數(shù)據(jù)選用孤山站的水文監(jiān)測數(shù)據(jù)。研究區(qū)的地理位置、雨量站和水文站的分布如圖1所示。
圖1 官山小流域地理位置Fig.1 Location of Guanshan small watershed
TOPMODEL模型是一個基于地形的半分布式水文模型,由Beven和Kirkby于1979年提出[9]。模型自提出以來,已在水文領(lǐng)域得到了廣泛應(yīng)用,其核心在于利用地形指數(shù)ln(α/tanβ)來反映以徑流運(yùn)動分布規(guī)律為主的流域水文現(xiàn)象[10]。其中,α為流域面積,β為流域坡度。
TOPMODEL模型主要采用飽和產(chǎn)流機(jī)制,即地下水達(dá)到地表時,土壤含水量達(dá)到或超過飽和含水量,流域產(chǎn)生飽和坡面流。模型具有3層結(jié)構(gòu),分別為植被根系層(Srz)、土壤非飽和層(Suz)和飽和地下水層(D);將目標(biāo)流域劃分為多個單元網(wǎng)格,降水首先按一定速率下滲進(jìn)入植被根系層補(bǔ)給該部分的土壤缺水量,植被根系層的水分損失為蒸發(fā)和下滲,當(dāng)降水滿足該層的土層儲水能力和蒸發(fā)損失后,剩余部分水分下滲進(jìn)入土壤非飽和層。土壤非飽和層中的水分按一定重力排水速率VG下滲進(jìn)入飽和地下水層,并最終經(jīng)由側(cè)向徑流形成壤中流Qb。飽和地下水層受到土壤非飽和層的補(bǔ)給后,隨著地下水位的上升,流域下游河道會產(chǎn)生飽和坡面流Qs。流域總徑流Q等于壤中流和坡面流的總和,即Q=Qb+Qs。
新安江模型是集總式水文模型,由趙人俊于1984年提出[11]。該模型結(jié)構(gòu)簡單,模型參數(shù)概念明確,且已在我國南方濕潤流域的水文模擬研究得到廣泛應(yīng)用。模型核心為蓄水容量曲線,首先計(jì)算各蓄滿單元流域的產(chǎn)匯流,得到各單元流域的出口流量過程后,進(jìn)行出口以下的河道洪水演算,疊加各單元流域的出流過程,進(jìn)而獲得流域總徑流過程。模型的蒸散發(fā)結(jié)構(gòu)為3層;產(chǎn)流計(jì)算過程采用嚴(yán)格的蓄滿產(chǎn)流模式,以田間持水量為控制條件,流域單元土壤含水量達(dá)到田間持水量,流域單元內(nèi)凈雨全部轉(zhuǎn)化為徑流;總徑流為地表徑流、壤中流和地下徑流三者之和;地表徑流匯流計(jì)算采用單位線法;壤中流和地下徑流的匯流計(jì)算采用線性水庫法。
TUWMODEL模型是維也納工業(yè)大學(xué)開發(fā)的集總式水文模型。模型的核心是產(chǎn)流隨土壤濕度呈指數(shù)變化。該模型包含融雪模塊、土壤水模塊和徑流模塊。融雪模塊通過簡單的度日因子法計(jì)算融雪過程;土壤水模塊表示產(chǎn)流和流域土壤水的變化。產(chǎn)流機(jī)制與新安江模型不同,模型采用接近蓄滿產(chǎn)流的概念,產(chǎn)流量主要受土壤濕度控制,是土壤濕度的指數(shù)函數(shù),與降雨強(qiáng)度無關(guān),土壤越濕,越容易產(chǎn)流;徑流模塊包含上下層水庫,土壤水模塊的水量可以下滲到上層水庫,上層水庫的水分會產(chǎn)生地表徑流和壤中流并可以通過下滲進(jìn)入下層水庫,下層水庫以一個較小的速率產(chǎn)生徑流。上下層水庫產(chǎn)生的徑流三角函數(shù)進(jìn)行匯流計(jì)算。
3.4.1 BMA方法
BMA方法是通過加權(quán)平均多個模型的模擬結(jié)果以期獲得更為可靠的綜合模擬結(jié)果[12]。基于BMA方法原理,可以根據(jù)k個模型(f1,f2,f3,…,fk)的獨(dú)立分布,計(jì)算綜合各模型模擬流量Q的總概率分布P(Q|OBD),進(jìn)而計(jì)算綜合各模型模擬結(jié)果的模擬流量Q。其計(jì)算過程如下。
總概率分布為
式中:OBD為觀測流量樣本數(shù)據(jù);fi為所選模型;pi(Q|fi,OBD)為所選模型的后驗(yàn)分布。
綜合模擬流量Q通過計(jì)算總概率分布的均值獲得,即
式中:g(Q|fi,σ2)為模型fi同觀測流量做線性回歸獲得的分布;σ為所得分布的標(biāo)準(zhǔn)差。令p(fi|OBD)=ωi,則可得
(3)
3.4.2 流量區(qū)間劃分
各模型結(jié)構(gòu)和產(chǎn)匯流機(jī)制不同,導(dǎo)致了不同模型在模擬不同量級的洪水流量過程時,往往具有不同的表現(xiàn)能力和相對優(yōu)勢。相較于傳統(tǒng)的BMA方法,分流量區(qū)間BMA方法考慮不同模型模擬不同量級洪水過程時相對優(yōu)勢的差異,為模型分配一個隨洪水流量量級發(fā)生變化的模型權(quán)重,可以更好地發(fā)揮不同模型模擬不同量級洪水過程的優(yōu)勢,提高洪水模擬精度[8]。本研究首先將率定期內(nèi)所有場次的雨洪過程拼接成長序列的雨洪過程,以該長序列的降雨流量數(shù)據(jù)分別率定新安江、TUWMODEL和TOPMODEL水文模型。按分位數(shù)(10%、25%、50%、75%、90%、100%)將長序列流量過程劃分為6個子區(qū)間,分別提取出3個單一模型在對應(yīng)子區(qū)間模擬的流量過程。在6個子區(qū)間內(nèi),以劃分出的實(shí)際流量過程為BMA方法中的觀測樣本數(shù)據(jù)(OBD,見3.4.1節(jié)),以模型模擬的流量過程為模型后驗(yàn)分布(pi(Q|fi,OBD),見3.4.1節(jié)),基于BMA方法分別計(jì)算各子區(qū)間內(nèi)新安江、TUWMODEL和TOPMODEL各水文模型的模型權(quán)重。檢驗(yàn)期洪水計(jì)算則根據(jù)洪水流量所處區(qū)間,采用分區(qū)間BMA方法計(jì)算的模型權(quán)重和檢驗(yàn)期模型模擬的流量過程,加權(quán)平均各模型優(yōu)勢的模擬流量過程。
模型誤差評價(jià)指標(biāo)如下:
(1)納什效率系數(shù)NSE計(jì)算公式為
(4)
(2)水量平衡系數(shù)VE計(jì)算公式為
(5)
(3)洪峰流量相對誤差Qe計(jì)算公式為
(6)
式中:Qe為洪峰流量相對誤差;Qp,sim和Qp,obs分別為洪峰流量模擬值和觀測值。本研究中流域洪峰流量平均相對誤差為所有場次洪峰流量相對誤差的均值,許可誤差為實(shí)測值的±20%[15]。
(4)洪量相對誤差Ve計(jì)算公式為
(7)
式中:Ve為洪量相對誤差;Vsim和Vobs分別為洪量的模擬值和觀測值。本研究中流域洪量平均相對誤差為所有場次洪量相對誤差的均值,許可誤差為實(shí)測值的±20%[15]。
(5)峰現(xiàn)時間誤差Te計(jì)算公式為
Te=Tsim-Tobs。
(8)
式中:Te為峰現(xiàn)時間誤差;Tsim和Tobs分別為峰現(xiàn)時間的模擬值和觀測值。本研究中流域平均峰現(xiàn)時間誤差為所有場次峰現(xiàn)時間誤差絕對值的均值,許可誤差為3 h[15]。
(6)合格率QR計(jì)算公式為
(9)
式中:QR為合格率;N為合格的模擬次數(shù);M為總的模擬次數(shù)。
針對以上各指標(biāo),納什效率系數(shù)越接近于1,說明該模型模擬洪水流量的整體過程越好;VE越接近于0,說明該模型模擬洪水過程中水量平衡模擬得越好;洪峰流量、洪量相對誤差和峰現(xiàn)時間相對誤差越接近于0,說明該模型模擬的洪峰流量、洪量和峰現(xiàn)時間越準(zhǔn)確;洪峰流量、洪量和峰現(xiàn)時間合格率越接近于100%,說明使用該模型模擬的洪峰流量、洪量和峰現(xiàn)時間的可靠性越高。
圖2對比分析了分流量區(qū)間BMA方法和3個水文模型的流量模擬結(jié)果在率定期和檢驗(yàn)期的水量平衡系數(shù)和納什效率系數(shù)之間的差異。結(jié)果表明:分流量區(qū)間BMA方法在檢驗(yàn)期的納什效率系數(shù)顯著高于3個單一水文模型,率定期和檢驗(yàn)期的中位數(shù)分別為0.87和0.79;水量平衡系數(shù)中位數(shù)分別為-8.8×10-11和0.02,比3個單一水文模型的水量平衡系數(shù)的中位數(shù)更接近0。
圖2 率定期及檢驗(yàn)期各模型評價(jià)指標(biāo)對比Fig.2 Comparison of evaluation indices in calibration period and verification period among different models
圖3 20100825號洪水期間降雨徑流過程 Fig.3 Rainfall and runoff processes during flood event 20100825
圖4 檢驗(yàn)期洪峰、洪量及峰現(xiàn)時間平均相對誤差和合格率Fig.4 Averaged relative errors and pass rate of flood peak, flood volume, and time of peak in verification period
分流量區(qū)間BMA方法可以計(jì)算出模型在各流量區(qū)間的模型權(quán)重,分析模型權(quán)重在不同區(qū)間的相對大小,可以為官山小流域挑選適用水文模型的提供參考。圖5顯示了分區(qū)間BMA方法估計(jì)的新安江模型、TUWMODEL和TOPMODEL在不同流量區(qū)間的模型權(quán)重。由圖5可以看出,新安江模型僅在流量區(qū)間為90%~100%時具有較高的模型權(quán)重。TUWMODEL在各區(qū)間的模型權(quán)重的中位數(shù)(除90%~100%外)均大于新安江模型和TOPMODEL,且在50%~75%和75%~90%區(qū)間內(nèi)高達(dá)0.85和0.82。TOPMODEL在各區(qū)間的模型權(quán)重中位數(shù)均較小,其模型權(quán)重中位數(shù)僅在25%~50%區(qū)間內(nèi)略高于新安江模型。其原因可能在于不同區(qū)間洪水的產(chǎn)匯流機(jī)制不同,各模型僅在符合自身產(chǎn)匯流假設(shè)的條件下得到更好的模擬效果,從而被分配到較高的權(quán)重。新安江模型采用蓄滿產(chǎn)流機(jī)制[16-17],故僅在降雨量較高,流域發(fā)生大洪水時,官山小流域的土壤含水量趨近于飽和時具有更高的模型權(quán)重。TUWMODEL采用的產(chǎn)流機(jī)制不是嚴(yán)格的蓄滿產(chǎn)流,其產(chǎn)流量是土壤含水量的指數(shù)函數(shù),受降雨強(qiáng)度約束較小,主要受土壤濕度控制,土壤越濕流域越容易產(chǎn)流[16, 18-19]。該產(chǎn)流機(jī)制可能比較符合官山小流域的暴雨山洪產(chǎn)流特性。雖然TOPMODEL考慮了地形因素,但因其模型的表現(xiàn)受洪峰影響較大,而所選的10場洪水,洪水量級差異明顯,故在模型率定過程中難以確定最優(yōu)參數(shù),致使模型表現(xiàn)性能較差[7, 20]。
圖5 新安江模型、TUWMODEL和TOPMODEL權(quán)重Fig.5 Weights of Xin’anjiang model, TUWMODEL, and TOPMODEL
針對小流域山區(qū)洪水模擬精度較低及缺乏水文資料帶來的洪水模擬可靠性低的問題,本研究以秦巴地區(qū)內(nèi)的典型小流域官山小流域?yàn)槔?,選取了2009—2015年10場降雨徑流關(guān)系較好的場次洪水,將其中任意7場作為率定期,剩余3場作為檢驗(yàn)期,分別率定了新安江模型、TUWMODEL和TOPMODEL。基于以上3種模型的模擬結(jié)果,利用分流量區(qū)間BMA方法估計(jì)了綜合各模型優(yōu)勢的模擬洪水流量過程,并對比了該方法的模擬結(jié)果和其他3種單一模型的差異。結(jié)果顯示,分流量區(qū)間BMA方法相比于其他3種單一水文模型更好地模擬了該區(qū)域的洪水過程。其檢驗(yàn)期納什效率系數(shù)的中位數(shù)為0.8,遠(yuǎn)高于其他3種單一水文模型。其檢驗(yàn)期的洪峰流量、洪量和峰現(xiàn)時間合格率分別為30.6%、70.8%和82.2%,高于其他3種單一水文模型,較于最優(yōu)水文模型分別高出4.5%、39.7%和48.9%。模型的權(quán)重結(jié)果顯示,TUWMODEL相較于其他2種水文模型在官山小流域更為適用,其原因在于TUWMODEL的產(chǎn)流機(jī)制更符合官山小流域。
本研究僅選取了3個水文模型運(yùn)行分流量區(qū)間BMA方法,未來的研究中可以增加水文模型的數(shù)量,選取更多的具有不同產(chǎn)匯流機(jī)制的水文模型運(yùn)行分區(qū)間BMA方法;也可使用不同目標(biāo)函數(shù)率定的水文模型運(yùn)行分流量區(qū)間BMA方法,來全面分析不同目標(biāo)函數(shù)對不同流量區(qū)間官山小流域流量過程的影響。此外,分流量區(qū)間BMA方法在產(chǎn)匯流機(jī)理相比于產(chǎn)匯流假設(shè)固定的單一水文模型,在復(fù)雜的山區(qū)小流域具有較好的適用性,對于雨洪資料匱乏、產(chǎn)匯流機(jī)制研究不完善的山區(qū)小流域,使用分流量區(qū)間BMA方法可以綜合各模型模擬優(yōu)勢,提高山區(qū)小流域的洪水模擬精度。對于資料比較豐富、能夠有效進(jìn)行分析和推斷產(chǎn)流機(jī)制和山洪形成過程的流域,應(yīng)該更多從物理機(jī)制入手,建立合理的物理過程模型。