楊瑞祥，梁　川，*，景　楠,楊來(lái)順

( 1.四川大學(xué) 水利水電學(xué)院，成都 610065; 2. 水力學(xué)與山區(qū)河流開(kāi)發(fā)保護(hù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，成都 610065;3.三明水文水資源勘測(cè)分局，福建 三明 365000)

?

基于粒子濾波同化方法的實(shí)時(shí)洪水預(yù)報(bào)

楊瑞祥1，2，梁川1，2，*，景楠1，2,楊來(lái)順3

( 1.四川大學(xué) 水利水電學(xué)院，成都 610065; 2. 水力學(xué)與山區(qū)河流開(kāi)發(fā)保護(hù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，成都 610065;3.三明水文水資源勘測(cè)分局，福建 三明 365000)

為了充分利用實(shí)時(shí)水文觀測(cè)數(shù)據(jù)，提高預(yù)報(bào)的可靠性，將馬斯京根分段連續(xù)流量演算法與粒子濾波同化方法進(jìn)行耦合，構(gòu)建了基于粒子濾波同化方法的實(shí)時(shí)洪水預(yù)報(bào)模型。采用沙溪干流的永安站至沙縣站區(qū)間6場(chǎng)典型洪水過(guò)程對(duì)上述實(shí)時(shí)洪水預(yù)報(bào)模型進(jìn)行了驗(yàn)證。結(jié)果表明，粒子濾波同化方法對(duì)馬斯京根分段連續(xù)流量演算模型具有較好的適應(yīng)性，由此得到的校正參數(shù)方案的確定性系數(shù)為0.85,與無(wú)實(shí)時(shí)校正結(jié)果方案和校正方案相比均有一定提高，且預(yù)報(bào)精度為甲級(jí)的洪水場(chǎng)次分別增加了2場(chǎng)和1場(chǎng)，使洪水預(yù)報(bào)精度得到了一定的提高。

馬斯京根法；粒子濾波；實(shí)時(shí)洪水預(yù)報(bào)；數(shù)據(jù)同化

實(shí)時(shí)洪水預(yù)報(bào)能利用水情測(cè)報(bào)系統(tǒng)實(shí)時(shí)測(cè)量的水文數(shù)據(jù)，對(duì)原有洪水預(yù)報(bào)方案進(jìn)行校正，以達(dá)到有效提高預(yù)報(bào)精度和防災(zāi)減災(zāi)決策水平的目的[1]。在實(shí)際工作中，常用的實(shí)時(shí)校正方法是最鄰近誤差校正法[2]和數(shù)據(jù)同化方法。后者能將觀測(cè)數(shù)據(jù)融合到模型擬合結(jié)果中，對(duì)系統(tǒng)狀態(tài)和參數(shù)不斷進(jìn)行更新，主要有變分方法和順序同化方法。變分方法的代表是在氣象領(lǐng)域常用的4D-Var算法[3]，順序同化方法的代表是在水文領(lǐng)域得到應(yīng)用的卡爾曼濾波及其優(yōu)化算法[4]。但水文過(guò)程的非線(xiàn)性、非高斯特點(diǎn)使卡爾曼濾波的應(yīng)用受到一定的限制。近年來(lái)得到進(jìn)一步發(fā)展的粒子濾波算法適用于能用狀態(tài)空間模型表示的非線(xiàn)性、非高斯系統(tǒng)，因此在水文領(lǐng)域的應(yīng)用前景非常廣闊。Moradkhani 等在評(píng)估水文模型參數(shù)和狀態(tài)方面，應(yīng)用粒子濾波算法并取得成功[5]。畢海蕓等將粒子濾波算法應(yīng)用于VIC模型中，對(duì)土壤水分進(jìn)行估算[6]。徐興亞等將粒子濾波算法與圣維南方程組結(jié)合，建立了河道洪水實(shí)時(shí)概率預(yù)報(bào)模型[7]。

本文在充分理解標(biāo)準(zhǔn)粒子濾波算法的基礎(chǔ)上，將其與馬斯京根分段連續(xù)流量演算模型進(jìn)行耦合，建立實(shí)時(shí)洪水預(yù)報(bào)模型，并應(yīng)用于福建省沙溪流域干流永安站至沙縣站區(qū)間2002—2010年共6場(chǎng)典型洪水過(guò)程的檢驗(yàn)及實(shí)時(shí)校正。

1　研究方法

粒子濾波同化方法是一種采用蒙特卡羅算法實(shí)現(xiàn)貝葉斯估計(jì)理論的算法，能夠達(dá)到最優(yōu)貝葉斯估計(jì)的效果[8]。其基本思想是利用從狀態(tài)空間選取一組加權(quán)的隨機(jī)樣本粒子，來(lái)實(shí)現(xiàn)對(duì)狀態(tài)的概率密度分布的逼近，然后用樣本均值代替積分運(yùn)算，以獲得狀態(tài)的最小方差估計(jì)。實(shí)際中一般采用序貫重要性采樣方法來(lái)選取隨機(jī)樣本粒子。

1.1序貫重要性采樣

序貫重要性采樣是從采用的重要性密度函數(shù)中生成隨機(jī)的粒子集合，利用最新的觀測(cè)值遞推更新得到當(dāng)前的權(quán)值，將粒子歸一化處理，再將每一個(gè)權(quán)值對(duì)應(yīng)一個(gè)粒子，可得到相應(yīng)的分布。

(1)

重要性密度函數(shù)q(xk|xk-1,z1：k)的選擇會(huì)對(duì)粒子濾波的性能造成很大的影響，抑制退化問(wèn)題的有效方法是選取好的重要性密度函數(shù)，從而減小需要的粒子數(shù)目，提高運(yùn)行速度。常用的重要性函數(shù)有兩種：①最優(yōu)的重要性密度函數(shù)，性能更好但更難實(shí)現(xiàn)；②先驗(yàn)密度函數(shù)，雖然沒(méi)有融入最新的觀測(cè)值，但實(shí)現(xiàn)簡(jiǎn)單[9]。重要性密度函數(shù):

(2)

式中p(xk|xk-1)為先驗(yàn)密度函數(shù)。

則粒子權(quán)值的表達(dá)式為:

(3)

式中p(xk|z1:k)為該時(shí)刻狀態(tài)變量的后驗(yàn)概率分布。

(4)

然而，上述方法的主要缺點(diǎn)是出現(xiàn)粒子匱乏現(xiàn)象。原因是在經(jīng)過(guò)數(shù)次遞推計(jì)算后，只有少數(shù)粒子的權(quán)值較大，其余粒子的權(quán)值太小以至于可忽略不計(jì)，選擇任何重要性密度函數(shù)都會(huì)出現(xiàn)這種問(wèn)題。因此，Gordan等[10]增加了重采樣的步驟。

1.2重采樣

重采樣是在保持過(guò)程中粒子總數(shù)不變的情況下，從狀態(tài)的后驗(yàn)概率密度中重新采樣，將權(quán)重較小的粒子舍棄，保留或復(fù)制權(quán)重較大的粒子，將原來(lái)帶權(quán)重的粒子集映射為新的等權(quán)重的粒子集。

通常采用粒子有效個(gè)數(shù)Neff來(lái)衡量粒子有效的數(shù)量[11]，近似為：

(5)

在序貫重要性采樣時(shí)，如果Neff小于某個(gè)值，就應(yīng)當(dāng)進(jìn)行重采樣。

1.3算法實(shí)現(xiàn)

標(biāo)準(zhǔn)粒子濾波算法的基本公式是式(2)和式(3)[12]。實(shí)現(xiàn)步驟如下：

1)粒子權(quán)值更新。假設(shè)N個(gè)權(quán)值均等的預(yù)測(cè)粒子可近似表示t-1時(shí)刻狀態(tài)變量的先驗(yàn)概率密度。在獲得當(dāng)前時(shí)刻的觀測(cè)值z(mì)k-1后，重新計(jì)算每個(gè)粒子的權(quán)重。與觀測(cè)值比較接近的粒子被賦予較大的權(quán)重，而與觀測(cè)值相隔較遠(yuǎn)的粒子被賦予較小的權(quán)重，之后將粒子權(quán)重歸一化。

2)粒子重采樣。對(duì)粒子進(jìn)行復(fù)制，權(quán)值越大則復(fù)制的次數(shù)越多，權(quán)值越小則復(fù)制的次數(shù)越少，權(quán)值過(guò)小的粒子直接被舍棄。經(jīng)過(guò)重采樣之后的粒子權(quán)重重新被均等地設(shè)置為1/N。

