李變，屈俐俐，陳永奇

毫秒脈沖星計(jì)時(shí)噪聲處理方法研究

李變1,2，屈俐俐1,2，陳永奇1,2

（1. 中國科學(xué)院 國家授時(shí)中心，西安 710600；2. 中國科學(xué)院 時(shí)間頻率基準(zhǔn)重點(diǎn)實(shí)驗(yàn)室，西安 710600）

毫秒脈沖星；計(jì)時(shí)噪聲；經(jīng)驗(yàn)?zāi)B(tài)分解

0 引言

自1982年D. C. Backer等人發(fā)現(xiàn)第一顆毫秒脈沖星PSR B1937+21以來，便開始了毫秒脈沖星計(jì)時(shí)觀測。因毫秒脈沖星自轉(zhuǎn)周期非常穩(wěn)定，被譽(yù)為自然界中最穩(wěn)定的時(shí)鐘。長期的計(jì)時(shí)觀測表明，毫秒脈沖星自轉(zhuǎn)周期的變化率很小，文獻(xiàn)[1]給出毫秒脈沖星PSR B1855+09和PSR B1937+21一年以上的頻率穩(wěn)定度優(yōu)于10-14，這與現(xiàn)有守時(shí)用主流商品原子鐘的穩(wěn)定度水平相當(dāng)。因此，通過對毫秒脈沖星計(jì)時(shí)觀測數(shù)據(jù)進(jìn)行分析處理，利用合適的算法就可以建立僅依賴于自然存在的天體的時(shí)間尺度——脈沖星時(shí)間尺度。該時(shí)間尺度的建立可以用于原子時(shí)或其他時(shí)間尺度長期性能的檢測與驗(yàn)證，對于深空導(dǎo)航、定位及科學(xué)研究具有重大而深遠(yuǎn)的意義，這使得脈沖星計(jì)時(shí)陣觀測項(xiàng)目成為目前多國科學(xué)家研究的熱點(diǎn)。

脈沖星計(jì)時(shí)殘差是脈沖預(yù)報(bào)到達(dá)時(shí)間與實(shí)際到達(dá)時(shí)間之差，如果用于預(yù)報(bào)的脈沖星鐘模型非常準(zhǔn)確，那么計(jì)時(shí)殘差將主要受測量誤差的影響。然而受轉(zhuǎn)換誤差、軌道伴星及自轉(zhuǎn)不穩(wěn)定等不確定因素的影響，計(jì)時(shí)殘差中表現(xiàn)出難以模型化的不規(guī)則特性。計(jì)時(shí)噪聲就是這種不規(guī)則性表現(xiàn)的主要類型之一，它由低頻信號組成，是計(jì)時(shí)殘差中不可預(yù)報(bào)的長期變化趨勢。由于計(jì)時(shí)殘差的統(tǒng)計(jì)量符合非平穩(wěn)過程的定義，即波動(dòng)性和均值具有趨向性，因此對于計(jì)時(shí)殘差預(yù)測的關(guān)鍵是如何提取其時(shí)間序列中不同的時(shí)頻特性，構(gòu)建時(shí)間序列中的趨勢和波動(dòng)模型。計(jì)時(shí)殘差的趨勢由低頻信號組成，是不可預(yù)報(bào)的長期變化趨勢，波動(dòng)則屬于高頻信號，其正態(tài)分布的特性符合時(shí)間序列平穩(wěn)性的要求。

脈沖星鐘模型參數(shù)是利用最小二乘法方法，對計(jì)時(shí)殘差進(jìn)行多次迭代擬合得到，如果是雙星系統(tǒng)，脈沖預(yù)報(bào)到達(dá)SSB處的時(shí)間還應(yīng)考慮雙星系統(tǒng)的影響。

計(jì)時(shí)噪聲產(chǎn)生的原因主要有脈沖星自身的自轉(zhuǎn)不穩(wěn)定、軌道伴星，以及轉(zhuǎn)換誤差、歷表誤差及引力波等[3-4]。計(jì)時(shí)噪聲對脈沖星鐘模型參數(shù)的精確測定有重要影響，具體包括：計(jì)時(shí)噪聲的線性項(xiàng)將影響脈沖星自轉(zhuǎn)頻率參數(shù)，二次項(xiàng)變化會被脈沖星自轉(zhuǎn)頻率一階導(dǎo)數(shù)吸收，高階項(xiàng)會影響脈沖星自轉(zhuǎn)頻率高階導(dǎo)數(shù)參數(shù)，某種周期性變化也可能會影響天體測量參數(shù)的正確擬合。對脈沖星時(shí)間尺度建立來說，計(jì)時(shí)噪聲不僅降低計(jì)時(shí)殘差的精度，而且影響脈沖星時(shí)間尺度的準(zhǔn)確度和穩(wěn)定度。

