全球导航卫星系统(global navigation satellite system, GNSS)旨在为全球用户提供全天候、全时段、高精度的导航、定位与授时服务^[1]。北斗卫星作为北斗卫星导航系统(BeiDou navigation satellite system，BDS)空间端的主要组成部分，为整个系统提供重要空间基准，而星载原子钟身为北斗卫星的重要载荷之一，维持着导航系统星上时间基准，决定了北斗卫星在轨寿命，其性能不仅会影响北斗卫星自主导航能力，还会对其服务能力产生较大影响^[2-3]。

GNSS星载原子钟在轨运行期间，会受到外部环境变化和设备老化等因素影响，导致其输出的钟差序列中存在粗差、钟跳等异常情况，而干净的卫星钟差数据是后续卫星钟性能研究与预报的前提^[4-6]，采用预处理后的钟差数据分析BDS卫星钟差特性可以取得较为可信的结果。常用的钟差数据质量控制方法主要有中位数法(median absolute deviation, MAD)和Baarda粗差探测法等。其中中位数法对粗差大小不敏感，简单有效，但很难探测出小的误差；Baarda法则是通过对验后残差进行处理，能够很好地识别隐藏在钟差数据中的小误差，但在进行多项式拟合钟差数据时，容易受到钟跳等异常值的影响无法准确获得拟合系数。针对常用钟差异常值探测方法的缺点，本文采用改进后的Baarda法对频率数据进行处理，能够准确探测出频率数据中的异常点，并采用多项式拟合值进行补充处理。

2015年, 我国发射了第一颗北斗全球卫星导航系统(简称北斗三号系统，BDS-3)试验卫星。截止到2018年4月，BDS-3卫星共在轨5颗试验卫星(C31~C35)和8颗组网卫星(C19、C20、C21、C22、C27、C28、C29、C30)。文献[7]指出我国北斗星载原子钟天稳定度要比GPS和GLONASS系统的星载原子钟低一个数量级，因此提高我国原子钟的物理性能是提高BDS定位精度一个亟待解决的问题。文献[8]利用北斗精密钟差评估了北斗频率稳定度，发现北斗原子钟的千秒稳定性、万秒稳定性和日稳定性分别能达到(3~4)×10^-13、(1~2)×10^-13、(1~2)×10^-14量级。文献[9-11]分析了BDS-3试验卫星的频率稳定度，结果表明BDS-3试验卫星频率稳定度比北斗二号系统(BDS-2)高一个量级。为分析BDS-3卫星钟差的性能，本文首先设计了3 d弧段的BDS-3精密卫星钟差解算实验；然后，根据北斗卫星钟差的噪声特性改进了Baarda法，为计算BDS-3卫星周期特性和频率稳定性指标提供了干净的钟差数据；最后，详细探讨了BDS-2/BDS-3之间、不同BDS卫星类型之间和不同星载原子钟类型之间的频率稳定性的差异。

1 BDS观测数据采集

BDS-2卫星精密钟差可以从国际GNSS服务组织(International GNSS Service, IGS)和国际GNSS监测评估系统(International GNSS Monitoring & Assessment System, iGMAS)等途径获得，但是现有的途径均不提供BDS-3精密钟差产品。

实验采用了iGMAS测站估计的3月1日—31日3 d弧段的BDS-2/BDS-3卫星钟差, 各卫星参数如表 1所示。

目前，已经建设好的iGMAS测站共有24个，其中能够接收到BDS-3信号的测站有17个，选择hmns和bjf1两个iGMAS测站分析跟踪BDS-3信号，结果如图 1所示。由图 1可以看出，iGMAS测站跟踪倾斜地球同步轨道(inclined geosynchronous orbit，IGSO)卫星的弧段长度要略长于跟踪中圆地球轨道(medium earth orbit, MEO)卫星的弧段长度，这也是估计得到的IGSO卫星钟差精度优于MEO卫星钟差精度的主要原因。图 2统计了iGMAS站接收BDS-3卫星信号的数据有效率，基本在50%左右。

图 3和表 2统计了BDS-3卫星钟差重叠弧段的均方根误差(root mean square, RMS)和标准差(standard deviation，STD)。从图 3和表 2中可以发现，STD和RMS的精度相近，这说明连续两次定轨过程中不存在明显的系统性误差，而且大部分钟差精度优于1 ns。其中IGSO卫星的钟差精度优于MEO卫星，主要原因可能是iGMAS测站分布不均匀，MEO卫星的观测弧段较少。

2 卫星钟差数据预处理

