基于集合卡爾曼濾波的實時校正方法

顧爐華,賴錫軍

(1.中國科學院南京地理與湖泊研究所，江蘇 南京 210008; 2.中國科學院大學資源與環(huán)境學院，北京 101407)

基于集合卡爾曼濾波的實時校正方法

顧爐華1,2,賴錫軍1

(1.中國科學院南京地理與湖泊研究所，江蘇 南京 210008; 2.中國科學院大學資源與環(huán)境學院，北京 101407)

為減少非恒定水流計算中的不確定性,基于集合卡爾曼濾波提出多變量交替校正的方法。該方法通過交替校正水位和流量,避開了濾波過程中的大矩陣計算,實現(xiàn)了利用觀測信息直接校正非恒定流狀態(tài)的目的;同時,應用尺度轉換方法提高水位濾波精度。數(shù)值試驗重點考察了觀測誤差和水位變換系數(shù)對模型計算精度的影響。結果表明：觀測誤差越小,模型的計算精度越高;水位尺度變換系數(shù)能顯著增強多變量交替校正方法的效果,變換系數(shù)越大,計算精度越高;基于集合卡爾曼濾波的多變量交替校正方法具有良好的校正性能,能顯著提高河道水流的預報精度。

水動力學模型;集合卡爾曼濾波;非恒定流;實時校正技術;洪水預報

如何將水動力學模型與實時校正技術結合起來提高河網(wǎng)水情的仿真與預報精度,是當前洪水預報領域的重要課題。賴錫軍等[1]基于確定性理論框架給出了非恒定流反問題計算的變分法;而基于統(tǒng)計估值理論的卡爾曼濾波方法(Kalman filter,KF)憑借在線估計時間短、存儲量小、適合實時處理和計算機計算等優(yōu)點成為該領域的研究熱點。葛守西等[2]考察了不同狀態(tài)變量對校正效果的影響,確定了誤差協(xié)方差陣的最優(yōu)結構形式和濾波器最優(yōu)參數(shù)值;周全等[3-5]提出了交替校正方法,分別建立水位與流量的狀態(tài)方程,進行交替濾波計算;吳曉玲等[6]基于交替校正方法,通過濾波增益信息研究了模型誤差的空間分布,改善了動態(tài)噪聲干擾矩陣;Madsen等[7]同樣基于增益信息的空間分布,提出適用于觀測延遲或缺失的優(yōu)化方法?；贙F,Evensen[8]提出了集合卡爾曼濾波方法(ensemble Kalman filter,EnKF),通過對集合成員的統(tǒng)計來計算誤差方差陣,此方法具有求解過程簡單、預報精度高、適用于強非線性系統(tǒng)等優(yōu)點,被廣泛地應用于大氣、海洋、水文等領域的數(shù)據(jù)同化[9-12]。如賴錫軍[13]將EnKF引入到水量預測和洪水預報中來,提出了水動力模型與EnKF耦合的實時多變量分析方法;陳一帆等[14]將EnKF運用于河網(wǎng)的水情數(shù)據(jù)同化。

本文擬建立基于集合卡爾曼濾波的多變量交替校正方法。通過交替校正水位和流量,充分挖掘隱藏在觀測數(shù)據(jù)中的有效信息,避免濾波過程中的大矩陣計算,并引入水位尺度變換系數(shù),探究其對濾波效果的影響。

1 一維非恒定流基本方程

一維非恒定流的基本方程,又稱為圣維南方程,可寫成:

(1)

(2)

式中:Z為水位,m;Q為流量,m3/s;K為流量模數(shù)，m3/s;q為單位河長旁側入流,m2/s;A為過水斷面面積,m2;α為動量校正系數(shù);g為重力加速度,m/s2;x為沿水流方向距離,m;t為時間,s;B為水面寬度,m。

采用Preissmann四點加權隱格式離散上述兩式得到:

A1jΔQj+B1jΔZj+C1jΔQj+1+D1jΔZj+1=E1j(3)

