基于殘差卡方檢驗(yàn)和動(dòng)態(tài)Allan方差的INS/GPS故障檢測(cè)與定位算法研究

張春熹，盧　鑫，高　爽，王　璐

(北京航空航天大學(xué)儀器科學(xué)與光學(xué)工程學(xué)院，北京 100191)

0　引言

一般的慣性導(dǎo)航系統(tǒng)通常由兩種傳感器組成，即沿正交方向的3軸陀螺和3軸加速度計(jì)，并通過與衛(wèi)星、磁強(qiáng)計(jì)等外部信息的組合，提供精確可靠的導(dǎo)航數(shù)據(jù)。而慣性傳感器作為INS的核心器件極易受到外界的干擾，造成INS性能的下降，包括零偏和標(biāo)度因數(shù)等參數(shù)都會(huì)受到影響，可能短期內(nèi)就會(huì)在姿態(tài)、速度、位置中累積巨大的誤差[1]。因此,需要對(duì)INS的工作狀態(tài)進(jìn)行實(shí)時(shí)的監(jiān)控，當(dāng)其出現(xiàn)故障時(shí)及時(shí)報(bào)警。

一般來(lái)說,故障可分為硬故障和軟故障。硬故障主要通過BIT(Built In Test)來(lái)檢測(cè)[2]；軟故障可以認(rèn)為是由于環(huán)境等因素造成的慣性器件的性能下降即傳感器信號(hào)中的未知改變，及由于衛(wèi)星信號(hào)受到干擾、遮擋或載體高動(dòng)態(tài)運(yùn)動(dòng)造成的衛(wèi)星數(shù)據(jù)的不確定性。這類故障通常表現(xiàn)為陀螺或加速度計(jì)出現(xiàn)異常的偏值，輸出噪聲特性發(fā)生變化或者衛(wèi)星信號(hào)中夾雜有額外地噪聲等[3]。

在系統(tǒng)層面進(jìn)行故障檢測(cè)時(shí)，通常使用卡方檢驗(yàn)法作為故障檢測(cè)算法，這種方法一般是使用狀態(tài)評(píng)價(jià)器并構(gòu)建統(tǒng)計(jì)學(xué)方程來(lái)進(jìn)行系統(tǒng)層面的狀態(tài)監(jiān)測(cè)[4]。這時(shí)由環(huán)境干擾等造成的傳感器性能的改變及用于組合導(dǎo)航的外部信息的故障都會(huì)耦合進(jìn)入統(tǒng)計(jì)特性中，因此使用卡方檢驗(yàn)法可以對(duì)整個(gè)慣導(dǎo)系統(tǒng)的工作狀態(tài)進(jìn)行實(shí)時(shí)的監(jiān)控，當(dāng)出現(xiàn)故障時(shí)及時(shí)報(bào)警，然而該方法并不能夠準(zhǔn)確地識(shí)別出故障發(fā)生的位置。

要想實(shí)現(xiàn)對(duì)故障的準(zhǔn)確定位，需要一種能夠在傳感器層面進(jìn)行故障檢測(cè)的算法。一般來(lái)說，為了分辨和評(píng)價(jià)慣性傳感器的噪聲特性，通常使用功率譜密度(PSD)、均方根(RMS)、Allan方差等方法，這些方法都是以統(tǒng)計(jì)學(xué)的方法分析當(dāng)前數(shù)據(jù)，判斷傳感器的工作狀態(tài)，并對(duì)未來(lái)可能出現(xiàn)的特性改變做出預(yù)測(cè)和辨識(shí)[5]。其中以Allan方差法最具代表性，作為IEEE協(xié)會(huì)認(rèn)可的一種評(píng)價(jià)陀螺性能的方法，被廣泛應(yīng)用于各種慣性傳感器性能的評(píng)估。由于實(shí)際的傳感器故障信號(hào)可能會(huì)在短時(shí)間內(nèi)表現(xiàn)出不穩(wěn)定性，且長(zhǎng)時(shí)間內(nèi)被測(cè)信號(hào)可能會(huì)受到由于環(huán)境變化、自身老化或突然故障造成的自身性能下降的影響，傳統(tǒng)Allan方差法的使用受到了極大的限制。動(dòng)態(tài)Allan方差法的出現(xiàn)使得跟蹤和描述信號(hào)隨時(shí)間變化的特性和各種噪聲系數(shù)成為可能[6]，然而這種方法的問題在于計(jì)算量很大，這對(duì)于小型化的導(dǎo)航計(jì)算機(jī)來(lái)說是一個(gè)巨大的負(fù)擔(dān)，無(wú)法保證故障檢測(cè)的實(shí)時(shí)性。

