楊瑞祥，梁　川，*，景　楠,楊來順

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

?

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

楊瑞祥1，2，梁川1，2，*，景楠1，2,楊來順3

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

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

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

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

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

1　研究方法

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

1.1序貫重要性采樣

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

(1)

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

(2)

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

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

(3)

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

(4)

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

1.2重采樣

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

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

(5)

在序貫重要性采樣時，如果Neff小于某個值，就應(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個權(quán)值均等的預(yù)測粒子可近似表示t-1時刻狀態(tài)變量的先驗(yàn)概率密度。在獲得當(dāng)前時刻的觀測值zk-1后，重新計算每個粒子的權(quán)重。與觀測值比較接近的粒子被賦予較大的權(quán)重，而與觀測值相隔較遠(yuǎn)的粒子被賦予較小的權(quán)重，之后將粒子權(quán)重歸一化。

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

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

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

2　模型應(yīng)用

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

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

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

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

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

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

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

Table 1Parameters of the Xinanjiang Model (3 components)

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

2.3實(shí)時校正方案

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

2.3.1校正結(jié)果方案

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

(6)

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

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

2.3.2校正參數(shù)方案

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

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

2.4評價指標(biāo)

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

(7)

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

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

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

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

(8)

(9)

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

3　效果評價

3.1預(yù)報結(jié)果

各場洪水預(yù)報流量過程線見圖1，其中用作率定模型參數(shù)的是前5場洪水。

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

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

3.2結(jié)果評價

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

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

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

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

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

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

4　結(jié)　論

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

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

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

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

[2]安全,范瑞琪,朱禎,等. 構(gòu)建式流域?qū)崟r水文預(yù)報系統(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é)報:自然科學(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]. 中國科學(xué):地球科學(xué),2014(5):1002-1016.

[7]徐興亞,方紅衛(wèi),張?jiān)婪?等. 河道洪水實(shí)時概率預(yù)報模型與應(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]魏國華. 粒子濾波算法研究與實(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] 中華人民共和國國家質(zhì)量監(jiān)督檢驗(yàn)檢疫總局，中國國家標(biāo)準(zhǔn)化管理委員會.水文情報預(yù)報規(guī)范：GB/T 22482-2008[S].北京：中國標(biāo)準(zhǔn)出版社，2009.

[14] 包為民. 水文預(yù)報：第4版[M].北京:中國水利水電出版社,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

國家自然科學(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