一种不受GLONASS码频间偏差影响的实时精密卫星钟差估计方法

摘要

本发明公开了一种不受GLONASS码频间偏差影响的实时精密卫星钟差估计方法，利用GNSS参考站网络进行GPS、GLONASS实时精密卫星钟差估计，通过在函数模型中设置多个独立“时频偏差”参数，吸收参与解算测站接收机的GLONASS码频间偏差，并通过选择参考钟及附加测站ISFB约束条件，实现时频偏差ISFB参数与接收机钟差、GLONASS卫星钟差钟差的有效分离，从而避免不同测站码频间偏差被GLONASS卫星钟差估值吸收所带来的不利影响；与现有技术相比，在不影响GPS实时卫星钟差估值精度的情况下，能明显提高GLONASS实时卫星钟差估值的精度。

说明书

技术领域

本发明涉及卫星大地测量及导航定位技术领域，尤其涉及一种不受GLONASS码频间偏差 影响的实时精密卫星钟差估计方法。

背景技术

高精度的实时卫星钟差改正是支持实时精密单点定位(Real-Time Precise Point  Positioning,RT-PPP)技术的关键。由于星载原子钟易受频率稳定性、环境温度变化等因素 影响，IGS及其分析中心提供的超快速卫星钟差产品在数小时内的预报精度可能下降至数个 ns水平，无法满足RT-PPP的要求。一般利用全球或局域跟踪站网络提供的实时非差载波相 位和伪距观测值，在固定测站位置和IGS超快星历提供预报轨道的前提下进行卫星钟差的实 时估计。基于GPS/GLONASS组合的RT-PPP需同时提供高精度的实时GPS和GLONASS卫星钟差 产品。2013年4月1日，IGS正式推出可提供实时卫星轨道和卫星钟差产品的RTS(Real-time  Service)服务；该项服务目前正处于测试阶段，且仅能提供实时GPS卫星钟差改正产品。迄 今为止，IGS尚未提供官方的GLONASS精密卫星钟差产品，仅有ESOC、IAC、GFZ及NRCan四 个分析中心提供事后的精密卫星钟差产品。随着GLONASS系统的逐步恢复，全球GLONASS跟 踪站网分布较以往有显著增加，已具备提供实时GLONASS精密卫星钟差产品的基础。

与GPS不同，当前的GLONASS系统使用频分多址(frequency division multiple access, FDMA)区分来自不同卫星的信号，由此在接收机内部会产生因通道频率而异的硬件延迟偏差 项(inter-channel bias,ICB)，又称为频间偏差(inter-frequnency bias,IFB)。已有研 究表明，GLONASS码频间偏差在不同频率通道间的差异可达数米，且与接收机类型、固件版 本、天线类型等因素有关。在联合GPS/GLONASS进行卫星钟差估计时，需估计两个接收机钟 差项，即GPS接收机钟差和GLONASS接收机钟差；若GPS与GLONASS卫星钟差采用的时间参 考基准不同，则还应考虑二者的系统时间偏差。现有技术一般通过引入“系统时差”参数， 将GLONASS接收机钟差表达成GPS接收机钟差与系统时差之和的形式。对于GPS接收机码硬 件延迟，由于所有GPS卫星信号采用相同固定频率，该延迟将由GPS接收机钟差吸收；对于 GLONASS接收机码硬件延迟，码频间偏差因频率不同而存在差异，因而引入的系统时差将仅 吸收接收机GLONASS码硬件延迟的公共部分。因此，“系统时差”参数将由三部分组成，即： GPS与GLONASS的系统时间偏差、GLONASS接收机码平均延迟及需扣除的GPS接收机码延迟。 显然，现有技术在进行卫星钟差估计函数建模中并未考虑GLONASS码频间偏差的影响； GLONASS码频间偏差一部分会被估计的GLONASS卫星钟差吸收，残余部分将体现在码伪距观 测值残差中。在进行卫星钟差估计时，由于相位模糊度参数的存在，高精度的相位观测值仅 决定卫星钟差在历元间的精确变化，而伪距观测值对钟差解的主要贡献在于提供所估计钟差 的时间参考基准。GPS/GLONASS精密卫星钟差的联合实时估计，一般利用来自全球或区域分 布的GNSS参考站网络实时观测数据流；对于每颗GLONASS卫星的钟差估值，在不同历元间所 参与求解的测站可能并不相同。现有技术并未考虑GLONASS码伪距观测值中的频间偏差项， 这将导致不同参考站码伪距确定的卫星钟差参考基准很难收敛至稳定值。在扣除卫星钟常数 偏差后，各参考站码频间偏差的综合影响将会使GLONASS卫星钟差估值出现较大范围的波动， 从而导致GLONASS实时卫星钟差估计精度的降低。