經(jīng)驗(yàn)?zāi)B(tài)分解（EMD）方法就是從復(fù)雜信號分離出本征模態(tài)分量（IMF）的篩選過程，從信號處理角度來講，EMD是一個(gè)不斷從高頻到低頻的濾波過程，體現(xiàn)了多分辨分析的特性[10-12]。對于均值具有趨向性的非線性、非平穩(wěn)信號的趨勢項(xiàng)的提取，EMD方法是非常適合的。脈沖星計(jì)時(shí)噪聲是均值具有趨向性的非平穩(wěn)信號，本文將EMD方法用于毫秒脈沖星計(jì)時(shí)噪聲的提取和白化處理。

1 計(jì)時(shí)噪聲的穩(wěn)定度評估

1.2 計(jì)時(shí)噪聲表現(xiàn)形式及分析

以毫秒脈沖星PSR J0437-4715和PSR J1939+2134為例進(jìn)行實(shí)驗(yàn)分析，數(shù)據(jù)采用澳大利亞Parkes計(jì)時(shí)陣（PTA）4年多的計(jì)時(shí)觀測資料。計(jì)時(shí)殘差由Tempo2軟件得到，與早期軟件相比，該軟件具有更完善、更精確的計(jì)時(shí)模型，并且增加了大氣延遲改正、由行星引起的Shapiro延遲改正、計(jì)時(shí)噪聲白化以及同時(shí)擬合多顆脈沖星計(jì)時(shí)殘差等功能。圖1和圖2是PSR J0437-4715和PSR J1939+2134的計(jì)時(shí)觀測數(shù)據(jù)經(jīng)過Tempo2軟件處理后得到計(jì)時(shí)殘差。為了分析計(jì)時(shí)殘差的分布情況，給出了這兩顆脈沖星的計(jì)時(shí)殘差柱形分布圖，具體如圖3和圖4所示。

圖1 PSRJ0437-4715的計(jì)時(shí)殘差（RMS=0.109μs）

圖2 PSR J1939+2134的計(jì)時(shí)殘差（RMS=1.399μs）

圖3 PSR J0437-4715的計(jì)時(shí)殘差分布圖

圖4 PSR J1939+2134的計(jì)時(shí)殘差分布圖

圖5 PSR J0437-4715的穩(wěn)定度

圖6 PSR J1939+2134穩(wěn)定度

2 計(jì)時(shí)噪聲處理方法

觀測數(shù)據(jù)不均勻、數(shù)據(jù)量少是毫秒脈沖星計(jì)時(shí)觀測的特點(diǎn)，而且一部分毫秒脈沖星中含有顯著的計(jì)時(shí)噪聲。為了提高計(jì)時(shí)殘差精度，進(jìn)而充分利用脈沖星計(jì)時(shí)觀測資源，有必要研究計(jì)時(shí)噪聲的處理方法。針對計(jì)時(shí)噪聲的特點(diǎn)，將EMD方法用于毫秒脈沖星計(jì)時(shí)噪聲的白化處理[15]。

2.1 經(jīng)驗(yàn)?zāi)B(tài)分解方法

EMD是一種能夠用于非線性、非平穩(wěn)信號處理的，新的時(shí)間序列信號分析方法。從本質(zhì)上講，該方法是對信號進(jìn)行平穩(wěn)化處理的過程。復(fù)雜信號通過EMD方法可以被分解成一系列本征模態(tài)函數(shù)（IMF）的單分量信號。

經(jīng)驗(yàn)?zāi)B(tài)分解方法假設(shè)信號由不同的IMF組成，這些IMF可能是線性的，也可能是非線性的。但是每個(gè)IMF分量都必須滿足兩個(gè)條件：① 極值點(diǎn)個(gè)數(shù)和過零點(diǎn)數(shù)相同或最多相差一個(gè)；② 任意時(shí)刻由極大值點(diǎn)定義的上包絡(luò)線和由極小指點(diǎn)定義的下包絡(luò)線的均值為零，也就是說，信號的上下包絡(luò)線關(guān)于時(shí)間軸對稱。

2.2 數(shù)據(jù)分析

圖7 PSR J0437-4715白化后的計(jì)時(shí)殘差（RMS = 0.075 μs）

