高 峰高玉平童明雷楊廷高趙成仕

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

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

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

(4西安科技大學(xué)應(yīng)用物理系西安710054)

1 引言

脈沖星計(jì)時(shí)分析是整個(gè)脈沖星領(lǐng)域研究的基石,在諸如強(qiáng)引力場中的廣義相對論效應(yīng)、檢驗(yàn)引力波、探測太陽系外行星以及提供計(jì)時(shí)標(biāo)準(zhǔn)等研究中[1?4],都需要高精度的脈沖星計(jì)時(shí)分析的結(jié)果.當(dāng)前,由于計(jì)時(shí)模型不完善及一些未能消除的誤差源的存在,限制了脈沖星計(jì)時(shí)分析精度進(jìn)一步提高,其影響主要體現(xiàn)在脈沖星計(jì)時(shí)殘差的波動性上.脈沖星計(jì)時(shí)殘差是脈沖星時(shí)與原子時(shí)比對的產(chǎn)物,其波動性主要來源于脈沖星計(jì)時(shí)觀測誤差、脈沖到達(dá)時(shí)間(Time of Arrival,TOA)轉(zhuǎn)換誤差、脈沖星本身內(nèi)部噪聲等不確定因素的影響.不同的誤差源對計(jì)時(shí)殘差的影響各不相同[5?7],分析計(jì)時(shí)殘差諸多波動特征,將有利于解讀上述各誤差的來源.參考文獻(xiàn)[8?10]已報(bào)道了利用指數(shù)模型、頻域方法、功率譜方法等實(shí)現(xiàn)計(jì)時(shí)殘差中低頻噪聲的度量及分類.頻域分析的最大優(yōu)點(diǎn)是通過傅立葉或小波變換能將波動信號的某些特征清楚地顯現(xiàn)出來.然而,上述變換所用的基底函數(shù)系和“窗口”分析都具有明顯的主觀選擇性.Huang等人于1998年提出了經(jīng)驗(yàn)?zāi)B(tài)分解(Empirical Mode Decomposition,EMD)方法用于分析非線性非平穩(wěn)時(shí)間序列[11].EMD方法具有良好的功能:對數(shù)據(jù)具有自適應(yīng)性,能用幾個(gè)內(nèi)在的本征模態(tài)函數(shù)(Intrinsic Mode Function,IMF)和一個(gè)剩余趨勢揭示時(shí)間序列的振蕩結(jié)構(gòu)和非平穩(wěn)性.EMD方法一經(jīng)提出,便在振動故障診斷、光學(xué)信號、濾波及噪聲消除等多個(gè)領(lǐng)域得到了廣泛應(yīng)用.

本文基于EMD方法分析時(shí)間序列的良好功能,將脈沖星計(jì)時(shí)殘差看作時(shí)間序列并采用EMD方法進(jìn)行分析,旨在揭示計(jì)時(shí)殘差內(nèi)在的振蕩結(jié)構(gòu),這有助于理解不同誤差源對脈沖星計(jì)時(shí)的影響;其次,脈沖星計(jì)時(shí)殘差的預(yù)報(bào)能反推各類誤差源影響的趨勢,本文進(jìn)一步利用EMD方法分解后的IMF分量或其組合,探討EMD方法在脈沖星計(jì)時(shí)殘差預(yù)報(bào)方面所具有的優(yōu)勢.

2 EMD方法的數(shù)學(xué)原理

EMD方法是分解信號序列的一種有效工具,它把時(shí)間序列分解為有限個(gè)IMF分量和一個(gè)剩余趨勢項(xiàng)[12].各分量應(yīng)滿足如下約束條件:(1)在整個(gè)資料集里,極值點(diǎn)的數(shù)目與穿零點(diǎn)的數(shù)目必須相等或者最多相差1個(gè);(2)由局部極大值所構(gòu)成的包絡(luò)線以及由局部極小值所構(gòu)成的包絡(luò)線的平均值為零.

