基于双谱滤波的综合脉冲星时算法研究∗

周庆勇 魏子卿 张 华 刘思伟 孙鹏飞 张 奋姜 坤 刘晓刚 李 奎

(1 地理信息工程国家重点实验室西安710054)(2 西安测绘研究所西安710054)(3 信息工程大学地理空间信息学院郑州450051)(4 西安电子科技大学空间科学与技术学院西安710126)(5 中国科学院国家授时中心西安710600)(6 北京通信与跟踪技术研究所北京100090)

1 引言

脉冲星是一类高速自转的中子星, 具有稳定度极高的自转频率, 部分毫秒脉冲星自转周期长期稳定性与原子钟相当甚至更优[1], 这表明脉冲星能够提供一种基于遥远自然天体并持续数百万乃至数十亿年的独立时间频率, 并且长期稳定性优于原子钟. 基于脉冲星高稳定度自转频率可构建一种新的时间尺度, 称为脉冲星时. 脉冲星时具有以下优势: (1)能够对地球时(Terrestrial time, TT)提供独立外部检核, 检测现有原子时的长期稳定性; (2)与原子时的工作原理不同, 脉冲星时是建立在恒星质量天体的运动过程上,不易受干扰; (3)原子钟工作寿命有限，随着原子钟器件的老化，其守时性能逐渐变差，而脉冲星时可稳定持续工作数百万乃至数十亿年, 且作用范围广. 脉冲星时具有高稳定性、全自主性和全宇宙性的特点, 在时空基准建设方面有良好的应用前景. 利用脉冲星计时数据可改进原子时的长期稳定性, 同时结合原子钟的短期稳定性和脉冲星的长期可用性、稳定性, 可构建一个新的综合时间尺度.

国外学者在脉冲星发现不久就意识到脉冲星自转频率在时间基准方面的应用价值. 1982年, Backer等人发现极其稳定的毫秒脉冲星[2]. 1984年II’in等人提出了脉冲时间尺度的概念, 并对脉冲星时间建立原理进行了初步讨论[3]. 1988年, Guinot指出脉冲星计时有利于对地球时的认识[4]. 1991年, Guniot等人分析发现毫秒脉冲星可改善原子时间的长期稳定性义,可用于建立新的时间基准[5]. 1996年Petit等人提出了综合脉冲星时的概念, 研究了综合脉冲星时经典加权方法[6]. 2008年以来, Rodin等人研究了基于维纳滤波的综合脉冲星时加权平均算法, 并利用帕克斯脉冲星计时阵(Parkes Pulsar Timing Array, PPTA)数据进行了分析[7–8]. 2012年, Hobbs等人提出将钟信号包含在计时模型中, 采用全局最小二乘拟合来构建综合脉冲星时, 并利用PPTA的19颗毫秒脉冲星观测数据构建了第1个具有与国际原子时精度相当的脉冲时间尺度TT(PPTA11)[9], 数字11表示采用截至2011年PPTA数据. 2019年, Hobbs等人采用谱估计和贝叶斯估计方法分别对国际脉冲星计时阵(International Pulsar Timing Array,IPTA)观测数据噪声进行建模, 并构建了更加稳定的脉冲星时间尺度TT (IPTA16),发现TT (IPTA16)与TT (BIPM17)具有较好的一致性, 并将TT (IPTA16)的功率谱与TT (BIPM2017)功率谱进行比较, 发现脉冲星时可检查原子时间尺度的不稳定性[10].BIPM是国际计量局(International Bureau of Weights and Measures)的简称, 负责国际原子时(International Atomic Time, TAI)的统一和发布. 近年来, 我国学者跟踪国际上脉冲星时相关技术的发展动态, 在综合脉冲星时算法方面取得较好的进展[11–14]. 2007年仲崇霞[15]将小波分解方法引入综合脉冲星时算法构建中, 提高了综合脉冲星时的稳定性. 2011年, 陈鼎等人介绍了脉冲星时的应用, 并分析了基于PPTA项目4颗毫秒脉冲星的综合脉冲星时[16]. 2016年, 尹东山等人使用Vondrak滤波处理了北美纳米赫兹天文台(North American Nanohertz Observatory for Gravitational Waves, NANOGrav) 9 yr的观测数据, 构建的脉冲星时长期稳定性为3.4×10−15[17–18]. 2017年, 张彩红等人利用PPTA数据构建综合脉冲星时, 比较了脉冲星时与原子时的长期稳定性[19].

