毫秒脉冲星计时噪声处理方法研究

李变，屈俐俐，陈永奇

毫秒脉冲星计时噪声处理方法研究

李变1,2，屈俐俐1,2，陈永奇1,2

（1. 中国科学院 国家授时中心，西安 710600；2. 中国科学院 时间频率基准重点实验室，西安 710600）

毫秒脉冲星；计时噪声；经验模态分解

0 引言

自1982年D. C. Backer等人发现第一颗毫秒脉冲星PSR B1937+21以来，便开始了毫秒脉冲星计时观测。因毫秒脉冲星自转周期非常稳定，被誉为自然界中最稳定的时钟。长期的计时观测表明，毫秒脉冲星自转周期的变化率很小，文献[1]给出毫秒脉冲星PSR B1855+09和PSR B1937+21一年以上的频率稳定度优于10-14，这与现有守时用主流商品原子钟的稳定度水平相当。因此，通过对毫秒脉冲星计时观测数据进行分析处理，利用合适的算法就可以建立仅依赖于自然存在的天体的时间尺度——脉冲星时间尺度。该时间尺度的建立可以用于原子时或其他时间尺度长期性能的检测与验证，对于深空导航、定位及科学研究具有重大而深远的意义，这使得脉冲星计时阵观测项目成为目前多国科学家研究的热点。

脉冲星计时残差是脉冲预报到达时间与实际到达时间之差，如果用于预报的脉冲星钟模型非常准确，那么计时残差将主要受测量误差的影响。然而受转换误差、轨道伴星及自转不稳定等不确定因素的影响，计时残差中表现出难以模型化的不规则特性。计时噪声就是这种不规则性表现的主要类型之一，它由低频信号组成，是计时残差中不可预报的长期变化趋势。由于计时残差的统计量符合非平稳过程的定义，即波动性和均值具有趋向性，因此对于计时残差预测的关键是如何提取其时间序列中不同的时频特性，构建时间序列中的趋势和波动模型。计时残差的趋势由低频信号组成，是不可预报的长期变化趋势，波动则属于高频信号，其正态分布的特性符合时间序列平稳性的要求。

脉冲星钟模型参数是利用最小二乘法方法，对计时残差进行多次迭代拟合得到，如果是双星系统，脉冲预报到达SSB处的时间还应考虑双星系统的影响。

计时噪声产生的原因主要有脉冲星自身的自转不稳定、轨道伴星，以及转换误差、历表误差及引力波等[3-4]。计时噪声对脉冲星钟模型参数的精确测定有重要影响，具体包括：计时噪声的线性项将影响脉冲星自转频率参数，二次项变化会被脉冲星自转频率一阶导数吸收，高阶项会影响脉冲星自转频率高阶导数参数，某种周期性变化也可能会影响天体测量参数的正确拟合。对脉冲星时间尺度建立来说，计时噪声不仅降低计时残差的精度，而且影响脉冲星时间尺度的准确度和稳定度。

经验模态分解（EMD）方法就是从复杂信号分离出本征模态分量（IMF）的筛选过程，从信号处理角度来讲，EMD是一个不断从高频到低频的滤波过程，体现了多分辨分析的特性[10-12]。对于均值具有趋向性的非线性、非平稳信号的趋势项的提取，EMD方法是非常适合的。脉冲星计时噪声是均值具有趋向性的非平稳信号，本文将EMD方法用于毫秒脉冲星计时噪声的提取和白化处理。

1 计时噪声的稳定度评估

1.2 计时噪声表现形式及分析

以毫秒脉冲星PSR J0437-4715和PSR J1939+2134为例进行实验分析，数据采用澳大利亚Parkes计时阵（PTA）4年多的计时观测资料。计时残差由Tempo2软件得到，与早期软件相比，该软件具有更完善、更精确的计时模型，并且增加了大气延迟改正、由行星引起的Shapiro延迟改正、计时噪声白化以及同时拟合多颗脉冲星计时残差等功能。图1和图2是PSR J0437-4715和PSR J1939+2134的计时观测数据经过Tempo2软件处理后得到计时残差。为了分析计时残差的分布情况，给出了这两颗脉冲星的计时残差柱形分布图，具体如图3和图4所示。