3)預(yù)測(cè)下一時(shí)刻狀態(tài)變量。應(yīng)用系統(tǒng)的狀態(tài)轉(zhuǎn)移方程對(duì)t時(shí)刻每個(gè)粒子的狀態(tài)進(jìn)行預(yù)測(cè)，得到t時(shí)刻的觀測(cè)數(shù)據(jù)，據(jù)此更新粒子權(quán)重，具體做法同1)。

4)重復(fù)執(zhí)行以上步驟，直到所有時(shí)刻運(yùn)行完畢。

2　模型應(yīng)用

2.1區(qū)域概況及洪水分析

永安(二)水文站與沙縣(石橋)水文站兩站間流域面積2 619 km2，河段全長(zhǎng)78 km，洪水傳播時(shí)間為6 h。本次共收集兩站區(qū)間2002—2010年6場(chǎng)典型洪水過(guò)程，第4場(chǎng)洪水用作檢驗(yàn)，其余5場(chǎng)洪水用作率定模型參數(shù)。原始資料觀測(cè)時(shí)間間隔為1～5 h，為了達(dá)到預(yù)報(bào)要求，將其插值得到時(shí)間間隔為1 h的流量觀測(cè)數(shù)據(jù)。經(jīng)過(guò)分析可得，6場(chǎng)洪水中第5、6兩場(chǎng)為雙峰洪水，第5場(chǎng)洪水的第2個(gè)洪峰、第6場(chǎng)洪水的第1個(gè)洪峰峰現(xiàn)時(shí)間上游永安站不早于下游沙縣站，且在同一場(chǎng)洪水中，2個(gè)洪峰相對(duì)大小相反。這些實(shí)際情況客觀上造成預(yù)報(bào)較為困難。