EMD方法的關(guān)鍵技術(shù)問題是如何把一個(gè)非線性非平穩(wěn)時(shí)間序列(假定初始時(shí)間序列為X(t))分解為有限個(gè)IMF分量和一個(gè)趨勢項(xiàng).解決問題的方法是通過篩選,即:首先,通過3次樣條曲線聯(lián)結(jié)所在局部極大值得到上包絡(luò)線,用同樣的方法連接局部極小值形成下包絡(luò)線,上下兩條包絡(luò)線應(yīng)當(dāng)分別涵蓋全部的極大值和極小值;其次,取上下包絡(luò)線的均值線為m1(t).原序列X(t)與均值m1(t)之差為h1(t),即

這樣得到的h1(t)未必滿足上述IMF分量的約束條件.把h1(t)當(dāng)作初始時(shí)間序列,重復(fù)篩選過程

直到第k次篩選的結(jié)果符合要求,這樣就得到了第1個(gè)IMF分量C1.C1應(yīng)是信號的特征時(shí)間尺度最小的高頻IMF分量.接下來,把C1從X(t)中分離出去,得到剩余時(shí)間序列.剩余時(shí)間序列重復(fù)(1)–(2)式的篩選過程,直到剩余項(xiàng)成為一個(gè)單調(diào)函數(shù).于是得到一系列的IMF分量,即

(3)式中,C1,···,Cn是IMF的n個(gè)分量,r是剩余趨勢項(xiàng).

3 數(shù)據(jù)來源及分析

本文選取4顆毫秒脈沖星B1937+21、J0613–0200、J1012+5307、J1643–1224進(jìn)行分析,其相關(guān)TOA數(shù)據(jù)和脈沖星參數(shù)來源于Arzoumanian等[13]2015年公布的北美計(jì)時(shí)陣(North American Nanohertz Observatory for Gravitational Waves,NANOGrav)觀測的37顆毫秒脈沖星數(shù)據(jù).上述數(shù)據(jù)經(jīng)標(biāo)準(zhǔn)脈沖星數(shù)據(jù)處理流程擬合脈沖星天文參數(shù),隨后獲得的計(jì)時(shí)殘差數(shù)據(jù)作為本文分析的原始數(shù)據(jù).4顆毫秒脈沖星觀測跨度最長為9.2 yr,釋放的TOA數(shù)據(jù)計(jì)數(shù)多達(dá)11995點(diǎn)(部分基本參數(shù)見表1,P為脈沖星自轉(zhuǎn)周期).

表1 4顆脈沖星基本參數(shù)Table 1 The basic parameters of 4 pulsars

