GNSS观测值域增强PPP定位中综合改正数质量控制算法

洪菊，涂锐,3，张世旋，张鹏飞，王思遥，李芳馨，刘明玥

(1.中国科学院国家授时中心,西安 710600；2.中国科学院大学,北京 100049；3.中国科学院精密导航定位与定时技术重点实验室,西安 710600)

0 引 言

实际观测中，全球卫星导航系统(GNSS)观测值由于受到观测信号、传播路径误差、接收机等影响不可避免的存在粗差.若不对异常值进行剔除或控制，势必会造成参数估值的扭曲[1].因此，异常值的探测与消除是GNSS 质量控制的重要工作.同时，对于未准确探测出来的小周跳也会影响GNSS精密定位的结果.社会生产和生活对高精度导航定位的需求促进了多系统GNSS 以及多种高精度、高可靠性定位技术的发展[2-6].观测值域增强精密单点定位(PPP) 技术利用参考站提取的综合改正数实现快速高精度的PPP 定位，具有收敛快、精度高等特点而受到广泛关注.该项技术中，参考站的综合改正数作为虚拟观测值直接参与用户端PPP 解算，综合改正数中的异常对用户端定位性能会产生直接影响.

当前，对粗差的质量控制，主要有两种方法，第一种是将粗差归入函数模型，将含粗差的观测值看作与其他同类观测值具有相同的方差，但数学期望发生改变，通过构造检验统计量来对粗差进行判断和剔除[7]；第二种是将粗差归入随机模型的方差膨胀模型，将含粗差的观测值看作与其他同类观测值具有相同的数学期望，但其方差发生变化，即抗差估计[8].通过等价权来对粗差观测值进行控制，以消除或减弱异常误差的影响，如IGG 方案、双因子等价权以及等价方差-协方差函数等[9-11]，但抗差估计处理对于海量观测数据会增加计算负担[12].

基于这些前期的理论基础和观测值域增强PPP 技术中综合改正数特性，本文在数据预处理阶段，对综合改正数进行频率间和二阶历元间差分，对得到的差分组合值采用中位数法探测可能存在的粗差或者周跳，并对使用异常综合改正数的卫星进行模糊度重初始化、降权或剔除控制处理，从而有效控制综合改正数异常对定位结果的影响.

1 观测值域增强PPP 和综合改正数质量控制方法

1.1 观测值域增强PPP 算法

单频相位观测值可表达为

式中：下标b和i分别为测站和频率；上标 s 为卫星；Rsb为几何距离；Lsb,i为相位观测值；dts为卫星钟差；c为光速；λi为频率i的波长；dtb为接收机钟差；为i频点的模糊度；为i频点电离层延迟误差；ZWD为对流层天顶湿延迟；db,i和dis分别为接收机和卫星i频点相位硬件延迟偏差；为观测值噪声以及多路径效应和未模型化误差.对流层干延迟、相位缠绕误差、天线相位误差以及潮汐误差等均使用相应模型改正，未在公式中体现.精密卫星钟差产品包含消电离层组合形式的卫星伪距硬件延迟偏差，即

已知参考站坐标，引入卫星星历与卫星钟差产品，消除与测站相关的可模型化误差和站星几何距离后，以载波观测值为例，综合改正数可表达为

式中：orbsb和sclksb分别为卫星轨道和卫星钟差残余量；omcsb,i为观测值减去计算值(observed minus computed，omc)，称之为观测值域综合改正数.

流动站上的改正数可以根据其周围参考站估计的非差omc线性插值得到，当选择三个参考站(b1、b2、b3)时非差omc内插值omcsL,i,b可表达为[13]

式中，A、B、C为内插系数，计算方法在此不在赘述.

在站网范围不大、大气变化不太剧烈的情况下可以认为内插后的改正数与流动站真实值差距很小，因此内插后的综合改正数应用在流动站且忽略观测噪声和多径效应后，单频载波观测值可表达为

其中：

1.2 观测值域综合改正数质量控制方法

相位电离层残差组合二阶历元间差分受电离层残差等各项误差影响较小，因此常被用来探测周跳[14].同样，不同频率相位观测值域改正值之差可以表示为(以第1、2 频率为例)

