高精度的时间同步是各个GNSS卫星定位导航系统的关键技术，是定位导航精度的重要保证。时间比对算法是实现时间同步技术的核心，是确定不同地面接收机时间差的基础算法，时间比对算法精度直接决定时间同步精度^[1]。

现有时间比对算法有共视时间比对、全视时间比对和载波相位时间比对等算法^[2-5]。全视时间比对需要多个接收机同时观测多颗卫星，对接收机和比对算法都有较高要求。载波相位法虽然精度较高，但观测设备的成本也相对较高，且算法复杂，不易实现。基于伪距观测的共视时间比对，只需被测量的接收机同时观测一颗卫星即可实现时间比对，是易于实现且可以满足精度要求的时间比对方式。

测量的伪距受定轨精度、传输时延等影响，需要结合GPS、Galileo、COMPASS和GLONASS系统特点进行补偿和修正，通常通过引入滤波器来实现。为了减弱或消除随机误差对时间比对结果的影响进而提高时间比对精度，文献[6]将自适应卡尔曼滤波算法应用在共视时间比对算法中, 有效的剔除野值，使视察数据更接近真实钟差，提高了卫星共视系统的比对精度；文献[7]提出了Vondrak平滑和三次样条插值处理，打破了传统的曲线拟合和滑动平均方法中多项式项数需人为选取和数据必须为等间距的限制。文献[8]提出了基于离散小波的多尺度卡尔曼滤波算法，用理论方程代替经验方程，有效的减小了人为误差。

现有滤波和平滑方法多针对当前解算的伪距进行修正，并未较好利用整个观测段的所有伪距信息。因此，本文将借助卡尔曼滤波算法消除计算过程中随机误差，对观测数据进行前处理，在卡尔曼滤波算法后加入RTS(rauch-tung-striebel)事后平滑算法，据此在卡尔曼滤波的基础上，进一步提高时间比对精度。

1 基于伪距观测的共识时间比对算法 1.1 建立伪距观测方程

伪距测量是基于伪距观测的GNSS共视时间比对算法的基础，其过程是位于地球上不同位置的多个地面接收机同时对同一颗本系统的GNSS卫星导航信号中的时标进行观测，并将本地时钟钟面时与卫星时标信号解算出的卫星钟面时作差再乘以信号传播速度，进而计算出地面接收机与被观测卫星的伪距值，再将这些伪距值通过互联网进行传递^[9]。伪距测量过程如图 1所示。

伪距观测方程为

