顧爐華,賴(lài)錫軍

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

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

顧爐華1,2,賴(lài)錫軍1

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

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

水動(dòng)力學(xué)模型;集合卡爾曼濾波;非恒定流;實(shí)時(shí)校正技術(shù);洪水預(yù)報(bào)

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

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

1 一維非恒定流基本方程

一維非恒定流的基本方程,又稱(chēng)為圣維南方程,可寫(xiě)成:

(1)

(2)

式中:Z為水位,m;Q為流量,m3/s;K為流量模數(shù)，m3/s;q為單位河長(zhǎng)旁側(cè)入流,m2/s;A為過(guò)水?dāng)嗝婷娣e,m2;α為動(dòng)量校正系數(shù);g為重力加速度,m/s2;x為沿水流方向距離,m;t為時(shí)間,s;B為水面寬度,m。

采用Preissmann四點(diǎn)加權(quán)隱格式離散上述兩式得到:

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

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

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

2 集合卡爾曼濾波

2.1 濾波原理

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

2.2 計(jì)算過(guò)程

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

Xk=AXk-1+BUk+Wk

(5)

Zk=Hk+Vk

(6)

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

狀態(tài)變量預(yù)測(cè):

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時(shí)刻觀測(cè)得到的分析解；DK為k時(shí)刻觀測(cè)值。

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

(11)

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

3 多變量交替校正方法

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

圖1 多變量交替校正計(jì)算流程

4 數(shù)值試驗(yàn)

4.1 案例選擇

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

4.2 試驗(yàn)設(shè)計(jì)

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

表1 不同變換系數(shù)和誤差標(biāo)準(zhǔn)差的試驗(yàn)方案

初始集合是在恒定流條件基礎(chǔ)上疊加相應(yīng)的擾動(dòng)得到,表達(dá)式如下:

Uk=Uk(1+Nkσu)

(12)

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

4.3 試驗(yàn)結(jié)果及分析

4.3.1 觀測(cè)誤差對(duì)EnKF的影響

選取距河道起點(diǎn)11 km和47 km處的斷面作為研究對(duì)象,計(jì)算不同條件下的分析值、預(yù)報(bào)值和真值,得到流量隨時(shí)間的變化曲線(xiàn)如圖2所示。

圖2 流量過(guò)程比較

圖2(a)顯示隨著觀測(cè)誤差的增加,分析值與真值差距變大,分析值的精度越來(lái)越低,這是因?yàn)樵谏嫌芜吔鐥l件不準(zhǔn)確時(shí),模型的精度主要取決于觀測(cè)的精度,觀測(cè)精度越高,模型精度也越高。圖2(b)中各種精度下的分析值與真值已經(jīng)相差無(wú)幾,體現(xiàn)了很好的校正效果。對(duì)比兩圖能夠發(fā)現(xiàn)距河道起點(diǎn)11 km的斷面的校正效果明顯差于距河道起點(diǎn)47 km的斷面,這是因?yàn)槟Ｐ椭饕牟淮_定性來(lái)自上游來(lái)流,接近上游的斷面在濾波時(shí)受上游邊界條件影響較大;而下游邊界是準(zhǔn)確的,對(duì)濾波中的下游斷面施加影響從而使得校正精度很高。

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

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

圖3 不同觀測(cè)誤差下的均方誤差

圖4 觀測(cè)誤差與殘差關(guān)系曲線(xiàn)

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

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

圖5 中間斷面流量過(guò)程

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

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

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

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

5 結(jié) 語(yǔ)

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

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

[1] 賴(lài)錫軍,傅國(guó)圣,孫波，等.非恒定水流計(jì)算的最優(yōu)控制問(wèn)題及其變分求解[J].水科學(xué)進(jìn)展,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] 葛守西,程海云,李玉榮，等.水動(dòng)力學(xué)模型卡爾曼濾波實(shí)時(shí)校正技術(shù)[J].水利學(xué)報(bào),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] 周全.洪水預(yù)報(bào)實(shí)時(shí)校正方法研究[D].南京:河海大學(xué),2005.

[4] 王船海,吳曉玲,周全，等.卡爾曼濾波校正技術(shù)在水動(dòng)力學(xué)模型實(shí)時(shí)洪水預(yù)報(bào)中的應(yīng)用[J].河海大學(xué)學(xué)報(bào)(自然科學(xué)版),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] 李大勇,董增川,劉凌,等.一維非恒定流模型與卡爾曼濾波耦合的實(shí)時(shí)交替校正方法研究[J].水利學(xué)報(bào),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] 吳曉玲,王船海,向小華,等.實(shí)時(shí)校正中系統(tǒng)噪聲均值的空間分布[J].河海大學(xué)學(xué)報(bào)(自然科學(xué)版),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] 賴(lài)錫軍.水動(dòng)力學(xué)模型與集合卡爾曼濾波耦合的實(shí)時(shí)校正多變量分析方法[J].水科學(xué)進(jìn)展,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] 陳一帆,程偉平,錢(qián)鏡林,等.基于集合卡爾曼濾波的河網(wǎng)水情數(shù)據(jù)同化[J].四川大學(xué)學(xué)報(bào)(工程科學(xué)版),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] 賴(lài)錫軍,姜加虎,黃群,等.漫灘河道洪水演算的水動(dòng)力學(xué)模型[J].水利水運(yùn)工程學(xué)報(bào),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] 汪德?tīng)?計(jì)算水力學(xué)理論與應(yīng)用[M].南京: 河海大學(xué)出版社,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

水體污染控制與治理科技重大專(zhuān)項(xiàng)(2012ZX07101-011);國(guó)家自然科學(xué)基金(51379059, 51279047)

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

賴(lài)錫軍(1977—), 男,副研究員,主要從事環(huán)境水力學(xué)研究。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 編輯：駱超)