脉冲星时最新发展方向之一是改善卫星导航系统时间基准的稳定性. 2018年12月,欧空局宣布运行了一个基于脉冲星的时钟项目(Pulsar Chronos, PulChron), 其目的是通过观测毫秒脉冲星监测和改善伽利略卫星导航系统时间的稳定性[20]. 报道称PulChron已经运行了较长时间, 且其结果令人鼓舞, 但未见更多详细文献, 也可预见将来脉冲星时在我国北斗导航系统中的应用.

单颗毫秒脉冲星的自转频率可定义一个脉冲星时(Pulsar Time, PT), 然而毫秒脉冲星观测及数据处理中受到各种误差的影响, 如原子钟误差、行星历表精度、星际介质不确定性、脉冲星自转不稳性(周期跃变及计时噪声)、设备观测误差等, 这些误差会使脉冲星观测数据中包含高频噪声、高斯白噪声、红噪声等, 往往这些噪声影响是耦合的.除脉冲星自转不稳定性和太阳系行星历表误差外, 其他误差源对不同脉冲星观测的影响是独立的, 可通过多颗毫秒脉冲星定义的综合脉冲星时, 削弱这些噪声源的影响, 提高脉冲星时的稳定性. 此外, 也可使用谱分析方法, 进一步降低脉冲星计时残差中的非高斯噪声影响, 提高综合脉冲星时的长期稳定度. 本文研究了一种基于双谱滤波的综合脉冲星时构建算法, 并利用IPTA 4颗观测时间超过10 yr的毫秒脉冲星数据, 来验证该方法的有效性.

2 综合脉冲星时计时及双谱滤波基本原理

毫秒脉冲星自转极其稳定, 其周期变化率在10−20s/s量级, 1 d周期变化累积量为10−14s. 脉冲星自转具有稳定性、可观测性, 故毫秒脉冲星计时观测可构建基于脉冲星自转频率的时间尺度, 称为脉冲星时. 脉冲星时是一种相对的时间尺度, 主要通过脉冲星时与参考原子时间的钟差测量而实现. 对毫秒脉冲星进行高精度计时观测, 可得到脉冲星的计时残差, 计时残差是太阳系质心处脉冲到达时间预报值和观测值之差. 预报值是基于反映脉冲星自转规律的脉冲相位模型(也称钟模型), 代表了脉冲星时, 而通过射电观测得到脉冲到达时间是以原子时(AT)为参考的, 故计时残差含有脉冲星时与参考原子时之间的钟差. 如果脉冲星计时处理中, 各种效应得到完全修正, 那么计时残差即为PT与AT之差. 脉冲信号到达测站的时刻是通过测站氢钟记录的, 之后通过钟差修正溯源到TAI. 当前, BIPM提供了TAI的两种实现方式—准实时TT (TAI)和事后TT(BIPM).

随着大量高精度原子钟的加入、新型基准钟的涌现以及原子时算法的改进, TAI的准确度和稳定度不断提升.然而考虑实时性的需求,BIPM每月计算并发表一个版本TAI,称为TT (TAI). 由于兼顾修正模型优化和计算效率, TT (TAI)精度会受到一定程度上的限制, 故其不是一个最高精度的时间尺度. 同时, BIPM提供了另一种更高精度的时间尺度, 即TT (BIPM). TT (BIPM)每年计算一个版本, 利用了全球所有基准钟数据计算重新得到TT的精确实现, 并在计算时段内对所采用的基准钟加入了黑体辐射、相对论效应等更多改正, 且对频率进行平滑和插值后得到, 使得TT (BIPM)比TT (TAI)更准确和更稳定. 本文选择了TT (BIPM2015)作为脉冲星观测数据处理的时间基准.

综合脉冲星时定义为多个脉冲星计时残差的加权平均:

式中PTi为第i个脉冲星时,ωi为权重, 以脉冲星观测误差σi作为计算依据.

