基于残差统计特性的中国区域格网电离层GIVE算法优化

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

1. 中国科学院上海天文台，上海 200235； 2. 32021部队，北京 100000； 3. 北京师范大学，北京 100036； 4. 中国科学院空天信息创新研究院，北京 100036

星基增强系统(SBAS)通过GEO卫星播发GNSS导航卫星轨道、钟差和格网电离层改正数以及UDRE和GIVE等完好性参数，用以提升区域用户GNSS导航服务精度与可靠性。1992年美国联邦航空管理局(FAA)提出了广域差分(WAAS)构想后，各国开始建设自己的星基增强系统。目前主要的星基增强系统有美国的WAAS、欧盟的EGNOS、日本的MSAS和印度的GAGAN。美国航空无线电技术委员会(Radio Technical Commission for Aeronautics,RTCA)2016年颁布的《全球定位系统/广域增强系统机载设备最低运行性能标准》(DO-229E)[1]详细描述了WAAS报文数据格式与内容。RTCA协议针对GPS L1 C/A频点单频用户进行增强。

北斗二号采用区域网定轨确定卫星轨道，基于星地双向时间同步确定卫星钟差[2-5]，利用北斗Klobuchar 8参数广播电离层模型修正电离层延迟以提供基本导航服务，北斗二号SBAS系统播发等效钟差改正数用于修正卫星钟差与轨道径向误差[6-7]，并以5°×2.5°空间分辨率播发格网电离层修正参数[8]。北斗三号采用区域网加星间链路定轨确定卫星轨道，利用星地星间双向时间同步确定卫星钟差[9-11]，采用北斗BDGIM 9参数广播电离层模型修正电离层延迟以提供基本导航服务[12]。2018年12月27日《北斗卫星导航系统公开服务性能规范》[13]发布，北斗三号全球系统开始提供基本导航服务[14]，随着北斗三号全球系统的逐步建成与完善，BDSBAS星基增强系统将成为日后工作重心。BDSBAS需要遵照RTCA协议增强GPS L1 C/A用户，同时需要增强BDS B1C用户。

电离层延迟是影响GNSS单频用户服务精度的主要误差源之一，除各GNSS系统播发的广播电离层参数外，星基增强系统也采用薄壳电离层模型提供更高精度的电离层延迟误差改正[15-16]。电离层延迟可通过双频伪距或精密单点定位提取[17]。WAAS播发电离层格网报文以修正用户电离层延迟[1，18]。WAAS利用基于梯度变化的克里金插值(KT)算法计算格网点垂直延迟，并基于卡方因子算法实现电离层异常状态的检验及GIVE参数的计算[19-20]。相较于平面拟合算法，基于克里金插值算法处理得到的格网电离层改正精度可提升3%～6%[21]。此后，基于最差用户几何位置的格网点延迟方差估值算法解决了由电离层风暴引起区域电离层梯度剧烈变化导致的服务完好性变差问题[22]；基于电离层层析技术的克里金插值后处理算法用以测量美国与巴西上空电子密度总含量[23]。印度GAGAN采用平面拟合算法计算格网点垂直延迟，采用与WAAS一致的算法计算GIVE，并且给出了适合赤道地区电离层活动的先验参考分布[24]。基于印度上空实测数据表明，相较于平面拟合算法，采用克里金插值计算格网点垂直延迟更为稳定[25]。适合中国区域的电离层球壳高度和格网电离层模型在后续研究中得到了详细论证[26-28]。陆态网实测数据表明，中国区域格网电离层采用平面拟合算法估算格网点垂直延迟RMS略低于克里金插值，但存在严重的边际效应[29]。

本文参考WAAS格网电离层计算算法，综合平面拟合与反比距离加权方法实现BDSBAS格网电离层格网点垂直延迟估值的计算。为避免粗差导致平面拟合估计值偏离实际值，根据克里金插值平稳假设，电离层活动正常状态下，利用小范围(5°×5°)内的穿刺点中值估计样本均值，利用样本方差序列中值估计样本方差剔除粗差。考虑到由电离层活动不规则导致的残差分布的非正态性，GIVE的计算考虑了残差分布的峰度系数与偏度系数，最终结合基于历史信息的卡方因子，能够实现99.9%以上的GIVE包络率并有效避免过包络。本文首先介绍了BDSBAS格网电离层容错算法与格网改正数解算算法；接着重点介绍了基于残差分布偏度与峰度系数的GIVE算法；最终用2020年1月BDSBAS监测接收机实测数据验证了算法的有效性并简单评估了格网电离层服务精度。

1 粗差剔除与残差统计

采用监测接收机双频相位平滑伪距观测数据提取电离层延迟，不可避免地存在粗差。粗差的存在会导致格网点垂直延迟估值失真，且显著增加平差迭代次数，导致格网电离层解算耗时增加。

