尹東山高玉平趙書紅

(1中國科學(xué)院國家授時中心西安710600)

(2中國科學(xué)院時間頻率基準(zhǔn)重點實驗室西安710600)

(3中國科學(xué)院大學(xué)北京100049)

綜合脈沖星時間尺度?

尹東山1,2,3?高玉平1,2趙書紅1,2

(1中國科學(xué)院國家授時中心西安710600)

(2中國科學(xué)院時間頻率基準(zhǔn)重點實驗室西安710600)

(3中國科學(xué)院大學(xué)北京100049)

由于截然不同的物理機(jī)制,毫秒脈沖星可產(chǎn)生與原子時尺度完全不同的另一類時間尺度.利用長期的脈沖星計時觀測數(shù)據(jù)可建立綜合脈沖星時間尺度,但通常脈沖星計時觀測是非等間隔采樣,并且有些脈沖星的觀測間隔過長等,這些因素均對綜合脈沖星時間尺度的產(chǎn)生造成很大的困難.因此,提出了一種新的算法來計算綜合脈沖星時間尺度,即:首先選用三次樣條插值方法進(jìn)行內(nèi)插,并使觀測數(shù)據(jù)間隔均勻;其次采用Vondrak濾波方法對觀測序列進(jìn)行平滑處理并去除高頻噪聲;最后再利用加權(quán)算法產(chǎn)生綜合脈沖星時間尺度.利用NANOGRAV(North American Nanohertz Observatory for Gravitational Waves)最新發(fā)布的長達(dá)9 yr的觀測數(shù)據(jù),計算得到了綜合脈沖星時間尺度,它的長期穩(wěn)定度(取樣間隔大于1 yr)優(yōu)于3.4×10?15.研究結(jié)果表明:這種方法可有效降低脈沖星計時殘差中的噪聲影響,同時可提高綜合脈沖星時間尺度的長期穩(wěn)定度.

時間,脈沖星,方法:觀測,數(shù)據(jù)分析

1 引言

目前,國際原子時TAI(International Atomic Time)是準(zhǔn)實時計算的時間尺度,隨著大量的高精度原子鐘和基準(zhǔn)鐘的涌現(xiàn),以及原子時算法的改進(jìn),TAI的精度得到了不斷提高.但由于考慮到實時性的影響,TAI必須每月進(jìn)行一次計算,且在計算中受到很多約束條件限制,精度會在一定程度上受到影響,因此TAI不是TT(Terrestrial Time)的最優(yōu)實現(xiàn).為了提供更完美的原子時時間尺度,BIPM(International Bureau of Weights and Measures)提供了另一種事后計算的時間尺度,即TT(BIPM).TT(BIPM)每年計算一個新的版本,利用了所有可用的基準(zhǔn)鐘的數(shù)據(jù)計算得到TT的事后實現(xiàn),在計算時間區(qū)間內(nèi)對所采用的基準(zhǔn)鐘加入了黑體輻射以及其他改正,并對頻率進(jìn)行平滑和插值后得到TT(BIPM),這使得TT(BIPM)比TAI更準(zhǔn)確也更穩(wěn)定,因此TT(BIPM)可用作綜合脈沖星時間尺度的參考時間.

由于與原子鐘截然不同的物理機(jī)制,毫秒脈沖星可產(chǎn)生與原子時尺度完全不同的另一類時間尺度,也是理想時間尺度的一種實現(xiàn)方式,同樣具有穩(wěn)定、準(zhǔn)確、可測量的特性[1?2].在2012年,國際天文學(xué)聯(lián)合會時間委員會已經(jīng)成立了一個脈沖星時間尺度工作組,協(xié)調(diào)該領(lǐng)域的國際合作,即說明了國際上開始重視綜合脈沖星時間尺度研究.但一般情況下,由于脈沖星計時觀測是非等間隔采樣、采樣率遠(yuǎn)遠(yuǎn)低于原子時的比對測量、脈沖星鐘模型與計時噪聲不同于原子鐘的模型等,這些因素均對綜合脈沖星時間尺度的計算造成很大的困難,因此脈沖星時間尺度算法是研究的重點.