$ \left\{ \begin{array}{l} \rho _i^{\left( s \right)} = {r_i} + \delta {t_{ui}} - \delta {t_i}^{\left( s \right)} + {I_i} + {T_i} + {\varepsilon _{{\rho _i}}}\\ \rho _j^{\left( s \right)} = {r_j} + \delta {t_{uj}} - \delta {t_j}^{\left( s \right)} + {I_j} + {T_j} + {\varepsilon _{{\rho _j}}} \end{array} \right. $

(1)

式中：s为卫星，u为接收机，i和j为不同的接收机。ρ^(s)是接收机的伪距观测值，δt^(s)为卫星钟差，I为电离层延时，T为对流层延时，r为接收机到卫星的真实距离，ε_ρ为伪距观测噪声，在设计算法时可以忽略。

1.2 电离层延时修正值的获取

双频接收机电离层延时估算算法可以适用于GPS、GLONASS、Galileo和COMPASS等系统^[10]，这种算法的优势在于不需要数学模型，双频接收机通过伪距观测和计算遍可获得实时的电离层时延值^[11]。公式如下

$ \left\{ \begin{array}{l} {I_1} = \frac{{f_2^2}}{{f_1^2 - f_2^2}}\left( {\rho _{{L_2}}^{\left( s \right)} - \rho _{{L_1}}^{\left( s \right)}} \right) = \frac{1}{{1 - {\gamma _{12}}}}\left( {\rho _{{L_2}}^{\left( s \right)} - \rho _{{L_1}}^{\left( s \right)}} \right)\\ {I_2} = \frac{{{\gamma _{12}}}}{{1 - {\gamma _{{1_2}}}}}\left( {\rho _{{L_2}}^{\left( s \right)} - \rho _{{L_1}}^{\left( s \right)}} \right) \end{array} \right. $

(2)

式中：I₁和I₂分别为双频信号传播过程中的电离层延迟修正值，ρ_L1^(s)和ρ_L2^(s)分别为接收机通过L₁信道和L₂信道观测的伪距值，f₁和f₂分别不同信道的频率，γ₁₂为f₁和f₂比值的平方。

1.3 对流层延时修正值的获取

在估算GPS、Galileo、GLONASS和COMPASS测量值中的对流层时延时采用统一的Hopfield模型，可分为干分量时延和湿分量时延两种情况^[12]。

对流层延时天顶方向干分量T_d的估算公式为

$ {T_d} = 1.552 \times {10^{ - 5}}\frac{{{P_0}}}{{{T_{{k_0}}}}}{H_d} $

(3)

式中：P₀与T_k0分别代表在地面上高度为零处的大气总压力与热力学温度，H_d代表地面到对流层干分量边界的距离。

天顶向对流层延时湿分量T_w的估算公式为

$ {T_w} = 0.074\;6\frac{{{e_{00}}}}{{T_{k0}^2}}{H_w} $

(4)

式中：e₀₀=11.691 MPa为地面零高度处的水汽分压，H_w为地面到对流层湿分量边界的距离。

在信号传播方向上的对流层延时T，即

$ T = {T_d}{F_d} + {T_w}{F_w} $

(5)

干分量倾斜率F_d的估算模型为

$ {F_d} = \frac{1}{{\sin \sqrt {{\theta ^2} + {{\left( {\frac{{2.5{\rm{ \mathsf{ π} }}}}{{180}}} \right)}^2}} }} $

(6)

湿分量倾斜率F_w的估算模型为

$ {F_w} = \frac{1}{{\sin \sqrt {{\theta ^2} + {{\left( {\frac{{1.5{\rm{ \mathsf{ π} }}}}{{180}}} \right)}^2}} }} $

(7)

式中θ为卫星与地面接收机之间形成的高度角。

1.4 地球自转效应修正

利用卫星星历参数进行卫星轨道计算，其计算结果是卫星在信号发射t时刻时的GPS地心地固坐标系中的坐标值。尽管接收机测量的是在这一信号发射时刻t时的卫星位置至此信号被接收时(即信号接收时刻t+τ，其中τ为信号从卫星到接收机的传播时间)的接收机位置之间的几何距离，但是由于地球自转以及依附在地球上的地心地固坐标系统在τ时段内的相应旋转，必须对卫星在信号发射时刻t时的地心地固坐标系中的这一位置坐标值进行地球自转校正，将它转变成在信号接收时刻t+τ时的地心地固坐标系中的位置坐标值。

1.5 确定卫星到接收机的真实距离

式(1)中r为接收机至卫星的几何距离，卫星位置可依据卫星星历计算获得，接收机位置也是已知的，如下

$ \left\{ \begin{array}{l} {{\rm{r}}_i} = \sqrt {{{\left( {x_i^s - {x_{{\rm{u}}i}}} \right)}^2} + {{\left( {y_i^s - {y_{ui}}} \right)}^2} + {{\left( {z_i^s - {z_{ui}}} \right)}^2}} \\ {{\rm{r}}_j} = \sqrt {{{\left( {x_j^s - {x_{{\rm{u}}j}}} \right)}^2} + {{\left( {y_j^s - {y_{uj}}} \right)}^2} + {{\left( {z_j^s - {z_{uj}}} \right)}^2}} \end{array} \right. $

(8)

式中：d为接收机，x^s、y^s、z^s为卫星在地心地固坐标系中的坐标值，x_u、y_u、z_u为接收机的位置坐标。GPS、COMPASS和Galileo系统通过解算卫星星历电文可以实时获得卫星轨道位置^[11]，GLONASS可以采用卫星轨道推算的方法。。

1.6 单星双频时间差观测方程

由1.1节给出的伪距观测方程推导出两地面接收机时间差的表达式为

$ \begin{array}{*{20}{c}} {{\delta _u} = \frac{{{\rho _{{\rm{i}}\_{L_2}}} - {\rho _{j\_{L_2}}} + {\gamma _{12}}{\rho _{j\_{L_1}}} - {\gamma _{12}}{\rho _{i\_{L_1}}}}}{{1 - {\gamma _{12}}}} + }\\ {{r_j} - {r_i} + {T_j} - {T_i} + {\varepsilon _{ij}}} \end{array} $

(9)

式中：δ_u为两接收机时间差，ρ为观测伪距，γ₁₂为两信号频率比值的平方，T为对流层延时修正值，ε_ij为观测噪声。

1.7 基于伪距观测的共视时间比对算法流程

1) 将导航电文、广播星历和接收机坐标位置作为算法的输入量，解码导航电文或卫星星历并计算卫星轨道位置；2)进行地球自转效应修正；3)使用接收机位置和修正后的卫星位置计算星地几何距离，并进行对流层延时计算；4)将上述计算结果通过互联网实现数据交换，最终利用GNSS双频接收机时间差观测方程计算实时时间差值; 5)对实时时间差值进行卡尔曼滤波和RTS事后处理获得事后时间比对结果。

