利用Crab脈沖星X射線觀測校準星載原子鐘頻率

童明雷，韓孟納，楊廷高，趙成仕，朱幸芝

1．中國科學院 國家授時中心，西安 710600 2．中國科學院大學，北京 100049

脈沖星是大質量恒星超新星爆發或雙星吸積塌縮后形成的帶強磁場的中子星，其旋轉產生的輻射束可被地面射電望遠鏡或空間高能探測器接收。因其超高溫、超高壓、超高密度、超強磁場、超強電場、超強引力場等極端物理條件，脈沖星成為天文和物理研究的天然實驗室。脈沖星是20世紀60年代的四大天文發現之一。脈沖星的觀測和理論研究自發現以來已持續了50多年，是天體物理領域重要的研究對象。

國內外對脈沖星的射電觀測已長達數十年之久，隨著各種地面大型射電望遠鏡的陸續建成，天文學家對脈沖星的形成、輻射機制、內部結構及旋轉特征等進行了深入研究［1］。值得一提的是，2020年12月1日，因饋源平臺墜落導致受損嚴重的阿雷西博射電望遠鏡（Arecibo）在其服役期間對脈沖星的觀測取得了眾多劃時代成果［2］：1974年Hulst和Taylor發現了第1個脈沖星-中子星雙星系統（PSR B1913+16），間接證明了引力波的存在［3］；1982年發現了第1顆毫秒脈沖星（PSR B1937+21）［4］；1992年Wolszczan和Frail發現第1顆伴星為行星的脈沖星（PSR B1257+12）［5］，這是人類首次發現系外行星。隨著Are?cibo的落幕，被譽為“中國天眼”的500 m口徑球面射電望遠鏡（FAST）于2021年3月31日對全球開放，接受來自全世界的觀測申請。國際上，中國參與了平方公里陣（SKA）項目，并將脈沖星搜尋、脈沖星測時及引力檢驗作為重點研究方向之一。

深空探測任務具有距離遠、延遲長、信號弱等特點，在此條件下高精度導航始終是深空測控技術需要重點解決的問題。X射線脈沖星導航作為一種新的天文導航方式，通過觀測脈沖星輻射的X射線光子，便可獲得高精度的測距信息與時間信息［6］。脈沖星輻射的高能X射線集中了絕大部分的輻射能量，易于被小型化探測器接收，有利于減小航天器有效載荷的尺寸。同時，脈沖星導航彌補了全球衛星導航系統（GNSS）導航精度隨測控距離的增加而降低的缺陷，是一種真正意義上的自主導航，適用于廣闊的太陽系空間的深空探測甚至星際飛行任務的航天器自主導航，近年來逐漸成為國內外研究的熱點。

脈沖星導航雖然未來應用前景廣泛，但由于X射線探測器的靈敏度較低，探測流量極低的毫秒脈沖星還比較困難，因此目前并沒有實質的應用。而利用脈沖星進行空間自主守時有望更早地獲得應用。由于脈沖星時的長期穩定性好，基于射電波段的長期脈沖星計時觀測資料可以建立脈沖星時［7-8］。有些觀測精度較高的脈沖星，利用其2～3年計時數據就可以建立穩定脈沖星時［9］。對于既有射電輻射又有X射線輻射的脈沖星，在利用射電波段構建脈沖星計時模型參數后可以作為數據庫，聯合X射線波段的計時觀測用于空間守時系統。由于脈沖星的自轉頻率變化率通常很小，因此短時間內脈沖星的自轉頻率非常穩定，構成一種頻率基準，基于此可以校準原子鐘的頻率，提高其準確度。本文將結合Crab脈沖星射電波段數據庫和X射線計時觀測進行空間自主守時研究，利用脈沖星改正空間原子鐘的頻率偏差而不依賴于地面。