采用上述粗差剔除算法，可显著降低迭代次数并提高格网电离层格网点垂直延迟估值的精度与可靠性。图1给出了采用中值容错前后的格网电离层残差分布：图1(a)为未采用中值容错，仅采用多次迭代方式剔除粗差的格网电离层残差频率分布直方图；图1(b)为采用中值容错后的格网电离层残差频率分布直方图；图1(c)为对应残差分布的Q-Q图，图中十字点线越靠近参考虚线，则表明样本分布为正态分布的可能性越大；图1(d)为对应残差分布的Q-Q图。

图1 采用中值容错前后的格网电离层残差分布Fig.1 Grid ionospheric residual distribution of median fault tolerance

由图1可知，采用中值容错剔除粗差后，Q-Q图十字点线更靠近参考虚线，样本分布偏度与峰度值更接近0；格网电离层残差分布会更加靠近正态分布。这说明了粗差剔除合理有效。

2 格网点延迟估计

(1)

(2)

同一卫星或者同一监测接收机的观测量间存在相关性，故令协方差矩阵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]反映了残差分布的非对称性，正态分布偏度系数为0。采样于未知分布的样本偏度值绝对值越大，表明该未知分布是正态分布的可能性越小。

偏度(skewness)α的计算方法为

(7)

式中，μ为未知分布期望，可用样本均值近似；ki为未知分布的i阶矩,可用样本i阶矩近似，计算公式为

(8)

式中，xj为采样于未知分布的样本；N为样本总数。

峰度(kurtosis)[30]值β反映了残差分布的拖尾性质，正态分布峰度系数为0。采样于未知分布的样本峰度值绝对值越大，表明该未知分布是正态分布的可能性越小，当峰度(kurtosis)值β小于0时，表示样本分布为厚尾分布。

峰度(kurtosis)β的计算方法为

(9)

采用平面拟合得到格网点垂直延迟后，根据RTCA协议[1]提供的格网电离层用户算法插值计算穿刺点垂直延迟，与双频提取电离层垂直延迟做差得到残差序列，根据残差序列的统计特性确定完好性参数GIVE。

图2为2019年11月23日—2019年11月25日双频提取与格网计算电离层延迟互差的频率分布直方图，期间未发生电离层暴，计算残差分布偏度值为0.116、峰度值为0.585。

图2 未发生电离层暴时电离层格网残差分布的偏度系数与峰度系数Fig.2 Skewness coefficient and kurtosis coefficient of the ionospheric grid residual distribution under ideal conditions

统计2019年10月20日—2020年3月18日的格网电离层残差分布偏度值与峰度值得到表1。从表中可以看出，99.9%的格网点残差分布偏度值绝对值在2.507以下、峰度值绝对值在9.459以下。当残差分布偏度值与峰度值偏离0时，说明残差分布为正态分布的可能性变小，基于正态分布假设给出的0.999概率置信区间已经不在可靠，需要扩大置信区间。

表1 格网电离层残差分布峰度与偏度值

3.2 边界方差的估计

为反映电离层格网点垂直延迟估值的可靠性，为用户提供最差包络，必须同时考虑残差分布特性与电离层历史信息，WAAS采用卡方因子来实现电离层历史信息的引入。格网点边界方差由式(10)[20]计算

(10)

(11)

卡方统计量计算公式为

(12)

去相关参数δdecorr被解释为平面拟合算法模型表达误差，采用去相关参数δdecorr可以有效提高包络率，但可能导致过包络问题。由于地磁赤道穿过中国南部地区，电离层梯度变化更为复杂，当采用与WAAS相同的去相关参数时，服务区南部地区部分格网点包络率低于99.9%，若扩大去相关参数会导致服务区内北部地区部分格网点过包络，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 卡方检验因子对格网电离层残差分布偏离正态分布时的响应情况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监测站数据入站率稳定，监测接收机IFB稳定，能够满足区域格网电离层解算需求。监测站均匀分布在中国境内，每个监测站均配备3台BDS监测接收机、3台GNSS多模监测接收机，可以同时解码获取GPS、GLONASS、Galileo 3系统观测数据与导航电文。目前仅利用北斗三号B1C-B2A、北斗二号B1I-B3I、GPS L1 C/A-L2P双频相位平滑伪距(CNMC)[31]提取电离层延迟参与格网解算。

按照本文所述算法剔除数据粗差后，再根据本文所述算法计算电离层格网点垂直延迟和GIVE，统计格网电离层改正残差RMS、改正百分比和GIVE包络率，以验证算法精度与可靠性。同时为考察BDSBAS格网电离层增强服务精度，利用本文算法计算得到的格网电离层产品，选取服务范围内的4个监测站伪距单点定位，并统计定位精度提升情况。统计结果时间系统均采用北斗时；格网电离层延迟单位为m。