圖8 PSR J1939+2134白化后的計(jì)時(shí)殘差（RMS = 0.256 μs）

如圖7所示，PSR J0437-4715白化后計(jì)時(shí)殘差精度由0.109 μs提高到0.075 μs，波動(dòng)范圍也略有改善，這是由于PSR J0437-4715的計(jì)時(shí)噪聲較小，主要受測量噪聲影響。圖9是該毫秒脈沖星白化后的計(jì)時(shí)殘差分布圖。圖8中PSR J1939+2134白化后計(jì)時(shí)殘差的波動(dòng)范圍和精度均有明顯改善，波動(dòng)范圍為（-0.8 ~+0.8 μs），精度從1.399 μs提高到0.256 μs，其白化后的分布也基本符合正態(tài)分布，具體如圖10所示。

圖9 PSR J0437-4715白化后的計(jì)時(shí)殘差分布

圖10 PSR J1939+2134白化后的計(jì)時(shí)殘差分布

圖11 PSR J0437-4715計(jì)時(shí)殘差白化后的穩(wěn)定度

圖12 PSR J1939+2134計(jì)時(shí)殘差白化后的穩(wěn)定度

3 結(jié)論

受軌道運(yùn)動(dòng)的擾動(dòng)、星際介質(zhì)色散量的變化、星際閃爍效應(yīng)、太陽系行星歷表的誤差、地球無線電干擾及鐘的誤差等外部因素和脈沖星本身自轉(zhuǎn)不穩(wěn)定性的影響，脈沖星計(jì)時(shí)結(jié)果中會出現(xiàn)一些無法用模型擬合的、非白色噪聲譜的計(jì)時(shí)噪聲，并且對不同的脈沖星，該噪聲具有不同的表現(xiàn)形式。計(jì)時(shí)噪聲由低頻信號組成，是計(jì)時(shí)殘差中不可預(yù)報(bào)的長期變化趨勢。

趨勢項(xiàng)是信號中緩慢變化的、控制信號變化趨勢的低頻成份。EMD方法根據(jù)特征時(shí)間尺度對信號進(jìn)行分解，得到從高頻到低頻的一系列信號分量，適用于均值具有趨向性的非平穩(wěn)信號趨勢項(xiàng)的提取。本文以PSR J0437-4715和PSR J1939+2134兩顆具有代表性的毫秒脈沖星為例，研究了基于EMD的脈沖星計(jì)時(shí)噪聲提取方法，結(jié)果表明：該方法能夠用于脈沖星計(jì)時(shí)噪聲的白化，對計(jì)時(shí)噪聲占主導(dǎo)的毫秒脈沖星，基于EMD方法白化后的計(jì)時(shí)殘差精度有明顯提高，而且具有方便高效、完全自適應(yīng)的特點(diǎn)。

毫秒脈沖星計(jì)時(shí)觀測數(shù)據(jù)不均勻、數(shù)據(jù)量少、計(jì)時(shí)精度相差懸殊，本文提出的計(jì)時(shí)噪聲處理方法首先提高了可利用的毫秒脈沖星的數(shù)量，使可用觀測數(shù)據(jù)增加；其次，大幅提高計(jì)時(shí)噪聲占主導(dǎo)的毫秒脈沖星計(jì)時(shí)殘差精度，并對計(jì)時(shí)噪聲不明顯的毫秒脈沖星的計(jì)時(shí)殘差精度也有一定提高。

[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ì)時(shí)噪聲估計(jì)[J]. 時(shí)間頻率學(xué)報(bào), 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. 希爾伯特黃變換及其應(yīng)用[M]. 張海勇, 韓東, 王芳, 等, 譯. 2版.北京: 國防工業(yè)出版社, 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] 丁永恒, 童明雷, 趙成仕. 影響脈沖星時(shí)穩(wěn)定度的因素分析[J]. 時(shí)間頻率學(xué)報(bào), 2017, 40(4): 260-267.

[15] 高峰, 高玉平, 童明雷, 等. 經(jīng)驗(yàn)?zāi)B(tài)分解方法在脈沖星計(jì)時(shí)殘差分析及預(yù)報(bào)中的應(yīng)用[J]. 天文學(xué)報(bào), 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ì)時(shí)噪聲處理方法研究[J]. 時(shí)間頻率學(xué)報(bào), 2020, 43(4): 310-317.

2020-01-25；

2020-04-28

國家自然科學(xué)基金資助項(xiàng)目（11473029；11873049；11973046；91736207；42030105）