附加IRI模型约束的全球电离层建模及定位精度分析

杨 玲，周春元，苏小宁，李博峰

（1.同济大学测绘与地理信息学院，上海 200092；2.兰州交通大学测绘与地理信息学院，甘肃兰州 730070）

电离层延迟是全球导航卫星系统（Global Navigation Satellite Systems，GNSS）进行导航定位的一种重要误差源，也是致使一般差分GNSS系统的定位精度随用户和基准站间的距离增加而迅速降低的主要原因［1-2］。由电离层引起的GNSS导航信号时延可达数米甚至百米级，会严重削弱导航定位的精度和可靠性。目前GNSS可以长期连续监测电离层活动，基于双频观测数据利用载波相位平滑伪距的方法进行电离层总电子含量（total electron content，TEC）的提取被证明是一种精度较高的手段［3-5］。基于提取的高精度电离层TEC观测数据，利用球谐函数可建立全球电离层模型（Global Ionosphere Map，GIM），反映电离层TEC的时空分布特征和电离层活动规律，并为单频GNSS接收机用户提供电离层延迟改正，从而提高导航定位的精度［6-7］。

从1998年起，IGS（International GNSS Service）电离层工作组的各分支中心基于GNSS双频观测数据开始独立解算全球电离层产品，IGS提供的电离层产品为各分支中心产品的加权平均值［8］。由于在海洋和南极区域测站数目较少，在这些区域利用球谐函数进行全球电离层建模时存在异常值现象，已有相关研究利用不等式约束最小二乘来消除建模区域的TEC负值［9］。为弥补观测数据的不足，利用前一日的GIM产品或经验模型在数据缺失区域添加虚拟观测也被证明是有效的方法［10-11］。本文利用最新的国际参考电离层（International Reference Ionosphere，IRI）模型2016版本在观测数据不足的区域添加虚拟观测值进行约束，分析模型的改善效果，并利用单点定位结果评估添加IRI模型约束的有效性。同时，对于只存在C1、P2观测值的测站，在C1、P1码间偏差改正效果不佳的情况下，考虑将C1观测值作为低精度的P1观测值，增加全球的观测数据，分析建模效果并在定位中评估数据选取策略的必要性。

1 全球电离层建模基本原理

1.1 求解STEC的方法

GNSS原始观测量包含伪距和载波相位等，伪距观测方程可表示为

式中：Pi为载波i（i=1，2）上的伪距观测值；ρ0为卫星和接收机之间的真实距离；dion，i为电离层延迟量；dtrop为对流层延迟量；、δr，t为卫星和接收机钟差；、d r，i分别为卫星和接收机的码延迟量，s和r分别为卫星和接收机编号；c为光速；ε为随机噪声。

在进行电离层延迟和差分码偏差（differential code bias，DCB）估计时，通常对P1、P2观测值进行组合，得到组合的P4观测值为

式中：dion，1、dion，2分别 为L1、L2载波上的电离层延迟；Ds、Dr分别为卫星和接收机差分码偏差。为获取更精确的DCB估计，需利用载波L4（L4=L1−L2）观测值对P4观测值进行逐弧段的平滑处理。在进行载波相位平滑伪距之前，应利用MW组合或电离层残差组合检测历元中的周跳和粗差，并进行剔除［1］。

根据电磁波信号穿过电离层介质时群速、相速的关系，GPS测距码伪距和载波相位观测值所受到的电离层延迟误差大小相等、符号相反［12］。仅顾及f2项时，电磁波在电离层中传播时所受到的电离层延迟可表示为

式中：fi表示载波i的频率表示沿着信号传播路径s对电子密度Ne进行积分，并将其称为倾斜路径上的总电子含量I，I以TECU（total electron content unit）为单位，1 TECU表示每平方米含有1016个电子。将式（3）代入式（2）中，经过变换可得到

1.2 球谐函数模型

电离层分布在距离地面60～2 000km高度范围内，信号传播路径的高度角和方位角各不相同，为将倾斜路径的总电子含量转化到天顶方向，利用单层模型将路径上的所有自由电子都集中在F2层高度处。转换关系可描述为［1-2］

式中：R为地球半径，取6 371km；H为近似单层高度；α为常数，取0.978 2。

为适应电离层的区域特征，通常采用球谐函数模型进行全球电离层建模。球谐函数模型可表示为［13］

式中：V(β，s)为地面上空某一高度的球壳上纬度为β、经度为s处的垂直总电子含量；nmax为拟合时采用的最大阶数；β为电离层穿刺点地心纬度；s=λ−λ0是日固坐标系中穿刺点经度，λ、λ0分别为穿刺点的经度和真太阳时；A nm、Bnm为全球电离层模型系数；是标准化的勒让德多项式。

将式（4）、（5）和（6）进行结合并变换，可得到组合观测方程为

