王旭瀅 包為民 王玉麗 游 洋 鐘 華

(1.河海大學(xué) 水文水資源學(xué)院,南京 210098;2.中國(guó)電建集團(tuán)中南勘測(cè)設(shè)計(jì)研究院有限公司,長(zhǎng)沙410014;3.江蘇省水利廳水資源處,南京 210029;4.南京水利科學(xué)研究院,南京 210029)

洪澇災(zāi)害是威脅人類(lèi)生命財(cái)產(chǎn)安全的十大隱患之一,隨著社會(huì)經(jīng)濟(jì)規(guī)模的擴(kuò)大與發(fā)展,提高防洪能力勢(shì)在必行.修建水庫(kù)是防洪防汛工程的典型措施之一,水庫(kù)入庫(kù)洪水的精確預(yù)報(bào)是水庫(kù)優(yōu)化調(diào)度、保證水庫(kù)下游安全的重要前提,因此研究水庫(kù)洪水預(yù)報(bào)方案意義重大.水庫(kù)流域洪水預(yù)報(bào)不同于一般的流域預(yù)報(bào),水庫(kù)大壩的建設(shè)會(huì)改變天然水文規(guī)律,從而給預(yù)報(bào)帶來(lái)一定的復(fù)雜性和挑戰(zhàn)性.受到影響的水文要素有:庫(kù)區(qū)下墊面特性、河道洪水波的演進(jìn)、壩址流量的獲取方式、水文資料序列的一致性發(fā)生變化等.加之洪水預(yù)報(bào)自身的不確定性、預(yù)報(bào)方法的選擇等其它公共因素,使得入庫(kù)洪水預(yù)報(bào)結(jié)果誤差較大,可信度較低[1].這些不確定因素會(huì)反映在反推得到的入庫(kù)流量中,從而對(duì)實(shí)時(shí)校正結(jié)果造成“污染”.

觀測(cè)誤差按照大小可以分為三類(lèi):小誤差、粗差、極值誤差,后兩者對(duì)估值影響最為明顯.傳統(tǒng)誤差修正方法的基礎(chǔ)是最小二乘法,它要求誤差呈正態(tài)分布,然而實(shí)際中誤差分布是未知的,難以概述,這會(huì)造成估值的不準(zhǔn)確.抗差方法在未知誤差具體分布的情況下,針對(duì)實(shí)時(shí)校正結(jié)果中的“污染”,提出了誤差“污染分布”的概念,利用等價(jià)權(quán),構(gòu)造出特定的目標(biāo)函數(shù),從而抵御粗差和極值誤差對(duì)估值產(chǎn)生的影響[2].流量資料的抗差方法已在理想模型中得到驗(yàn)證[3],本文針對(duì)水庫(kù)預(yù)報(bào)問(wèn)題,結(jié)合流量資料的抗差方法,對(duì)水庫(kù)洪水預(yù)報(bào)實(shí)時(shí)校正進(jìn)行改進(jìn).

1 水文模型與實(shí)時(shí)校正方法

根據(jù)本文研究區(qū)域的水文氣象條件,采用新安江模型作為其入庫(kù)洪水預(yù)報(bào)模型.

新安江模型屬于概念性模型,其計(jì)算流程如圖1所示.模型最大特點(diǎn)是“三分”:分單元、分水源、分階段.分單元即根據(jù)流域下墊面的空間分布差異性對(duì)其進(jìn)行分塊;分水源即徑流分水源,分為地表徑流、地下徑流、壤中流,在匯流速度上地表徑流＞壤中流＞地下徑流;分階段即匯流過(guò)程分為坡面匯流、河網(wǎng)匯流.馬斯京根河道匯流演算是通過(guò)聯(lián)立水量平衡方程和槽蓄方程得到的,已知入流過(guò)程,通過(guò)上一時(shí)段出流量推算出下一時(shí)段的出流[4].

圖1 三水源新安江模型計(jì)算流程圖

洪水預(yù)報(bào)實(shí)時(shí)校正的對(duì)象可以是預(yù)報(bào)模型參數(shù)、預(yù)報(bào)誤差以及狀態(tài)變量等.常見(jiàn)的實(shí)時(shí)校正方法有自回歸修正(AR模型)、遞推最小二乘法、卡爾曼濾波等.這些方法能夠?qū)崟r(shí)處理預(yù)報(bào)系統(tǒng)中出現(xiàn)的誤差,并且以這些誤差作為修正的依據(jù),在預(yù)報(bào)系統(tǒng)中作出及時(shí)的糾正反應(yīng)[5].較為常用的是誤差自回歸修正法,即認(rèn)為當(dāng)前時(shí)刻的誤差與此時(shí)刻之前的誤差系列緊密相關(guān),通過(guò)實(shí)測(cè)流量序列與預(yù)報(bào)流量序列建立殘差序列,從而對(duì)殘差進(jìn)行實(shí)時(shí)校正,最后疊加到預(yù)報(bào)序列上.但該方法對(duì)自回歸系數(shù)的估計(jì)是時(shí)不變的,而對(duì)于一個(gè)流域系統(tǒng)來(lái)說(shuō),實(shí)際應(yīng)為時(shí)變非線性系統(tǒng).因此,本文采用遞推最小二乘法,根據(jù)實(shí)時(shí)輸入的信息修正預(yù)報(bào)誤差.