由式(6)可知，该组合值消去了星历误差、对流层延迟等与频率无关的误差，只包含频率间电离层延迟、模糊度参数、多路径效应和观测噪声、卫星端未校准相位延迟(UPD)以及接收机端UPD.硬件延迟和多径随时间变化比较缓慢[15]，因此，对频间差分改正数进行历元间二次差后得到的二阶历元间差分组合值 Δ∇omcsL,GF,b在连续弧段下主要受观测值精度和电离层残差的影响，在零值附近波动.因此，通过对上述差分组合值序列检验识别综合改正数中存在的异常.考虑到观测噪声和电离层残余量会随着历元间隔变化而不同，采用中位数(MAD)探测方法，定义为[16]

其中，因子0.674 5 使得MAD 等同于正态分布数据的标准偏差.当数据在 [m-n·MAD,m+n·MAD] 范围之外时，被识别为异常值，整数n为常量.若用过于严格的判别标准导致的错判会影响解的精度和效率，因此本文取n为10.

根据式(4)和误差传播定律得

根据经验值[17]，载波测量中误差均为mφ=±0.01周，以4倍检测量中误差为限差作为检验异常的最小阈值判断条件，综合中位数判别条件在数据预处理阶段对存在的异常值进行判别.

结合上述理论方法，本文采用的异常值识别、定位与控制流程如图1所示.首先计算参考站第j历元综合改正数omcsL,GF,b(j)，并得到二阶历元间差分组合值Δ∇omcsL,GF,b(j).利用公式(7)～(9)确定异常识别阈值，若未超出阈值且上一历元无异常，则正常解算，否则根据上一历元组合值特性分为两种情况控制异常值：(a)当上一历元 Δ∇omcsL,GF,b(j-1) 为异常值时，若Δ∇omcsL,GF,b(j) 与 Δ∇omcsL,GF,b(j-1) 大小相近、符号相反，判断第j-1 个历元存在未探测出的周跳，参数估计时对该卫星模糊度进行重新初始化，同时从第j个历元开始重新计算差分组合值序列，反之，采用将Δ∇omcsL,GF,b(j-1)设为未存在异常的相邻历元差分组合值等方法反算 Δ∇omcsL,GF,b(j)，利用公式(7)～(9)确定异常识别阈值，若超出阈值则标记为异常，对卫星进行降权或者剔除处理；(b)当上一历元Δ∇omcsL,GF,b(j-1)无异常时，对卫星进行降权或剔除处理.

图1 综合改正数异常值识别、定位与控制流程

2 试验分析

为分析综合改正数异常对定位结果影响以及本文所提方法对异常值的探测情况，作者选用香港CORS 测站2021年第141 天HKST、HKSL、HKOH和HKSC 测站数据，HKST、HKSL 和HKOH 作为参考站构成参考网，平均边长约为26 km，HKSC 作为流动站，测站分布如图2所示.用户站与参考站平均距离为15 km，利用多个参考站内插得到的用户端大气延迟与真实大气延迟空间相关性强，实验中忽略大气延迟残余误差的影响.

图2 测站分布图

对不同采样率下的相位omc二阶历元间差分组合值特性分析以探究异常值探测水平.图3展示了该组数据采样间隔为1 s、15 s、30 s 时相位omc二次历元间差分组合值序列.1 s、15 s 和30 s 采样间隔下各卫星组合值始终在零值附近波动，其平均标准差分别为0.009 周、0.021 周和0.025 周.受观测噪声和电离层残差等影响，1 s、15 s、30 s 采样率下组合值的标准差有所不同，随着采样间隔的变长历元间差分组合值变大.

图3 1 s、15 s、30 s 采样间隔下相位omc 二阶历元间差分组合值时间序列

以采样间隔30 s 数据为例，分析不同粗差大小对仿动态增强PPP 定位结果的影响.选用单天数据，每8 h 重新解算.随机选取每个时段第960 个历元某颗卫星按以下方案进行测试.方案一，人为在第一频点L1 改正数上加入0.5 周粗差，分析较小粗差对GPS动态增强PPP 的影响情况；方案二，人为在第一频点L1 改正数上加入2 周粗差，分析较大粗差对GPS 动态增强PPP 的影响情况；方案三，利用无异常值的综合改正数进行增强PPP 解算.图4展示了参与解算的卫星以及位置精度因子(PDOP)值，卫星数量平均为7.4 颗，PDOP 值平均为2.04，卫星几何分布较好.

图4 仿动态PPP 增强解算中PDOP 值和卫星数量

