基于卡爾曼濾波的MEMS 陀螺噪聲抑制算法設(shè)計(jì)與實(shí)現(xiàn)

韓穎

（中國空空導(dǎo)彈研究院，河南洛陽，471000）

0 引言

將MEMS 器件應(yīng)用在飛航導(dǎo)彈系統(tǒng)中是未來微型化、小型化的發(fā)展趨勢。因技術(shù)受限，MEMS 陀螺儀目前多為商業(yè)級(jí)[1]。對于精度要求高的穩(wěn)像系統(tǒng)，MEMS 陀螺儀噪聲在彈載飛行過程中引起運(yùn)動(dòng)載體輸出軸的抖動(dòng)以及角度的漂移，會(huì)導(dǎo)致視頻圖像質(zhì)量的下降。為了使穩(wěn)像系統(tǒng)有高的跟蹤精度，在本文中設(shè)計(jì)使用卡爾曼濾波算法，對MEMS 陀螺儀的噪聲進(jìn)行抑制或去除。

1 基于時(shí)間序列ARMA 模型的隨機(jī)誤差建模

MEMS 陀螺儀隨機(jī)噪聲信號(hào)屬于非平穩(wěn)的時(shí)間序列，對于載體機(jī)動(dòng)頻繁，不確定因素較多的系統(tǒng)，為避免將真實(shí)數(shù)據(jù)當(dāng)作噪聲過濾掉，在卡爾曼濾波前，應(yīng)對陀螺儀噪聲信號(hào)進(jìn)行預(yù)處理，保證其為零均值、平穩(wěn)、正態(tài)時(shí)間序列后，建立ARMA 模型[2,3]。

1.1 去趨勢項(xiàng)

趨勢項(xiàng)是是MEMS 陀螺儀線性漂移的主要影響因素，一般常采用最小二乘法將其剔除。

在本研究中使用Matlab 運(yùn)算工具中polyfit()函數(shù)和polyval()函數(shù)代替一般代數(shù)方法求得趨勢項(xiàng)系數(shù)，在對陀螺儀數(shù)據(jù)進(jìn)行3 階、4 階、5 階、6 階、7 階多項(xiàng)式趨勢項(xiàng)去除后發(fā)現(xiàn)，去除5 階、6 階和7 階趨勢項(xiàng)后的數(shù)據(jù)方差很接近，見表1。

表1 去趨勢項(xiàng)后數(shù)據(jù)方差

考慮到在實(shí)際應(yīng)用中，模型階數(shù)大，計(jì)算量會(huì)隨之增大。因此，本文選擇采用剔除5 階趨勢項(xiàng)方法使陀螺的隨機(jī)漂移滿足平穩(wěn)性要求。

1.2 去除野值

采集信號(hào)中偶然出現(xiàn)的異常數(shù)據(jù)，被稱為野值（outl ier）?？紤]到計(jì)算量的問題，設(shè)計(jì)中使用3σ 準(zhǔn)則循環(huán)3 次去除野值項(xiàng)，如圖1。

圖1 野值點(diǎn)與3 次去野值后的波形

1.3 零均值檢驗(yàn)

利用Matlab 工具箱中dtrend()函數(shù)應(yīng)對去除野值之后的陀螺儀隨機(jī)噪聲信號(hào)進(jìn)行平穩(wěn)化，零化處理，使噪聲信號(hào)為零均值序列，處理前后均值見表2。

表2 零化處理前后的均值

1.4 平穩(wěn)性檢驗(yàn)

平穩(wěn)性檢驗(yàn)是對時(shí)間序列進(jìn)行ARMA 建模的前提條件?？梢杂脕頇z驗(yàn)陀螺儀漂移噪聲序列是否不隨時(shí)間推移而變化。本設(shè)計(jì)中采用平穩(wěn)性的Daniel 檢驗(yàn)方法對陀螺隨機(jī)噪聲序列進(jìn)行檢驗(yàn)。結(jié)果如表3：

表3 Daniel 平穩(wěn)性檢驗(yàn)方法檢驗(yàn)結(jié)果

1.5 正態(tài)性檢驗(yàn)

總體分布的正態(tài)性檢驗(yàn)一般采取Jarque-Bera 檢驗(yàn)，觀察偏態(tài)系數(shù)g3與峰態(tài)系數(shù)g4是否滿足正太隨機(jī)變量的特性。對兩組去除趨勢項(xiàng)和野值點(diǎn)的陀螺隨機(jī)噪聲信號(hào)做正態(tài)性檢驗(yàn)，其結(jié)果如表4。

表4 正態(tài)性檢驗(yàn)結(jié)果

由上表可以看出，g3、g4都約等于0，因此可以認(rèn)此時(shí)MEMS 陀螺隨機(jī)信號(hào)序列均滿足正態(tài)性要求。

1.6ARMA 模型階數(shù)的確定

將MEMS 陀螺漂移信號(hào)轉(zhuǎn)化為零均值、正態(tài)、平穩(wěn)時(shí)間序列，符合建模要求后，選用AIC 準(zhǔn)則和BIC 準(zhǔn)則對陀螺隨機(jī)信號(hào)序列進(jìn)行ARMA(p，q)模型階數(shù)的確定。