2.2無(wú)實(shí)時(shí)校正方案

區(qū)間入流采用三水源新安江模型，模型參數(shù)率定結(jié)果見(jiàn)表1。

考慮河道平均傳播時(shí)間作為預(yù)見(jiàn)期，取6 h，區(qū)間雨量站降雨過(guò)程及永安站入流過(guò)程由相關(guān)方案展延獲得。

馬斯京根分段連續(xù)演算法參數(shù)率定成果為：蓄量流量關(guān)系曲線(xiàn)的坡度K=6，流量比重系數(shù)x=0.4，時(shí)間間隔Δt=1h。則馬斯京根分段連續(xù)流量演算法的參數(shù)Kl=Δt=1h，河段數(shù)n=6，x1=-0.1,計(jì)算得C0=0.375,C1=0.250,C2=0.375。

表1新安江模型(三水源)參數(shù)

Table 1Parameters of the Xinanjiang Model (3 components)

序號(hào)參數(shù)參數(shù)值序號(hào)參數(shù)參數(shù)值1WM1209IM0.452WUM1510SM50.453WLM8511EX1.154K0.3312KG0.115B0.5013KI0.796C0.2014F26197CI0.8815LAG58CG0.3316CS0.95

2.3實(shí)時(shí)校正方案

實(shí)時(shí)校正方案采用校正結(jié)果和校正參數(shù)兩種方案進(jìn)行比較，均實(shí)現(xiàn)了滾動(dòng)預(yù)報(bào)。

2.3.1校正結(jié)果方案

校正結(jié)果方案是最直接的實(shí)時(shí)校正方案，通過(guò)最鄰近誤差校正方法實(shí)現(xiàn)。最鄰近誤差指的是與當(dāng)前預(yù)報(bào)所依據(jù)的最鄰近時(shí)刻的實(shí)測(cè)值與預(yù)報(bào)值之差，本文所依據(jù)的最鄰近時(shí)刻與當(dāng)前預(yù)報(bào)時(shí)刻之差就是預(yù)見(jiàn)期。校正后的預(yù)報(bào)值采用下式計(jì)算：

(6)

式中Qu,t+1為校正后的預(yù)報(bào)值；Qc,t+6為校正前的預(yù)報(bào)值；εt為最臨近誤差。

