一種雙星故障條件下RAIM可用性評(píng)估的改進(jìn)方法

馬下平 余科根 賀小星 嚴(yán) 麗 趙立都

1 西安科技大學(xué)測(cè)繪科學(xué)與技術(shù)學(xué)院，西安市雁塔路58號(hào)，710054 2 中國礦業(yè)大學(xué)環(huán)境與測(cè)繪學(xué)院，江蘇省徐州市大學(xué)路1號(hào)，221116 3 江西理工大學(xué)土木與測(cè)繪工程學(xué)院，江西省贛州市紅旗大道86號(hào)，341000 4 江西省防震減災(zāi)與工程地質(zhì)災(zāi)害探測(cè)工程研究中心，南昌市廣蘭大道418號(hào)，330013 5 重慶交通大學(xué)土木工程學(xué)院，重慶市學(xué)府大道66號(hào)，400074

BDS是我國自主建設(shè)、獨(dú)立運(yùn)行的全球衛(wèi)星導(dǎo)航系統(tǒng)[1-4]，可為用戶提供更多高定位精度和可靠性的導(dǎo)航信息，同時(shí)也可為民用航空接收機(jī)自主完好性監(jiān)測(cè)(receiver autonomous integrity monitoring，RAIM)提供有力支撐。RAIM對(duì)衛(wèi)星進(jìn)行故障監(jiān)測(cè)會(huì)受到衛(wèi)星個(gè)數(shù)和衛(wèi)星幾何分布的影響。當(dāng)可見衛(wèi)星數(shù)過少或幾何分布不理想時(shí)，由衛(wèi)星幾何結(jié)構(gòu)引起的定位誤差會(huì)遮掩衛(wèi)星故障引起的定位誤差，造成完好性監(jiān)測(cè)結(jié)果不可信，漏檢率增加。因此在進(jìn)行故障監(jiān)測(cè)之前必須進(jìn)行可用性判斷，以保證不會(huì)影響故障監(jiān)測(cè)性能[5]。RAIM可用性計(jì)算中涉及到漏檢率和誤警率參數(shù)，王爾申等[6]研究發(fā)現(xiàn)，RAIM可用性隨漏檢率和誤警率的降低而降低。

目前，基于單星故障的RAIM可用性評(píng)估方法主要有水平保護(hù)法(horizontal protection level method，HPLM)、最大水平精度因子法(δHDOPmax)和近似徑向誤差保護(hù)法(approximate radial-error protection method，ARPM)，理論上這3種方法具有等價(jià)性[7-8]；考慮到2顆衛(wèi)星同時(shí)發(fā)生故障的可能性，倪育德等[9]從1顆衛(wèi)星故障的可用性出發(fā)，推導(dǎo)出2顆衛(wèi)星故障的可用性；針對(duì)多故障衛(wèi)星的可用性，陳金平等[10]提出基于漏檢概率和圓概率誤差的可用性分析方法。隨著GNSS的發(fā)展與應(yīng)用以及飛機(jī)更高階段精密進(jìn)近的需求，具備垂直導(dǎo)航能力的高級(jí)接收機(jī)自主完好性監(jiān)測(cè)(advanced RAIM，ARAIM)被引入。ARAIM旨在全球范圍內(nèi)提供垂直導(dǎo)向-200 ft定標(biāo)性能(LPV-200)級(jí)別的航空導(dǎo)航服務(wù)，眾多學(xué)者對(duì)此進(jìn)行了相關(guān)算法改進(jìn)和定位性能方面研究[11-17]。

現(xiàn)有的雙星故障條件下RAIM可用性評(píng)估大多沿用已有的矩陣最大特征值方法(matrix maximum eigenvalue method，MMEM)[18]，該方法可解決2個(gè)粗差比值的計(jì)算，但該值還可采用高等數(shù)學(xué)中求極值方法得到。為此，本文提出針對(duì)RAIM可用性評(píng)估的極大值方法(maxima method，MM)，并基于完好性風(fēng)險(xiǎn)參數(shù)和IGMAS中8個(gè)GNSS站的實(shí)測(cè)數(shù)據(jù)對(duì)這2種方法的RAIM可用性進(jìn)行驗(yàn)證，得到最新的中國境內(nèi)BDS可用性性能。

1 RAIM方法原理

1.1 RAIM觀測(cè)方程建立