2 时间差数据滤波与平滑 2.1 卡尔曼滤波

上述算法最终可以得到不同接收机时间差数据，但由于观测噪声的影响，这组数据的精度很难达到要求。卡尔曼滤波是一个对信号或数据进行处理和变换的过程，其最主要目的是去掉或削弱不想要的成分对估计值的影响，并增强想要成分的权重^[13]。卡尔曼滤波采用递推处理，用上个采样时刻的异地钟差估计值和当前时刻的异地钟差观测值，来估计当前时刻的异地钟差估计值，当前时刻以后的观测值不会对当前时刻的估计值产生任何影响，因而适合于实时的共视观测资料处理^[14]。因此对时间差数据进行卡尔曼滤波是提高时间比对精度的常用方法。对含噪声的钟差数据(观测量)进行卡尔曼滤波，可估计出准确的钟差。假设在k时刻钟差真值用x_k表示，它构成了状态变量X_k:

$ {\mathit{\boldsymbol{X}}_k} = \left( {{\mathit{\boldsymbol{x}}_k}} \right) $

(10)

卡尔曼滤波的状态方程为

$ {\mathit{\boldsymbol{X}}_k} = {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{k,k - 1}}{\mathit{\boldsymbol{X}}_{k - 1}} + {\mathit{\boldsymbol{W}}_{k - 1}} $

(11)

量测方程用1.5节中得到的时间差矩阵代替。在式(11)中，Φ_{k, k-1}为状态转移矩阵，取值为[1]，W_k-1为模型噪声。

卡尔曼滤波器的动态系统维数n、观测系统维数m均为1^[15-16]。

考虑接收机未能按照共视表规定时刻及时锁定卫星，甚至在共视表所规定的整个跟踪时间段内，一直未能锁定卫星，致使共视数据中缺少该记录。采用等间隔的卡尔曼滤波器，对于共视间歇和未能锁星成功的时间段，用前3个时刻点上钟差估计值$ \mathit{\boldsymbol{\hat x}}\left( {k - 3} \right)、\mathit{\boldsymbol{\hat x}}\left( {k - 2} \right)、\mathit{\boldsymbol{\hat x}}\left( {k - 1} \right) $进行二阶多项式外插，作为当前时刻的观测值x(k)，使卡尔曼滤波继续进行。二阶多项式外插考虑到了两站原子钟的频差和相对漂移的影响，可以对间歇时间段内两站钟差变化量进行较为准确的估计，从而保证了卡尔曼滤波器性能。

在共视资料预处理后，对含噪声的异地钟差数据序列进行卡尔曼滤波，其算法过程为：

1) 状态变量X_k(只含钟差真值一个分量)与其卡尔曼估计$ {\mathit{\boldsymbol{\hat X}}_k} $之间的均方误差矩阵称为估计误差的协方差矩阵，用C_k表示。给定C₀一个初值，根据下式

$ {\mathit{\boldsymbol{P}}_{k + 1}} = {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{k + 1,k}}{\mathit{\boldsymbol{C}}_k}\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{k + 1,k}^{\rm{T}} + {\mathit{\boldsymbol{Q}}_k} $

(12)