脈沖星導航試驗衛星（XPNAV-1）是由中國航天科技集團第五研究院研制的中國首顆X射線脈沖星導航試驗衛星，于北京時間2016年11月10日在酒泉衛星發射中心由長征11號運載火箭發射升空［10］。衛星軌道高度500 km，軌道周期94 min。因在執行觀測任務時為降低對探測器的損害，衛星在觀測時避開了南大西洋異常區，每個觀測弧段的持續時間通常只有30～50 min［11-12］。衛星采用整星零動量三軸穩定姿態控制方式，搭載了掠入射Wolter-I聚焦型探測器和準直型微通道板探測器［13］。Crab脈沖星（PSR B0531+21）作為天空中最明亮的X射線源之一，成為XPNAV-1衛星進行在軌標定的標準源。本文主要對XPNAV-1衛星的Crab脈沖星觀測數據進行處理，分析研究了星載原子鐘的頻率偏差對脈沖到達時間的影響，據此進一步給出可對頻率偏差進行修正的方法。

1 X射線脈沖星計時

與射電脈沖星計時不同，X射線探測器記錄的是X射線光子的到達時刻。XPNAV-1衛星的公開發布數據［14］包括光子事件文件和軌道文件2類數據。由于光子到達探測器的時刻與衛星軌道遙測時刻并不一致，因此需要在光子到達時刻內插出航天器的位置與速度。衛星在繞地飛行時，其所處的引力場是在不斷變化的，“多普勒效應”“引力紅移”以及視差等效應會使探測器探測到的光子到達時刻失去周期性，需將光子到達探測器的時刻轉換為到達太陽系質心（SSB）的時刻。到達時間的轉換模型包括：衛星繞地運動和地球公轉引起的幾何延遲、引力場中光線彎曲引起的Shap?iro延遲、引力場中時間尺度變換的Einstein延遲等，具體過程可見文獻［15-16］中的描述。

1.1 脈沖輪廓折疊

要得到脈沖到達SSB的時刻（TOA），需要將光子到達SSB的時間序列先進行歷元折疊得到積分脈沖輪廓與標準脈沖輪廓，再由二者的相關獲得脈沖到達時間。歷元折疊方法是將脈沖星的自轉周期均勻劃分為N等份（bin數），統計落入每個bin的光子數，從而獲得脈沖輪廓。折疊過程使用Jodrell Bank天文臺射電觀測的Crab脈沖星自轉模型參數［17］，使用的2段Crab脈沖星自轉參數如表1所示。在公開發布數據時間跨度內，Crab脈沖星自轉參數更新過一次。表1中P1代表第1段自轉參數，P2代表第2段自轉參數。當光子到達時刻（以約化儒略日表示）<57727時，使用P1，否則使用P2。

表1 Crab脈沖星自轉參數Table 1 Spin parameters of Crab pulsar

在進行脈沖輪廓折疊時，若子相位間隔數N取值太大，則折疊的脈沖輪廓會損失信噪比；若N取值太小，脈沖輪廓將過于平滑，丟失頻域信息。圖1給出了采用第24組數據、按不同N（分別取64、128、256、512、1024、4096）折疊的積分脈沖輪廓，顯然脈沖輪廓的信噪比明顯地依賴于N值的選取。表2給出了各脈沖輪廓的信噪比（SNR），計算公式［18］為

圖1 積分脈沖輪廓隨子相位間隔數的變化Fig. 1 Variations of integrated pulse profile with number of sub-phase intervals

式中：Npeak為脈沖輪廓主峰的光子數；Ntotal為觀測的總光子數；p為曝光時間。表2中p均為44735 s。此外，N的取值也與探測器的靈敏度以及每個觀測文件的數據量有關。脈沖星輪廓的SNR只是作計時數據處理的一個參考，不能作為唯一標準。盡管表2中N=64的信噪比最高，但其脈沖星輪廓的形狀可能已經失真。后面的分析發現N?。?4，512］的范圍內，計時結果很穩定。因此如無特別說明，不失一般性地均取N=256［19］。

表2 按不同子相位間隔數折疊的積分脈沖輪廓SNRTable 2 Signal-to-noise ratio of integrated pulse profile folded by different number of sub-phase intervals

