Vondrak-Cepek组合滤波在北斗共视和卫星双向时间比对融合中的应用*

王威雄，董绍武，武文俊，郭 栋，王 翔，高 喆

(1. 中国科学院国家授时中心， 陕西 西安 710600； 2. 中国科学院时间频率基准重点实验室， 陕西 西安 710600； 3. 中国科学院大学 天文与空间科学学院， 北京 100049)

高精度时间比对技术是国际标准时间协调世界时(Coordinated Universal Time, UTC)产生的必要手段。基于全球卫星导航系统(Global Navigation Satellite System, GNSS)的共视(Common View, CV)、全视和精密单点定位(Precise Point Positioning, PPP)时间比对技术以及另一种完全独立的卫星双向时间比对(Two-Way Satellite Time and Frequency Transfer, TWSTFT)技术是目前国际上主要采用的几种时间比对手段[1]。TWSTFT从1999年就成为UTC计算的常规比对手段，且目前传递了超过60%的原子钟数据以及大部分用于UTC计算的一级频标的数据[2]。随着光钟、喷泉钟等更高精度的频率基准的发展，对远距离时间比对的精度也提出了更高的要求。TWSTFT平均时间为1 d内的时间稳定度优于1 ns且其长期稳定度较好，但其采样率较低且在大部分链路中会受到周日效应的影响，导致其没有优良的短期稳定度[3]。从2016年2月起，国际权度局(Bureau International des Poids et Mesures, BIPM)组织在欧亚、欧美以及欧洲内部的守时实验室间开展基于软件接收机(Software-Defined Radio, SDR)的TWSTFT试点研究，降低了硬件TWSTFT中的周日效应，并于2020年3月首次正式将SDR TWSTFT用于UTC的计算，但SDR TWSTFT结果中仍有残余的周日效应存在[4-5]。

随着我国北斗等多个全球卫星导航系统的发展，时间比对技术呈现多元化的趋势。为进一步提高UTC时间比对链路的准确性、稳定性和鲁棒性，在2017年召开的第21届CCTF大会上提出充分利用冗余链路来提高当前时间比对链路稳定度并改善TWSTFT中周日效应问题的议案。Vondrak-Cepek(简称V-C)组合滤波是基于最小二乘原理的Vondrak平滑方法的进一步变形，该方法刚被提出时主要用于测量地球自转参数时GPS法和甚长基线干涉测量法(Very Long Baseline Interferometry, VLBI)的融合[6]。有关学者对该方法在时频领域的应用也进行了研究，Jiang等通过该方法将利用载波相位数据的GPS PPP时间比对和TWSTFT数据进行了融合，提高了时间比对链路的短期稳定度[7]；姜萌等采用V-C法对短期稳定度较好的氢原子钟数据和长期稳定度较好的铯原子钟数据在时间尺度方面进行了融合，获得了理想的氢铯融合时间尺度[8]。因此，本文基于目前中科院国家授时中心(National Time Service Center, NTSC)参与建立的卫星双向时间比对网和北斗链路，分别在长短基线上利用Vondrak-Cepek方法对没有周日效应现象的北斗CV和TWSTFT数据进行有效融合，通过时间偏差、频谱分析以及另一种独立的GPS PPP为参考对融合结果作内外符合评估。融合结果表明该方法可以明显改善TWSTFT的周日效应。

1 时间比对原理

1.1 SATRE和SDR TWSTFT原理

TWSTFT被广泛应用于UTC守时实验室的国际时间比对，其原理如图1所示。地面站1的本地主钟产生的信号输入到调制解调器与经过二进制相移键控方式调制的伪随机码(PseudoRandom Noise, PRN)信号同步，生成的70 MHz中频信号经上变频器变频为上行频率后，经功率放大器放大，由天线系统将信号发射到卫星；位于地面站2的天线系统接收到卫星转发的时间信号后，经下变频器将频率转换为70 MHz中频信号输入本地卫星时间和距离测量设备(SAtellite Time and Ranging Equipment, SATRE)的接收通道；将接收到的信号进行解调后得到站1和站2的信号传输时间，两站按照比对计划表进行常规比对，各自得到对方的比对文件后就可得到双向比对结果。SDR TWSTFT由软件编程来实现接收机功能，目前在卫星双向地面站，SDR TWSTFT的发射信号采用与硬件SATRE一致，接收通道与硬件SATRE相互独立，天线系统接收的信号经下变频器处理后输入SDR接收通道，通过信号转换器的模数转换、现场可编程门阵列(Field Programmable Gate Array, FPGA)的解调以及滤波后得到两地钟差。SATRE和SDR TWSTFT都可获得发射站与接收站之间的高精度钟差，两者相互冗余[9-10]。

图1 卫星双向时间比对原理Fig.1 Principle of TWSTFT

