对于低轨卫星星载GPS双频接收机，采用GPS双频无电离层（ionospheric free，IF）组合伪距和载波相位观测值消除电离层延迟误差，且不考虑对流层延迟等误差，其观测方程可简化为：

式中，$P_{\mathrm{I}\mathrm{F}}$和$L_{\mathrm{I}\mathrm{F}}$分别为GPS双频伪距和载波相位无电离层组合观测值，单位：m；$\rho$为星载接收机与GPS卫星之间的几何距离；$\mathrm{\delta }{t}_r$和$\mathrm{\delta }{t}^s$分别为星载接收机和GPS卫星的钟差参数；c为真空中的光速；$B_{\mathrm{I}\mathrm{F}}$为组合载波相位的模糊度；$M_P$和$M_L$分别为组合伪距和载波相位的多路径误差；${\epsilon }_P$和${\epsilon }_L$分别为组合伪距和载波相位的观测噪声。

使用与GPS广播星历参数中星历/时钟数据期号相符且最近的实时改正数，用户可得到GPS卫星的精密轨道和钟差。假设当前历元信号发射t时刻由GPS广播星历计算的轨道和钟差分别为${\mathit{r}}_{\mathrm{b}\mathrm{r}\mathrm{d}\mathrm{c}}^s$和$\mathrm{d}{t}_{\mathrm{b}\mathrm{r}\mathrm{d}\mathrm{c}}^s$，则GPS卫星精密轨道${\mathit{r}}^s$和精密钟差$\mathrm{d}{t}^s$的计算如下：

式中，$\mathit{r}\mathbf{、}\mathit{v}$分别为由GPS广播星历计算的在当前历元信号发射t时刻WGS84坐标系下GPS卫星位置、速度；${\mathit{e}}_r\mathbf{、}{\mathit{e}}_{\mathit{a}}\mathrm{、}{\mathit{e}}_c$分别为在GPS卫星轨道径向（radial，R）、切向（along，A）、法向（cross，C）上的单位向量；$\mathrm{\delta }\mathit{o}$为GPS卫星在轨道坐标系下RAC 3个方向上的改正量；$\mathrm{\delta }c$为GPS卫星钟差改正量；c为真空中的光速；$(\mathrm{\delta }{o}_r,\mathrm{\delta }{o}_a,\mathrm{\delta }{o}_c)$和$(\mathrm{d}\mathrm{\delta }{o}_r,\mathrm{d}\mathrm{\delta }{o}_a,\mathrm{d}\mathrm{\delta }{o}_c)$分别为$t_{\mathrm{0}}$参考时刻的GPS卫星在轨道坐标系下RAC 3个方向的轨道改正值及其变化率；$(c_{\mathrm{0}},c_{\mathrm{1}},c_{\mathrm{2}})$为GPS卫星钟差改正的系数。这些轨道与钟差改正数信息可直接从实时产品中得到，在短时接收中断情况下，用最近时刻的改正数来外推当前时刻的改正数。

由于实时改正数产品是由IGS分析中心等机构根据精密定轨得到的轨道与钟差参数经短时间预报产生。当前，由于轨道动力学模型精度较高，短时轨道预报精度可达到厘米级。但是GPS星座使用不同类型的原子钟，卫星钟差的预报精度存在一定差异。另外，在极区以及接收中断等情况下，直接使用最近时刻的改正数来外推当前时刻的改正数，易导致改正数的误差特性发生较大的变化。由实时改正数及其短时外推改正数计算的GPS卫星轨道与钟差，与真实的卫星轨道和钟差仍存在一定差异，这种差异在GPS信号传播路径（location of signal，LOS）上引起的星历综合误差$E_{\mathrm{L}\mathrm{O}\mathrm{S}}$可表示为^[5，12]：

式中，${\rho }^{\mathrm{*}}$和$c\mathrm{\delta }{t^s}^{\mathrm{*}}$分别为真实的低轨卫星与GPS卫星之间的几何距离和GPS钟差；$\rho$和$c\mathrm{\delta }{t}^s$分别为由GPS广播星历和实时改正数计算的几何距离和GPS钟差；c为真空中的光速；$\mathrm{d}\rho$和$c·\mathrm{d}\mathrm{\delta }{t}^s$分别为由星历误差引起的在GPS信号传播路径上的几何距离误差与GPS钟误差。