本文選取脈沖星數(shù)據(jù)的標(biāo)準(zhǔn)是基于以下考慮:首先,脈沖星計(jì)時(shí)模型參數(shù)精度的提高雖然有賴于多種因素,但觀測TOA計(jì)數(shù)越多,越有利于脈沖星計(jì)時(shí)模型參數(shù)的擬合,從而降低擬合誤差的引入;其次,采樣時(shí)間跨度較長,將有利于脈沖星計(jì)時(shí)殘差中紅噪聲成分的測量.所選的4顆脈沖星在文獻(xiàn)[13]中皆給出了殘差的紅噪聲模型,在本文分析中,其約9 yr觀測跨度將非常有利于解析低頻噪聲對計(jì)時(shí)殘差的影響.目前,脈沖星計(jì)時(shí)觀測并不是規(guī)則采樣,為了避免或減小因數(shù)據(jù)不規(guī)則采樣引入新的統(tǒng)計(jì)誤差,本文所選脈沖星計(jì)時(shí)殘差數(shù)據(jù)均做了規(guī)則化預(yù)處理.其處理過程以脈沖星J0613–0200為例進(jìn)行說明:根據(jù)表1第3列和第5列可知,脈沖星J0613–0200在約8.6 yr的觀測跨度內(nèi)有7651個(gè)計(jì)數(shù)點(diǎn),對應(yīng)的約化儒略日(Modified Julian Day,MJD)有較高的分辨率,即同一天可能對應(yīng)數(shù)個(gè)觀測點(diǎn).令同一天有觀測點(diǎn)時(shí)取樣1個(gè),其值為同一天內(nèi)數(shù)個(gè)點(diǎn)求平均,對應(yīng)MJD數(shù)值向下取整.這樣處理過程等同于1 d之內(nèi)脈沖星時(shí)與原子時(shí)比對1次取計(jì)時(shí)殘差,經(jīng)過處理后脈沖星J0613–0200計(jì)時(shí)殘差為171個(gè)點(diǎn).隨后進(jìn)行線性插值取樣,獲得均勻分布的同樣數(shù)目計(jì)時(shí)殘差點(diǎn)數(shù),即171個(gè)點(diǎn).對其他3顆脈沖星殘差做相同規(guī)則化預(yù)處理,預(yù)處理后的4顆脈沖星計(jì)時(shí)殘差作為本文框架內(nèi)分析的初始數(shù)據(jù),其殘差分布如圖1所示.

將脈沖星計(jì)時(shí)殘差數(shù)據(jù)作為信號序列應(yīng)用EMD方法進(jìn)行分解,將討論以下問題:

(1)各脈沖星計(jì)時(shí)殘差的波動由哪些時(shí)間尺度的振蕩模態(tài)變化構(gòu)成;

(2)哪些振蕩模態(tài)對整個(gè)計(jì)時(shí)殘差的變化起主要作用,不同脈沖星計(jì)時(shí)殘差波動有何共性及個(gè)性;

(3)不同振蕩模態(tài)時(shí)間尺度的變化與脈沖星計(jì)時(shí)測量誤差來源有何聯(lián)系.

圖2(a)–(h)分別表示脈沖星計(jì)時(shí)殘差數(shù)據(jù)X(t)、應(yīng)用EMD方法將X(t)進(jìn)行分解后得到的6個(gè)IMF分量C1–C6和剩余趨勢項(xiàng)r的分布.實(shí)線、虛線、點(diǎn)線、點(diǎn)虛線分別表示數(shù)據(jù)來源于脈沖星B1937+21、J0613–0200、J1012+5307、J1643–1224的殘差.整體來看,對比X(t)可以發(fā)現(xiàn):圖2(a)中4顆脈沖星殘差數(shù)據(jù)X(t)幅值大小漲落很不規(guī)則且有明顯的低頻特征.圖2(b)–(g)中分解之后C1–C6分量呈現(xiàn)出圍繞零均值線振蕩且局部極大值和極小值基本對稱的振蕩形式.它們的均值都為零,不隨時(shí)間變化,頻率變化較小,波形具有比X(t)規(guī)整、簡單、非平穩(wěn)性減弱等特點(diǎn).圖2(h)中剩余趨勢項(xiàng)r是近似單調(diào)變化的曲線.除此之外,分解之后通過不同IMF分量可以清晰地表明X(t)的多尺度變化,C1–C6分量時(shí)間尺度逐漸增大,亦即從高頻到低頻.上述特點(diǎn)都非常符合EMD方法分解時(shí)序信號的特征.

圖1 4顆脈沖星計(jì)時(shí)殘差分布.縱坐標(biāo)是計(jì)時(shí)殘差,橫坐標(biāo)是MJD的跨度.實(shí)線、虛線、點(diǎn)線、點(diǎn)虛線分別代表脈沖星B1937+21、J0613–0200、J1012+5307、J1643–1224的殘差數(shù)據(jù).Fig.1 The timing residuals of 4 pulsars.The vertical axis stands for the timing residual,and the horizontal axis presents span of MJD.The solid,dashed,dotted,and dot-dashed lines represent B1937+21,J0613–0200,J1012+5307,and J1643–1224,respectively.