本文基于NANOGRAV最新發(fā)布的長達(dá)9 yr的觀測數(shù)據(jù),利用新的算法計算綜合脈沖星時間尺度,分析了脈沖星時間殘差中的噪聲影響,并評定了綜合脈沖星時間尺度的穩(wěn)定度.

2 算法原理

2.1 Vondrak濾波

Vondrak濾波方法不需要了解觀測資料的變化規(guī)律,就能夠?qū)Φ乳g隔或非等間隔的觀測資料進(jìn)行有效平滑.Vondrak濾波實質(zhì)上類似于單邊濾波器,通過選擇一個合適的截斷頻率,以保證盡可能多地保存有用信號,在數(shù)據(jù)處理中可以作為數(shù)字濾波器使用.

設(shè)觀測資料x(ti)(i=1,2,···,N),Vondrak濾波的基本準(zhǔn)則為[3]:

2.2 經(jīng)典加權(quán)算法

設(shè)鐘組中有n臺脈沖星鐘,其讀數(shù)為hi(t),i=1,2,···,n,依據(jù)經(jīng)典加權(quán)算法,計算得到一個時間尺度TA(t),一般形式為:

ωi(t)為鐘clocki的權(quán)重.時間尺度TA(t)的噪聲是鐘組中所有鐘的噪聲加權(quán)和:

為了使綜合鐘的噪聲δS(t)最小,通常以如下(7)式?jīng)Q定權(quán)重: )

利用每臺鐘相對于TA(t)的數(shù)據(jù),計算阿倫方差或標(biāo)準(zhǔn)方差σ2i,根據(jù)(7)式計算得到每臺鐘的權(quán)重.

2.3 σz(τ)估計方法

由于脈沖星的自轉(zhuǎn)相位、頻率、自轉(zhuǎn)變慢速率和相應(yīng)天體測量參數(shù)應(yīng)用計時模型已經(jīng)精確測定,所以計時殘差的三次差是脈沖星計時殘差中剩余的最低階偏差,因此σz(τ)是適用于脈沖星時間穩(wěn)定度分析研究的估計方法[1].σz(τ)估計方法的計算流程如下:

(1)按觀測時間順序進(jìn)行先后排序,將觀測時間、計時殘差以及誤差分別記為ti、xi、σi,i=1,2,···,n,所取數(shù)據(jù)的總時間跨度為T=tn?t1.

(2)求出σz(τ),將這些數(shù)據(jù)按觀測時間順序連續(xù)分成時間間隔為τ的子序列.設(shè)t0為任意參考時刻,在每個子序列里應(yīng)用三次多項式,如下式所示:

式中ci(i=0,1,2,3)為待擬合的系數(shù).對(8)式進(jìn)行最小二乘擬合,使得[(xi?X(ti))/σi]2為最小.定義σz(τ)為:

其中??是指所有子序列的加權(quán)平均,權(quán)重反比于c3的誤差平方值.

3 綜合脈沖星時間尺度的算法

由于脈沖星計時觀測是非等間隔采樣,采樣率也遠(yuǎn)遠(yuǎn)低于原子時的比對測量等,這些因素均對綜合脈沖星時間尺度的產(chǎn)生造成很大的困難,因此脈沖星時間尺度算法是研究的重點.綜合脈沖星時間尺度的新算法流程為:

(1)選用三次樣條插值方法來加密觀測數(shù)據(jù),并使觀測數(shù)據(jù)間隔均勻.

(2)采用Vondrak濾波方法來對觀測序列進(jìn)行平滑處理并去除高頻噪聲.

(3)利用加權(quán)算法產(chǎn)生綜合脈沖星時間尺度,并最終對綜合脈沖星時間尺度進(jìn)行穩(wěn)定度分析.具體的計算框圖如圖1所示.

圖1 綜合脈沖星時間尺度算法框圖Fig.1 The schematic of the algorithm of the ensemble pulsar time scale

4 試驗結(jié)果及分析

4.1 試驗數(shù)據(jù)介紹