可以计算出P₁。其中，P_k为状态变量X_k与其在无观测噪声与模型噪声条件下的估计$ {\mathit{\boldsymbol{\hat X'}}_k} $之间的均方误差阵，Q_k为1×1阶模型噪声叭的协方差阵。

2) 得到P₁后，根据卡尔曼增益矩阵G_k的表示式：

$ {\mathit{\boldsymbol{G}}_k} = {\mathit{\boldsymbol{P}}_k}\mathit{\boldsymbol{H}}_k^{\rm{T}}{\left[ {{\mathit{\boldsymbol{H}}_k}{\mathit{\boldsymbol{P}}_k}\mathit{\boldsymbol{H}}_k^{\rm{T}} + {\mathit{\boldsymbol{R}}_k}} \right]^{ - 1}} $

(13)

求得G₁，其中R_k为1×1阶观测噪声V_k的协方差阵。

3) 根据下式

$ {{\mathit{\boldsymbol{\hat X}}}_k} = {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{k,k - 1}}{{\mathit{\boldsymbol{\hat X}}}_{k - 1}} + {\mathit{\boldsymbol{G}}_k}\left[ {{\mathit{\boldsymbol{Y}}_k} - {\mathit{\boldsymbol{H}}_k}{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{k,k - 1}}{{\mathit{\boldsymbol{\hat X}}}_{k - 1}}} \right] $

(14)

可得k=1时刻的状态变量估计值$ {\mathit{\boldsymbol{\hat X}}_1} $，(即k=1时刻钟差卡尔曼估计值)。

4) 将P₁代入下式

$ {\mathit{\boldsymbol{C}}_k} = \left( {{\bf{I}} - {\mathit{\boldsymbol{G}}_k}{\mathit{\boldsymbol{H}}_k}} \right){\mathit{\boldsymbol{P}}_k} $

(15)

可求得k=1时刻的1×1阶估计误差的协方差阵C₁。然后，进入下次循环。

在近实时共视中，应该对异地钟差X₀有初步的估计，用这个值对滤波器进行初始化，会使滤波器的收敛速度加快。此时，估计误差的协方差矩阵初值取为

$ {\mathit{\boldsymbol{C}}_0} = {\rm{E}}\left[ {\left( {{\mathit{\boldsymbol{X}}_0} - {{\mathit{\boldsymbol{\dot X}}}_0}} \right){{\left( {{\mathit{\boldsymbol{X}}_0}{{\mathit{\boldsymbol{\hat X}}}_0}} \right)}^{\rm{T}}}} \right] = {\mathop{\rm var}} \left( {{\mathit{\boldsymbol{X}}_0}} \right) $

(16)

2.2 RTS事后算法

RTS固定区间最优平滑算法是在卡尔曼滤波的基础上，利用整个时间间隔内所有量测数据得到状态的最小方差估计，可以获得比卡尔曼滤波精度更高的融合结果。平滑解算过程相对于滤波过程是逆向的。因此，RTS固定区间平滑算法在传递对准精度评估等侧重于初始状态获取的应用中，最终平滑值的读取方式与普通前向滤波器估计值的读取方式相反。平滑过程首先是进行前向卡尔曼滤波，获得滤波估计值，然后经过一个反向平滑过程，进而得到平滑估计值。因此，平滑解算需要在滤波过程中实时存储数据，所存储的数据为4个矩阵，分别为估计值$ \mathit{\boldsymbol{\hat X}}_{\rm{K}}^f $、一步转移阵ϕ_{k, k-1}^f、估计的均方差阵P_k^f、以及一步预测均方差阵P_{k, k-1}^f。

平滑公式为

