楊成偉，鄭建華

（1. 北京理工大學(xué) 機(jī)電學(xué)院，北京 100081；2. 中國科學(xué)院 國家空間科學(xué)中心，北京 100190；3. 中國科學(xué)院大學(xué)，北京 100049）

0 引 言

脈沖星自主導(dǎo)航技術(shù)經(jīng)歷多年的發(fā)展正由最初的理論探索轉(zhuǎn)入工程應(yīng)用，正由概念逐步變?yōu)楝F(xiàn)實(shí)。國內(nèi)外的大量學(xué)者為脈沖星導(dǎo)航的發(fā)展付出了艱辛的努力與探索。

脈沖星之所以能夠用于導(dǎo)航是因?yàn)槠湓谖锢肀举|(zhì)上具有極為出色的時間特性。脈沖星輻射的物理信號從射電到X射線再到射線都有一定的時間特征，其中最適合于航天器導(dǎo)航的是X射線波段的物理信號。

本文就作者曾從事的關(guān)于脈沖星X射線數(shù)據(jù)處理和分析的基本技術(shù)進(jìn)行綜述，為從事X射線脈沖星導(dǎo)航和脈沖星天文研究提供參考。

1 FITS文件

航天器觀測到的X射線數(shù)據(jù)一般以FITS（Flexible Image Transport System）格式文件進(jìn)行存儲，F(xiàn)ITS文件在天文數(shù)據(jù)處理中有廣泛應(yīng)用。天文觀測中的FITS文件一般由文件頭和數(shù)據(jù)組成文件頭含有數(shù)據(jù)觀測的時間系統(tǒng)、起止時刻、觀測設(shè)備、觀測對象信息等大量內(nèi)容數(shù)據(jù)則是觀測設(shè)備記錄下的基本觀測量。在X射線觀測中往往記錄的是光子到達(dá)探測器的時刻、位置和能量等信息，這些信息在文件頭（Header）中都有詳細(xì)的描述。

Heasoft軟件包是X射線數(shù)據(jù)處理中的一款流行軟件，能夠完成多種處理功能?，F(xiàn)僅簡要介紹一些基本功能的實(shí)現(xiàn)過程。在完成具體功能前，首先要進(jìn)行Heasoft軟件包的初始化，由于主要針對NuSTAR（The Nuclear Spectroscopic Telescope Array）任務(wù)進(jìn)行數(shù)據(jù)分析，因此參考作者的使用每次運(yùn)行Heasoft軟件包前都要運(yùn)行以下語句：

其中：第1句初始化CALDB校準(zhǔn)數(shù)據(jù)庫；第2句針對使用XSELECT功能的初始化；第3句對Heasoft軟件包的初始化。

Heasoft軟件包中的FV（Fits Viewer）命令能夠通過圖像或列表的形式對FITS文件進(jìn)行查看、修改和選擇。如要讀取脈沖星PSR B1509-58在觀測ID 400240 01002下的數(shù)據(jù)時，在Linux系統(tǒng)終端輸入：

此時會打開fv界面顯示相關(guān)數(shù)據(jù)如圖1所示。從該界面可以看到該文件由3部分組成，分別是EVENTS、GTI和REG部分。EVENTS部分包含19 067個X射線光子到達(dá)時間，每個光子特有屬性共有14列的說明，除時間外還包括能量、坐標(biāo)、通道等信息；GTI（Good Time Interval）表示可用光子數(shù)據(jù)的時間范圍；REG為光子數(shù)據(jù)的位置范圍利用圖像模式的X射線探測器獲得的數(shù)據(jù)需要先根據(jù)觀測源的位置，選定數(shù)據(jù)位置范圍才能獲得較高信噪比的觀測源數(shù)據(jù)。

圖1 利用FV讀取B1509-58脈沖星的X射線觀測數(shù)據(jù)Fig. 1 Using FV to read the X-ray observation data of B1509-58 pulsar

2 X射線數(shù)據(jù)初步處理

