附加历元间约束的滑动窗单频实时精密单点定位算法
杨凯淳1
1. 信息工程大学地理空间信息学院, 河南 郑州 450001;
2. 哈尔滨工业大学(深圳)空间科学与应用技术研究院, 广东 深圳 518055;
3. 中国科学院精密测量科学与技术创新研究院大地测量与地球动力学国家重点实验室, 湖北 武汉 430077;
4. 西安科技大学测绘科学与技术学院, 陕西 西安 710054
基金项目：国家自然科学基金(42104033);深圳市科技计划(KQTD20180410161218820);广东省基础与应用基础研究基金(2021A1515012600);河南省重点研发与推广专项(科技攻关)(212102310428)
关键词：多历元递推算法 单频精密单点定位 约束条件 病态性 电离层延迟
杨凯淳, 吕志平, 李林阳, 等. 附加历元间约束的滑动窗单频实时精密单点定位算法[J]. 测绘学 ，2022，51(5)：648-657. DOI:
10.11947/j.AGCS.2022.20210207
YANG Kaichun, Lü Zhiping, LI Linyang, et al. Sliding window single-frequency real time precise point positioning algorithm with epoch constraints[J]. Acta Geodaetica et Cartographica Sinica, 2022, 51(5): 648-657. DOI: 10.11947/j.AGCS.2022.20210207
阅读全文：
http://xb.sinomaps.com/article/2022/1001-1595/20220503.htm
引 言
单频接收机由于价格低廉、操作方便，已广泛应用于变形监测、低轨卫星定轨及地面 解算等领域[1-3]，然而单频精密单点定位(precise point positioning, PPP)的定位精度与多频PPP相比，还存在一定差距，最大的难点在于电离层延迟的改正[4-6]。不少学者对其开展了很多研究，得到了有效的处理方法[7-9]。
目前，主要采取以下4种方法改正单频PPP的电离层延迟：①采用电离层模型，但利用GPS导航电文中的Klobuchar模型，只能消除50%~60%的电离层延迟，而IGS格 数据精度为2 TECU，在电离层活跃及IGS站很少的区域其精度更低[9-10]；②构成半合观测值，该方法只消除了电离层延迟的一阶项，且由于模糊度与接收机钟差高度相关，无法获得真实的钟差估值[11-12]；③将电离层斜延迟作为未知参数估计[13-16]，目前普遍采用的是文献[13]中的三参数法，同时估计对天顶延迟、水平梯度、映射函数，但模型过于简单，无法较好地描述电离层的活动，特别是在电离层延迟剧烈变化的区域；④通过多频接收机建立区域电离层模型，对单频接收机进行改正[17-18]，该方法在距离较远且条件不好时，不方便使用。若不使用电离层模型也不拟合电离层斜延迟，而是直接将其作为参数求解，不仅可以获取高精度的电离层延迟信息与测站坐标，还能在一定程度上提高收敛速度[14, 19]。
1 单历元单频精密单点定位算法
参数合并后的单频PPP模型为
(1)
式中，Pi、Φi对应卫星i的码和相位观测量；ρi为接收机到卫星的几何距离；c为光速；δtr为接收机钟差；Mi和dzpd分别为对流层延迟映射函数和天顶延迟；dioni为卫星i的电离层延迟；dr、dsi为接收机和卫星i的码硬件延迟；δr、δsi分别为接收机和卫星i的相位硬件延迟；λ为载波波长；Ni为卫星i的整周模糊度；εPi、εΦi分别为伪距和相位测量噪声；φr与φsi分别为接收机和卫星的初始相位。钟差参数δt吸收了接收机端硬件延迟，即δt=δtr+dr；电离层路径延迟参数Ii吸收了卫星端硬件延迟，即Ii=dioni+dsi；bi为模糊度参数，它包含了整数模糊度，卫星和接收机的硬件延迟以及初始相位，即bi=dsi－dr+δr+δsi+λ(Ni+φr－φsi)。
除了三参数法[13]和附加约束的单频PPP算法[14]，还可以添加模糊度约束[15]。文献[15]提出了同时估计电离层延迟的方法，将式(1)的未知参数分为两类
(2)
(3)
式中，k表示第k历元。
可将观测方程写为
(4)
式中，Qk为对角阵；Hg,k和Hy,k是对应未知向量gk和yk的偏导数矩阵。
该方法通过参数约化减少解算时间，但定位精度仍受限于约束方程的准确性。
2 多历元单频定位模型及解算方法
相对于用模型改正电离层延迟，直接对电离层延迟进行估计定位精度更高。约束条件中，先验信息的准确性及约束建立的客观性都会导致精度下降及收敛时间的增加。多频非差非组合PPP对电离层延迟的处理方式分为3种：①仅通过电离层模型对其进行修复[20]；②直接将电离层斜延迟作为参数进行估计[21]；③添加电离层延迟约束，并对电离层斜延迟进行参数估计[22]。
方法①与方法③已在单频PPP中得到应用，但秩亏导致无法直接采用方法②。在多频非差非组合PPP中，方法③通过添加电离层延迟约束，分离硬件延迟与电离层延迟参数，定位精度比方法①与方法②更高，而单频PPP中电离层延迟参数只含有卫星端硬件延迟，可通过DCB(differential code biases)文件进行改正，不添加电离层延迟约束就达到较高的定位精度，避免外部先验约束的影响。
对于方法②造成的秩亏问题，假设可用卫星数为n，则单历元观测的方程个数为2×n个，待定参数包括三维坐标、接收机钟差、对流层天顶湿延迟、n个站星斜径分量电离层延迟以及n个模糊度参数，总共5+2×n个参数，秩亏数为5。此时，可通过联合多个历元的观测值，在解决秩亏问题的同时不引入额外的误差。
2.1 观测方程
假定有m个历元，n颗卫星，线性化后的多历元单频PPP模型可表示为
(5)
滑动窗内的观测方程可简化为
(6)
(7)
式中，Y为联合m个历元的观测向量；A为联合m个历元线性化后的系数矩阵；Pmn为第m个历元第n颗可用卫星的伪距；φmn为第m个历元第n颗可用卫星的相位；x为联合m个历元的参数向量；ε为联合m个历元的观测误差向量；Ym为第m个历元n×1维观测向量；Am为第m个历元线性化后的系数矩阵；xm为第m个历元的参数向量；εm为第m个历元的观测误差向量。
未发生周跳时，整周模糊度数值不变，因此，在多历元联合解算时，每颗卫星只需解算一个模糊度参数，此时通过联合多个历元，即可得到唯一解。假设有n颗卫星，联合m个历元，则共有2×n×m个观测值，静态情况下，需要估计3个坐标参数、1个对流层湿延迟参数、m个钟差参数、n×m个电离层参数及n个模糊度参数，待求参数有4+m+n×m+n个，其自由度为n×m－1)－m－4；动态情况下有4×m+1+n×m+n个参数，自由度为n×(m－1)－4×m－1。此时，除坐标参数需要估计3×m个以外，其余参数设置均与静态情况相同。
若静、动态情况下解唯一，则可用卫星数与联合的历元数需要满足的条件为
(8)
式中，nstatic为静态情况下的可用卫星数；nkinematic为动态情况下的可用卫星数；INT(·)为取整函数。
2.2 随机模型
利用卫星高度角和观测值噪声来确定伪距和载波观测值的权，其中σP=0.1 m，σφ=0.001 m。
(9)
式中，QYmYm为第m个历元观测值之间的方差-协方差阵；QY1Ym指第1个历元的观测值与第m个历元观测值之间的方差-协方差阵。
(10)
因此，自相关函数中元素的自相关系数可以导出为
(11)
式中，1≤τ≤m－1;
是第i历元的偏差；en=[1 1 … 1]n×1T；σi与σi+τ分别为第i与第i+τ历元观测值的中误差；ri与ri+τ分别为第i与第i+τ历元观测值的残差，并且满足ri=ri+τ=n－1。
对于参数的随机模型，静态情况下，坐标固定且先验约束为100 m，动态情况下，坐标采用白噪声模型且先验约束也为100 m；接收机钟差采用白噪声模型；对流层干延迟采用Saastamoinen模型改正，其湿延迟采用随机游走模型估计；电离层延迟采用随机游走模型；模糊度当作常数处理，且为了克服接收机钟差、电离层延迟和相位模糊度之间高度相关导致的列秩亏问题(秩亏数为1)，选择第1个滑动窗的第1组模糊度作为S基准，并将其约束到近似值，该近似值通过单频的相位观测值减去伪距观测值得到[24-25]。
对于不同历元间参数的先验方差-协方差阵，不同测站确定的先验方差-协方差阵不同，这里以SHAO站DOY 91的观测数据为例，首先通过同时估计电离层单频PPP算法100个历元的计算结果，确定相邻历元间钟差及电离层参数的相关系数，结果见表 1、表 2；然后根据其确定多历元递推单频PPP算法的先验方差-协方差阵。
表 1 历元间钟差的相关系数
Tab. 1 Correlation coefficient of clock correction between epochs
表 2 历元间电离层延迟的相关系数
Tab. 2 Correlation coefficient of ionospheric delay between epochs
在确定先验方差-协方差阵后，进行多历元递推单频PPP的解算，图 1给出了不同历元间各个钟差及电离层延迟参数之间的相关系数。将图 1与表 1、表 2进行对比，可以看出多历元单频PPP比单历元单频PPP参数之间的时间相关性要强，在进行多历元联合解算的过程中必须要考虑。
图 1 历元间电离层延迟与接收机钟差的相关系数
Fig. 1Correlation coefficient of ionospheric delay and receiver clock correction between epochs
2.3 参数估计方法
将观测方程式(6)中的参数按被约束与不被约束分组，写成误差方程形式
(12)
式中，a为不被约束的m个历元的参数；d为被约束的参数；B与C分别为a与d的系数行满秩矩阵；l为误差向量。
约束条件方程写为
(13)
式中，D为秩为r的r×r维设计矩阵，r为被约束的参数个数；E为r×1维的常量向量。然后分别对
及
进行量测更新与时间更新[15]。
考虑到计算效率以及计算结果的精度，预先对不同约束组合进行试验，并根据试验结果分别在动静态情况下对各参数进行不同的约束，见表 3。
表 3 各参数约束情况Tab. 3 Constraints of each parameter
表 3共联合了5个历元的观测值，静态情况下只对接收机钟差与电离层延迟进行约束，接收机钟差与电离层延迟的约束历元数都为4个，分别表示对前4个历元的接收机钟差与电离层延迟进行约束；动态情况下只对坐标与接收机钟差进行约束，其约束历元数都为3个。
图 2 多历元递推单频PPP算法流程
Fig. 2The flow chart of multi-epoch recursive single-frequency PPP algorithm
(1) 预处理，通过多普勒积分法与相位伪距组合法探测并修复周跳[26]，利用数据探测法探测并修复粗差，探测与修复钟跳，并剔除卫星数据不好的时段。
(2) 联合预处理后的多历元观测信息，建立观测方程并将其线性化。
(3) 根据自相关系数法确定观测值之间的方差-协方差阵，通过同时估计电离层单频PPP算法100个历元的计算结果确定参数的先验方差-协方差阵。
(4) 进行参数解算，选择第1个滑动窗的第一组模糊度作为S基准，并使用带约束的SRIF算法进行解算。
3 试验分析
图 3 测站分布Fig. 3Stations distribution
3.1 静态定位
图 4 多历元与单历元算法的PDOP值及可见卫星数对比
Fig. 4Comparison of PDOP value and visible satellites between multi-epoch and single-epoch algorithm
图 4为多历元算法与单历元算法在单天内的PDOP值及可用卫星数的变化。在多数情况下，多历元比单历元的可用卫星数多且卫星几何分布更好，但由于多历元算法对数据质量要求较高，需要剔除数据质量不好的卫星，在部分历元，多历元算法反而比单历元算法的PDOP值低且可见卫星数少。另外，对未经预处理的原始数据进行试验，发现多历元递推算法更易受到周跳及粗差的影响。
经过以上分析可知，多历元算法比单历元算法更好。使用传统的间接平差对多历元联合算法进行数据处理，可以得到图 5中的定位结果。由图 5可知，定位结果较差，E分量最大达到了1 m，N分量在0.5 m左右波动，U分量大部分在1 m左右波动。根据2.3节的分析可知，间接平差只处理当前存储的数据，并没有充分利用上一次解算的结果，若当前解算的m个历元的数据质量较差就会导致偏差很大，在U分量上更为明显。
图 5 SHAO站多历元联合的单频PPP(间接平差)
Fig. 5SHAO station single-frequency PPP based on multi-epoch (indirect adjustment)
图 6 3种方案下SHAO站的静态定位精度Fig. 6Static positioning accuracy of SHAO station under three different schemes
对比图 6(b)与图 5可知，相比使用间接平差的多历元联合单频PPP算法，采用多历元递推单频PPP算法拥有更高的定位精度，水平分量的精度能够达到厘米级，且在高程分量上有较大改善。
剔除部分存在跳变的历元，综合15个测站14 d的定位结果，统计了静态定位的RMS、STD及平均收敛历元数，见表 4。方案(b)中多历元递推算法的E、N和U分量上的RMS分别为0.017、0.008和0.028 m，优于方案(a)。方案(b)的收敛时间大大缩短，与方案(c)双频无电离层PPP的收敛速度基本保持一致。
表 4 静态定位平均RMS、STD和收敛时间
Tab. 4 Average RMS, STD and the average convergence time of static positioning
声明：本站部分文章及图片源自用户投稿，如本站任何资料有侵权请您尽早请联系jinwei@zod.com.cn进行处理,非常感谢！