发明内容

本发明所要解决的技术问题是针对背景技术中所涉及的缺陷，提供一种不受GLONASS码 频间偏差影响的实时精密卫星钟差估计方法。

本发明为解决上述技术问题采用以下技术方案:

一种不受GLONASS码频间偏差影响的实时精密卫星钟差估计方法，包括以下步骤：

步骤1)，获取来自GNSS参考站网络的实时GPS、GLONASS观测数据；

步骤2)，对步骤1中获取的GPS、GLONASS观测数据进行实时周跳探测，并对发生周跳 的历元进行标记；

步骤3)，对步骤1)中获取的每个GNSS参考站的GPS、GLONASS观测数据，分别组成消 电离层组合观测值；

步骤4)，对步骤1)中获取的每个GNSS参考站的GPS、GLONASS观测数据，分别建立GPS、 GLONASS卫星观测方程，其中将GLONASS接收机码频间偏差参数与系统时差参数进行合并， 为每颗观测的GLONASS卫星均设置一个独立时频偏差参数；

步骤5)，从IGS最新公布的SINEX文件中提取各GNSS参考站位置坐标；从IGS超快星 历提供的预报轨道提取当前观测历元的卫星位置坐标；对观测方程中的各项误差源进行建模 改正；选定一个配备稳定频标的GNSS参考站，将其接收机钟作为参考钟；

步骤6)，将所有参与解算GNSS参考站的GPS/GLONASS卫星观测方程线性化，并表示成 矩阵形式；

步骤7)，引入附加约束条件，即假定每个测站所有待估“时频偏差”之和为0；

步骤8)，采用扩展卡尔曼滤波方法进行参数估计；其中，对标记为发生周跳的历元，将 模糊度参数重置后再进行扩展卡尔曼滤波估计；

步骤9)，输出实时卫星钟差结果。

作为本发明一种不受GLONASS码频间偏差影响的实时精密卫星钟差估计方法进一步的优 化方案，步骤3)中对于测站r、GPS卫星i组成消电离层组合观测值的计算公式为：

$\left(\begin{array}{c}{P}_{\mathrm{IF},i,r}^G=\frac{{f_1}^2{P_1}^G-{f_2}^2{P_2}^G}{{f_1}^2-{f_2}^2}\\ {\Phi }_{\mathrm{IF},i,r}^G=\frac{{f_1}^2{{\Phi }_1}^G-{f_2}^2{{\Phi }_2}^G}{{f_1}^2-{f_2}^2}\end{array}\right)$

式中，G表示GPS系统；分别为码伪距和载波相位的消电离层组合观测值；f₁、 f₂分别表示L₁和L₂载波的频率；P₁^G、P₂^G与Φ₁^G、Φ₂^G分别为测站r相应频率的码伪距和载波 相位原始观测值；

对于测站r、GLONASS卫星j组成消电离层组合观测值的计算公式为：

$\left(\begin{array}{c}{P}_{\mathrm{IF},j,r}^R=\frac{f_{K,1}^2{P_1}^R-f_{K,2}^2{P_2}^R}{f_{K,1}^2-f_{K,2}^2}\\ {\Phi }_{\mathrm{IF},j,r}^R=\frac{f_{K,1}^2{{\Phi }_1}^R-f_{K,2}^2{{\Phi }_2}^R}{f_{K,1}^2-f_{K,2}^2}\end{array}\right)$

式中，R表示GLONASS系统，分别为码伪距和载波相位的消电离层组合观测 值；K表示卫星j所对应的GLONASS频率通道号；f_K,1、f_K,2分别表示卫星j所对应的GLONASS 频率通道K的L₁和L₂载波的频率；P₁^R、P₂^R与Φ₁^R、Φ₂^R分别为测站r相应频率的码伪距和载 波相位原始观测值。

作为本发明一种不受GLONASS码频间偏差影响的实时精密卫星钟差估计方法进一步的优 化方案，步骤4)中对于测站r、GPS卫星i，其观测方程表示为：

$\left(\begin{array}{c}{P}_{\mathrm{IF},i,r}^G={\rho }_r^G+c({\mathrm{dt}}_r^G-{\mathrm{dt}}_i^{G/s})+{\delta }_{{\mathrm{trop}}_{i,r}}+{\delta }_{\mathrm{tide}}+{\delta }_{\mathrm{rel}}+{ϵ}_P^G\\ {\Phi }_{\mathrm{IF},i,r}^G={\rho }_r^G+c({\mathrm{dt}}_r^G-{\mathrm{dt}}_i^{G/s})+{\lambda }^G{N}_{i,r}^G+{\delta }_{{\mathrm{trop}}_{i,r}}+{\delta }_{\mathrm{tide}}+{\delta }_{\mathrm{rel}}+{\delta }_{\mathrm{phw}}+{ϵ}_{\Phi }^G\end{array}\right)$