經(jīng)過空間觀測，能夠得到X射線光子到達(dá)時間的原始數(shù)據(jù)，但由于像素壞點(diǎn)、坐標(biāo)系轉(zhuǎn)換、虛光、像素位置誤差等因素的影響必須要對數(shù)據(jù)進(jìn)行校準(zhǔn)。美國加州理工大學(xué)的科學(xué)數(shù)據(jù)運(yùn)行中心為NuSTAR任務(wù)[1]編制了用于數(shù)據(jù)校準(zhǔn)的軟件NuSTARDAS（NuSTAR Data Analysis Software），利用該軟件可以完成NuSTAR數(shù)據(jù)的處理和分析等操作。

NuSTAR任務(wù)配備了兩個相同的X射線光子焦平面收集組件（Focal Plane Module，F(xiàn)PM），分別為FPMA和FPMB，有效探測能級為3.0～79.0 keV。在NuSTARDAS軟件中重要的處理命令為nupipeline，主要完成坐標(biāo)系轉(zhuǎn)換、姿態(tài)校準(zhǔn)、壞像素剔除、能量校準(zhǔn)、天球位置計算、根據(jù)校準(zhǔn)文件庫的數(shù)據(jù)剔除、數(shù)據(jù)可用時間段計算、由時間分辨率引起的無效時間的校準(zhǔn)、曝光時間計算等處理，調(diào)用格式如下：

其中：indir為觀測數(shù)據(jù)的位置；steminputs為觀測數(shù)據(jù)的ID；outdir為經(jīng)處理后的數(shù)據(jù)存放位置。

在完成數(shù)據(jù)校準(zhǔn)后便得到了可用于科學(xué)分析的數(shù)據(jù)，由于空間背景的存在使得信噪比水平較低。因此，需要進(jìn)行脈沖星源的選取，需要去掉背景輻射光子才能進(jìn)行脈沖星的時間分析和譜分析。首先，需要用到DS9軟件選取星源區(qū)域，再用Xselect命令選取指定區(qū)域的光子序列，獲得的指定區(qū)域的X射線脈沖星光子序列是光子到達(dá)航天器的時間（Local Time）。由于航天器高速飛行產(chǎn)生的多普勒效應(yīng)會使折疊出的脈沖輪廓并不精確，為了進(jìn)行精確的時間分析必須將光子到達(dá)航天器的時間轉(zhuǎn)換到光子到達(dá)太陽系質(zhì)心的時間，并修正星載時鐘的時差。Heasoft軟件包集成了該時間轉(zhuǎn)換命令barycorr調(diào)用格式為

其中：infile為需要進(jìn)行太陽系質(zhì)心時間轉(zhuǎn)換的光子序列文件；outfile為時間轉(zhuǎn)換的輸出文件；orbitfiles為衛(wèi)星軌道文件；clockfile為星載時鐘的鐘差文件；衛(wèi)星任務(wù)組會定期發(fā)布，RA為脈沖星的赤經(jīng)；DEC為脈沖星的赤緯。

barycorr軟件中進(jìn)行的時間轉(zhuǎn)換為

其中：tb為光子到達(dá)太陽系質(zhì)心的時間即脈沖星導(dǎo)航中利用的光子到達(dá)時間；tobs為光子到達(dá)航天器的觀測時間；（clock）為將航天器的本地時間系統(tǒng)轉(zhuǎn)換為地球質(zhì)心力學(xué)時，包括星載時鐘的修正，該項與具體的X射線觀測望遠(yuǎn)鏡有關(guān)。

在RXTE任務(wù)中，只要進(jìn)行適當(dāng)?shù)臅r鐘修正，星載記錄時間就可以直接表達(dá)為地球質(zhì)心力學(xué)時。（dispersion）為色散改正，對于X射線來講色，散改正量為0，（geometric）為幾何學(xué)時間延遲量，（"Einstein"）為由相對論效應(yīng)引起的愛因斯坦延遲，（"Shapiro"）為沙皮諾延遲是由太陽系內(nèi)引力場的光線彎曲引起的。

經(jīng)上述一系列的數(shù)據(jù)處理后就可以進(jìn)行時間分析和譜分析。時間分析主要進(jìn)行周期檢測、周期的一階和二階導(dǎo)數(shù)擬合、光子計數(shù)及計數(shù)率統(tǒng)計、基于歷元折疊的輪廓分析、周期躍變、雙星系統(tǒng)軌道參數(shù)分析等。工作譜分析主要進(jìn)行星源能量譜的模型擬合、背景譜分析、物理輻射模型建立與驗(yàn)證等工作，往往需要結(jié)合具體的X射線望遠(yuǎn)鏡的相關(guān)參數(shù)和文件進(jìn)行分析。

