EMD濾波在捷聯(lián)慣導(dǎo)晃動基座初始對準(zhǔn)中的應(yīng)用

于曉雪，張志鑫，夏金橋

(1.北京航天時(shí)代光電科技有限公司，北京 100094；2.北京航天控制儀器研究所，北京 100039)

0 引言

初始對準(zhǔn)是捷聯(lián)慣導(dǎo)系統(tǒng)進(jìn)行導(dǎo)航的前提和基礎(chǔ)。初始對準(zhǔn)一般分為2個(gè)階段，先進(jìn)行粗對準(zhǔn)再進(jìn)行精對準(zhǔn)。在實(shí)際環(huán)境中，影響初始對準(zhǔn)性能的因素，一方面是對準(zhǔn)算法；另一方面是外界噪聲，如車載基座經(jīng)常受到的影響有：車輛的振動、風(fēng)吹導(dǎo)致車輛的晃動、人員上下車、開關(guān)車門等，都會造成較大的噪聲。

傳統(tǒng)的解析粗對準(zhǔn)方法[1]抗晃動能力差，對準(zhǔn)精度低，原因是晃動干擾角速度遠(yuǎn)大于地球自轉(zhuǎn)角速度，無法直接根據(jù)陀螺和加速度計(jì)的輸出計(jì)算姿態(tài)陣?；谀虘T性系的粗對準(zhǔn)方案[2]可以有效地抑制角運(yùn)動干擾對對準(zhǔn)精度的影響，但是無法抑制線運(yùn)動的干擾。胡華峰提出了對慣性坐標(biāo)系下的參考矢量進(jìn)行多項(xiàng)式擬合的抗干擾對準(zhǔn)算法[3]，嚴(yán)恭敏等[4]提出了基于重力加速度積分的解析粗對準(zhǔn)方法，但這兩種方法抑制線運(yùn)動干擾對對準(zhǔn)精度的影響有限；趙長山等[5]提出了抗干擾重力加速度積分解析粗對準(zhǔn)算法，將重力加速度矢量投影到慣性坐標(biāo)系下來隔離角運(yùn)動的干擾，對包含線運(yùn)動干擾的矢量進(jìn)行最小二乘擬合，以抑制線運(yùn)動干擾；練軍想[6]提出了利用多級FIR數(shù)字濾波器進(jìn)行濾波，抑制噪聲干擾，在晃動基座下進(jìn)行初始對準(zhǔn)，但是根據(jù)實(shí)際的應(yīng)用環(huán)境需要設(shè)置相應(yīng)的參數(shù)，適應(yīng)性差。頻域?yàn)V波方案主要通過IIR或FIR濾波[7-9],對高頻線運(yùn)動干擾的抑制效果較好，但是無法抑制異常噪聲的干擾，沒有考慮階躍響應(yīng)的時(shí)間問題，不利于初始對準(zhǔn)方案的快速性應(yīng)用。

近年來小波變換因其良好的時(shí)頻特性和多分辨率能力在慣性技術(shù)領(lǐng)域應(yīng)用越來越多[10-12]。但是小波降噪方法不能夠抑制異常噪聲的干擾，去噪的結(jié)果受閾值和閾值函數(shù)的影響較大，并且不能滿足實(shí)時(shí)、在線降噪處理的要求。文獻(xiàn)[13]提出了一種自適應(yīng)的信號分析方法——經(jīng)驗(yàn)?zāi)B(tài)分解法，該分解方法的時(shí)間尺度為極值點(diǎn)間隔時(shí)間尺度，可以直接從待分解信號的波形當(dāng)中得到，不需要設(shè)置先驗(yàn)分解函數(shù)基，不需要借助窗函數(shù)，具有自適應(yīng)性，無需遞歸運(yùn)算，適合對非平穩(wěn)信號進(jìn)行濾波；并且本質(zhì)上為時(shí)域處理過程，不需要考慮時(shí)間延遲的問題，相對于小波濾波計(jì)算量大幅降低,被廣泛地應(yīng)用在非平穩(wěn)信號去噪中。

因此，本文利用EMD法進(jìn)行預(yù)濾波以消除線運(yùn)動干擾的影響，采用基于凝固慣性系的粗對準(zhǔn)方案抑制角運(yùn)動干擾對對準(zhǔn)精度的影響。

1 坐標(biāo)系定義

1)地心地球坐標(biāo)系(e系)：原點(diǎn)位于地心，xe軸位于赤道平面內(nèi)，從地心指向載體所在點(diǎn)的子午線，ze軸沿地球自轉(zhuǎn)角速度方向，ye軸與xe軸、ze軸構(gòu)成右手坐標(biāo)系。

