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

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

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

0 引 言

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

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

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

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

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

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

(1)

(2)

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

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

C0ΔZ1+D0ΔQ1=G0

(3)

CNΔZN+DNΔQN=GN

(4)

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

1.2 觀測變量為狀態向量的Kalman濾波自適應模型

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

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

(5)

則

(6)

結合式(3)～式(6)，加上噪聲得到狀態方程

(7)

觀測方程為

yt=xt+vt

(8)

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

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

當狀態變量維數n等于量測變量維數m時，可以估計出Rk和Qk值；當狀態變量維數n大于量測變量維數m，雖然能夠估計出Rk的統計值，但難以估計出Qk統計值。考慮到實際預報中，預報斷面數目往往多于實測斷面數，理論上難遇估計出Qk值。因此，本研究假設Qk已知，建立以觀測變量為狀態向量的Kalman濾波自適應遞推模型。

2 基于Kalman濾波的實時校正模型應用分析

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

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

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

表1 Kalman濾波模型實時校正結果比較

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

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

圖3 正陽關水位Kalman濾波模型外推6 h預報過程

3 結論與討論

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