Jodrell Bank天文臺發布的Crab脈沖星星歷參數是基于美國JPL太陽系行星歷表DE200得到的，而JPL歷表參考的時間尺度是太陽系質心力學時（TDB），因此表1中t0給出的是TDB時刻。折疊的標準脈沖輪廓如圖2所示。

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

1.2 脈沖輪廓平滑

當脈沖輪廓折疊所用的光子數過少時，得到的脈沖輪廓信噪比會很低，從而影響脈沖到達時間的確定，通常需將輪廓進行平滑降噪處理。核回歸作為一種非參數估計方法，其對數據的分布特征不附加任何假定，是一種從數據本身出發研究數據分布特征的方法［20］。定義x處的觀測量為

式中：yi為x處的已知觀測量；wi(x)為x處的權函數。權函數的定義為

其中：h為核寬度，h的大小決定了平滑程度；κ(m)為核函數。本文比較了3種常用的核平滑方法（Gaussian核平滑、Epanech?nikov核平滑以及Tri-cube核平滑）對脈沖輪廓的平滑效果。3種核函數的定義為

式中，κ1(m)、κ2(m)、κ3(m)分別為Gaussian核函數、Epanechnikov核函數以及Tri-cube核函數。

平滑的具體方法是選取其中一組原始積分脈沖輪廓，采用上述3種方法進行降噪處理，處理后的結果如圖3所示?？梢园l現，當使用相同的核寬度時，相較于其他2種核回歸法，Gaussian核回歸可以有效去除脈沖輪廓中的隨機噪聲并保留信號的峰值特征，平滑效果最好。且由于高斯核函數靈活性比較高，當特征數較小時，一般采用高斯核回歸平滑含噪聲的信號。因此下文均采用Gaussian核回歸平滑法對脈沖輪廓降噪處理。

圖3 脈沖輪廓經3種核回歸處理前后的結果比較Fig. 3 Results of pulse profile before and after being processed by three kernel regressions

1.3 計時殘差獲取

計時殘差是衡量脈沖星計時水平的一個重要因素，要得到計時殘差需要首先獲得TOA數據。本文的TOA是采用積分脈沖輪廓與標準脈沖輪廓經DFT之后在頻域互相關［21］的方法得到的，頻域互相關相較于時域互相關可以獲得更高的TOA測量精度。計時殘差為測量的TOA與模型預報的TOA之差［22］，即

式中：R(t)表示計時殘差；?(t)為脈沖到達時刻t對應的相位，N(t)為與?(t)最接近的整數，ν為所觀測脈沖星的自轉頻率。

將XPNAV-1衛星的35組觀測數據，結合Jodrell Bank發布的Crab脈沖星的射電計時模型參數，并擬合掉2個波段的零點相位差后得到的擬合前計時殘差如圖4所示，圖中脈沖到達時刻用約化儒略日（MJD）表示。零點相位差不存在長期變化趨勢，在短時間內變化很小，可近似視為常數［23-24］。圖4的結果利用了全部數據構建的標準脈沖星輪廓，得到的計時殘差均方根（RMS）值為29.4853 μs。在這35組觀測數據的時間跨度內，Crab脈沖星更新了一次星歷參數，可以嘗試使用不同星歷的觀測數據段分別構建各自的標準脈沖輪廓，本文暫不討論這種處理方法。圖4中并沒有出現線性項，說明Crab脈沖星射電計時模型參數對于一個月的X射線觀測數據比較準確，計時殘差主要由觀測白噪聲引起。當然，可以對計時殘差進行二次多項式擬合，得到更新后更加準確一點的自轉模型參數［16］。

圖4 擬合前計時殘差Fig. 4 Pre-fit timing residuals

2 X射線脈沖星校準星載原子鐘

2.1 參考鐘對脈沖星計時的影響