該法的優(yōu)勢(shì)是簡(jiǎn)便易行，便于預(yù)報(bào)人員根據(jù)實(shí)際情況具體操作，缺點(diǎn)是容易造成校正后的預(yù)報(bào)值發(fā)生劇烈振蕩，且可能降低預(yù)報(bào)精度。

2.3.2校正參數(shù)方案

校正狀態(tài)方案是通過(guò)粒子濾波方法，根據(jù)預(yù)報(bào)起始時(shí)刻到當(dāng)前時(shí)刻的實(shí)測(cè)流量，優(yōu)化模型的參數(shù)，據(jù)此向前預(yù)報(bào)。當(dāng)有新的實(shí)測(cè)數(shù)據(jù)加入，則重新優(yōu)化參數(shù)，繼續(xù)滾動(dòng)向前預(yù)報(bào)。該方案需要設(shè)定粒子數(shù)量、狀態(tài)方差和觀測(cè)方差3個(gè)參數(shù)。

粒子數(shù)量越多，計(jì)算效果理論上越好，但是計(jì)算時(shí)間增加。狀態(tài)方差和觀測(cè)方差越大，校正后的流量過(guò)程線(xiàn)偏離校正前的過(guò)程線(xiàn)的幅度越大，反之，則與校正前的過(guò)程線(xiàn)接近。經(jīng)過(guò)多次率定，粒子數(shù)量取200，在保證一定精度的情況下不至于導(dǎo)致計(jì)算時(shí)間過(guò)長(zhǎng)；該方案有兩個(gè)變化的參數(shù)Kl和xl，因此狀態(tài)方差取0.03和0.015，觀測(cè)方差取100，使得能從流量過(guò)程線(xiàn)上看出有比較明顯的校正，又不至于出現(xiàn)大的震蕩現(xiàn)象。若采用校正參數(shù)C0、C1、C2的方法，反算回Kl、xl會(huì)發(fā)現(xiàn)取值很不合理，失去物理意義，因此不予采用。

2.4評(píng)價(jià)指標(biāo)

根據(jù)《水文情報(bào)預(yù)報(bào)規(guī)范》(GB/T 22482-2008)[13]規(guī)定的洪水預(yù)報(bào)確定性系數(shù)和許可誤差標(biāo)準(zhǔn)對(duì)上述3種方案進(jìn)行評(píng)價(jià)。確定性系數(shù)DC按下式計(jì)算：

(7)

預(yù)報(bào)項(xiàng)目的精度按確定性系數(shù)的大小分為3個(gè)等級(jí)，精度等級(jí)按表2確定。

降雨徑流預(yù)報(bào)以實(shí)測(cè)洪峰流量的15%作為洪峰預(yù)報(bào)許可誤差；過(guò)程預(yù)報(bào)以90%的預(yù)報(bào)區(qū)間表示，精度指標(biāo)選取觀測(cè)值落入預(yù)報(bào)90%區(qū)間內(nèi)的頻率。

表2　預(yù)報(bào)項(xiàng)目精度等級(jí)

校正后模型效果DCt就是模型計(jì)算加上實(shí)時(shí)信息進(jìn)行誤差校正的總預(yù)報(bào)效果，校正效果DCu就是相對(duì)于原模型誤差的效果[14]。效果定量評(píng)價(jià)系數(shù)分別為：

(8)

(9)

式中yci,u為實(shí)時(shí)信息進(jìn)行誤差校正的預(yù)報(bào)總流量。

3　效果評(píng)價(jià)

3.1預(yù)報(bào)結(jié)果

各場(chǎng)洪水預(yù)報(bào)流量過(guò)程線(xiàn)見(jiàn)圖1，其中用作率定模型參數(shù)的是前5場(chǎng)洪水。

圖1　沙縣站洪水預(yù)報(bào)結(jié)果

