王文博, 徐 穎

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

0 引 言

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

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

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

1 REKF RAIM方法

1.1 REKF RAIM模型

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

Yk=HkXk+εk

(1)

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

系統(tǒng)的狀態(tài)轉(zhuǎn)移方程為

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

(2)

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

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

預測:

(3)

(4)

(5)

修正:

(6)

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

(7)

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

(8)

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

rk～N(0,Dk)

(9)

式中,Dk為新息的協(xié)方差矩陣

(10)

根據(jù)新息構造檢驗統(tǒng)計量

(11)

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

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

(12)

1.2 M估計及IGG III等價權函數(shù)

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

圖1 REKF RAIM算法流程圖

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

(13)

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

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

(14)

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

(15)

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

2 兩步抗差MM估計

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

2.1 LT S估計

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

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

(16)

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

LTS估計的步驟如下:

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

h≤n-q

(17)

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

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

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

(18)

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

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

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

3 仿真驗證

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

圖2 初始時刻衛(wèi)星分布

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

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

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

統(tǒng)計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估計的穩(wěn)健性受到嚴重破壞,對故障的檢測效果較差。而從圖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種方法全部故障星組合的定位誤差均值統(tǒng)計如表1所示,本文方法對全部故障組合均有很好的抗差效果,相比于基于M估計的REKF RAIM,定位誤差分別降低了10.807%、12.122%、18.076%和21.718%。

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

4 結 論

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