對于X射線脈沖星觀測而言，星載原子鐘的穩定度與準確度會影響光子到達探測器時刻的測量精度。在理想情況下，原子鐘頻率源的輸出信號為頻率恒定的正弦波。但由于系統噪聲和隨機噪聲的存在，頻率源的實際輸出頻率會偏離其標稱頻率，瞬時相對頻率偏差y(t)的定義為

式中：f(t)為頻率源的輸出頻率；f0為標稱頻率；x(t)為原子鐘實際輸出信號相對于理想時間信號的時間偏差，可以表示為［25］

式中：x0為初始時刻時間偏差；y0為初始時刻相對頻率偏差；a為原子鐘頻率漂移；xr(t)為隨機噪聲項。由于頻率源輸出頻率的不穩定以及航天器所處空間磁場環境的變化，星載原子鐘的輸出頻率與標稱頻率之間會存在一個頻率偏差，該偏差值會影響航天器本地時間的準確度。

XPNAV-1衛星直接接收GPS衛星發播的系統時間。GPS時是GPS信號的時間基準，由地面監控站和衛星上的原子鐘經加權算法得到［26］。GPS衛星軌道高度約為20200 km，星載鐘主要受廣義相對論的“引力紅移”效應，若衛星發射之前星載鐘未作頻率調整，GPS衛星鐘的走時將比地球時（TT）快。為使GPS鐘的走時速率與TT一致，GPS衛星發射前將對原子鐘進行頻率調整，扣除相對論效應，不再以SI秒為單位測量原時。這種方法的優點是無需考慮地球的引力場及衛星相對于地心的速度對星載原子鐘走時速度的影響。但這種時間接收方法只適用于地球表面的近地空間，對于更廣闊的太陽系以及星際空間，由于距離問題，航天器無法接收到GPS信號，只能通過攜帶的星載鐘來計時。對于XPNAV-1衛星，若用星載原子鐘記錄光子到達時刻，則需要考慮廣義相對論的時空線元和度規，將原時轉化為坐標時，并最終轉為TDB。對于地球軌道衛星，原時到坐標時轉換公式為［27］

式中：U為地球在衛星處產生的引力勢大小，衛星所處的引力勢與衛星的位置有關；v為衛星相對于地心的速度；c為光速。TT與地心坐標時（TCG）的關系為［28-30］

式中：LG≡6.969290134×10?10為一常數。若衛星飛行軌道近似為圓軌道，則衛星所處的引力場可視為常數。由式（7）與式（8）可得原時與地球時的關系：

式中：G為引力常量；M為地球質量；r為地球半徑；H為衛星的軌道高度。若XPNAV-1衛星使用自己攜帶的原子鐘，經計算，星載原子鐘相對于TT走時的變化率為?2.6941×10?10。

對于各種波段的脈沖星計時觀測，TOA是由參考鐘記錄的，進而通過時間轉換溯源到更高精度的時間尺度。脈沖星的計時殘差方程［2］：

原子鐘的鐘差、頻率偏差、頻率漂移均會使記錄的光子到達時刻不準確，最終分別導致擬合前計時殘差存在常數偏差、線性項以及二次項。式中：R為t時刻的計時殘差；R0為t0時刻的殘差；?v0和?v?分別為v0、v?的修正量；?α、?δ分別為對脈沖星赤經、赤緯的修正；μα、μδ分別為赤經方向上的自行與赤緯方向上的自行；A、B為脈沖星位置修正項的系數。

為研究星載時鐘的頻率偏差對X射線脈沖到達時間的影響，本文仿真了參考時鐘存在頻率偏差時脈沖星光子到達探測器的時刻，這里鐘的相對頻率偏差取為y0=2.6941×10?10，即衛星星載鐘原時與TT速率的相對偏差。此時，第i個光子的到達時間變為

式中：ti表示參考時沒有頻率偏差時第i個光子的達到時間。