式中，为GPS卫地距；分别为待估的GPS接收机钟差和卫星钟差；为 GPS卫星i对应的消电离层相位模糊度；为对流层斜延迟；δ_tide为潮汐效应改正；δ_rel为 相对论效应改正；δ_phw为天线相位缠绕改正；为观测噪声；c为真空中光速，λ^G为GPS 消电离层相位观测值所对应的波长；

对于测站r、GLONASS卫星j，其观测方程表示为：

$\left(\begin{array}{c}{P}_{\mathrm{IF},j,r}^R={\rho }_r^R+c({\mathrm{dt}}_r^G-{\mathrm{dt}}_j^{R/s})+{\beta }_{j,r}+{\delta }_{\mathrm{tro}{p}_{j,r}}+{\delta }_{\mathrm{tide}}+{\delta }_{\mathrm{rel}}++{ϵ}_P^R\\ {\Phi }_{\mathrm{IF},j,r}^R={\rho }_r^R+c({\mathrm{dt}}_r^G-{\mathrm{dt}}_j^{R/s})+{\beta }_{j,r}+{\lambda }_j^R{N}_{j,r}^R+{\delta }_{{\mathrm{trop}}_{j,r}}+{\delta }_{\mathrm{tide}}+{\delta }_{\mathrm{rel}}+{ϵ}_{\Phi }^R\end{array}\right)$

式中，为GLONASS卫地距；为待估的GLONASS卫星钟差；为GLONASS卫星j 对应的消电离层相位模糊度；为卫星j消电离层相位观测值所对应的波长；β_j,r为待估的 GLONASS卫星j的时频偏差参数，具体表示为：

${\beta }_{j,r}={\mathrm{dt}}_{\mathrm{SYS}}+M_r^R-D_r^G+{\bar{B}}_{j,r}^R,(j=\mathrm{1,2},...,n)$

式中，dt_SYS为系统时间偏差；为GLONASS接收机码平均延迟偏差；为GPS接收机 码延迟偏差；为GLONASS接收机码频间偏差，具体表示为：

${\bar{B}}_{j,r}^R=\frac{f_{K,1}^2{\hat{B}}_{j,r}^{P1}-f_{K,2}^2{\hat{B}}_{j,r}^{P2}}{f_{K,1}^2-f_{K,2}^2}$

式中，为原始观测值对应的GLONASS接收机码频间偏差。

作为本发明一种不受GLONASS码频间偏差影响的实时精密卫星钟差估计方法进一步的优 化方案，步骤6)中观测方程线性化后写成矩阵形式表示为：

$\left(\begin{array}{c}\underset{p\times 1}{{\tilde{P}}_{\mathrm{IF}}^G\left(t\right)}\\ \underset{p\times 1}{{\tilde{\Phi }}_{\mathrm{IF}}^G\left(t\right)}\\ \underset{q\times 1}{{\tilde{P}}_{\mathrm{IF}}^R\left(t\right)}\\ \underset{q\times 1}{{\tilde{\Phi }}_{\mathrm{IF}}^R\left(t\right)}\end{array}\right)=\left(\begin{array}{ccccccc}\underset{2p\times (k-1)}{A_r^G}& \underset{2p\times k}{A_{\mathrm{zwd}}^G}& \underset{2p\times m}{A_s^G}& \underset{2p\times n}{O}& \underset{2p\times p}{A_N^G}& \underset{2p\times q}{O}& \underset{2p\times q}{O}\\ \underset{2q\times (k-1)}{A_r^R}& \underset{2q\times k}{A_{\mathrm{zwd}}^R}& \underset{2q\times m}{O}& \underset{2q\times n}{A_s^R}& \underset{2q\times n}{A_s^R}& \underset{2q\times q}{A_N^R}& \underset{2q\times q}{A_{\mathrm{ISFB}}}\end{array}\right)\left(\begin{array}{c}\underset{(k-1)\times 1}{{\mathrm{\delta t}}_r^G}\\ \underset{k\times 1}{{\delta }_{\mathrm{zwd}}}\\ \underset{m\times 1}{{\mathrm{\delta t}}^{G/s}}\\ \underset{n\times 1}{{\mathrm{\delta t}}^{R/s}}\\ \underset{p\times 1}{N_{\mathrm{IF}}^G}\\ \underset{q\times 1}{N_{\mathrm{IF}}^R}\\ \underset{q\times 1}{{\mathrm{\delta t}}_{\mathrm{ISFB}}}\end{array}\right)+\left(\begin{array}{c}{ϵ}^G\\ {ϵ}^R\end{array}\right)$

