非高斯環(huán)境下船變測量最大熵卡爾曼濾波方法

龍子旋, 周 琪, 彭俠夫, 張霄力,*

(1.廈門大學(xué)航空航天學(xué)院, 福建 廈門 361101;2.中國航空工業(yè)集團(tuán)公司西安飛行自動(dòng)控制研究所, 陜西 西安 710076)

0 引 言

慣性量匹配方法因其精度高、實(shí)時(shí)性好而在船體變形測量中得到廣泛應(yīng)用,其本質(zhì)就是建立船變測量過程的卡爾曼濾波方程。目前主流的慣性匹配算法包括角速度匹配法[1-3]、姿態(tài)陣匹配法[4-5]、姿態(tài)角匹配法[6]和姿態(tài)四元數(shù)匹配法[7-8]。

文獻(xiàn)[9]指出,當(dāng)噪聲建模失準(zhǔn)或測量系統(tǒng)受到未知噪聲干擾時(shí),傳統(tǒng)卡爾曼濾波或無法準(zhǔn)確跟蹤角形變。在船體變形測量中,為簡化建模過程,往往將量測噪聲理想化為模型準(zhǔn)確已知的高斯白噪聲,對(duì)模型失準(zhǔn)系統(tǒng)受到有色噪聲干擾下的濾波研究較少。文獻(xiàn)[10]將系統(tǒng)噪聲建模為白噪聲驅(qū)動(dòng)的有色噪聲,然后修改了卡爾曼濾波器的迭代方程來完成非高斯濾波。文獻(xiàn)[11]提出了一種強(qiáng)跟蹤的最大互相關(guān)熵(簡稱為最大熵)卡爾曼濾波,即將自適應(yīng)因子引入傳統(tǒng)的最大熵濾波器改進(jìn)姿態(tài)陣匹配法,減小模型失配誤差。

針對(duì)噪聲統(tǒng)計(jì)特性建模失準(zhǔn)問題,學(xué)者們提出了一種自適應(yīng)卡爾曼濾波,其核心思想是,利用量測數(shù)據(jù)自動(dòng)進(jìn)行相關(guān)迭代矩陣的在線估計(jì)和修正。文獻(xiàn)[12-13]介紹了一種Sage-Husa自適應(yīng)卡爾曼濾波,該濾波器引入單一的遺忘因子來減小歷史數(shù)據(jù)的影響,直接估計(jì)噪聲的一階和二階統(tǒng)計(jì)特性。此后,單自適應(yīng)漸消因子的概念被提出,其本質(zhì)是強(qiáng)迫新息協(xié)方差的理論值與估計(jì)值相等以穩(wěn)定濾波過程[14]。文獻(xiàn)[15-18]引入多通道自適應(yīng)漸消因子對(duì)狀態(tài)一步預(yù)測誤差矩陣進(jìn)行實(shí)時(shí)更新調(diào)整以實(shí)現(xiàn)對(duì)多通道時(shí)變?cè)肼暤淖赃m應(yīng)。此外,文獻(xiàn)[19-20]還引入了χ2檢驗(yàn)來提高自適應(yīng)漸消因子引入時(shí)機(jī),抑制可能出現(xiàn)的狀態(tài)跳變,成功應(yīng)用于捷聯(lián)慣性導(dǎo)航系統(tǒng)(strapdown inertial navigation system, SINS)初始對(duì)準(zhǔn)問題,并推廣到非線性系統(tǒng)中[19]。

