蔡保杰, 邵 雷

(空軍工程大學防空反導學院, 陜西 西安 710051)

0 引 言

導航系統(tǒng)與日常生活關系越來越密切,無論是車載導航還是飛機、航天器等都需要精確的定位。當行駛到隧道或高樓林立的街道,或是在現(xiàn)代戰(zhàn)爭中衛(wèi)星信號被大功率雜波干擾,都會使衛(wèi)星導航觀測值出現(xiàn)故障,對導航精度造成嚴重影響[1-2],抗差濾波便是在這樣的背景下被提出的,Krarup[3]等人提出理論,Caspary[4]完善了這一理論。

抗差濾波大概分為兩種:等價權函數(shù)法[5]和卡方檢驗法[6]。等價權函數(shù)法通過權函數(shù)來構造相應的等價權,對觀測噪聲方差陣進行降權處理,經(jīng)常與自適應濾波[7-9]配合使用;卡方檢驗法的原理是服從高斯分布的狀態(tài)量的馬氏距離的平方服從卡方分布,通過設定閾值,對故障進行檢測。在這兩種抗差濾波方法上也有很多拓展:文獻[10]中分析了由于觀測值之間的相關性而導致的粗差誤判問題,判斷整體模型是否存在異常。文獻[11]提出了錯誤預警率和故障探測率的概念,對故障的定位更加精確。文獻[12]根據(jù)殘差向量高斯分布統(tǒng)計特性對故障進行檢測,對增益矩陣進行加權處理來抵消異常觀測值影響。

上述抗差濾波算法都可以較好地解決觀測值故障的問題,但也有可改進之處:兩段判別域的故障檢測方法中,被判定為無故障的區(qū)域解算出的導航信息也會出現(xiàn)偏差;當觀測值出現(xiàn)超差情況時,當前時刻的觀測值已經(jīng)沒有利用價值,簡單的降權處理并不能保證觀測值處在合理的范圍內(nèi),應作廢并構造一個新的觀測值進行替換。據(jù)此,本文提出基于三段判斷域的故障檢測方法和最小二乘擬合的抗差濾波算法(以下簡稱擬合法),將檢測區(qū)間分為三段:無故障、偏差和超差,當衛(wèi)星導航觀測值出現(xiàn)偏差時,對故障觀測值進行降權處理;當檢測到當前時刻衛(wèi)星導航觀測值有超差情況時,利用前n個時刻的衛(wèi)星導航觀測值構造的擬合函數(shù)進行一個時刻的外延,替換當前時刻的故障觀測值。相對于傳統(tǒng)的兩段判別域的抗差濾波算法,三段判別域和基于最小二乘擬合的抗差濾波算法在導航精度和穩(wěn)定性方面均有所提高。

1 以組合導航為模型的卡爾曼濾波算法

1.1 卡爾曼濾波

卡爾曼濾波[13]是20世紀60年代由數(shù)學家Kalman等人提出的基于線性最小方差估計的一種遞推算法。

在許多的工程實踐中,準確得到所需狀態(tài)變量的真值是不可能的,只能通過觀測信號和運動模型對狀態(tài)變量進行估計和預測,來無限地接近真值,這也是卡爾曼濾波的基本思想[14-16]?？柭鼮V波迭代流程如圖1所示。

圖1 經(jīng)典卡爾曼濾波流程圖

1.2 松組合模型

組合導航一般為慣性導航系統(tǒng)(inertial navigation system,INS)和全球衛(wèi)星導航系統(tǒng)(global navigation satellite system,GNSS)[17-18]進行信息融合組成的導航系統(tǒng)[19-20],其組合方式大概分為3種:松組合、緊組合和深組合,基于算法層面的是松組合和緊組合。緊組合與松組合相比,狀態(tài)量和觀測量多了偽距和偽距率信息。松組合狀態(tài)向量為

(1)

