基于穩健MM估計的REKF RAIM算法

王文博, 徐 穎

(1. 中國科學院空天信息創新研究院, 北京 100094;2. 中國科學院大學電子電氣與通信工程學院, 北京 100049)

0 引 言

接收機自主完好性監測(receiver autonomous integrity monitoring, RAIM)利用接收機觀測信息的冗余性進行一致性檢驗來實現對衛星故障的檢測與識別[1-3],是導航系統用戶端完好性監測的主要手段之一。典型的RAIM算法包括僅利用當前時刻觀測信息進行一致性檢驗的“快照式”算法[4-6]和利用當前觀測信息與歷史觀測信息的“濾波式”算法[7-9]。

傳統的RAIM算法主要是基于殘差或新息構造卡方檢驗的檢驗統計量,因此檢驗統計量也是參數估計值的函數,在故障模式下,同樣會受異常觀測值的影響而產生漏檢或誤檢[10]。針對這一問題,有學者提出用抗差估計取代最小二乘估計[11-14],其中基于極大似然估計(maximum likelihood estimate, MLE)的穩健M估計[15]既具備一定的穩健性,又保留了最小二乘估計效率高的優點,是最常用的抗差估計方法。M估計通過設計等價權函數對異常觀測值進行降權或除權處理,以此來削弱異常觀測值對RAIM故障檢測與識別的影響,國內外研究學者提出了大量可供選擇的等價權函數[16-19]。文獻[20]基于M估計進行了RAIM粗差檢測實驗,證明了基于M估計的RAIM具備更高的故障檢測率。文獻[21]基于IGG III方案設計了魯棒擴展卡爾曼濾波(robust extended Kalman filter, REKF)RAIM方法,驗證了基于REKF的RAIM對微小緩變故障有更好的檢測效果。文獻[22]提出了基于魯棒擴展卡爾曼粒子濾波的RAIM算法,有效縮短了對故障衛星的檢測延遲時間。上述文獻多是針對單星故障模式,但M估計的穩健性依賴于迭代初值的可靠性,對于多星故障模式,由于故障對殘差的影響相互耦合,尤其當故障矢量具有較高的空間一致性時,M估計的迭代初值會受到嚴重影響,此時得到的參數估值穩健性較差,對故障的檢測與識別效果也將受到影響。

本文提出基于穩健MM估計[23]的REKF RAIM算法,穩健MM估計器分為兩部分,首先采用具備高崩潰污染率的最小截斷二乘(least trimmed square, LTS)估計[24-25]得到高穩健性的估值初值,然后采用具備高抗差效率的IGG III函數輸出抗差估計結果。另外,針對LTS估計的計算量隨觀測信息數量增加而顯著增大的問題,設計了基于特征斜率的快速選星算法來降低LTS估計器的計算量。仿真結果表明,本文提出的方法能有效提高雙星故障模式下M估計迭代初值的穩健性,對雙星故障有更好的檢測效果和更高的定位精度。

1 REKF RAIM方法

1.1 REKF RAIM模型

線性化后接收機的觀測方程為

Yk=HkXk+εk

(1)

式中,Yk為時刻k的觀測向量;Hk為n×m階觀測矩陣，n為衛星數，m為待估狀態數;Xk為待估計狀態向量；εk為觀測噪聲向量,且εk～N(0,Rk),其中Rk為觀測噪聲的協方差矩陣。

系統的狀態轉移方程為

Xk=Φk,k-1Xk-1+wk-1

(2)

式中,Φk,k-1為狀態轉移矩陣;wk為過程噪聲向量,且wk～N(0,Qk),其中Qk為過程噪聲的協方差矩陣。

REKF估計器包括預測和修正兩個步驟。

預測:

(3)

(4)

(5)

修正:

(6)

Pk=(In-KkHk)Pk,k-1

(7)

新息為觀測值與預測觀測值之差,表示為

(8)

在無故障模式下,新息服從高斯分布,有

rk～N(0,Dk)

(9)

式中,Dk為新息的協方差矩陣

(10)

根據新息構造檢驗統計量

(11)

Tk服從自由度為(n-m)的卡方分布,Tk～χ2(n-m)。在故障模式下,Tk服從非中心卡方分布,Tk～χ2(n-m,λ),λ為非中心化參數。

根據給定的虛警概率Pfa計算檢驗門限TD,即

(12)