3 X射線時間分析

3.1 周期檢測

X射線時間分析的重要內(nèi)容是進(jìn)行星源的周期檢測，包括脈沖星自身的旋轉(zhuǎn)周期檢測和雙星系統(tǒng)中軌道周期的檢測。這能夠幫助我們了解脈沖星的年齡等物理參數(shù)，進(jìn)而了解演化進(jìn)程。在X射線脈沖星導(dǎo)航數(shù)據(jù)庫的構(gòu)建中，周期檢測是基礎(chǔ)，只有首先確定脈沖星星源的信號周期性以及輪廓特征才能得到導(dǎo)航脈沖星數(shù)據(jù)庫的精確參數(shù)。

離散傅里葉周期檢測適用于在周期未知情況下的檢測，但計算量較大。因此，在已知脈沖星周期范圍時往往采用周期檢測或周期檢測，這2種方法都在一定的周期范圍內(nèi)進(jìn)行搜索極大值點(diǎn)，即與脈沖星周期相比離散傅里葉變換方法精度更高。周期檢測方法較周期檢測方法計算量大但有著更高的搜尋精度。

在周期范圍了解的情況下，可以采用X射線天文學(xué)中廣泛使用的檢測標(biāo)準(zhǔn)（卡方檢測）進(jìn)行周期檢測檢驗(yàn)。按照某一周期折疊的脈沖輪廓是否具有足夠好的分布，即：脈沖輪廓各個Bin中的光子數(shù)目距平均值的離散程度。在指定范圍內(nèi)獲得離散程度最高的輪廓，其對應(yīng)的周期值即為脈沖星的自轉(zhuǎn)周期。

用于脈沖星輪廓檢測的卡方計算標(biāo)準(zhǔn)定義為

圖2 卡方周期檢測Fig. 2 Period detection using Chi-square test

雙星系統(tǒng)的軌道周期較長，每個Bin的時間間隔也相對較長。由于光子探測無效時間的存在，如果簡單地利用各個Bin中的光子數(shù)目進(jìn)行檢測往往會造成很大的誤差。因此，在雙星系統(tǒng)軌道周期檢測中采用光子計數(shù)率（單位時間內(nèi)的光子數(shù)目）來進(jìn)行卡方分布統(tǒng)計。定義式為

在多個搜索周期下，進(jìn)行折疊將得到一系列的統(tǒng)計量S。由于光子到達(dá)時間為隨機(jī)量，在背景噪聲的影響下，卡方檢測的最大值不一定是最優(yōu)的結(jié)果。但大量S值能夠表現(xiàn)出一定的統(tǒng)計特性。利用含有周期項的統(tǒng)計函數(shù)進(jìn)行擬合，將使結(jié)果更具有可信性。Leahy提出的Epoch Folding方法[2]能夠?qū)ǚ椒植嫉腟量進(jìn)行擬合，使結(jié)果更為精確。在光子較少的情況下Epoch Folding方法同樣具有很高的精度，Epoch Folding方法的擬合函數(shù)為

利用Epoch Folding方法進(jìn)行4U 0142+61脈沖星的周期搜索，該脈沖星位于仙后座中，是非常規(guī)X射線脈沖星的一員，距地球約1.3萬光年、自轉(zhuǎn)周期約為8.69 s、周期一階導(dǎo)數(shù)為0.2 × 10–11s s–1、兩級磁場強(qiáng)度為

在2～10 keV能段，在非常規(guī)X射線脈沖星中具有較高的光度，約為1 × 1035erg/s。從紅外波段到硬X射線波段均有輻射。在觀測中也發(fā)現(xiàn)了該脈沖星的脈沖輪廓有緩慢的變化，并且存在X射線暴現(xiàn)象。NuSTAR任務(wù)在2014年3月對該脈沖星進(jìn)行了觀測。Epoch Folding方法在進(jìn)行周期擬合后得到自轉(zhuǎn)周期為P= 8.689 157 5 ± 3.799 ×10－6s，同其他文獻(xiàn)中得到的自轉(zhuǎn)周期一致，如圖3所示。

