基于殘差統(tǒng)計(jì)特性的中國區(qū)域格網(wǎng)電離層GIVE算法優(yōu)化

馬岳鑫，唐成盼，胡小工，常志巧，蒲俊宇，邢 楠，曹月玲，王寧波

1. 中國科學(xué)院上海天文臺，上海 200235； 2. 32021部隊(duì)，北京 100000； 3. 北京師范大學(xué)，北京 100036； 4. 中國科學(xué)院空天信息創(chuàng)新研究院，北京 100036

星基增強(qiáng)系統(tǒng)(SBAS)通過GEO衛(wèi)星播發(fā)GNSS導(dǎo)航衛(wèi)星軌道、鐘差和格網(wǎng)電離層改正數(shù)以及UDRE和GIVE等完好性參數(shù)，用以提升區(qū)域用戶GNSS導(dǎo)航服務(wù)精度與可靠性。1992年美國聯(lián)邦航空管理局(FAA)提出了廣域差分(WAAS)構(gòu)想后，各國開始建設(shè)自己的星基增強(qiáng)系統(tǒng)。目前主要的星基增強(qiáng)系統(tǒng)有美國的WAAS、歐盟的EGNOS、日本的MSAS和印度的GAGAN。美國航空無線電技術(shù)委員會(Radio Technical Commission for Aeronautics,RTCA)2016年頒布的《全球定位系統(tǒng)/廣域增強(qiáng)系統(tǒng)機(jī)載設(shè)備最低運(yùn)行性能標(biāo)準(zhǔn)》(DO-229E)[1]詳細(xì)描述了WAAS報(bào)文數(shù)據(jù)格式與內(nèi)容。RTCA協(xié)議針對GPS L1 C/A頻點(diǎn)單頻用戶進(jìn)行增強(qiáng)。

北斗二號采用區(qū)域網(wǎng)定軌確定衛(wèi)星軌道，基于星地雙向時間同步確定衛(wèi)星鐘差[2-5]，利用北斗Klobuchar 8參數(shù)廣播電離層模型修正電離層延遲以提供基本導(dǎo)航服務(wù)，北斗二號SBAS系統(tǒng)播發(fā)等效鐘差改正數(shù)用于修正衛(wèi)星鐘差與軌道徑向誤差[6-7]，并以5°×2.5°空間分辨率播發(fā)格網(wǎng)電離層修正參數(shù)[8]。北斗三號采用區(qū)域網(wǎng)加星間鏈路定軌確定衛(wèi)星軌道，利用星地星間雙向時間同步確定衛(wèi)星鐘差[9-11]，采用北斗BDGIM 9參數(shù)廣播電離層模型修正電離層延遲以提供基本導(dǎo)航服務(wù)[12]。2018年12月27日《北斗衛(wèi)星導(dǎo)航系統(tǒng)公開服務(wù)性能規(guī)范》[13]發(fā)布，北斗三號全球系統(tǒng)開始提供基本導(dǎo)航服務(wù)[14]，隨著北斗三號全球系統(tǒng)的逐步建成與完善，BDSBAS星基增強(qiáng)系統(tǒng)將成為日后工作重心。BDSBAS需要遵照RTCA協(xié)議增強(qiáng)GPS L1 C/A用戶，同時需要增強(qiáng)BDS B1C用戶。