2 抗差修正方法

抗差修正來(lái)源于1960年Tukey[6]提出的兩部分組合而成的污染分布:主體分布、污染分布.他認(rèn)為主體分布是服從正態(tài)分布的,而污染分布則是極值誤差和粗差聚集而成的,其分布無(wú)規(guī)律可循,難以用數(shù)學(xué)語(yǔ)言描述.Huber[7]通過(guò)把α鄰域的分布函數(shù)轉(zhuǎn)換為相應(yīng)的極值函數(shù)ρ進(jìn)行抗差估計(jì).可以說(shuō),抗差估計(jì)是參數(shù)估計(jì)與非參數(shù)估計(jì)的結(jié)合,其目的就是削弱變方差對(duì)估值造成的不良影響.抗差估計(jì)的有效性主要取決于抗差權(quán)函數(shù),目前常用的抗差估計(jì)方法有Huber法[7]、Tukey 雙 權(quán) 法[6]、Hampel估 計(jì)[8]、IGG法[9].

對(duì)于水庫(kù)流量修正,宜采用雙權(quán)法的抗差權(quán)函數(shù).由于入庫(kù)流量反推帶來(lái)粗差,流量過(guò)程線呈上下波動(dòng)的鋸齒狀.為了將其修正成光滑的曲線,首先進(jìn)行二次曲線擬合,對(duì)場(chǎng)次洪水實(shí)測(cè)值進(jìn)行分段擬合處理得到光滑流量Qgh.二次曲線擬合的原則是最小二乘法,即最小平方法.分段的節(jié)點(diǎn)選擇非常重要,關(guān)鍵是要結(jié)合計(jì)算流量的漲落特點(diǎn)來(lái)考慮[10].分段段數(shù)需盡可能多,否則會(huì)造成誤差的增大.但是對(duì)于水庫(kù)實(shí)測(cè)入庫(kù)流量來(lái)說(shuō),資料信息較少,因此不能完全剔除某個(gè)實(shí)測(cè)值.但是為了限制異常值的影響,采用如下方法:建立降權(quán)區(qū),在降權(quán)區(qū)內(nèi)做抗差處理,之外則保留原信息[11].相較于理想模型中人為對(duì)粗差進(jìn)行控制,實(shí)測(cè)信息中波動(dòng)更為隨機(jī),規(guī)律性不明顯,故要特別注意不能“抗錯(cuò)”和“抗過(guò)”.在保留實(shí)測(cè)有效信息的基礎(chǔ)上,對(duì)殘差較小的系列統(tǒng)一賦權(quán)重,只對(duì)殘差較大系列進(jìn)行抗差處理.因此,根據(jù)實(shí)測(cè)流量和光滑流量之間殘差ζ(t)采用如下公式進(jìn)行計(jì)算抗差權(quán)函數(shù)w(t)和權(quán)函數(shù)特征量:

為了最大可能保留有效實(shí)測(cè)信息,選擇k=1.0.式(2)中m為數(shù)據(jù)點(diǎn)個(gè)數(shù),可以認(rèn)為數(shù)據(jù)點(diǎn)等權(quán)重,wi=1,因此每一場(chǎng)洪水σ為定值.

最終修正得到的流量即是光滑流量和實(shí)測(cè)流量按照權(quán)重進(jìn)行分配,計(jì)算公式如下:

其中,Qsc(t)為實(shí)測(cè)流量,Qxz(t)為修正后的光滑抗差流量.

水庫(kù)原入庫(kù)實(shí)測(cè)流量過(guò)程線的波動(dòng)大小可以用波動(dòng)系數(shù)α表示.計(jì)算公式如下:

3 實(shí)例應(yīng)用

3.1 流域概況

