王玉麗 ,包為民,沈丹丹,周俊偉,孫逸群
(河海大學(xué)水文水資源學(xué)院,南京 210098)
實(shí)時(shí)洪水預(yù)報(bào)誤差修正主要利用實(shí)時(shí)系統(tǒng)能獲得的觀測(cè)信息和一切能利用的其他信息對(duì)預(yù)報(bào)誤差進(jìn)行實(shí)時(shí)校正,以彌補(bǔ)流域水文模型的不足[1]。實(shí)時(shí)修正的方法有很多,如抗差修正[2]、神經(jīng)網(wǎng)絡(luò)修正[3]、卡爾曼濾波法以及自回歸修正等。
當(dāng)水文預(yù)報(bào)模型采用概念性水文模型,如新安江模型時(shí),實(shí)時(shí)校正模型可采用誤差自回歸( AR)模型[4]。一般實(shí)測(cè)資料誤差服從正態(tài)分布,但若實(shí)測(cè)資料遭到污染出現(xiàn)異常值,在誤差自回歸模型中就轉(zhuǎn)化為誤差中存在異常值,即存在不滿(mǎn)足正態(tài)分布的誤差甚至錯(cuò)誤,再利用傳統(tǒng)的估計(jì)方法得到的結(jié)果就不可信[5],本文把抗差估計(jì)方法引入到AR模型的參數(shù)估計(jì)中,由抗差能夠抵御異常值的特性來(lái)減小異常值對(duì)估計(jì)方法的影響,并與傳統(tǒng)方法進(jìn)行比較分析。
n階的AR模型:
ε(t)=XTtC+δtt=1,2,…,N
(1)
式中:ε(t)為t時(shí)刻的實(shí)測(cè)與計(jì)算流量之間的差值;Xt=ε(t-1),ε(t-2),…,ε(t-N);C為AR模型的參數(shù)矩陣。
當(dāng)流量資料系列存在異常誤差時(shí),就轉(zhuǎn)化為 系列存在著異常值。t時(shí)刻的參數(shù)估計(jì)為:
(2)
Φt=[X1,X2,…,Xt];Yt=[ε(1),ε(2),…,ε(t)]
根據(jù)矩陣的分塊乘法,得到遞推公式為:
(4)
由以上可知,遞推最小二乘算法的思想就是:新的參數(shù)估計(jì)值=舊的參數(shù)估計(jì)值+修正項(xiàng),即新的遞推參數(shù)估值就是在舊的遞推估值的基礎(chǔ)上修正而成的。
給AR模型的誤差系列的每個(gè)值加上初始權(quán)重,初始權(quán)重為單位陣,則式(2)就轉(zhuǎn)換為:
(5)
式中:Wt=diag[w(1),w(2),…,w(t)],W(t)是t時(shí)刻的權(quán)重值。
此處設(shè)P(t)=(ΦTtWtΦt)-1,得到抗差遞推最小二乘法的公式為:
(7)
針對(duì)所研究的問(wèn)題,此處采用以下抗差權(quán)函數(shù)形式:
(8)
式中:σ為加權(quán)殘差均方差;δ(t)為殘差;k1和k2為常數(shù),k1=1.5,k2=2.5。
異常誤差生成模式:
(9)
閩江流域處于亞熱帶季風(fēng)氣候區(qū),流域多年平均降水量 1 724 mm,年內(nèi)分布不均,年際變化大,其中60% 左右的年降水量主要集中在4-6月,豐水年降雨量是枯水年的2~3倍。受氣候的影響,流域內(nèi)多年平均蒸發(fā)量為915.0 mm,夏季蒸發(fā)量較大,冬季蒸發(fā)量較小。水吉流域地處閩江上游地區(qū),流域面積3 767 km2,多年平均氣溫15~20 ℃。本文選用水吉流域1988-1999年之間的18場(chǎng)洪水資料進(jìn)行分析,流域年平均降水量為1820 mm,年平均徑流深972.0 mm,徑流系數(shù)大于0.4,屬于典型的濕潤(rùn)地區(qū),因此使用三水源新安江模型進(jìn)行洪水模擬預(yù)報(bào)是較為合理的[6]。
2.2.1 無(wú)異常值情況
當(dāng)實(shí)測(cè)流量資料中無(wú)異常誤差時(shí),分別運(yùn)用遞推最小二乘與抗差遞推最小二乘法進(jìn)行計(jì)算,將遞推最小二乘的估計(jì)結(jié)果作為近似真值,記為C0。表1是無(wú)異常值情況下抗差遞推最小二乘法計(jì)算的參數(shù)估值相對(duì)于C0的偏離程度。并用兩種算的參數(shù)估計(jì)結(jié)果進(jìn)行實(shí)時(shí)修正,用確定性系數(shù)(DC)作為衡量指標(biāo),結(jié)果見(jiàn)表2。
表1 無(wú)異常值情況下抗差遞推最小二乘法結(jié)果的偏離程度Tab.1 Parameter dispersion degree of robust recursive least-squares result without abnormal error
表2 無(wú)異常值情況下兩種算法的確定性系數(shù)Tab.2 Deterministic coefficient of the two methods without abnormal error
分析表1,參數(shù)估值的均方差能夠反映參數(shù)估計(jì)的穩(wěn)定程度,估計(jì)參數(shù)的波動(dòng)越大,所對(duì)應(yīng)的參數(shù)均方差也越大,因此,可以采用參數(shù)均方差來(lái)比較分析兩種算法的好壞。由表1可以看出,當(dāng)資料中不含異常值時(shí),抗差遞推最小二乘法相對(duì)于真值的離散程度很小,說(shuō)明與真值非常接近,即在觀測(cè)值無(wú)異常值的情況下,抗差遞推最小二乘法能夠獲得最優(yōu)估值。
分析表2,把兩種算法得到的參數(shù)估值進(jìn)行實(shí)時(shí)修正,以洪水過(guò)程確定性系數(shù)進(jìn)行衡量,可以看出,兩種算法的結(jié)果非常接近,并且能獲得較好的實(shí)時(shí)校正效果。
2.2.2 有異常值情況
采用式(9)生成異常誤差,添加到正常的流量資料中,以此產(chǎn)生含有異常誤差的數(shù)據(jù)。通過(guò)不同的p和T值情況下的計(jì)算,來(lái)分析不同量級(jí)和頻率的異常誤差對(duì)抗差估計(jì)的影響,本文選取p=2、3、4,T=5、8、10,共9種組合情況。表3表示在這9種情況下兩種算法計(jì)算的參數(shù)估值相對(duì)于真值的偏離程度,表4則顯示了這兩種算法所獲得的參數(shù)值進(jìn)行實(shí)時(shí)洪水校正的效果,以確定性系數(shù)(DC)作為衡量指標(biāo)。
表3 不同量級(jí)和頻率的異常誤差情況下兩種算法的參數(shù)均方差Tab.3 Parameter Mean Square Error of the tow methods with different magnitude and frequency of the abnormal error
注:V為未抗差,VR為進(jìn)行了抗差。
表4 不同量級(jí)和頻率的異常誤差情況下兩種算法的確定性系數(shù)Tab.4 Deterministic coefficient of the tow methods with different magnitude and frequency of the abnormal error
這5場(chǎng)洪水的計(jì)算結(jié)果相似,現(xiàn)選取970605次洪水為例來(lái)比較兩種算法在不同量級(jí)和頻率的異常誤差情況下的結(jié)果,選取p為3、4時(shí)在頻率為5的情況下的兩種算法的效果,其結(jié)果如圖1所示(其中RRLS為抗差遞推最小二乘法,RLS為遞推最小二乘法),然后再比較兩種量級(jí)(p=3、4)下不同頻率(T=5、10)的抗差效果,結(jié)果如圖2所示。
圖1 970605號(hào)洪水不同量級(jí)異常誤差下的兩種算法的效果(以p=3、4,T=5為例)Fig.1 Tow methods Result of 970605 flood with different magnitude of the abnormal error(take p=3、4,T=5、10 as example )
圖2 970605號(hào)洪水不同頻率異常誤差下的抗差效果(以p=3、4,T=5、10為例)Fig.2 Result of 970605 flood with different frequency of the abnormal error(take p=3、4,T=5、10 as example )
分析表3可以看出,在不同量級(jí)和頻率的異常誤差情況下,遞推最小二乘算法的參數(shù)估計(jì)值嚴(yán)重偏離真值,隨著異常誤差量級(jí)的增大和頻率的減小,參數(shù)均方差越大,即估計(jì)效果越差。而抗差遞推最小二乘法依據(jù)抗差能夠抵御異常誤差對(duì)參數(shù)估計(jì)的干擾這一特性,獲得參數(shù)均方差比遞推最小二乘的小。
分析表4同樣可以發(fā)現(xiàn),遞推最小二乘算法的校正效果明顯受異常值得影響,且隨著異常誤量級(jí)的增大,校正精度越差,而抗差遞推最小二乘的結(jié)果相對(duì)較穩(wěn)定,且效果比遞推最小二乘的結(jié)果要好很多。
分析圖1可以看出,遞推最小二乘的參數(shù)估計(jì)結(jié)果波動(dòng)很大,且嚴(yán)重偏離真值。而抗差遞推最小二乘的結(jié)果和真值比較貼近,但是當(dāng)p=4的結(jié)果比p=3的結(jié)果要好些,這是因?yàn)楫惓U`差的量級(jí)越大,則越容易被識(shí)別出來(lái),抗差的效果隨著異常值量級(jí)的增大而越好。
分析圖2顯示,同一量級(jí)的異常誤差情況下,T=10的效果要比T=5的效果要好,這是因?yàn)殡S著異常值的個(gè)數(shù)增多,抗差準(zhǔn)確探測(cè)異常的概率就會(huì)降低,相應(yīng)的修正效果也降低。由圖2也可以看出,同一發(fā)生頻率的異常誤差下,誤差量級(jí)越大則修正效果越好。
(1)當(dāng)觀測(cè)值中無(wú)異常誤差時(shí),對(duì)于AR模型的參數(shù)估計(jì),利用遞推最小二乘法與抗差遞推最小二乘法都可以獲得理想的值,且兩種算法的結(jié)果相近,將兩種算法獲得的參數(shù)估值進(jìn)行實(shí)時(shí)計(jì)算,可提高實(shí)時(shí)洪水預(yù)報(bào)的精度。
(2)當(dāng)觀測(cè)值中含有異常誤差時(shí),遞推最小二乘法的結(jié)果嚴(yán)重偏離真值,抗差遞推最小二乘法依據(jù)抗差理論的特性,在不同量級(jí)和頻率的誤差情況下都能夠抵御異常值的干擾,獲得比較穩(wěn)定可靠的參數(shù),由此得到的實(shí)時(shí)修正效果也非常理想。
(3)無(wú)論觀測(cè)數(shù)據(jù)中是否含有異常誤差,抗差遞推最小二乘法獲得的參數(shù)結(jié)果都非常接近優(yōu)值,也能得到精度相對(duì)較高的實(shí)時(shí)校正結(jié)果。
(4)可以看出異常值的數(shù)目對(duì)抗差估計(jì)的效果是有影響的,對(duì)于這一數(shù)目的確定有待進(jìn)一步研究。
□
[1] 包為民.水文預(yù)報(bào)[M].北京:中國(guó)水利水電出版,2006.
[2] 包為民,嵇海祥,胡其美,等.抗差理論及在水文學(xué)中的應(yīng)用[J]. 水科學(xué)進(jìn)展,2003,14(4):528-532.
[3] 覃光華,丁 晶.基于人工神經(jīng)網(wǎng)絡(luò)的卡爾曼濾波實(shí)時(shí)校正技術(shù)[J].水力發(fā)電,2002,(11): 9-12.
[4] 包為民,王 浩,趙 超,等.AR模型參數(shù)的抗差估計(jì)研究[J].河海大學(xué)學(xué)報(bào)(自然科學(xué)版),2006,34(3):258-261.
[5] 趙 超,包為民,王葉琴,等.河段匯流參數(shù)抗差估計(jì)研究[J].河海大學(xué)學(xué)報(bào)(自然科學(xué)版),2006,36(1):15-17.
[6] 張章新.閩江流域水文特性分析[J].水文,2000,20(6):55-58.
[7] 趙人俊.流域水文模擬-新安江模型與陜北模型[M].北京:水利電力出版社,1984:110-112.
[8] 包為民,瞿思敏,黃賢慶,等.水文系統(tǒng)抗差權(quán)函數(shù)分析與檢驗(yàn)[J].清華大學(xué)學(xué)報(bào)(自然科學(xué)版),2003,43(8):1 127-1 129.
[10] Li Qian, Bao Weimin,Qian jinglin. An error updating system for real-time flood forecasting based on robust procedure[J].KSCE Journal of Civil Engineering, 2015,19(3):796-803.
[11] 瞿思敏,包為民.實(shí)時(shí)洪水預(yù)報(bào)綜合修正方法初探[J].水科學(xué)進(jìn)展,2003,14(2):167-171.
[12] 瞿思敏.抗差理論在洪水預(yù)報(bào)中的應(yīng)用研究[D].南京:河海大學(xué),2004.
[13] 包為民,李榮容,王 濤,等.波動(dòng)系數(shù)與洪水預(yù)報(bào)抗差效果關(guān)系分析[J]. 水力發(fā)電學(xué)報(bào),2009,3(28):57-61.
[14] 李榮容,吳國(guó)堯.水庫(kù)入庫(kù)流量抗差修正研究[J].中國(guó)農(nóng)村水利水電,2008,(11):12-14.
[15] 沈丹丹,包為民,劉可新,等. 馬斯京根匯流參數(shù)抗差估計(jì)研究[J]. 中國(guó)農(nóng)村水利水電,2016,(7):72-74.