式中，m、n分别为在轨GPS、GLONASS卫星数量，k为历元t参与解算测站的数量，p、 q为所有测站GPS、GLONASS可见卫星数，为观测向量；δ_zwd、δt^G/s、δt^R/s、及δt_ISFB为待估参数向量，分别对应GPS接收机钟差、对流层 天顶湿延迟、GPS卫星钟差、GLONASS卫星钟差、GPS相位模糊度、GLONASS相位模糊度及时 频偏差；A表示设计阵，下标r、zwd、s、N、ISFB分别表示接收机钟差、对流层天顶湿 延迟、卫星钟差、相位模糊度及时频偏差，上标G、R分别表示GPS与GLONASS系统；其中 I_q为q维单位阵，表示克罗内克积；O表示零矩阵；ε^G、 ε^R为测量噪声向量，服从零均值的整体分布，其对应的协方差矩阵Ω^G、Ω^R均为对角阵， 按下式计算矩阵的对角线元素：

${\sigma }^2=\frac{{\sigma }_0^2}{\mathrm{si}{n}^{2}E}$

式中，σ为观测值先验中误差；E为测站处的卫星高度角；σ₀为天顶方向观测值的先验 中误差。

作为本发明一种不受GLONASS码频间偏差影响的实时精密卫星钟差估计方法进一步的优 化方案，步骤7)中假定每个测站所有待估“时频偏差”之和为0的公式为：

$\underset{j=1}{\overset{n}{\Sigma }}{\mathrm{dt}}_{\mathrm{ISFB},j,r}=0$

将附加约束条件视为虚拟观测方程，其矩阵形式表示为：

$O_{\mathrm{ISFB}}=A_c·{\mathrm{\delta t}}_{\mathrm{ISF}{B}_r}$

式中，为测站r对应的时频偏差ISFB参数向量；A_c为对应的设计阵，其中 e_k为各元素均为1的k维列向量，I_u为单位阵，下标u表示该测站可见GLONASS 卫星数；O_ISFB为零矩阵。

作为本发明一种不受GLONASS码频间偏差影响的实时精密卫星钟差估计方法进一步的优 化方案，步骤8)中所述扩展卡尔曼滤波的状态方程表示如下：

X(t)＝Ψ_t,t-1X(t-1)+w_t-1,w_t-1～(0,Ω_w)

式中，X(t)、X(t-1)为相邻历元的系统状态向量；Ψ_t,t-1为状态转移矩阵；w_t-1为状态 噪声向量，且服从均值为0、协方差矩阵为Ω_w的正态分布。

状态转移矩阵Ψ_t,t-1表示为：

${\Psi }_{t,t-1}=\left(\begin{array}{ccccc}\underset{(k-1)\times (k-1)}{O}& & & & \\ & \underset{k\times k}{I}& & & \\ & & \underset{(m+n)\times (m+n)}{I}& & \\ & & & \underset{(p+q)\times (p+q)}{I}& \\ & & & & \underset{q\times q}{I}\end{array}\right)$

式中，O表示零矩阵；I表示单位阵；

协方差矩阵Ω_w表示为：

${\Omega }_w=\left(\begin{array}{ccccc}{q}_{{\mathrm{\delta t}}_r^G}\Delta & & & & \\ & q_{{\delta }_{\mathrm{zwd}}}\mathrm{\Delta t}& & & \\ & & \underset{(m+n)\times (m+n)}{I}& & \\ & & & \underset{(p+q)\times (p+q)}{O}& \\ & & & & \underset{q\times q}{q_{\beta }}\mathrm{\Delta t}\end{array}\right)$

式中，Δt为相邻历元的时间间隔，q_β分别表示待求GPS接收机钟差、对流 层天顶湿延迟及时频偏差参数。

作为本发明一种不受GLONASS码频间偏差影响的实时精密卫星钟差估计方法进一步的优 化方案，所述步骤4)的GLONASS观测方程中接收机钟差参数与GPS观测方程中的相同，且 引入的时频偏差参数由GPS与GLONASS的系统时间偏差、GLONASS接收机码平均延迟、需扣 除的GPS接收机码延迟及GLONASS接收机码频间偏差四部分组成。

作为本发明一种不受GLONASS码频间偏差影响的实时精密卫星钟差估计方法进一步的优 化方案，将时频偏差参数变化视为随机游走过程处理，其先验中误差取值

本发明采用以上技术方案与现有技术相比，具有以下技术效果：