白水坑水庫(kù)流域位于浙江省西南部江山港峽口水庫(kù)上游,東經(jīng)118°22'37″,北緯28°15'26″.白水坑水庫(kù)的東邊是本省的衢縣、遂昌縣,向南與福建省浦城縣相鄰,西與江西省玉山縣、廣豐縣接壤,北靠本省的常山縣.白水坑水庫(kù)流域發(fā)源于浙閩贛三省交界的山區(qū),地勢(shì)南高北低.江山港是該流域的干流,屬于錢(qián)塘江水系,江山港與常山港在衢州市雙港口匯合后稱(chēng)為衢江.江山港主流兩岸多山,少平原,河谷兩岸多為山丘.山區(qū)的沉積物土質(zhì)較薄,以砂礫質(zhì)土為主.江山港屬山溪性河流,洪水暴漲暴落,集流時(shí)間短.洪水的主要誘因是臺(tái)風(fēng)和梅雨,因此汛期在5～7月.洪水特征是峰高量大,呈“高胖”狀.洪水過(guò)程一般歷時(shí)3～4 d.白水坑流域?qū)賮啛釒Ъ撅L(fēng)氣候區(qū),流域控制面積是330 km2.多年平均氣溫17.1℃左右,多年平均降水量1 754 mm,但分配不均,3～6月降雨占全年57%以上.

白水坑水庫(kù)為江山港梯級(jí)開(kāi)發(fā)的龍頭水庫(kù),是集防洪、發(fā)電、灌溉、生態(tài)效益于一體的控制型大Ⅱ型水利樞紐.水庫(kù)壩址位于大巒口鄉(xiāng)雪花淤橋附近,距江山市區(qū)65 km,庫(kù)區(qū)集雨面積330 km2,主流長(zhǎng)26 km,多年平均徑流總量4億m3[13].

3.2 參數(shù)率定

選擇白水坑水庫(kù)流域2004～2009年日模資料以及21場(chǎng)次洪資料對(duì)新安江模型以及馬斯京根流量演算參數(shù)進(jìn)行率定,另外選擇2010～2012年9場(chǎng)洪水進(jìn)行驗(yàn)證.日模資料包括逐日蒸發(fā)、逐日降雨以及逐日入庫(kù)流量,次洪模型資料包括次洪流量過(guò)程以及該場(chǎng)洪水所對(duì)應(yīng)的降雨過(guò)程.參數(shù)率定的結(jié)果見(jiàn)表1.

表1 參數(shù)率定結(jié)果

率定期共有21場(chǎng)洪水,四項(xiàng)指標(biāo)全部合格的共有21場(chǎng)洪水,合格率為100%.進(jìn)一步分析模型模擬精度,將次洪模型實(shí)測(cè)徑流深與計(jì)算徑流深、實(shí)測(cè)洪峰流量與計(jì)算洪峰流量建立相關(guān)關(guān)系,如圖2所示,點(diǎn)據(jù)均分布在45°線附近,相關(guān)關(guān)系良好.因此將此次率定的參數(shù)作為模型參數(shù)代入模型檢驗(yàn).

圖2 相關(guān)關(guān)系圖

經(jīng)過(guò)檢驗(yàn)期9場(chǎng)洪水的驗(yàn)證,均模擬合格,表明該套參數(shù)能夠反映白水坑流域水文特性,可以作為預(yù)報(bào)參數(shù).

3.3 抗差修正

選擇2010～2012年波動(dòng)強(qiáng)烈的8場(chǎng)洪水進(jìn)行修正.具體步驟如下:首先采用分段二次曲線擬合得到光滑流量,每一段通過(guò)最小二乘法得到該段的擬合方程,得到連續(xù)的光滑流量過(guò)程線Qgh.其次利用式(1)、(2)計(jì)算抗差權(quán)重w(t)、權(quán)函數(shù)特征量σ,通過(guò)式(4)計(jì)算其波動(dòng)系數(shù)α,估計(jì)其波動(dòng)誤差的影響情況.最終通過(guò)式(3)得到8場(chǎng)洪水的光滑抗差流量Qxz(t),將此作為實(shí)測(cè)流量來(lái)對(duì)預(yù)報(bào)流量進(jìn)行實(shí)時(shí)校正.實(shí)時(shí)校正采用遞推最小二乘法,計(jì)算修正后的殘差從而得到計(jì)算模擬的校正值.

流量殘差均方差從統(tǒng)計(jì)意義的角度反映了預(yù)報(bào)流量相對(duì)于實(shí)測(cè)值的偏離程度,但是原先未經(jīng)抗差的實(shí)測(cè)流量存在著異常值的干擾.因此,采用流量殘差抗差均方差來(lái)檢驗(yàn)抗差效果,進(jìn)一步計(jì)算殘差減小率得到抗差的表現(xiàn)如何.計(jì)算公式如下:

式中,Vkc為流量殘差抗差均方差;Vsx為流量殘差均方差;Qo、Qc分別為光滑抗差流量和抗差后實(shí)時(shí)校正的計(jì)算流量;Q'o,Q'c分別為不帶抗差的實(shí)測(cè)值和其對(duì)應(yīng)的實(shí)時(shí)校正計(jì)算值;η為殘差減小率.