图1 PSRJ0437-4715的计时残差（RMS=0.109µs）

图2 PSR J1939+2134的计时残差（RMS=1.399µs）

图3 PSR J0437-4715的计时残差分布图

图4 PSR J1939+2134的计时残差分布图

图5 PSR J0437-4715的稳定度

图6 PSR J1939+2134稳定度

2 计时噪声处理方法

观测数据不均匀、数据量少是毫秒脉冲星计时观测的特点，而且一部分毫秒脉冲星中含有显著的计时噪声。为了提高计时残差精度，进而充分利用脉冲星计时观测资源，有必要研究计时噪声的处理方法。针对计时噪声的特点，将EMD方法用于毫秒脉冲星计时噪声的白化处理[15]。

2.1 经验模态分解方法

EMD是一种能够用于非线性、非平稳信号处理的，新的时间序列信号分析方法。从本质上讲，该方法是对信号进行平稳化处理的过程。复杂信号通过EMD方法可以被分解成一系列本征模态函数（IMF）的单分量信号。

经验模态分解方法假设信号由不同的IMF组成，这些IMF可能是线性的，也可能是非线性的。但是每个IMF分量都必须满足两个条件：① 极值点个数和过零点数相同或最多相差一个；② 任意时刻由极大值点定义的上包络线和由极小指点定义的下包络线的均值为零，也就是说，信号的上下包络线关于时间轴对称。

2.2 数据分析

图7 PSR J0437-4715白化后的计时残差（RMS = 0.075 µs）

图8 PSR J1939+2134白化后的计时残差（RMS = 0.256 µs）

如图7所示，PSR J0437-4715白化后计时残差精度由0.109 µs提高到0.075 µs，波动范围也略有改善，这是由于PSR J0437-4715的计时噪声较小，主要受测量噪声影响。图9是该毫秒脉冲星白化后的计时残差分布图。图8中PSR J1939+2134白化后计时残差的波动范围和精度均有明显改善，波动范围为（-0.8 ~+0.8 µs），精度从1.399 µs提高到0.256 µs，其白化后的分布也基本符合正态分布，具体如图10所示。

图9 PSR J0437-4715白化后的计时残差分布

图10 PSR J1939+2134白化后的计时残差分布

图11 PSR J0437-4715计时残差白化后的稳定度

图12 PSR J1939+2134计时残差白化后的稳定度

3 结论

受轨道运动的扰动、星际介质色散量的变化、星际闪烁效应、太阳系行星历表的误差、地球无线电干扰及钟的误差等外部因素和脉冲星本身自转不稳定性的影响，脉冲星计时结果中会出现一些无法用模型拟合的、非白色噪声谱的计时噪声，并且对不同的脉冲星，该噪声具有不同的表现形式。计时噪声由低频信号组成，是计时残差中不可预报的长期变化趋势。

趋势项是信号中缓慢变化的、控制信号变化趋势的低频成份。EMD方法根据特征时间尺度对信号进行分解，得到从高频到低频的一系列信号分量，适用于均值具有趋向性的非平稳信号趋势项的提取。本文以PSR J0437-4715和PSR J1939+2134两颗具有代表性的毫秒脉冲星为例，研究了基于EMD的脉冲星计时噪声提取方法，结果表明：该方法能够用于脉冲星计时噪声的白化，对计时噪声占主导的毫秒脉冲星，基于EMD方法白化后的计时残差精度有明显提高，而且具有方便高效、完全自适应的特点。

毫秒脉冲星计时观测数据不均匀、数据量少、计时精度相差悬殊，本文提出的计时噪声处理方法首先提高了可利用的毫秒脉冲星的数量，使可用观测数据增加；其次，大幅提高计时噪声占主导的毫秒脉冲星计时残差精度，并对计时噪声不明显的毫秒脉冲星的计时残差精度也有一定提高。

