北斗三号卫星共视时间比对性能分析

王威雄，董绍武，武文俊，王 翔，郭 栋，广 伟

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

0 引 言

时间是一个国家的重要战略参数资源，在经济建设，军事国防和各类科学研究都有广泛应用。国际标准时间协调世界时(Coordinated universal time, UTC)是在国际地球自转服务组织(International earth rotation service, IERS)协助下由国际权度局(Bureau international des poids et mesures, BIPM)利用分布在全球80多个守时实验室的500多台原子钟统一归算产生[1]，这些原子钟通过国际时间比对网参与UTC的计算。因此，远距离时间比对是UTC产生中的重要环节。现国际上主要采用基于全球导航卫星系统(Global navigation satellite system, GNSS)的时间比对和地球同步轨道卫星(Geosy-nchronous earth orbit, GEO)的双向时间比对(Two-way satellite time and frequency transfer, TWSTFT)两种手段，但由于TWSTFT设备昂贵且受基线长度限制，只有少数十几个实验室才有TWSTFT链路，而GNSS时间比对仍是国际上广泛采用的技术[2]。GPS共视法(Common view, CV)从20世纪80年代就被BIPM用于UTC的计算。为提高时间比对的可靠性、稳定性及准确性，BIPM鼓励采用多手段冗余时间比对技术[3]。2009年俄罗斯GLONASS时间比对第一次被引入UTC的计算，2017年国际时间频率咨询委员会(Consultative committee for time and frequency, CCTF)国际原子时实验室工作组会议也提出应积极推动北斗参与UTC的计算[4]。

北斗卫星导航系统是我国独立自主开发的卫星导航系统。近年来，国内外学者针对北斗二号卫星开展了诸多时间比对研究。例如，Huang等[5]在欧洲内部对北斗GEO、倾斜地球同步轨道(Inclined Geosynchronous orbit, IGSO)以及中圆地球轨道(Medium earth orbit, MEO)不同星座的共视时间比对进行了分析[5]；Guang等[6]进一步开展了亚欧长基线北斗共视时间比对试验；而Andreas等[7]也基于德国物理技术研究院(Physikalisch-technische bundesanstalt, PTB)的国家时间基准系统对北斗信号以及其亚欧共视比对结果进行了有效评估。上述诸多研究表明：北斗二号卫星时间比对性能稳定，但由于欧洲可视卫星数较少，比对性能仍有提升空间。文献[6-7]中也都提到期望对北斗三号卫星进行更深入的研究。

目前北斗三号卫星正处于全球系统的组网阶段，现已发射了20颗组网卫星，能够向“一带一路”国家和地区提供基本导航服务。全球系统预计在2020年底建成，届时整个星座将由3颗GEO卫星、3颗IGSO卫星和24颗MEO卫星组成，为全球提供导航、测速及授时服务。北斗三号卫星配备了高性能的铷原子钟和氢原子钟，铷原子钟天稳定度为E-14量级，氢原子钟天稳定度为E-15量级，比北斗二号卫星星载钟的稳定度提高了一个数量级[8-9]。我国标准时间由中科院国家授时中心(National time service center, NTSC)负责产生并对外发播，而北斗时(BDS time, BDT)通过UTC(NTSC)与UTC取得联系[10]。本文依据CCTF发布的全球卫星导航时间比对标准(Common GNSS generic time transfer standard, CGGTTS)，利用目前NTSC和捷克光电研究院(Institute of photonics and electronics, TP)可观测的北斗三号卫星对亚欧长基线共视时间比对性能进行了初步计算评估。

1 北斗共视比对算法实现

地面上任意两个守时实验室进行北斗共视时间比对是由分别位于两地的多频多模接收机在同一时刻同时观测一颗或多颗北斗卫星，通过伪距测量得到本地时与北斗系统时BDT之间的偏差，该偏差包括卫星、信号传播路径和接收机端引入的各项误差，两地通过数据交换扣除各项误差后求差，就可以获得两实验室之间的时差[11-13]，比对原理如图1所示。

图1 北斗共视时间比对原理Fig.1 The principle of common view time comparison by BeiDou

设A地的本地时间为tA，B地的本地时间为tB，ΔtABDT是A地与BDT的差，ΔtBBDT是B地与BDT的差，有：

ΔtABDT=tA-tBDT

(1)

ΔtBBDT=tB-tBDT

(2)

式(1)减去式(2)可得两地的时间差：

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

tBDT)=tA-tB=ΔtAB

(3)