由图5可知，相比于平面方向，加入0.5 周小粗差对天顶(U)方向影响较大，但是同样大小粗差对定位结果的影响不同，第1 和第3 时段对U 方向影响在1～2 cm，在第2 个时段对U 方向影响约为1 dm，这与观测值的权有关.一般认为高度角较大时，观测值多路径效应、对流层延迟及噪声等误差较小，GNSS 数据处理中通常采用卫星高度角相关的随机模型[18-20].第1～3 时段所选卫星分别为G28、G13 和G29，加入误差时对应高度角分别为17°、75°和20°，所以加入相同误差时第2 个时段对定位结果影响较大.当加入2 周粗差后，结果如图6所示，对东(E)和U 方向定位精度均有影响，对U 方向定位精度影响最大，第2 个时段影响最大可达0.3 m.因此选择合适的质量控制方法是十分必要的.

图5 加入0.5 周较小异常时PPP 增强定位结果序列

图6 加入2 周较大异常时PPP 增强定位结果序列

中位数法中常量取值不同对异常探测能力不同，本文放宽检验标准，n取10.图7展示了不同卫星在不同时段的探测阈值，图7(a)、(b)分别代表阈值上限和下限，若omc二阶历元间差分组合值超过阈值被判定为异常，即内插后的综合改正数存在的异常引起差分组合值超过阈值时可被识别.图中不同颜色代表不同卫星，由图可知，不同卫星在不同时段异常探测能力不同，当n取10 时可探测出差分组合值中大部分1 周以上的较大粗差和部分1 周以内的小粗差.表1结果表明利用该方法能够有效探测出加入的0.5 周和2 周粗差.但是根据图7所示，某些时段阈值较大，不能完全探测出0.5 周小粗差，对于未探测出的异常可考虑调整n大小或者采用抗差估计消除或减弱异常误差的影响.

图7 异常探测方法阈值(不同颜色代表不同卫星)

表1 异常值探测情况

为了进一步检验上述方法的有效性，笔者选用科廷大学CORS 网2018年第182 天CUT0-CUT2 零基线数据进行仿动态实验.以CUT2 为参考站，图8展示了CUT0 测站未进行异常值探测时PPP 增强解算的结果.可以看出，后半段定位结果存在分米量级的异常抖动.经分析，G21 卫星综合改正数存在异常，较大异常值处的历元二次历元间差分组合值为-0.284 周，下一历元其值为0.289 周，因此判断异常处发生周跳.其对当前历元以及后面所有历元定位结果造成高达3 dm 的误差.因此，本文对该异常值处卫星降权处理并对下一历元模糊度重新初始化.图9展示的是未进行异常值控制与进行异常值控制后PPP 增强定位结果的对比，由图9可知，对综合改正数进行质量控制后，PPP 增强定位结果精度有了明显提升.

图8 未进行异常值识别与控制时PPP 增强定位结果序列

图9 异常值控制前后PPP 增强定位结果比对

3 结论及展望

综合改正数质量直接影响PPP 增强定位精度，因此，对综合改正数中异常值的探测与控制在观测值域增强PPP 定位中具有重要作用.本文对观测值域增强改正数质量控制方法展开研究，针对综合改正数特性，提出利用经频间和二阶历元间差分后的差分组合值，作为检验量进行异常值识别，并对使用该异常值的卫星采用模糊度重新初始化、降权或剔除方法进行处理.差分组合值受观测值精度的影响，因此可以通过计算组合观测值精度探测综合改正数中可能存在的异常，但随着历元间隔变化观测噪声和电离层残差会变化同时避免过大的异常观测值对检验所产生的破坏性影响，选用中位数法设置探测阈值.中位数法常量设置决定了探测水平，为适当放宽检验标准，本文常量取值为10.经实验验证该方法能够有效探测出差分组合值中大部分1 周以上的较大异常和部分1 周以内的异常，并针对实验中出现的异常通过判断异常类型采用模糊度重新初始化方法抑制了异常对定位的影响.在GNSS 定位解算中，通常采用高度角定权模型，同时经过实验分析可知，不同高度角卫星的综合改正数存在的异常对定位结果影响不同.因此，针对不同高度角卫星存在的异常值处理方法仍需进一步探究.上述方法存在一定的局限性，比如无法识别不同频点改正数存在的相同异常，下一步继续研究结合MW 组合用于异常识别的方法.

致谢：感谢科廷大学GNSS-SPAN 和香港特别行政区政府地政总署测绘处开源的CORS 数据，本文的研究得到了国家自然科学基金项目(41974032)的支持.