電離層延遲是影響GNSS單頻用戶服務(wù)精度的主要誤差源之一，除各GNSS系統(tǒng)播發(fā)的廣播電離層參數(shù)外，星基增強(qiáng)系統(tǒng)也采用薄殼電離層模型提供更高精度的電離層延遲誤差改正[15-16]。電離層延遲可通過雙頻偽距或精密單點(diǎn)定位提取[17]。WAAS播發(fā)電離層格網(wǎng)報(bào)文以修正用戶電離層延遲[1，18]。WAAS利用基于梯度變化的克里金插值(KT)算法計(jì)算格網(wǎng)點(diǎn)垂直延遲，并基于卡方因子算法實(shí)現(xiàn)電離層異常狀態(tài)的檢驗(yàn)及GIVE參數(shù)的計(jì)算[19-20]。相較于平面擬合算法，基于克里金插值算法處理得到的格網(wǎng)電離層改正精度可提升3%～6%[21]。此后，基于最差用戶幾何位置的格網(wǎng)點(diǎn)延遲方差估值算法解決了由電離層風(fēng)暴引起區(qū)域電離層梯度劇烈變化導(dǎo)致的服務(wù)完好性變差問題[22]；基于電離層層析技術(shù)的克里金插值后處理算法用以測量美國與巴西上空電子密度總含量[23]。印度GAGAN采用平面擬合算法計(jì)算格網(wǎng)點(diǎn)垂直延遲，采用與WAAS一致的算法計(jì)算GIVE，并且給出了適合赤道地區(qū)電離層活動的先驗(yàn)參考分布[24]。基于印度上空實(shí)測數(shù)據(jù)表明，相較于平面擬合算法，采用克里金插值計(jì)算格網(wǎng)點(diǎn)垂直延遲更為穩(wěn)定[25]。適合中國區(qū)域的電離層球殼高度和格網(wǎng)電離層模型在后續(xù)研究中得到了詳細(xì)論證[26-28]。陸態(tài)網(wǎng)實(shí)測數(shù)據(jù)表明，中國區(qū)域格網(wǎng)電離層采用平面擬合算法估算格網(wǎng)點(diǎn)垂直延遲RMS略低于克里金插值，但存在嚴(yán)重的邊際效應(yīng)[29]。

本文參考WAAS格網(wǎng)電離層計(jì)算算法，綜合平面擬合與反比距離加權(quán)方法實(shí)現(xiàn)BDSBAS格網(wǎng)電離層格網(wǎng)點(diǎn)垂直延遲估值的計(jì)算。為避免粗差導(dǎo)致平面擬合估計(jì)值偏離實(shí)際值，根據(jù)克里金插值平穩(wěn)假設(shè)，電離層活動正常狀態(tài)下，利用小范圍(5°×5°)內(nèi)的穿刺點(diǎn)中值估計(jì)樣本均值，利用樣本方差序列中值估計(jì)樣本方差剔除粗差?？紤]到由電離層活動不規(guī)則導(dǎo)致的殘差分布的非正態(tài)性，GIVE的計(jì)算考慮了殘差分布的峰度系數(shù)與偏度系數(shù)，最終結(jié)合基于歷史信息的卡方因子，能夠?qū)崿F(xiàn)99.9%以上的GIVE包絡(luò)率并有效避免過包絡(luò)。本文首先介紹了BDSBAS格網(wǎng)電離層容錯算法與格網(wǎng)改正數(shù)解算算法；接著重點(diǎn)介紹了基于殘差分布偏度與峰度系數(shù)的GIVE算法；最終用2020年1月BDSBAS監(jiān)測接收機(jī)實(shí)測數(shù)據(jù)驗(yàn)證了算法的有效性并簡單評估了格網(wǎng)電離層服務(wù)精度。

1 粗差剔除與殘差統(tǒng)計(jì)

采用監(jiān)測接收機(jī)雙頻相位平滑偽距觀測數(shù)據(jù)提取電離層延遲，不可避免地存在粗差。粗差的存在會導(dǎo)致格網(wǎng)點(diǎn)垂直延遲估值失真，且顯著增加平差迭代次數(shù)，導(dǎo)致格網(wǎng)電離層解算耗時增加。

采用上述粗差剔除算法，可顯著降低迭代次數(shù)并提高格網(wǎng)電離層格網(wǎng)點(diǎn)垂直延遲估值的精度與可靠性。圖1給出了采用中值容錯前后的格網(wǎng)電離層殘差分布：圖1(a)為未采用中值容錯，僅采用多次迭代方式剔除粗差的格網(wǎng)電離層殘差頻率分布直方圖；圖1(b)為采用中值容錯后的格網(wǎng)電離層殘差頻率分布直方圖；圖1(c)為對應(yīng)殘差分布的Q-Q圖，圖中十字點(diǎn)線越靠近參考虛線，則表明樣本分布為正態(tài)分布的可能性越大；圖1(d)為對應(yīng)殘差分布的Q-Q圖。