我們選用了NANOGRAV最新發(fā)表的9 yr來所觀測的37顆毫秒脈沖星的計時觀測數(shù)據(jù),這些數(shù)據(jù)由GBT(Green Bank Telescope)和ARECIBO望遠(yuǎn)鏡觀測得到,因為PSR B1937+21所展現(xiàn)出的計時噪聲無法建立準(zhǔn)確模型并加以改正,而且其噪聲水平明顯高于白噪聲,所以剔除了這顆脈沖星的數(shù)據(jù).其余36顆脈沖星觀測數(shù)據(jù)的詳細(xì)信息列于圖2和表1中[4],表1列出了名稱、自轉(zhuǎn)周期(Period)、色散延遲量(Dispersion Measurement,DM)、擬合殘差(Post- fit residual)、數(shù)據(jù)跨度(Span)、觀測數(shù)據(jù)個數(shù)(Nobs)、初始?xì)v元(First epoch)與結(jié)束歷元(Last epoch)等基本參數(shù).其中所有脈沖星擬合后的計時殘差均由TEMPO2處理得到,歷表選用DE421太陽系行星質(zhì)心歷表[5?7].

4.2 計算結(jié)果及分析

為了不降低脈沖星時間尺度的穩(wěn)定度,我們首先利用三次樣條插值方法將原始觀測數(shù)據(jù)插值為等間隔的數(shù)據(jù),其次采用Vondrak濾波方法對處理后的等間隔數(shù)據(jù)進(jìn)行濾波,其目的是盡量消除高頻噪聲,最后再利用新的算法計算綜合脈沖星時間尺度,并對其穩(wěn)定度進(jìn)行分析.

4.2.1 三次樣條插值

從圖2和表1中可以看出,NANOGRAV的觀測數(shù)據(jù)一般比較稀疏,并且數(shù)據(jù)點之間的間隔不均勻.為此我們嘗試了多種不同的插值方法來增加數(shù)據(jù)點的個數(shù),使數(shù)據(jù)點的采樣間隔均勻化.在嘗試了線性插值、二次項插值、三次樣條插值等方法之后,結(jié)果顯示三次樣條插值能提供最優(yōu)的數(shù)據(jù)序列,此方法是數(shù)據(jù)處理中常用的方法,具體過程略去.

需要特別注意的是,由于不同觀測數(shù)據(jù)序列的時間跨度不同,并且觀測數(shù)據(jù)個數(shù)差異也比較大,因此我們只對觀測數(shù)據(jù)在其觀測時間跨度內(nèi)進(jìn)行三次樣條插值,目的是保證觀測數(shù)據(jù)的可靠性和準(zhǔn)確性,如果進(jìn)行長時間的外推,勢必影響到脈沖星觀測數(shù)據(jù)的精度,進(jìn)而會降低生成的綜合時間尺度的穩(wěn)定性.

圖2 NANOGRAV 36顆源的擬合殘差,觀測歷元從53216(MJD)至56599(MJD)Fig.2 The post- fit residuals of 36 pulsars included in NANOGRAV,within an observation duration from 53216(MJD)to 56599(MJD)

4.2.2 Vondrak濾波降噪分析

由于脈沖星觀測受到宇宙背景輻射、射電頻段的人為干擾、星際介質(zhì)、設(shè)備的熱噪聲等因素影響,會導(dǎo)致脈沖星觀測數(shù)據(jù)中包含高頻噪聲,因此在計算綜合脈沖星時間尺度之前,對脈沖星觀測數(shù)據(jù)進(jìn)行降噪處理.本文選用Vondrak濾波方法來降低觀測數(shù)據(jù)中的高頻噪聲,以脈沖星J1012+5307為例,圖3給出了濾波前后的計時殘差曲線圖. Vondrak平滑因子選為1×10?9,從圖3、圖4和圖5中可以看出,濾波后的曲線能夠反映出脈沖星數(shù)據(jù)的細(xì)節(jié),實測的觀測數(shù)據(jù)相對于濾波后曲線值的殘差在零值附近上下波動,且殘差的分布服從高斯分布,符合相位白噪聲即高頻噪聲的性質(zhì).

表1 文中用到的36顆毫秒脈沖星的詳細(xì)參數(shù)Table 1 The detailed parameters of 36 millisecond pulsars included in this paper

圖3 Vondrak濾波前(點線)后(點虛線)PSR J1012+5307的計時殘差Fig.3 The post- fit residuals of PSR J1012+5307 before(dotted line)and after(dot-dashed line)the Vondrak filter