所以必須綜合利用系統(tǒng)級(jí)和傳感器級(jí)的故障檢測(cè)方法，充分發(fā)揮各方法的優(yōu)勢(shì)，實(shí)現(xiàn)對(duì)故障的及時(shí)檢測(cè)和準(zhǔn)確定位。

1　原理

1.1　基于殘差卡方檢驗(yàn)法的系統(tǒng)故障檢測(cè)

卡方檢驗(yàn)是判斷統(tǒng)計(jì)樣本的實(shí)際觀測(cè)值與理論值之間偏離程度的一種方法，實(shí)際觀測(cè)值與理論值之間的偏離程度越大，則卡方值越大。

根據(jù)所構(gòu)造的隨機(jī)向量的不同，卡方檢驗(yàn)法又可分為殘差卡方檢驗(yàn)法和狀態(tài)卡方檢驗(yàn)法[7]。狀態(tài)卡方檢驗(yàn)法在報(bào)警期間沒有漏檢，但因其計(jì)算量大，報(bào)警延遲高，而且靈敏度隨著濾波的不斷進(jìn)行有所下降；而殘差卡方檢驗(yàn)法報(bào)警延遲小，對(duì)量測(cè)故障比較敏感，可以直接使用Kalman濾波器的計(jì)算結(jié)果從而運(yùn)算量較低，是最為廣泛使用的故障檢測(cè)算法。

殘差卡方檢驗(yàn)法一般有兩個(gè)步驟，首先是通過Kalman濾波的新息迭代來(lái)提供待檢測(cè)信號(hào)，然后再使用閾值函數(shù)來(lái)界定故障和干擾。故障檢測(cè)算法的設(shè)計(jì)就是要盡量將干擾的影響降至最低，并盡最大可能分辨出故障的發(fā)生[8]。

一般基于松組合的INS/GPS系統(tǒng)的狀態(tài)變量可以表示為：

Xk=[θTδvTδpTεTT]T

(1)

狀態(tài)方程可以寫為:

Xk+1=FkXk+Gkωk

(2)

Fk為系統(tǒng)的狀態(tài)轉(zhuǎn)移矩陣，ωk為系統(tǒng)的噪聲過程，Gk視作白噪聲且其方差陣為一常值，為系統(tǒng)的噪聲轉(zhuǎn)移矩陣。

量測(cè)向量及觀測(cè)方程為：

Yk=[(v-vGPS)T(p-pGPS)T]T

(3)

Yk=HkXk+vk

(4)

Yk為INS和GPS的速度和位置差，Hk為觀測(cè)系數(shù)矩陣，vk為觀測(cè)噪聲并假設(shè)其為白噪聲過程。

在k時(shí)刻得到的Kalman濾波狀態(tài)估計(jì)方程為：

(5)

一旦傳感器或外部信息即觀測(cè)值出現(xiàn)故障，都會(huì)與之前的假設(shè)ωk和vk為相互獨(dú)立的白噪聲相矛盾。當(dāng)系統(tǒng)正常工作未出現(xiàn)故障時(shí)，rk應(yīng)該是Gauss白噪聲服從0均值方差vk的正態(tài)分布，且不同時(shí)刻的新息向量相互獨(dú)立：

rk=Yk-HkXk,k-1～N(0,Vk)

(6)

Vk=HkPkHkT+Rk

(7)

為了檢測(cè)出這時(shí)噪聲方差陣的改變，需要構(gòu)建一個(gè)合適的統(tǒng)計(jì)學(xué)方程，選擇數(shù)據(jù)窗口長(zhǎng)度為p，即截取p個(gè)新息的樣本數(shù)據(jù)作為一組，則其總體服從自由度為6p的卡方分布。觀察系統(tǒng)模型方程和Kalman濾波的迭代方程，評(píng)價(jià)函數(shù)可以設(shè)計(jì)為：

(8)