針對(duì)非高斯白噪聲的濾波問題,在信息論(information theoretical learning, ITL)中,有關(guān)學(xué)者提出了另一種最優(yōu)理論,即最大熵準(zhǔn)則(maximum correntropy criterion, MCC)[21-22]。其核心思想是,在卡爾曼濾波迭代過程中引入MCC構(gòu)造核函數(shù),以此作為卡爾曼濾波器的目標(biāo)函數(shù),重新求解濾波器的增益矩陣,充分利用量測值的高階信息,這便是用來處理非高斯噪聲的最大相關(guān)熵卡爾曼濾波器。文獻(xiàn)[23]在進(jìn)行增益矩陣的求解時(shí),近似認(rèn)為狀態(tài)值與一步預(yù)測值相等,得到基于MCC的卡爾曼濾波迭代(MCC Kalman filter, MCKF)算法,由于其近似處理較多,故而其算法簡單但濾波器精度較差。文獻(xiàn)[24]通過改進(jìn)近似處理過程并融合加權(quán)最小二乘(weighted least squares, WLS)理論修改高斯核函數(shù)得到性能更優(yōu)的非高斯濾波器,其濾波收斂性更好、精度更高、穩(wěn)定性更好。此外,基于擴(kuò)展卡爾曼濾波器(extend Kalman filter, EKF)和無跡卡爾曼濾波器(unscented Kalman filter, UKF)的最大熵濾波器也相繼被提出以解決工程遇到的實(shí)際問題[25-26],將MCC推廣到非線性系統(tǒng),進(jìn)一步說明了MCC在非線性環(huán)境中同樣適用。

受相關(guān)論文啟發(fā),本文介紹一種基于權(quán)重矩陣的MCCKF算法解決船體變形角估計(jì)中遇到的實(shí)際工程問題。

首先，介紹了最大熵概念,推導(dǎo)了傳統(tǒng)MCKF和改進(jìn)MCCKF迭代公式。然后,引入兩種主流自適應(yīng)濾波器:多重自適應(yīng)漸消因子濾波器(簡稱為MAKF)和基于χ2檢驗(yàn)的多重自適應(yīng)漸消因子濾波器(簡稱為TAKF),并與MCKF一同作為對(duì)比算法，基于船體變形角估計(jì)應(yīng)用場景,在兩種典型非高斯噪聲和兩種船體靜態(tài)變形建模情形下,針對(duì)4種濾波器開展多組對(duì)比實(shí)驗(yàn)。此外,還探討了MAKF、TAKF實(shí)際使用受限問題。最后，通過實(shí)驗(yàn)結(jié)果分析,驗(yàn)證了本文介紹的濾波算法在收斂性、變形估計(jì)精度、實(shí)用性等綜合性能上相對(duì)于主流自適應(yīng)和最大熵濾波器具備一定的優(yōu)越性。

1 MCKF

1.1 傳統(tǒng)MCKF

MCC準(zhǔn)則可充分利用信息的高階統(tǒng)計(jì)特性來衡量兩個(gè)隨機(jī)變量之間的相關(guān)性[27-28],在非高斯環(huán)境中的濾波有廣泛的應(yīng)用[29-30]。最大相關(guān)熵濾波器源于對(duì)經(jīng)典卡爾曼濾波迭代公式中的代價(jià)函數(shù)(又稱目標(biāo)函數(shù))進(jìn)行替換繼而推導(dǎo)得到的。首先,采用如下形式的高斯核函數(shù)替代目標(biāo)函數(shù):

(1)

式(1)表明新的高斯核函數(shù)由兩個(gè)核函數(shù)組成,其中,

(2)

類似于標(biāo)準(zhǔn)卡爾曼濾波器的最小均方誤差(minimum square error, MSE)理論求解最優(yōu)解的方式,求解式(2)的最大熵,令

(3)

得

(4)

即

(5)

為求狀態(tài)變量更新公式,通常在較短濾波更新周期內(nèi)認(rèn)為

(6)

則,式(5)中可有以下相關(guān)近似項(xiàng):

(7)

式中:I為n維的單位矩陣。

聯(lián)立式(5)與式(7)得到狀態(tài)量更新公式:

(8)

式中：增益矩陣可表示為

(9)

(10)

(11)

(12)

最后求得改進(jìn)后的MCKF狀態(tài)量更新公式:

(13)

式(13)還可以得到狀態(tài)預(yù)測的增益矩陣:

(14)

用式(13)和式(14)替換掉經(jīng)典卡爾曼濾波狀態(tài)更新公式便可得到MCKF的迭代算法。算法流程如下所示。

