J0613-0200驾驭铯钟的方法

2023-03-12李变屈俐俐高玉平

航空学报 2023年3期

李变，屈俐俐，高玉平，3

1．中国科学院国家授时中心，西安 710600 2．中国科学院时间频率基准重点实验室，西安 710600 3．中国科学院大学，北京 100049

目前，精度最高的时间尺度国际原子时（TAI）、协调世界时（UTC）和地球时的另一个实现TT（BIPMXX）都是由国际权度局（BIPM）基于原子钟建立并保持的。BIPM首先收集分布在全球80余个守时实验室约500台原子钟的数据资料，采用ALGOS算法加权平均算法得到自由原子时（EAL），然后利用10余台基准频标对EAL进行频率校准后获得TAI。UTC是TAI与世界时相互协调的产物，与TAI只差整数秒，截止目前二者相差37 s［1］。

参与EAL计算的原子钟主要有铯钟和氢钟2种类型，分别约占总钟数的59%和40%，TT（BIPMXX）是BIPM每年年初利用所有可用的基准频标和原子钟的数据资料，采用事后处理的方式获得的滞后一年的时间尺度，也是目前全球最准确、最稳定的原子时尺度［2］。然而TT（BIPMXX）、TAI和UTC是滞后的纸面时间，无法满足时间/频率作为参考实时性的要求。因此，各个国家或地区的守时实验室都建立与保持着具有实时、连续、稳定的物理信号UTC（k），作为UTC的物理实现用于实时测量比对，k为实验室名称缩写。

毫秒脉冲星通常被认为是由双星系统演化而来的，双星系统中的正常脉冲星通过吸积伴星质量来增加角动量，将脉冲星自转周期加速到毫秒量级。与年轻正常脉冲星相比，毫秒脉冲星具有辐射流量弱、自转稳定度高、脉冲形状陡、计时观测精度高、且很少发生自转频率的突变（glitch）等特征，有些自转周期变化率可达到10−19～10−21，被誉为自然界最稳定的“时钟”。1982年Backer等［3］发现第一颗毫秒脉冲星PSR B1937+21之后，人们便开始了毫秒脉冲星的计时观测，关于脉冲星计时资料分析及脉冲星时算法的研究也随之兴起［4-7］，2012年澳大利亚Parkes天文台的Hobbs等得到了可以与地球时TT（BIPMXX）相媲美的脉冲星时［8-9］。

受观测设备、观测技术、星际介质传播延迟等因素的影响，目前毫秒脉冲星最好的计时精度仅为100 ns，远低于原子钟测量比对精度。但原子时的长期频率稳定度受原子钟性能、季节性变化、周年效应、以及原子时算法等因素的影响，长期稳定度不如脉冲星时。由于脉冲星时没有独立定义时刻起点，因此脉冲星只能被视为长期稳定度极高的频率源，用于其他时频设备（例如：原子钟）输出频率长期性能的验证或校准。

近年来随着脉冲星观测设备性能提升，观测到的毫秒脉冲星的数量有所增加、观测精度不断提高。尤其是我国FAST的投入运行，将积累大量的、优质的毫秒脉冲星观测资料［10］，同时中国SKA先导项目已顺利开展，这些为改变目前时频体系单一依赖原子钟的格局，进而提高当前原子时体系的可靠性和稳定性提供有力的数据及技术支持。

脉冲星时与原子时具有优势互补的特点，将脉冲星时用于原子时的检测、校准与驾驭，利用二者的优势，可以建立一种能够在保持原子时短期稳定度优势的基础上，又具有脉冲星时长期稳定度的优势综合时间尺度。

1 驾驭铯钟的参考毫秒脉冲星选取

国际脉冲星计时阵（IPTA）2019年发布了第二批脉冲星计时观测数据（IPTA dr2）。这些观测数据分别来自欧洲脉冲星计时阵（EPTA）、北美脉冲星计时阵（NANOGrav）和澳大利亚Parkes脉冲星计时阵（PPTA）。与IPTA发布的第一批数据（IPTA dr1）相比，IPTA dr2中包含65颗毫秒脉冲星，比IPTA dr1多了16颗［11-12］。

1.1 IPTA dr2概况

IPTA dr2中共有65颗毫秒脉冲星，其中EPTA包含42颗，观测时间跨度7～18 a，NANO⁃Grav包含37颗，观测时间跨度0.6～9.2 a，PPTA包含20颗，观测时间跨度5.6～29.4 a。