接下來(lái)選擇相應(yīng)的顯著性水平為a(0

(9)

這里根據(jù)自由度6p和系統(tǒng)要求的虛警率參照卡方分布表來(lái)選取合適的閾值ε，當(dāng)評(píng)價(jià)函數(shù)大于閾值時(shí)，說明系統(tǒng)有故障發(fā)生，否則系統(tǒng)正常工作，即：

通過不斷地對(duì)評(píng)價(jià)函數(shù)值與閾值進(jìn)行比較，可以實(shí)現(xiàn)實(shí)時(shí)的故障檢測(cè)。

1.2　基于動(dòng)態(tài)Allan方差法的傳感器故障檢測(cè)

Allan方差法是一種時(shí)域的分析方法，可以將傳感器隨機(jī)模型中的噪聲項(xiàng)量化成偏值不穩(wěn)定性、角度隨機(jī)游走、量化噪聲等[9]。傳統(tǒng)的Allan方差是以采樣周期τ0對(duì)陀螺輸出的角增量進(jìn)行采樣，在連續(xù)采樣N個(gè)數(shù)據(jù)點(diǎn)后，將其分為K組，每組包含M(M<(N-1/2))個(gè)采樣點(diǎn)，每組數(shù)據(jù)的持續(xù)時(shí)間為τ=Mτ0，即相關(guān)時(shí)間。

不同的分組方式對(duì)應(yīng)每一組的平均值為：

(k=1,2,…,K)

(10)

則Allan方差的計(jì)算公式為：

(11)

根據(jù)不同相關(guān)時(shí)間對(duì)應(yīng)的不同Allan方差值，可以繪制出它們的雙對(duì)數(shù)曲線圖，并通過其與原始數(shù)據(jù)中噪聲功率譜對(duì)應(yīng)的關(guān)系，借由曲線擬合得到不同的噪聲系數(shù)，此時(shí)Allan方差可以表示為幾種誤差源方差的平方和：

(12)

相應(yīng)的噪聲系數(shù)N、K、B、Q、R分別代表角度隨機(jī)游走、速率隨機(jī)游走、零偏穩(wěn)定性、量化噪聲、速率斜坡5種噪聲源的系數(shù)，其計(jì)算公式為[10]：

(13)

動(dòng)態(tài)Allan方差是經(jīng)典Allan方差的擴(kuò)展，它的基本思路是分別計(jì)算信號(hào)在不同時(shí)段內(nèi)的Allan方差，并將計(jì)算結(jié)果繪制在同一幅3維圖中。其計(jì)算方式為定義窗口長(zhǎng)度為M的窗函數(shù)PL(t′)，步進(jìn)長(zhǎng)度N(N

y(t,t0)=x(t0)PL(t-t0)

(14)

假設(shè)待處理的樣本總長(zhǎng)度為L(zhǎng)，則y(t,t1)由長(zhǎng)度為M的采樣信號(hào)組成，對(duì)這段信號(hào)進(jìn)行Allan方差分析，隨后依次步進(jìn)至t0+kN(k=1,2,…)點(diǎn)作為下一組待處理樣本的中點(diǎn),繼續(xù)進(jìn)行Allan方差分析，最后將所有的分析結(jié)果以x軸、y軸和z軸分別代表時(shí)間、相關(guān)時(shí)間、Allan方差值的關(guān)系繪制在3維坐標(biāo)系中，就可以得到這段樣本的動(dòng)態(tài)Allan方差分析圖，并可以分別計(jì)算出每次步進(jìn)后該組樣本的噪聲系數(shù)。

2　仿真驗(yàn)證

從圖1、圖2可知，排除500s～600s的GPS故障和1000s～1200s的傳感器故障外，系統(tǒng)運(yùn)行穩(wěn)定，3軸速度誤差幾乎為0，且評(píng)價(jià)函數(shù)的值一直維持在一個(gè)較低的水平即小于20；而在500s～600s時(shí)由于GPS信號(hào)受到干擾，速度誤差值顯著增加，評(píng)價(jià)函數(shù)值也隨之出現(xiàn)了劇烈的跳變；而在1000s～1200s時(shí)，由于慣性傳感器輸出的噪聲參數(shù)發(fā)生改變，影響了濾波新息的卡方分布，造成了評(píng)價(jià)函數(shù)的突變，而且由于Kalman濾波器所建立的噪聲模型與傳感器真實(shí)的噪聲模型存在巨大偏差，故濾波器一直無(wú)法收斂，導(dǎo)致評(píng)價(jià)函數(shù)的值一直無(wú)法達(dá)到一個(gè)比較小的水平，即系統(tǒng)一直處于故障的狀態(tài)，速度誤差圖反映的情況基本與故障發(fā)生的情況相符，1200s后傳感器和GPS恢復(fù)正常狀態(tài)，系統(tǒng)又回到正常工作狀態(tài)。

隨后對(duì)實(shí)驗(yàn)2系統(tǒng)故障時(shí)的運(yùn)行數(shù)據(jù)進(jìn)行動(dòng)態(tài)Allan方差分析，計(jì)算條件選擇矩形窗，并將窗口長(zhǎng)度設(shè)置為50s，窗口中點(diǎn)的步進(jìn)長(zhǎng)度為16s，擬合方式選擇最小二乘擬合,具體對(duì)比如圖3、圖4所示。

由于卡方檢驗(yàn)已經(jīng)將故障出現(xiàn)的時(shí)間段確定在600s～700s和1000s～1200s之間，接下來(lái)需要確定故障出現(xiàn)的位置，僅需在故障發(fā)生前后的局部區(qū)域?qū)Σ煌瑐鞲衅鬟M(jìn)行動(dòng)態(tài)Allan方差分析即可，故選擇500s～800s及950s～1250s的傳感器數(shù)據(jù)進(jìn)行動(dòng)態(tài)Allan方差分析，即圖3中的0s～300s。由圖3(a)可知，950s～1250s的Allan方差分析圖出現(xiàn)了明顯的波動(dòng)，這是其中一段數(shù)據(jù)加速度計(jì)出現(xiàn)了故障導(dǎo)致的，故障時(shí)的噪聲系數(shù)也明顯增大；而由圖3(b)可知，正常工作的傳感器的Allan方差分析圖非常平穩(wěn)，且噪聲系數(shù)保持在一個(gè)平穩(wěn)的范圍內(nèi)，由此可以判斷此次故障應(yīng)該是由GPS信號(hào)的異常造成的。

3　實(shí)驗(yàn)驗(yàn)證

由于慣性傳感器可能受到復(fù)雜外界環(huán)境的干擾，不失一般性選用實(shí)際采集的振動(dòng)條件下的INS的數(shù)據(jù)進(jìn)行故障檢測(cè)算法分析，卡方檢驗(yàn)的結(jié)果如圖5所示。

圖5中可見在0s～300s評(píng)價(jià)函數(shù)值非常小，說明系統(tǒng)在該時(shí)間段工作正常未出現(xiàn)故障，而在300s～600s評(píng)價(jià)函數(shù)出現(xiàn)了不平穩(wěn)的現(xiàn)象，且在600s～700s時(shí)出現(xiàn)劇烈波動(dòng)，并在1000s后恢復(fù)正常，說明INS系統(tǒng)在300s～1000s的時(shí)間段內(nèi)非正常工作。

隨后對(duì)故障出現(xiàn)的200s～1100s的任意一軸陀螺數(shù)據(jù)進(jìn)行動(dòng)態(tài)Allan方差分析，數(shù)據(jù)采樣間隔2.5ms，數(shù)據(jù)總長(zhǎng)度1250s， Allan方差選擇矩形窗并將窗口長(zhǎng)度設(shè)置為50s，窗口中點(diǎn)的步進(jìn)長(zhǎng)度16s，擬合方式選擇最小二乘擬合，具體如圖6、圖7所示。

由圖6、圖7可知，前300s的噪聲系數(shù)非常小，動(dòng)態(tài)Allan方差分析圖也非常穩(wěn)定，說明INS處于正常工作中狀態(tài)良好；而300s后圖像出現(xiàn)明顯的改變，相應(yīng)傳感器的噪聲系數(shù)也出現(xiàn)了輕微的增加并在600s處出現(xiàn)了顯著的增加，說明受到振動(dòng)環(huán)境的干擾傳感器的輸出特性發(fā)生了改變，而且由圖中可見出現(xiàn)了兩種不同分布的噪聲，這是由于振動(dòng)臺(tái)進(jìn)入開機(jī)的狀態(tài)后對(duì)INS有所干擾，隨后距離振動(dòng)真正的開始大約有300s的時(shí)間間隔，無(wú)論在Allan方差分析圖或是在噪聲系數(shù)上都能夠清晰地反映出這幾種狀態(tài)的變化，完全符合實(shí)驗(yàn)條件，也從側(cè)面佐證了殘差卡方檢驗(yàn)配合動(dòng)態(tài)Allan方差法在慣性傳感器故障檢測(cè)方面的有效性。

4　結(jié)論

INS/GPS系統(tǒng)作為一種精密的測(cè)量系統(tǒng)，其內(nèi)部傳感器極易受到外部溫度、濕度、振動(dòng)、輻射等環(huán)境因素的影響及其本身器件老化造成的性能下降，而GPS的信號(hào)也容易受到干擾并不完全可靠，因此提出了綜合運(yùn)用殘差卡方檢驗(yàn)法和動(dòng)態(tài)Allan方差法，分別從系統(tǒng)層面和傳感器層面對(duì)系統(tǒng)故障進(jìn)行檢測(cè)及定位的方法。

其中，殘差卡方檢驗(yàn)法借助Kalman濾波的中間變量構(gòu)建評(píng)價(jià)方程，通過不斷的與閾值進(jìn)行比較，可以實(shí)現(xiàn)實(shí)時(shí)地系統(tǒng)故障診斷，能夠從系統(tǒng)層面實(shí)時(shí)地對(duì)系統(tǒng)性能進(jìn)行監(jiān)測(cè)。當(dāng)傳感器的性能下降或者外部信息出現(xiàn)故障時(shí)，能夠及時(shí)檢測(cè)出來(lái)并準(zhǔn)確判斷故障發(fā)生的時(shí)間，為后續(xù)Allan方差的使用劃出范圍，大大優(yōu)化了Allan方差的計(jì)算量。然而傳統(tǒng)的卡方檢驗(yàn)法是建立在噪聲為Gauss分布的條件下對(duì)系統(tǒng)進(jìn)行檢測(cè)，在實(shí)際工程應(yīng)用中噪聲很難完全滿足理想的Gauss分布，因此判斷故障與否的閾值應(yīng)根據(jù)實(shí)際情況適當(dāng)調(diào)整。隨后的動(dòng)態(tài)Allan方差法在已知故障發(fā)生時(shí)間的基礎(chǔ)上，彌補(bǔ)了殘差卡方檢驗(yàn)法在出現(xiàn)故障位置識(shí)別方面的劣勢(shì)，但是由于該方法本身計(jì)算量較大且受計(jì)算方式的影響需要一定的數(shù)據(jù)長(zhǎng)度來(lái)保證結(jié)果的準(zhǔn)確性，因此實(shí)時(shí)性較差。

最后，借助這兩種方法的綜合使用，充分發(fā)揮各方法的優(yōu)勢(shì)，可以在INS出現(xiàn)故障時(shí)及時(shí)準(zhǔn)確地檢測(cè)出來(lái)，并準(zhǔn)確定位至故障發(fā)生的位置，為INS故障后的修復(fù)提供依據(jù)，確保INS正常工作。

[1]Zhong M Y, Guo J, Guo D F, et al. An extended HiH∞optimization approach to fault detection of INS/GPS-integrated system[J]. IEEE Transactions on Instrumentation and Measurement, 2016, 65(11) :2495-2504.

[2]Pecht M, Dube M, Natishan M, et al. Evaluation of built-in test[J]. IEEE Transactions on Aerospace & Electronic Systems,2001, 37(1):266-271.

[3]Vitanov I, Aouf N. Fault diagnosis for MEMS INS using unscented Kalman filter enhanced by Gaussian process adaptation[C]. 2014 NASA/ESA Conference on Adaptive Hardware and Systems(AHS), 2014：120-126.

[4]Zhong M Y, Guo J, Yang Z H. On real time performance evaluation of the inertial sensors for INS/GPS integrated systems[J]. IEEE Sensors Journal, 2016, 16(17):6652-6661.

[5]Evaluation[EB/OL].http://en.wikipedia. org/wiki/evaluation, 2015.

[6]張娜, 李緒友. 動(dòng)態(tài)Allan方差的理論改進(jìn)及其應(yīng)用研究[J].光學(xué)學(xué)報(bào), 2011, 31(11):70-75.

ZHANG Na，LI Xu-you. Research on theoretical improvement of dynamic Allan variance and itsapplication[J]. Acta Optica Sinica, 2011, 31(11):70-75.

[7]翁浚, 成研, 秦永元,等. 車輛運(yùn)動(dòng)約束在SINS/OD 系統(tǒng)故障檢測(cè)中的應(yīng)用[J].中國(guó)慣性技術(shù)學(xué)報(bào), 2013, 21(3): 406-410.

WENG Jun, CHENG Yan, QIN Yong-yuan, et al. Application of vehicle constraints in SINS/OD system's fault detection[J]. Journal of Chinese Inertial Technology, 2013,

21(3): 406-410.

[8]Hwang D H, Sang H O, Sang J L, et al. Design of a low-cost attitude determination GPS/INS integrated navigation system[J]. GPS Solutions, 2005, 9(4):294-311.

[9]李冀辰, 高鳳岐, 王廣龍, 等. 光纖陀螺振動(dòng)和變溫條件下的DAVAR分析[J].中國(guó)激光, 2013, 40(9):184-190.

LI Ji-chen, GAO Feng-qi, WANG Guang-long, et al. Analysis of dynamic Allan variance for fiber optic gyro under vibration and variable temperature conditions[J]. Chinese Journal of Lasers, 2013, 40(9):184-190.

[10]Galleani L. The dynamic Allan variance III: confidence and detection surfaces[J]. IEEE Transactions on Ultrasonics, Ferroelectric, and Frequency Control, 2011, 58(8): 1550-1558.