組合導航系統(tǒng)是一個復雜的系統(tǒng),最簡單的松組合系統(tǒng)也需要多信息融合才能運行,信號傳遞和數(shù)據(jù)處理過程中難免產(chǎn)生一些無用的噪聲,所以導航系統(tǒng)一般使用卡爾曼濾波進行降噪處理[22]。以誤差信息為狀態(tài)量的松組合系統(tǒng)是一個線性系統(tǒng),利用基本卡爾曼濾波便可以得到較好的效果。

2 基于卡方檢驗的故障檢測方法

2.1 基于馬氏距離的卡方檢驗法

卡爾曼濾波是一個迭代過程,殘差向量包含有當前時刻全部的觀測信息。殘差向量表達式為

(2)

(3)

式中,μ取零向量;DM的平方服從卡方分布:

(4)

根據(jù)假設檢驗[24-26]的方法,設原假設為H0,假設導航系統(tǒng)觀測值沒有發(fā)生故障,也就是假設殘差向量服從均值為零的高斯分布:

1.4 統(tǒng)計學方法 采用SPSS 22.0統(tǒng)計軟件分析數(shù)據(jù)。計量資料以x±s表示，采用t檢驗；計數(shù)資料以百分比表示，采用χ2檢驗。以P<0.05為差異有統(tǒng)計學意義。

H0:DMk～χ(m)

(5)

式中，DMk為R時刻殘差向量馬氏距離的平方值。設顯著性水平為α,即認為導航系統(tǒng)沒有發(fā)生故障的概率為α,設

(6)

式中,χa(·)是顯著水平為α的卡方分布。T(m)為檢測量容許的最大值,也稱閾值。當檢測量DM

若當前時刻檢測到故障時,可以將當前時刻的所有觀測量當作故障值,但有可能會將其他正常的觀測值當作故障值處理,造成的誤差會在迭代的過程中逐漸積累,造成濾波發(fā)散。為防止這種現(xiàn)象的發(fā)生,本文采用單維度的故障檢測方法[27],將設定的閾值變?yōu)?/p>

(7)

檢測量的值變?yōu)?/p>

(8)

式中,i表示狀態(tài)向量的第i維。

當檢測量DM(i)

2.2 閾值的選取

DM(i)是根據(jù)當前時刻殘差向量、估計誤差方差陣等計算所得,而閾值T(1)需要根據(jù)卡方統(tǒng)計表選取,如表1所示。閾值選取的準確性決定了抗差濾波的準確性,表征對故障的容忍程度。閾值選取過大,含有較大粗差的觀測值不容易被檢測到,使用抗差濾波前后幾乎無差別;閾值選取過小,一些相對正常的觀測值會被當做故障觀測值來處理,產(chǎn)生的誤差經(jīng)過迭代也會造成濾波發(fā)散。

如表1所示是顯著性水平α和χ2(1)的分布表，從表中可以看出，顯著性水平α越小，認為GNSS觀測值發(fā)生故障的概率就越大，χ2(1)的值就越大，判斷條件就越寬松。本文通過大量的實驗，由檢測量DM(i)的分布圖分析得，檢測量DM(i)的大小與觀測噪聲方差陣Rk和當前時刻粗差大小有關。如圖2～圖6所示為經(jīng)度檢測量分布圖，這里仿真時長共500 s，設定每隔10 s加入粗差，提取這些時刻的卡方檢測量的值作圖。從圖2～圖6中可以看出，當誤差小于5 m時，檢測量數(shù)值有90%的概率小于2.706，當誤差大于25 m時，檢測量數(shù)值有96%的概率大于6.635。當誤差為5～25 m時，檢測量的值位于2.706～6.635之間，這個區(qū)間內(nèi)有一段模糊區(qū)，檢測量DM(i)穿插在卡方分布分位數(shù)的幾個數(shù)值當中，僅用一個閾值判定故障與否很容易出現(xiàn)誤判、漏判的情況。

圖2 誤差為5 m時卡方檢測量分布圖

圖3 誤差為10 m時卡方檢測量分布圖

圖4 誤差為15 m時卡方檢測量分布圖

圖5 誤差為20 m時卡方檢測量分布圖

圖6 誤差為25 m時卡方檢測量分布圖

3 三段判別域和最小二乘擬合的抗差濾波算法