Baarda法能够有效地探测出钟差序列中的小粗差并将其剔除^[12]。为了避免粗差、钟跳等异常值对二次多项式模型参数求解的影响，并且考虑到频率数据较钟差相位数据更能突出异常点，本文采用抗差最小二乘估计模型参数对频率数据进行建模，对探测出的异常点采用多项式拟合值进行填充，并将频率数据恢复成钟差相位数据。具体模型如下：

$ x(t) = {a_0} + {a_1}t + {a_2}{t^2} + {e_t} $

(1)

$ {x^\prime }(t) = {a_1} + 2{a_2}t + {w_t} $

(2)

$ y(t)=x^{\prime}(t)=\frac{x(t+\tau)-x(t)}{\tau} $

(3)

式中，x(t)为钟差相位数据；y(t)为频率数据；a₀、a₁、a₂分别表示相位、钟速和钟漂；τ为采样间隔；e_t和w_t分别表示钟差和频率时间序列的拟合残差。

为了确定频率数据的权，本文选用中国矿业大学分析中心解算的C19和C34两颗卫星在2018年年积日第60~90天共31 d的BDS-2/BDS-3精密钟差数据，来分析BDS星载原子钟的噪声水平，采样间隔为300 s。采用直方图和分位(quantile-quantile, Q-Q)图对多项式拟合残差进行分析，检验拟合残差是否符合标准正态分布，结果分别如图 4、图 5所示。图 4显示了拟合残差的均值和标准差，从中可以看出拟合残差序列基本符合正态分布；图 5表明拟合残差序列在95%的置信区间内可以很好的拟合。从以上分析中可以得出：钟差残差序列基本符合正态分布。

为了直观地展示本文改进的钟差数据异常值探测方法，图 6给出了本文进行数据预处理的详细流程图, 首先根据式(3)，将读取的钟差数据转变为频率数据后，用线性模型对频率数据进行拟合得到钟速、钟漂和拟合残差，然后根据残差的标准差，归一化拟合残差得：

$ \tilde{v}_{i}=\frac{v_{i}}{\sigma_{v_{i}}} $

(4)

式中，$ {{\sigma _{{v_i}}}}$为残差v_i对应的标准差。若归一化残差大于阈值k(一般取3)，则剔除最大残差所对应的频率数据，并重新构建法方程，估计钟差参数得到残差, 直至没有残差大于阈值为止。

根据改进的异常值探测方法，对C34号卫星钟差序列中存在的异常值进行探测，探测结果如图 7所示。根据图 7(a)中的频率数据可以看出，在钟差序列中存在些许较小的粗差，而图 7(b)中的频率数据没有发现显著的跳变，这表明本文所改进的Baarda粗差探测法能够有效探测出钟差序列中的小粗差。尤其是针对频漂特性比较明显的钟差序列，采用传统的Baarda粗差探测法和中位数法均不能得到很好的处理，而采用本文所改进的Baarda粗差探测法效果十分明显。另外，由于钟跳的出现会对钟差性能指标的计算产生较大的影响，有必要对钟差序列中存在的钟跳进行处理，将钟跳当作频率数据中的粗差进行处理，能够较好地处理钟差序列中存在的钟跳情况。

3 钟差特性分析 3.1 周期特性分析

卫星轨道和钟差高度相关，尤其是轨道径向的误差会被卫星钟差所吸收，因此在钟差数据中会有一定的周期特性^[13]。本实验采用§2中处理好的卫星钟差数据，在去除二次趋势项之后，得到相应的钟差残差数据，再用快速傅里叶变换的方法获得卫星钟差的主要频率，以此来分析北斗卫星钟差的周期特性^[14-15]。在采用傅里叶变换将时间域信号转换到频率域信号的过程中，原始的复杂信号可以分解为多个简单的正弦和余弦信号的叠加，在这些简单的正弦和余弦信号中, 分别有其各自的频率, 并且都可以通过功率谱体现出来，功率谱越大, 则表明对应频率的信号在原始信号中的作用越大, 否则越小^[16-19]。因此，可以根据频谱分析结果中振幅的大小，选择钟差数据中的显著周期项。BDS卫星钟差主频率分析结果如图 8、图 9所示，图中横轴给出的是一天中出现某个周期的次数。

从图 8、图 9可以看出:(1)BDS-2和BDS-3的IGSO卫星以及BDS-2的C14卫星有明显的12 h和24 h周期项，另外除了C08，其他卫星还存在明显的8 h或6 h周期项，并且随着周期的缩短，振幅有逐渐降低的趋势，与文献[16]所得的结论一致。在BDS-2和BDS-3的MEO卫星中，C11、C28主要周期项不明显。(2)造成周期性不一致的原因可能是卫星钟差对外界环境变化比较敏感，导致估计的卫星钟差除了跟轨道周期耦合外，还会受到昼夜变化的影响。卫星钟差本身也存在一定的误差，导致拟合残差中也存在误差。另外能够接收到BDS-3信号的测站少、分布不均匀，导致钟差数据存在部分历元的缺失，也可能会造成周期性不一致。(3)不同的卫星主要周期项不同，钟差数据中的周期项将会对预报精度产生较大的影响。在对钟差数据进行建模预报时，必须考虑到周期特性的影响，根据文献[17]所得的研究结果，每颗卫星选择两个主要周期项加到式(1)中，构建附加周期项的钟差模型如下：