式中n为脉冲星的数量. 对于脉冲星时间稳定度, Matsakis等人提出了σz估计方法进行评估[1]. 假设脉冲星的天体测量参数、脉冲相位模型参数及双星模型参数都得到精确测定, 暂不考虑星际介质延迟、计时噪声、行星历表、参考原子时等误差, 那么扣除二次多项式后的残差就是脉冲星时与参考原子时剩余的最低阶偏差, 因此σz适用于脉冲星时间稳定度的估计[1].σz估计方法计算基本思路如下, 其将时间残差测量序列分成m个子序列(m=1,2,4,8,···), 利用每一个子序列的残差测量数据进行最小二乘拟合, 使

式中χ(tj)为观测时间tj的计时残差,ψj为其误差,j为观测的子序列号,m为子序列个数.X(tj)的表达式为:

式中ck(k=0,1,2,3)为其系数,t0为参考历元.σz定义如下:

式中τ为子序列连续采样时间,⟨⟩表示对所有子序列的c3进行加权平均, 权重类似于(2)式的计算.σz估计的上下限采用统计方法估计, 详见参考文献[1].

双谱滤波是通过计算计时残差的双谱, 利用双谱频域滤波器对双谱幅值和相位进行滤波处理, 消除信号噪声, 再对信号进行重构获取滤波后的信号. 非参数直接法双谱估计的方差较小, 方法简单且容易实现, 所以在信号处理中得到了较好的应用, 直接法双谱估计的原理请参考文献[21]. 双谱对高斯过程是“盲”的, 因此通过双谱滤波可以很好地抑制信号中的高斯噪声, 同时在双谱域结合阈值滤波器, 也可有效地抑制加性噪声、乘性噪声等一些非高斯噪声成分.

3 算例分析

射电望远镜都是多科学任务的实体, 脉冲星计时观测无法做到等间隔观测. 从IPTA发布的最新数据来看, 不同脉冲星观测数据起始观测时间不一样, 观测间隔非常不均匀,有些脉冲星间隔1 yr才有观测, 特别是早期数据, 这种现象尤为明显. 随着技术的发展,脉冲星观测精度提高, 计时残差误差减小. 综合脉冲星时的计算及估计流程如下:

(1)计算原始数据的综合脉冲星时. 对所有观测序列去重, 并按照时间排序. 以天为单位, 对位于该天所有观测系列进行加权平均, 权重计算见(2)式. 如果该天没有观测, 不均匀化处理, 设置为无数据, 避免引入未知插值误差. 得到新的计时残差系列, 即综合脉冲星时与参考原子时之差;

(2)计算脉冲星时的σz. 按照(3)–(5)式计算综合脉冲星时的稳定度, 同时考虑两个约束条件: 每个计时残差子序列中观测量数大于4; 子序列中时间间隔要大于子序列时间长度的同时对于m=1的情况, 采用最新观测序列的计时残差;

(3)使用双谱滤波对每个脉冲星的观测序列进行处理, 消除高斯噪声的影响, 然后重复(1)–(2)步, 得到滤波后综合脉冲星时的稳定度.

2019年IPTA发布了第2批共65颗脉冲星的观测数据[22], 包括两种类型的数据.第1种类型数据继承了第1批IPTA数据的形式, 拟合了色散量(Dispersion Measure, DM)和白噪声参数化. 第2种数据在第1类基础上进行红噪声拟合. 本文选择了比较稳定的PSR J0437−4715、J0613−0200、J1713+0747和J1909−37444颗脉冲星的第1类数据, 没有处理红噪声. 4颗脉冲星的基本情况见表1, 拟合后残差见图1, 计时残差的时间参考为TT (BIPM2015), 历表采用DE436行星历表. 表1中PSR为脉冲星名称, Period为毫秒脉冲星的自转频率, RMS为计时残差(residual, Res)的均方根值, Span为观测周期长度, NTOA为参与解算的观测量数量, Start和Finish为观测数据的起始、结束时刻, 用简化儒略日(Modified Julian Date, MJD). 由图1可知, 由于观测设备更新和时间系统提升, 脉冲星观测精度也相应地提高.

表1 IPTA 4颗脉冲星的基本信息Table 1 The basic information of four IPTA pulsars

图1 IPTA 4颗脉冲星的计时残差(时间基准TT (BIPM2015))Fig.1 Timing residuals of four IPTA Pulsars (time reference is TT (BIPM2015))

