抗差Kalman濾波算法研究及其在GPS監(jiān)測數據中的應用

焦雄風 陳 錚 楊興旺 索廣建 張獻州

(1.西南交通大學地球科學與環(huán)境工程學院,成都 611756; 2.西南交通大學高速鐵路運營安全空間信息技術國家地方聯合工程實驗室,成都 611756; 3.上海鐵路北斗測量工程技術有限公司,上海 200071; 4.中國鐵路上海局集團有限公司,上海 200071)

1 概述

Kalman濾波在許多領域應用廣泛,但受制于其動態(tài)系統模型中隨機噪聲設定(必須為白噪聲),對觀測手段和環(huán)境要求較高?？共頚alman濾波的提出可有效解決上述問題。趙長勝針對GPS精密定位問題,提出有色噪聲作用下的抗差Kalman濾波,并對其公式進行詳細推導[1];劉華夏探討一種抗差自適應Kalman濾波方法,并將其應用于高速鐵路變形監(jiān)測數據的分析中[2];張海平針對BDS-3星座的中長基線實時動態(tài)定位,等提出一種基于組合觀測值的實時動態(tài)(RTK)卡爾曼濾波定位算法,并獲得厘米級的定位精度[3-4];賀晗等將抗差Kalman濾波應用于塌陷區(qū)監(jiān)測數據的處理中[5];余航等推導了針對動態(tài)EIV模型的總體卡爾曼濾波方法及其近似精度評定公式[6]。

然而,基于不同準則下的抗差Kalman濾波性能差異較大。以下對抗差Kalman濾波、經典Kalman濾波以及小波閾值去噪方法應用于實際GPS單歷元解算進行研究,以期獲得最優(yōu)方法。

2 抗差Kalman濾波

2.1 經典Kalman濾波

Kalman濾波是一種能夠對包含噪聲的觀測數據進行降噪處理,從而得到觀測目標狀態(tài)最優(yōu)估計的算法[7],其離散線性系統下的Kalman濾波基本函數模型方程組由狀態(tài)方程和觀測方程構成,即

式中,xk為(n×1)初始狀態(tài)參數矩陣;zk為(m×1)測量參數矩陣;wk為(n×1)動態(tài)噪聲;vk為(m×1)測量噪聲;Φk/k-1為(n×n)狀態(tài)轉移矩陣;Hk為(m×n)觀測矩陣;Bk/k-1為(n×r)控制參數的增益矩陣;uk-1為(r×1)控制參數矩陣,下標k表示第k時刻。

上述系統的隨機模型設定為

式中,Qk為動態(tài)噪聲wk的非負定方差矩陣;Rk為系統觀測噪聲vk的正定方差矩陣;δkj為克羅內克函數,即

上述^X可按如下幾步公式遞推演算得到。

狀態(tài)值及估計協方差矩陣一步預測為

觀測噪聲協方差及最優(yōu)卡爾曼增益矩陣一步更新為

更新的狀態(tài)估及協方差估計為

通過上述公式遞推流程,若給定初值狀態(tài)參數以及在k時刻下的觀測值Zk,就能得到在k時刻下目標狀態(tài)的一步預報值^xk/k-1和估值xk(k=1,2,3,…,N)。

2.2 Huber法抗差Kalman濾波

經典Kalman濾波迭代方式主要由預測加更新的線性回歸方式構成,其主要結合觀測量和增益矩陣K進行估計[8],可將其構造為

令

由此,卡爾曼濾波迭代式可等價于式(10)的線形回歸方程式,而對于求解上述線形回歸方程式,可以利用穩(wěn)健估計中的“M估計”來解決[9],即

式中,εk,i為εk的第i個元素。

通過上述穩(wěn)健估計方法,可得Kalman更新方程為

2.3 最大相關熵法抗差Kalman濾波

如2.2節(jié)推導類似,由經典Kalman濾波迭代公式可以構造得到式(8),再令

式中,K(p)、K(r)是用平方根法對P、R進行分解得到[10],對式(8)同時乘,可得

由此,Kalman濾波迭代式同樣可等價于式(13)的線形回歸方程式,而對于求解上述線形回歸方程式,可以采用如下方法[11],即

式中,N為已知信息個數,c(i)k為Ck的第i個元素,d(i)k為Dk的第i行元素。根據最大相關熵準則可知xk的最優(yōu)估值為

求取上述最優(yōu)估值等價于求取xk的一個固定解方程,即

求解上述固定解方程,再結合固定解最終形式及經典Kalman濾波的更新方程[12],有

3 精度評價指標

引入均方根誤差RMSE(Root Mean Square Error)、平均絕對誤差MAE(Mean Absolute Error)、信噪比dB(signal-to-noise ratio)、平滑度r(smoothing of signal)來對濾波去噪效果進行評價[13]。其中,RMSE能反映實際信號與濾波信號之間的偏差,其值越小表明其濾波效果精度越好;MAE能反映濾波信號誤差的實際情況,其值越小表明模型擬合程度越高,dB能衡量實際信號和濾波信號的相似度,其值越高表明濾波效果越好;r用來評價濾波信號的整體誤差情況,其平滑度越小表明濾波整體誤差偏移量越小,具體公式為

式中,n均為采樣次數;νk為濾波值與實際值的差值。

