一種基于軟卡方檢測的自適應Kalman濾波方法

嚴恭敏, 劉 璠,, 李梓陽, 周 琪

(1. 西北工業大學自動化學院, 西安 710072;2. 中國航空工業集團公司西安飛行自動控制研究所, 西安 710076)

0 引言

卡爾曼濾波(Kalman filter, KF)是一種利用隨機線性系統狀態空間模型,通過系統輸出的觀測數據,對系統狀態進行最優估計的方法,在導航制導與控制、通信以及圖像處理等眾多領域有著廣泛的應用。

理論上,只有在隨機線性系統的結構參數和噪聲統計特性參數都準確已知的條件下,KF才能獲得狀態的最優估計。如果實際量測的統計特性與系統建模參數不匹配,會導致KF狀態估計精度下降,嚴重時還可能會引起濾波發散。卡方檢測(Chi-square test, C2T)方法是一種傳統的量測故障檢測方法[1],在狀態空間建模準確的前提下,通過KF的新息及其均方差陣構造卡方統計量,可檢測出量測是否存在異常,如果量測正常則進行量測更新,反之,如果量測異常則放棄量測更新。然而,實際應用中的系統噪聲參數往往難以準確確定,從而使得卡方檢測統計量的閾值設置成為一大難題。如果卡方檢測閾值設置太大,則可能將異常量測引入濾波器,降低濾波估計精度;如果閾值設置太小,則可能排除了一些正常的量測,使得量測利用率下降,同樣也會降低狀態估計精度。造成這一問題的根源在于傳統卡方檢測的二值化量測異常判斷原則,量測要么被判斷為正常、要么被判斷為異常,再沒有其他任何中間情形。

借鑒Sage-Husa自適應濾波方法[2-5](Sage-Husa adaptive KF, SHAKF)或抗差自適應濾波方法[6-7],至少可以更精細地將量測狀況劃分為“正常/信任”“可疑/部分信任”和“異常/丟棄”這三種類別,針對不同的量測類別修改量測噪聲方差陣的大小或賦予不同的量測利用權重,或者直接根據KF新息大小計算權重,實現權重大小的無縫連續過渡。但與上述兩種自適應濾波的思路不完全相同,本文通過分析傳統卡方檢測的新息權重二值化利用特點,給出了新息權重的連續化利用方法,將權重的二值化取值改進為連續化取值,并將其推廣,從而獲得了多維量測的多分量卡方檢測方法。文中將所提新方法稱為軟卡方檢測KF方法(soft Chi-square test KF,SC2TKF),傳統方法可被稱為硬卡方檢測KF方法。顯然,新方法更具普遍性,傳統硬卡方檢測可以看作是新軟卡方檢測方法的一個特例。最后,利用慣導/衛導組合進行了軟卡方檢測方法的KF仿真驗證,結果顯示新方法具有比Sage-Husa自適應濾波更好的濾波精度和穩定性。

1 傳統的量測故障卡方檢測原理

隨機線性系統的狀態空間模型記為

(1)

式中,Xk是n×1維的狀態矢量,Zk是m×1維的量測矢量;Φk/k-1,Γk-1,Hk是已知的系統結構參數,分別稱為n×n維的狀態一步轉移矩陣、n×l維的系統噪聲分配矩陣、m×n維的量測矩陣;Wk-1是l×1維的系統噪聲矢量,Vk是m×1維的量測噪聲矢量,兩者都是零均值的高斯白噪聲矢量序列,且它們之間互不相關,即滿足如下Kalman濾波關于噪聲的基本假設條件

(2)

其中,Qk是半正定的,Rk是正定的,δkj為克羅內克函數。

針對系統式(1)的Kalman濾波公式為

(3)

由量測新息及其均方差陣可構造統計量[1]

(4)

其中,λk服從自由度為m的卡方分布,即λk～χ2(m)。在量測正常情況下,統計量λk的數值應當比較小;而如果量測出現異常,λk將變得較大,量測正常與否一般以某選定的閾值TDm作為判斷門限,即

(5)