按照X射線脈沖星計時流程重新做處理，得到擬合前計時殘差，注意，這里依然用了參考時沒有偏差時構建的標準脈沖輪廓，這是為了避免參考鐘偏差引起的標準脈沖輪廓變形和信噪比的降低，而這在實際情況中是比較容易做到的，只需利用Crab脈沖星時間校準過的歷史數據構建標準脈沖輪廓即可。同時，為了避免復雜性，后面的討論和計算不考慮分段構建標準脈沖輪廓的情況。對計時殘差進行線性擬合，擬合結果如圖5所示，擬合斜率為2.6858×10?10，擬合斜率的不確定度為4.2842×10?14。從擬合結果可以看出，通過計時殘差的擬合斜率值可估計原子鐘的相對頻率偏差。因此，可通過X射線脈沖星觀測自主校準星載鐘的頻率偏差，在一定程度上保證星載原子鐘的準確度。而通過相對密集的頻率駕馭，還可以提高星載鐘的穩定度。需要指出的是，本文是基于真實觀測數據仿真了記錄光子到達探測器時刻的原子鐘存在的頻率偏差，這只改變了光子到達時刻的時間溯源，并沒有仿真脈沖星信號的本身，因此完全保留了實測數據的噪聲特征。

圖5 星載鐘存在頻率偏差時的計時殘差及線性擬合結果Fig. 5 Timing residuals and linear fitting results of spaceborne clock with frequency deviation

2.2 不同N值確定的相對頻率偏差

考慮到不同bin數對脈沖星計時殘差的影響較大，計算了當參考鐘存在頻率偏差，N取不同值時相對頻率偏差估計值yc的變化情況，結果如表3所示，同時給出了相對頻率偏差的相對誤差δyc：δyc=（yc?y0）/y0。由表3可知，N取值不同，頻率偏差估計值與實際設定值之間的偏離程度不同。若N取的太大，則脈沖輪廓信噪比太低，會影響脈沖到達時間的測量精度；若N取的太小，則脈沖輪廓會過于平滑，將間接掩蓋頻率偏差引起的計時殘差的系統性趨勢。當N取的過大時，隨著N的增加，頻率偏差估計值與實際設定值之間的相對誤差也隨之變大。因此，對于鐘差漂移的修正，N要取合適的值。也同時給出了脈沖輪廓經高斯平滑后得到的相對頻率偏差估計值ys及其與實際設定值的相對誤差δys。比較δyc和δys可知，脈沖輪廓經核回歸處理后可有效改正脈沖輪廓低信噪比對頻率偏差估計精度的影響。為了更加直觀地反映相對頻率偏差估計值隨N值的變化，圖6給出了相對頻率偏差估計值與實際設定值的差值（即絕對誤差）隨N的變化趨勢。如圖6所示，當N>512時，yc誤差的絕對值隨N增大而增大；而ys基本不隨N值的變化而變化，即對N的取值不敏感。因此平滑脈沖輪廓有利于減小結果對N值的依賴性。由表3給出的相對頻率偏差估計值與實際設定值之間的相對誤差可知，如果不平滑脈沖輪廓，則N最好≤512。

表3 不同子相位間隔數對應的頻率偏差估計值及其相對誤差Table 3 Estimates of frequency deviations correspond?ing to different sub-phase intervals and their relative errors

圖6 相對頻率偏差估計值的絕對誤差隨N的變化Fig. 6 Variations of absolute error of estimated relative frequency deviation along with values of N

2.3 駕馭不同水平的星載原子鐘