脉冲星观测及数据处理中, 不可避免地受到各种误差源的影响, 如各种残余的微小射电干扰、氢钟的不稳定性及设备的热噪声、数据处理中DM偏差、行星历表原点误差、计时噪声难以白化处理等. 这些影响会在计时残差中产生高斯及非高斯噪声, 而双谱滤波可抑制其中部分噪声. 每颗脉冲星的计时残差包含的噪声成分不尽相同, 因此表现出不同残差形态, 双谱滤波方法可能对不同脉冲星的处理效果不尽相同. 作为对比,按照上节计算过程, 采用原始和滤波后残差数据分别构建4颗脉冲星综合时, 见图2.

脉冲星观测具有独立性, 不具有时间关联性, 故在数据处理中不建议均匀插值采样,那样做尽管可填补空白时间段的观测数据, 但插值没有真实观测量对应, 容易白化和平滑各种噪声, 导致综合脉冲星时评估精度过高. 由图1和图2可知, 基于原始数据得到综合脉冲星时, 实际上是4颗脉冲星时的加权平均, 综合脉冲星时精度主要受观测精度略差的PSR J0613−0200、J1713+0747的影响. 基于双谱滤波数据的综合脉冲星时精度得到明显的提升, 计时残差的平均值从5.69×10−7s减小到1.12×10−7s, 方差从1.25×10−6s降低到9.58×10−8s, 精度有大幅度提高. 使用σz评估两种方法的综合脉冲星时稳定性, 分别见图3、图4.

图2 IPTA 4颗脉冲星构建的综合脉冲星时(时间基准参考于TT (BIPM2015), 上图基于原始数据, 下图基于滤波后数据).Fig.2 Ensemble pulsar time of four IPTA Pulsars (time reference is TT (BIPM2015), the top panel is based on the original data, and the bottom panel is based on the filtered data).

图3 基于IPTA 4颗脉冲星原始数据构建的综合脉冲星时的稳定度Fig.3 The stability of ensemble pulsar time established by the original data of four IPTA Pulsars

由图3可知, 4颗脉冲星J0437−4715、J0613−0200、J1713+0747、J1909−3744和综合脉冲星时的月稳定度分别为2.26×10−13、2.37×10−11、1.03×10−9、2.55×10−13、3.77×10−13; 年稳定度分别为3.30×10−14、1.94×10−13、7.32×10−14、1.65×10−14、7.77×10−14; 10 yr的稳定度为1.23×10−15、2.25×10−14、2.86×10−15、8.75×10−16、8.56×10−16. 基于所有滤波数据实现的综合脉冲星时稳定度为2.03×10−16, 可见4颗脉冲星时与综合脉冲星时的稳定度随着观测时间增加而更加稳定, 稳定度较好的两颗脉冲星时是观测精度高的J0437−4715和J1909−3744. 当观测时间超过10 yr,综合脉冲星时的稳定性更优, 原因在于观测数据充分, 且有大批高质量的新观测数据加入. 由图4可知, 基于双谱滤波算法构建4颗脉冲星时和综合脉冲星时, 其稳定度的基本规律与基于原始数据构建的脉冲星时一样, 但稳定度更高. 4颗脉冲星J0437−4715、J0613−0200、J1713+0747、J1909−3744和综合脉冲星时的月稳定度分别为3.77×10−14、1.60×10−12、6.75×10−11、4.76×10−14、7.40×10−14; 年稳定度分别为6.94×10−15、4.39×10−14、2.03×10−14、3.19×10−15、1.50×10−14; 10 yr的稳定度为2.91×10−15、1.74×10−14、6.53×10−16、7.36×10−17、3.50×10−16. 基于所有滤波数据实现的综合脉冲星时稳定度为8.21×10−17. 相比于经典加权算法, 基于双谱滤波算法的综合脉冲星时构建算法, 能够较好地提升单脉冲星时及综合脉冲星时的稳定性. 需要说明的是, 当时间尺度大于数据观测周期的一半时, 即m=1, 会优先使用最新的观测数据, 这样的数据处理策略会导致脉冲星时的稳定度较好.

图4 基于IPTA 4颗脉冲星滤波数据构建的综合脉冲星时的稳定度Fig.4 The stability of the ensemble pulsar time established by the filtering data of four IPTA Pulsars