为标准化数据处理流程，国际时间频率咨询委员会在1993年规定了GPS时间比对的格式标准GGTTS V01，1998年该标准又加入了俄罗斯GLONASS卫星导航系统，并命名为CGGTTS V02，随着欧洲Galileo和中国北斗等卫星导航系统的不断发展，2015年CCTF的GNSS时间传递组再次对时间比对标准进行了扩展，并更名为CGGTTS V2E，该版本与旧版本数据格式兼容，可在不同版本间完成共视时间比对[13]，其中几项重要的参数信息如表1所示。

表1 CGGTTS标准文件的参数信息Table 1 Parameter information of CGGTTS file

在CGGTTS时间比对标准框架下，其处理算法是连续跟踪一个或多个卫星并采集13 min的伪码观测数据，观测值1 s一个共780个，每15个测量值作一个二次多项式拟合，然后对拟合结果作各项误差修正，取其拟合的中间值。13 min的跟踪周期内有52个拟合中间值，最后将52个拟合中间值进行线性一次拟合，取中间值是本地时间尺度和卫星系统时的偏差[14]。在误差项处理中，电离层时延和对流层时延分别采用双频无电离层组合模型和标准的NATO模型来修正，地球自转效应使用Sagnac效应算法来修正，设备时延经校准后和测量的线缆时延和参考时延一起写入文件表头用于共视时间比对的计算[15-16]。

2 数据平滑滤波与异常值剔除

在共视比对中，由于卫星钟的不确定性，以及环境干扰及测量噪声因素的存在，比对结果中会出现一些异常值和随机噪声，本文采用3σ准则进行异常值剔除和Vondrak平滑法进行滤波处理。

2.1 Vondrak平滑滤波

Vondrak滤波是由捷克天文学家J.Vondrak在whittaker修匀的基础上提出的一种滤波方法，其既适用于等间隔的观测数据又适用于非等间隔的观测数据，且可以在观测数据变化规律未知的情况下对观测数据进行有效的平滑，其基本准则是[17]：

(4)

式中：pi为观测值的权，yi为观测值，y′i为待求的平滑值，ε=1/λ是一个给定的无量纲正数，定义为平滑因子。F为平滑值与观测值的拟合度，S为待求平滑曲线的平滑度。Vondrak滤波实质上就是在观测值的绝对拟合和绝对平滑之间选择一条折衷的曲线。ε作为调节折衷程度的参数，ε越小，曲线平滑程度越强，反之，平滑程度越弱[18]。

但在实际应用中，对于相同精度的观测数据，若采用同一平滑因子对不同采样间隔和不同长度的数据进行平滑，则得到的平滑曲线具有不同的平滑程度，必须选用不同的平滑因子才能得到同样平滑程度的曲线，这就给Vondrak的实际运用带来了不便。为此，Vondrak在1976年对其基本准则进行了改进，将式(4)中的拟合度和平滑度分别用平均值来代替：

(5)

根据式(5)得到的平滑方法，对于相同精度的观测数据，无论取样间隔和数据长度如何，选用同一平滑因子就可得到具有同样平滑程度的结果。

Vondrak滤波的关键是平滑因子的选取，目前常用的选取方法有观测误差法，频率响应法和交叉证认法等[19]。本文根据原始共视比对数据的幅值频谱图，采用频率响应法计算恰当的平滑因子作为最后确定的最佳平滑因子。频率响应函数的解析表达式为：

(6)

式中：f为频率，ε为平滑因子，A为在不同频率和平滑因子下的频率响应函数。

对于滤波结果的评价标准采用互相关系数(R)和滤波后噪声的均方根误差(NRMSE)来评估。

1)互相关系数R的计算公式为：

(7)

式中：cov(yi,y′i)是观测值yi与滤波值y′i的协方差，σy为观测值yi的标准差，σy′为滤波值y′i的标准差。R越大，说明滤波值与原始观测值相关性越大，即滤波值与原始值越接近。

2)滤波后噪声的均方根误差NRMSE计算公式为：

(8)

式中：yi为观测值，y′i为滤波值，n为观测值个数。NRMSE越小，说明滤波值与原始观测值拟合性越好，测量精度越高。

2.2 异常值剔除

3σ准则是随机误差服从正态分布情况下的异常值判别方法。设服从正态分布的一维随机变量X的概率密度为：

(9)

式中：x∈(-∞,+∞)，μ为X的数学期望，σ2为X的方差。

正态随机变量X出现在给定区间(μ-kσ，μ-kσ)内的概率(k为正数)为：

(10)

(11)

查概率积分表可得，当k=3时，正态随机变量X出现在区间μ±3σ范围内的概率约为99.7%。利用3σ准则剔除共视比对结果中异常值的流程图如图2所示。