2)導(dǎo)航坐標(biāo)系(n系)：用于導(dǎo)航解算的坐標(biāo)系，選東北天。

3)載體坐標(biāo)系(b系)：原點(diǎn)位于載體質(zhì)心，xb軸、yb軸、zb軸分別指向載體的右前上。

4)慣性坐標(biāo)系(i系)：xi軸在赤道平面內(nèi)且指向春分點(diǎn)，zi軸指向地球自轉(zhuǎn)，三軸構(gòu)成右手坐標(biāo)系。

5)載體慣性坐標(biāo)系(b0系)：在初始對準(zhǔn)時(shí)刻，將b系慣性凝固后的坐標(biāo)系。

6)導(dǎo)航慣性坐標(biāo)系(n0系)：在初始對準(zhǔn)時(shí)刻，將n系慣性凝固后的坐標(biāo)系。

2 抗干擾初始對準(zhǔn)總體方案設(shè)計(jì)

為了減少線運(yùn)動干擾的影響，設(shè)計(jì)基于EMD法的對準(zhǔn)預(yù)濾波方案；為了減少角運(yùn)動干擾的影響，采用基于凝固慣性系的粗對準(zhǔn)方案，精對準(zhǔn)選擇Kalman濾波精對準(zhǔn)方法，初始對準(zhǔn)流程如圖1所示。

圖1 初始對準(zhǔn)流程Fig.1 Process of initial alignment

2.1 粗對準(zhǔn)方案

(1)

(2)

(3)

構(gòu)造如下2個(gè)參考矢量

(4)

(5)

根據(jù)TRIAD法求解姿態(tài)矩陣，計(jì)算方式如下所示

(6)

其中，t1和t2為選取的2個(gè)時(shí)刻，一般為中間時(shí)刻和最后時(shí)刻。

2.2 Kalman濾波精對準(zhǔn)

在晃動基座條件下，載體質(zhì)心平均位置幾乎不變，載體的速度也保持在零值附近搖擺，唯一有明顯變化的是載體的姿態(tài)，會隨著基座的晃動而變化。因此，在晃動基座條件下進(jìn)行精對準(zhǔn)時(shí)，捷聯(lián)慣導(dǎo)的力學(xué)編排要做相應(yīng)的簡化，只考慮速度和姿態(tài)的更新解算，位置不做更新，這與實(shí)際情況相符。

可推導(dǎo)出SINS的誤差方程為

1) 姿態(tài)誤差方程

(7)

其中，δKG為陀螺組件的刻度系數(shù)誤差，δG為陀螺組件的安裝誤差角，具體如下

在慣性導(dǎo)航系統(tǒng)精確標(biāo)定的情況下，δKG和δG的量級都相當(dāng)?shù)?，對φ的影響都可以忽略，因此，姿態(tài)誤差方程可以簡化為

(8)

2) 速度誤差方程

(9)

其中，δKA為加速度計(jì)組件的刻度系數(shù)誤差，δA為加速度計(jì)組件的安裝誤差角，具體如下

在慣性導(dǎo)航系統(tǒng)精確標(biāo)定的情況下，δKA和δA的量級都相當(dāng)?shù)?。因此，速度誤差方程可以簡化為

(10)

(11)

其中，ωg和ωa分別為陀螺組件和加速度計(jì)組件的輸出噪聲的標(biāo)準(zhǔn)差陣。

由于載體的真實(shí)速度為0附近的小值，而速度解算結(jié)果因失準(zhǔn)角φ的影響會呈發(fā)散趨勢，因此解算得到的速度Vn可以近似等價(jià)于速度誤差δVn，即：δVn≈Vn。

因此量測量取為速度誤差δVn，構(gòu)造量測方程為

Zk=HkXk+vk

(12)

其中，vk為量測噪聲，在晃動基座條件下主要受線擾動幅度的影響。

至此，完整的Kalman濾波模型可以表示為

(13)

3 EMD濾波方案

EMD濾波過程如下：設(shè)信號為x(t)，m1為上下包絡(luò)線的平均值，則函數(shù)h10=x(t)-m1，h1為待選擇的本征模態(tài)函數(shù)。當(dāng)h1滿足一定的篩選條件，則為固有的模態(tài)分量，第二次篩選時(shí)，將第一次篩選得到的h10=x(t)-m11作為原始信號。則有

(14)

由式(14)分解得到h1i，i為篩分的次數(shù)，直到第i次的h1i為本征函數(shù)為止，此時(shí)，令c1=h1i，c1為信號的第一個(gè)本征模態(tài)函數(shù)。剩余信號為