BIPM定期发布各实验室保持的地方原子时TAI (XX)与国际原子时TAI的偏差文件, XX代表具体实验室, 用英文名称简写表示. 同样使用σz方法估算了中科院国家授时中心(National Time Service Center, NTSC)、德国技术物理研究所(Physikalisch-Technische Bundesanstalt, PTB)、俄罗斯国家技术物理及无线电工程研究院(VNIIFTRI, 简写SU)、美国国家标准技术研究院(National Institute of Standards and Technology, NIST) 4家授时单位原子时的稳定性, 时间从MJD50814—MJD58904. PTB是目前国际上所有授时单位中原子钟精度最高的单位之一, 在国际综合原子时建立中的权重较大, TAI (PTB)的数据由于基准钟升级存在明显的2次跳变, 且第2段数据较短, 故本文分别处理了第1段和第3段数据, 第3段数据精度明显优于第1段, 下面分别标记为TAI (PTB1)、TAI (PTB2), 第1段数据周期为MJD50814—MJD55439, 第2段数据周期为MJD56079—MJD58904. 4家单位原子钟稳定性见图5.

原子时的稳定度一般使用Allan方差表示, 由于方便与综合脉冲星(Ensemble Pulsar Time, EPT)的比较, 也使用σz方法评估其稳定度. 由图5可知,σz方法表示的原子时稳定度基本上保持稳定, 首先随着观测时间原子时稳定度缓慢提高, 约6 yr后稳定性下降.TAI (NIST)、TAI (SU)、TAI (NTSC)、TAI (PTB1)和TAI (PTB2)的年稳定度分别为1.19×10−12、1.38×10−13、3.53×10−14、2.31×10−13、1.46×10−15; 5 yr稳定度分别为1.05×10−13、1.12×10−14、2.51×10−15、5.51×10−14、2.22×10−16. 综合脉冲星时的年稳定度优于TA(NIST)、TA(SU)、TAI(NTSC)、TAI(PTB1), 逊色于TAI(PTB2),TAI (PTB2)代表当前原子钟的最高水平. 下面以TAI (PTB2)为例进行分析, 随着观测精度的提高, 综合脉冲星时稳定性提升, 而原子时稳定性恒定, 5 yr时间两者稳定性相当,然而5 yr后综合脉冲星时优于TAI (PTB2), 可用于改善原子时稳定性. 脉冲星时短稳主要受频率闪烁噪声和随机游走频率噪声等系统噪声影响, 主要来源包括脉冲星自转的不规则性、观测系统噪声或是星际闪烁噪声. 脉冲星时长稳主要受脉冲星自转不规律性的影响, 主要包括计时噪声, 故精确的计时噪声模型、优秀的观测设备及更合理的观测方案有利于脉冲星时稳定度的提升.

图5 国际上4家授时单位原子时的稳定度Fig.5 Stability of atomic time of four international time service units

4 结论与分析

本文研究了一种基于双谱滤波的综合脉冲星时算法, 并利用IPTA的4颗毫秒脉冲星观测数据展示其良好性能, 比较分析了脉冲星时与原子时的稳定性. 由于脉冲星的动能会因自身的辐射和周围吸积的旋转物质同磁层间的相互作用这两个因素逐渐消耗, 使脉冲星旋转角动量减少, 所以脉冲星自转周期都是随着时间逐渐变长, 然而也存在一类自转不稳定性现象, 如计时噪声. 计时噪声是指脉冲星自转参数发生的连续、时标较长(通常为几个月或几年)的扰动. 一般来说, 计时噪声表现为无规律的随机信号, 是一种均值为零的高斯分布, 该类噪声可使用双谱滤波较好地消除. 然而不少脉冲星的计时噪声表现为准周期性、红噪声, 其功率在低频端更强. 当前计时噪声产生原因不明确, 普遍认为其与中子星内部的超流过程、内部温度的变化以及磁层中的物理过程有关. 下一步可对红噪声进行建模拟合, 如采用Cholesky方法和贝叶斯估计方法对计时噪声进行估计,减弱红噪声的影响.

本文的研究也可为我国脉冲星时间系统的建设及应用提供参考. 我国新疆天文台南山射电望远镜于1996年率先在国内开展脉冲星观测, 对脉冲星已经连续观测20多年, 同时我国于2016年建成的全球最大的射电望远镜500 m口径球面射电望远镜(Fivehundred-meter Aperture Spherical radio Telescope, FAST)能实现当前毫秒脉冲星最高精度的观测, 它们都将为我国脉冲星时间尺度的研究提供宝贵的资料.