在联合GPS/GLONASS进行实时精密卫星钟差估计过程中，传统方法在函数建模中未考虑 GLONASS码伪距观测值中的频间偏差项，所估计的GLONASS卫星钟差将吸收来自不同GNSS参 考站的部分GLONASS接收机码频间偏差。在不同历元间所参与求解的测站可能并不相同，各 参考站码频间偏差的综合影响将会使GLONASS卫星钟差估值出现较大范围的波动，从而造成 GLONASS实时卫星钟差估计精度的降低。本发明通过在函数模型中为每个参与解算测站设置 多个独立“时频偏差”参数，吸收该测站接收机的GLONASS码频间偏差，并通过选择参考钟 及附加测站ISFB约束条件，实现时频偏差ISFB参数与接收机钟差、GLONASS卫星钟差钟差 的有效分离，避免了不同测站码频间偏差被GLONASS卫星钟差估值吸收所带来的不利影响； 与现有技术相比，在不影响GPS实时卫星钟差估值精度的情况下，能明显提高GLONASS实时 卫星钟差估值的精度。

附图说明

图1是本发明方法的流程图；

图2是本发明方法得到的GPS卫星钟差实时估计结果偏差图；

图3是传统方法得到的GLONASS卫星钟差实时估计结果偏差图；

图4是本发明方法得到的GLONASS卫星钟差实时估计结果偏差图；

图5是两种方法得到的GPS、GLONASS卫星钟差实时估计结果精度对比图；

图6是利用本发明方法得到的GPS、GLONASS实时卫星钟差进行实时精密单点定位的结果 偏差RMS值统计图。

具体实施方式

下面结合附图对本发明的技术方案做进一步的详细说明，但是本发明的保护范围不局限 于所述实施例。

如图1所示，本发明提供的一种不受GLONASS码频间偏差影响的实时精密卫星钟差估 计方法，具体步骤如下：

步骤1：获取来自GNSS参考站网络的实时GPS、GLONASS观测数据；为检验本发明所 提供方法的有效性，选取全球均匀分布的30个GPS/GLONASS跟踪站用于卫星钟差的实时 估计，测站分布如图1所示；从IGS下载所选测站30s间隔观测数据RINEX文件用于模拟 实时数据流，观测时间为2014年年积日(Day of Year,DOY)41；将GPS/GLONASS卫星钟 差实时估计结果与欧洲航天局(European Space Agency,ESA)提供的事后精密卫星钟差产品 进行比较，用于评估本发明方法用于GPS/GLONASS卫星钟差实时估计所能达到的精度；

步骤2：对步骤1中获取的GPS、GLONASS实时观测数据流进行周跳探测，并对发生周 跳的历元进行标记；

步骤3：对步骤1获取的每个GNSS参考站的GPS、GLONASS观测数据，分别组成消电 离层组合观测值，通过以下公式计算得到：

对于测站r、GPS卫星i：

$\left(\begin{array}{c}{P}_{\mathrm{IF},i,r}^G=\frac{{f_1}^2{P_1}^G-{f_2}^2{P_2}^G}{{f_1}^2-{f_2}^2}\\ {\Phi }_{\mathrm{IF},i,r}^G=\frac{{f_1}^2{{\Phi }_1}^G-{f_2}^2{{\Phi }_2}^G}{{f_1}^2-{f_2}^2}\end{array}\right)---\left(1\right),$

式中，分别为码伪距和载波相位的消电离层组合观测值；f₁、f₂分别表示L₁和L₂载波的频率；P₁^G、P₂^G与Φ₁^G、Φ₂^G分别为测站r相应频率的码伪距和载波相位原始观测 值；

对于测站r、GLONASS卫星j：

$\left(\begin{array}{c}{P}_{\mathrm{IF},j,r}^R=\frac{f_{K,1}^2{P_1}^R-f_{K,2}^2{P_2}^R}{f_{K,1}^2-f_{K,2}^2}\\ {\Phi }_{\mathrm{IF},j,r}^R=\frac{f_{K,1}^2{{\Phi }_1}^R-f_{K,2}^2{{\Phi }_2}^R}{f_{K,1}^2-f_{K,2}^2}\end{array}\right)---\left(2\right),$

式中，与GPS计算公式符号相同的参数其含义是相同的，此处用上标R表示GLONASS 系统，下标K表示卫星j所对应的GLONASS频率通道号；

步骤4：对步骤1获取的每个GNSS参考站的GPS、GLONASS观测数据，分别建立GPS、 GLONASS卫星观测方程，其中将GLONASS接收机码“频间偏差”与“系统时差”参数进 行合并，为每颗观测的GLONASS卫星均设置一个独立“时频偏差”(Inter-System and Inter -Frequency Bias,ISFB)参数，具体表示如下：