圖1 采用中值容錯前后的格網(wǎng)電離層殘差分布Fig.1 Grid ionospheric residual distribution of median fault tolerance

由圖1可知，采用中值容錯剔除粗差后，Q-Q圖十字點(diǎn)線更靠近參考虛線，樣本分布偏度與峰度值更接近0；格網(wǎng)電離層殘差分布會更加靠近正態(tài)分布。這說明了粗差剔除合理有效。

2 格網(wǎng)點(diǎn)延遲估計(jì)

(1)

(2)

同一衛(wèi)星或者同一監(jiān)測接收機(jī)的觀測量間存在相關(guān)性，故令協(xié)方差矩陣W-1為

(3)

(4)

GT·W·Iv,IPP]

(5)

式中，Iv,IPP=[Iv,IPP1Iv,IPP2…Iv,IPPn]T為格網(wǎng)點(diǎn)附近的穿刺點(diǎn)(IPP)處電離層垂直延遲。

(6)

3 GIVE的計(jì)算

3.1 偏度與峰度

偏度(skewness)[20]反映了殘差分布的非對稱性，正態(tài)分布偏度系數(shù)為0。采樣于未知分布的樣本偏度值絕對值越大，表明該未知分布是正態(tài)分布的可能性越小。

偏度(skewness)α的計(jì)算方法為

(7)

式中，μ為未知分布期望，可用樣本均值近似；ki為未知分布的i階矩,可用樣本i階矩近似，計(jì)算公式為

(8)

式中，xj為采樣于未知分布的樣本；N為樣本總數(shù)。

峰度(kurtosis)[30]值β反映了殘差分布的拖尾性質(zhì)，正態(tài)分布峰度系數(shù)為0。采樣于未知分布的樣本峰度值絕對值越大，表明該未知分布是正態(tài)分布的可能性越小，當(dāng)峰度(kurtosis)值β小于0時，表示樣本分布為厚尾分布。

峰度(kurtosis)β的計(jì)算方法為

(9)

采用平面擬合得到格網(wǎng)點(diǎn)垂直延遲后，根據(jù)RTCA協(xié)議[1]提供的格網(wǎng)電離層用戶算法插值計(jì)算穿刺點(diǎn)垂直延遲，與雙頻提取電離層垂直延遲做差得到殘差序列，根據(jù)殘差序列的統(tǒng)計(jì)特性確定完好性參數(shù)GIVE。

圖2為2019年11月23日—2019年11月25日雙頻提取與格網(wǎng)計(jì)算電離層延遲互差的頻率分布直方圖，期間未發(fā)生電離層暴，計(jì)算殘差分布偏度值為0.116、峰度值為0.585。

圖2 未發(fā)生電離層暴時電離層格網(wǎng)殘差分布的偏度系數(shù)與峰度系數(shù)Fig.2 Skewness coefficient and kurtosis coefficient of the ionospheric grid residual distribution under ideal conditions

統(tǒng)計(jì)2019年10月20日—2020年3月18日的格網(wǎng)電離層殘差分布偏度值與峰度值得到表1。從表中可以看出，99.9%的格網(wǎng)點(diǎn)殘差分布偏度值絕對值在2.507以下、峰度值絕對值在9.459以下。當(dāng)殘差分布偏度值與峰度值偏離0時，說明殘差分布為正態(tài)分布的可能性變小，基于正態(tài)分布假設(shè)給出的0.999概率置信區(qū)間已經(jīng)不在可靠，需要擴(kuò)大置信區(qū)間。

表1 格網(wǎng)電離層殘差分布峰度與偏度值

3.2 邊界方差的估計(jì)

為反映電離層格網(wǎng)點(diǎn)垂直延遲估值的可靠性，為用戶提供最差包絡(luò)，必須同時考慮殘差分布特性與電離層歷史信息，WAAS采用卡方因子來實(shí)現(xiàn)電離層歷史信息的引入。格網(wǎng)點(diǎn)邊界方差由式(10)[20]計(jì)算

(10)

(11)

卡方統(tǒng)計(jì)量計(jì)算公式為

(12)