由于$E_{\mathrm{L}\mathrm{O}\mathrm{S}}$参数和式（1）中模糊度参数$B_{\mathrm{I}\mathrm{F}}$的系数一致且均未知，无法准确分离，同时，为了减少后续实时定轨滤波模型中待估参数的数量，提高实时定轨计算效率，本文将这两个参数合并为一个参数，称为伪模糊度参数^[5，12]，即$A=B_{\mathrm{I}\mathrm{F}}+E_{\mathrm{L}\mathrm{O}\mathrm{S}}$。因此，可将观测方程改写为：

低轨卫星围绕地球运转的动力学方程可表示为：

式中，${\mathit{a}}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}$为低轨卫星所受的总加速度；${\mathit{a}}_g$为由保守力（包括地球中心引力和非球形引力、N体引力、地球固体潮汐和海洋潮汐摄动力等）引起的加速度；${\mathit{a}}_{ng}$为由非保守力（包括大气阻力、太阳光压力等）引起的加速度；${\mathit{a}}_w$为人为引入的经验加速度，用于补偿无法模型化或错误模型的微小摄动力的影响，采用一阶高斯-马尔可夫随机模型对RAC 3个方向进行动力学模型补偿^[5]。

以GPS双频无电离层组合伪距和载波相位作为主要观测值，结合低轨卫星简化动力学模型，用扩展卡尔曼滤波估计卫星的位置、速度等状态参数，实时定轨滤波模型中的状态方程和观测方程如下：

式中，${\mathit{X}}_k$和${\mathit{Z}}_k^{}$分别为扩展卡尔曼滤波状态量和观测量；${\mathit{\Phi }}_{k,k-\mathrm{1}}$和${\mathit{H}}_k$分别为状态转移矩阵和观测矩阵；${\mathit{W}}_{k-\mathrm{1}}$和${\mathit{V}}_k$为系统噪声和观测噪声。其中，待估滤波状态量选取为：

式中，${\mathit{Y}}_{\mathrm{6}\times \mathrm{1}}=\left[\begin{array}{cc}\mathit{r}& \mathit{v}\end{array}\right]$为$k$时刻低轨卫星轨道参数（包括三维位置$\mathit{r}$和速度v）；${\mathit{T}}_{\mathrm{2}\times \mathrm{1}}=\left[\begin{array}{cc}\mathrm{\delta }{t}_r& \mathrm{d}\mathrm{\delta }{t}_r\end{array}\right]$为GPS接收机钟差$\mathrm{\delta }{t}_r$和钟漂$\mathrm{d}\mathrm{\delta }{t}_r$参数；${\mathit{P}}_{\mathrm{5}\times \mathrm{1}}=\left[\begin{array}{ccc}{c}_d& c_r& {\mathit{a}}_w\end{array}\right]$为动力学模型相关参数（包括大气阻力系数$c_d$、太阳光压系数$c_r$、RAC 3个方向的经验加速度${\mathit{a}}_w$）；${\mathit{A}}_{n\times \mathrm{1}}=\left[\begin{array}{cc}{A}_{\mathrm{1}}& \begin{array}{ccc}{A}_{\mathrm{2}}& \cdots & A_n\end{array}\end{array}\right]$为所观测到的n颗GPS卫星待估的伪模糊度参数。

在参数估计时，模糊度$B_{\mathrm{I}\mathrm{F}}$为常量，伪模糊度参数$A$的随机模型与$E_{\mathrm{L}\mathrm{O}\mathrm{S}}$参数紧密相关，需要根据$E_{\mathrm{L}\mathrm{O}\mathrm{S}}$的变化特性来合理建模。通过伪模糊度参数的实时准确估计，减小$E_{\mathrm{L}\mathrm{O}\mathrm{S}}$对实时定轨的影响，可提高实时定轨精度。

本文以IGS最终精密星历和精密钟差（30 s）产品为参考，分析IGS实时改正数产品的精度及其变化特性。利用BNC（BKG Ntrip Client）软件接收了CLK93（采样率为5 s）挂载点在2019-04-08—2019-04-14，即年积日（day of year，DOY）第98~104天，的IGS-RTS数据流。使用式（2）~式（7），以30 s为间隔计算CLK93实时精密轨道与钟差结果，并与IGS精密星历和精密钟差产品计算结果作比较，在计算过程中，要将IGS精密轨道从卫星质心改正到天线相位中心。

将两者的轨道差值转换到R、A、C方向，并统计各方向的均方根误差（root mean square error，RMSE）。为消除卫星钟差的系统偏差影响，本文选取G01卫星作为基准星，采用求二次差的方法评估实时钟差的精度^[13]，统计钟差二次差的均值（MEAN）、RMSE、标准差（STD）。