站1和站2时间尺度间的差值计算公式为：

T1-T2=0.5×(TTW1+TCALR1+TESDVAR1-TTW2-TCALR2-TESDVAR2)+TREFDLY1-TREFDLY2

(1)

式中：Tk为实验室1和实验室2的本地时间；TTWk为地面站1和地面站2调制解调器得到的信号传输时延值；TCALRk为经移动校准器标定的地面站校准时延值；TESDVARk为链路校准后的本地站设备变化时延值；TREFDLYk为实验室UTC(k)参考点到调制解调器发射端的时延值。

1.2 北斗共视比对原理

北斗共视时间比对是位于地球任何两个位置的守时实验室同时观测一颗或多颗北斗导航卫星，通过连接本地参考UTC(k)的接收机输出伪距测量值，建立观测方程分别得到本地时间与北斗系统时(BeiDou system Time, BDT)的时差，共视法可消除卫星钟误差和卫星轨道误差的影响，消除传播路径和接收机端的各项误差后就可以间接得到两地之间的时差，北斗共视比对原理如图2所示[11-12]。

图2 北斗共视时间比对原理Fig.2 Principle of BDS CV time transfer

假设A、B两地的本地时间分别为tA和tB，ΔtABDT、ΔtBBDT是A地、B地与BDT的差，有：

ΔtABDT=tA-tBDT

(2)

ΔtBBDT=tB-tBDT

(3)

式(2)减去式(3)可得两地的时间差ΔtAB：

ΔtAB=ΔtABDT-ΔtBBDT=(tA-tBDT)-(tB-tBDT)

=tA-tB

(4)

2 北斗共视和卫星双向融合算法

由于GNSS CV结果没有周日效应，因此采用Vondrak-Cepek组合滤波来对北斗共视和TWSTFT结果进行融合，用共视结果的一阶差分值对双向时间比对结果进行干预，来减小TWSTFT的周日效应，改善卫星双向时间比对的性能。

2.1 BDS CV和TWSTFT融合模型

假设TWSTFT在时间ti的观测量为y(t)，该时刻相对应的导数y′(t)定义为[y(ti+1)-y(ti)]/(ti+1-ti)，假设另一组独立的北斗共视观测量为Y(t)，其导数Y′(t)=[Y(tj+1)-Y(tj)]/(tj+1-tj)，这两组导数的物理意义是相同的，均表示两地钟差的变化率。当ti=tj时，y′(t)=Y′(t)，且

(5)

(6)

其中，MJD为简化儒略日。因为所选的双向时间比对间隔和北斗共视的时间间隔相同，所以式(6)下一时刻的TWSTFT观测值可表示为本时刻的TWSTFT观测值与该时段内北斗共视比对结果增量的和，总观测误差即TWSTFT实际观测值与通过共视表示的TWSTFT观测值在积分时间段内的差异。采用Vondrak-Cepek组合滤波方法建立平滑曲线的组合约束条件来实现对该观测误差的最小化。

2.2 Vondrak-Cepek组合滤波原理

1)曲线平滑度。

(7)

2)平滑曲线对观测值的逼真度。

(8)

3)平滑曲线对观测值一阶导数的逼真度。

(9)

Vondrak-Cepek组合滤波对数据融合的目的是找出输出的平滑序列yi在下列三种条件下的折中：

1) 曲线是平滑的(最小化S)。

2) 平滑序列应该接近观测函数值(最小化F)。

对三种条件进行最小化约束，有：

(10)

选用频率响应法来确定平滑因子，其公式为：

(11)

表1 频率响应值、平滑因子和周期间的关系

图3 P=1时频率响应值与平滑因子的关系Fig.3 Relationship between frequency response value and smoothing factor when P=1

3 结果分析

由于BIPM计算国际原子时(International Atomic Time, TAI)采用一个月的数据，因此分别选取几何基线长度约为7 700 km的NTSC和德国联邦物理技术研究院(Physikalisch-Technische Bundesanstalt, PTB)链路2019年10月至11月(MJD58775-MJD58804)以及几何基线长度约为900 km的NTSC和中国计量科学研究院(National Institute of Metrology, NIM)链路2020年1月至2月(MJD58852-MJD58881)的SATRE和SDR双向以及北斗共视数据进行分析。

3.1 BDS CV与TWSTFT融合结果

参与比对的卫星双向系统和北斗共视比对系统的参考同源且都为UTC(k)(k指不同的实验室)，实验室和设备信息见表2。

表 2 参与比对的守时实验室和设备信息

本次共视比对选择16 min无电离层组合的标准共视程序，组合信号频点为B1I和B2I。长基线NTSC-PTB链路的北斗共视结果及比对期间可共视卫星数统计如图4(a)和4(b)所示，短基线NTSC-NIM链路如图4(c)和4(d)所示。