去相關(guān)參數(shù)δdecorr被解釋為平面擬合算法模型表達(dá)誤差，采用去相關(guān)參數(shù)δdecorr可以有效提高包絡(luò)率，但可能導(dǎo)致過包絡(luò)問題。由于地磁赤道穿過中國南部地區(qū)，電離層梯度變化更為復(fù)雜，當(dāng)采用與WAAS相同的去相關(guān)參數(shù)時，服務(wù)區(qū)南部地區(qū)部分格網(wǎng)點(diǎn)包絡(luò)率低于99.9%，若擴(kuò)大去相關(guān)參數(shù)會導(dǎo)致服務(wù)區(qū)內(nèi)北部地區(qū)部分格網(wǎng)點(diǎn)過包絡(luò)，BDSBAS格網(wǎng)電離層不宜采用固定的去相關(guān)參數(shù)δdecorr辦法。

(13)

式中

(14)

式中，α和β為利用格網(wǎng)點(diǎn)附近穿刺點(diǎn)殘差計(jì)算的樣本偏度與峰度系數(shù)，根據(jù)表1，α0與β0取其99.9%分位數(shù)。

3.3 GIVE的計(jì)算

(15)

圖3 卡方檢驗(yàn)因子對電離層異常天氣的響應(yīng)情況圖Fig.3 Response of chi-square test factor to abnormal space weather in ionosphere

圖4 卡方檢驗(yàn)因子對格網(wǎng)電離層殘差分布偏離正態(tài)分布時的響應(yīng)情況Fig.4 Response of chi-square test factor to grid ionospheric residual distribution deviation from the normal distribution

格網(wǎng)點(diǎn)完好性參數(shù)GIVE采用式(16)計(jì)算

(16)

4 BDSBAS格網(wǎng)電離層服務(wù)精度評估

截至2020年3月，中國境內(nèi)32個BDSBAS監(jiān)測站數(shù)據(jù)入站率穩(wěn)定，監(jiān)測接收機(jī)IFB穩(wěn)定，能夠滿足區(qū)域格網(wǎng)電離層解算需求。監(jiān)測站均勻分布在中國境內(nèi)，每個監(jiān)測站均配備3臺BDS監(jiān)測接收機(jī)、3臺GNSS多模監(jiān)測接收機(jī)，可以同時解碼獲取GPS、GLONASS、Galileo 3系統(tǒng)觀測數(shù)據(jù)與導(dǎo)航電文。目前僅利用北斗三號B1C-B2A、北斗二號B1I-B3I、GPS L1 C/A-L2P雙頻相位平滑偽距(CNMC)[31]提取電離層延遲參與格網(wǎng)解算。

按照本文所述算法剔除數(shù)據(jù)粗差后，再根據(jù)本文所述算法計(jì)算電離層格網(wǎng)點(diǎn)垂直延遲和GIVE，統(tǒng)計(jì)格網(wǎng)電離層改正殘差RMS、改正百分比和GIVE包絡(luò)率，以驗(yàn)證算法精度與可靠性。同時為考察BDSBAS格網(wǎng)電離層增強(qiáng)服務(wù)精度，利用本文算法計(jì)算得到的格網(wǎng)電離層產(chǎn)品，選取服務(wù)范圍內(nèi)的4個監(jiān)測站偽距單點(diǎn)定位，并統(tǒng)計(jì)定位精度提升情況。統(tǒng)計(jì)結(jié)果時間系統(tǒng)均采用北斗時；格網(wǎng)電離層延遲單位為m。

采用2019年12月12日32個監(jiān)測站雙頻提取電離層延遲，截止高度角設(shè)置為15°，利用三角投影函數(shù)[32]計(jì)算穿刺點(diǎn)垂直延遲。利用上文所述算法計(jì)算格網(wǎng)點(diǎn)垂直延遲，計(jì)算GIVE時不考慮去相關(guān)參數(shù)，統(tǒng)計(jì)GIVE一天包絡(luò)率并繪制圖5。

圖5 2019年12月12日BDSBAS格網(wǎng)電離層未增加因子和去相關(guān)參數(shù)的GIVE包絡(luò)情況Fig.5 GIVE envelope diagram of the BDSBAS grid ionosphere without factor and decorrelation parameters on December 12, 2019

