王玉麗 ，包為民，沈丹丹，周俊偉，孫逸群

(河海大學(xué)水文水資源學(xué)院，南京 210098)

實(shí)時(shí)洪水預(yù)報(bào)誤差修正主要利用實(shí)時(shí)系統(tǒng)能獲得的觀測信息和一切能利用的其他信息對預(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í)測資料誤差服從正態(tài)分布，但若實(shí)測資料遭到污染出現(xiàn)異常值，在誤差自回歸模型中就轉(zhuǎn)化為誤差中存在異常值，即存在不滿足正態(tài)分布的誤差甚至錯(cuò)誤，再利用傳統(tǒng)的估計(jì)方法得到的結(jié)果就不可信[5]，本文把抗差估計(jì)方法引入到AR模型的參數(shù)估計(jì)中，由抗差能夠抵御異常值的特性來減小異常值對估計(jì)方法的影響，并與傳統(tǒng)方法進(jìn)行比較分析。

1 方法介紹

1.1 遞推最小二乘法

n階的AR模型：

ε(t)=XTtC+δtt=1,2,…,N

(1)

式中：ε(t)為t時(shí)刻的實(shí)測與計(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ǔ)上修正而成的。

1.2 抗差遞推最小二乘法

給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)

針對所研究的問題，此處采用以下抗差權(quán)函數(shù)形式：

(8)

式中：σ為加權(quán)殘差均方差；δ(t)為殘差；k1和k2為常數(shù)，k1=1.5，k2=2.5。

1.3 誤差生成模式

異常誤差生成模式：

(9)

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

2.1 流域概況