图 1为DOY第98~104天CLK93实时产品中所有卫星轨道和钟差每天的精度统计。由图 1可以看出，CLK93实时产品在RAC方向的轨道精度优于5 cm，钟差精度为0.2~0.5 ns。

图  1  CLK93产品的轨道与钟差精度

Figure 1.  Orbit and Clock Accuracies of CLK93 Products

图 2给出了DOY第101天每颗GPS卫星轨道和钟差的统计结果，其中G04没有轨道和钟差改正信息。由图 2可以看出，不同卫星之间的轨道误差和卫星钟同步误差都存在一定的不一致性，其中，G08卫星轨道误差超过了10 cm，G03、G18、G21等卫星钟同步误差都超过了1 ns。实时产品误差将直接影响低轨卫星实时精密定轨的轨道精度，在实时定轨算法中，需要加以建模来减小其影响。

图  2  各GPS卫星的轨道与钟差精度

Figure 2.  Orbit and Clock Accuracies for Each GPS Satellite

考虑到GEO覆盖盲区，本文假定低轨卫星星下点纬度大于65°的区域为极区范围，对于轨道高度约500 km的极轨卫星，通过极区所需时间大约为13 min。在极区改正数接收中断时，直接使用最后可用时刻的改正数来外推得到。

以IGS精密星历和精密钟差（30 s）产品的计算结果代替式（4）的真实轨道与钟差，使用GPS广播星历和CLK93实时产品，计算GRACE C卫星观测的所有GPS卫星的$\mathrm{d}\rho$、$c·\mathrm{d}\mathrm{\delta }{t}^s$、$E_{\mathrm{L}\mathrm{O}\mathrm{S}}$和$\mathrm{d}{E}_{\mathrm{L}\mathrm{O}\mathrm{S}}$。图 3给出了DOY第101天连续4次过极区期间的与$E_{\mathrm{L}\mathrm{O}\mathrm{S}}$相关参数的变化情况，每次选取1颗GPS卫星的连续跟踪弧段为例，其中各曲线黑色部分对应为极区弧段。

图  3  与$E_{\mathrm{L}\mathrm{O}\mathrm{S}}$相关参数的变化曲线

Figure 3.  Variation Curves of Parameters Related to $E_{\mathrm{L}\mathrm{O}\mathrm{S}}$

由图 3可知，在非极区弧段内，$\mathrm{d}\rho$和$c·\mathrm{d}\mathrm{\delta }{t}^s$的曲线变化平缓且光滑，但是由变化率$\mathrm{d}{E}_{\mathrm{L}\mathrm{O}\mathrm{S}}$曲线可知，$E_{\mathrm{L}\mathrm{O}\mathrm{S}}$仍受到小量级的随机误差的影响。在极区弧段内，因受到外推期间广播星历误差的影响，$E_{\mathrm{L}\mathrm{O}\mathrm{S}}$的波动明显大于非极区，其变化特性与广播星历相近。进入极区时，$E_{\mathrm{L}\mathrm{O}\mathrm{S}}$变化不明显，但在出极区时刻，$E_{\mathrm{L}\mathrm{O}\mathrm{S}}$出现显著的台阶变化，这是由于出极区后重新接收到了新的改正数，$E_{\mathrm{L}\mathrm{O}\mathrm{S}}$迅速减小。

针对$E_{\mathrm{L}\mathrm{O}\mathrm{S}}$在极区和非极区的不同变化特性，本文采用分段随机游走过程对伪模糊度A进行建模：

当低轨卫星进入极区弧段时，A不用重新初始化，由于直接使用最后可用时刻的改正数来外推当前时刻的改正数，相当于在广播星历上进行了常值改正，其变化特性与广播星历相近，所以过程噪声要设置为广播星历对应的数值（广播星历$E_{\mathrm{L}\mathrm{O}\mathrm{S}}$变化特性与随机模型设置详见文献[4-5，12]）；在极区弧段进入非极区弧段时，因为$E_{\mathrm{L}\mathrm{O}\mathrm{S}}$出现显著的台阶变化，所以A要重新初始化，其过程噪声设置与实时改正数的精度有关。本文将A参数的过程噪声方差设置为（1 mm/s·dt）²，其中，dt为滤波中前后历元的时间间隔。

