基于Kalman濾波的實時校正模型研究

王莉莉，包紅軍，2，李致家

(1.國家氣象中心，北京 100081；2.中國氣象局-河海大學(xué)水文氣象研究聯(lián)合實驗室，北京 100081；3.河海大學(xué)水文水資源學(xué)院，江蘇 南京 210098)

0 引 言

洪水實時預(yù)報以動態(tài)預(yù)報洪水變化過程為目的，由洪水預(yù)報模型和實時校正方法兩部分組成[1]。洪水預(yù)報模型反映水文系統(tǒng)的內(nèi)在變化規(guī)律，從根本上決定流域水文物理過程的模擬精度[2]；實時校正方法根據(jù)當(dāng)前實時信息，應(yīng)用現(xiàn)代系統(tǒng)理論和方法，對預(yù)報模型參數(shù)或狀態(tài)變量或輸入向量或預(yù)報值進(jìn)行實時訂正，建立對系統(tǒng)模型與預(yù)報的現(xiàn)時校正的回饋機制[3]。洪水預(yù)報實時校正的依據(jù)來自于預(yù)報誤差，主要包括模型結(jié)構(gòu)誤差、模型參數(shù)誤差、模型輸入誤差、初始條件誤差[4-7]。

根據(jù)實時校正所依據(jù)理論區(qū)分，實時預(yù)報方法主要可劃分為[3，8]：①以現(xiàn)代系統(tǒng)理論為基礎(chǔ)的系統(tǒng)類方法，以最小二乘法、Kalman濾波等為最為典型[9]；②以現(xiàn)代時間序列分析為基礎(chǔ)的統(tǒng)計類方法，對預(yù)報模型的參差序列進(jìn)行預(yù)測，與模型本身分離，包括誤差自回歸模型、自回歸滑動平均模型等；③簡易經(jīng)驗校正方法及利用計算機圖形功能的人機交互校正方法。

目前，水動力學(xué)模型在洪水預(yù)報中的研究表明，單純求解水動力學(xué)方程的差分技術(shù)已經(jīng)比較充分，進(jìn)一步提高精度的潛力不大[9]。為了探索從實時校正角度提高洪水預(yù)報精度，本研究以淮河中游河道為例，將水動力學(xué)模型動態(tài)化，建立以觀測變量為狀態(tài)向量的水動力學(xué)模型Kalman濾波實時校正模型，并在試驗流域典型洪水預(yù)報中進(jìn)行驗證。

1 基于Kalman濾波的實時校正模型

1.1 圣維南方程的離散化和線性化

擴展的圣維南方程由連續(xù)方程和動量方程組成

(1)

(2)

可采用Preissmann四點隱式差分格式對連續(xù)方程和動量方程進(jìn)行時間(i)和空間(j)的離散化。

形成的兩個方程，對N個斷面的N-1個河段可寫出2N-2個方程，加上兩個邊界條件方程，總共是2N個代數(shù)方程。其中，含有2N個未知量(N個斷面的ΔZj、ΔQj)；方程組有唯一解。在河道實時洪水預(yù)報中，其上下邊界條件分別為上下游水位或者流量隨時間的變化過程或者水位～流量關(guān)系曲線。上邊界選用流量過程，下邊界選用水位過程。則上下邊界方程分別為

C0ΔZ1+D0ΔQ1=G0

(3)

CNΔZN+DNΔQN=GN

(4)

如果將每個斷面時段初時刻的水位、流量看作輸入，末時刻的水位、流量看作輸出，描述的正是一個長河道洪水變化的線性系統(tǒng)。

1.2 觀測變量為狀態(tài)向量的Kalman濾波自適應(yīng)模型

選擇1～N個斷面的水位、流量(系統(tǒng)輸入輸出變量)為狀態(tài)向量xt并與觀測向量yt同形，即

yt=xt=[Z′1Q′1Z′2Q′2……Z′NQ′N]T

(5)

則

(6)

結(jié)合式(3)～式(6)，加上噪聲得到狀態(tài)方程

(7)

觀測方程為

yt=xt+vt

(8)

這就構(gòu)成了非馬爾可夫形式的狀態(tài)空間方程。其中，狀態(tài)轉(zhuǎn)移陣、觀測矩陣皆為單位陣，ωi、vi分別為模型誤差向量和觀測誤差向量。使用式(7)、(8)即可調(diào)用Kalman濾波器，進(jìn)行濾波的結(jié)果是校正模型和觀測兩種噪聲(誤差)，使xt逐步逼近真值，達(dá)到提高精度的目的。