圖3 0142+61自轉(zhuǎn)周期檢測Fig. 3 0142+61 rotation cycle period detection

根據(jù)脈沖星時間相位模型每個光子到達(dá)時間都對應(yīng)一個脈沖相位即

其中：N為光子總數(shù)；n為諧波分量的個數(shù)。

3.2 脈沖比例

在完成周期檢測后，需要進(jìn)行脈沖比例（Pulsed Fraction）的分析，用以評價折疊出的脈沖輪廓的強(qiáng)度及好壞。該分析結(jié)果在X射線時間分析中有重要作用。比較不同輪廓間的差別來分析脈沖輻射相對于背景輻射的強(qiáng)度。若要在脈沖星導(dǎo)航數(shù)據(jù)庫建立中起到參考作用更傾向于選擇脈沖比例大的脈沖星。脈沖比例的計算方法一般有以下幾種。

圖4 方法進(jìn)行周期檢測Fig. 4 Period detection using test

1）峰值脈沖比例（Peak-to-Peak Pulsed Fraction）

峰值脈沖比例以脈沖輪廓中的最小值和最大值進(jìn)行脈沖輪廓的評估。

圖5 正弦波形脈沖輪廓Fig. 5 Sinusoidal pulse profile

則峰值脈沖比例為

2）均方根脈沖比例（Root Mean Square Pulsed Fraction）

均方根脈沖比例，基于傅里葉變換通過傅里葉能譜評估脈沖輪廓的強(qiáng)度。由于該評價方法利用了輪廓中所有的bin值，相比峰值脈沖比例，有著更好的統(tǒng)計性能。

其中：N為脈沖輪廓中相位間隔的數(shù)目；pi為每個相位間隔中的光子數(shù)；k為傅里葉變換的諧波數(shù)目；和分別為ak和bk的不確定度。

3）面積脈沖比例（Area Pulsed Fraction）

從表達(dá)式上很容易理解面積脈沖比例計算的是脈沖振幅的面積與整個脈沖輪廓面積的比值反映了脈沖輪廓相對于噪聲的強(qiáng)度。

在不同能級下，同一顆脈沖星的脈沖輪廓存在差異，這時便可以利用脈沖比例進(jìn)行脈沖強(qiáng)度的分析。例如在不同能級下，4U 0142+61具有不同的累積脈沖輪廓，按照3～5、5～8、8～20、20～35、35～50、50～79 keV劃分能量段可得到不同能級下的脈沖輪廓，如圖6所示，從中可以看出在不同能級上脈沖輪廓的形態(tài)有很大的差異。

按照上述的能段劃分利用Area方法進(jìn)行脈沖比例的計算，可以得到不同能段的脈沖比例。Area脈沖比例分別為：12.7% ± 1.2%（3～5 keV）、15.4% ±2%（5～8 keV）、14.9% ± 3%（8～20 keV）、30.8% ±7.2%（20～35 keV）、33.9% ± 12.6%（35～50 keV）以及52.7% ± 20%（50～79 keV）?？芍狝rea脈沖比例隨著能級的增加持續(xù)增長，這種趨勢同Hartog[4]得出的結(jié)論一致。

4 X射線譜分析

XSPEC軟件是美國航空航天局（National Aeronautics and Space Administration，NASA）管理和發(fā)布的關(guān)于X射線譜分析的軟件，能夠處理多個X射線探測任務(wù)的數(shù)據(jù)、完成數(shù)種輻射模型的建立和擬合，并對擬合結(jié)果進(jìn)行統(tǒng)計評估。通過圖形顯示和數(shù)值分析，幫助人們推斷天體輻射的物理特性。XSPEC也能夠模擬特定X射線探測任務(wù)的能譜，在科學(xué)任務(wù)分析和規(guī)劃中起到重要作用，在高能X射線譜分析中應(yīng)用十分普遍。在脈沖星導(dǎo)航中，可以用于之前X射線探測任務(wù)的結(jié)果，進(jìn)行對比分析以確認(rèn)脈沖星導(dǎo)航探測器是否存在問題，也可以在導(dǎo)航的同時兼顧科學(xué)分析。