与IPTA dr1不同，IPTA dr2分别采用脉冲星计时软件TEMPO2和TEMPOnest拟合计时模型得到A、B两个版本的脉冲星计时残差。IPTA dr2具有以下特点：① 对于同一颗脉冲星多个观测带宽的情况，根据脉冲到达时间（TOA）的不确定度计算权重，并将权重最大的数据序列作为计时模型的参考，进而得到其他时间序列相对于该参考的系统差；② 处理方法更全面，计时残差精度更高。例如：对低偏心率脉冲双星，采用ELL1模型拟合拉普拉斯-拉格朗日参数。对低偏心率中-高倾角的脉冲双星，采用DDH模型对轨道周期的3次谐波振幅和连续谐波的振幅比进行拟合等；③ 计时观测的时间跨度相差大。毫秒脉冲星J1939+2134的观测间隔最长，约29.4 a，J0931-1902的观测间隔最短，只有0.6 a；④ 计时残差数据质量差异大，不仅表现在不同的脉冲星之间，而且同一颗脉冲星不同观测带宽之间的数据质量差异也非常明显。

1.2 参考毫秒脉冲星选取

用于铯钟驾驭的脉冲星参考必须满足测量精度高、观测间隔尽可能长、长期稳定度高，并且计时观测数据的缺失尽可能少等条件。为此，首先在IPTA dr2中选取满足观测间隔τ>15 a、TOA测量精度σTOA<1.5 μs条件，并且计时观测数据缺失尽可能少的6颗毫秒脉冲星作为候选参考（图1）。

图1是对初步选取的编号分别为J0437-4715、J1640+2224、J1713+0747、J0613-0200、J1744-1134和J1939+2134的6颗毫秒脉冲星不同观测波段计时观测数据，采用TEMPO2软件处理后的计时残差。图中横坐标是约化儒略日（MJD）表示的时间序列，纵坐标计时残差。

图1 6颗毫秒脉冲星的计时残差Fig. 1 Timing residuals of 6 millisecond pulsars

利用毫秒脉冲星驾驭铯钟的宗旨是在保持原子钟短期稳定度的基础上，发挥脉冲星长期稳定度的优势。因此，稳定度是必须考虑的关键因素。脉冲星数据采样具有不规则性，其稳定度估计采用与Allan方差类似的σz(τ)估计方法［13］：

式中：τ为采样间隔；c为τ间隔上多项式拟合的三次项系数；表示统计平均。

对选取的6颗候选毫秒脉冲星，采用σz(τ)分别计算稳定度，结果如图2所示。

图2 6颗毫秒脉冲星的稳定度Fig. 2 Stabilities of 6 millisecond pulsars

计时噪声也是决定毫秒脉冲星参考选取的一个重要考虑因素。计时噪声由低频信号组成，是计时残差中不可预报的长期变化趋势［14-17］。表1是6颗毫秒脉冲星的观测间隔、TOA测量精度、计时残差精度及是否含有计时噪声等情况。

表1 6颗毫秒脉冲星比较Table 1 Comparison of 6 millisecond pulsars

根据表1、图1及图2中6颗候选毫秒脉冲星的稳定度，选取J0613-0200作为铯钟驾驭的参考，开展毫秒脉冲星驾驭铯钟的实验研究。

2 J0613-0200驾驭铯钟实验

2.1 实验系统

利用中国科学院国家授时中心（NTSC）的一台5071A铯钟（编号：Cs2928）、相位微调仪、计数器和相位频率分配放大器等相关时频设备以及IPTA dr2的计时观测资料，建立单颗毫秒脉冲星驾驭铯钟实验系统（图3）。

图3中：f为10 MHz频率信号；p为1 PPS信号。时间基准为NTSC建立和保持的中国标准时间UTC（NTSC），作为计数器的开门信号。计数器测量Cs2928及其被驾驭后输出的实验系统时间（APTCs）与UTC（NTSC）的相位差，用于实验数据的分析处理。TEMPO2软件用于分析处理IPTA dr2中的J0613-0200计时观测资料，并基于J0613-0200的参数文件，定期模拟用于Cs2928驾驭的计时观测资料。

图3 驾驭实验系统Fig. 3 Experimental system for steering

2.2 J0613-0200模拟数据

IPTA dr2发布的J0613-0200计时观测数据截止到2014年5月18日（MJD=56795），为了开展毫秒脉冲星驾驭铯钟的实验研究，根据J0613-0200的自转参数、天体测量参数及计时残差的统计特性，利用TEMPO2软件中的fake插件，每月定期模拟用于Cs2928驾驭的J0613-0200计时观测资料。图4是利用J0613-0200已有观测资料和模拟数据计算得到的计时残差。图中横坐标是MJD表示的时间序列，纵坐标是J0613-0200的计时残差。图中黑色曲线（Obs）是根据IPTA dr2中J0613-0200（MJD=50931～56795）的计时观测数据，得到的计时残差，灰色曲线（Sim）是根据J0613-0200模拟数据得到的计时残差。

图4 J0613-0200计时残差Fig. 4 Timing residuals of J0613-0200

J0613-0200模拟数据中计时红噪声的功率谱密度P(f)为［18］

式中：P0为振幅；fc为截止频率；α为谱指数。

2.3 实验系统铯钟性能分析