如果Qk、Rk已知，則可以得出一般的Kalman濾波遞推方程，但是將Kalman濾波遞推模型用于實際的洪水預(yù)報存在幾個問題：①Q(mào)k與Rk已知；②Qk與Rk偏差大，E{ωk}≠0。增益矩陣主要基于Qk與Rk，但是在洪水預(yù)報前得到準(zhǔn)備的Qk與Rk值是很困難的[10]。如何在洪水預(yù)報實時的估計出Qk與Rk值，稱之為自適應(yīng)洪水預(yù)報濾波模型。如果濾波系統(tǒng)是定常的，wk和vk屬于平穩(wěn)隨機過程，在濾波系統(tǒng)過度過程結(jié)束后，Pk與Gk將都為定常。實際河道中水位或者流量都是非平穩(wěn)的隨機過程；因此，系統(tǒng)噪聲wk和量測vk也是屬于非平穩(wěn)的隨機過程。由于洪水預(yù)報實際作業(yè)中的時間步長與樣本長度，到洪水退水期，濾波系統(tǒng)才能達(dá)到平穩(wěn)。

當(dāng)狀態(tài)變量維數(shù)n等于量測變量維數(shù)m時，可以估計出Rk和Qk值；當(dāng)狀態(tài)變量維數(shù)n大于量測變量維數(shù)m，雖然能夠估計出Rk的統(tǒng)計值，但難以估計出Qk統(tǒng)計值?？紤]到實際預(yù)報中，預(yù)報斷面數(shù)目往往多于實測斷面數(shù)，理論上難遇估計出Qk值。因此，本研究假設(shè)Qk已知，建立以觀測變量為狀態(tài)向量的Kalman濾波自適應(yīng)遞推模型。

2 基于Kalman濾波的實時校正模型應(yīng)用分析

本次研究以王家壩水文站至魯臺子水文站的淮河干流為試驗河段(見圖1)。河段右岸有史河、淠河2個支流匯入，并且還有兩個蓄洪區(qū)(城西湖蓄洪區(qū)、城東湖蓄洪區(qū))；左岸有潁河1個主要支流匯入，另外還有洪河分洪道、谷河和潤河匯入和蒙洼蓄洪區(qū)、姜唐聯(lián)湖蓄洪區(qū)、南潤段、邱家湖、潤趙段3個行洪區(qū)。王家壩至魯臺子河段共有3個水文站與4個水位站，在整個流域防汛中，一直處于洪水預(yù)報關(guān)鍵之處。

選取淮河王家壩至魯臺子河段作為試驗河段，基于一維水動力學(xué)模型，以南照集水位、潤河集水位和流量、汪集水位、正陽關(guān)水位作為觀測變量進(jìn)行建立狀態(tài)方程。在半自適應(yīng)濾波模型中，Qk預(yù)先給出，Rk是自動實時估計出來的；Qk是wk的協(xié)方差矩陣。由于wk是非平穩(wěn)的隨機過程，Qk也是時變的，在是實際應(yīng)用中只能預(yù)估個均值。本次研究中取為1 000。Kalman濾波模型中決定校正量的是增益矩陣Gk，洪水預(yù)報中Qk的取值不當(dāng)，濾波器會通過調(diào)整Rk來式增益比保持一定的數(shù)值，促使濾波器保持最優(yōu)。

模型應(yīng)用結(jié)果如圖1～圖3和表1所示。結(jié)果表明：預(yù)見期為6 h的實時外延預(yù)報中，取得了提高預(yù)報精度的較好效果。2007年洪水中各個水位及流量站的預(yù)報結(jié)果均有不同幅度的提高：南照集水位誤差降低了0.23 mm，潤河集、汪集、正陽關(guān)、魯臺子的水位誤差分別降低了0.10、0.09、0.08 mm和0.36 mm，流量的確定系數(shù)均有一定的提高，說明校正模型的建模是合理的，對洪水預(yù)報的精度有一定的提高。

表1 Kalman濾波模型實時校正結(jié)果比較

圖1 潤河集站Kalman濾波模型外推6 h預(yù)報過程

圖2 魯臺子站流量Kalman濾波模型外推6 h預(yù)報過程

圖3 正陽關(guān)水位Kalman濾波模型外推6 h預(yù)報過程

3 結(jié)論與討論

本文建立了以一維水動力學(xué)模型為基礎(chǔ)，以水位作為狀態(tài)量的濾波器校正模型。對于實際河道中要預(yù)報的斷面總是多于有實測資料的斷面。所以從理論上講，無法把Qk的統(tǒng)計特性全部估計出來，但狀態(tài)變量的維數(shù)n大于量測變量的維數(shù)m，Rk的統(tǒng)計特性可以全部實時地估計出來，建立的濾波模型為自適應(yīng)濾波模型。將Kalman濾波應(yīng)用到一維水動力學(xué)模型的實時校正中，使用觀測變量作為狀態(tài)向量，狀態(tài)方程和觀測方程很容易從線性化后的水動力學(xué)差分方程組改建，從而使用自適應(yīng)Kalman濾波實現(xiàn)實時校正，達(dá)到提高精度的目的。在淮河的應(yīng)用表明，它比原水力學(xué)模型預(yù)報精度得到了提高。