$ \begin{array}{*{20}{c}} {\left\{ \begin{array}{l} {\mathit{\boldsymbol{K}}_{s,k}} = {\mathit{\boldsymbol{P}}_{f,k}}\phi _{k + 1,k}^{\rm{T}}\mathit{\boldsymbol{P}}_{f,k + 1/k}^{ - 1}\\ {{\mathit{\boldsymbol{\hat X}}}_{s,k}} = {{\mathit{\boldsymbol{\hat X}}}_{f,k}} + {\mathit{\boldsymbol{K}}_{s,k}}\left( {{{\mathit{\boldsymbol{\hat X}}}_{s,k + 1}} - {{\mathit{\boldsymbol{\hat X}}}_{f,k + 1/k}}} \right)\\ {\mathit{\boldsymbol{P}}_{s,k}} = {\mathit{\boldsymbol{P}}_{f,k}} - {\mathit{\boldsymbol{K}}_{s,k}}\left( {{\mathit{\boldsymbol{P}}_{f,k + 1/k}} - {\mathit{\boldsymbol{P}}_{s,k + 1}}} \right)\mathit{\boldsymbol{K}}_{s,k}^{\rm{T}} \end{array} \right.,}\\ {k = N - 1,N - 2, \cdots ,2,1,0} \end{array} $

(17)

式中：$ {\mathit{\boldsymbol{P}}_{s, N}} = {\mathit{\boldsymbol{P}}_{f, N}}, {\mathit{\boldsymbol{\hat X}}_{s, N}} = {\mathit{\boldsymbol{\hat X}}_{f, N}}。{\mathit{\boldsymbol{\hat X}}_{s, k}} $为RTS平滑值，K_{s, k}为滤波增益，P_{s, k}为协方差阵。

3 数值仿真与结果分析

本文以位于北京和天津的两个接收机的时间比对过程为数学仿真算例，采用基于伪距观测的单星双频共视时间比对算法，同时观测同一颗COMPASS卫星播发的B₁信道，使用Matlab软件对时间比对过程进行编程仿真。接收机在大地坐标系中的坐标分别为(40, 116, 0)和(39, 117, 0)，表格1为仿真参数。

图 2依次为状态方程系数和观测方程系数为1、0.5、2、0.3、3时的仿真结果。当卡尔曼滤波中的状态方程系数和观测方程系数均为1时，滤波效果最好。从图中可以看出，滤波前的时间差观测值的精度为12 ns，滤波后的精度为3 ns，RTS处理后的精度为1 ns。图 2(a)中可以看出卡尔曼滤波结果收敛到0.4 ns以内，由于图像波动过小，无法正确反映出观测值的变化趋势。从图 2(b)中可以看出卡尔曼滤波结果波动幅度明显增大，甚至某些滤波值已经超过了观测值，因此导致滤波效果不理想。从图 2(c)、(d)中可以看出，状态方程系数影响卡尔曼滤波的震动幅度，当系数增大时，震动幅度迅速增大，反之亦然；当系数取值为1时，滤波结果既可以正确的反映观测值的变化趋势，也可以取得较好精度。图 2(d)中可以看出卡尔曼滤波结果图像振幅变大，RTS处理结果值变得更为平滑。但滤波精度只能达到3 ns，因此0.3不是最理想的观测方程系数。从图 2(e)中可以看出卡尔曼滤波结果图像振幅变小，RTS处理结果不是很理想，虽然滤波精度达到1.5 ns，但RTS处理结果与卡尔曼滤波结果区别不够明显，曲线不够平滑。因此3也不是最理想的观测方程系数。从图 2(d)、(e)中可以看出观测方程系数的变化对卡尔曼滤波结果和RTS处理结果都有影响，当观测方程系数增大时，卡尔曼滤波结果的振幅减小，精度提高，但RTS处理的结果平滑度降低，精度提升不明显。反之，当观测方程系数减小时，卡尔曼滤波结果振幅增加，精度下降，但RTS处理后的图像平滑度增加。

综上所述，当状态方程和观测方程的系数均为1时，滤波结果最为理想。能够正确的表现出时间差观测值的变化趋势，同时，精度可以达到1 ns以内。因此，本文所述的这种基于伪距观测的单星双频时间比对算法具有成本较低，易于实现，精度较高三个特点，适合在可以接收卫星双频导航电文的情况下，确定两个或多个相距较远的接收机之间的时间差时使用。

4 结论

1) 对流层、电离层和地球自转效应等链路误差对共视时间比对精度影响较大，为了提高共视精度，在实时比对算法和事后比对算法中均需要考虑上述链路误差。

2) 卡尔曼滤波和RTS事后平滑相结合的算法可以有效的提高事后共视室间比对精度，可以使GPS、COMPASS、Galileo和GLONASS等GNSS系统双频接收机的事后共视时间比对精度达到2 ns。在研究过程中发现，链路误差模型采用的是经验模型，精度还有提高的空间。

未来可以在本文算法的基础上，继续研究精度更高的链路误差，事后滤波及平滑等算法，进一步挺高共视时间比对精度。