由圖1可見(jiàn)，從洪水過(guò)程線(xiàn)來(lái)看，校正結(jié)果方案得到的洪水過(guò)程線(xiàn)震蕩較大，校正參數(shù)方案和無(wú)校正方案得到的洪水過(guò)程線(xiàn)較為接近，但前者與實(shí)測(cè)流量點(diǎn)據(jù)的配合較后者好。從洪峰流量和峰現(xiàn)時(shí)間上看，校正結(jié)果方案與實(shí)測(cè)值偏離最大，校正參數(shù)方案和無(wú)校正方案對(duì)于不同場(chǎng)次的洪水來(lái)說(shuō)各有優(yōu)劣。出現(xiàn)上述現(xiàn)象的原因是預(yù)報(bào)有6 h的預(yù)見(jiàn)期，如果在這6 h實(shí)測(cè)洪水過(guò)程線(xiàn)突然變陡或變緩，則會(huì)導(dǎo)致預(yù)報(bào)模型來(lái)不及響應(yīng)，影響預(yù)報(bào)精度。

3.2結(jié)果評(píng)價(jià)

3種預(yù)報(bào)方案率定期、檢驗(yàn)期和總的預(yù)報(bào)精度分別見(jiàn)表3、表4、表5。

表3　率定期洪水預(yù)報(bào)精度

表4　檢驗(yàn)期洪水預(yù)報(bào)精度

表5　全部洪水的總預(yù)報(bào)精度

由表5可見(jiàn)，校正結(jié)果方案在確定性系數(shù)和洪峰合格率指標(biāo)上均不如無(wú)校正方案，綜合來(lái)看降低了無(wú)校正方案的精度。而相比于無(wú)校正方案，校正參數(shù)方案在各項(xiàng)評(píng)價(jià)指標(biāo)上均有提高，說(shuō)明粒子濾波同化方法的應(yīng)用能夠達(dá)到提高預(yù)報(bào)精度的目的。值得注意的是，預(yù)報(bào)值落在90%的區(qū)間內(nèi)的頻率僅有60%左右，說(shuō)明過(guò)程預(yù)報(bào)精度有待進(jìn)一步提高。

總之，本次預(yù)報(bào)效果從好到壞的排序是校正參數(shù)方案>無(wú)校正方案>校正結(jié)果方案。

4　結(jié)　論

1)首次提出粒子濾波同化方法與馬斯京根分段連續(xù)流量演算法的耦合模型，該模型能夠根據(jù)實(shí)時(shí)輸入的實(shí)測(cè)流量值，校正預(yù)報(bào)誤差，實(shí)現(xiàn)了對(duì)預(yù)報(bào)參數(shù)的時(shí)段動(dòng)態(tài)更新，達(dá)到提高實(shí)時(shí)校正精度的目的。

2)以沙溪干流永安站至沙縣站區(qū)間為例，根據(jù)流域典型洪水的資料，率定預(yù)報(bào)模型參數(shù)并采用校正結(jié)果和校正參數(shù)兩種實(shí)時(shí)校正方案進(jìn)行校正。計(jì)算結(jié)果表明，無(wú)校正方案的確定性系數(shù)為0.83，采用校正結(jié)果方案的確定性系數(shù)為0.81，采用校正參數(shù)方案的確定性系數(shù)為0.85，且預(yù)報(bào)精度為甲級(jí)的洪水場(chǎng)次較無(wú)校正方案增加了2場(chǎng)。說(shuō)明所采用的粒子濾波同化方法的實(shí)時(shí)洪水預(yù)報(bào)模型能在一定程度上提高預(yù)報(bào)精度。

3)由于標(biāo)準(zhǔn)粒子濾波存在的粒子退化和重要性密度函數(shù)的選擇問(wèn)題，采用更有效的重采樣算法和重要性密度函數(shù)是今后研究的重點(diǎn)。

[1]陸波. 流域水文模型與卡爾曼濾波耦合實(shí)時(shí)洪水預(yù)報(bào)研究[D].河海大學(xué),2006.

[2]安全,范瑞琪,朱禎,等. 構(gòu)建式流域?qū)崟r(shí)水文預(yù)報(bào)系統(tǒng)關(guān)鍵技術(shù)研究[J]. 水文,2011(4):6-11.

[3]Dimet F X L, Talagr and O. Variational algorithms for analysis and assimilation o1 meteorological observations; theoretical aspects[J]. Tellus A,1986,38(2):97-110.