为了验证本文提出的实时精密定轨方法的有效性，根据使用GPS星历的不同，本文设计了5种实时精密定轨实验方案，方案1：仅采用GPS广播星历；方案2：采用GPS广播星历和CLK93实时产品，不考虑极区改正数接收中断情况，即全弧段使用CLK93实时产品，伪模糊度参数采用随机常量模型，过程噪声方差设为0，认为改正数没有误差；方案3：采用GPS广播星历和CLK93实时产品，考虑极区改正数接收中断情况，伪模糊度的随机模型采用本文提出的分段随机游走过程；方案4：使用的数据与方案2相同，只是伪模糊度的随机模型选取不同，这里采用随机游走过程；方案5：采用IGS精密轨道和精密钟差（30 s）产品。5种方案定轨策略及设置均相同，如表 1所示。

根据上述5种实验方案设置，利用自主研制的实时精密定轨软件SATODS对2019年DOY第98~104天共7 d的GRACE C卫星GPS双频实测数据模拟连续实时精密定轨数据处理，计算的轨道结果与JPL精密参考轨道求差，并统计位置与速度在R、A、C方向和三维（3D）的精度。

图 4给出了5种实验方案每天的3D位置和速度精度统计。表 2给出了5种实验方案连续7天实时定轨的位置与速度精度。除了方案1仅采用广播星历外，方案2~4都使用CLK93实时产品后，实时定轨精度从30.9 cm提高到优于10 cm，实现了厘米级实时定轨。其中，方案4全弧段使用CLK93实时产品且考虑伪模糊度随机模型与方案5采用IGS精密产品的实时定轨精度相当。

图  4  5种方案的每天位置和速度精度统计

Figure 4.  Daily Position and Velocity Accuracy Statistics of Five Schemes

图 5给出了DOY第101天方案2、3、5分别与方案4的3D位置误差对比。由图 5可以看出，全弧段都使用CLK93实时产品时，方案2的大部分弧段轨道误差要大于方案4，说明不同GPS卫星不同时段的实时改正数精度并不均匀，对实时定轨精度的影响较大。方案4采用随机游走过程的伪模糊度模型，大幅减弱了实时改正数精度不均匀的影响，定轨精度接近于方案5。与方案4相比，方案3考虑了极区实时改正数接收中断情况，采用分段随机游走过程来描述伪模糊度参数，可以消除大部分极区弧段的改正数外推误差的影响，定轨精度可以达到7.04 cm，但部分弧段较差，甚至达到了20 cm。

图  5  4种方案的三维位置误差比较

Figure 5.  3D Position Error Comparison of Schemes 2, 3, 4 and 5

图 6给出了方案3与方案4部分精度较差弧段（对应图 5中虚线框内所示范围）在RAC方向和3D位置误差对比。由图 6可以看出，方案3在R、C方向受极区改正数外推的影响较小，而A方向受影响较大，位置误差曲线的波峰出现在LEO卫星出极区弧段，可能与部分极区弧段改正数的外推值不能准确消除广播星历的误差有关，极区弧段伪模糊度的随机模型切换到广播星历的随机

图  6  方案3和方案4定轨位置误差对比

Figure 6.  Comparison of Orbit Errors for Schemes 3 and 4

模型参数设置，如果该弧段广播星历的误差不显著，定轨误差会得到较好保持，如果该弧段广播星历的误差较大，定轨误差将会出现逐渐增加的现象，定轨精度将会退化到方案1的分米级。但是在出极区后，由于重新初始化了伪模糊度参数，随着新改正数的使用，伪模糊度参数将逐渐重新收敛，非极区弧段的定轨误差又达到了厘米级。

本文利用IGS播发的GPS卫星轨道与钟差改正数实时数据流，提出用分段随机游走过程来构建伪模糊度参数随机模型，减小星历综合误差对实时精密定轨精度的影响，从而实现厘米级星载GPS

实时定轨。基于CLK93实时产品，设计了多种实时定轨实验方案，对GRACE C卫星GPS双频实测数据模拟在轨实时精密定轨处理，得出以下结论：（1）不考虑实时改正数中断的情况下，基于伪模糊度随机模型的实时精密定轨算法，定轨精度优于6 cm（3D RMSE），接近于IGS精密轨道与钟差产品的实时定轨精度；（2）针对极区弧段实时改正数中断情况下，基于本文提出的伪模糊度分段随机游走模型的实时精密定轨算法，定轨精度约为7 cm（3D RMSE），大幅减小实时改正数外推误差对定轨精度的影响。

综上所述，本文提出基于IGS-RTS的厘米级星载GPS实时定轨的数学模型及其处理策略具有可行性和有效性，较好地解决了极区实时改正数中断时外推误差对实时定轨精度的影响，对于部分极区弧段实时定轨精度较差的影响因素，还需要进一步的研究。