图2 3σ准则异常值剔除流程图Fig.2 Flow chart of 3σ criterion outlier elimination

3 实验与结果分析

选取2019年3月28日(MJD 58570)2019年4月27日(MJD 58600)间NTSC和TP的多频多模接收机观测获得的BeiDou和GPS CGGTTS数据。该试验中，每台接收机的时间频率参考都为国家时间基准UTC(k)(k指不同的守时实验室)，参与此次时间比对实验的接收机配置信息如表2所示。

表2 接收机配置信息Table 2 Receiver configuration information

3.1 零基线共钟比对

为测试接收机的稳定性和对北斗三号卫星测距码的噪声水平进行评估，需要对接收机进行零基线共钟时间比对(Common clock difference, CCD)实验。利用共视法对捷克光电研究院的两台同类型接收机TP01和TP02开展了北斗三号卫星CCD实验，并和GPS CCD的结果进行了比较。利用比对结果的标准差(Standard deviation, STD)对噪声水平进行评估，比对结果如图3和图4所示。由于NTSC目前暂未有两台可接收北斗三号卫星信号的接收机，所以未进行CCD实验，但参与本次比对的接收机与TP相同，故认为TP的CCD结果也适用于NTSC。

图3 北斗三号卫星零基线共钟比对结果Fig.3 The CCD results of BDS-3 satellite

图4 GPS零基线共钟比对结果Fig.4 The CCD results of GPS

由图3和图4可以看出，TP01和TP02接收机基于北斗三号卫星的零基线时间比对结果稳定，其STD为0.74 ns，与GPS零基线共钟比对结果的STD处于同一量级，但由于北斗三号卫星正处于全球组网阶段，观测到的BDS-3卫星较GPS卫星数少，所以BDS-3卫星的测量噪声水平略差于GPS，但可以用于远距离的高精度时间比对。

3.2 单站观测数据分析

利用CGGTTS标准共视文件在比对时间段内对NTSC和TP的观测数据分别进行分析。两个实验室观测到不同卫星数的观测量占整个观测量的百分比统计如图5～6所示。

图5 NTSC的可视卫星数统计Fig.5 The visual satellite number statistics of NTSC

图6 TP的可视卫星数统计Fig.6 The visual satellite number statistics of TP

从图5～6可以看出，由于北斗二号卫星是区域卫星导航系统，亚太地区覆盖范围比欧洲区域好，因此NTSC观测到的北斗二号卫星数目比TP多，北斗三号卫星是全球卫星导航系统，目前在轨服务的卫星和GPS星座卫星相似都为MEO卫星，所以北斗三号卫星和GPS在亚太和欧洲地区都有相似的覆盖情况，在NTSC和TP观测到的北斗三号卫星数量大多都分布在2～5颗，GPS卫星数量大多为5～9颗。

相关研究发现测量信号的码伪距多路径噪声与卫星高度角有关，且高度角越高，多路径噪声越小[6]。图7和图8分别为观测期间(MJD 58570-58571)在NTSC和TP不同高度角区间的卫星观测量统计情况，由于北斗二号卫星的GEO卫星(C01～C05)在NTSC的可观测数多于TP，因此在频率分布直方图中未统计GEO卫星。

图7 NTSC观测卫星的高度角统计情况Fig.7 Altitude statistics of NTSC visual satellites

图8 TP观测卫星的高度角统计情况Fig.8 Altitude statistics of TP visual satellites

从图7、8可以得出，北斗二号卫星在NTSC的观测高度角约有70%都在50°以上，信号覆盖情况优于TP，数据的观测质量也更好。而NTSC和TP的北斗三号卫星和GPS卫星在不同高度角区间内的观测数量在本地和异地都比较相近。

3.3 共视比对结果分析

在GNSS共视时间比对时，较多的可用卫星数可以平均出更好的时间比对结果。图9给出了NTSC和TP在相同时刻能够共视到的卫星数。

图9 NTSC和TP共视卫星数Fig.9 The satellites number of CV between NTSC and TP

从图9得知，NTSC和TP共视到的GPS卫星数大多为2～3颗，北斗二号卫星数大多为2～4颗，北斗三号卫星数最少，大多为1～2颗。

基于北斗三号卫星，依据GNSS标准共视数据处理规范并利用式(3)得到了几何基线为7500 km上的NTSC-TP共视时间比对结果。由于卫星或者观测条件的限制，GNSS时间比对还存在不同程度的粗差。因此，根据图2中的流程先后对原始计算结果利用频率响应法进行了Vondrak滤波平滑和3σ准则粗差剔除。由于Vondrak滤波实质是低通滤波，这里选用的频率响应值A=0.2，原始共视比对结果频谱图低频区域中较高频率幅值最大的频率为f，依据式(6)可得平滑因子ε=3493000。对原始比对结果进行滤波平滑和粗差剔除后的三种链路的共视比对结果如图10所示。为便于分析，图10中对GPS和北斗二号卫星共视结果进行了常数平移。图11～13给出了三种链路时间比对结果滤波前后的幅值频谱变化情况。