对于测站r、GPS卫星i，其观测方程表示为：

$\left(\begin{array}{c}{P}_{\mathrm{IF},i,r}^G={\rho }_r^G+c({\mathrm{dt}}_r^G-{\mathrm{dt}}_i^{G/s})+{\delta }_{{\mathrm{trop}}_{i,r}}+{\delta }_{\mathrm{tide}}+{\delta }_{\mathrm{rel}}+{ϵ}_P^G\\ {\Phi }_{\mathrm{IF},i,r}^G={\rho }_r^G+c({\mathrm{dt}}_r^G-{\mathrm{dt}}_i^{G/s})+{\lambda }^G{N}_{i,r}^G+{\delta }_{{\mathrm{trop}}_{i,r}}+{\delta }_{\mathrm{tide}}+{\delta }_{\mathrm{rel}}+{\delta }_{\mathrm{phw}}+{ϵ}_{\Phi }^G\end{array}\right)---\left(3\right),$

式中，为GPS卫地距；分别为待估的GPS接收机钟差和卫星钟差；为 GPS卫星i对应的消电离层相位模糊度；为对流层斜延迟；δ_tide为潮汐效应改正；δ_rel为 相对论效应改正；δ_phw为天线相位缠绕改正；为观测噪声；

对于测站r、GLONASS卫星j，其观测方程表示为：

$\left(\begin{array}{c}{P}_{\mathrm{IF},j,r}^R={\rho }_r^R+c({\mathrm{dt}}_r^G-{\mathrm{dt}}_j^{R/s})+{\beta }_{j,r}+{\delta }_{\mathrm{tro}{p}_{j,r}}+{\delta }_{\mathrm{tide}}+{\delta }_{\mathrm{rel}}++{ϵ}_P^R\\ {\Phi }_{\mathrm{IF},j,r}^R={\rho }_r^R+c({\mathrm{dt}}_r^G-{\mathrm{dt}}_j^{R/s})+{\beta }_{j,r}+{\lambda }_j^R{N}_{j,r}^R+{\delta }_{{\mathrm{trop}}_{j,r}}+{\delta }_{\mathrm{tide}}+{\delta }_{\mathrm{rel}}+{ϵ}_{\Phi }^R\end{array}\right)---\left(4\right),$

式中，与GPS观测方程符号相同的参数其含义是相同的，此处用上标R表示GLONASS 系统，β_j,r为待估的GLONASS卫星j的时频偏差ISFB参数，具体表示为：

${\beta }_{j,r}={\mathrm{dt}}_{\mathrm{SYS}}+M_r^R-D_r^G+{\bar{B}}_{j,r}^R,(j=\mathrm{1,2},...,n)---\left(5\right),$

式中，dt_SYS为系统时间偏差；为GLONASS接收机码平均延迟偏差；为GPS接收 机码延迟偏差；为GLONASS接收机码频间偏差，具体表示为：

${\bar{B}}_{j,r}^R=\frac{f_{K,1}^2{\hat{B}}_{j,r}^{P1}-f_{K,2}^2{\hat{B}}_{j,r}^{P2}}{f_{K,1}^2-f_{K,2}^2}---\left(6\right),$

式中，为原始观测值对应的GLONASS接收机码频间偏差，f_K,1、f_K,2为相应 的信号频率，对应的GLONASS频率通道号为K；

注意本步骤中GLONASS观测方程中接收机钟差参数与GPS观测方程中的相同，GPS与 GLONASS系统时间差及GLONASS码频间偏差均由时频偏差(ISFB)参数吸收，并作为待 估参数进行估计。

步骤5：从IGS最新公布的SINEX文件中提取各GNSS参考站位置坐标；从ESA超快星 历提供的预报轨道提取当前观测历元的卫星位置坐标；选择测站MATE的接收机钟作为参考 钟；对观测方程中的各项误差源进行建模改正，具体策略如表1所示；

表1：实时卫星钟差估计观测模型及误差源改正策略

步骤6：将所有参与解算GNSS参考站的GPS/GLONASS卫星观测方程线性化，并表示 成矩阵形式；假设在轨GPS、GLONASS卫星分别为m、n颗，在历元t共有k个测站参与解 算，所有测站GPS、GLONASS可见卫星数分别为p、q颗，观测方程线性化后写成矩阵形 式表示为：