圖4 原始觀測數(shù)據(jù)與濾波曲線的相位差Fig.4 The phase di ff erence between the raw data and the filtered data

圖5 原始數(shù)據(jù)與平滑曲線值的殘差分布圖Fig.5 The residual distribution between the raw data and the filtered data

4.2.3 綜合脈沖星時間尺度

對36顆毫秒脈沖星的觀測數(shù)據(jù)采用了三次樣條插值以及Vondrak濾波算法后,采用加權(quán)算法對處理后的觀測數(shù)據(jù)計算,最終獲得綜合脈沖星時間尺度.由于各脈沖星的觀測時間跨度是不同的,因此在相同的觀測時間上毫秒脈沖星的個數(shù)是不同的,這也是其與慣用的加權(quán)算法的不同之處.最終計算的脈沖星時間尺度如圖6所示[8].

圖6 綜合脈沖星時間尺度Fig.6 The ensemble pulsar time scale

圖6 所示為綜合脈沖星時間尺度,參考時間為TDB(Barycentric Dynamic Time),均方根為0.056μs.在最初的4～5 yr內(nèi)(2009年之前),數(shù)據(jù)的波動幅度較大,而在2009年之后數(shù)據(jù)的波動幅度便明顯減小,這主要由兩個原因引起:(1)觀測設(shè)備更新,主要是數(shù)據(jù)處理終端的更新,用GUPPI和PUPPI分別取代了GASP和ASP;(2)后期加入大批觀測精度更高、計時殘差更小的毫秒脈沖星的觀測數(shù)據(jù).

選用σz(τ)估計方法來分析最終生成的脈沖星時間尺度的穩(wěn)定性,圖7給出了綜合脈沖星時間尺度的σz(τ)結(jié)果[9],當(dāng)采樣間隔小于100 d時,其不穩(wěn)定度大于1×10?14,而當(dāng)采樣間隔大于100 d時,其不穩(wěn)定度小于1×10?14.從結(jié)果可以看出,脈沖星時間尺度有很好的長期頻率穩(wěn)定度,但短期頻率穩(wěn)定度不佳.分析其原因,當(dāng)取樣間隔小于80 d時,影響頻率穩(wěn)定度的主要原因是頻率閃爍噪聲和隨機(jī)游走頻率噪聲,這兩類噪聲應(yīng)該是由系統(tǒng)噪聲引起的,例如脈沖星自轉(zhuǎn)的不規(guī)則性、觀測系統(tǒng)噪聲或星際閃爍噪聲.而當(dāng)采樣間隔大于80 d時,噪聲的主要成分是頻率白噪聲,通常是由觀測噪聲引起的.因此為進(jìn)一步提高脈沖星時間尺度的穩(wěn)定度和準(zhǔn)確度,還需要引進(jìn)更精確的Timing噪聲模型、更好的觀測設(shè)備以及更長的觀測數(shù)據(jù)序列,這將是脈沖星時間尺度后續(xù)發(fā)展的方向.

分析數(shù)據(jù)表明有多顆脈沖星的計時殘差表現(xiàn)出明顯的紅噪聲特征,目前關(guān)于紅噪聲的產(chǎn)生機(jī)制尚不明確,只是在完成脈沖星的各項參數(shù)擬合之后,尚殘留一些低頻的觀測噪聲[10?11].經(jīng)過前期的研究可以推斷,能夠引起紅噪聲的因素有很多,例如,觀測采用的地球質(zhì)心歷表誤差、TT與TCB(Barycentric Coordinate Time)之間轉(zhuǎn)換時所采用的理論表達(dá)式的參數(shù)誤差、脈沖星自身的物理機(jī)制,以及某些常數(shù)采用值誤差等,均是計時殘差中紅噪聲的成因.在后續(xù)的工作中,我們計劃對擬合模型進(jìn)行進(jìn)一步優(yōu)化,逐步剝離紅噪聲中的系統(tǒng)性組成部分,通過對紅噪聲的分析和解析,可以提高各項參數(shù)的擬合精度.

圖7 綜合脈沖星時間尺度的頻率穩(wěn)定度Fig.7 The frequency stability of the ensemble pulsar time scale

5 結(jié)語