$ \begin{array}{*{20}{c}} {x(t) = {a_0} + {a_1}t + {a_2}{t^2} + \sum\limits_{i = 1}^P {\left( {{A_i}\sin \left( {\frac{{2\pi }}{T}t} \right) + } \right.} }\\ {\left. {{B_i}\cos \left( {\frac{{2\pi }}{T}t} \right)} \right) + {\mathit{\Delta} _i}} \end{array} $

(5)

式中，P为周期项的个数; A_i、B_i、T分别为正、余弦函数的振幅和周期；${{\mathit{\Delta} _i}} $表示附加周期项拟合残差。

图 10为C14、C19、C34 3颗卫星不加周期项和加入两个主要周期项之后的拟合残差序列图，可以看出，加入两个周期项后，拟合残差序列更加稳定，拟合精度得到了很大的提升，这也证明了不同卫星钟差序列中存在几何周期项^[17]。

3.2 频率稳定性分析

本实验仍采用§2处理后的数据，将钟差序列中的粗差、钟跳、天跳变以及频率跳变进行细密的处理后，计算卫星钟Allan方差和Hadamard方差(下文简称Had)，结果如图 11、图 12所示。表 3中给出了BDS-2的IGSO、MEO卫星和BDS-3卫星的千秒稳定性、万秒稳定性和日稳定性。可以看出，采用Allan方差和Had方差计算的频率稳定性指标基本一致。表 4中还统计了两种算法计算的频率稳定性指标之间的相似度，结果发现，相似度达到95%以上，且平滑时间在600~2 000 s左右时，Allan方差和Had方差序列均接近$ \tau^{-\frac{1}{2}}$，主要受到频率白噪声的影响^{[13, 16, 18-19]}。超过2 000 s时，序列表现出复杂的、非幂率变化，并在10 000 s附近出现一个“凸起”，这是由于周期特性引起的，可见钟差序列中的周期特性极大地影响了钟差的频率稳定性。

表 5给出了BDS-3和BDS-2卫星的频率稳定性指标的对比结果，统计结果表明，BDS-3稳定性整体优于BDS-2卫星，BDS-3的频率稳定性相对于BDS-2频率稳定性最高提升了75.09%。BDS-3卫星性能的提升主要得益于载荷有新型国产铷钟和被动型氢原子钟。

表 6给出了BDS系统IGSO和MEO卫星频率稳定性的统计结果。从表 6中可以看出，MEO卫星的频率稳定性整体优于IGSO卫星，而BDS-3 IGSO卫星的日稳定性优于MEO卫星的日稳定性，主要由于BDS-3 C32卫星载荷有氢原子钟，侧面反映出氢原子钟的长期稳定性要优于铷原子钟。

表 7给出了BDS系统搭载的不同类型原子钟的频率稳定性的统计结果(实验期间仅C32卫星启用了氢原子钟)，结果表明，氢原子钟千秒稳定性、万秒稳定性和日稳定性分别高于铷原子钟40.16%、47.10%、65.60%，并且氢原子钟的长期稳定性要明显优于铷原子钟。

4 结语

本文详细分析了BDS-2 IGSO、MEO卫星和BDS-3卫星的周期特性和频率稳定性，采用不同算法讨论了BDS系统不同卫星类型之间、不同星载原子钟类型之间的频率稳定性的差异。分析实验结果得到以下结论：(1)采用Allan方差和Had方差估计的BDS卫星钟差频率稳定性基本一致，相似度达到95%以上。(2)BDS-3卫星载荷了新一代的国产铷钟和被动型氢原子钟，其性能相对于BDS-2系统来说具有较大的提升，BDS-3卫星的频率稳定性相对于BDS-2卫星提升了20.85%~75.09%。(3)BDS卫星中MEO卫星的频率稳定性要优于IGSO卫星。(4)氢原子钟千秒稳定性、万秒稳定性和日稳定性分别高于铷原子钟40.16%、47.10%、65.60%，且氢原子钟的长期稳定性要明显优于铷原子钟。

致谢: 感谢江苏省双创团队项目支持和国际GNSS监测评估系统(iGMAS)提供的数据产品。