將其與抗差前、后的校正流量進(jìn)行比較,分別計(jì)算其殘差均方差,并計(jì)算殘差減小率,結(jié)果見(jiàn)表2.

表2 抗差計(jì)算結(jié)果

以上結(jié)果表明,對(duì)8場(chǎng)歷史洪水的實(shí)測(cè)流量抗差之后再進(jìn)行校正,均能夠有效減小殘差,抵抗粗差、極值帶來(lái)的不良影響,進(jìn)而有效提高模擬精度.由于水庫(kù)其調(diào)蓄特性,易出現(xiàn)陡漲陡落型洪水,波動(dòng)系數(shù)較大,波動(dòng)系數(shù)越大反映其產(chǎn)生粗差的概率越高.

4 結(jié) 論

本文采用新安江模型和馬斯京根河道匯流模型對(duì)白水坑水庫(kù)歷史洪水進(jìn)行模擬,率定相關(guān)參數(shù),以水文站的實(shí)測(cè)流量以及歷史雨情資料對(duì)新安江模型的參數(shù)進(jìn)行率定和檢驗(yàn),得到適合該流域的一套參數(shù)優(yōu)值,將其作為白水坑水庫(kù)洪水預(yù)報(bào)方案中的默認(rèn)參數(shù)值.針對(duì)入庫(kù)洪水反推時(shí)造成的波動(dòng)誤差,利用抗差原理對(duì)8場(chǎng)“鋸齒”狀洪水進(jìn)行修正之后再進(jìn)行實(shí)時(shí)校正,以消除非正常因素對(duì)校正過(guò)程的不良影響.最終效果明顯,可以消除大部分粗差,大大提高洪水預(yù)報(bào)精度.相較于其它修正方法,抗差計(jì)算過(guò)程簡(jiǎn)便,能夠?yàn)樗畮?kù)預(yù)報(bào)系統(tǒng)提供較為準(zhǔn)確的流量資料,可將其應(yīng)用于水庫(kù)實(shí)時(shí)預(yù)報(bào)系統(tǒng)中.

由于水庫(kù)實(shí)時(shí)洪水的波動(dòng)更為隨機(jī),因此,在建立抗差降權(quán)域、選擇參值數(shù)時(shí)要著重考慮水庫(kù)歷史入庫(kù)洪水出現(xiàn)粗差的頻率以及幅度.目前抗差估計(jì)方法所針對(duì)的主要是線性估值問(wèn)題,對(duì)于非線性函數(shù)的參數(shù)如何估計(jì)還需進(jìn)一步研究.

[1] 李樹(shù)平,張懷生.水庫(kù)洪水預(yù)報(bào)誤差分析[J].東北水利水電,2006,24(262):42-43.

[2] 包為民,嵇海祥,胡其美,等.抗差理論及在水文學(xué)中的應(yīng)用[J].水科學(xué)進(jìn)展,2003,14(4):428-432.

[3] 李佳佳,包為民,劉可新,等.流量資料抗差方法[J].水力發(fā)電,2016,42(4):27-30.

[4] 包為民.水文預(yù)報(bào)[M].3版.北京:中國(guó)水利水電出版社,2005.

[5] 田 雨,雷曉輝,蔣云鐘,等.洪水預(yù)報(bào)實(shí)時(shí)校正技術(shù)研究綜述[J].人民黃河,2011,33(3):25-26,28.

[6] Tukey J W.A Survey of Sampling from Contaminated Distribution in Contributions to Probability and Statistics[M].Stanford,Calif.:Stanford University Press,1960.

[7] Huber P J.Robust Estimation of a Location Parameter[J].Ann Math Statist,1964,35(1):73-101.

[8] Hampel F R.A General Qualitative Definition of Robustness[J].Ann Math Statist,1971,42(6):1887-1896.

[9] 周江文,黃幼才,楊元喜,等.抗差最小二乘法[M].武漢:華中理工大學(xué)出版社,1997.

[10]包為民,李榮容,王 濤,等.波動(dòng)系數(shù)與洪水預(yù)報(bào)抗差效果關(guān)系分析[J].水力發(fā)電學(xué)報(bào),2009,28(3):57-61.

[11]李榮容,吳國(guó)堯.水庫(kù)入庫(kù)流量抗差修正研究[J].中國(guó)農(nóng)村水利水電,2008(11):12-14.

[12]Li Qian,Bao Weimin,et al.An Error Updating System for Real-time Flood Forecasting based on Robust Procedure[J].KSCE JOURNAL OF CIVIL ENGINEERING,2015,19(3):796-803.

[13]郭浩亮,吳士偉,姜新忠.白水坑水庫(kù)泄洪洞設(shè)計(jì)[J].浙江水利科技,2006,143(1):56-58,77.