-
随着全球导航卫星系统(Global Navigation Satellite System，GNSS)系统稳步建设，已有GPS、GLONASS、Galileo、BDS系统提供导航定位服务。2003年，国际GNSS服务(international GNSS service，IGS)组织成立了多GNSS工作组，用于多GNSS系统数据收集及产品生成[1]。目前，IGS组织提供可供事后多GNSS精密定轨的全球分布测站已经超过200个。另外，国际GNSS监测与评估系统(international GNSS continuous monitoring and assessment system, iGMAS)[2]、中国大陆构造环境监测网络工程、澳洲连续运行参考网也提供了充足的精密定轨区域GNSS观测站，为多系统GNSS实时精密轨道产品解算提供了支撑。多GNSS系统融合数据处理一方面能够有效减弱系统误差，降低多路径误差影响，提高解的精度、时空分辨率及可靠性，另一方面，大量增加的地面观测值及新兴的卫星系统的加入增加了计算模型的复杂度和参数解算的计算时间，降低了多GNSS系统精密定轨解算效率[3]。随着GNSS技术的快速发展，全球实时高精度应用成为了重要的研究方向，对多GNSS系统实时精密轨道解算精度和更新效率提出了较高要求。
GNSS卫星轨道受轨道动力学约束较为稳定，短期内轨道预报精度损失较少，目前广泛采用数值积分方法对卫星轨道进行预报。Beutler等[4]首次提出合并过去连续3个单天定轨弧段法方程方法，通过综合轨道参数的方法达到与事后3天弧段定轨精度一致的效果; 随后，Brockmann[5]借鉴该思想，将该方法应用到卫星轨道、地球自转参数及参考框架的研究中，提高了相关参数的解算精度; Dong等[6]研究了定轨弧长对地球形变参数的影响; Andersen[7]对多个定轨弧段合并中随机参数先验信息兼容问题进行了研究，提出一种合并最小二乘信息滤波与平滑的方法，实现了常量参数、白噪声参数及有色参数的前后一致，该方法在保证参数解算精度的前提下也提升了解算效率; 姚宜斌[8]通过法方程综合算法将单天解算轨道综合成三天解，进一步采用动力平滑算法进行轨道平滑，可以提高解的精度和可靠性; 楼益栋等[9]采用滑动窗口多弧段法方程叠加方法，通过单天滑动窗口解算的GPS轨道精度与IGS快速轨道产品精度相当; Lutz等[10]研究了单天解与多天合并解的轨道精度，并着重分析了其对地球自转参数的解算精度的影响; 刘伟平等[11]研究表明，在弧段连接处引入伪随机参数可以有效改善轨道跳变的现象。上述研究主要集中于事后产品精度、实时轨道产品服务方面，法国国家太空中心于2016年开始提供实时轨道钟差改正服务[12]，其实时轨道基于滤波算法解算，轨道更新时效性较好，但是受滤波算法影响，参数估计只能顾及“过去”时刻，精度略逊于最小二乘方法。Li等[13]基于最小二乘批处理方法研究了超快速实时轨道计算方法，但仅从定轨弧长、数据采样率、提升计算机硬件性能方面进行了优化，未能从解算方法上改进原有算法。针对实时轨道解算问题，本文提出一种基于分块递推最小二乘配置的多系统卫星精密轨道解算方法，该方法可同时兼容传统的最小二乘批处理方法[14]与实时滤波解方法[15]，实现四系统轨道超快速更新，有效解决了实时轨道产品更新精度与时效性的问题，相比于传统事后批处理定轨预报方法有显著提升。得益于高效的解算方法和策略，目前多GNSS系统轨道可实现1 h更新，延迟0.5 h发布，并已提供公开服务，服务产品包括0.5 h、3 h、6 h更新四系统实时轨道、钟差及ERP产品，相关产品可通过武汉大学IGS数据中心FTP服务获取(ftp://igs.gnsswhu.cn/pub/whu/MGEX)。
-
GNSS卫星精密定轨的函数模型可概括为运动方程和观测方程两部分[16]，在惯性系下, 运动方程为：
$$ \ddot{\bar{r}}={{\bar{a}}_{g}}+{{\bar{a}}_{ng}}+{{\bar{a}}_{\text{emp}}} $$ (1) 式中，r为卫星位置向量；ag为卫星所受保守力之和；ang为卫星所受非保守力之和；aemp为经验加速度力。其中力模型及偏导数计算可通过相关模型获取[16-18]。进一步将运动方程转换成一阶微分方程为：
$$ \left\{ \begin{array}{l} \dot r = v\\ \dot v = {{\bar a}_g} + {{\bar a}_{ng}} + {{\bar a}_{{\rm{emp}}}}\\ \dot p = 0 \end{array} \right. $$ (2) 式中，v表示卫星速度；p表示动力学模型中的待估参数，如太阳光压参数等。式(2)中，给定卫星初始参数(卫星初始位置、速度及动力学待估参数)，通过数值积分即可求得初始状态参数与指定时刻卫星状态的转换关系[19]。
假设具有先验参数约束的GNSS观测方程为：
$$ \left\{ \begin{array}{l} \mathit{\boldsymbol{L}} = \mathit{\boldsymbol{BX}} + \mathit{\boldsymbol{ \boldsymbol{\varDelta} }}\\ {{\mathit{\boldsymbol{X'}}}_{i + 1}} = \mathit{\boldsymbol{G}}{{\mathit{\boldsymbol{X'}}}_i} \end{array} \right. $$ (3) 式中，L为观测向量，即GNSS载波相位和伪距观测值；X为倾向参数，如卫星轨道参数、测站坐标参数等；X′为随机参数，如卫星钟差参数、对流层参数等，一般可用高斯-马尔科夫过程描述其变化；B为设计矩阵；G为参数配置数学模型矩阵；Δ为随机误差。对式(2)中轨道初始参数积分求得卫星位置并代入到式(3)中，对式(3)进行线性化，通过状态转移矩阵建立轨道初始状态参数与指定历元的联系，即可得到可用于线性估计理论求解的观测方程，进一步选择适当的线性估计方法进行参数解算[20]。精密定轨解算示意图如图 1所示, 通过状态转移矩阵逐个将卫星位置信息与初始状态参数建立联系，同时剔除个别粗差轨道位置数据，估计初始状态参数。由于定轨过程中综合考虑了轨道力学约束和卫星几何观测约束，因此称这种定轨方式为简化动力学定轨[20]。中高轨导航卫星受轨道力学约束较为稳定，短时间内轨道衰减较小，因此一般通过对卫星状态参数积分即可获得实时轨道。
图 1 GNSS卫星精密定轨原理示意图
Figure 1. Diagram of GNSS Satellites Precise Orbit Determination
-
本文提出的分块递推最小二乘配置方法可分为3个步骤：累积实时数据形成法方程与历史数据法方程叠加、轨道相关参数松弛、状态参数递推，具体实现过程如下。
考虑式(3)并顾及待估参数的先验精度，最小二乘估计准则[15]表示为：
$$ \left\{ \begin{array}{l} J\left( \mathit{\boldsymbol{x}} \right) = {\mathit{\boldsymbol{v}}^{\rm{T}}}\mathit{\boldsymbol{Pv}} + {\mathit{\boldsymbol{\bar x}}^{\rm{T}}}\mathit{\boldsymbol{ \boldsymbol{\varLambda} \bar x}} = \min \\ \mathit{\boldsymbol{v}} = \mathit{\boldsymbol{Bx}} - \mathit{\boldsymbol{l}},\mathit{\boldsymbol{\bar x}} = \mathit{\boldsymbol{x}} - \mathit{\boldsymbol{\tilde x}} \end{array} \right. $$ (4) 式中，P为观测方程先验权；x和${\mathit{\boldsymbol{\tilde x}}}$分别为待估参数及其先验无偏估计量；Λ为先验参数信息矩阵；B为观测方程关于待估参数的偏导数矩阵。为使二次型表达式(4)最小，令其一阶导数为0，整理后可得：
$$ \left( {\mathit{\boldsymbol{ \boldsymbol{\varLambda} }} + {\mathit{\boldsymbol{B}}^{\rm{T}}}\mathit{\boldsymbol{PB}}} \right)\mathit{\boldsymbol{\hat x}} = \mathit{\boldsymbol{ \boldsymbol{\varLambda} \bar x}} - {\mathit{\boldsymbol{B}}^{\rm{T}}}\mathit{\boldsymbol{l}} $$ (5) 在进行参数估计时，可利用事后获取的先验信息矩阵对后续待估参数求解进行配置。实时数据处理时，将观测数据按时间累计的方式分块，现将第一块数据采用最小二乘批处理形成的法方程表示为：
$$ \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{N}}_{11}}}&{{\mathit{\boldsymbol{N}}_{12}}}\\ {{\mathit{\boldsymbol{N}}_{21}}}&{{\mathit{\boldsymbol{N}}_{22}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{X}}_1}}\\ \mathit{\boldsymbol{Y}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{W}}_1}}\\ {{\mathit{\boldsymbol{W}}_2}} \end{array}} \right] $$ (6) 式中，X1为轨道状态参数；Y为其他参数；N为法方程系数；W为观测误差项。采用分块递推最小二乘配置方法时，需将新观测数据生成的法方程与之前的法方程叠加。当累计观测数据弧长达到指定定轨弧长时，考虑到陈旧数据累计会造成观测数据饱和，并引起新观测数据对解的贡献降低问题，需对轨道状态参数进行松弛处理，公式为：
$$ {\mathit{\boldsymbol{X}}_2} - {\mathit{\boldsymbol{X}}_1} = 0,\mathit{\boldsymbol{P}} $$ (7) 式中，X1为历史参数；X2为松弛后参数；P为松弛因子。
将式(7)作为虚拟观测方程，通过附有限制条件的最小二乘配置方法形成法方程，并叠加到式(6)中，采用参数消除方法[21]消除历史参数，可得参数松弛后的法方程：
$$ \left[ {\begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\tilde N}}}_{11}}}&{{{\mathit{\boldsymbol{\tilde N}}}_{12}}}\\ {{{\mathit{\boldsymbol{\tilde N}}}_{21}}}&{{{\mathit{\boldsymbol{\tilde N}}}_{22}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{X}}_2}}\\ \mathit{\boldsymbol{Y}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\tilde W}}}_1}}\\ {{{\mathit{\boldsymbol{\tilde W}}}_2}} \end{array}} \right] $$ (8) 式中，
$$ \left\{ \begin{array}{l} {{\mathit{\boldsymbol{\tilde N}}}_{11}} = \mathit{\boldsymbol{P}} + \mathit{\boldsymbol{P}}{\left( {{\mathit{\boldsymbol{N}}_{11}} + \mathit{\boldsymbol{P}}} \right)^{ - 1}}\mathit{\boldsymbol{P}}\\ {{\mathit{\boldsymbol{\tilde N}}}_{12}} = \mathit{\boldsymbol{P}}{\left( {{\mathit{\boldsymbol{N}}_{11}} + \mathit{\boldsymbol{P}}} \right)^{ - 1}}{\mathit{\boldsymbol{N}}_{12}}\\ {{\mathit{\boldsymbol{\tilde N}}}_{21}} = {{\mathit{\boldsymbol{\tilde N}}}_{12}}\\ {{\mathit{\boldsymbol{\tilde N}}}_{22}} = {\mathit{\boldsymbol{N}}_{22}} - {\mathit{\boldsymbol{N}}_{21}}{\left( {{\mathit{\boldsymbol{N}}_{11}} + \mathit{\boldsymbol{P}}} \right)^{ - 1}}{\mathit{\boldsymbol{N}}_{12}}\\ {{\mathit{\boldsymbol{\tilde W}}}_1} = - \mathit{\boldsymbol{P}}{\left( {{\mathit{\boldsymbol{N}}_{11}} + \mathit{\boldsymbol{P}}} \right)^{ - 1}}{\mathit{\boldsymbol{W}}_1}\\ {{\mathit{\boldsymbol{\tilde W}}}_2} = {\mathit{\boldsymbol{W}}_2} - {\mathit{\boldsymbol{N}}_{21}}{\left( {{\mathit{\boldsymbol{N}}_{11}} + \mathit{\boldsymbol{P}}} \right)^{ - 1}}{\mathit{\boldsymbol{W}}_1} \end{array} \right. $$ (9) 需要指出的是，批处理所用的卫星初始状态参数对轨道的约束力会随轨道积分长度增加而逐渐减弱，进而引起轨道精度衰减，因此采用分块递推最小二乘配置方法时，需保持卫星初始状态参数强度，控制轨道解算弧长，定时对卫星状态参数进行向后递推转移，该项工作主要通过轨道状态参数与状态转移矩阵实现。利用式(2)微分方程求解卫星位置时，可通过状态转移矩阵将所求时刻卫星位置转化为求解卫星初始状态参数问题，具体形式为：
$$ \mathit{\boldsymbol{X}} = \mathit{\boldsymbol{\psi }}\left( {t,{t_0}} \right){\mathit{\boldsymbol{X}}_0} $$ (10) 式中，X0为卫星初始状态参数；X为t时刻卫星状态参数；ψ(t, t0)为从t时刻到t0时刻的状态转移矩阵，可通过轨道积分求解。当定轨弧段累计到设计长度时，需将卫星初始状态参数转移递推到下一时刻，叠加新的观测数据以保持卫星定轨弧长不变。轨道状态参数递推如图 2所示。
图 2 轨道初始状态转移递推示意图
Figure 2. Diagram of Transferring and Recursion for Initial Orbit Status
值得注意的是，本文提出的分块递推最小二乘配置方法兼容事后批处理轨道解算方法，既可以实现与滤波方法一致的快速更新，同时又具备了批处理高精度和滤波解快速更新的优点。在实际处理中，当有稳定实时数据源时，采用分块递推最小二乘配置方法，当数据源不能实时获取时，可通过对实时数据累计，分批解算实时轨道。
-
本文采用简化动力学定轨方法，在精密定轨解算过程中，既需要解算与轨道相关的参数，同时也要解算一些其他参数。按参数性质可分为状态参数、过程参数及模糊度参数。状态参数是描述特定时刻状态的参数，一般在一定期间内保持不变；而过程参数是伴随参数解算过程，只在特定时间段内有效的参数；模糊度参数是根据测站卫星实际观测情况在一段时间内有效的参数。在分块递推最小二乘算法实现过程中，实时数据添加时需进行参数初始化，在配置历史先验参数时，要谨慎地连接全局参数，而对于在一定时间内有效的参数也需严格判断，保证实时数据初始化与历史参数一致。
-
卫星初始状态参数包括卫星初始位置参数、初始速度参数及力模型参数；状态参数描述卫星初始参考时刻的状态，在解算弧段内不随时间改变，因此与叠加的实时数据一致，只需将历史数据解算得到的参数和法方程信息配置到实时弧段中。但是需要注意的是，在指定时刻，轨道状态参数具有连续性，实际轨道解算时，人为地将定轨弧段分成两段，而由于受初始轨道精度影响及观测误差的存在，两弧段轨道连接处动力学参数不连续，在进行轨道初始状态参数配置时，需考虑轨道连接处的状态参数不一致问题。本文采用第一个定轨与第二弧段估计得到的轨道连接处轨道状态参数相等作为约束，建立误差方程为：
$$ \left[ {\begin{array}{*{20}{c}} {\Delta {x^{{t_0}}}}\\ {\Delta {v^{{t_0}}}}\\ {\Delta {p^{{t_0}}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {x_0^{{t_k},{t_0}} - x_0^{{t_0},{t_0}}}\\ {v_0^{{t_k},{t_0}} - v_0^{{t_0},{t_0}}}\\ {p_0^{{t_k},{t_0}} - p_0^{{t_0},{t_0}}} \end{array}} \right] + \mathit{\boldsymbol{ \boldsymbol{\varPsi} }}\left( {{t_0},{t_k}} \right)\left[ {\begin{array}{*{20}{c}} {\Delta {x^{{t_k}}}}\\ {\Delta {v^{{t_k}}}}\\ {\Delta {p^{{t_k}}}} \end{array}} \right] $$ (11) 式中，Δxt0、Δvt0、Δpt0分别为两弧段连接引起的轨道初始状态参数误差量；x0tk, t0、v0tk, t0、p0tk, t0和x0t0, t0、v0t0, t0、p0t0, t0分别为第二弧段轨道状态参数初值；Δxtk、Δvtk、Δptk分别为第二弧段计算得到的轨道状态参数改正数；Ψ(t0, tk)为第二弧段到第一弧段状态转移矩阵，可通过轨道积分获取。将式(11)所示的观测方程采用§2.1方法配置到第一弧段法方程中，即可实现轨道状态参数连接。本文保留的历史弧段长度与实时数据弧段长度差异较大，为了避免过多的历史“陈旧”数据累积造成数据过饱和问题[22-23]，同时考虑到新的观测数据对实时轨道解算作用较大，在参数配置过程中，本文对历史弧段状态参数做了适当的松弛处理，以改善方程状态，方法见式(6)至式(8)。
-
定轨过程采用双频无电离层组合观测值，因此过程参数一般包括卫星与接收机钟差参数和对流层延迟参数。随机参数可通过马尔科夫随机过程表示：
$$ {\mathit{\boldsymbol{P}}_{k + 1}} = {M_{k + 1,k}}{\mathit{\boldsymbol{P}}_k} + {\mathit{\boldsymbol{w}}_k} $$ (12) 式中，wk为随机噪声量；当参数满足白噪声特性时，可设Mk+1, k=0，而当参数按随机游走模型估计时，可设Mk+1, k=1。通过式(12)可将轨道连接处的过程参数进行连接，同样可将式(12)作为虚拟观测方程叠加到第一弧段法方程中进行整个弧段综合求解。对于钟差参数采用白噪声逐历元解算，参数只对当前历元有效，失效参数随即消除，因此实时数据解算过程参数重新初始化即可。天顶对流层延迟参数一般采用固定时长的分段常数估计或随机游走方法估计，对于采用分段常数估计时，处理方法与§2.2.1状态参数一致，而对于采用随机游走方法计算对流层参数时，需将第一弧段参数协方差信息作为先验值配置到第二弧段，以保持两段参数的一致性。
-
模糊度配置是关键问题，正确连接历史弧段与实时弧段模糊度参数不仅可以形成完整的法方程，而且能够减少总体解算参数，增加模糊度参数观测数量，提供更多的模糊度固定备选方案。由于模糊度参数与卫星及测站个数有关，甚至数据处理策略不同也会影响模糊度参数的估值及作用时间。针对前后两段模糊度参数配置问题，首先根据弧段连接处观测值的数据特性及模糊度有效时间判断能否进行模糊度参数连接。判断模糊度连接特性时，本文采用Melbourne-Wübbena(MW)组合和无几何距离组合进行分析。
假设两个待连接模糊度观测弧段的MW组合平均值和方差分别为N1、N2和σ1、σ2，则给定判断条件：
$$ \left\{ \begin{array}{l} \left| {{{\bar N}_1} - {{\bar N}_2}} \right| < 3\sigma \\ \sigma = \sqrt {\sigma _1^2 + \sigma _2^2} \end{array} \right. $$ (13) 若满足上述条件，则认为两弧段连接处不存在周跳或粗差，可以进行配置连接。但当两个载波的观测值发生大小相等的周跳时，采用MW组合不能检测出该类周跳[24-25]，需要增加无距离组合的连续性检验。
首先对历史弧段中待连接模糊度观测值的无几何距离组合观测值Lg进行多项式拟合，并外推到第二个弧段的观测历元，然后将外推值与第二个弧段的真实观测值求差。由于无几何距离组合观测值主要受电离层延迟误差的影响[26]，因此所求的差值反映了电离层和周跳的变化。由于短时间内电离层变化缓慢[27]，对多个差值求平均，给定以下条件，即可判断在连接处是否发生周跳或较大粗差：
$$ \left| {{{\bar L}_g} - {{\bar P}_g}} \right| < \varepsilon $$ (14) 式中，Lg为一段时间无几何距离组合平均值；Pg为相应时间外推值的平均值。当同时满足式(13)和式(14)的条件时，则通过第一弧段模糊度参数对第二弧段模糊度参数进行配置，否则需新建立一个模糊度参数，并重新初始化。由于历史弧段参数解算精度较高，能够较准确地给出模糊度参数估值及先验精度，因此可通过历史弧段获取的模糊度参数协方差信息对实时弧段模糊度参数配置提供较好的先验信息约束，进一步能够提高实时弧段的模糊度固定效率与正确率。精密定轨时采用无电离层组合观测值，可将模糊度分解为宽巷及窄巷组合形式，但是由于接收机端和卫星端均有未校准的硬件延迟影响，解算的模糊度参数并不具有整数特性[28]。在进行第二段数据处理时，首先对观测数据在测站之间作单差，用以消除卫星端的硬件延迟，然后对卫星间数据作差分，可进一步消除接收机端硬件延迟，形成双差观测值后的模糊度理论上具备整数特性，将载波相位无电离层组合模糊度分解为宽巷与窄巷组合形式：
$$ {N_{{\rm{LC}}}} = \frac{{{f_1}{f_2}}}{{f_1^2 - f_2^2}}\left[ {\left( {{N_1} - {N_2}} \right) + \frac{{{f_1} - {f_2}}}{{{f_2}}}{N_1}} \right] $$ (15) 式中，NLC为LC组合模糊度；f1、f2为载波相位频率；N1-N2为宽巷组合模糊度；N1为窄巷组合模糊度。由于宽巷模糊度波长较长，一般可通过对多历元MW组合观测值平滑后取整求得，宽巷模糊度固定后，利用解算得到的LC组合模糊度可推导出窄巷模糊度估值，再通过判别函数取整法依次进行固定，具体方法可参见文献[29-31]；待双差模糊度参数固定后，将得到的双差模糊度通过构建虚拟观测方程方法映射到非差数据处理法方程中，以增加解的强度。
-
为了提高参数解算的效率，本文采用GPS/ GLONASS和GPS/Galileo/BDS两个线程进行四系统精密轨道解算。坐标参考框架强约束到IGS发布的IGb14框架系统。通过递推最小二乘配置法将历史弧段解算的轨道状态参数、模糊度参数约束到实时弧段参数解算中，具体解算流程如图 3所示。
图 3 定轨流程示意图
Figure 3. Workflow of Orbit Determination
实时轨道解算时，考虑到较长的解算弧长可以改善轨道预报精度，本文采用的历史累计弧段长度为33 h，实时弧段部分设置为3 h，轨道解算部分长度保持为36 h。首次累计解算33 h弧长历史数据耗时约为90 min，计算3 h实时轨道弧长部分耗时约为10 min，考虑到实时数据流及事后多系统小时观测文件传到本地所需时间，可预留20 min用于实时数据等待时间，本文解算四系统实时轨道产品为30 min延迟后更新。解算延迟与发布如图 4所示。
图 4 实时轨道解算延迟与发布示意图
Figure 4. Illustration of Delay and Updating of Real-Time Orbit Solutions
目前已实现整点之后20 min开始轨道解算，10 min即可完成实时数据部分轨道计算并发布产品，四系统产品每小时更新并延迟30 min发布，用户可用部分为预报轨道的0.5~1.5 h部分。本文定轨采用的模型改正项及定轨参数设置如表 1所示。
表 1 模型改正项及定轨参数设置
Table 1. Model Correction Term and Orbit Parameters
项目 模型 解算系统 GPS、GLONASS、BDS、Galileo 解算弧段 历史弧长33 h，实时弧长3 h 观测值组合 非差无电离层组合 信号频率 L1&L2、G1&G2、B1&B2、E1&E5a 采样间隔 300 s 截止高度角 7°(GPS & GLONASS); 10°(BDS & Galileo) 光压模型 ECOM 5 轨道积分步长 60 s 轨道摄动项 非球形地心引力、多体摄动、固体潮、 海潮、地球返照辐射、太阳光压 天线相位中心 GPS & GLONASS: igs14.atx； BDS & Galileo: WHU[14] -
目前IGS组织提供全球均匀分布的测站约有450个，其中每小时数据更新测站超过300个，且同时提供多系统服务的测站约200个，其中大部分测站包含Galileo数据，但可用于每小时更新测站且包含北斗数据的测站数量约为40个，测站全球分布不均匀，因此本文加入了iGMAS、陆态网及澳洲连续运行参考网数据，用以在局部区域进行观测数据的补充。本文定轨选取120个多系统实时测站，测站分布图如图 5所示，其中所有测站均可进行GPS和GLONASS卫星观测，BDS和Galileo可用测站个数分别为85个和81个。
图 5 四系统定轨测站分布图
Figure 5. Distribution of Quad-Constellation Orbit Determination Stations
笔者通过对武汉大学的精密定轨软件PANDA进行二次开发，增加了实时轨道解算分块最小二乘递推功能。目前已经分别提供了GPS、GLONASS、Galileo、BDS四系统0.5 h更新轨道、3 h更新轨道及6 h更新轨道，随着多GNSS系统的建设发展，个别IGS MGEX(Multi-GNSS Experiment)分析中心也提供了多系统实时轨道。为了比较相关算法及性能的改进，本文对武汉大学发布的3种轨道、GFZ发布的多系统超快速实时轨道及IGS超快速实时轨道用户可用部分精度进行了分析比较。由于实时数据等待及数据解算耗时，实时用户可用部分往往要推迟一段时间，而由于实时轨道精度随预报时长推移而衰减，因此缩短轨道更新时间成为了提高轨道精度的一个重要因素。表 2列出了武汉大学、IGS及GFZ目前已经发布的产品更新、延迟及解算方案，同时还列出了实时用户可用部分轨道预报时间。其中WHU 1 h更新轨道采用递推最小二乘配置方法解算，WHU 3 h、6 h及GFZ发布轨道均为最小二乘批处理解算，而IGS发布轨道为通过对多个分析中心发布轨道综合而得。本文在进行轨道精度统计时，GPS轨道以IGS发布的最终轨道产品作为参考，而GLONASS、Galileo及北斗轨道以CODE发布的事后多系统轨道作为参考，轨道精度统计时段为2018年年积日231~265天，共35天数据。为了表示方便，对各类实时轨道产品定义如下：案例1:WHU 1 h更新延迟30 min发布四系统轨道产品；案例2:WHU 3 h更新延迟2 h发布四系统轨道产品；案例3:WHU 6 h更新延迟2 h发布四系统轨道产品；案例4:GFZ 3 h更新延迟2 h发布四系统轨道产品；案例5:IGS 6 h更新延迟3 h发布四系统轨道产品。
表 2 实时轨道延迟与更新方案
Table 2. Schemes of Ultra-Rapid Orbit Delay and Updating
机构 更新方
案/h延迟时
间/h用户可用部
分/h实时轨道解
算方法WHU 1 0.5 0.5~1.5 递推配置 WHU 3 2 3~6 批处理解算 WHU 6 2 2~8 批处理解算 GFZ 3 2 3~6 批处理解算 IGS 6 3 3~9 ACs综合解 5种案例实时轨道的精度统计如表 3所示，其中CODE提供的多系统参考轨道未提供北斗GEO卫星，故本文只给出北斗IGSO和MEO卫星实时用户可用部分精度。另外，由于IGS超快速综合轨道仅发布GPS系统产品，故案例5中其他3个系统未有比较结果。本文统计的比较结果为相关轨道与参考轨道间的三维均方根误差(root mean square, RMS)，其中1 h、3 h、6 h更新轨道精度统计量分别为840、280和140个。
表 3 实时轨道产品精度统计/mm
Table 3. Precisions of Real-Time Orbit Products /mm
系统 案例1 案例2 案例3 案例4 案例5 GPS 28 42 46 39 42 GLONASS 85 110 117 101 - Galileo 50 67 74 62 - BDS 115 155 198 153 - 图 6为上述5种案例对应的每个系统中所有卫星参与统计的3D RMS值，可以看出，与其他案例相比，使用本文方法计算的案例1实时轨道用户可用部分精度均有较大提高，GPS、GLONASS、Galileo及BDS系统实时用户可用部分三维精度分别为2.8 cm、8.5 cm、5.0 cm及11.5 cm。其中，针对WHU发布的GPS精度，相比于3 h和6 h更新轨道，本文方法分别提高33.3%和39.1%；同样地，GLONASS精度分别提高22.7%和27.4%；Galileo精度分别提高25.3%和32.4%；BDS IGSO/MEO精度分别提高25.8%和41.9%。案例1实时轨道精度也分别高于GFZ和IGS发布的超快速产品精度。
图 6 四系统轨道5种案例精度统计
Figure 6. Statistical Results of Five Cases for Quad-Constellation Orbit Determination
图 7为各系统每颗卫星用户可用部分精度的统计结果，其中对于GPS，分别统计了武汉大学实时轨道、GFZ超快速实时轨道和IGS超快速实时轨道用户可用部分精度，而对于其他3个系统，则分别比较了武汉大学实时轨道和GFZ超快速实时轨道用户可用部分精度。考虑到实时轨道精度随轨道预报时间的增加而逐渐衰减，武汉大学实时轨道解算时采用了分块递推最小二乘配置算法，提高了轨道解算效率，可使轨道更新延迟时间小于0.5 h。GFZ实时轨道解算方法为最小二乘批处理，待估参数较多，导致轨道解算效率低，其轨道更新延迟时间约为2 h。而IGS超快速实时轨道需要综合各个超快速分析中心的产品，其产品延迟更新时间约为3 h。在实时应用中，用户关心的是轨道更新后的可用部分精度，因此在先确保轨道确定精度的前提下，缩短轨道更新延迟显得尤为重要。从图 7中统计结果可以看出，通过高效的参数估计方法和策略，极大地缩短了实时轨道预报时间，武汉大学解算的实时轨道用户可用部分精度在4个系统中均为最优。不难发现，本文提出的方法可在一定程度上提高各系统各类卫星实时用户可用部分精度。GPS卫星轨道3D RMS大部分优于30 mm，但G05、G10、G22卫星精度稍差，可能是由于在数据统计期间，这3颗卫星恰好处于地影期，导致轨道预报精度降低。不过值得注意的是, 得益于较高的轨道更新频率，武汉大学实时轨道仍可保持较高的轨道精度。GLONASS卫星轨道3D RMS一般优于90 mm。BDS卫星由于特殊的轨道设计，IGSO卫星轨道高度和轨道周期与其他3个系统均有较大差异，相关的卫星动力学模型、姿态控制及积分理论仍在研究中，尚没有一套成熟稳健的算法供各个分析中心统一采用，因此目前轨道解算及预报精度稍差，轨道3D RMS约为120 mm。Galileo卫星在轨正常运行卫星较多，可以较稳健地进行模糊度固定，且卫星未有较大的动力学模型及姿态模型误差，实时轨道3D RMS约为50 mm。
图 7 2018年231~237天各系统卫星平均RMS统计
Figure 7. Average RMS for Quad-Constellation Satellites During DOY 231 to 237 in 2018
目前, 武汉大学IGS数据中心FTP服务器已提供计算的1 h更新延迟30 min发布四系统实时精密轨道产品，同时，也已提供基于该产品搭建的实时产品服务系统。实时产品服务系统采用计算的实时轨道，利用滤波方法计算四系统卫星实时钟差，并通过Ntrip协议发布精密轨道、钟差产品相对于广播星历产品的改正数，可为广域用户提供实时精密轨道钟差改正，目前相关服务可从实时产品服务系统获取，IP地址及端口为59.175.223.165 :2101。
-
本文提出一种基于分块递推最小二乘配置的实时多GNSS系统精密轨道解算方法，该方法原理上可兼容最小二乘批处理和实时滤波两种方法，降低了实时轨道解算对实时数据的依赖，同时，在进行数据质量控制时采用了严格的后处理方法，较滤波方法数据处理更加有效可靠。经统计，该方法解算的GPS轨道3D RMS为2.8 cm，精度优于IGS超快速产品；GLONASS轨道3D RMS为8.5 cm，Galileo、BDS轨道3D RMS分别为5.0 cm和11.5 cm，轨道精度均优于GFZ发布的超快速产品。相较于武汉大学之前发布的3 h和6 h轨道，GPS轨道精度提高了33.3%和39.1%；GLONASS轨道精度提高22.7%和27.4%；BDS IGSO/MEO轨道精度提高25.8%和41.9%；Galileo轨道精度提高25.3%和32.4%。实验证明，本文方法能够有效地提高参数解算效率，改善轨道更新频率。目前，采用该方法解算的多GNSS系统每小时更新延迟30 min发布的轨道产品已经发布到武汉大学IGS数据中心，实时用户可通过FTP服务获取相关产品。同时，基于Ntrip协议的实时产品服务系统也正式上线，可为实时用户提供广域精密单点定位服务。
-
摘要: 多系统全球导航卫星系统（Global Navigation Satellite System，GNSS）精密轨道确定及其预报是实现高精度实时精密定位的前提。针对多GNSS系统超快速轨道解算时效性及轨道预报精度随时间下降的问题，提出一种基于分块递推最小二乘配置方法，该方法通过对动力学和几何学待估参数松弛、连接以及轨道状态参数转移递推，能够同时兼容事后及实时滤波定轨方法。该方法能够有效地提高多GNSS系统轨道解算效率，缩短实时轨道更新时间。基于全球实测数据验证了该方法的可靠性和有效性，轨道精度优于国际GNSS服务组织发布的GPS超快速轨道及德国地学研究中心发布的超快速轨道，实验结果表明，采用该方法，GPS/GLONASS/Galileo/BDS四系统120个地面测站精密定轨可以实现1 h更新，延迟30 min发布，统计GPS/GLONASS/Galileo/BDS实时轨道可用部分3D均方根分别为2.8 cm、8.5 cm、5.0 cm及11.5 cm（IGSO/MEO）。目前，1 h更新多GNSS系统轨道及实时产品服务系统已业务化发布，较之前发布的3 h更新及6 h更新轨道分别有20%~40%的精度提升。Abstract: Multi-GNSS (Global Navigation Satellite System) precise orbit determination and prediction is the prerequisite of high-precision real-time precise positioning. To improve the efficiency of ultra-rapid orbit solution and prediction precision in multi-GNSS, this paper proposes a collocation method based on block recursive least squares. The proposed method can combine post estimation with real-time filtering orbit determination by ① relaxing and connecting kinetics and geometry parameters, and ② transferring and recurring orbit status parameters. In addition, this method increases the efficiency of GNSS orbit solution and shortens the orbit updating time. Furthermore, experimental result verifies the reliability and effectiveness of our proposed method. Our orbit precisions are superior to the precisions of ultra-rapid GPS orbits provided by international GNSS service (IGS) and German Research Centre for Geosciences. The results show that, for quad-constellation (GPS/GLONASS/Galileo/BDS) precise orbit determination by 120 ground observation stations, our method is able to realize one-hour updating and provides service with 30 minutes delay. The 3D root mean squares of available parts in GPS/GLONASS/Galileo/BDS real-time orbit are 2.8 cm, 8.5 cm, 5.0 cm and 11.5 cm(IGSO/MEO), respectively. At present, one-hour updating for multi-GNSS orbits and real-time service system have been provided, and compared with previous three-hour and six-hour updating orbit, the precision has been improved by 20% to 40%.
-
图 2 轨道初始状态转移递推示意图
Figure 2. Diagram of Transferring and Recursion for Initial Orbit Status
图 4 实时轨道解算延迟与发布示意图
Figure 4. Illustration of Delay and Updating of Real-Time Orbit Solutions
图 5 四系统定轨测站分布图
Figure 5. Distribution of Quad-Constellation Orbit Determination Stations
图 6 四系统轨道5种案例精度统计
Figure 6. Statistical Results of Five Cases for Quad-Constellation Orbit Determination
图 7 2018年231~237天各系统卫星平均RMS统计
Figure 7. Average RMS for Quad-Constellation Satellites During DOY 231 to 237 in 2018
表 1 模型改正项及定轨参数设置
Table 1. Model Correction Term and Orbit Parameters
项目 模型 解算系统 GPS、GLONASS、BDS、Galileo 解算弧段 历史弧长33 h，实时弧长3 h 观测值组合 非差无电离层组合 信号频率 L1&L2、G1&G2、B1&B2、E1&E5a 采样间隔 300 s 截止高度角 7°(GPS & GLONASS); 10°(BDS & Galileo) 光压模型 ECOM 5 轨道积分步长 60 s 轨道摄动项 非球形地心引力、多体摄动、固体潮、 海潮、地球返照辐射、太阳光压 天线相位中心 GPS & GLONASS: igs14.atx； BDS & Galileo: WHU[14] 
下载: 导出CSV
表 2 实时轨道延迟与更新方案
Table 2. Schemes of Ultra-Rapid Orbit Delay and Updating
机构 更新方
案/h延迟时
间/h用户可用部
分/h实时轨道解
算方法WHU 1 0.5 0.5~1.5 递推配置 WHU 3 2 3~6 批处理解算 WHU 6 2 2~8 批处理解算 GFZ 3 2 3~6 批处理解算 IGS 6 3 3~9 ACs综合解
下载: 导出CSV
表 3 实时轨道产品精度统计/mm
Table 3. Precisions of Real-Time Orbit Products /mm
系统 案例1 案例2 案例3 案例4 案例5 GPS 28 42 46 39 42 GLONASS 85 110 117 101 - Galileo 50 67 74 62 - BDS 115 155 198 153 -
下载: 导出CSV
-
[1] Dow J M, Neilan R E, Rizos C. The International GNSS Service in a Changing Landscape of Global Navigation Satellite Systems[J]. Journal of Geodesy, 2009, 83(3-4):191-198 doi: 10.1007/s00190-008-0300-3 [2] Cai H, Chen K, Xu T, et al. The iGMAS Combined Products and the Analysis of Their Consistency[C]//China Satellite Navigation Conference (CSNC) 2015 Proceedings: Volume Ⅲ. Berlin, Heidelberg: Springer, 2015, 342: 213-226 [3] Montenbruck O, Steigenberger P, Prange L, et al. The Multi-GNSS Experiment (MGEX) of the International GNSS Service (IGS) -Achievements, Prospects and Challenges[J]. Advances in Space Research, 2017, 59(7):1671-1697 doi: 10.1016/j.asr.2017.01.011 [4] Beutler G, Rothacher M, Schaer S, et al. The International GPS Service (IGS):An Interdisciplinary Service in Support of Earth Sciences[J]. Advances in Space Research, 1999, 23(4):631-653 doi: 10.1016/S0273-1177(99)00160-X [5] Brockmann E. Combination of Solutions for Geode-tic and Geodynamic Applications of the Global Positioning System (GPS)[J]. Geod-Geophys Arb Schweiz, 1997, 55:55 [6] Dong D, Herring T A, King R W. Estimating Regional Deformation from a Combination of Space and Terrestrial Geodetic Data[J]. Journal of Geodesy, 1998, 72(4):200-214 doi: 10.1007/s001900050161 [7] Andersen P H. Multi-level Arc Combination with Stochastic Parameters[J]. Journal of Geodesy, 2000, 74(7):531-551 doi: 10.1007/s001900000115 [8] 姚宜斌. GPS轨道合成原理及其实现[J].武汉大学学报·信息科学版, 2007, 32(6):510-514 http://ch.whu.edu.cn/CN/abstract/abstract1906.shtml Yao Yibin. Theory and Realization of GPS Orbit Integration[J]. Geomatics and Information Science of Wuhan University, 2007, 32(6):510-514 http://ch.whu.edu.cn/CN/abstract/abstract1906.shtml [9] 楼益栋, 施闯, 葛茂荣, 等. GPS卫星实时精密定轨及初步结果分析[J].武汉大学学报·信息科学版, 2008, 33(8):815-817 http://ch.whu.edu.cn/CN/abstract/abstract1672.shtml Lou Yidong, Shi Chuang, Ge Maorong, et al. GPS Real Time Orbit Determination and Initial Results Analysis[J].Geomatics and Information Science of Wuhan University, 2008, 33(8):815-817 http://ch.whu.edu.cn/CN/abstract/abstract1672.shtml [10] Lutz S, Meindl M, Steigenberger P, et al. Impact of the Arc Length on GNSS Analysis Results[J]. Journal of Geodesy, 2016, 90(4):1-14 doi: 10.1007/s00190-015-0878-1 [11] 刘伟平, 郝金明, 于合理, 等.北斗卫星多天轨道合成方法及其精度分析[J].测绘学报, 2016, 45(10):1157-1164 doi: 10.11947/j.AGCS.2016.20160055 Liu Weiping, Hao Jinming, Yu Heli, et al. Solution Method and Precision Analysis of Multi-days Orbit Combination of BeiDou Satellites[J].Acta Geodaetica et Cartographica Sinica, 2016, 45(10):1157-1164 doi: 10.11947/j.AGCS.2016.20160055 [12] Kazmierski K, Soŝnica K, Hadas T. Quality Assessment of Multi-GNSS Orbits and Clocks for Real-Time Precise Point Positioning[J]. GPS Solutions, 2018, 22(1):11 doi: 10.1007/s10291-017-0678-6 [13] Li X, Chen X, Ge M, et al. Improving Multi-GNSS Ultra-Rapid Orbit Determination for Real-Time Precise Point Positioning[J]. Journal of Geodesy, 2018, 92(3):1-20 doi: 10.1007/s00190-018-1138-y [14] Guo J, Xu X, Zhao Q, et al. Precise Orbit Determination for Quad-Constellation Satellites at Wuhan University:Strategy, Result Validation, and Comparison[J]. Journal of Geodesy, 2016, 90(2):1-17 doi: 10.1007/s00190-015-0862-9 [15] 戴小蕾.基于平方根信息滤波的GNSS导航卫星实时精密定轨理论与方法[D].武汉: 武汉大学, 2016 http://cdmd.cnki.com.cn/Article/CDMD-10486-1016124266.htm Dai Xiaolei. Real-Time Precise GNSS Satellite Orbit Determination Using the SRIF Method: Theory and Implementation[D]. Wuhan: Wuhan University, 2016 http://cdmd.cnki.com.cn/Article/CDMD-10486-1016124266.htm [16] Pavlis N K, Holmes S A, Kenyon S C, et al. The Development and Evaluation of Earth Gravitational Model (EGM2008)[J]. Journal of Geophysical Research:Solid Earth, 2013, 118(5):2633 doi: 10.1002/jgrb.50167 [17] Beutler G, Brockmann E, Gurtner W, et al. Extended Orbit Modeling Techniques at the CODE Processing Center of the International GPS Service for Geodynamics (IGS):Theory and Initial Results[J]. Manuscripta Geodaetica, 1994, 19:367-386 [18] Petit G, Luzum B. IERS Conventions (2010)[R]. IERS Technical Note No. 36, IERS Convention Centre, 2010 [19] Montenbruck O, Gill E, Lutze F. Satellite Orbits:Models, Methods, and Applications[J]. Applied Mechanics Reviews, 2002, 55(2):2504-2510 http://d.old.wanfangdata.com.cn/Periodical/kzllyyy201004003 [20] 赵齐乐. GPS导航星座及低轨卫星的精密定轨理论和软件研究[D].武汉: 武汉大学, 2005 http://www.cqvip.com/QK/71135X/201107/15333621.html Zhao Qile. Research on Precision Orbit Determination Theory and Software of Both GPS Naviagation Constellation and LEO Satellites[D]. Wuhan: Wuhan University, 2005 http://www.cqvip.com/QK/71135X/201107/15333621.html [21] 葛茂荣. GPS卫星精密定轨理论及软件研究[D].武汉: 武汉大学, 1995 http://www.wanfangdata.com.cn/details/detail.do?_type=degree&id=Y219152 Ge Maorong. Study on Theory and Software of Precise Oribt Determination for GPS Satellite[D]. Wuhan: Wuhan University, 1995 http://www.wanfangdata.com.cn/details/detail.do?_type=degree&id=Y219152 [22] Duan Bingbing, Chen Junping. Extended Filter for Real-Time Multi-GNSS Orbit Determination[C]. IGS Workshop, Sydney, NSW, Australia, 2016 [23] Mervart L, Weber G. Real-Time Combination of GNSS Orbit and Clock Correction Streams Using a Kalman Filter Approach[C]. The 24th International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GNSS 2011), Portland, USA, 2011 [24] Cai C, Liu Z, Xia P, et al. Cycle Slip Detection and Repair for Undifferenced GPS Observations Under High Ionospheric Activity[J]. GPS Solutions, 2013, 17(2):247-260 doi: 10.1007/s10291-012-0275-7 [25] Zhao Q, Sun B, Dai Z, et al. Real-Time Detection and Repair of Cycle Slips in Triple-Frequency GNSS Measurements[J]. GPS Solutions, 2015, 19(3):381-391 doi: 10.1007/s10291-014-0396-2 [26] Feng Y. GNSS Three Carrier Ambiguity Resolution Using Ionosphere-Reduced Virtual Signals[J]. Journal of Geodesy, 2008, 82(12):847-862 doi: 10.1007/s00190-008-0209-x [27] Fritsche M, Dietrich R, Knöfel C, et al. Impact of Higher-Order Ionospheric Terms on GPS Estimates[J]. Geophysical Research Letters, 2005, 32(23):113-133 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=0e51aded8bdc164d04c84676d1de3686 [28] Li X, Ge M, Zhang H, et al. A Method for Improving Uncalibrated Phase Delay Estimation and Ambiguity-Fixing in Real-Time Precise Point Positioning[J]. Journal of Geodesy, 2013, 87(5):405-416 doi: 10.1007/s00190-013-0611-x [29] Dong D N, Bock Y. Global Positioning System Network Analysis with Phase Ambiguity Resolution Applied to Crustal Deformation Studies in California[J]. Journal of Geophysical Research:Atmospheres, 1989, 94(B4):3949-3966 doi: 10.1029/JB094iB04p03949 [30] Ge M, Gendt G, Dick G, et al. Improving Carrier-Phase Ambiguity Resolution in Global GPS Network Solutions[J]. Journal of Geodesy, 2005, 79(1):103-110 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=5202bb66ea9980502ebd35f7358148c4 [31] Ge M, Gendt G, Rothacher M, et al. Resolution of GPS Carrier-Phase Ambiguities in Precise Point Positioning (PPP) with Daily Observations[J]. Journal of Geodesy, 2008, 82(7):389-399 doi: 10.1007/s00190-007-0187-4 -
点击查看大图
图(7) / 表(3)
计量
- 文章访问数: 1558
- HTML全文浏览量: 109
- PDF下载量: 324
- 被引次数: 0