式中：A nm、Bnm为待求电离层模型参数；Ds、Dr即为待求卫星和接收机的差分码偏差。由于组合观测方程（7）是奇异的，通常添加所有GPS卫星DCB和为零的基准约束，实现观测方程的满秩，从而分离卫星和接收机的DCB。

1.3 IRI模型约束

IRI模型是根据大量电离层探空数据得到的电离层经验模型，可计算海拔50～2 000 km范围内的电子密度等参数，2 000 km以上至GPS卫星轨道高度范围内的电子密度可根据IRI模型外推获得［14-15］。2 000 km以下电子密度的方程式可描述为

式中：f(h)为高度h对应的电子密度；yj为高度hj对应的电子密度；w为待求系数。根据最小二乘解算系数w后可外推2 000 km高度以上的电子密度，累加50 km至GPS卫星轨道高度范围内的电子密度即可得到总电子含量。

由于IGS测站在大西洋、南半球部分海洋地区分布较为稀疏，导致该地区无穿刺点覆盖。在穿刺点空白区域逐历元利用IRI模型计算虚拟穿刺点处的VTEC作为虚拟观测值参与球谐函数建模，实现对全球电离层模型的约束。由于IRI模型计算的VTEC精度较低，通常降低虚拟观测值的权重，而平滑后的P4计算所得VTEC权重默认为单位权［16］。虚拟观测方程为

式中：βIPP、sIPP为所添加虚拟观测点的经纬度坐标。在进行模型约束时，通常需要经过一轮的解算，求得卫星和接收机的DCB参数，将式（7）转化为式（6）。考虑到卫星和接收机DCB短期内的稳定性，亦可利用前一日的卫星和接收机DCB参数代入式（7），得到式（6）中纯净的VTEC观测方程。最终联立式（6）和式（9）进行解算，便可得到附加IRI模型约束的全球电离层模型。

2 实验及分析

实验选取2019年10月2日（年积日为275）全球420个测站的观测数据，数据采样率为30 s，其中205个测站包含P1、P2伪距双频观测数据，其余215个测站P1观测值缺失，但存在C1、P2双频观测数据。实验中利用GPS单系统进行建模，采用4种建模方案：①基于P1、P2观测值的全球电离层建模，相关结果记为P1-P2；②在方案1的基础上，在RMS值超限区域利用IRI模型进行约束，相关结果记为P1-P2-Limited；③将C1观测值作为低精度数据取代缺失的P1观测值，相关结果记为C1-P2；④在方案3的基础上，在RMS超限区域利用IRI模型进行约束，相关结果记为C1-P2-Limited。

通过数据预处理，可获得方案1和方案3测站上空的电离层穿刺点分布情况如图1a和1b所示。可见，由于东太平洋和南半球海洋区域测站分布较少，其大气上方无电离层穿刺点分布。经过首轮解算，RMS值超限区域基本位于此范围内，因此方案2重点考虑在A、B、C、D 4个区域添加虚拟观测值，如图1a所示；方案4重点考虑在E、F、G 3个区域添加虚拟观测值，如图1b所示。

图1 电离层穿刺点分布Fig.1 Distribution of IPP

采用如表1中的设置及处理策略进行解算，得到卫星、接收机DCB和电离层参数，并利用电离层参数建立全球的电离层模型。

表1 设置及处理策略Tab.1 Strategies for setting and processing

2.1 模型解算精度分析

实验中，方案2和方案4加入了虚拟观测值，但方案1与方案2、方案3与方案4的DCB解算结果相同，仅在添加虚拟观测值后电离层参数的求解中存在差异。将解算的DCB结果与IGS提供的DCB产品进行对比，评定方案1与方案3中模型参数的解算精度。统计2种方案下卫星DCB的解算偏差及偏差分布如图2a和2b所示，接收机DCB的解算偏差及偏差分布如图3a和3b所示。

图2 方案1、3卫星DCB解算偏差对比Fig.2 Comparison of satellite DCB solution deviations of schemes 1 and 3

图3 方案1、3接收机DCB解算偏差对比Fig.3 Comparison of receiver DCB solution deviation of schemes 1 and 3

可以看出，利用P1、P2观测值时，由于精码伪距具有较高的精度，虽然参与解算的测站数较少，但卫星DCB的解算偏差均在0.1 ns之内，90%的接收机DCB解算偏差在0.6 ns之内；P1观测值缺失时，将C1观测值作为低精度的P1观测值后，所有测站的观测数据均可参与解算，90%的卫星DCB解算偏差在1 ns之内，90%的接收机DCB解算偏差在0.8 ns之内。对比C1作为低精度P1观测值前后，DCB解算精度有所下降，分析原因可能为：C1码本身测距精度较低，且C1码与P1码之间存在码间偏差，将C1作为低精度的P1观测值后，虽然观测数据增加，但由于添加的观测值精度较低，导致参数解算精度整体略有下降，但大部分的DCB解算偏差仍在1 ns之内，这也验证了模型参数解算的有效性。