RAIM主要是基于GNSS中的偽距觀測(cè)量。若在歷元t時(shí)接收機(jī)r與衛(wèi)星s進(jìn)行同步觀測(cè)，并觀測(cè)m顆衛(wèi)星，則觀測(cè)方程可表示為：

y=Hx+ε

(1)

式中，y為m維偽距觀測(cè)值與近似衛(wèi)地距的差向量；H為m行4列設(shè)計(jì)矩陣，表示各衛(wèi)星對(duì)用戶的投影向量；x為測(cè)站坐標(biāo)和接收機(jī)鐘差構(gòu)成的四維向量，ε為m維觀測(cè)偽距的誤差向量。

給定觀測(cè)量權(quán)陣為P，則式(1)的最小二乘解為：

=x+(HTPH)-1HTPε

(2)

=(I-H(HTPH)-1HTP)ε

(3)

令

(4)

(5)

1.2 檢驗(yàn)統(tǒng)計(jì)量構(gòu)建

(6)

(7)

一般取檢驗(yàn)統(tǒng)計(jì)量(T)為：

(8)

在導(dǎo)航過程中，將實(shí)時(shí)解算的T值與監(jiān)測(cè)門限(TD)進(jìn)行比較，TD計(jì)算可參考文獻(xiàn)[19]。若T

2 RAIM可用性分析方法

在進(jìn)行RAIM之前，首先根據(jù)性能指標(biāo)對(duì)當(dāng)前可見星的幾何分布進(jìn)行判斷，分析是否適合進(jìn)行完好性監(jiān)測(cè)，即判斷RAIM算法是否可用。

2.1 監(jiān)測(cè)統(tǒng)計(jì)量與定位誤差的關(guān)系

RAIM可用性判斷時(shí)需要構(gòu)建T與定位誤差之間的數(shù)學(xué)關(guān)系。由于RAIM不具備垂直導(dǎo)航能力，因此本文只討論水平保護(hù)級(jí)別(horizontal protection level，HPL)。

圖1中橫軸表示T，縱軸表示水平徑向誤差(horizontal radial error，HRE)，將這2個(gè)變量取比值后得到圖中斜線斜率(SLOPE)，即

(9)

圖1 T與HPL之間的關(guān)系Fig.1 The relationship between T and HRE

從圖1可以看出，T與HRE呈線性關(guān)系，沿此傾斜線，當(dāng)偽距觀測(cè)量中存在偏差時(shí)，T與HRE線性增加。偏差較小時(shí)，HRE小于HPL且T小于TD，系統(tǒng)處于正常狀態(tài)(圖1左下角)，PMD較低；偏差較大時(shí)，HRE大于HPL且T大于TD，系統(tǒng)能正確監(jiān)測(cè)偏差(圖1右上角)，PMD也較低；當(dāng)偏差處于兩者之間時(shí)，PMD可能達(dá)到最大(圖1左上角)。對(duì)應(yīng)相同T值時(shí)，SLOPE越大，衛(wèi)星發(fā)生故障時(shí)MD概率越高。對(duì)應(yīng)衛(wèi)星的SLOPEs為：

(10)

在m個(gè)觀測(cè)量中，HPL可表示為：

HPL=max(SLOPEs)·T=SLOPEmax·T

(11)

RAIM可用性分析就是利用SLOPEs在各歷元的最大值SLOPEmax進(jìn)行判定，依據(jù)SLOPEmax求出HPL，并與水平告警限值(horizontal alarm limit，HAL)進(jìn)行比較，如果HPL

2.2 雙星故障時(shí)MM求解HPL

在忽略觀測(cè)偽距隨機(jī)誤差的情況下，當(dāng)觀測(cè)誤差ε在第i個(gè)和第j個(gè)觀測(cè)值中存在粗差bi和bj時(shí)，則有：

(12)

根據(jù)式(5)可以得到定位解的誤差向量為：

(13)

則有：

(14)

將式(12)代入式(5)可得：

(15)

(16)

以往對(duì)k值主要采用矩陣最大特征值方法(MMEM)[9,18,20]進(jìn)行求解，但其也可采用高等數(shù)學(xué)中求極值的方法獲取，具體如下：

為求式(16)的極大值，兩邊取對(duì)數(shù)，并對(duì)式中k求偏導(dǎo)，可得一元二次方程式：

ak2+bk+c=0

(17)

設(shè)k1和k2分別為式(17)兩根，將k1、k2分別代入式(16)，則當(dāng)前歷元的(SLOPEij)max可表示為：