算法1 MCKF算法濾波過程迭代公式初始條件x0=x0,P0=P0一步預(yù)測^xk|k-1=Φk|k-1^xk-1Pk|k-1=Φk|k-1Pk-1ΦTk|k-1+Γk-1Qk-1ΓTk-1濾波增益 Kk=G1(zk-Hk^xk|k-1)HTkI+G1(zk-Hk^xk|k-1)HTkHk濾波更新^xk=^xk|k-1+ KkεkPk=(I- KkHk)Pk|k-1

1.2 改進(jìn)型MCCKF

(15)

類似式(3)方式求其最優(yōu)解,得

(16)

將式(16)改寫為

(17)

其中,

(18)

類似式(11)的處理,得

(19)

最后求得MCCKF狀態(tài)量更新公式:

(20)

類似式(14),得到濾波過程的增益矩陣為

(21)

同樣, MCCKF迭代算法流程如下所示。

算法2 MCCKF算法濾波過程迭代公式初始條件x0=x0,P0=P0一步預(yù)測^xk|k-1=Φk|k-1^xk-1Pk|k-1=Φk|k-1Pk-1ΦTk|k-1+Γk-1Qk-1ΓTk-1濾波增益^Kk=GkHTkR-1kP-1k|k-1+GkHTkR-1kHk濾波更新^xk=^xk|k-1+^KkεkPk=(I-^KkHk)Pk|k-1

2 實(shí)驗(yàn)與分析

目前角速度匹配法在慣性匹配方式體系中應(yīng)用較廣、技術(shù)成熟,因此以下變形濾波驗(yàn)證實(shí)驗(yàn)亦基于此。另外,因試驗(yàn)無法提供基準(zhǔn)及真實(shí)數(shù)據(jù),以下驗(yàn)證實(shí)驗(yàn)均以計(jì)算機(jī)仿真平臺(tái)為主。

2.1 角速度匹配法和實(shí)驗(yàn)參數(shù)設(shè)置

船體變形角包括靜態(tài)變形角和動(dòng)態(tài)變形角。靜態(tài)變形視為“準(zhǔn)靜態(tài)”模型,用隨機(jī)游走模型近似代替其長周期的緩變過程:

(22)

動(dòng)態(tài)變形θi視為白噪聲信號(hào)通過二階濾波器得到,其模型如下:

(23)

式中：i=x,y,z。光纖陀螺器件漂移模型分為兩類,即常值漂移εc和偏置穩(wěn)定性漂移εr,前者視為常量,后者視作一階Markov過程:

(24)

式中：j=1,2。式(23)和式(24)相關(guān)參數(shù)含義如表1所示。

表1 船體變形模型關(guān)鍵參數(shù)

綜合以上變形模型、陀螺漂移模型,將角速度匹配方法的狀態(tài)方程組表示如下:

(25)

另外,量測方程為

(26)

(27)

其中,各個(gè)矩陣的含義由于篇幅有限這里不作展開。實(shí)驗(yàn)給定的船變參數(shù)如表1所示。對(duì)搖擺譜的選擇需要符合船舶航行的海浪實(shí)況,實(shí)測數(shù)據(jù)表明搖擺激勵(lì)越小,觀測量越小,變形結(jié)果耦合的地速就越大,精度越低。表1中船體靜、動(dòng)態(tài)變形角示意圖如圖1所示。

圖1 船體變形角示意圖

2.2 實(shí)驗(yàn)與分析

本節(jié)通過兩種典型的非高斯噪聲來對(duì)比驗(yàn)證MCKF、MAKF、TAKF、MCCKF等4種濾波器的船體變形角估計(jì)性能:沖擊噪聲和高斯混合噪聲,又稱強(qiáng)拖尾噪聲。此外,為確保角度估計(jì)長時(shí)間的收斂穩(wěn)定性,整個(gè)濾波時(shí)間較長;且為避免角度估計(jì)結(jié)果的偶然性,所有實(shí)驗(yàn)的RMSE值均采用50次蒙特卡羅獨(dú)立重復(fù)實(shí)驗(yàn)計(jì)算其均值,并統(tǒng)一從濾波穩(wěn)定之后開始計(jì)算(5 min)。

