XPNAV-1衛星先期發布數據的計時分析

朱鴻旭，童明雷，楊廷高，趙成仕，李 琳，高玉平

(1. 中國科學院國家授時中心，西安 710600；2. 中國科學院時間頻率基準重點實驗室，西安 710600；3. 中國科學院大學，北京 100049)

0 引言

脈沖星是自轉非常穩定的自然天體，通常具有從射電一直到Gamma射線的多波段電磁輻射。而毫秒脈沖星由于自轉穩定性更高，更值得我們關注。對毫秒脈沖星的計時觀測，在科學和工程上有諸多應用。利用地面上的射電望遠鏡對脈沖星做長期計時觀測，可以開展脈沖星時間尺度建立[1-3]、探測低頻引力波[4-7]、測量太陽系行星質量[8]等方面的研究工作。利用X射線脈沖星可以實現各類空間飛行器自主導航，尤其適用于深空探測飛行器導航。理論上，若星載原子鐘精度足夠高，同時觀測三顆脈沖星即可對航天器展開絕對導航。在觀測多顆脈沖星相對困難的情況下，還可以利用脈沖星進行相對導航[9]。目前，深空探測飛行器導航主要依賴地面跟蹤測量技術，如美國的深空跟蹤網(Deep space network,DSN)或甚長基線干涉測量(Very long baseline interferometry,VLBI)技術。由于深空探測飛行器距離地球遙遠，信號延遲大(如1AU距離，信號延遲約8.5 min，探測冥王星飛行器信號延遲達到約4 h)，不能做到實時定位測量。另外，地面跟蹤測量精度隨著飛行器距離增大而降低。脈沖星作為宇宙中的自然天體，當其自轉參數(自轉頻率及其導數)與天體測量參數(方向，距離與自行)經過精確測定后，可以為太陽系范圍內各類飛行器進行導航，其導航精度與飛行器距離無關。另外，脈沖星也有可能用作衛星導航的備用或輔助導航手段以及某些特殊情況下的自主導航[10-11]。2004年Sheikh[12]、Sala等[13]基于現代衛星導航系統體制的思想，各自論證了基于X射線脈沖星航天器自主導航的可行性。自此，各國脈沖星自主導航的相關研究進入快軌道[14-19]。尤其隨著美國NICER/SEXTANT項目成功運行，脈沖星導航的可實現性與可操作性受到越來越廣泛的認可[20]。我國在X射線脈沖星導航領域近幾年也發展迅速。硬X射線調制望遠鏡(HXMT/慧眼衛星)、天宮二號搭載的POLAR、X射線脈沖星導航試驗衛星(XPNAV-1)等平臺均在進行X射線觀測，開展了脈沖星導航試驗，并取得了一系列進展和成果[21-24]。對于脈沖星導航而言，Crab脈沖星是這些X射線觀測任務的一個主要目標源。

Crab脈沖星是公元1054年超新星爆發的產物[25]，它的X射線波段輻射在所有脈沖星中最強。由于Crab脈沖星比較年輕，內部物理過程仍頻繁發生，導致其自轉不穩定，需要每月甚至更頻繁地更新自轉參數。但由于Crab脈沖星具有全波段的輻射，且流量很強，便于觀測、分析。因此，利用衛星對Crab脈沖星展開持續觀測，進行脈沖星導航試驗仍是很有必要的。