考察IMF分量對于殘差數(shù)據(jù)X(t)的相對重要性時(shí),可根據(jù)各分量振幅變化的量級(如圖2的縱坐標(biāo)數(shù)值)比較;也可以參照諧波分析,用方差貢獻(xiàn)表現(xiàn)各分量對原序列的重要性.考慮到標(biāo)準(zhǔn)差與數(shù)據(jù)的量綱一致,在描述數(shù)據(jù)波動性時(shí)更明顯,且方差與標(biāo)準(zhǔn)差很容易轉(zhuǎn)換,本文采用標(biāo)準(zhǔn)差貢獻(xiàn)來表征各分量對原序列的重要性.表2第2列是各脈沖星初始?xì)埐钚蛄械臉?biāo)準(zhǔn)差值,其中脈沖星J1643–1224標(biāo)準(zhǔn)差最大,為2.5658μs. 脈沖星J0613–0200標(biāo)準(zhǔn)差值最小,為0.7575 μs. 第3–9列是IMF各分量C1–C6及剩余趨勢項(xiàng)r對各脈沖星殘差序列波動性的標(biāo)準(zhǔn)差貢獻(xiàn).相對于后3顆脈沖星各分量貢獻(xiàn)而言,脈沖星B1937+21各分量中標(biāo)準(zhǔn)差最大的3個(gè)依次為C6,r和C4,它們的貢獻(xiàn)分別為0.7482μs,0.5321μs和0.4748μs. 尤其是C6,比C1–C3大了數(shù)倍,說明C6的貢獻(xiàn)對X(t)振蕩起決定性作用,即X(t)中最顯著的低頻特征由C6反映出來.脈沖星J0613–0200除C3–C4分量偏小外,其他分量及趨勢項(xiàng)標(biāo)準(zhǔn)差值相當(dāng).脈沖星J1012+5307中C2最大,為1.0623μs,除C6外,其他各項(xiàng)隨時(shí)間尺度的增大,標(biāo)準(zhǔn)差貢獻(xiàn)都是下降的.而脈沖星J1643–1224各IMF分量和趨勢項(xiàng)標(biāo)準(zhǔn)差值相差不大.

圖2 圖(a)–(h)分別是4顆脈沖星計(jì)時(shí)殘差X(t)、IMF分量C1–C6和剩余趨勢項(xiàng)r的分布.實(shí)線、虛線、點(diǎn)線、點(diǎn)虛線分別表示數(shù)據(jù)來源于脈沖星B1937+21、J0613–0200、J1012+5307、J1643–1224的殘差.Fig.2 The panels(a)–(h)indicate pulsar timing residuals X(t),C1–C6of the intrinsic mode function,and trend term r,respectively.The solid,dashed,dotted,and dot-dashed lines represent B1937+21,J0613–0200,J1012+5307,and J1643–1224,respectively.

表2 各IMF分量對4顆脈沖星計(jì)時(shí)殘差波動性的標(biāo)準(zhǔn)差貢獻(xiàn)Table 2 The contributions of IMF components to the standard deviations of the volatility for 4 pulsars timing residuals

