焦雄風(fēng) 陳 錚 楊興旺 索廣建 張獻(xiàn)州

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

1 概述

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

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

2 抗差Kalman濾波

2.1 經(jīng)典Kalman濾波

Kalman濾波是一種能夠?qū)Π肼暤挠^測(cè)數(shù)據(jù)進(jìn)行降噪處理,從而得到觀測(cè)目標(biāo)狀態(tài)最優(yōu)估計(jì)的算法[7],其離散線(xiàn)性系統(tǒng)下的Kalman濾波基本函數(shù)模型方程組由狀態(tài)方程和觀測(cè)方程構(gòu)成,即

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

上述系統(tǒng)的隨機(jī)模型設(shè)定為

式中,Qk為動(dòng)態(tài)噪聲wk的非負(fù)定方差矩陣;Rk為系統(tǒng)觀測(cè)噪聲vk的正定方差矩陣;δkj為克羅內(nèi)克函數(shù),即

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

狀態(tài)值及估計(jì)協(xié)方差矩陣一步預(yù)測(cè)為

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

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

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

2.2 Huber法抗差Kalman濾波

經(jīng)典Kalman濾波迭代方式主要由預(yù)測(cè)加更新的線(xiàn)性回歸方式構(gòu)成,其主要結(jié)合觀測(cè)量和增益矩陣K進(jìn)行估計(jì)[8],可將其構(gòu)造為

令

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

式中,εk,i為εk的第i個(gè)元素。

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

2.3 最大相關(guān)熵法抗差Kalman濾波

如2.2節(jié)推導(dǎo)類(lèi)似,由經(jīng)典Kalman濾波迭代公式可以構(gòu)造得到式(8),再令

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

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

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

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

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

3 精度評(píng)價(jià)指標(biāo)

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

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

4 實(shí)例分析

4.1 工程概況

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

4.2 濾波效果比較及精度評(píng)價(jià)

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

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

圖1 實(shí)例數(shù)據(jù)

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

圖2 de濾波結(jié)果及殘差

圖3 dn濾波結(jié)果及殘差

圖4 de小波閾值去噪結(jié)果及殘差

圖5 dn小波閾值去噪結(jié)果及殘差

由圖4、圖5可知,小波閾值分析方法只能消弱偶然誤差的影響,無(wú)法對(duì)粗差進(jìn)行有效探測(cè),達(dá)不到有效消除粗差的要求。為進(jìn)一步定量分析上述幾種方法的去噪效果,利用第3節(jié)提到的精度指標(biāo)來(lái)評(píng)價(jià)四者的去噪能力,相關(guān)統(tǒng)計(jì)數(shù)值如表1所示。

表1 三種濾波去噪能力評(píng)價(jià)指標(biāo)

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

5 結(jié)論

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