采用2019年12月12日32个监测站双频提取电离层延迟，截止高度角设置为15°，利用三角投影函数[32]计算穿刺点垂直延迟。利用上文所述算法计算格网点垂直延迟，计算GIVE时不考虑去相关参数，统计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并统计一天包络率得到图7。图7中实线为GIVE包络线，散点为格网电离层残差序列点。图7(a)为30°N—55°N范围内的穿刺点包络图，图7(b)为0°N—30°N范围内的穿刺点包络图。由图7可以看出，若采用固定的去相关参数，对于中国区域0°N—30°N范围内的穿刺点GIVE过包络现象明显；对于中国区域低纬度地区仍存在不少的漏警点，包络率仅为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日发生一次小型电离层暴(http:∥www.nsmc.org.cn/NSMC)，图9反映了电离层异常天气时的包络性能。图9(a)采用伯尔尼大学(AIUB)事后处理获取的全球电离层延迟图(GIM)绘制，从图中可以看出当日电离层暴导致了GIM异常双峰现象；图9(b)为当日境内BDSABS监测站所有穿刺点格网电离层残差频率分布直方图，其偏度值与峰度值靠近公式(14)给出的阈值；图9(c)为当天GIVE包络图。从图9可以看出，残差分布的偏度值峰度值会及时响应电离层暴等异常空间天气，放大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日发生电离层暴时GIVE包络情况Fig.9 The GIVE envelope diagram during Ionospheric storm on January 20,2020

为进一步评价本文提出的GIVE算法包络性和电离层格网点垂直延迟解算精度，利用2020年1月1日—2020年1月31日连续一个月32个BDSBAS境内监测站实测数据，按照本文算法计算格网点垂直延迟和GIVE，从格网电离层残差RMS(分别与双频实测数据和CODE GIM产品对比)、改正百分比、GIVE包络率3个角度评价BDSBAS格网电离层服务精度。选取华中、华北、华东、华南、西北、东北、西南地区7个格网点，统计2020年1月整月RMS、改正百分比和GIVE包络率得到表2。

表2 BDSBAS格网电离层服务精度统计评价

为简单评估BDSBAS格网电离层服务精度，利用2020年1月1日—2020年1月31日连续一个月32个BDSBAS境内监测站实测数据，按照本文算法计算格网点垂直延迟和GIVE，选取位于华中(Sta1)、华北(Sta2)、华东(Sta3)、华南(Sta4)的4个监测站统计增强定位提升百分比得到表3。伪距单点定位采用本团队编写的性能监视程序，其双频相位平滑伪距基本导航定位精度在5 m(95%)以内，将事后精密单点定位(PPP)解算坐标作为参考真值考察增强定位提升百分比。

表3 BDSBAS格网电离层服务定位精度统计评价

由表2和表3可知，BDSBAS格网电离层改正残差RMS平均在0.4 m左右，大约在2～3 TECU量级；改正百分比均在75%以上，平均在78%左右。利用本文提出的GIVE算法，GIVE包络率均在99.9%以上，由于考虑了残差分布的偏度特性与峰度特性，GIVE包络在降低漏警率的同时抑制了过包络现象。利用BDSBAS监测接收机进行伪距单点定位，在格网电离层增强的状态下，对于BDS BIC频点用户，定位精度增强不明显；对于GPS L1 C/A频点用户，可增强定位精度20%～40%。GPS L1 C/A频点用户增强效果明显，是由于GPS基本导航用户采用了Klobuchar 八参数广播电离层模型，该模型改正百分比大约在50%左右；而BDS BIC频点用户采用BDGIM广播电离层模型修正电离层延迟，该模型改正百分比全球平均75%[27]，中国境内改正百分比更高，故定位精度提升不明显。考虑到格网电离层电文更新频率比广播电离层模型高，相比于精度提升，电离层异常活动及时预警对于BDS BIC频点用户更为重要。

5 结 论

本文提出了基于残差分布峰度系数与偏度系数的格网电离层GIVE算法。相比于采用固定的去相关参数计算GIVE，基于残差分布偏度值与峰度值自适应调整去相关参数计算的GIVE可以实现在降低漏警率的同时有效抑制过包络现象。新算法计算的完好性参数GIVE可以做到对服务区内任何纬度范围内的格网点实现99.9%以上的包络率。除此之外，本文同时介绍了BDSBAS容错算法与格网改正数计算算法。基于BDSBAS实测数据的统计表明，格网电离层修正RMS约在2～3 TECU，改正百分比达到75%～79%；修正格网电离层后可提升GPS定位精度20%～40%。