(SLOPEij)max=

max((SLOPEij)k1,(SLOPEij)k2)

(18)

即可得到雙星故障時(shí)的HPL為[19]：

(19)

2.3 雙星故障時(shí)MMEM求解HPL

目前求解HPL主要采用MMEM方法，該方法主要是構(gòu)建關(guān)于粗差向量b的二次型矩陣，根據(jù)矩陣的最大值為矩陣最大特征值這一定律求解HPL，具體計(jì)算過程可參考文獻(xiàn)[18]。

3 數(shù)據(jù)處理

3.1 實(shí)驗(yàn)數(shù)據(jù)和處理策略

實(shí)驗(yàn)數(shù)據(jù)采用國際GNSS監(jiān)測(cè)評(píng)估系統(tǒng)(international GNSS monitoring and assessment system，IGMAS)提供的中國境內(nèi)BJF1、CHU1、GUA1、KUN1、LHA1、SHA1、WUH1、XIA1共8個(gè)連續(xù)跟蹤站2020-09-06的觀測(cè)數(shù)據(jù)和星歷文件。通過自編程序統(tǒng)計(jì)這些測(cè)站在截止高度角5°以上的BDS衛(wèi)星可見性、DOP值和基于2種方法的雙星故障HPL值，比較并分析中國境內(nèi)BDS可用性和2種雙星故障RAIM可用性的差異。具體數(shù)據(jù)情況及處理策略為：1)BJF1和CHU1站觀測(cè)時(shí)段為00:00:00～09:59:30；GUA1、KUN1、LHA1、SHA1站觀測(cè)時(shí)段為00:00:00～15:59:30；WUH1和XIA1站觀測(cè)時(shí)段為00:00:00～16:59:30；2)考慮到每個(gè)測(cè)站不完全包含雙頻或三頻觀測(cè)值，因此采用BDS B1I頻率的單頻觀測(cè)值；3)全球星歷文件中每隔1 h給出用于計(jì)算每個(gè)歷元的BDS衛(wèi)星軌道參數(shù)、衛(wèi)星鐘差參數(shù)和衛(wèi)星狀態(tài)信息，其1 d的全球星歷文件可從MGEX網(wǎng)站下載；4)電離層延遲采用Klobuchar模型進(jìn)行改正，對(duì)流層延遲采用Saastamoinen模型進(jìn)行改正；5)式(2)中P采用已知的高度角定權(quán)法[21]，定權(quán)時(shí)考慮3種異構(gòu)衛(wèi)星之間的差異，將GEO和IGSO+MEO衛(wèi)星的權(quán)比關(guān)系設(shè)置為1∶5[22]。

3.2 衛(wèi)星可見數(shù)和DOP值

利用8個(gè)IGMAS站的BDS觀測(cè)數(shù)據(jù)和全球星歷文件，計(jì)算出每個(gè)測(cè)站每個(gè)歷元對(duì)應(yīng)衛(wèi)星的高度角，并統(tǒng)計(jì)可見衛(wèi)星數(shù)、平均衛(wèi)星數(shù)、GDOP值及其平均值。結(jié)果表明：1)中國境內(nèi)8個(gè)IGMAS站的BDS可見衛(wèi)星數(shù)最小值為17顆，完全能夠滿足RAIM可用性評(píng)估；2)GDOP最小值、平均值和最大值的均值分別為1.117、1.437、2.186，說明BDS衛(wèi)星在中國境內(nèi)空間分布情況很好。

3.3 2種方法計(jì)算雙星故障時(shí)HPL比較

目前BDS在中國區(qū)域的定位精度為5m左右[4]。保守起見，將σ0設(shè)置為8 m，PMD設(shè)置為10-3，PFA設(shè)置為10-5[19]。利用2種方法分別計(jì)算BDS雙星故障時(shí)HPL值，圖2～5分別為8個(gè)IGMAS站的HPL坐標(biāo)時(shí)間序列，表1(單位m)

圖2 雙星故障條件下BJF1和CHU1站HPL時(shí)間序列Fig.2 HPL time series of BJF1 and CHU1 stations under double-satellite faults conditions

圖3 雙星故障條件下GUA1和KUN1站HPL時(shí)間序列Fig.3 HPL time series of GUA1 and KUN1 stations under double-satellite faults conditions