$\left(\begin{array}{c}\underset{p\times 1}{{\tilde{P}}_{\mathrm{IF}}^G\left(t\right)}\\ \underset{p\times 1}{{\tilde{\Phi }}_{\mathrm{IF}}^G\left(t\right)}\\ \underset{q\times 1}{{\tilde{P}}_{\mathrm{IF}}^R\left(t\right)}\\ \underset{q\times 1}{{\tilde{\Phi }}_{\mathrm{IF}}^R\left(t\right)}\end{array}\right)=\left(\begin{array}{ccccccc}\underset{2p\times (k-1)}{A_r^G}& \underset{2p\times k}{A_{\mathrm{zwd}}^G}& \underset{2p\times m}{A_s^G}& \underset{2p\times n}{O}& \underset{2p\times p}{A_N^G}& \underset{2p\times q}{O}& \underset{2p\times q}{O}\\ \underset{2q\times (k-1)}{A_r^R}& \underset{2q\times k}{A_{\mathrm{zwd}}^R}& \underset{2q\times m}{O}& \underset{2q\times n}{A_s^R}& \underset{2q\times n}{A_s^R}& \underset{2q\times q}{A_N^R}& \underset{2q\times q}{A_{\mathrm{ISFB}}}\end{array}\right)\left(\begin{array}{c}\underset{(k-1)\times 1}{{\mathrm{\delta t}}_r^G}\\ \underset{k\times 1}{{\delta }_{\mathrm{zwd}}}\\ \underset{m\times 1}{{\mathrm{\delta t}}^{G/s}}\\ \underset{n\times 1}{{\mathrm{\delta t}}^{R/s}}\\ \underset{p\times 1}{N_{\mathrm{IF}}^G}\\ \underset{q\times 1}{N_{\mathrm{IF}}^R}\\ \underset{q\times 1}{{\mathrm{\delta t}}_{\mathrm{ISFB}}}\end{array}\right)+\left(\begin{array}{c}{ϵ}^G\\ {ϵ}^R\end{array}\right)---\left(7\right),$

式中，为观测向量；δ_zwd、δt^G/s、δt^R/s、及δt_ISFB为待估参数向量，分别对应GPS接收机钟差、对流层天顶湿延迟、GPS卫星钟差、 GLONASS卫星钟差、GPS相位模糊度、GLONASS相位模糊度及时频偏差；A表示设计阵， 其中I_q为q维单位阵，表示克罗内克积；O表示零矩 阵；ε^G、ε^R为测量噪声向量，服从零均值的整体分布，其对应的协方差矩阵Ω^G、Ω^R均为 对角阵，按下式计算矩阵的对角线元素：

${\sigma }^2=\frac{{\sigma }_0^2}{\mathrm{si}{n}^{2}E}---\left(8\right),$

式中，σ为观测值先验中误差；E为测站处的卫星高度角；σ₀为天顶方向观测值的先验 中误差；

步骤7：引入附加约束条件，即假定每个测站所有待估“时频偏差”之和为0，即：

$\underset{j=1}{\overset{n}{\Sigma }}{\mathrm{dt}}_{\mathrm{ISFB},j,r}=0---\left(9\right),$

将附加约束条件视为虚拟观测方程，其矩阵形式表示为：

$O_{\mathrm{ISFB}}=A_c·{\mathrm{\delta t}}_{\mathrm{ISF}{B}_r}---\left(10\right),$

式中，为测站r对应的时频偏差ISFB参数向量；A_c为对应的设计阵，其中 为各元素均为1的k维列向量，I_u为单位阵，下标u表示该测站可见GLONASS 卫星数；O_ISFB为零矩阵。

步骤8：采用扩展卡尔曼滤波(Extended Kalman Filter,EKF)方法进行参数估计，其中对 标记为发生周跳历元的模糊度参数进行重置。EKF状态方程表示如下：

X(t)＝Ψ_t,t-1X(t-1)+w_t-1,w_t-1～(0,Ω_w)   (11)，

式中，X(t)、X(t-1)为相邻历元的系统状态向量；Ψ_t,t-1为状态转移矩阵；w_t-1为状态噪 声向量，且服从均值为0、协方差矩阵为Ω_w的正态分布。

状态转移矩阵Ψ_t,t-1表示为：

${\Psi }_{t,t-1}=\left(\begin{array}{ccccc}\underset{(k-1)\times (k-1)}{O}& & & & \\ & \underset{k\times k}{I}& & & \\ & & \underset{(m+n)\times (m+n)}{I}& & \\ & & & \underset{(p+q)\times (p+q)}{I}& \\ & & & & \underset{q\times q}{I}\end{array}\right)---\left(12\right),$

式中，O表示零矩阵；I表示单位阵；

协方差矩阵Ω_w表示为：