基于脈沖星觀測估計的星載鐘頻率偏差的不確定度依賴于脈沖星的TOA測量精度。Crab脈沖星屬于正常年輕的脈沖星，TOA測量誤差比較大，且其存在周期躍變現象，自轉穩定性遠低于毫秒脈沖星。遺憾的是XPNAV-1測量不到毫秒脈沖星的周期信號。下面，基于Crab脈沖星的觀測，討論一下其駕馭星載原子鐘頻率偏差的性能。假設星載鐘存在10?10、10?11和10?12這3個不同量級的相對頻率偏差，N=256，并利用高斯核回歸平滑脈沖輪廓，所得結果如表4所示。由表4可知，隨著頻率偏差量級的減小，該方法得到的頻率偏差估計值與真值的偏離程度變大。這是由于Crab脈沖星TOA測量精度不高，導致頻率偏差較小時對其估計精度偏低。在實際應用過程中，可通過觀測自轉更加穩定的毫秒脈沖星，獲得更高精度的星載鐘頻率偏差值。另一方面，選用的Crab脈沖星的數據長度只有1個月左右，如果數據時間跨度更長，則其檢測星載鐘相對頻率偏差的性能會更高。當然，實際情況中星載鐘還可能存在頻率漂移，那就需要做二次多項式擬合，并通過擬合的多項式系數改正星載原子鐘的鐘差、頻率偏差及漂移項。由于所用數據比較短，在模擬鐘差時只考慮了線性的頻率偏差。

表4 相對頻率偏差取值不同對應的改正結果Table 4 Correction results corresponding to different values of relative frequency deviation

在得到星載鐘相對頻率偏差的估計值之后，據此可以駕馭星載鐘的頻率，改正其頻率偏差，使之輸出更準確的時間。具體過程是：根據脈沖星擬合前計時殘差的擬合斜率，即星載鐘頻率偏差估計值ys，將星載鐘的相對頻率偏差反向補償一個因子(1+ys)；駕馭后星載鐘的等效輸出頻率為f′=f(1+ys)，其中f是未經駕馭的星載鐘的輸出頻率。這樣就在一定程度上校準了星載鐘的頻率偏差，提高了其準確度。圖7給出了星載鐘無頻率偏差時擬合前計時殘差和星載鐘有頻率偏差又被脈沖星駕馭頻率修正后的計時殘差的比較。這個結果用了N=256時的相對頻率偏差估計值2.6858×10?10。從圖7可以看出，計時殘差不含明顯的趨勢項，與參考鐘沒有頻率偏差時得到的計時殘差相比，二者的變化趨勢基本一致。

圖7 參考鐘經頻率偏差修正后的計時殘差Fig. 7 Timing residuals after correction of reference clock frequency deviation

3 結 論

針對XPNAV-1衛星公開的數據，分析了利用Crab脈沖星駕馭空間原子鐘頻率的可行性。通過比較3種核回歸方法對脈沖輪廓的平滑能力，發現Gaussian核回歸可以更加有效地降低脈沖輪廓中的隨機噪聲并保留信號的峰值特征，提高了計時精度。基于此，仿真了星載原子鐘存在頻率偏差時的光子到達時刻數據，研究了鐘的頻率偏差對TOA及計時殘差的影響。結果表明鐘的頻率偏差會使計時殘差產生線性變化趨勢，對其線性擬合獲得了參考原子鐘的相對頻率偏差。我們計算了不同bin數的結果，發現脈沖輪廓經過高斯核回歸法平滑后，得到的參考鐘相對頻率偏差擬合值基本與bin數無關，這減小了因N的取值不同帶來的不確定性，即提高了脈沖星頻率駕馭的穩定性。通過對參考原子鐘的頻率駕馭，可校準其頻率偏差，提高其準確度。通過分別仿真星載鐘10?10、10?11和10?12水平的相對頻率偏差，一月數據長度的Crab脈沖星校準星載鐘頻率的相對誤差分別為0.3%、42%和113%。

下一步，如果能夠利用更長時間跨度的Crab脈沖星計時數據，有利于降低或消除脈沖星本身的噪聲或其他具有周期性的影響因素，進而提高相對頻率偏差的估計精度。目前星載原子鐘頻率準確度已達到10?12量級甚至更高，為此必須提高脈沖星檢驗星載鐘頻率偏差的水平，這需要觀測自轉更加穩定、計時精度更高的X射線毫秒脈沖星，預期可以大大提高星載鐘的頻率校準精度。當然，這需要靈敏度更高的X射線探測器。

通過密集的原子鐘頻率駕馭，還可以提高其長期穩定度。因此，本文的研究有利于提高空間時間系統的長期自主保持能力。