A2jΔQj+B2jΔZj+C2jΔQj+1+D2jΔZj+1=E2j(4)

式中符號含義見文獻[16]。對于緩流,在已知上下游邊界條件的情況下,可用追趕法求出各系數(shù),進而求得河道各斷面每一時刻的水位值和流量值。

2 集合卡爾曼濾波

2.1 濾波原理

EnKF的基本思想是利用蒙特卡羅方法構造一個背景集合,集合元素以抽樣的方式從系統(tǒng)狀態(tài)中獲取,背景集合的樣本協(xié)方差視為預報誤差方差,運行濾波器對背景集合進行更新,得到分析集合,分析集合的均值視為狀態(tài)的最佳估計,分析集合的樣本協(xié)方差視為校正誤差方差。分析集合通過模型傳遞,得到下一時刻的背景集合。EnKF用集合統(tǒng)計的方法估算真實值,具有精度高、計算量小、適用范圍廣等優(yōu)點。

2.2 計算過程

首先,建立狀態(tài)方程和量測方程:

Xk=AXk-1+BUk+Wk

(5)

Zk=Hk+Vk

(6)

式中:Xk、Xk-1分別為k時刻、k-1時刻系統(tǒng)的狀態(tài)量;Uk為k時刻系統(tǒng)的控制量;A、B為狀態(tài)系統(tǒng)參數(shù)矩陣;Zk為k時刻的觀測量;Hk為量測系統(tǒng)的觀測矩陣;Wk、Vk分別為過程和測量的噪聲,它們各自的協(xié)方差為Q和R。

狀態(tài)變量預測:

Xk|k-1=AXk-1|k-1+BUk

(7)

集合誤差協(xié)方差:

(8)

采用下式求增益矩陣：

(9)

狀態(tài)變量更新：

Xk|k=Xk|k-1+K(DK-HkXk|k-1)

(10)

式中：Xk|k為根據(jù)k時刻觀測得到的分析解；DK為k時刻觀測值。

集合誤差協(xié)方差更新：

(11)

式中：Pa為更新后的集合誤差協(xié)方差；Pf為預測集合誤差協(xié)方差；Re為觀測誤差協(xié)方差。

3 多變量交替校正方法

交替校正方法是指針對不同類型狀態(tài)量分別建立濾波器方程,進行交替濾波計算以達到不同狀態(tài)量之間協(xié)同變化的方法。由周全等[3-5]首先提出并應用到一維水動力模型中,驗證了方法的可行性。吳曉玲等[6]基于卡爾曼濾波交替校正方法,從模型誤差空間分布的角度對原方法進行了改進,使得非實測斷面也有很高的校正精度;并將卡爾曼濾波多變量交替校正方法和多變量綜合校正方法進行對比,得出交替校正的效果略優(yōu)于綜合校正的結論。在前人研究的基礎上,本文將交替校正的思路拓展應用于EnKF,探究基于EnKF的一維非恒定流多變量交替校正方法的可行性。該方法的主要思想是構建兩組平行計算的集合卡爾曼濾波器。預報階段兩組狀態(tài)變量獨立計算,在有觀測融入的時刻分別運行水位濾波器和流量濾波器得到分析解,再將分析解賦給水位和流量作為下一步計算的初始條件,重新給集合矩陣賦值并增加擾動,循環(huán)計算。鑒于水位數(shù)量級較小,不能充分體現(xiàn)EnKF的性能,引入尺度變換系數(shù)M,在濾波前對水位進行尺度變換。多變量交替校正計算流程如圖1所示。

圖1 多變量交替校正計算流程

4 數(shù)值試驗

4.1 案例選擇