2.2.1 實(shí)驗(yàn)一: 沖擊噪聲

假設(shè)量測量Z每個(gè)通道均受到如下非高斯分布的劇烈沖擊噪聲(單位:rad/s):

Δvshot=5×10-4,t≥100 s

(28)

該沖擊信號(hào)次數(shù)設(shè)置為1 000次,且隨機(jī)施加給系統(tǒng)。單通道噪聲示意圖如圖2所示,4種濾波器各軸的變形角跟蹤曲線和誤差曲線分別如圖3和圖4所示。

圖2 沖擊噪聲

圖3 變形角跟蹤曲線(沖擊)

圖4 變形角估計(jì)誤差(沖擊)

各軸變形角估計(jì)誤差的RMSE值如表2所示。

表2 船體變形估計(jì)RMSE(沖擊)

從圖3和圖4可以看出,在系統(tǒng)受到隨機(jī)的沖擊噪聲下,傳統(tǒng)的MCKF濾波效果最差,在規(guī)定時(shí)間內(nèi)甚至無法收斂,這是因?yàn)殚L時(shí)間的濾波誤差積累導(dǎo)致發(fā)散;MAKF可以收斂,但由于自適應(yīng)因子的引入時(shí)機(jī)可能有誤,導(dǎo)致個(gè)別通道(橫扭角)狀態(tài)估計(jì)誤差跳動(dòng)范圍大,穩(wěn)定性較差;TAKF解決了這一問題,但濾波開始3 min內(nèi)角度跟蹤誤差曲線出現(xiàn)巨大跳變;而MCCKF在穩(wěn)定性上顯著優(yōu)于前者,但快速收斂之后抗隨機(jī)噪聲干擾的能力要稍弱于自適應(yīng)濾波器。根據(jù)表2顯示,TAKF在精度上具有微弱優(yōu)勢,穩(wěn)定后的RMSE值要略大于前者0.4″,這在更為強(qiáng)調(diào)穩(wěn)定性的工程上幾乎忽略不計(jì)?？梢哉J(rèn)為,在收斂精度滿足要求時(shí),可以選擇MCCKF以期實(shí)現(xiàn)更穩(wěn)定、更快的變形測量。

2.2.2 實(shí)驗(yàn)二: 混合噪聲

高斯混合噪聲,即將兩種非零均值、非高斯分布的噪聲混合,同樣施加給量測信息(均值單位:rad/s;方差單位:(rad/s)2),表達(dá)式為

(29)

式中：N(·，·)為正態(tài)分布。為方便起見,其中相關(guān)參數(shù)設(shè)置下:

(30)

(31)

不同于實(shí)驗(yàn)一隨機(jī)施加的沖擊噪聲,混合噪聲存在于整個(gè)實(shí)驗(yàn)過程,其單通道示意圖如圖5所示,各濾波器的變形角跟蹤曲線如圖6所示。

圖5 混合噪聲

圖6 變形角跟蹤曲線(混合)

相應(yīng)的誤差收斂曲線也給出,如圖7所示。

圖7 變形角估計(jì)誤差(混合)

收斂曲線表明自適應(yīng)濾波器和MCCKF均可快速達(dá)到穩(wěn)定狀態(tài),但明顯可見,MCCKF誤差更小,波動(dòng)閾值要小于MAKF和TAKF。

各軸變形估計(jì)的RMSE結(jié)果如表3所示。

表3 船體變形估計(jì)RMSE(混合)