4 實例分析

4.1 工程概況

某跨長江公鐵兩用特大橋全長約10km,位于江蘇省境內,其主航道橋段采用雙塔三索面箱桁組合梁斜拉橋方案,主跨長約1km,橋上鐵路列車設計速度為200km/h,在其主跨段雙塔頂端、跨中及東西兩側各設1個GNSS監(jiān)測點。

4.2 濾波效果比較及精度評價

上述監(jiān)測點獲取的GPS單歷元解算數據涉及到監(jiān)測點位置、變化速度以及時間的計算,這些參數主要受多路徑效應、觀測值噪聲及其他誤差源的影響,在上述誤差的影響下,可能會導致解算的數據出現較大偏差[14]。一般情況下,偶然誤差可以通過經典Kalman濾波以及其他一些去噪方法進行去除,但當誤差超出偶然誤差的范圍時,上述方法就不適用了。而抗差Kalman濾波可以基于穩(wěn)健估計準則或通過最大化相關熵,對異常粗差進行有效探測及去除,以保證GPS單歷元解算數據的準確性[15]。

為比較上述抗差Kalman濾波及經典Kalman濾波效果,選取該橋跨中監(jiān)測點某個觀測時段的GPS單歷元解算數據,得到該監(jiān)測點在北(N)、東(E)方向上的實時位置(de、dn),本次采樣個數共3587個,如圖1所示。在第100、250、500、900、1300、1800、2900、3300處,人為添加不同數值的粗差,再采用Huber抗差Kalman濾波(HKF)、最大相關熵法抗差Kalman濾波(MKF)、經典Kalman濾波(CKF)進行分析,為測試三者的去噪能力,設置同樣的初始參數,狀態(tài)向量參數設置為位置和速度,初始位置設為不同方向上第一個時間點采樣的位置,初始速度設為0,初始估計協方差初始動態(tài)噪聲的協方差矩陣Q0=初始觀測噪聲協方差矩陣R0=狀態(tài)轉移矩陣觀測向量H=具體濾波結果見圖2、圖3。

圖1 實例數據

由圖2、圖3可知,在無異常干擾的觀測數據處,三者濾波效果基本上都達到優(yōu)良范圍;而在人為添加的粗差處,最大相關熵法抗差Kalman濾波表現出強大的去噪能力,有效排除了異常噪聲的干擾,其殘差范圍控制在厘米級;Huber法抗差Kalman濾波也能探測到部分粗差,但其殘差范圍僅為分米級;經典Kalman濾波在異常噪聲處的去噪效果最差,基本上無法濾掉粗差,其在異常噪聲處的殘差最大為20m。為進一步對上述抗差Kalman方法的性能進行綜合評價,采用小波閾值分析的方法,對上述添加粗差的GPS單歷元差分解算數據進行去噪處理。小波閾值也是目前在GPS單歷元解算數據去噪方面較為常用的方法[16-17],閾值類型可選擇rigrsure(無偏估計),閾值使用方式可為軟閾值,且閾值處理不隨噪聲水平變化而變化,具體處理結果如圖4、圖5所示。

圖2 de濾波結果及殘差

圖3 dn濾波結果及殘差

圖4 de小波閾值去噪結果及殘差

圖5 dn小波閾值去噪結果及殘差

由圖4、圖5可知,小波閾值分析方法只能消弱偶然誤差的影響,無法對粗差進行有效探測,達不到有效消除粗差的要求。為進一步定量分析上述幾種方法的去噪效果,利用第3節(jié)提到的精度指標來評價四者的去噪能力,相關統計數值如表1所示。

表1 三種濾波去噪能力評價指標

由表1可知,最大相關熵法抗差Kalman濾波、Huber抗差Kalman濾波的精度指標總體上優(yōu)于經典Kalman濾波和小波閾值去噪方法。通過對比,表明最大相關熵法抗差Kalman濾波去噪效果最好,Huber法抗差Kalman濾波次之,經典Kalman濾波和小波閾值去噪效果最差。綜合可知,最大相關熵法抗差Kalman濾波能有效探測粗差并改善估值精度,去噪能力為四者中最優(yōu)。

5 結論

通過對含有粗差的GPS單歷元解算數據進行分析,比較最大相關熵法抗差Kalman濾波、Huber法抗差Kalman濾波以及經典Kalman濾波在異常噪聲下的去噪效果。實驗結果表明,最大相關熵法抗差Kalman濾波具有最優(yōu)降噪能力。其原因為,最大相關熵法抗差Kalman濾波從估計誤差和殘差的相關熵方面考慮,通過最大化相關熵來提高狀態(tài)估計的精度,該方法在異常噪聲的干擾下仍然適用;Huber法抗差Kalman濾波利用最小二乘的原理進行估計,這就導致即使觀測值為粗差時,仍然會對其賦予一定的權重,從而導致濾波精度下降;經典Kalman濾波也是利用加權平均的思想,但迭代方式決定了其對異常噪聲基本上無法有效去除;而小波閾值去噪的核心是根據臨界閾值來判斷噪聲和真實信息,臨界閾值則是根據小波系數的大小來設置,而粗差的出現會導致小波系數出現嚴重的不均勻性,繼而可能產生錯誤的閾值設定,導致去噪效果出現嚴重下滑。