綜上所述,我們可以看到:采用EMD方法分析能較好地用幾個(gè)內(nèi)在IMF分量和一個(gè)剩余趨勢揭示各脈沖星殘差序列的振蕩結(jié)構(gòu)和非平穩(wěn)趨勢.尤其是結(jié)合各分量的標(biāo)準(zhǔn)差貢獻(xiàn)可以清楚地看出某個(gè)或幾個(gè)IMF分量在原殘差序列中所起的作用.本文所選的4顆脈沖星其共性都具有明顯的低頻波動(紅噪聲),這與文獻(xiàn)[13]能用紅噪聲模型描述殘差中的紅噪聲成分是一致的.但通過EMD方法分解后分析可以看出,4顆脈沖星殘差的波動特征仍有較大的差異.以脈沖星B1937+21為例,其低頻成份貢獻(xiàn)非常明顯,這表明其殘差序列的頻譜分析在低頻端必定有明顯的能量包存在.而其他3顆脈沖星殘差低頻成分則沒有這么顯著.另外,4顆脈沖星的差異也表現(xiàn)在各自的剩余趨勢項(xiàng)r上,結(jié)合圖2(h)可以看出:脈沖星B1937+21和J1643–1224近10 yr來殘差序列總體上是上升趨勢,而脈沖星J0613–0200和J1012+5307是下降趨勢.最后,如果簡單地將脈沖星計(jì)時(shí)殘差區(qū)分為白噪聲和紅噪聲兩種成分,白噪聲主要來源于計(jì)時(shí)觀測的隨機(jī)誤差,具有高頻、隨機(jī)性強(qiáng)等特征.而紅噪聲來源廣泛,脈沖星本身內(nèi)部噪聲、傳播誤差、脈沖星到達(dá)時(shí)間轉(zhuǎn)換誤差等都可累積形成低頻噪聲.單個(gè)或多個(gè)因素組合形成周期性特征并且對殘差波動貢獻(xiàn)較強(qiáng)時(shí),理論上EMD方法都能在時(shí)域分析中做到較好的辨識.所以,EMD方法是研究不同誤差來源對脈沖星計(jì)時(shí)影響的一種重要的工具,可與計(jì)時(shí)殘差的頻譜分析形成交叉證認(rèn).

4 ARMA模型及其實(shí)現(xiàn)

自回歸滑動平均模型(Auto-Regressive Moving Average Model,ARMA)是時(shí)間序列分析與前推預(yù)測領(lǐng)域中一類常用的重要方法[14].模型形式ARMA(p,q)中包含了p階自回歸項(xiàng)和q階移動平均項(xiàng).ARMA(p,q)模型可以表示為

其中,Xt是平穩(wěn)或準(zhǔn)平穩(wěn)的時(shí)間序列,φi和θj分別是自回歸參數(shù)和滑動平均參數(shù),εt為均值為零、方差為σ2的正態(tài)白噪聲,可認(rèn)為是隨機(jī)干擾項(xiàng).

(4)式中,若Xt僅與t時(shí)刻的白噪聲εt和以前的變化Xt?i有關(guān),即

此時(shí),(5)式稱為p階自回歸(Auto-Regressive,AR)模型,可記作AR(p)或者ARMA(p,0),其中p為模型階數(shù).

若Xt僅與t時(shí)刻的白噪聲εt及以前的白噪聲εt?j存在關(guān)聯(lián),即

則(6)式稱為q階滑動平均(Moving Average,MA)模型,記作MA(q)或是ARMA(0,q).所以,ARMA(p,q)模型可看作是AR模型與MA模型的有機(jī)綜合.

利用ARMA模型進(jìn)行時(shí)間序列的預(yù)報(bào)通常由以下4步確認(rèn)或計(jì)算:

(1)模型識別 對于一個(gè)平穩(wěn)時(shí)間序列,可以通過該序列的自相關(guān)函數(shù)(Auto-Correlation Function,ACF)和偏相關(guān)函數(shù)(Partial Correlation Function,PCF)的“拖尾”(Trailing Lag,相關(guān)函數(shù)隨延遲量增加逐漸歸于零)或“截尾”(Truncation Lag,延遲量大于某值后相關(guān)函數(shù)恒等于零)特性來合理判斷序列更符合哪個(gè)模型,其判斷準(zhǔn)則可見表3.

(2)模型定階模型類別確立后,模型階次可由信息論準(zhǔn)則、傳遞函數(shù)準(zhǔn)則或是最終預(yù)測誤差準(zhǔn)則等多種方式確定模型的最優(yōu)階次[15].本文采用Akaike信息論法則(Akaike Information Criterion,AIC)來選擇模型的最佳階次.