本文可以采用的時(shí)間序列模型一共有五種[4]，即AR(1)、AR(2)、AR(3)、ARMA(1,1)、ARMA(2,1)。對這五種模型同時(shí)進(jìn)行AIC 準(zhǔn)則和BIC 準(zhǔn)則的適用性檢驗(yàn)，檢驗(yàn)結(jié)果如表5。

表5 AIC 準(zhǔn)則和BIC 準(zhǔn)則的檢驗(yàn)結(jié)果

通過比較后，選取AIC 和BIC 值相對較小的ARMA(2,1)模型作為陀螺儀的隨機(jī)漂移誤差模型。其表達(dá)式為：

在Matlab 中使用最小二乘法即可解得到相應(yīng)的ARMA(2,1)的模型參數(shù)，如表6 所示，將此數(shù)據(jù)作為卡爾曼濾波算法的初值。

表6 陀螺儀參數(shù)模型

2 卡爾曼濾波算法設(shè)計(jì)

2.1 卡爾曼濾波原理

卡爾曼濾波是根據(jù)前一個(gè)估計(jì)值和最近一個(gè)觀測值，按一定的遞推方式求得新的狀態(tài)估計(jì)值，然后不斷循環(huán)的過程?；A(chǔ)公式(1)之外，其狀態(tài)方程和測量方程分別如下[5,6]：

其中，R k是測量噪聲協(xié)方差矩陣，為正定矩陣；Q k是過程噪聲協(xié)方差矩陣，為非負(fù)正定矩陣，δ kj為克羅尼克（Kroneker）δ函數(shù)，其特性是：

如果被估計(jì)狀態(tài)和觀測量是滿足式(2)、(3)，則k 時(shí)刻的觀測的估計(jì)可按下述方程求解。

根據(jù)以上五條基本公式，給定初值X0和P0，根據(jù)k時(shí)刻的觀測值Zk，就可以遞推計(jì)算得k 時(shí)刻的狀態(tài)估計(jì)(k=1，2…,N)。

3 卡爾曼濾波軟件設(shè)計(jì)

卡爾曼濾波算法框圖如圖2 所示。

圖2 一步預(yù)測卡爾曼濾波算法

在試驗(yàn)中ARMA(2,1)模型參數(shù)如表6 所示。R k是觀測噪聲序列方差陣，取濾波數(shù)據(jù)的方差varX；Qk是系統(tǒng)噪聲序列方差陣，為正定矩陣，取其中σα為ARMA(2,1)模型系統(tǒng)激勵(lì)白噪聲方差。

在本設(shè)計(jì)中假定濾波估計(jì)從開始就是無偏的，且估計(jì)的誤差陣最小，初始值另有，Φk,k?1和Γk,k?1為tk,k?1時(shí)刻至tk時(shí)刻的狀態(tài)轉(zhuǎn)移矩陣，取

4 Allan 方差隨機(jī)噪聲分析方法

理想情況下，當(dāng)MEMS 陀螺輸入為零時(shí)，輸出也應(yīng)該為零，但是MEMS 陀螺存其它形式的隨機(jī)干擾和漂移[7]。在本研究中，選用IEEE 公認(rèn)的陀螺儀參數(shù)分析評(píng)價(jià)的標(biāo)準(zhǔn)方法——Allan 方差，對MEMS 陀螺儀濾波結(jié)果進(jìn)行評(píng)價(jià)。

5 實(shí)驗(yàn)

在25 ℃室溫靜態(tài)條件下，以500Hz 的采樣頻率對MEMS 陀螺樣品進(jìn)行靜態(tài)數(shù)據(jù)采集，敏感軸指天，陀螺標(biāo)度因數(shù)為7.0mV/°/sec，并假設(shè)其不變。采樣時(shí)間1min。利用Matlab2007a 對采集到的數(shù)據(jù)進(jìn)行預(yù)處理后，方可對陀螺儀輸出的隨機(jī)噪聲信號(hào)序列進(jìn)行濾。濾波前后對比，如圖3 所示。

圖3 卡爾曼濾波前后數(shù)據(jù)以及濾除噪聲

卡爾曼濾波前后Allan 方差對比如圖4 所示。

圖4 經(jīng)典卡爾曼濾波前后Allan 方差對比圖

卡爾曼濾波前后Allan方差分析各噪聲分量的結(jié)果如7。

卡爾曼濾波前后陀螺隨機(jī)誤差標(biāo)準(zhǔn)差如表8 所示。

表8 卡爾曼濾波前后陀螺隨機(jī)誤差標(biāo)準(zhǔn)差

6 結(jié)論

從圖3 及表7 中可以看出，在經(jīng)過卡爾曼濾波處理后，量化噪聲約為濾波前的50%，角速度隨機(jī)游走、零偏不穩(wěn)定性、速率隨機(jī)游走、速率斜坡約為濾波前的60%。在表8 中可以看出，通過卡爾曼濾波之后，陀螺儀隨機(jī)誤差標(biāo)準(zhǔn)差降為濾波前的50%。說明卡爾曼濾波對該MEMS 陀螺儀有明顯的濾波效果。

表7 卡爾曼濾波前后Allan 方差結(jié)果

將卡爾曼濾波算法與MEMS 陀螺儀相結(jié)合，運(yùn)用在彈載系統(tǒng)中，既節(jié)省了存儲(chǔ)空間，又易于工程實(shí)現(xiàn)，使濾波過程更加直接有效，非常適合應(yīng)用在穩(wěn)像系統(tǒng)中對MEMS 陀螺儀的噪聲處理中。