目前我國的天基脈沖星導航試驗均基于人造地球衛星開展，今后可能會基于其他天體衛星(如：人造火星衛星)展開試驗，并應用于深空探測飛行器。對于人造衛星而言，若利用脈沖星對其進行自主導航，目前比較成熟的是基于軌道動力學[21]與X射線脈沖星計時觀測相結合的Kalman濾波[26-27]的方法。Kalman濾波的狀態方程包括軌道動力學模型及星載原子鐘的狀態方程；Kalman濾波的測量方程即為脈沖星計時觀測[28]。更重要的是建立衛星的軌道動力學模型，一般包括地球(火星)中心引力以及對應非球形攝動力模型。若對地球衛星進行導航，第三體主要考慮太陽與月亮的影響。對于低軌道衛星，還應該建立大氣阻力模型。由此可見，若對人造衛星進行脈沖星導航，對地球(火星)重力場展開精確測量顯得尤為重要[29-30]。高精度動力學模型有利于衛星軌道預報。但動力學模型的誤差具有累積性，會導致衛星軌道位置偏離越來越大。因而需要依賴外部觀測手段，如脈沖星計時觀測，實時測量衛星的軌道位置。設計合適的Kalman濾波器，在濾波器內將動力學預報與脈沖星計時觀測結合起來，從而實現具有一定精度的衛星軌道位置的預報與更新。如果知道衛星的多普勒頻移信息，可以將其與前述信息進行融合以提高導航精度[31]。根據XPNAV-1衛星軌道參數特征，理論上，利用軌道動力學模型與Crab脈沖星觀測，可以驗證Kalman濾波方法精度。但XPNAV-1觀測Crab脈沖星35組數據時間跨度為31天，如果由每組觀測，計算得到一個TOA(TimE-of-artival)，在刪除錯誤觀測后，實際上平均每天約只有一個有效TOA測量值。1天內衛星繞衛星軌道運動大約近20圈。利用如此稀疏的測量數據，估計難以用Kalman濾波方法直接檢驗XPNAV-1導航精度。而在測量方程中，TOA的測量精度和脈沖星計時模型精度是影響脈沖星導航精度的兩個重要因素。如果脈沖到達時刻的測量和脈沖星計時模型的精度都很高，那么二者之差(計時殘差)也會很小。因此，計時殘差的大小對脈沖星導航的精確度有著根本的影響：計時殘差的標準差σ每增大10 μs，導航精度至少下降3 km。不僅如此，獲得X射線脈沖星的計時殘差，是構建天基授時系統的基礎，對建立綜合脈沖星時同樣有著深遠的意義。

鑒于此，本文將對XPNAV-1先期發布的35組Crab脈沖星觀測數據[32]進行計時精度分析，每組數據觀測積分時間為50 min左右。XPNAV-1搭載了兩種不同體制的X射線探測器，工作能譜為0.5～10 keV[22]。首先將觀測數據時間尺度由UTC(Universal time coordinated)時間尺度轉換為質子力學時(TDB)時間尺度，同時采用英國Jordrell Bank臺發布的Crab脈沖星星歷[33-34]與美國JPL發布的DE200歷表[35]將光子到達衛星時刻轉換至太陽系質心(Solar system barycenter,SSB)。隨后利用35組數據建立積分脈沖輪廓與標準脈沖輪廓并進行互相關處理，得到觀測的脈沖到達時刻(TOA)。再與Crab脈沖星模型參數預報的TOA作差得到計時殘差。

1 光子到達時刻轉換與脈沖輪廓建立

對于X射線脈沖星計時觀測而言，要獲得TOA，通常需要先將光子達到航天器時刻轉換到達到SSB的時刻，然后進行周期折疊建立積分脈沖輪廓，再與標準脈沖輪廓進行互相關，得到TOA。當然，也可以采用美國SEXTANT脈沖星導航試驗觀測脈沖星的處理方法。該方法利用最大似然估計，直接由觀測記錄的光子到達時刻序列數據，計算得到脈沖星脈沖到達探測器時刻與多普勒頻移兩個參數[36]。鑒于目前XPNAV-1公開發布的數據脈沖信噪比有限，因此本文仍采取將光子到達時刻轉換到SSB并進行輪廓折疊的方法。

先期發布的數據中，每組觀測包括衛星接收到Crab光子的時刻文件和遙測的軌道文件。由于光子到達衛星時刻與進行軌道遙測的時刻并不重合，因此需要將軌道位置與速度內插到光子到達衛星時刻。利用觀測期間對應的Crab脈沖星星歷與DE歷表，將光子到達時刻由衛星轉換到SSB。轉換方程為：

Δ☉=ΔE☉+ΔR☉+ΔS☉

(1)