這便是Kalman濾波的量測故障卡方檢測原理。在濾波過程中可以計算統計量λk,根據其大小實時監測量測是否異常,進而決定是否進行量測更新,式(3)中的量測及其均方差陣更新方程可改寫為

(6)

顯然,當χk=1即λk≤TDm時,進行正常的Kalman濾波量測更新;而當χk=0即λk>TDm時,則放棄量測更新,達到隔離異常量測的效果。

2 新的軟卡方檢測方法

傳統卡方檢測方法有效的前提條件是濾波系統模型準確無誤。在慣導/衛導組合導航等實際應用中,一般Kalman濾波的系統結構參數建模會相對比較準確且穩定,而噪聲參數會由于運動狀態、建模不完善、系統老化或運行環境等原因而變化,難以完全精確建模,這樣就會使得卡方檢測閾值TDm難以嚴格地按理論方法準確確定。閾值設置過大會造成故障檢測概率降低,存在將較多的異常值引入濾波量測的風險,從而造成濾波誤差增大;閾值設置過小又會造成故障虛警概率增大,經常性的虛警降低了量測利用率,量測修正作用減小也會降低濾波估計精度。當然,對于高精度慣導系統而言,隨機性棄用衛導定位量測50%以上,比如衛導量測即使從1 Hz降為0.1 Hz,對組合導航性能影響也不大;但對于低精度慣導系統,衛導利用從1 Hz變為0.5 Hz對組合精度的影響可能就比較顯著了。

傳統卡方檢測方法的結果是二值化的,非0即1,再無任何中間狀態。該方法在雷達目標探測等領域中主要以有/無為指標,其應用是非常合理的,但是,針對組合導航等場合,在量測信息正常(信任)與異常(丟棄)之間還可能存在大量的中間狀態(可疑),僅簡單地使用卡方檢測二值化結果就不太合適了,這也有悖于目前常用的量測自適應濾波和抗差濾波等方法的理念。比如,在Sage-Husa量測自適應Kalman濾波中,通過量測新息的大小自動調整量測噪聲方差陣的大小,相當于是將所有量測信息都進行了綜合利用,與正常量測相比,量測新息中等大小(可疑)時降低量測權重,新息明顯異常時量測權重很小,若權重趨于0則等效于棄用。論文根據自適應濾波的新息利用特點,將式(6)改進為

(7)

(8)

值得特別指出的是,如果式(1)中量測是多維的,傳統卡方檢測將多維量測視為整體,只要有任何一個分量出現異常,卡方檢測方法也會同時降低其他正常量測分量的新息利用率,這種處理方式不是很合理,為了避免該缺陷,需對各量測分量逐一做卡方檢測。如果在Kalman濾波模型中,各量測分量之間是不相關的,即量測噪聲均方差陣Rk為對角線矩陣,則采用序貫濾波處理后,各序貫量測更新的卡方檢測就是相當于對單個量測分量的逐一處理。如果Rk為非對角線矩陣,不妨記rk(i)為新息rk的第i(i=1,2,…,m)個分量,Ak(ii)為新息均方差陣Ak的第i行i列對角線元素,且記統計量

(9)

不難理解,各統計量λk(i)均服從自由度為1的卡方分布,即有λk(i)～χ2(1)。

參考式(7),推廣和建立由λk(i)構造的Kalman濾波量測更新方法,如下

(10)

(11)

(12)

在式(10)中,修改了濾波增益Kk的計算方式,其目的是保證之后均方差陣Pk計算的對稱性,顯然,式(7)是式(10)中ρk各分量都相等時的特殊情形。理論上,ρk可取對稱正定陣,但為了簡便,一般取為如式(11)所示的對角陣。在式(12)中,參數TD1是自由度為1的卡方統計量λk(i)的故障檢測設置閾值,當取顯著性水平為0.05時有TD1≈3.8。

3 仿真驗證

3.1 慣導/衛導組合導航模型

以慣導/衛導松組合導航為例進行仿真實驗驗證,系統狀態15維、量測3維,具體建模如下

(13)

其中

3.2 仿真條件