(a) NTSC-PTB共视结果(a) NTSC-PTB CV results

(b) NTSC-PTB共视卫星数统计(b) NTSC-PTB CV satellite statistics

(c) NTSC-NIM共视结果(c) NTSC-NIM CV results

(d) NTSC-NIM共视卫星数统计(d) NTSC-NIM CV satellite statistics图4 不同基线北斗共视时间比对结果Fig.4 Results of BDS CV over different baselines

从图4可以看出，两条链路的北斗共视比对结果中都没有明显的周日效应现象，NTSC-PTB链路噪声较NTSC-NIM链路大。这是因为NTSC-PTB链路基线长，在同一时刻共视到的卫星数较短基线链路NTSC-NIM少，在共视比对中，较多的共视卫星数可以平均出更好的时间比对结果，但该结果不会影响双向和北斗共视的融合计算。图5为两条链路SATRE和SDR TWSTFT的原始结果。

(a) NTSC-PTB链路(a) NTSC-PTB link

(b) NTSC-NIM链路(b) NTSC-NIM link图5 SATRE和SDR双向比对结果Fig.5 Results of SATRE and SDR TWSTFT

由图5可知，无论是长基线链路还是短基线链路，SATRE和SDR TWSTFT的比对结果中都会看到明显的周日效应现象。为对双向结果进行改善，利用第2节的V-C组合滤波方法分别对NTSC-PTB长基线以及NTSC-NIM短基线上的SATRE和SDR TWSTFT比对结果与北斗共视结果进行融合，由于GPS PPP是与TWSTFT相独立的时间比对技术且其基本不受周日效应影响，因此以GPS PPP时间比对结果为参考，结果如图6和图7所示。

(a) SATRE双向链路不同方法时间比对结果(a) Time transfer solutions of different methods for SATRE TW link

(b) SDR双向链路不同方法时间比对结果(b) Time transfer solutions of different methods for SDR TW link图6 NTSC-PTB链路融合结果Fig.6 Combined smoothing results of NTSC-PTB link

(a) SATRE双向链路不同方法时间比对结果(a) Time transfer solutions of different methods for SATRE TW link

(b) SDR双向链路不同方法时间比对结果(b) Time transfer solutions of different methods for SDR TW link图7 NTSC-NIM链路融合结果Fig.7 Combined smoothing results of NTSC-NIM link

从图6和图7的融合结果可知，SATRE和SDR TWSTFT经V-C组合滤波后周日效应明显改善，且融合结果的噪声水平也明显减小，这与V-C本身的滤波特性有关。另外，融合结果与GPS PPP参考链路的结果基本一致，保证了融合结果的可靠性。

3.2 内符合评估

为对融合结果的稳定度及周日效应改善程度进行量化，采用时间偏差σx评估融合前后时间比对链路的稳定度，σx越小，链路稳定度越高。用式(12)的增益因子c表示改善程度的大小。

(12)

式中，σx(y1)和σx(y2)表示不同时间比对链路的时间偏差。这里用SATRE & BDS CV表示SATRE双向和北斗共视融合，用SDR & BDS CV表示SDR双向和北斗共视融合，下同。融合前后的时间偏差如图8所示。

(a) NTSC-PTB不同链路的时间偏差(a) Time deviation of different links in NTSC-PTB

(b) NTSC-NIM不同链路的时间偏差(b) Time deviation of different links in NTSC-NIM图8 融合前后时间偏差比较Fig.8 Comparison of time deviation before and after combined smoothing

从图8可以看出，在时间偏差指标上，长短基线的SATRE TWSTFT和SDR TWSTFT在平均时间为1 d时有明显的“凸起”现象，说明其周日效应明显，经过V-C组合滤波的融合结果周日效应基本消失，且链路的短期稳定度大大提高，这一方面与不同链路的融合有关，另一方面也与V-C滤波方法自身的平滑性有关。用平均时间1 d的时间偏差稳定度来评估融合方法对原始双向的改善程度，结果见表3。从表3可知：在长基线NTSC-PTB链路的天稳指标上，融合结果对SATRE TWSTFT的增益因子为1.85，对SDR TWSTFT的增益因子为1.81；在短基线NTSC-NIM链路的天稳指标上，融合结果对SATRE TWSTFT的增益因子为1.69，对SDR TWSTFT的增益因子为1.59。

采用快速傅里叶变换后频率分量的幅值对融合前后双向链路的周日效应改善进行标定。用式(13)的增益因子d表示改善程度的大小。

(13)

式中，A(y1)和A(y2)表示融合前后时间比对链路的24h频率分量幅值。NTSC-PTB双向链路和北斗共视融合前后的频谱图如图9所示, 采取与NTSC-PTB相同的方法求得NTSC-NIM的幅值频谱，各链路融合前后的24 h分量幅值(即周日效应量)统计结果见表4。