圖6 2019年12月12日BDSBAS格網(wǎng)電離層GIVE未包絡(luò)部分殘差分布Fig.6 Distribution of residuals in the unenveloped part of the ionospheric GIVE on December 12, 2019

圖6為GIVE未包絡(luò)部分殘差分布圖。由圖6可知，不考慮去相關(guān)參數(shù)和殘差偏度系數(shù)與峰度系數(shù)時，GIVE包絡(luò)存在大量的漏警點(diǎn)；將漏警率高的時段單獨(dú)取出，繪制該時段殘差分布圖，得到圖6(b)和圖6(c)，可以看出，此時段殘差分布峰度系數(shù)與偏度系數(shù)均明顯偏離0。

仿照WAAS的做法，固定的去相關(guān)參數(shù)為0.35 m，計(jì)算GIVE并統(tǒng)計(jì)一天包絡(luò)率得到圖7。圖7中實(shí)線為GIVE包絡(luò)線，散點(diǎn)為格網(wǎng)電離層殘差序列點(diǎn)。圖7(a)為30°N—55°N范圍內(nèi)的穿刺點(diǎn)包絡(luò)圖，圖7(b)為0°N—30°N范圍內(nèi)的穿刺點(diǎn)包絡(luò)圖。由圖7可以看出，若采用固定的去相關(guān)參數(shù)，對于中國區(qū)域0°N—30°N范圍內(nèi)的穿刺點(diǎn)GIVE過包絡(luò)現(xiàn)象明顯；對于中國區(qū)域低緯度地區(qū)仍存在不少的漏警點(diǎn)，包絡(luò)率僅為99.8%左右。若擴(kuò)大去相關(guān)參數(shù)，可以進(jìn)一步減少漏警率并提高GIVE包絡(luò)率，但難免會引入過包絡(luò)的問題，使得GIVE估值較實(shí)際偏大。

圖7 2019年12月12日BDSBAS格網(wǎng)電離層采用固定的去相關(guān)參數(shù)的GIVE包絡(luò)情況Fig.7 GIVE envelope diagram of the BDSBAS ionosphere grid with fixed decorrelation parameters on December 12,2019

2020年1月20日發(fā)生一次小型電離層暴(http:∥www.nsmc.org.cn/NSMC)，圖9反映了電離層異常天氣時的包絡(luò)性能。圖9(a)采用伯爾尼大學(xué)(AIUB)事后處理獲取的全球電離層延遲圖(GIM)繪制，從圖中可以看出當(dāng)日電離層暴導(dǎo)致了GIM異常雙峰現(xiàn)象；圖9(b)為當(dāng)日境內(nèi)BDSABS監(jiān)測站所有穿刺點(diǎn)格網(wǎng)電離層殘差頻率分布直方圖，其偏度值與峰度值靠近公式(14)給出的閾值；圖9(c)為當(dāng)天GIVE包絡(luò)圖。從圖9可以看出，殘差分布的偏度值峰度值會及時響應(yīng)電離層暴等異?？臻g天氣，放大GIVE并保證包絡(luò)率仍在99.9%以上。

圖8 2019年12月12日BDSBAS格網(wǎng)電離層增加因子后GIVE包絡(luò)情況Fig.8 GIVE envelope diagram of the BDSBAS ionosphere grid based on December 12, 2019

圖9 2020年1月20日發(fā)生電離層暴時GIVE包絡(luò)情況Fig.9 The GIVE envelope diagram during Ionospheric storm on January 20,2020

為進(jìn)一步評價本文提出的GIVE算法包絡(luò)性和電離層格網(wǎng)點(diǎn)垂直延遲解算精度，利用2020年1月1日—2020年1月31日連續(xù)一個月32個BDSBAS境內(nèi)監(jiān)測站實(shí)測數(shù)據(jù)，按照本文算法計(jì)算格網(wǎng)點(diǎn)垂直延遲和GIVE，從格網(wǎng)電離層殘差RMS(分別與雙頻實(shí)測數(shù)據(jù)和CODE GIM產(chǎn)品對比)、改正百分比、GIVE包絡(luò)率3個角度評價BDSBAS格網(wǎng)電離層服務(wù)精度。選取華中、華北、華東、華南、西北、東北、西南地區(qū)7個格網(wǎng)點(diǎn)，統(tǒng)計(jì)2020年1月整月RMS、改正百分比和GIVE包絡(luò)率得到表2。

