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

蔡保杰, 邵 雷

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

0 引 言

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

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

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

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

1.1 卡爾曼濾波

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

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

圖1 經典卡爾曼濾波流程圖

1.2 松組合模型

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

(1)

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

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

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

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

(2)

(3)

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

(4)

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

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

H0:DMk～χ(m)

(5)

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

(6)

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

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

(7)

檢測量的值變為

(8)

式中,i表示狀態向量的第i維。

當檢測量DM(i)

2.2 閾值的選取

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

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

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

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

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

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

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

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

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

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

4 仿真分析

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

圖8 模擬飛機飛行軌跡圖

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

5 結 論

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