選擇一全長為60 km的河道進行測試,河道斷面為梯形,上游是流量邊界,流量-時間關系為Q=160+320texp[-0.5(t-3)-exp[-0.5(t-3)],下游為給定的水位邊界。計算河道共61個斷面,空間步長取1 km,初始時刻各斷面流量為160 m3/s。觀測站點分別距離上游11 km、23 km、35 km和47 km,觀測間隔為0.5 h。

4.2 試驗設計

案例通過對比用水量模型和基于集合卡爾曼濾波的多變量交替校正方法的計算結果,來驗證集合卡爾曼濾波模型的準確性,并通過改變觀測對象的誤差標準差和水位尺度變換系數(shù)來探究該方法的校正性能。由于很多河流上游監(jiān)測點很少乃至沒有,造成上游流量過程誤差較大。假定上游流量過程為主要誤差源,取其為真實值的80%作為上游邊界條件,在給定不同觀測精度和水位尺度變換系數(shù)條件下,用EnKF多變量交替校正方法進行分析。試驗選取1%、3%、5%、9%這4種觀測誤差和5、10、20、30和50這5種變換系數(shù)作為控制變量,試驗方案如表1所示。

表1 不同變換系數(shù)和誤差標準差的試驗方案

初始集合是在恒定流條件基礎上疊加相應的擾動得到,表達式如下:

Uk=Uk(1+Nkσu)

(12)

式中:Nk為第k個集合的偽光滑隨機擾動場,每個斷面的集合服從均值為0、差為1的正態(tài)分布,生成方式參考文獻[8];σu為模型誤差的相對標準差。在設置初始集合時取σu=0.1,在設置每個計算步的模型隨機誤差時取σu=2.5×10-4。每組試驗集合數(shù)均取100,試驗中明渠水流運動具體計算模型可見文獻[15]。

4.3 試驗結果及分析

4.3.1 觀測誤差對EnKF的影響

選取距河道起點11 km和47 km處的斷面作為研究對象,計算不同條件下的分析值、預報值和真值,得到流量隨時間的變化曲線如圖2所示。

圖2 流量過程比較

圖2(a)顯示隨著觀測誤差的增加,分析值與真值差距變大,分析值的精度越來越低,這是因為在上游邊界條件不準確時,模型的精度主要取決于觀測的精度,觀測精度越高,模型精度也越高。圖2(b)中各種精度下的分析值與真值已經(jīng)相差無幾,體現(xiàn)了很好的校正效果。對比兩圖能夠發(fā)現(xiàn)距河道起點11 km的斷面的校正效果明顯差于距河道起點47 km的斷面,這是因為模型主要的不確定性來自上游來流,接近上游的斷面在濾波時受上游邊界條件影響較大;而下游邊界是準確的,對濾波中的下游斷面施加影響從而使得校正精度很高。

為更好地比較不同觀測精度下模型的校正性能,引入均方誤差和殘差,其計算公式見文獻[13]。

圖3顯示了不同觀測誤差下流量過程均方誤差隨時間變化的趨勢,總體來說觀測精度越高,均方誤差越小,即校正效果越好,這與上面的結論是一致的。同時可以看出,均方誤差的變化與上游來流量總體一致,而且隨時間呈鋸齒狀變化。表明在每一個有觀測融入的時刻,模型能夠調整運行軌跡,使分析值接近真值。圖4是觀測誤差與所有計算時刻的均方誤差的平均值即殘差的關系曲線,可以看到殘差與觀測誤差值正相關,這也從另一個角度驗證了觀測精度越高,校正效果越好的結論。

圖3 不同觀測誤差下的均方誤差

圖4 觀測誤差與殘差關系曲線

4.3.2 水位尺度變換系數(shù)對EnKF的影響

因為水位量級較小,難以充分發(fā)揮EnKF的性能,本文設置了水位尺度變換系數(shù)M,考察其對EnKF校正效果的影響。取中間斷面即距河道起點30 km處流量過程作為研究對象,在5%觀測誤差條件下進行數(shù)值試驗,結果見圖5。

圖5 中間斷面流量過程

圖5顯示隨著水位尺度變換系數(shù)M的增大,分析值越來越接近真值，說明M對EnKF有重要影響,它的大小在很大程度上影響了EnKF的校正精度。這一點可以進一步通過均方誤差和殘差得到佐證。

圖6和圖7都表明水位尺度變換系數(shù)M的增大會提高EnKF的校正性能,圖7還顯示殘差隨變換系數(shù)呈對數(shù)分布,在M從5增加到20的過程中,殘差有快速下降的過程,隨后趨于平緩,這再次證明了上面的結論。M的這一變化特性值得上升到理論層面進行討論,筆者將進行后續(xù)研究。

圖6 不同水位尺度變換系數(shù)下的流量均方誤差

圖7 水位尺度變換系數(shù)與殘差關系曲線

5 結 語

本文建立了基于集合卡爾曼濾波的一維非恒定流多變量交替校正方法,通過分別構建水位和流量濾波器,實現(xiàn)了對水位和流量的單獨校正,再將校正結果統(tǒng)一作為下一個預報集合的初始場,增加擾動,循環(huán)計算。水位校正和流量校正交替進行,不僅充分挖掘了隱藏在觀測數(shù)據(jù)中的有效信息,而且避免了濾波過程中的大矩陣計算,顯著提高了計算效率。在濾波前,還對水位進行了尺度變換,改善了水位濾波效果。

數(shù)值試驗重點考察了觀測誤差和水位尺度變換系數(shù)對模型計算精度的影響,并以均方誤差為指標,檢驗了校正的效果。試驗結果表明,在集合卡爾曼濾波中校正效果主要取決于觀測精度,觀測精度越高,則均方誤差越小,均方誤差與觀測誤差之間大致呈線性關系;水位尺度變換系數(shù)能顯著增強多變量交替校正方法的效果,且變換系數(shù)越大,濾波精度越高。但濾波精度在變換系數(shù)增大到一定程度后提高的幅度越來越小。

[1] 賴錫軍,傅國圣,孫波，等.非恒定水流計算的最優(yōu)控制問題及其變分求解[J].水科學進展,2008,19(4):537-545.(LAI Xijun,FU Guosheng,SUN Bo,et al.Optimal control problems in unsteady flow computation and their variational solutions[J].Advances in Water Science,2008,19(4):537-545.(in Chinese))

[2] 葛守西,程海云,李玉榮，等.水動力學模型卡爾曼濾波實時校正技術[J].水利學報,2005,36(6):687-693.(GE Shouxi,CHENG Haiyun,LI Yurong,et al.Real time updating of hydrodynamic model by using Kalman filter[J].Journal of Hydraulic Engineering,2005,36(6):678-693.(in Chinese))

[3] 周全.洪水預報實時校正方法研究[D].南京:河海大學,2005.

[4] 王船海,吳曉玲,周全，等.卡爾曼濾波校正技術在水動力學模型實時洪水預報中的應用[J].河海大學學報(自然科學版),2008,36(3):300-305.(WANG Chuanhai,WU Xiaoling,ZHOU Quan,et al.Application of Kalman filter technique in real-time flood forecasting[J].Journal of Hohai University(Natural Sciences),2008,36(3):300-305.(in Chinese))

[5] 李大勇,董增川,劉凌,等.一維非恒定流模型與卡爾曼濾波耦合的實時交替校正方法研究[J].水利學報,2007,38(3):330-336.(LI Dayong,DONG Zengchuan,LIU Ling,et al.Real-time alternative updating model coupling 1-D unsteadyflow computation with Kalman filter[J].Journal of Hydraulic Engineering,2007,38(3):330-336.(in Chinese))

[6] 吳曉玲,王船海,向小華,等.實時校正中系統(tǒng)噪聲均值的空間分布[J].河海大學學報(自然科學版),2008,36(4):448-451.(WU Xiaoling,WANG Chuanhai,XIANG Xiaohua,et al.Spatial distribution of mean value of system noise in real-time updating[J].Journal of Hohai University(Natural Sciences),2008,36(4):448-451.(in Chinese))

[7] MADSEN H,SKOTNER C.Adaptive state updating in real-time river flow forecasting-a combined filtering and error forecasting procedure [J].Journal of Hydrology,2005,308(1/2/3/4):302-312.

[8] EVENSEN G.Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods to forecast error statistics[J].Journal of Geophysical Research,1994,99(5):10143-10162.

[9] HOUTEKAMER P,MITCHELL H.Data assimilation using an ensemble Kalman filter technique[J].Monthly Weather Review,1998,126(3):796-811.

[10] HAUGEN V,EVENSEN G.Assimilation of SLA and SST data into an OGCM for the Indian Ocean[J].Ocean Dynamics,2002,52(3):133-151.

[11] MClAUGHLIN D.An integrated approach to hydrologic data assimilation: interpolation,smoothing,and filtering[J].Advances in Water Resources,2002,25(8):1275-1286.

[12] MORADKHANI H,SOROOSHIAN S,GUPTA H,et al.Dual stateparameter estimation of hydrological models using ensemble Kalman filter[J].Advances in Water Resources,2005,28:135-147.

[13] 賴錫軍.水動力學模型與集合卡爾曼濾波耦合的實時校正多變量分析方法[J].水科學進展,2009,20(2):241-248.(LAI Xijun.Real-updating multivariate analysis for unsteady flows with ensemble Kalman filter [J].Advances in Water Science,2009,20(2): 241-248.(in Chinese))

[14] 陳一帆,程偉平,錢鏡林,等.基于集合卡爾曼濾波的河網(wǎng)水情數(shù)據(jù)同化[J].四川大學學報(工程科學版),2014,46(4):26-32.(CHEN Yifan,CHENG Weiping,QIAN Jinglin,et al.Data assimilation of river networks using ensemble kalman filtering algorithm[J].Journal of Sichuan University(Engineering Science Edition),2014,46(4):26-32.(in Chinese))

[15] 賴錫軍,姜加虎,黃群,等.漫灘河道洪水演算的水動力學模型[J].水利水運工程學報,2005(4):29-35.(LAI Xijun,JIANG Jiahu,HUANG Qun,et al.Hydrodynamicmodel for flood routing in a flood plain river[J].Hydro-Science and Engineering,2005(4):29-35.(in Chinese))

[16] 汪德爟.計算水力學理論與應用[M].南京: 河海大學出版社,1989.

A real-time alternating updating method based on ensemble Kalman filter

GU Luhua1,2, LAI Xijun1

(1.NanjingInstituteofGeographyandLimnology,ChineseAcademyofSciences,Nanjing210008,China; 2.CollegeofResourcesandEnvironment,UniversityofChineseAcademyofSciences,Beijing101407,China)

To reduce the uncertainty in calculation of unsteady flows, a multivariate alternate updating method is proposed based on the ensemble Kalman filter. This method updates water stage and discharge data alternately to calibrate unsteady flow, using the observed information without the large matrix calculating; meanwhile, scaling transformation is used in order to improve the water level filter precision. Numerical experiments emphatically investigate the effects of measurement accuracy and water level transformation coefficient on forecast precision of the method. The results show that the forecast error increases as the measurement accuracy decreases; the water level transformation coefficient can obviously improve the effect of the multivariate alternate updating method, the larger the water level transformation coefficient is, the higher the forecast precision will be; the multivariate alternate updating method has good calibrating performance and can improve forecast accuracy of unsteady flows in open channel.

hydrodynamic model; ensemble Kalman filter; unsteady flow; real-time updating; flood forecasting

水體污染控制與治理科技重大專項(2012ZX07101-011);國家自然科學基金(51379059, 51279047)

顧爐華(1993—), 男,碩士研究生,主要從事環(huán)境水力學研究。E-mail:lhgu21@163.com

賴錫軍(1977—), 男,副研究員,主要從事環(huán)境水力學研究。E-mail:xjlai@niglas.ac.cn

10.3880/j.issn.1006-7647.2017.02.013

TV131.2

：A

：1006-7647(2017)02-0073-05

2016-01-28 編輯：駱超)