(3)模型參數(shù)估計(jì)模型參數(shù)估計(jì)方法有多種,本文采用Yule-Walker方程估計(jì)法.

(4)預(yù)報(bào)模式選擇預(yù)報(bào)模式同樣多樣且各具優(yōu)點(diǎn),本文采用遞推模式進(jìn)行預(yù)報(bào),其特點(diǎn)是只需一次建模計(jì)算、計(jì)算量小、可減少建模的復(fù)雜性可能引入的多種誤差.

表3 ARMA模型識別準(zhǔn)則Table 3 ARMA identification criterion

5 脈沖星計(jì)時(shí)殘差的預(yù)報(bào)

在第4節(jié)中,我們討論了ARMA模型相關(guān)原理及其對平穩(wěn)時(shí)間序列進(jìn)行預(yù)報(bào)的步驟.本節(jié)中將進(jìn)一步討論應(yīng)用EMD方法對脈沖星計(jì)時(shí)殘差預(yù)報(bào)的優(yōu)勢.需要說明的是,時(shí)序數(shù)據(jù)的回歸預(yù)報(bào)對數(shù)據(jù)的平穩(wěn)性有較高的要求,脈沖星J0613–0200殘差數(shù)據(jù)相較于其他3顆脈沖星有更好的準(zhǔn)平穩(wěn)性.所以,本節(jié)中僅對脈沖星J0613–0200應(yīng)用ARMA模型進(jìn)行計(jì)時(shí)殘差序列的預(yù)報(bào),各參與預(yù)報(bào)的數(shù)據(jù)均通過了用于單位根檢驗(yàn)的擴(kuò)展迪基-福勒檢驗(yàn)(Augmented Dickey-Fuller test)[16?17].

開展計(jì)時(shí)殘差序列預(yù)報(bào)的思路是:首先截取初始數(shù)據(jù)前部采用直接建模預(yù)報(bào)和對截取數(shù)據(jù)應(yīng)用EMD方法分解后建模間接預(yù)報(bào)兩種模式;隨后兩種預(yù)報(bào)結(jié)果與初始數(shù)據(jù)后部(真實(shí)觀測數(shù)據(jù))進(jìn)行比對,以檢驗(yàn)兩種預(yù)報(bào)方式的優(yōu)劣及其精度.其具體步驟如下:

(1)直接預(yù)報(bào) 取脈沖星J0613-0200殘差數(shù)據(jù)(共171個(gè)點(diǎn))的前161個(gè)點(diǎn)建模并預(yù)報(bào)后10個(gè)點(diǎn);

(2)間接預(yù)報(bào) 取殘差數(shù)據(jù)的前161個(gè)點(diǎn)并應(yīng)用EMD方法分解,利用分解后的IMF分量或其組合建模并預(yù)報(bào)后10個(gè)點(diǎn).

采用直接預(yù)報(bào)模式時(shí),根據(jù)表3中ARMA模型識別準(zhǔn)則,發(fā)現(xiàn)建模數(shù)據(jù)適合AR(p)模型.隨后逐次進(jìn)行模型定階、參數(shù)估計(jì),并采用遞推模式進(jìn)行預(yù)報(bào).采用EMD方法分解進(jìn)行間接預(yù)報(bào)有兩點(diǎn)考慮:一是在比對預(yù)報(bào)結(jié)果時(shí)為體現(xiàn)兩種預(yù)報(bào)模式的優(yōu)劣,本文中將控制間接預(yù)報(bào)模型、參數(shù)估計(jì)和預(yù)報(bào)模式與直接預(yù)報(bào)流程中選擇盡量相同或類似;二是EMD方法分解后的IMF高頻分量會出現(xiàn)隨機(jī)性太強(qiáng)不能提取有益信息進(jìn)行預(yù)報(bào),或是低頻分量趨勢性太強(qiáng)不符合平穩(wěn)時(shí)間序列的要求.本文采用各分量組合重構(gòu)的方式解決上述問題.如圖3所示,上圖為C1與C3相加重構(gòu),中圖為C2加C5重構(gòu),下圖為C4與趨勢項(xiàng)r重構(gòu).這類組合方式并沒有檢驗(yàn)是否最優(yōu),但各重構(gòu)時(shí)序均通過了單位根檢驗(yàn),并根據(jù)表3中ARMA模型識別準(zhǔn)則,均適合AR(p)模型.類似于直接預(yù)報(bào)流程,3個(gè)重構(gòu)時(shí)序分別給出各自的預(yù)報(bào)結(jié)果,最后疊加給出間接預(yù)報(bào)的結(jié)果.

