王文博, 徐 穎

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

0 引 言

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

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

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

1 REKF RAIM方法

1.1 REKF RAIM模型

線(xiàn)性化后接收機(jī)的觀(guān)測(cè)方程為

Yk=HkXk+εk

(1)

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

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

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

(2)

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

REKF估計(jì)器包括預(yù)測(cè)和修正兩個(gè)步驟。

預(yù)測(cè):

(3)

(4)

(5)

修正:

(6)

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

(7)

新息為觀(guān)測(cè)值與預(yù)測(cè)觀(guān)測(cè)值之差,表示為

(8)

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

rk～N(0,Dk)

(9)

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

(10)

根據(jù)新息構(gòu)造檢驗(yàn)統(tǒng)計(jì)量

(11)

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

根據(jù)給定的虛警概率Pfa計(jì)算檢驗(yàn)門(mén)限TD,即

(12)

1.2 M估計(jì)及IGG III等價(jià)權(quán)函數(shù)

當(dāng)觀(guān)測(cè)信息中存在異常值時(shí),擴(kuò)展卡爾曼濾波估計(jì)值可能會(huì)失真甚至不收斂,因此引入抗差估計(jì)理論來(lái)抑制異常值的影響是有必要的。REKF RAIM方法的流程如圖1所示。

圖1 REKF RAIM算法流程圖

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

(13)

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

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

(14)

通過(guò)設(shè)計(jì)合適的權(quán)函數(shù)ρ來(lái)抑制異常觀(guān)測(cè)值的影響,常用的等價(jià)權(quán)函數(shù)有丹麥法、Huber法和IGG法等,本文采用三段式的IGG III等價(jià)權(quán)函數(shù)

(15)

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

2 兩步抗差MM估計(jì)

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

2.1 LT S估計(jì)

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

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

(16)

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

LTS估計(jì)的步驟如下:

步驟 3取殘差平方和最小的方程組的解作為L(zhǎng)TS的估計(jì)解。

h≤n-q

(17)

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

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

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

(18)

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

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

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

3 仿真驗(yàn)證

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

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

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

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

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

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

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

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

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

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

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

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

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

圖6 不同偽距偏差下的故障檢測(cè)與識(shí)別率

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

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

4 結(jié) 論

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