擬合就是從一些數(shù)據(jù)點集(xi,yi)中找出規(guī)律，構造一個函數(shù)關系y=P(x)將這種規(guī)律表現(xiàn)出來，這個函數(shù)關系曲線要求盡可能靠近數(shù)據(jù)點，相比于插值[28]，擬合曲線不嚴格要求通過數(shù)據(jù)點[29-30]。Matlab軟件提供了冪函數(shù)的最小二乘曲線擬合的功能，利用polyfit函數(shù)可以很快地求解出各個項的系數(shù)。

圖7 三段判別域最小二乘擬合抗差濾波算法流程圖

4 仿真分析

本文模擬一條飛機飛行軌跡進行仿真分析，飛行軌跡如圖8所示，將飛行軌跡數(shù)據(jù)送入導航數(shù)據(jù)生成器，得到帶有誤差的模擬真實飛行環(huán)境的導航數(shù)據(jù)，仿真步長為1 s，總時長為501 s。慣性器件誤差較小時的線性系統(tǒng)，一般僅利用位置觀測值便可達到較好的效果，為簡化模型，本文的仿真分析不使用速度觀測值。

圖8 模擬飛機飛行軌跡圖

4.1 最小二乘擬合法抗差效果分析

擬合法是利用前n個時刻的GNSS觀測值構造擬合曲線，進行一個時刻的外延，得到當前時刻的GNSS觀測值。一般來說，n值越大，擬合效果越好，但計算量也會隨之增大，本文經(jīng)過大量實驗分析得出n=4時擬合效果較好，計算量也不會很大。

為營造一個衛(wèi)星信號丟失的環(huán)境，隨機產(chǎn)生100個范圍在0～501之間的整數(shù)作為觀測值故障時刻，對這些時刻的GNSS東向、北向加入100 m的觀測噪聲，高度加入20 m觀測噪聲。

如圖9和圖10所示，一般的濾波方法不能消除觀測值粗差帶來的影響，在觀測值發(fā)生故障的時刻，位置誤差和速度誤差都有較大波動。如圖11和圖12所示，INS在短期內(nèi)的導航效果不錯，誤差較小，但由于INS積分誤差會逐漸積累，隨著時間的推移誤差會越來越大。相比之下，使用最小二乘擬合的抗差濾波算法可以很好地減小故障觀測值帶來的影響，導航精度保持在一個良好的范圍內(nèi)。但擬合法抗差濾波也存在誤差值高于INS的情況，這是由于基于擬合法的抗差濾波雖然可以修正故障觀測值，但比較依賴于構造的函數(shù)模型，在某個時刻可能會出現(xiàn)較大誤差?？柭鼮V波是一個迭代的過程，每一次迭代都需要上個時刻的估計值作為輸入，經(jīng)過觀測值的修正，曲線會逐漸趨于平穩(wěn)?？梢宰C明，最小二乘擬合的抗差濾波算法可以較為準確、穩(wěn)定地輸出導航信息。

圖9 擬合法抗差濾波與無抗差濾波位置誤差

圖10 擬合法抗差濾波與無抗差濾波速度誤差

圖11 擬合法抗差濾波與純慣導位置誤差

圖12 擬合法抗差濾波與純慣導速度誤差

4.2 三段判別域的故障檢測方法仿真分析

當檢測量的值小于2.706時，觀測值誤差小于5 m，認為觀測值無故障，輸出的導航信息也在可容忍的范圍內(nèi)。當檢測量的值在2.706～6.635之間時，在兩段判別域的故障檢測方法中判定為無故障，在三段判別域的故障檢測方法中判定為出現(xiàn)偏差，使用對觀測值降權的方法來減弱故障的影響，以加入的觀測值誤差為20 m為例，如圖13和圖14所示。從圖13和圖14中可以看出，當觀測值誤差為20 m時，如不使用抗差濾波算法解算出的導航信息還是存在一定誤差的，相對而言使用降權法抗差濾波時誤差較小，具體誤差范圍如表2所示。當觀測值誤差為20 m時，檢測量處在2.706～6.635這一區(qū)間的既有故障觀測值，也有無故障觀測值，但這些無故障觀測值經(jīng)過降權法抗差濾波之后的導航解算值仍處在合理范圍內(nèi)，所以在偏差情況使用降權法進行抗差濾波處理可以得到一個較好的結果。