圖4中給出了直接預(yù)報(bào)數(shù)據(jù)、間接預(yù)報(bào)數(shù)據(jù)與初始數(shù)據(jù)(后10個(gè)點(diǎn),下同)比對的結(jié)果,以此檢驗(yàn)兩種預(yù)報(bào)模式的優(yōu)劣及其精度.4個(gè)標(biāo)注參量X(t)、AR(X)、AR(C)和X?(t)分別表示初始數(shù)據(jù)、直接預(yù)報(bào)數(shù)據(jù)、間接預(yù)報(bào)數(shù)據(jù)和初始數(shù)據(jù)平滑后的值.

圖3 不同IMF分量的組合重構(gòu).上圖是C1和C3累加,中圖是C2和C5累加,下圖是C4和r累加.Fig.3 The reconstruction of different IMF components.The top panel stands for the sum of C1and C3,the middle panel is the sum of C2and C5,and the bottom panel presents the sum of C4and r.

圖4 脈沖星計(jì)時(shí)殘差的預(yù)報(bào)數(shù)據(jù)和初始數(shù)據(jù)的比對.虛線和點(diǎn)線分別表示直接預(yù)報(bào)數(shù)據(jù)和間接預(yù)報(bào)數(shù)據(jù).上下圖實(shí)線分別表示初始數(shù)據(jù)和平滑后的初始數(shù)據(jù).Fig.4 The comparison between the forecast data and the initial data of pulsar timing residuals.The dashed and dotted lines indicate the direct and indirect forecast data,respectively.The solid lines in the upper and lower panels represent the initial data and the smoothed initial data,respectively.

分析圖4可發(fā)現(xiàn):上圖顯示初始數(shù)據(jù)X(t)隨機(jī)起伏較大,直接預(yù)報(bào)數(shù)據(jù)AR(X)和間接預(yù)報(bào)數(shù)據(jù)AR(C)較為平滑.檢驗(yàn)預(yù)報(bào)結(jié)果好壞最直接的方法是比對預(yù)報(bào)數(shù)據(jù)與初始數(shù)據(jù)接近的程度.為消除初始數(shù)據(jù)隨機(jī)性因素的干擾,下圖將初始數(shù)據(jù)進(jìn)行平滑后與預(yù)報(bào)數(shù)據(jù)進(jìn)行比對.圖5中顯示的是兩種預(yù)報(bào)數(shù)據(jù)AR(X)、AR(C)分別與平滑數(shù)據(jù)X?(t)作差后的結(jié)果.從接近程度上看,AR(C)在短期(約4個(gè)月)范圍內(nèi)要明顯優(yōu)于AR(X).圖5中更長期(大于4個(gè)月)預(yù)報(bào)時(shí),直接預(yù)報(bào)數(shù)據(jù)AR(X)更接近平滑數(shù)據(jù)X?(t),但考慮初始數(shù)據(jù)末尾存在端部效應(yīng)等因素,上述結(jié)果不能形成明顯結(jié)論.除此之外,初始數(shù)據(jù)的1倍標(biāo)準(zhǔn)差約為0.37μs,短期預(yù)報(bào)時(shí),間接預(yù)報(bào)數(shù)據(jù)明顯落在初始數(shù)據(jù)1倍標(biāo)準(zhǔn)差之內(nèi),具有較高的精度.所以,綜上所述,采用EMD方法分解后對脈沖星J0613–0200計(jì)時(shí)殘差進(jìn)行短期預(yù)報(bào),其預(yù)報(bào)模式合理有效,且預(yù)報(bào)結(jié)果明顯優(yōu)于直接預(yù)報(bào).