r1=x(t)-c1

(15)

將r1重復(fù)上述步驟，得到c2，同理可以得到滿足條件的n個(gè)本征模態(tài)函數(shù)，從高頻到低頻有

(16)

直到rn不滿足分解的條件，即單調(diào)函數(shù)，EMD分解結(jié)束。常用的篩分標(biāo)準(zhǔn)函數(shù)為

(17)

式中，常用sd=0.2～0.3。

經(jīng)過上述過程，得到了n個(gè)本征模態(tài)函數(shù)，低頻段的主要成分是信號，高頻段的主要成分是噪聲。所以，基于EMD方法進(jìn)行濾波的關(guān)鍵就是求出分界點(diǎn)i，使得對于i之后的IMFi+1中信號為主要成分，而在i之前的IMF中噪聲為主要成分。

Boudraa[14]給出了一種基于連續(xù)均方誤差準(zhǔn)則的EMD濾波方法，即

(18)

其中，i=1,…,n-1。

通過式(18)，可以得到分界點(diǎn)i對應(yīng)的第i個(gè)本征模態(tài)函數(shù)，然后利用從i+1個(gè)之后的本征模態(tài)函數(shù)對信號進(jìn)行重構(gòu)

(19)

其中，分界點(diǎn)i定義為

(20)

即找到本征模態(tài)函數(shù)能量的全局極小值IMFi，i作為能量分布突變的分界點(diǎn)。

EMD濾波流程如圖2所示。

圖2 EMD濾波流程Fig.2 Process of EMD filtering

4 EMD分解存在問題及改進(jìn)

4.1 端點(diǎn)效應(yīng)問題

EMD分解過程中，在連接極值點(diǎn)得到包絡(luò)線時(shí)，通常端點(diǎn)處不會剛好是極值點(diǎn)，因此利用三次樣條曲線連接極值點(diǎn)得到上、下包絡(luò)線時(shí)，端點(diǎn)處的上下包絡(luò)會出現(xiàn)發(fā)散的情況，即端點(diǎn)處的包絡(luò)線會出現(xiàn)較大的擬合誤差。所以，隨著IMF不斷的篩分，誤差會不斷積累。

例如：信號的左端點(diǎn)剛好是極大值點(diǎn)，那么上包絡(luò)線以該點(diǎn)為左端起點(diǎn)，則包絡(luò)線擬合誤差??；此時(shí)，既然左端點(diǎn)是極大值點(diǎn)，則它不是極小值點(diǎn)，因此無法確定下包絡(luò)線的左端起點(diǎn)，將會產(chǎn)生較大的擬合誤差，誤差隨著差分過程不斷積累，導(dǎo)致所得分解結(jié)果嚴(yán)重失真。

具體解決：基于極值點(diǎn)的相關(guān)性來預(yù)測端點(diǎn)處的下一個(gè)極大值和下一個(gè)極小值。

1)求信號的每一個(gè)極大值和極小值，qmax(end)和qmin(end)為最右端的極值，求每個(gè)極大值與qmax(end)、每個(gè)極小值與qmin(end)之間的差，它們的絕對值之和為誤差集；

2)求最右端的極值的左邊的單調(diào)性，與所求誤差集中的每一點(diǎn)的左邊的單調(diào)性進(jìn)行比較，有四種情況(增增，增減，減增，減減)，如果單調(diào)性一致(即增增與減減)，則保留，否則刪除；

3)再計(jì)算保留的這些點(diǎn)的單調(diào)性與右邊相鄰的一個(gè)極值點(diǎn)的單調(diào)性的情況，比較單調(diào)性一致(即增增與減減)出現(xiàn)的頻率，選擇出現(xiàn)頻率較高的情況對應(yīng)的那些點(diǎn)；

構(gòu)造式(21)所示的方程組，求解系數(shù)A0、A1、…、Am

(21)

再根據(jù)式(22)來求解預(yù)測的右端點(diǎn)的下一個(gè)極值qmax(end+1)

Amqmax(end)+…+A0qmax(end-m)=qmax(end+1)

(22)

將信號倒置后，同理，可以求出qmin(end+1)、左端點(diǎn)的極值qmax(start-1)和qmin(start-1)。

改進(jìn)算法的關(guān)鍵是預(yù)測信號兩端的極值，所以，通過比較預(yù)測的極值與真實(shí)值間的誤差，就可以判斷出端點(diǎn)效應(yīng)的強(qiáng)弱。

由圖3中可以看出，改進(jìn)后的算法擬合得到的包絡(luò)線較為靠近真實(shí)的值，因此減弱了帶來的端點(diǎn)效應(yīng)。