式中：Δ☉表示光子到達時刻轉換到SSB產生的總延遲量；ΔE☉表示Einstein延遲，表示時間尺度轉換產生的延遲量；ΔR☉表示Roemer延遲，代表幾何真空延遲；ΔS☉表示Shapiro延遲，代表時空彎曲引起的延遲。Jordrell Bank臺80年代就開始發布Crab脈沖星星歷。為了保證前后數據歸算采用的參數一致，目前其給出的Crab脈沖星星歷仍基于DE200歷表與TDB時間尺度，且不包含視差與自行。因而在作時間轉換時也必須與Crab脈沖星星歷保持一致，使用DE200歷表與TDB時間尺度，同樣也不考慮視差與自行，以避免產生額外的誤差。觀測文件中光子達到時刻是用UTC時間尺度來記錄的，需轉換為TDB時間尺度。其他部分的展開式可參照文獻[37]，下文主要給出UTC到TDB的時間尺度轉換表達式：

ΔE☉=ΔTAI-UTC+ΔTT-TAI+ΔTCB-TT+ΔTDB-TCB

(2)

進一步展開，有[38-39]：

ΔTAI-UTC=l

(3)

ΔTT-TAI=32s.184+δ

(4)

(5)

ΔTDB-TCB=-LBΔtTCB

(6)

(7)

(8)

建立積分脈沖輪廓時，將每組觀測第一個光子到達SSB的時刻作為起點(即積分脈沖輪廓的相位零點)，利用式(8)計算該組觀測每個光子到達SSB時刻對應的Crab脈沖星自轉相位φ1i，φ2i，…，φni，…。其中i表示某組觀測，n表示該組觀測的第n個光子。將每組觀測的所有光子的相位都減去該組第一個光子相位φ1i再進行折疊。取歸一化自轉相位的小數值即為該光子所對應積分脈沖輪廓的相位值。將一個自轉周期等分為N份，即為N個子相位間隔。N個子相位間隔與對應的光子數構成的直方圖就是該組測量得到的積分脈沖輪廓。

標準脈沖輪廓由于要獲得比積分脈沖輪廓更高的信噪比，因此，需利用全部可利用的觀測資料進行折疊。折疊時，除統一以Crab脈沖星星歷參數定義的相位零點為起點外，其余與積分脈沖輪廓折疊方法相同，且必須同樣取N個相位間隔。折疊后的35組輪廓相加即得標準脈沖輪廓。圖1、圖2是當N=256時，得到的一組積分脈沖輪廓和利用全部35組觀測資料得到的標準脈沖輪廓。

圖1 某一組積分脈沖輪廓Fig.1 A set of integrated pulse profile

圖2 標準脈沖輪廓Fig.2 Standard pulse profile

一般定義脈沖輪廓的主峰頂點對應的時刻標記為TOA。Crab脈沖星星歷定義脈沖輪廓的主峰頂點為相位零點，即式(8)中的φ0=0。由圖2可見，標準脈沖輪廓主峰頂點并不在相位零點，而是略有偏差，該相位差Δφ對應的時間差Δtφ約1.38 ms。這是由于兩種不同的觀測頻段引起的零點相位差。因為Crab脈沖星在射電波段流量更強，觀測得到的脈沖信噪比更高，故Jodrell Bank提供的星歷基于射電波段給出。由于Crab脈沖星的輪廓隨觀測頻段變化而有所不同[41]，當獲得Crab脈沖星星歷所使用的觀測波段(射電)與衛星探測光子的頻段(X射線)不一致時，利用星歷折疊出的脈沖輪廓主峰頂點會產生相位差。一般認為，不同波段的光子是由脈沖星不同位置輻射出的[42]。因此，當用X射線脈沖星導航，而其星歷參數來自于射電波段觀測的情況下，相位零點差的精確測量是非常必要的。

2 獲得觀測TOA與計時殘差

每組積分脈沖輪廓分別與標準脈沖輪廓進行互相關計算可以獲得二者之間的相位差，進而獲得時間差，確定TOA。離散傅里葉變換(DFT)法是一種常用的相關方法[43]。由于DFT法不需要知道脈沖的統計信息，且計算精度只與脈沖輪廓的信噪比有關，而與子相位間隔無關。鑒于目前衛星觀測的信噪比較低，導致N難以取很高的值，脈沖直接包含的頻域信息也較少[44]。且由于未獲取大量觀測數據，較難得知脈沖的確切統計信息。因此，采用DFT法進行計算是合適的。