閩江流域處于亞熱帶季風(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場洪水資料進(jìn)行分析，流域年平均降水量為1820 mm，年平均徑流深972.0 mm，徑流系數(shù)大于0.4，屬于典型的濕潤地區(qū)，因此使用三水源新安江模型進(jìn)行洪水模擬預(yù)報(bào)是較為合理的[6]。

2.2 應(yīng)用檢驗(yàn)

2.2.1 無異常值情況

當(dāng)實(shí)測流量資料中無異常誤差時(shí)，分別運(yùn)用遞推最小二乘與抗差遞推最小二乘法進(jìn)行計(jì)算，將遞推最小二乘的估計(jì)結(jié)果作為近似真值，記為C0。表1是無異常值情況下抗差遞推最小二乘法計(jì)算的參數(shù)估值相對于C0的偏離程度。并用兩種算的參數(shù)估計(jì)結(jié)果進(jìn)行實(shí)時(shí)修正，用確定性系數(shù)(DC)作為衡量指標(biāo)，結(jié)果見表2。

表1 無異常值情況下抗差遞推最小二乘法結(jié)果的偏離程度Tab.1 Parameter dispersion degree of robust recursive least-squares result without abnormal error

表2 無異常值情況下兩種算法的確定性系數(shù)Tab.2 Deterministic coefficient of the two methods without abnormal error

分析表1，參數(shù)估值的均方差能夠反映參數(shù)估計(jì)的穩(wěn)定程度，估計(jì)參數(shù)的波動(dòng)越大，所對應(yīng)的參數(shù)均方差也越大，因此，可以采用參數(shù)均方差來比較分析兩種算法的好壞。由表1可以看出，當(dāng)資料中不含異常值時(shí)，抗差遞推最小二乘法相對于真值的離散程度很小，說明與真值非常接近，即在觀測值無異常值的情況下，抗差遞推最小二乘法能夠獲得最優(yōu)估值。

分析表2，把兩種算法得到的參數(shù)估值進(jìn)行實(shí)時(shí)修正，以洪水過程確定性系數(shù)進(jìn)行衡量，可以看出，兩種算法的結(jié)果非常接近，并且能獲得較好的實(shí)時(shí)校正效果。

2.2.2 有異常值情況

采用式(9)生成異常誤差，添加到正常的流量資料中，以此產(chǎn)生含有異常誤差的數(shù)據(jù)。通過不同的p和T值情況下的計(jì)算，來分析不同量級和頻率的異常誤差對抗差估計(jì)的影響，本文選取p=2、3、4，T=5、8、10，共9種組合情況。表3表示在這9種情況下兩種算法計(jì)算的參數(shù)估值相對于真值的偏離程度，表4則顯示了這兩種算法所獲得的參數(shù)值進(jìn)行實(shí)時(shí)洪水校正的效果，以確定性系數(shù)(DC)作為衡量指標(biāo)。

表3 不同量級和頻率的異常誤差情況下兩種算法的參數(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 不同量級和頻率的異常誤差情況下兩種算法的確定性系數(shù)Tab.4 Deterministic coefficient of the tow methods with different magnitude and frequency of the abnormal error

這5場洪水的計(jì)算結(jié)果相似，現(xiàn)選取970605次洪水為例來比較兩種算法在不同量級和頻率的異常誤差情況下的結(jié)果，選取p為3、4時(shí)在頻率為5的情況下的兩種算法的效果，其結(jié)果如圖1所示(其中RRLS為抗差遞推最小二乘法，RLS為遞推最小二乘法)，然后再比較兩種量級(p=3、4)下不同頻率(T=5、10)的抗差效果，結(jié)果如圖2所示。

圖1 970605號洪水不同量級異常誤差下的兩種算法的效果(以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號洪水不同頻率異常誤差下的抗差效果(以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可以看出，在不同量級和頻率的異常誤差情況下，遞推最小二乘算法的參數(shù)估計(jì)值嚴(yán)重偏離真值，隨著異常誤差量級的增大和頻率的減小，參數(shù)均方差越大，即估計(jì)效果越差。而抗差遞推最小二乘法依據(jù)抗差能夠抵御異常誤差對參數(shù)估計(jì)的干擾這一特性，獲得參數(shù)均方差比遞推最小二乘的小。

分析表4同樣可以發(fā)現(xiàn)，遞推最小二乘算法的校正效果明顯受異常值得影響，且隨著異常誤量級的增大，校正精度越差，而抗差遞推最小二乘的結(jié)果相對較穩(wěn)定，且效果比遞推最小二乘的結(jié)果要好很多。

分析圖1可以看出，遞推最小二乘的參數(shù)估計(jì)結(jié)果波動(dòng)很大，且嚴(yán)重偏離真值。而抗差遞推最小二乘的結(jié)果和真值比較貼近，但是當(dāng)p=4的結(jié)果比p=3的結(jié)果要好些，這是因?yàn)楫惓Ｕ`差的量級越大，則越容易被識別出來，抗差的效果隨著異常值量級的增大而越好。

分析圖2顯示，同一量級的異常誤差情況下，T=10的效果要比T=5的效果要好，這是因?yàn)殡S著異常值的個(gè)數(shù)增多，抗差準(zhǔn)確探測異常的概率就會降低，相應(yīng)的修正效果也降低。由圖2也可以看出，同一發(fā)生頻率的異常誤差下，誤差量級越大則修正效果越好。

3 結(jié) 語

(1)當(dāng)觀測值中無異常誤差時(shí)，對于AR模型的參數(shù)估計(jì)，利用遞推最小二乘法與抗差遞推最小二乘法都可以獲得理想的值，且兩種算法的結(jié)果相近，將兩種算法獲得的參數(shù)估值進(jìn)行實(shí)時(shí)計(jì)算，可提高實(shí)時(shí)洪水預(yù)報(bào)的精度。

(2)當(dāng)觀測值中含有異常誤差時(shí)，遞推最小二乘法的結(jié)果嚴(yán)重偏離真值，抗差遞推最小二乘法依據(jù)抗差理論的特性，在不同量級和頻率的誤差情況下都能夠抵御異常值的干擾，獲得比較穩(wěn)定可靠的參數(shù)，由此得到的實(shí)時(shí)修正效果也非常理想。

(3)無論觀測數(shù)據(jù)中是否含有異常誤差，抗差遞推最小二乘法獲得的參數(shù)結(jié)果都非常接近優(yōu)值，也能得到精度相對較高的實(shí)時(shí)校正結(jié)果。

(4)可以看出異常值的數(shù)目對抗差估計(jì)的效果是有影響的，對于這一數(shù)目的確定有待進(jìn)一步研究。

□

[1] 包為民.水文預(yù)報(bào)[M].北京:中國水利水電出版,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] 李榮容,吳國堯.水庫入庫流量抗差修正研究[J].中國農(nóng)村水利水電,2008,(11):12-14.

[15] 沈丹丹,包為民,劉可新,等. 馬斯京根匯流參數(shù)抗差估計(jì)研究[J]. 中國農(nóng)村水利水電,2016,(7):72-74.