圖3 新舊包絡(luò)對比圖Fig.3 Comparison of old and new envelopes

4.2 模態(tài)混疊問題

當(dāng)信號的信噪比較低且存在脈沖干擾或異常噪聲時(shí)，信號的極值點(diǎn)分布不均勻，從而產(chǎn)生模態(tài)混疊現(xiàn)象。

基于奇異值分解抑制模態(tài)混疊：

x={x1,x2,…,xL}為含噪聲的一維離散信號，構(gòu)造矩陣如下

(23)

對Dm矩陣進(jìn)行奇異值分解得到

D=USVT

(24)

其中，U為m×m正交矩陣；V為n×n正交矩陣；S為m×n矩陣。

(25)

式(25)中，Λ=diag(λ1,λ2,…,λγ)，從大到小排列，γ為矩陣Dm的秩。將矩陣Λ的前k(k<γ)個(gè)奇異值保留，令其余奇異值為0，再通過逆奇異值分解過程來重構(gòu)信號。

5 實(shí)驗(yàn)及分析

為驗(yàn)證初始對準(zhǔn)方案的有效性，將光纖捷聯(lián)慣導(dǎo)系統(tǒng)搭載在導(dǎo)航車上進(jìn)行車載初始對準(zhǔn)實(shí)驗(yàn)，對準(zhǔn)期間有人員上下車、開關(guān)車門以及導(dǎo)航車處于怠速狀態(tài)，如圖4所示。

圖4 車載初始對準(zhǔn)試驗(yàn)設(shè)備Fig.4 Test equipment of vehicle initial alignment

采集陀螺儀和加速度計(jì)的輸出數(shù)據(jù)，分別用EMD法和改進(jìn)后的EMD法對數(shù)據(jù)進(jìn)行處理，將實(shí)測數(shù)據(jù)在MATLAB中模擬降噪濾波過程，得到結(jié)果如圖5所示。

(a)z軸陀螺輸出濾波結(jié)果對比

(b)z軸加表輸出濾波結(jié)果對比圖5 輸出結(jié)果對比Fig.5 Comparison of filtering output results

表1 濾波結(jié)果對比Tab.1 Comparison of filtering results

從圖5和表1可以看出，濾波取得了明顯的效果，而且改進(jìn)后的EMD濾波法要優(yōu)于原來的EMD濾波法。對濾波算法的濾波時(shí)間進(jìn)行測試，使用原EMD濾波法用時(shí)0.0713s，使用改進(jìn)后的EMD濾波法用時(shí)0.0760s，數(shù)據(jù)處理速度僅稍有降低。綜合考慮，改進(jìn)后的EMD濾波法效果更好。

圖6 航向角初始對準(zhǔn)結(jié)果對比Fig.6 Comparison of heading angle of initial alignment results

再對初始對準(zhǔn)結(jié)果進(jìn)行分析，以航向角為例。從圖6可以看出，經(jīng)過EMD法濾波后，姿態(tài)角的振蕩明顯減小，收斂速度明顯加快。沒有經(jīng)過預(yù)濾波時(shí)，即僅基于慣性系粗對準(zhǔn)和Kalman濾波精對準(zhǔn)，對準(zhǔn)時(shí)間超過200s，航向角的均方差為6.52′；預(yù)濾波后，對準(zhǔn)150s即達(dá)到了較高精度，航向角的均方差為2.85′；未濾波的對準(zhǔn)算法與基于EMD濾波的對準(zhǔn)算法的航向角對準(zhǔn)結(jié)果均值差為4.372′。

6 結(jié)論

本文綜合考慮了晃動基座情況下進(jìn)行初始對準(zhǔn)受到的外界干擾——角運(yùn)動干擾和線運(yùn)動干擾，設(shè)計(jì)了一種基于EMD的預(yù)濾波方法，并進(jìn)行基于凝固慣性系的粗對準(zhǔn)和Kalman濾波精對準(zhǔn)的車載實(shí)驗(yàn)?？梢缘贸鲆韵陆Y(jié)論：

1)提出了一種基于EMD去噪的抗干擾初始對準(zhǔn)算法，針對EMD法存在的2個(gè)問題進(jìn)行了分析并改進(jìn)EMD法。

2)該抗干擾對準(zhǔn)方法無需分析應(yīng)用的具體環(huán)境噪聲特性，無需設(shè)置先驗(yàn)函數(shù)基。

3)能夠?qū)崿F(xiàn)車載晃動基座下的高精度抗干擾快速初始對準(zhǔn)，具有一定的工程應(yīng)用價(jià)值。