設s(t)，p(t)分別為按照上文方法求得的Crab脈沖星的標準脈沖輪廓和一組積分脈沖輪廓。其中0≤t≤P，P為Crab脈沖星的自轉周期。二者關系可表示為[43]：

p(t)=a+bs(t-τ)+g(t)

(9)

式中:a為幅值差，b為比例因子，τ為時間延遲，g(t)為高斯白噪聲。τ是最為關心的待求量。對上式兩邊分別作DFT，根據DFT的線性性質，有：

Pkexp(iθk)=aN+bSkexp[i(φk+2πνkτ)]+Gk

(10)

χ2(τ)=

(11)

(12)

由于待求量τ包含在三角函數內，難以直接求解。故一般用迭代方法求解。初值可以選擇0值。為提高迭代效率，最好選擇時域相關計算的τ。根據式(12)計算出的每組的τ，再加上每組觀測第一個光子到達時刻即代表每組數據觀測的TOA。

3 計時殘差分析

脈沖星觀測TOA與根據脈沖星星歷參數預報的TOA之差即為計時殘差。第2節計算得到了每組觀測得到的TOA，將觀測TOA代入式(8)中的t,利用Jodrell Bank給出的Crab自轉參數，可計算得到觀測TOA對應的相位殘差。將相位殘差除以觀測時刻Crab的自轉頻率就得到以時間單位表示的計時殘差，即自轉模型預報TOA與觀測TOA之差。

利用上述算法，歸算N=256時35組觀測的計時殘差，結果如圖3所示。由于脈沖的信噪比有限，有些組的觀測出現了迭代不收斂的現象，計算出的計時殘差很大。去掉在所有觀測的計時殘差中迭代不收斂的觀測組，再去掉超過標準差±3σ以上的觀測組構成有效觀測數據。即去掉序號為2，5，9，14，15，18，25的共7組觀測，有效觀測為剩余的28組觀測數據。減掉Δφ引起的Δtφ(1.38 ms)之后，計時殘差的結果如圖4所示。注意，圖3和圖4是擬合前計時殘差，即不對脈沖星計時模型參數做擬合。

圖3 所有觀測數據的計時殘差Fig.3 Timing residuals of all observed data

需要指出的是，當利用脈沖星數據處理軟件Tempo2[46]進行光子到達時刻轉換時，一般使用TCB時間尺度進行計算，與Crab脈沖星星歷采用的時間尺度不一致，這樣就產生了系統誤差[47]。而本文在做時間尺度轉換時，未利用TEMPO2軟件，而是采用了Python環境下的Astropy[47]等包，并采用了TDB時間尺度計算，避免了由時間尺度不一致導致的系統誤差。如式(6)所示，這種誤差以線性趨勢項為主。

基于Crab脈沖星的導航試驗還必須使用定期更新的星歷，以及與歸算星歷時一致的太陽系歷表。英國Jodrell Bank提供的Crab脈沖星星歷定期更新，且在目前精確度最高。然而，該星歷基于DE200歷表給出，其發布時間很早，誤差較大。所提供的各天體的位置與較新的歷表(如DE430)約有數百公里的差，且帶有明顯的趨勢項。若使用DE430歷表歸算，盡管其自身的精確度高于DE200歷表，但由于與歸算脈沖星星歷使用的歷表不同，采用同樣的數據處理方法，結果得到的計時殘差比采用DE200歷表的計時殘差反而彌散度更大。另外，英國Jodrell Bank天文臺提供的Crab脈沖星星歷沒有給出Crab脈沖星的自行與視差參數，在數據處理時又試驗性地增加了該脈沖星的自行與視差參數(取自于國際脈沖星數據庫)，結果也表明殘差彌散變大。這證明，在脈沖星導航數據處理時，采用的脈沖星星歷參數(包括自轉參數與位置等天體測量參數)必須是內部自洽的一套參數，不能人為地增加或減少某些個別參數，否則會降低精度。采用與脈沖星星歷相一致的DE歷表與質心坐標時也同樣是非常重要的，否則會產生系統誤差。對于利用Crab脈沖星的導航試驗研究，盡快利用中國的射電計時觀測，實時提供高精度星歷表是必要的。地面射電與空間X射線同步觀測研究是非常有意義的。