将每种方案解算的12组电离层参数利用球谐函数模型展开，加入历元、经纬度信息可得到一天内任意时刻、任意经纬度的电离层VTEC值。图4分别统计了2019年10月2日0：00、6：00、12：00和18：00时刻不同方案的电离层建模结果与IGS产品的差值，在默认IGS提供的电离层产品精度较高的情况下，分析不同方案建模结果的优劣。

图4 附加IRI模型约束前后全球电离层VTEC对比（单位：TECU）Fig.4 Comparison of global ionospheric VTEC after addition of IRI model constraints(unit:TECU)

该图显示，方案1中，由于含有P1、P2观测值的测站在海洋区域数量较少且分布不均匀，导致在太平洋东部、大西洋中部和印度洋南部海域建模结果存在异常值。方案2中添加虚拟观测值后，6：00和12：00时刻在太平洋南部和印度洋南部海域的异常值范围显著减小，建模效果显著提升；0：00时刻在印度洋西南部的异常值改善效果显著，而太平洋大部分异常值改善效果甚微；18：00时刻2种方式的建模效果无显著差异。结合对比其他时刻的建模结果表明，添加虚拟观测值后，观测数据在全球得到充分补充，建模效果在部分区域提升明显，但较为稀疏的虚拟观测数据对模型的约束力并不强，使得改善效果存在不确定性，并不能从根本上解决观测数据缺失的问题。

方案3中将C1观测值作为低精度的P1观测值后，观测数据明显增多，与方案1相比建模异常值的范围显著缩小，全球整体建模效果提升显著，但在东太平洋无测站的海域依然存在异常值现象。方案4中添加虚拟观测值后，由于添加虚拟观测的范围有限且精度较低，与方案3相比各历元的建模效果无显著差异。

为客观分析建模效果，统计了4种方案各历元时刻电离层VTEC值与IGS参考值差值的平均值（Mean）、均方根（RMS）和标准差（STD）如图5所示，4种方案所有历元各统计量的均值如表2所示。

表2 不同方案全部历元统计量均值统计Tab.2 Mean value of all epoch statistics in differ⁃ent schemes（单位：TECU）

图5 附加IRI模型约束前后各统计量对比Fig.5 Comparison of statistics after addition of IRI model constraints

可以看到，与方案1相比，方案2中添加IRI模型约束后，8：00和20：00时刻的各统计量绝对值都有所减小，说明模型的内外符合精度都有提升，模型改善效果显著；其他历元时刻的Mean的绝对值均有所增大，RMS值均有不同程度的减小。统计所有历元时刻的Mean、RMS和STD平均值后发现，Mean平均值的绝对值由0.41增大到0.69，RMS和STD的平均值分别由4.22、4.14减小为2.99、2.85，模型内符合精度降低而外符合精度有所提高，分析原因可能为添加虚拟观测值后在观测数据较多的区域模型精度略有降低，但在观测数据缺失区域建模效果提升显著。所有历元附加约束前后STD值均有不同程度的减小，说明附加IRI模型约束可显著降低模型中异常值的概率，改善了建模效果。

方案3中利用部分C1观测值作为低精度P1观测值后，与方案1相比各历元Mean值的绝对值均有所增大，而RMS值除14：00时刻外均减小，说明加入额外的低精度观测数据后，模型内符合精度降低而外符合精度提高。方案4中添加IRI模型约束后，整体建模效果略好于方案3，但并无显著差异。总体而言，利用IRI模型添加虚拟观测值和利用C1观测值作为低精度P1观测值的目的都是为了增加观测量，扩大电离层穿刺点的覆盖范围，提升全球整体建模的效果，但是由于添加的多余观测数据精度往往更低，因此会降低电离层模型的内符合精度。

2.2 定位精度分析

为分析不同方案电离层建模效果的优劣，采用伪距单点定位的方式评定不同方案所得电离层产品的定位效果。考虑到全球IGS测站的分布及赤道附近电离层双驼峰结构的影响，以经度0°、纬度30°N和30°S为分界线将全球划分为6大区域并标记为A、B、C、D、E和F，每一区域选取6个测站的观测数据用于定位精度分析，区域划分及参与测试的测站分布如图6所示。

图6 区域划分及测试测站分布Fig.6 Area division and distribution of test stations

