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

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

?

基于粒子濾波同化方法的實時洪水預報

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

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

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

馬斯京根法；粒子濾波；實時洪水預報；數據同化

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

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

1　研究方法

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

1.1序貫重要性采樣

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

(1)

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

(2)

式中p(xk|xk-1)為先驗密度函數。

則粒子權值的表達式為:

(3)

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

(4)

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

1.2重采樣

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

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

(5)

在序貫重要性采樣時，如果Neff小于某個值，就應當進行重采樣。

1.3算法實現

標準粒子濾波算法的基本公式是式(2)和式(3)[12]。實現步驟如下：

1)粒子權值更新。假設N個權值均等的預測粒子可近似表示t-1時刻狀態(tài)變量的先驗概率密度。在獲得當前時刻的觀測值zk-1后，重新計算每個粒子的權重。與觀測值比較接近的粒子被賦予較大的權重，而與觀測值相隔較遠的粒子被賦予較小的權重，之后將粒子權重歸一化。

2)粒子重采樣。對粒子進行復制，權值越大則復制的次數越多，權值越小則復制的次數越少，權值過小的粒子直接被舍棄。經過重采樣之后的粒子權重重新被均等地設置為1/N。

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

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

2　模型應用

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

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

2.2無實時校正方案

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

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

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

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

Table 1Parameters of the Xinanjiang Model (3 components)

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

2.3實時校正方案

實時校正方案采用校正結果和校正參數兩種方案進行比較，均實現了滾動預報。

2.3.1校正結果方案

校正結果方案是最直接的實時校正方案，通過最鄰近誤差校正方法實現。最鄰近誤差指的是與當前預報所依據的最鄰近時刻的實測值與預報值之差，本文所依據的最鄰近時刻與當前預報時刻之差就是預見期。校正后的預報值采用下式計算：

(6)

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

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

2.3.2校正參數方案

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

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

2.4評價指標

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

(7)

預報項目的精度按確定性系數的大小分為3個等級，精度等級按表2確定。

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

表2　預報項目精度等級

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

(8)

(9)

式中yci,u為實時信息進行誤差校正的預報總流量。

3　效果評價

3.1預報結果

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

圖1　沙縣站洪水預報結果

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

3.2結果評價

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

表3　率定期洪水預報精度

表4　檢驗期洪水預報精度

表5　全部洪水的總預報精度

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

總之，本次預報效果從好到壞的排序是校正參數方案>無校正方案>校正結果方案。

4　結　論

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

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

3)由于標準粒子濾波存在的粒子退化和重要性密度函數的選擇問題，采用更有效的重采樣算法和重要性密度函數是今后研究的重點。

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

[2]安全,范瑞琪,朱禎,等. 構建式流域實時水文預報系統(tǒng)關鍵技術研究[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]王文,寇小華. 水文數據同化方法及遙感數據在水文數據同化中的應用進展[J]. 河海大學學報:自然科學版,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]畢海蕓,馬建文,秦思嫻,等. 基于殘差重采樣粒子濾波的土壤水分估算和水力參數同步優(yōu)化[J]. 中國科學:地球科學,2014(5):1002-1016.

[7]徐興亞,方紅衛(wèi),張岳峰,等. 河道洪水實時概率預報模型與應用[J]. 水科學進展,2015(3):356-364.

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

[9]魏國華. 粒子濾波算法研究與實現[D].北京：電子科技大學,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] 朱志宇.粒子濾波算法及其應用[M].北京:科學出版社,2010

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

[14] 包為民. 水文預報：第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

國家自然科學基金資助項目(61006403)；水文學及水資源博士點基金資助項目(20130181110045)

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

P338

A

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