由于受到各脈沖星的觀測數(shù)據(jù)采樣時間非等間隔、總觀測時間跨度不同和多種噪聲源等因素的影響,降低了由單顆脈沖星產(chǎn)生的脈沖星時間尺度的穩(wěn)定度和準(zhǔn)確度.為消除這些因素的影響,本文提出了綜合脈沖星時間尺度的新算法,并利用實測數(shù)據(jù)進(jìn)行試驗驗證.試驗結(jié)果表明新的脈沖星時間尺度可有效降低脈沖星時間殘差中的噪聲影響,并提高了綜合脈沖星時間尺度的長期穩(wěn)定度.

致謝感謝美國國立射電天文臺Scott Ranson以及Paul Demorest在觀測數(shù)據(jù)方面所提供的幫助.

[1]仲崇霞.綜合脈沖星時算法及脈沖星時應(yīng)用.西安:中國科學(xué)院國家授時中心,2007

[2]Hobbs G,Coles W,Machester R N,et al.MNRAS,2012,427:2780

[3]Vondrak J.BAICz,1977,28:84

[4]Arzoumanian Z,Brazier A,Burke-spolaor S,et al.ApJ,2015,813:65

[5]Hobbs G,Edwards R T,Manchester R N,et al.MNRAS,2006,369:655

[6]Hobbs G,Jenet F,Lee K J,et al.MNRAS,2009,394:1945

[7]Edwards R T,Hobbs G,Manchester R N.MNRAS,2006,372:1549

[8]Hobbs G,Archibald A,Arzoumanian Z,et al.CQGra,2010,27:2006

[9]Riley W J.Handbook of Frequency Stability Analysis.Washington:U.S.Government printing office, 2008

[10]楊廷高,高玉平,童明雷,等.天文學(xué)報,2015,56:370

[11]Yang T G,Gao Y P,Tong M L,et al.ChA&A,2016,40:256

Ensemble Pulsar Time Scale

YIN Dong-shan1,2,3GAO Yu-ping1,2ZHAO Shu-hong1,2
((1 National Time Service Centre,Chinese Academy of Sciences,Xi’an 710600)
(2 Key Laboratory of Time and Frequency Primary Standards,National Time Service Center,Chinese Academy of Sciences,Xi’an 710600)
(3 University of Chinese Academy of Sciences,Beijing 100049)

Millisecond pulsars can generate another type of time scale that is totally independent of the atomic time scale,because the physical mechanisms of the pulsar time scale and the atomic time scale are quite di ff erent from each other.Usually the pulsar timing observational data are not evenly sampled,and the internals between data points range from several hours to more than half a month.What’s more,these data sets are sparse.And all these make it difficult to generate an ensemble pulsar time scale. Hence,a new algorithm to calculate the ensemble pulsar time scale is proposed.Firstly, we use cubic spline interpolation to densify the data set,and make the intervals between data points even.Then,we employ the Vondrak filter to smooth the data set,and get rid of high-frequency noise, finally adopt the weighted average method to generate the ensemble pulsar time scale.The pulsar timing residuals represent clock di ff erence between the pulsar time and atomic time,and the high precision pulsar timing data mean the clock di ff erence measurement between the pulsar time and atomic time with a high signal to noise ratio,which is fundamental to generate pulsar time.We use the latest released NANOGRAV(North American Nanohertz Observatory for Gravitational Waves)9-year data set to generate the ensemble pulsar time scale.This data set is from the newest NANOGRAV data release,which includes 9-year observational data of 37 millisecond pulsars using the 100-meter Green Bank telescope and 305-meter Arecibo telescope.We find that the algorithm used in this paper can lower the in fl uence caused by noises in timing residuals,and improve long-term stability of pulsar time.Results show that the long-term(>1 yr)frequency stability of the pulsar time is better than 3.4×10?15.

time,pulsars,methods:observational,data analysis

P127

:A

10.15940/j.cnki.0001-5245.2016.03.008

2015-10-26收到原稿,2015-12-21收到修改稿

?國家自然科學(xué)基金項目(11303032,11473029)及地理信息工程國家重點實驗室開放基金課題(SKLGIE2013-M-1-1)資助

?yds@ntsc.ac.cn