在定位中，分别加入4种方案所得到的电离层产品用于修正电离层延迟量，同时添加Klobuchar模型和IGS提供的GIM产品作为对比。为客观地分析不同方案电离层产品的定位效果，统计了各区域6个测站定位偏差的平均RMS分布如图7所示，各方向的平均定位误差（95%显著水平）如图8所示。可以看出，4种方案电离层产品的定位效果普遍优于Klobuchar模型，在观测数据较多的A、B、D和F区域，利用不同方案电离层产品的定位效果与IGS提供的GIM产品效果相当，而在缺少观测数据的C和E区域，不同方案电离层产品的定位效果表现出较大的差异。

图7 不同方案各区域平均RMS分布Fig.7 Average RMS distribution of different schemes in each area

可以发现，方案2中添加虚拟观测值之后，与方案1相比在测站数据较多的A、B、D和F区域各方向的RMS整体有略微增大的趋势，N、E、U方向RMS值增大的平均值为0.015 m、0.001 m和0.002 m。而在测站数据较少的C和E区域，各方向上的RMS显著减小，其中以E区域N、E方向最为明显，RMS分别减小0.109 m、0.073 m，U方向减小0.045 m。图8显示，添加虚拟观测值后，在A、B、D和F区域各方向定位误差整体略微增大，N、E、U方向定位误差平均增加0.016 m、0.003 m和0.012 m，在观测数据较少的C区域定位误差略微减小，而观测数据最少的E区域定位误差减小更为显著，N、E、U方向上分别减小0.260 m、0.146 m和0.103 m。说明在加入低精度虚拟观测值之后，对于观测数据较多的区域，重新解算电离层模型时参数精度有所下降，影响了该区域电离层模型的精确建立，最终导致定位效果略有下降。而在观测数据不足的南半球海洋区域，添加低精度虚拟观测值一定程度上弥补了观测数据不足的问题，改善了电离层参数中异常值的出现，增加了该区域电离层模型的整体符合度，并一定程度上提升了该地区的定位效果。

图8 不同方案各区域平均定位误差（95%显著水平）Fig.8 Average positioning error of different schemes in each area(95%significant level)

方案3中将C1观测值作为低精度P1观测值之后观测数据显著增多，与方案1相比，在观测数据较多的A、B、D和F区域RMS和定位误差有略微增大趋势，而在观测数据不足的C和E区域RMS和定位误差的增大趋势则更为显著，其中E区域的RMS在N、E、U方向上分别增加0.200 m、0.024 m和0.384 m，定位误差在N、E、U方向上分别增加0.524 m、0.111 m和1.549 m，方案3的定位效果整体不如方案1，说明利用C1观测值作为低精度的P1观测值后，虽然图4a和4c直观显示全球整体建模效果有所提升，但是一定程度上损失了模型的精度，导致在测试区域采用方案3的电离层模型的定位效果并不如方案1。

以上分析表明，4种方案中方案2整体上最优，与方案1相比，虽然在测站数目较多区域的定位效果有所下降，但下降幅度仅在厘米级，而在缺少观测数据的区域，定位效果的提升能达到分米级，且在严重缺少观测数据的E区域最为明显。因此，在利用载波相位平滑伪距的方法进行全球电离层建模时，选取利用包含P1、P2双频观测值的观测数据是非常有必要的，在观测数据严重不足的区域，利用IRI模型添加虚拟观测值能在保证定位精度的基础上提升模型整体的建模效果，而存在C1、P2观测值时，则不建议将C1观测值作为低精度的P1观测值参与电离层模型的建立。

3 结论

采用载波相位平滑伪距的全球电离层建模方法，分析对比4种电离层建模方案的效果差异，评定了不同方案所得到的电离层产品在单频伪距单点定位中的改善效果，最终确定了不同数据选取策略和添加IRI模型约束的必要性和有效性。

基于P1、P2原始双频观测数据添加IRI模型约束后，部分历元电离层建模整体效果改善显著，异常值出现范围显著缩小，但是该方法降低了模型的内符合精度，且无法从根本上解决观测数据缺失的问题。在定位效果上，添加IRI模型约束后，在观测数据不足的区域定位误差显著减小，N、E、U方向上分别减小0.260 m、0.146 m和0.103 m。因此在进行全球电离层建模时，利用IRI模型在观测数据不足区域进行约束是非常有效的。

将C1观测值作为低精度P1观测值后，电离层整体建模效果明显提升，但模型精度损失严重，数据不足区域的定位误差在N、E、U方向上分别增加0.524 m、0.111 m和1.549 m。因此进行全球电离层建模时，利用包含P1、P2双频观测值的观测数据是非常有必要的，而存在C1、P2观测值时，则不建议将C1观测值作为低精度的P1观测值参与电离层模型的建立。

作者贡献声明：

杨 玲：统筹论文的研究工作，指导论文的研究方向并修改论文。

周春元：查找文献、分析数据，负责论文的撰写。

苏小宁：提供部分测站数据及相关问题的咨询。

李博峰：参与部分数据分析，提供修改意见。