混合噪聲擾動(dòng)下,與實(shí)驗(yàn)一不同的是,此時(shí)MCCKF的收斂精度要高于TAKF,這是因?yàn)榛旌显肼曉谡麄€(gè)濾波過程均存在,自適應(yīng)濾波器在這種情況下不具備自適應(yīng)抗隨機(jī)擾動(dòng)優(yōu)勢,而MAKF不存在判斷失誤的情況因此并未出現(xiàn)誤差跳變的情形。從實(shí)驗(yàn)結(jié)果可以看出,MCCKF在收斂精度和穩(wěn)定性上也要比MAKF、TAKF優(yōu)1″左右。此外,當(dāng)逐步加大噪聲幅值,MCCKF的RMSE值要比TAKF小至數(shù)個(gè)角秒,表明強(qiáng)非高斯環(huán)境下,MCCKF具備更明顯的優(yōu)勢。

2.2.3 實(shí)驗(yàn)三: 常值靜態(tài)變形角

此外,針對(duì)船體靜態(tài)變形角建模為常值的情形,進(jìn)一步驗(yàn)證MCCKF濾波算法的魯棒性。靜態(tài)變形角大小均設(shè)置為0.1°,以沖擊噪聲為例,靜態(tài)變形角示意圖、變形角跟蹤曲線分別如圖8和圖9所示。

圖8 常值靜態(tài)變形角

圖9 變形角跟蹤曲線(常值靜態(tài)變形角+沖擊噪聲)

相應(yīng)的誤差收斂曲線、變形角估計(jì)誤差的RMSE值分別如圖10和表4所示。

圖10 變形角估計(jì)誤差(常值靜態(tài)變形角+沖擊噪聲)

表4 船體變形估計(jì)RMSE

當(dāng)靜態(tài)變形角建模為常值情形時(shí),實(shí)驗(yàn)結(jié)果再次驗(yàn)證了實(shí)驗(yàn)一的結(jié)論,但MCKF有收斂趨勢,延長濾波時(shí)間或可收斂。實(shí)驗(yàn)表明本文介紹的MCCKF收斂性能始終優(yōu)于MCKF,且對(duì)不同的變形測量模型也具有較好的魯棒性。

此外,當(dāng)加大實(shí)驗(yàn)三的沖擊噪聲頻率次數(shù),MCCKF的RMSE值與TAKF差距越來越小,再次驗(yàn)證MCCKF相對(duì)于自適應(yīng)濾波器,對(duì)強(qiáng)非高斯環(huán)境濾波更加具備優(yōu)勢。

最后,從算法本身來講,基于權(quán)重矩陣的MCCKF相較于兩種自適應(yīng)濾波器MAKF、TAKF,迭代算法簡潔和直觀,無需設(shè)置判斷語句,更方便編程。并且,后者涉及自適應(yīng)漸消因子縮放系數(shù)的選取,選取結(jié)果直接影響自適應(yīng)濾波性能?？紤]到最大熵濾波器對(duì)隨機(jī)外擾、強(qiáng)非高斯均具備不亞于自適應(yīng)濾波器的性能,本文介紹的MCCKF通用性也更強(qiáng)。

3 結(jié) 論

本文介紹了一種基于最大熵準(zhǔn)則和權(quán)重矩陣核函數(shù)的改進(jìn)卡爾曼濾波器,并成功應(yīng)用于非高斯白噪聲環(huán)境下的船體變形測量。重點(diǎn)介紹了基于最大熵的卡爾曼濾波器算法迭代公式,并對(duì)比兩種典型的自適應(yīng)濾波器和傳統(tǒng)最大熵濾波器來驗(yàn)證MCCKF濾波性能。各項(xiàng)實(shí)驗(yàn)的結(jié)果證實(shí),新算法在濾波精度和收斂速度上滿足相關(guān)收斂指標(biāo)前提下,還可更好地解決噪聲模型失配和非高斯噪聲環(huán)境下的變形角估計(jì)問題,且在不同實(shí)驗(yàn)條件下仍保持良好的魯棒性。此外,算法不僅可直接應(yīng)用于復(fù)雜干擾、噪聲建模失準(zhǔn)環(huán)境下的光纖陀螺組件船體變形測量問題,對(duì)于其他用到卡爾曼濾波算法的實(shí)際工程問題也具有一定的指導(dǎo)意義,如運(yùn)載體跟蹤、導(dǎo)航和定位,組合導(dǎo)航的組合濾波算法等。