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

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

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

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

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

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

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

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

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

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

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

2 格網點延遲估計

(1)

(2)

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

(3)

(4)

GT·W·Iv,IPP]

(5)

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

(6)

3 GIVE的計算

3.1 偏度與峰度

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

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

(7)

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

(8)

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

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

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

(9)

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

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

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

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

表1 格網電離層殘差分布峰度與偏度值

3.2 邊界方差的估計

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

(10)

(11)

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

(12)

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

(13)

式中

(14)

式中，α和β為利用格網點附近穿刺點殘差計算的樣本偏度與峰度系數，根據表1，α0與β0取其99.9%分位數。

3.3 GIVE的計算

(15)

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

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

格網點完好性參數GIVE采用式(16)計算

(16)

4 BDSBAS格網電離層服務精度評估

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

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

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

圖5 2019年12月12日BDSBAS格網電離層未增加因子和去相關參數的GIVE包絡情況Fig.5 GIVE envelope diagram of the BDSBAS grid ionosphere without factor and decorrelation parameters on December 12, 2019

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

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

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

圖7 2019年12月12日BDSBAS格網電離層采用固定的去相關參數的GIVE包絡情況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反映了電離層異常天氣時的包絡性能。圖9(a)采用伯爾尼大學(AIUB)事后處理獲取的全球電離層延遲圖(GIM)繪制，從圖中可以看出當日電離層暴導致了GIM異常雙峰現象；圖9(b)為當日境內BDSABS監(jiān)測站所有穿刺點格網電離層殘差頻率分布直方圖，其偏度值與峰度值靠近公式(14)給出的閾值；圖9(c)為當天GIVE包絡圖。從圖9可以看出，殘差分布的偏度值峰度值會及時響應電離層暴等異?？臻g天氣，放大GIVE并保證包絡率仍在99.9%以上。

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

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

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

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

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

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

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

5 結 論

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