守时原子钟组通常由铯钟和氢钟组成，铯钟拥有优良的长期稳定度（铯钟的闪烁底噪1×10−14～1×10−15，典型值5×10−15），氢钟拥有出色的短期稳定度（典型氢钟的闪烁底噪1×10−15～1×10−16，典型值5×10−16），但因存在频率漂移，其长期稳定度取决于频率漂移的变化与大小。本实验的目的是证明毫秒脉冲星驾驭原子钟的可行性，因此，选择没有频率漂移或频率漂移很小的铯钟作为被驾驭的原子钟。

Cs2928是NTSC一台连续工作近8年的铯钟，该钟的频率稳定性略高于同类铯钟平均水平，并且拥有完整的历史测量资料，噪声特性长期持续稳定。图5是Cs2928的稳定度曲线。由图5可知：Cs2928的噪声主要包括白色调频噪声（WFM）、频率随机游走调频噪声（RWFM）和闪烁调频噪声（FFM）。Cs2928的频率漂移量值很小（<1×10−17），被淹没在随机游走频率噪声中。因此，Cs2928的驾驭实验不考虑频率漂移的影响，只需根据历史测量资料，结合其噪声特点，准确预报其频率变化。

2.4 驾驭算法

实验系统旨在产生和保持既具有Cs2928短期稳定度，同时又具有J0613-0200长期稳定度特性的时间/频率信号。因此，铯钟驾驭实验系统采用双重驾驭的算法。首先，基于Cs2928的历史数据特征及噪声特点，以地方原子时TA（NTSC）为参考，计算用于每天驾驭的日驾驭量（offset1）。由2.3节内容可知，Cs2928的频率漂移量值很小，被随机游走噪声淹没。因此，每天的日驾驭量仅考虑Cs2928的频率估计；其次，每月根据fake插件模拟的J0613-0200计时观测资料，分析计算用于每月驾驭的月驾驭量（offset2）。

对于月驾驭，需要评估J0613-0200和实验系统产生的时间APTCs之间的偏差。以TAI为共同参考，评估J0613-0200与APTCs之间的偏差：

采用卡尔曼滤波方法，估计滤波后的相位和频率［19-20］，即

式中：B0和B1分别为滤波后的相位和频率；H为测量矩阵。

设X为状态向量；Φ为状态转移矩阵；Zk为测量输入数据(J0613-0200)−APTCs；P为状态向量的协方差矩阵；Q为滤波器的过程噪声矩阵；R为测量噪声矩阵，则

式中：

滤波器初始值由前一计算间隔的(J0613-0200)−APTCs计算得到，卡尔曼滤波估计的频率即为用于每月驾驭的月驾驭量offset2：

这样，每天的驾驭量是这2项之和：

3 实验结果

J0613-0200驾驭Cs2928的日驾驭量每天定时送入相位微调仪对Cs2928进行初次驾驭；月驾驭量每月根据J0613-0200的模拟数据计算，对Cs2928进行二次驾驭，输出实验系统时间APTCs。利用时间间隔计数器将APTCs时间信号与NTSC保持的我国标准时间UTC（NTSC）进行测量比对，用于实验系统时间APTCs的评估。

作为国际标准时间的UTC由TAI闰秒后得到，通常BIPM滞后15～45 d发布，采样间隔为5 d。快速协调世界时（UTCr）是BIPM于2012年开始，在40多个守时实验室和国际时间频率咨询委员会（CCTF）的大力支持下，正式启动的一项研究。近年来，UTCr和UTC的一致性很好，最大偏差优于2 ns，每周三BIPM公布最新结果，并且UTCr的采样间隔为1 d。因此，实验系统时间APTCs的评估以BIPM每周公布的快速协调世界时（UTCr）为参考。

铯钟驾驭实验从2021年2月下旬开始，图6为MJD=59262～59323的Cs2928未经驾驭的运行结果相对于UTCr的相位偏差UTCr－Cs2928和Cs2928经驾驭后的实验结果APTCs相对于UTCr的相位偏差UTCr－APTCs。

图6 J0613-0200未经驾驭和驾驭Cs2928实验结果Fig. 6 Experimental results of J0613-0200 Unsteering and steering Cs2928

由于APTCs的相位偏差是任意的，因此为了便于分析，UTCr−APTCs曲线在起点MJD=59262处被人为对齐到零。由图6（b）可知，APTCs与UTCr之间存在1.39×10−15的频偏。主要原因是：①J0613-0200与UTC/TAI之间存在2.82×10−15的频率偏差；②J0613-0200的模拟数据与观测数据存在1.27×10−15的频偏；③Cs2928的一次驾驭采用历史测量数据；④模型估计的不确定性也是一个影响因素。因此，APTCs与UTCr之间存在该量值的频偏是合理的，结果符合毫秒脉冲星驾驭铯钟的预期，证明脉冲星驾驭铯钟实现脉冲星时间尺度频率的可行性。

APTCs与UTCr之间固定的频率偏差，并不影响2个时间尺度的相对稳定性。截止目前，实验系统运行7个多月，APTCs相对于UTCr的频率稳定度曲线如图7所示。