在仿真中需要用到兩個重要文件，即RMF文件和ARF文件，描述X射線探測設(shè)備的性能。RMF文件是響應(yīng)矩陣文件（Response Matrix File），描述了X射線探測設(shè)備對X射線能量及收集通道的響應(yīng)。ARF是輔助響應(yīng)文件（Ancillary Response File），包含了不同能級。X射線光子的有效面積響應(yīng)及其參數(shù)影響因素，主要有X射線探測設(shè)備的設(shè)置、視場位置、有效面積、量子效率、濾波器轉(zhuǎn)換效能、暈影改正等因素。

在獲得兩個文件之后便可以建立輻射模型。XSPEC能根據(jù)輻射模型和設(shè)備的響應(yīng)參數(shù)進(jìn)行隨機(jī)數(shù)打靶，獲得較高置信度的模擬X射線能譜。對該能譜進(jìn)行擬合，便能給出該探測設(shè)備在指定的觀測時間內(nèi)，所能獲得的探測目標(biāo)物理參數(shù)的精度。其中，幾個比較重要的命令為model命令、fakeit命令和fit命令。

model命令可以根據(jù)設(shè)定的參數(shù)，進(jìn)行光子能量譜模型的建立，用于分析星源的輻射模型。常用的輻射模型為黑體輻射譜模型（blackbody）和冪律譜模型（powerlaw），黑體譜可由黑體輻射引起，冪律譜可由熱軔致輻射、同步加速輻射和逆康普頓散射引起。

XSPEC中的黑體輻射譜模型為

其中：k為玻爾茲曼常數(shù)；T為黑體絕對溫度，單位為keV；為星源光度，單位為1039erg/s，為到星源的距離，單位為10 kpc；E為輻射能量。

圖6 4U 0142+61不同能級的脈沖輪廓Fig. 6 4U 0142+61 pulse profiles with different energy level

XSPEC中的冪律譜模型為

fakeit命令能夠根據(jù)X射線探測器的RMF和ARF文件，仿真出該X射線探測器能捕捉到的X射線光子能量數(shù)據(jù)，這可以作為模擬數(shù)據(jù)進(jìn)行模型的擬合，具有很高的可信度，用于觀測任務(wù)論證等場合。調(diào)用格式為

其中，nicer.rmf和nicer.rmf為NICER探測器對應(yīng)的RMF和ARF文件，J1847.fak為仿真出的數(shù)據(jù)文件，200 000為仿真的觀測時長，單位為s、y和blah，為格式輸入。

fit命令根據(jù)建立的輻射模型，利用觀測或模擬生成的數(shù)據(jù)進(jìn)行模型擬合，并給出多種統(tǒng)計量以評判擬合結(jié)果的可信性，當(dāng)可信性較高時，便可以接受模型，推斷出星源的物理輻射形式。

仿真中用到的RMF文件和ARF文件，由美國宇航局戈達(dá)德航天中心提供。XSPEC中模型的參數(shù)需要參考相關(guān)探測目標(biāo)的研究文獻(xiàn)。但有些探測目標(biāo)由于沒有進(jìn)行觀測或觀測時間較短，或X射線望遠(yuǎn)鏡的觀測能力有限，無法從文獻(xiàn)中查找，此時可以簡單地通過黑體輻射熱力學(xué)和普朗克函數(shù)進(jìn)行估計。

普朗克黑體輻射定律為

如果已知星源的半徑、溫度和距離，便能通過普朗克黑體輻射定律的積分求出指定能量段的輻射流量，進(jìn)而得到黑體輻射模型的系數(shù)K。例如，磁星J1847-0310的半徑約為12 km、溫度約為0.13 keV、距離為8.4 kpc，那么對普朗克輻射定律求積分后，可得到在0.2～10 keV能量段的總輻射能量為7.597 91029erg/s，X射線探測器單位面積上接收的光子流量為9.032 810–17erg/s/cm2。

5 結(jié) 論

本文對脈沖星X射線的基本數(shù)據(jù)處理流程和基本分析方法進(jìn)行了綜述，使讀者能夠快速進(jìn)入脈沖星X射線數(shù)據(jù)處理的大門，為脈沖星導(dǎo)航和脈沖星天文研究提供了參考。