圖5 脈沖星計(jì)時(shí)殘差的預(yù)報(bào)數(shù)據(jù)和平滑后初始數(shù)據(jù)的比對.實(shí)線和短劃線分別來自直接預(yù)報(bào)數(shù)據(jù)和間接預(yù)報(bào)數(shù)據(jù)與平滑后初始數(shù)據(jù)作差的絕對值.Fig.5 The comparison between the forecast data and smoothed initial data of pulsar timing residuals.The solid line indicates the absolute value of the difference between direct forecast data and the smoothed initial data,and the dashed line represents the absolute value of the difference between the indirect forecast data and the smoothed initial data.

6 結(jié)論

分析脈沖星計(jì)時(shí)殘差的波動特征有利于研究各類誤差源對脈沖星計(jì)時(shí)精度的影響.EMD方法分析具有對時(shí)序數(shù)據(jù)較強(qiáng)的自適應(yīng)性,能用幾個(gè)內(nèi)在的本征模態(tài)函數(shù)和一個(gè)剩余趨勢揭示時(shí)間序列的振蕩結(jié)構(gòu)和非平穩(wěn)性等特點(diǎn),非常適合脈沖星計(jì)時(shí)殘差數(shù)據(jù)的分析.本文采用EMD方法對脈沖星B1937+21、J0613–0200、J1012+5307、J1643–1224的計(jì)時(shí)殘差振蕩模式進(jìn)行了分析,并進(jìn)一步探究了應(yīng)用EMD方法在計(jì)時(shí)殘差預(yù)報(bào)方面所具有的優(yōu)勢.通過以上研究,我們發(fā)現(xiàn):

(1)脈沖星的計(jì)時(shí)殘差的波動可以分解為多個(gè)振蕩模態(tài)和一個(gè)趨勢項(xiàng),不同的振蕩模態(tài)清晰地表明了計(jì)時(shí)殘差波動的多時(shí)間尺度的變化.通過比較4顆脈沖星各振蕩模態(tài)的方差貢獻(xiàn),發(fā)現(xiàn)脈沖星B1937+21低頻成分對殘差波動貢獻(xiàn)最大,其余3顆脈沖星則不是很顯著.這種脈沖星殘差序列波動的差異性也表現(xiàn)在分解后的剩余趨勢項(xiàng)上.

(2)EMD方法具有在時(shí)域上辨識殘差序列高頻成分與低頻成分的能力,在研究不同誤差來源對脈沖星計(jì)時(shí)殘差影響領(lǐng)域,可與殘差的頻譜分析形成交叉證認(rèn).

(3)通過分析單顆脈沖星J0613–0200計(jì)時(shí)殘差序列回歸預(yù)報(bào),揭示在短期預(yù)報(bào)范圍內(nèi)采用EMD方法將計(jì)時(shí)殘差序列分解后再預(yù)報(bào)要優(yōu)于直接預(yù)報(bào)結(jié)果.究其原因,可認(rèn)為EMD方法是將時(shí)間序列平穩(wěn)化的一種過程,這將有利于提高模型預(yù)報(bào)精度.

致謝此研究利用了北美計(jì)時(shí)陣的觀測數(shù)據(jù),并特別感謝中國科學(xué)院國家授時(shí)中心雷雨博士對自回歸模型的有益討論.