图10 NTSC-TP链路共视比对结果Fig.10 The CV results of NTSC-TP link

表3采用原始结果和滤波结果的互相关系数和滤波后噪声的均方根误差对滤波结果进行了评估。

表3 三种共视链路滤波结果对比Table 3 Comparison of filtering results of three CV links

从图10和表3可以看出，GPS共视链路、北斗二号卫星共视链路和北斗三号卫星共视链路的Vondrak平滑值与时间比对计算值有很好的吻合性，可以反映出原始共视比对数据的具体特征。二者具有很强的相关性，滤波值和原始共视结果的相关系数可以达到0.985，0.957和0.968。从图11～13可知，原始比对结果中含有很强的高频噪声，需对噪声加以抑制或消除。从Vondrak滤波后的频谱图可见，经滤波后高频噪声得到了很好的抑制。因此，Vondrak滤波可以对共视比对数据进行有效平滑。滤波后三种共视比对链路噪声的均方根误差分别为0.839 ns，1.417 ns和1.161 ns。由图9，图10及表3可知，在目前北斗三号卫星比北斗二号卫星共视可视卫星少的情况下，其共视比对精度较北斗二号卫星有了明显提高，提高幅度约19%。北斗三号卫星正处于全球组网阶段，目前其共视比对精度较GPS略低。

图14为原始共视结果和滤波值的残差统计图及正态拟合曲线。由图可以看出，三种比对链路实测的共视比对结果相对于平滑结果的残差分布服从正态分布。一方面证明本文通过3σ准则剔除异常值的恰当性，另一方面也说明利用Vondrak滤波主要滤掉的是比对链路中的随机噪声，平滑后的曲线能够保留原始比对结果的固有信息，具有很高的保真度。

三种远程共视比对链路的时间方差和修正阿伦方差如图15、16所示。

从图15、16可以看出，三种共视链路1 d的时间和频率稳定度都在1 ns量级和10-14量级。当平均时间小于10000 s时，GPS共视比对链路的稳定度最高，北斗三号卫星优于北斗二号卫星。当平均时间大于10000 s时，三者的稳定度指标基本相当。

图11 GPS共视比对结果的频谱图Fig.11 The spectrum of GPS CV results

图12 北斗二号卫星共视比对结果的频谱图Fig.12 The spectrum of BDS-2 satellite CV results

图13 北斗三号卫星共视比对结果的频谱图Fig.13 The spectrum of BDS-3 satellite CV results

图14 原始结果与滤波值残差分布Fig.14 Residual distribution of original results and filtering results

图15 不同链路的时间方差Fig.15 Time variance of different links

图16 不同链路的修正阿伦方差Fig.16 Modified Allan variance of different links

4 结 论

利用目前中捷两地可接收的北斗三号卫星完成了NTSC和TP之间亚欧长基线共视时间比对初步试验与评估，对于北斗系统早日正式纳入UTC的归算提供了技术参考。基于北斗三号卫星完成了零基线共钟时间比对，利用CGGTTS数据分析了单站可视卫星数及高度角情况，通过Vondrak滤波对共视比对结果进行滤波平滑与粗差剔除，并和北斗二号系统及GPS系统做了比较验证，并得出以下结论：

1)北斗三号卫星的零基线共钟时间比对结果连续稳定，其标准差为0.74 ns，可以用于高精度时间比对。

2)Vondrak滤波可以对北斗共视结果进行有效平滑，滤波值具有很高的保真度。

3)在当前中捷之间北斗三号卫星共视可视卫星数比北斗二号卫星共视可视卫星数少一半的情况下，北斗三号卫星的共视比对精度为1.16 ns，北斗二号卫星的共视比对精度为1.42 ns，其比对精度较北斗二号卫星有明显提升，提升幅度约19%。

4)北斗三号卫星，北斗二号卫星和GPS共视链路1 d的时间和频率稳定度都在1 ns量级和10-14量级。10000 s以内的链路稳定度GPS共视最高，北斗三号卫星优于北斗二号卫星。当平均时间大于10000 s时，三者的稳定度指标基本相当。

致 谢

感谢捷克科学院光电研究院提供的北斗和GPS RINE以及CGGTTS观测文件。