上文以N=256為例，進行了試驗，得到了上述結論。但N的取值未必是恰當的。依次取N=128～512，重復上述處理過程，可得到N與殘差的RMS，及無效觀測組數之間的關系。由圖5可以看出，當N逐漸增加時，殘差的RMS先線性降低，再基本保持平穩，后幅度劇烈變化。由于輪廓折疊等價于進行低通濾波，而觀測的信噪比較低，于是當N過大時，濾波效果不明顯，導致輪廓信噪比迅速降低；當N過小時，又會導致輪廓過平滑，導致頻域信息缺失過多，影響相關計算結果。這與由圖5得出的信息是一致的。因此，從實際應用的角度出發，對XPNAV-1的數據而言，取N=256是恰當的。

如圖6所示，無效觀測組數大致呈線性增加的趨勢。盡管每次觀測積分時間均為50 min左右，但由于衛星所在位置不同，受到的遮擋與干擾因素也不盡相同。故每次觀測時的脈沖信噪比不一致。當N持續增大時，脈沖信噪比較差的觀測濾波效果有限，從而導致相關計算過程中出現不收斂現象。對于XPNAV-1而言，需持續改進觀測方案，提升有效觀測時長，并盡可能減少背景光子的影響。同時，在積累了大量數據后，還可以利用先驗知識，對TOA測量方法加以改進。這些都是提高TOA測量精度，增加有效觀測數可行的手段。

圖5 子相位間隔數與殘差均方根(RMS)的關系Fig.5 The relationship between the number of sub-phase intervals and the RMS of timing residuals

圖6 相位間隔數與無效觀測組數的關系Fig.6 The relationship between the number of sub-phase intervals and the number of unavailable observed data set

圖7 兩種取權方法的一組輪廓比較Fig.7 Comparison of a set of pulse profiles for two weighting methods

圖8 兩組取權方法的標準脈沖輪廓比較Fig.8 Comparison of standard pulse profiles for two weighting methods

4 結 論

本文利用Python環境下的Astropy包處理了XPNAV-1先期公布的35組觀測數據，分別建立了積分脈沖輪廓和標準脈沖輪廓，采用DFT方法將二者做互相關處理獲得了測量TOA及其測量誤差σTOA。結合英國Jordrell Bank臺發布的Crab脈沖星星歷與美國JPL發布的DE200歷表，進一步得到了擬合前計時殘差。該殘差不存在系統趨勢項，從而驗證了XPNAV-1觀測數據的有效性。同時，我們分析了在目前觀測條件下，建立脈沖輪廓時子相位間隔N的取值與計時殘差的RMS之間的關系。我們發現N取256是合適的。N=256時，計時殘差的RMS為1.14 μs。同時也表明在現有觀測條件和數據量的情況下，適宜采用折疊-DFT相關法獲得測量TOA。從而證明利用中國自主觀測資料，采用合適數據處理方法，能夠獲得較高的X射線脈沖星計時精度。由于Crab脈沖星的星歷由射電波段的觀測給出，而測量的TOA來源于X射線波段的觀測，這二者存在零點相位差。在這35組的觀測時間段內，我們得到零點相位差約為0.04，對應的時間差約為1.38 ms。本文利用XPNAV-1發布的數據只討論了計時精度，旨在驗證觀測數據的有效性和數據質量。

在構建標準脈沖輪廓時，本文采用了全部觀測數據以提高信噪比。而Crab脈沖星的脈沖輪廓會隨時間推移產生變化[49]。先期發布的觀測數據時間跨度約為一個月，這種變化可以忽略。當獲得大量數據后，如果采用全部觀測數據構建標準脈沖輪廓，觀測數據的時間跨度可能會對其信噪比的提升造成影響。采用多長時間段的數據建立高信噪比的標準脈沖輪廓也是值得研究的一個問題。