表2 BDSBAS格網(wǎng)電離層服務(wù)精度統(tǒng)計(jì)評價

為簡單評估BDSBAS格網(wǎng)電離層服務(wù)精度，利用2020年1月1日—2020年1月31日連續(xù)一個月32個BDSBAS境內(nèi)監(jiān)測站實(shí)測數(shù)據(jù)，按照本文算法計(jì)算格網(wǎng)點(diǎn)垂直延遲和GIVE，選取位于華中(Sta1)、華北(Sta2)、華東(Sta3)、華南(Sta4)的4個監(jiān)測站統(tǒng)計(jì)增強(qiáng)定位提升百分比得到表3。偽距單點(diǎn)定位采用本團(tuán)隊(duì)編寫的性能監(jiān)視程序，其雙頻相位平滑偽距基本導(dǎo)航定位精度在5 m(95%)以內(nèi)，將事后精密單點(diǎn)定位(PPP)解算坐標(biāo)作為參考真值考察增強(qiáng)定位提升百分比。

表3 BDSBAS格網(wǎng)電離層服務(wù)定位精度統(tǒng)計(jì)評價

由表2和表3可知，BDSBAS格網(wǎng)電離層改正殘差RMS平均在0.4 m左右，大約在2～3 TECU量級；改正百分比均在75%以上，平均在78%左右。利用本文提出的GIVE算法，GIVE包絡(luò)率均在99.9%以上，由于考慮了殘差分布的偏度特性與峰度特性，GIVE包絡(luò)在降低漏警率的同時抑制了過包絡(luò)現(xiàn)象。利用BDSBAS監(jiān)測接收機(jī)進(jìn)行偽距單點(diǎn)定位，在格網(wǎng)電離層增強(qiáng)的狀態(tài)下，對于BDS BIC頻點(diǎn)用戶，定位精度增強(qiáng)不明顯；對于GPS L1 C/A頻點(diǎn)用戶，可增強(qiáng)定位精度20%～40%。GPS L1 C/A頻點(diǎn)用戶增強(qiáng)效果明顯，是由于GPS基本導(dǎo)航用戶采用了Klobuchar 八參數(shù)廣播電離層模型，該模型改正百分比大約在50%左右；而BDS BIC頻點(diǎn)用戶采用BDGIM廣播電離層模型修正電離層延遲，該模型改正百分比全球平均75%[27]，中國境內(nèi)改正百分比更高，故定位精度提升不明顯?？紤]到格網(wǎng)電離層電文更新頻率比廣播電離層模型高，相比于精度提升，電離層異?；顒蛹皶r預(yù)警對于BDS BIC頻點(diǎn)用戶更為重要。

5 結(jié) 論

本文提出了基于殘差分布峰度系數(shù)與偏度系數(shù)的格網(wǎng)電離層GIVE算法。相比于采用固定的去相關(guān)參數(shù)計(jì)算GIVE，基于殘差分布偏度值與峰度值自適應(yīng)調(diào)整去相關(guān)參數(shù)計(jì)算的GIVE可以實(shí)現(xiàn)在降低漏警率的同時有效抑制過包絡(luò)現(xiàn)象。新算法計(jì)算的完好性參數(shù)GIVE可以做到對服務(wù)區(qū)內(nèi)任何緯度范圍內(nèi)的格網(wǎng)點(diǎn)實(shí)現(xiàn)99.9%以上的包絡(luò)率。除此之外，本文同時介紹了BDSBAS容錯算法與格網(wǎng)改正數(shù)計(jì)算算法?；贐DSBAS實(shí)測數(shù)據(jù)的統(tǒng)計(jì)表明，格網(wǎng)電離層修正RMS約在2～3 TECU，改正百分比達(dá)到75%～79%；修正格網(wǎng)電離層后可提升GPS定位精度20%～40%。