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

蔡保杰, 邵 雷

(空軍工程大學(xué)防空反導(dǎo)學(xué)院, 陜西 西安 710051)

0 引 言

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

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

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

1 以組合導(dǎo)航為模型的卡爾曼濾波算法

1.1 卡爾曼濾波

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

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

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

1.2 松組合模型

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

(1)

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

2 基于卡方檢驗(yàn)的故障檢測(cè)方法

2.1 基于馬氏距離的卡方檢驗(yàn)法

卡爾曼濾波是一個(gè)迭代過程,殘差向量包含有當(dāng)前時(shí)刻全部的觀測(cè)信息。殘差向量表達(dá)式為

(2)

(3)

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

(4)

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

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

H0:DMk～χ(m)

(5)

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

(6)

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

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

(7)

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

(8)

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

當(dāng)檢測(cè)量DM(i)

2.2 閾值的選取

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

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

圖2 誤差為5 m時(shí)卡方檢測(cè)量分布圖

圖3 誤差為10 m時(shí)卡方檢測(cè)量分布圖

圖4 誤差為15 m時(shí)卡方檢測(cè)量分布圖

圖5 誤差為20 m時(shí)卡方檢測(cè)量分布圖

圖6 誤差為25 m時(shí)卡方檢測(cè)量分布圖

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

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

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

4 仿真分析

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

圖8 模擬飛機(jī)飛行軌跡圖

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

擬合法是利用前n個(gè)時(shí)刻的GNSS觀測(cè)值構(gòu)造擬合曲線，進(jìn)行一個(gè)時(shí)刻的外延，得到當(dāng)前時(shí)刻的GNSS觀測(cè)值。一般來(lái)說(shuō)，n值越大，擬合效果越好，但計(jì)算量也會(huì)隨之增大，本文經(jīng)過大量實(shí)驗(yàn)分析得出n=4時(shí)擬合效果較好，計(jì)算量也不會(huì)很大。

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

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

圖9 擬合法抗差濾波與無(wú)抗差濾波位置誤差

圖10 擬合法抗差濾波與無(wú)抗差濾波速度誤差

圖11 擬合法抗差濾波與純慣導(dǎo)位置誤差

圖12 擬合法抗差濾波與純慣導(dǎo)速度誤差

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

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

圖13 偏差情況降權(quán)法和無(wú)抗差濾波位置誤差比較

圖14 偏差情況降權(quán)法和無(wú)抗差濾波速度誤差比較

抗差濾波無(wú)抗差濾波東向位置誤差/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

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

表3 觀測(cè)量誤差為100 m時(shí)導(dǎo)航誤差范圍

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

圖15 單時(shí)刻故障擬合法與降權(quán)法位置誤差

圖16 單時(shí)刻故障擬合法與降權(quán)法速度誤差

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

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

5 結(jié) 論

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