${\Omega }_w=\left(\begin{array}{ccccc}{q}_{{\mathrm{\delta t}}_r^G}\Delta & & & & \\ & q_{{\delta }_{\mathrm{zwd}}}\mathrm{\Delta t}& & & \\ & & \underset{(m+n)\times (m+n)}{I}& & \\ & & & \underset{(p+q)\times (p+q)}{O}& \\ & & & & \underset{q\times q}{q_{\beta }}\mathrm{\Delta t}\end{array}\right)---\left(13\right),$

式中，Δt为相邻历元的时间间隔，q表示各状态参数的谱密度矩阵。在随机模型中，将 GPS接收机钟差视为白噪声处理，其先验中误差设置为10²m；天顶对流层湿延迟视为随机游 走过程处理，其先验中误差设置为由于相邻历元间隔Δt通常较小，在此较 短时间内可以认为ISFB参数变化是稳定的，因而同样将其视为随机游走过程处理，其先验中 误差设置为GPS与GLONASS模糊度参数均视为不随时间变化的量。

为分析本发明提供方法的效果，分别利用传统方法和本发明方法对实验数据进行计算。 传统方法：每个测站估计一个“系统时差”参数，函数模型中不考虑GLONASS码频间偏差； 本发明方法：每个测站每颗观测GLONASS卫星均估计一个“时频偏差”参数，函数模型中 考虑各测站的GLONASS码频间偏差。以上两种情况仅待估参数不同，其观测模型与参数估 计策略均设为相同，具体信息如表1、表2所示。

表2：实时卫星钟差估计参数估计策略

注：x₀、σ₀、σ分别为参数初值、先验中误差及中误差；为伪距单点定位求得GPS与 GLONASS接收机钟差；为ESA超快GPS与GLONASS卫星预报钟差；为所选择 参考站的GPS接收机钟差。

为说明方便起见，GPS卫星钟按星座及频标类型任选一颗卫星，GLONASS卫星钟按信 号频率选择一颗卫星，将本发明发法所估计卫星钟差的结果与对比传统方法进行对比；相关 卫星信息如表3所示。

表3：用于钟差估计结果分析的卫星钟信息

由于GLONASS码频间偏差对GPS卫星钟差估计影响甚微，两种方法得到的GPS钟差结 果RMS值差异非常小(小于±0.01ns)，在此仅给出顾及码频间偏差处理方案的GPS实时钟 差估计结果；2014年DOY 41全天的GPS卫星钟差实时估计结果偏差如图2所示。从图2可 以看出，在卫星钟差估计初始阶段，由于受伪距误差的影响，各卫星初始钟偏差的收敛需要 较长时间，GPS卫星钟差结果波动范围较大；在经过5～6h收敛之后，各卫星钟差估值偏差 均能稳定在±0.1m范围内。

图3、图4分别给出了采用传统方法与本发明方法所得到的GLONASS实时卫星钟差结 果。对比图2、图3可以发现，传统方法得到的GLONASS卫星钟差估值偏差要明显大于本 发明方法，经过数小时收敛后仍在±0.2m范围内波动；而采用本发明方法，多数卫星钟差估 值在±0.1m范围内波动。

在GPS/GLONASS实时卫星钟差估计收敛后，计算2014年DOY 41全天的卫星钟差估值 与ESA事后精密钟差产品互差的RMS值，如图5所示。可以进一步看出，两种方法所得到 的GPS卫星钟差精度统计值较为符合，且平均RMS值均小于0.2ns；对于GLONASS卫星 钟差，传统方法估值偏差的RMS值要明显大于本发明方法，少数卫星RMS值接近1ns，而 本发明方法多数卫星RMS值小于0.3ns。

图6为利用本发明方法估计得到的GPS/GLONASS卫星钟差进行实时精密单点定位 (Real-Time Precise Point Positioning,RT-PPP)的3D偏差RMS统计结果。用于测试的数据来 自欧洲区域的10个IGS参考站，观测时间为2014年DOY 41～68，按时段长2h、6h、12h、 24h进行划分，预报卫星轨道由ESA发布的超快星历获得，分别在GLONASS单系统、GPS 单系统及GPS/GLONASS组合系统3种模式下进行静态RT-PPP，并计算其与IGS公布测站 坐标互差RMS值。从图5可以看出，当观测时段长为2h时，GPS/GLONASS组合系统的定 位精度要优于GPS单系统，所有测站偏差RMS值均小于10cm，而GLONASS单系统3D定 位精度约为15cm左右；当观测时段长超过6h后，GPS单系统与GPS/GLONASS组合系统 的定位精度基本相当；对于24h实时单天解，三种模式下进行RT-PPP均能获得优于5cm的 三维定位精度。

本发明的具体实施方式中未涉及的说明属于本领域公知的技术，可参考公知技术加以实 施。