[4]王文,寇小華. 水文數(shù)據(jù)同化方法及遙感數(shù)據(jù)在水文數(shù)據(jù)同化中的應(yīng)用進(jìn)展[J]. 河海大學(xué)學(xué)報(bào):自然科學(xué)版,2009(5):556-562.

[5]Moradkhani H, Hsu K-L, Gupta H, et al. Uncertainty assessment of hydrologic model states and parameters; sequential data assimilation using the particle filter[J]. Water Resources Research, 2004,41(5):W05012. doi: 10. 1029/2004WR003604.

[6]畢海蕓,馬建文,秦思嫻,等. 基于殘差重采樣粒子濾波的土壤水分估算和水力參數(shù)同步優(yōu)化[J]. 中國(guó)科學(xué):地球科學(xué),2014(5):1002-1016.

[7]徐興亞,方紅衛(wèi),張?jiān)婪?等. 河道洪水實(shí)時(shí)概率預(yù)報(bào)模型與應(yīng)用[J]. 水科學(xué)進(jìn)展,2015(3):356-364.

[8]畢海蕓,馬建文. 粒子濾波算法在數(shù)據(jù)同化中的應(yīng)用研究進(jìn)展[J]. 遙感技術(shù)與應(yīng)用,2014(5):701-710.

[9]魏國(guó)華. 粒子濾波算法研究與實(shí)現(xiàn)[D].北京：電子科技大學(xué),2015.

[10] Gordon N J, Salmond D J, Smith A F M. Novel approach to nonlinear/non-gaussian bayesian state estimation[J]. IEEE Proceedings on Radar and Signal Processing, 1993,140(2):107-113.

[11] Bei Cao, Cai Wen Ma. A fine resampling algorithm for general particle filters[C]//4th International Congress on Image and Signal Processing,2011.

[12] 朱志宇.粒子濾波算法及其應(yīng)用[M].北京:科學(xué)出版社,2010

[13] 中華人民共和國(guó)國(guó)家質(zhì)量監(jiān)督檢驗(yàn)檢疫總局，中國(guó)國(guó)家標(biāo)準(zhǔn)化管理委員會(huì).水文情報(bào)預(yù)報(bào)規(guī)范：GB/T 22482-2008[S].北京：中國(guó)標(biāo)準(zhǔn)出版社，2009.

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

Real-time flood forecasting based on particle filter assimilation method

YANG Rui-Xiang1,2, LIANG Chuan1,2,*, JING Nan1,2, YANG Lai-Shun3

(1.CollegeofWaterResourceandHydropower,SichuanUniversity,Chengdu610065,China; 2.StateKeyLaboratoryofHydraulicsandMountainRiverEngineering,SichuanUniversity,Chengdu610065,China;3.HydrologicandWaterResourcesSurveyBureau,SanmingBranch，Sanming365000,Fujian，China)

To make full use of real-time hydrological data and improve reliability of forecasting, Muskingum piecewise continuous flow calculation model and particle filter assimilation method are coupled to build the real-time flood forecasting model based on particle filter data assimilation method. The above real-time flood forecasting model is verified by using interval six typical flood process from Yongan Station to Shaxian Station of Shaxi mainstream. The results show that the particle filter assimilation method has better adaptability to the Muskingum piecewise continuous flow calculation model. Based on this, the determined coefficient of the scheme of correcting parameter is 0.85. Compared with the scheme of non-real-time correction and correcting result, the determined coefficient is improved and the number of precision grade A of forecast are respectively increased by 2 and 1. The flood forecasting accuracy has been improved to some extent.

Muskingum method; particle filter; real-time flood forecasting; data assimilation

10.13524/j.2095-008x.2016.03.033

2016-07-29

國(guó)家自然科學(xué)基金資助項(xiàng)目(61006403)；水文學(xué)及水資源博士點(diǎn)基金資助項(xiàng)目(20130181110045)

楊瑞祥(1994-)，男，福建三明人，碩士研究生，研究方向: 水文水資源，E-mail: smsxfypy@qq.com; *通訊作者:梁川(1957-)，男，四川雅安人，教授，博士，博士研究生導(dǎo)師，研究方向: 水文水資源與水環(huán)境，E-mail: lshester@sohu.com。

P338

A

2095-008X(2016)03-0001-06