1.2 M估計及IGG III等價權函數

當觀測信息中存在異常值時,擴展卡爾曼濾波估計值可能會失真甚至不收斂,因此引入抗差估計理論來抑制異常值的影響是有必要的。REKF RAIM方法的流程如圖1所示。

圖1 REKF RAIM算法流程圖

基于MLE理論的M估計是較為常用的抗差估計器,其本質上是一種加權最小二乘估計,根據標準化殘差對各觀測信息賦予不同的權重,其目標函數可表示為

(13)

式中,si為標準化殘差,si=ri/MAD(r),其中尺度參數MAD(v)為殘差的中位絕對偏差,表示為

MAD(v)=1.483·med(|(v-med(v)|)

(14)

通過設計合適的權函數ρ來抑制異常觀測值的影響,常用的等價權函數有丹麥法、Huber法和IGG法等,本文采用三段式的IGG III等價權函數

(15)

式中,k0和k1為經驗性抗差參數,通常取k0∈[1.5,2.5],k1∈[3,4.5]。IGG III等價權函數包括保權區、降權區和除權區,通過擴大異常值的方差來抑制異常值對參數估計的負面影響。

2 兩步抗差MM估計

M估計通常采用最小二乘估計值作為迭代初值,由于最小二乘估計值對異常觀測值十分敏感,會在一定程度上削弱M估計的抗差能力,尤其是在多星故障模式下,當故障矢量具有較高的空間一致性時,最小二乘估計結果將向故障矢量產生明顯偏移,此時M估計的抗差能力將受到嚴重破壞。本文提出采用MM估計替代M估計的方法,MM估計由LTS估計和M估計兩部分組成,以LTS估計值作為M估計的迭代初值,并設計了基于特征斜率的快速選星方法來提高LTS估計的計算效率。

2.1 LT S估計

崩潰污染率是衡量參數估計方法穩健性的指標,表示該估計方法允許數據中存在的不致使估計值失效的異常值數量的最小比例,估計方法的崩潰污染率越高,其可容忍的異常觀測值越多,得到的估計值穩健性越好。

LTS估計是一種具有高崩潰污染率的穩健估計方法,其目標函數可表示為

(16)

LTS估計的核心思想是使目標函數升序排列的前h個殘差平方和達到最小值[25]。LTS估計在h=int[(n+m+1)/2]時崩潰污染率可達到理論上的最大值50%,因此LTS估計具備更高的穩健性。

LTS估計的步驟如下:

步驟 3取殘差平方和最小的方程組的解作為LTS的估計解。

h≤n-q

(17)

h的取值越小,LTS估計的崩潰污染率越高,但需要解算的方程組的數量也會隨之增大,因此LTS估計的效率并不高。根據完好性風險與故障衛星數量的關系[26],三星及三星以上故障模式的概率小于完好性指標對漏檢概率需求,因此多星故障模式僅需考慮雙星故障模式,故本文取h=n-2。

2.2 基于特征斜率的快速選星方法

特征斜率[27]是描述定位誤差與檢驗統計量之間關系的線性無噪聲模型,各衛星的特征斜率等價于各星對PDOP值的貢獻度

(18)

式中,A=(HTH)-1HT和S=In-H(HTH)-1HT均由觀測矩陣計算得到,因此特征斜率僅與衛星的空間位置有關。

相比于常用的基于精度因子(dilution of precision,DOP)值或幾何分布的選星算法[28-29],以特征斜率作為判斷依據僅需要簡單的數值計算,避免了大量的矩陣運算,具備更高的選星效率,給定終止判定閾值λ,衛星數量閾值nall,單星座衛星數量保護閾值nsin,具體的步驟如下。

算法 1 基于特征斜率的快速選星輸入 H,n輸出 Hsel1 while n>nall do2 根據H計算當前集合的PDOP3 根據式(18)計算各星的特征斜率slope4 排序并選出slopemin對應的衛星5 if 該星所屬星座為受保護星座then6 跳過該星7 else8 刪除H中該星所對應的行得到Hsel9 n=n-110end if11根據Hsel計算PDOPsel12δ=(PDOPsel-PDOP)/PDOP13if δ>λ then14 break15else16 if 該星座剩余衛星數小于保護閾值 then17 將該星座標記為保護星座18 end if19end if20if 全部星被標記為受保護星座 then21 break22end if23 end while

3 仿真驗證

為驗證本文算法的有效性,選取3種參照算法(EKF RAIM、基于M估計的REKF RAIM和基于MM估計的REKF RAIM)對比雙星故障的檢測識別能力。仿真中使用從CELESTRAK下載的全球定位系統(global positioning system,GPS)和北斗三號全球衛星導航系統(簡稱為BDS -3)在 2020年第11天的軌道數據,模擬的目標觀測點為(116°E, 40°N),先驗觀測噪聲均方根為4 m,仿真起始時間為北京時間12:00:01,仿真時長為150 s,仿真步長為1 s,在此期間觀測到8顆GPS衛星(編號1～編號8)和9顆BDS-3衛星(編號9～編號17),在仿真起始時刻全部衛星的位置分布如圖2所示。

圖2 初始時刻衛星分布

由圖2可以看出,4號星和13號星在空間中相對于接收機的位置具有較高的一致性,在加入相同的偽距偏差時,偏差矢量具有很高的空間一致性。

單星座衛星數量保護閾值nsin=4,終止判定閾值λ=0.01,圖2中紅色點為采用快速選星算法選出的10顆衛星,全部衛星集合的PDOP值為1.341 4,選星后的PDOP值為1.542 1,相比增長了14.96%,但LTS估計的組合數從136降至45,相比降低了66.91%。將4號星和13號星設為故障衛星,從1 s至50 s加入50 m階躍型偽距偏差,從101 s至150 s加入1 m/s的緩變型偽距偏差,采用各RAIM算法后的定位誤差如圖3所示。

由圖3(a)和圖3(b)可以看出,在存在故障的情況下EKF RAIM的定位誤差明顯增大,而基于M估計的REKF RAIM也僅在極個別時刻正確地檢測并剔除故障衛星。對比圖3(c)和圖3(d),基于MM估計的REKF RAIM和本文所提方法的定位誤差在階躍型故障下均未出現明顯變化,在緩變型故障達到20 m后也能完全正確地檢測并剔除故障星,結果表明MM估計相比于M估計具備更好的抵抗雙星故障的能力。

統計4種算法的計算時間分別為0.116 5 s、0.121 8 s、2.891 0 s和0.839 2 s,結果表明本文所提算法縮短了70.97%的計算時間,在保證抗差效果的前提下大幅度提升了MM估計的計算效率。

分別加入-100～100 m偽距偏差,步長為5 m,3種方案的定位誤差均值如圖4所示。

圖3 故障模式下的定位誤差

圖4 不同偽距偏差值下的定位誤差

由圖4(a)和圖4(b)可以看出,相比于EKF RAIM,REKF RAIM的定位誤差得到明顯改善,但當兩個故障星的偽距偏差具有較高的一致性時,M估計的穩健性受到嚴重破壞,對故障的檢測效果較差。而從圖4(c)可以看出,本文方法則不受影響,對故障星具備更好的檢測及識別效果。

分別向全部雙星故障組合加入50 m的偽距偏差,3種方法的定位誤差均值如圖5所示,REKF RAIM對6號星和7號星、6號星和17號星的故障組合的檢測識別效果也較差,而本文方法對全部故障組合均能正確地檢測與識別。

分別加入20～70 m的偽距偏差,3種方法對全部故障星組合的正確檢測與識別率如圖6所示,當偽距偏差達到45 m后,本文方法的故障檢測與識別率能達到100%,而其他兩種方法均存在漏檢或誤檢的情況。

圖5 不同故障組合下的定位誤差

圖6 不同偽距偏差下的故障檢測與識別率

在加入30 m、50 m、80 m和100 m的偽距偏差時,3種方法全部故障星組合的定位誤差均值統計如表1所示,本文方法對全部故障組合均有很好的抗差效果,相比于基于M估計的REKF RAIM,定位誤差分別降低了10.807%、12.122%、18.076%和21.718%。

表1 不同偽距偏差值下的定位誤差

4 結 論

針對雙星故障模式,本文對基于M估計的REKF RAIM進行了改進,采用崩潰污染率高的LTS估計獲得M估計的尺度參數和迭代初值,并設計了基于特征斜率的快速選星算法優化LTS估計的計算量。仿真結果表明,在雙星故障矢量具有較高的空間一致性時,M估計的穩健性受到嚴重破壞,本文方法則能有效地解決此問題,提高對雙星故障的正確檢測與識別率,降低定位誤差。