圖4 雙星故障條件下LHA1和SHA1站 HPL時(shí)間序列Fig.4 HPL time series of LHA1 and SHA1 stations under double-satellite faults conditions

圖5 雙星故障條件下WUH1和XIA1站HPL時(shí)間序列Fig.5 HPL time series of WUH1 and XIA1 stations under double-satellite faults conditions

表1 2種方法計(jì)算的HPL值統(tǒng)計(jì)

為8個(gè)IGMAS站的HPL最小值、最大值和平均值。通過分析可知：

1)從整體上看，MMEM計(jì)算的HPL小于MM計(jì)算的HPL，表明MMEM得到的RAIM可用性高。利用MM求解的8個(gè)IGMAS站的HPL最小值、平均值、最大值的均值分別為124.921 m、233.141m、3 020.005 m；MMEM求解的HPL最小值、平均值和最大值的均值分別為219.093 m、515.110 m、1 302.518 m。

2)2種方法得到的HPL時(shí)間序列中，大部分歷元的HPL值一致性很好，HPL平均值的差值約為281.969 m，主要原因可能為這2種方法的計(jì)算公式中所用參數(shù)個(gè)數(shù)和精度存在差異。

3)在中國區(qū)域，雙星故障條件下HPL值的空間分布存在差異，表明RAIM可用性與地理位置相關(guān)。

3.4 計(jì)算時(shí)間比較

表2(單位s)為2種方法的計(jì)算時(shí)間，從表中可以看出，MM和MMEM的平均計(jì)算時(shí)間分別為11.733 s和21.387 s，MM比MMEM計(jì)算時(shí)間短。

3.5 2種方法雙星故障時(shí)可用性統(tǒng)計(jì)

飛機(jī)在飛行階段包含非精密進(jìn)近(non-precision approach，NPA)、終端、本土航路和遠(yuǎn)洋航路4個(gè)飛行階段[23]，4個(gè)階段的HAL值分別為555.6 m、1 852 m、5 556 m和7 408 m?；?種方法得到雙星故障條件下HPL值，根據(jù)每個(gè)測(cè)站計(jì)算出觀測(cè)歷元中各階段的可用性百分比，再計(jì)算8個(gè)測(cè)站的可用性比值均值，得到每個(gè)階段的可用性結(jié)果(表3)。從表中可以看出：1)MM計(jì)算的可用性在航路和遠(yuǎn)洋階段達(dá)到100%，而MMEM計(jì)算的RAIM可用性除在NPA階段外，其他階段的可用性均達(dá)到100%；2)結(jié)合所計(jì)算的HPL值，雖然MMEM在NPA階段的可用性不及MM，但整體上MMEM可用性優(yōu)于MM。

表2 2種方法的計(jì)算時(shí)間比較

表3 單星和雙星故障條件下RAIM可用性統(tǒng)計(jì)

4 結(jié) 語

本文從雙星存在故障的前提條件出發(fā)，在總結(jié)已有的最大特征值計(jì)算RAIM可用性方法的基礎(chǔ)上，提出RAIM可用性評(píng)估的極大值方法，并以2020-09-06中國區(qū)域8個(gè)IGMAS連續(xù)跟蹤站的BDS觀測(cè)數(shù)據(jù)為例，比較和分析中國境內(nèi)HPL變化及RAIM可用性情況，得到以下結(jié)論：

1)中國境內(nèi)可見衛(wèi)星數(shù)和GDOP完全能夠滿足可用性計(jì)算需求，中國境內(nèi)BDS目前至少可觀測(cè)到17顆可見衛(wèi)星，GDOP值最大為2.186，可為RAIM可用性計(jì)算和故障監(jiān)測(cè)提供豐富的數(shù)據(jù)源。

2)雖然MM計(jì)算的RAIM可用性在NPA階段高于MMEM，且耗時(shí)較少，但總體來說，MMEM計(jì)算HPL的時(shí)間序列小于MM，即對(duì)應(yīng)的RAIM可用性更高。

3)2種方法計(jì)算的HPL存在差異，主要原因可能為各方法所用參數(shù)個(gè)數(shù)及本身精度存在差異，且RAIM可用性計(jì)算中涉及到完好性風(fēng)險(xiǎn)參數(shù)，不同的參數(shù)初值會(huì)對(duì)RAIM可用性評(píng)估產(chǎn)生不同影響，因此參數(shù)取值應(yīng)引起重視。