采用PSINS工具箱軟件模擬一段車載行駛軌跡,軌跡包含靜止、加減速、勻速、轉彎、爬坡等階段,總仿真時間約1 000 s,行駛里程約7.5 km。行車速度如圖1所示,軌跡相對位置變化的二維平面圖如圖2所示,圖中標識“★”為載車行駛起始點。

圖1 載車行駛速度Fig.1 The vehicular velocities

圖2 載車行駛平面軌跡Fig.2 The vehicular trajectory

根據軌跡仿真生成慣性傳感器數據和衛導定位數據,其中慣性數據采樣100 Hz、衛導測量1 Hz。再進行慣性導航解算和慣導/衛導組合導航Kalman濾波解算,仿真過程中注入的各種誤差如表1所列。

表1 仿真誤差設置

衛導的位置量測誤差(緯、經、高,3個分量)均設置為零均值的高斯白噪聲,正常時間段內方差為(1 m)2,以下兩個時間段除外:(1)200～400 s之間的噪聲方差設置為(10 m)2;(2)600～800 s之間設置為污染噪聲,以衛導高度測量誤差為例,其表示為

(14)

其中,“w.p.”表示“以…概率”(with probability)之意。式(14)表示噪聲δH在正常情況下以90%的大概率服從方差為(1 m)2的高斯分布,視為正常噪聲;而在異常情況下以10%的小概率服從方差為(100 m)2的高斯分布,視為野值噪聲。

一組高度量測誤差仿真的示例如圖3所示,由圖可見,在200～400 s之間高度誤差總體波動變大,模擬噪聲方差變化;而在600～800 s之間個別誤差幅值跳變很大,模擬小概率的野值跳變。

圖3 衛導高度測量誤差Fig.3 GNSS altitude measurement error

3.3 仿真結果

慣導/衛導組合導航的姿態失準角、速度誤差和位置誤差的結果如圖4所示,從圖中可以看出:在200～400 s之間隨著衛導的量測噪聲增大,組合導航的誤差也相應變大,該段誤差的量級與全程衛導誤差方差都設置為(10 m)2的結果是基本一致的;在600～800 s之間狀態估計基本不受衛導量測噪聲野值的影響,濾波器具有良好的野值隔離效果,與只存在小方差(1 m)2時的量測噪聲一樣,始終保持較高的濾波精度。

(a) 姿態誤差

(b) 速度誤差

(c) 位置誤差圖4 基于軟卡方檢測的估計誤差Fig.4 State estimation error based on SC2TKF

此外,還對本文新方法(SC2TKF)和傳統Sage-Husa量測噪聲自適應濾波(SHAKF)作了對比仿真,進行了50次的蒙特卡洛仿真并統計各導航參數誤差,結果表明SC2TKF更加穩定,其中緯度誤差的均方根(root mean square,RMS)誤差統計如圖5所示,其他誤差情況類似,不再一一列出。值得注意的是,在SHAKF中還有一個漸消因子需要精心設置(文中取b=0.5),若漸消因子設置不合適,對濾波結果有較大的負面影響。綜上,新方法無需設置任何參數,具有使用方便和適應性強等優點。

圖5 不同濾波方法的緯度誤差比較Fig.5 Comparison of latitude error by using SC2TKF vs SHAKF

4 結論

自適應Kalman濾波的方法很多,廣義上說,能根據實時量測自動調整濾波器參數并獲得良好狀態估計效果的方法均可稱為自適應濾波,比如Sage-Husa自適應濾波[2-3]、抗差自適應濾波[4,6-7]、強跟蹤濾波[8]以及近年來一系列文獻中根據變分貝葉斯估計原理給出的RSTKF、MCKF、SSMKF和CERKF等眾多方法[9-13]。在組合導航實際系統中,除基于量測噪聲自適應的經典Sage-Husa自適應濾波比較實用外,其他自適應方法理論推導結果和仿真效果雖然較好,但往往前提假設過于理想或者參數設置繁瑣,難以適用于復雜高維系統,因而少有實際應用的報道。繼Sage-Husa自適應濾波之后,論文提出的軟卡方檢測Kalman濾波方法有望在組合導航領域成為一種更加實用的自適應濾波方法。