[1] KASPI V M, TAYLOR J H, RYBA M F. High-precision timing of millisecond pulsars III: long term monitoring of PSRs B1855+09 and B1937+21[J]. The Astrophysical Journal, 1994(428): 713-728.

[2] BACKER D C, HELLINGS R W. Pulsar timing and general relativity[J]. Annual Review of Astronomy and Astrophysics, 1986(24): 537-575.

[3] LENTATI L, SHANNON R M, COLES W A, et al. From spin noise to systematics: stochastic processes in the first international pulsar timing array data release[J]. Monthly Notices of the Royal Astronomical Society, 2016, 458(2): 2161-2187.

[4] HOBBS G, LYNE A, KRAMER M. Pulsar timing noise[J]. Chinese Journal of Astronomy and Astrophysics, 2006, 6(52): 169-175.

[5] URAMA J, LINK B, WEISBERG J M. A strong correlation in radio pulsar with implications for torque variations[J]. Monthly Notices of the Royal Astronomical Society, 2006, 370(1): 76-79.

[6] ARZOUMANIAN Z, NICE D J, TAYLOR J H, et al. Timing behavior of 96 radio pulsars[J]. The Astrophysical Journal, 1994, 422(2): 671-680.

[7] 杨廷高, 童明雷, 高玉平. 毫秒脉冲星计时噪声估计[J]. 时间频率学报, 2014, 37(2): 80-88.

[8] SHANNON R M, CORDES J M. Assessing the role of spin noise in the precision timing of millisecond pulsars[J]. The Astrophysical Journal, 2010, 725(2): 1607-1619.

[9] REARDON D J, HOBBS G, COLES W, et al. Timing analysis for 20 millisecond pulsars in the parkes pulsar timing array[J]. Monthly Notices of the Royal Astronomical Society, 2016, 455(2): 1751-1769.

[10] HUANG N E, SHEN Z, LONG S R, et al. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[J]. Mathematical, Physical and Engineering Sciences, 1998, 454(1971): 903-995.

[11] HUANG N E, SHIN H, LONG S R. The ages of large amplitude coastal seiches on the caribbean coast of Puerto Rico[J]. Journal of Physical Oceanography, 2000, 30(8): 2001-2012.

[12] HUANG N E, SHEN S S P. 希尔伯特黄变换及其应用[M]. 张海勇, 韩东, 王芳, 等, 译. 2版.北京: 国防工业出版社, 2017.

[13] MATSAKIS D N, TAYLOR J H, EUBANKS T M. A statistic for describing pulsar and clock stabilities[J]. Astronomy and Astrophysics, 1997(326): 924-928.

[14] 丁永恒, 童明雷, 赵成仕. 影响脉冲星时稳定度的因素分析[J]. 时间频率学报, 2017, 40(4): 260-267.

[15] 高峰, 高玉平, 童明雷, 等. 经验模态分解方法在脉冲星计时残差分析及预报中的应用[J]. 天文学报, 2018, 59(4): 3-13.

Study on the timing noise analysis method of millisecond pulsar

LI Bian1,2, QU Li-li1,2, CHEN Yong-qi1,2

(1. National Time Service Center, Chinese Academy of Sciences, Xi’an 710600, China; 2. Key Laboratory of Time and Frequency Primary Standards, Chinese Academy of Sciences, Xi’an 710600, China)

millisecond pulsar; timing noise; empirical mode decomposition (EMD)

10.13875/j.issn.1674-0637.2020-04-0310-08

李变，屈俐俐，陈永奇. 毫秒脉冲星计时噪声处理方法研究[J]. 时间频率学报, 2020, 43(4): 310-317.

2020-01-25；

2020-04-28

国家自然科学基金资助项目（11473029；11873049；11973046；91736207；42030105）