表3 NTSC-PTB和NTSC-NIM链路时间偏差改善统计

从图9和表4融合前后24 h分量的幅值变化可以看出，经过V-C融合后，两条链路SATRE及SDR TWSTFT中的周日效应都基本消失。对于亚欧NTSC-PTB链路，V-C融合结果24 h分量对SATRE TWSTFT的增益因子为4.30，对SDR TWSTFT的增益因子为2.50；对于亚洲内部的NTSC-NIM链路，V-C融合结果24 h分量对SATRE TWSTFT的增益因子为2.77，对SDR TWSTFT的增益因子为4.85。通过以上分析可知，利用V-C对双向和北斗共视融合后，长短基线上的SATRE和SDR TWSTFT的周日效应基本消失，且其1 d内的短期稳定度也得到提高。

(a) SATRE TWSTFT频谱分析(a) Spectral analysis in SATRE TWSTFT

(b) SATRE & BDS CV融合频谱分析(b) Spectral analysis in fusion of SATRE & BDS CV

(c) SDR TWSTFT频谱分析(c) Spectral analysis in SDR TWSTFT

(d) SDR & BDS CV融合频谱分析(d) Spectral analysis in fusion of SDR & BDS CV图9 NTSC-PTB链路融合前后频谱分析Fig.9 Spectral analysis of NTSC-PTB link before and after combined smoothing

表4 NTSC-PTB和NTSC-NIM链路24 h频率 幅值改善统计Tab.4 Gain in 24 h frequency amplitude of NTSC-PTB and NTSC-NIM links

3.3 外符合评估

外符合评估是指将融合结果与采用第三种方法获得的结果进行比较，GPS PPP时间比对是与TWSTFT完全独立的时间比对技术，其由于利用载波相位观测量，短期稳定度高且不受周日效应影响。因此，以同一基线的GPS PPP时间比对链路为参考，将TWSTFT与BDS CV经V-C融合后的结果与GPS PPP作差，得到双差(Double Clock Difference, DCD)结果后对融合结果的准确度与可靠度进行评估。图10为NTSC-PTB和NTSC-NIM两条链路SATRE和SDR双向与北斗共视融合结果分别与参考链路GPS PPP的DCD结果，表5为相关统计情况。

结合图10、表5可以看出： NTSC-PTB长基线链路SATRE & BDS CV融合结果与GPS PPP的差值绝对值最大为1.31 ns， SDR & BDS CV融合结果与GPS PPP的差值绝对值最大为1.16 ns； NTSC-NIM短基线链路SATRE和SDR与BDS CV融合结果与GPS PPP的差值绝对值的最大值分别为1.65 ns和1.64 ns，且融合前后DCD结果均值基本无变化，这说明经过V-C融合的结果与GPS PPP结果保持一致，V-C融合并没有改变原始双向比对的校准信息，保证了融合结果的准确度。而原始SATRE和SDR双向结果由于周日效应及链路不稳定性的存在，与GPS PPP的DCD最大值明显超过了链路允许的最大不确定度。与式(12)时间偏差指标增益因子计算相类似，从表5融合前后DCD结果的标准差比值计算可得，NTSC-PTB链路SATRE和SDR TWSTFT与GPS PPP DCD结果融合前后的标准差增益因子分别为2.71和1.94，NTSC-NIM链路融合前后的标准差增益因子分别为1.67和1.82，进一步证明融合结果可以减小TWSTFT结果中的周日效应并提高链路稳定度，与上一节内符合评估中得出的结论一致。

(a) NTSC-PTB链路(a) NTSC-PTB link

(b) NTSC-NIM链路(b) NTSC-NIM link图10 融合结果与GPS PPP DCD结果比较Fig.10 DCD results between combined solutions and GPS PPP

表5 不同链路的DCD结果

4 结论

采用Vondrak-Cepek组合滤波方法对没有周日效应现象的BDS CV和TWSTFT结果在不同基线长度上进行融合处理，通过融合前后各链路的时间偏差，频谱24 h分量幅值以及与GPS PPP的DCD等指标的分析，证明了V-C方法对TWSTFT和BDS CV融合的有效性。结果表明：融合结果在保证时间比对准确度的同时，可以基本消除SATRE和SDR TWSTFT中的周日效应，且其短期稳定度大幅提高。对于长基线NTSC-PTB链路，融合结果1 d的时间偏差稳定度对SATRE双向的增益因子为1.85，对SDR 双向的增益因子为1.81；对于短基线NTSC-NIM链路，融合结果对SATRE双向的增益因子为1.69，对SDR双向的增益因子为1.59。随着UTC时间比对手段的不断丰富，多链路多手段时间比对融合将是未来国际时间比对工作的重要议题。