圖13 偏差情況降權法和無抗差濾波位置誤差比較

圖14 偏差情況降權法和無抗差濾波速度誤差比較

抗差濾波無抗差濾波東向位置誤差/m-11.8^6.06-25.8^32.6北向位置誤差/m-6.64^7.54-20.5^34.9天向位置誤差/m-2.62^2.57-7.9^13.9東向速度誤/(m/s)-0.3^0.3-1.16^1.11東向速度誤差/(m/s)-1.3^1.2-1.8^5.2東向速度誤差/(m/s)-0.2^0.3-0.8^1.11

當檢測量的值大于6.635時，觀測值誤差大于25 m，認為出現(xiàn)超差的情況，用基于最小二乘擬合的抗差濾波算法進行修正。當出現(xiàn)超差的情況時，理論上也可使用降權法對故障觀測值進行修正。如圖15和圖16所示為每隔10個時刻觀測值加入100 m的誤差，擬合法與降權法的位置誤差和速度誤差圖。從圖15、圖16及表3中可以看出，當故障時刻為單點時刻時，降權法和擬合法都可以減弱故障觀測值帶來的影響，由于擬合法構造擬合函數(shù)時的計算誤差較大，降權法的抗差效果略好于擬合法，適用于丟包的情況。

表3 觀測量誤差為100 m時導航誤差范圍

但降權法也存在以下問題:當發(fā)生故障的時刻為連續(xù)的一段時間時,降權法的抗差效果并不理想。如圖17和圖18所示為10～200 s這一連續(xù)時間內(nèi)觀測值加入100 m的誤差時,擬合法與降權法的位置誤差和速度誤差圖。從圖17和圖18中可以看出,在10～200 s這一區(qū)間內(nèi),使用降權法進行抗差濾波,誤差從0開始慢慢增大,隨后誤差值開始振蕩,200 s后觀測值正常,誤差趨于平穩(wěn)。相比之下,基于最小二乘擬合的抗差濾波算法在整個仿真區(qū)間位置誤差和速度誤差都比較穩(wěn)定,抗差效果較好。歸結原因主要為:擬合法在觀測值發(fā)生故障時,將故障的GNSS觀測值進行替換,使之處在相對合理的范圍內(nèi),是從根源上解決觀測值故障的一種方法,誤差相對固定。降權法通過增大觀測噪聲方差陣Rk的值來減小濾波增益Kk的值,而Kk的值決定著估計均方誤差陣Pk的大小,Pk經(jīng)過迭代又作用于下一時刻的濾波增益值Kk+1,在連續(xù)時間的迭代中,會使增益陣漸漸失去合適的加權作用,造成發(fā)散。在單時刻故障中,當前時刻計算出的異常的Pk、Kk值會隨著接下來幾個時刻的正常值趨于正常,故導航效果良好。

圖15 單時刻故障擬合法與降權法位置誤差

圖16 單時刻故障擬合法與降權法速度誤差

圖17 連續(xù)時刻故障擬合法與降權法位置誤差

圖18 連續(xù)時刻故障擬合法與降權法速度誤差

5 結 論

針對導航系統(tǒng)在工作過程中GNSS觀測值出現(xiàn)故障的問題,本文提出一種三段判別域的故障檢測方法和基于最小二乘擬合的抗差濾波算法。三段判別域?qū)⒐收戏譃闊o故障、偏差和超差3種情況,若判定觀測值故障為偏差,用對觀測值降權的方法進行抗差濾波處理,若故障為超差,用最小二乘擬合的方法進行抗差濾波。實驗結果表明,三段判別域中對偏差情況進行降權處理比兩段判別域中不做處理導航精度高。在超差情況下,單時刻故障時,使用降權法比擬合法精度略高,適用于丟包的情況。連續(xù)時間故障時使用擬合法比降權法導航精度高且更具穩(wěn)定性,適用于鏈路中斷的情況。