脉冲星X射线数据处理与分析

2018-11-07杨成伟郑建华

深空探测学报 2018年3期

关键词：脉冲星光子轮廓

杨成伟，郑建华

（1. 北京理工大学机电学院，北京 100081；2. 中国科学院国家空间科学中心，北京 100190；3. 中国科学院大学，北京 100049）

0 引言

脉冲星自主导航技术经历多年的发展正由最初的理论探索转入工程应用，正由概念逐步变为现实。国内外的大量学者为脉冲星导航的发展付出了艰辛的努力与探索。

脉冲星之所以能够用于导航是因为其在物理本质上具有极为出色的时间特性。脉冲星辐射的物理信号从射电到X射线再到射线都有一定的时间特征，其中最适合于航天器导航的是X射线波段的物理信号。

本文就作者曾从事的关于脉冲星X射线数据处理和分析的基本技术进行综述，为从事X射线脉冲星导航和脉冲星天文研究提供参考。

1 FITS文件

航天器观测到的X射线数据一般以FITS（Flexible Image Transport System）格式文件进行存储，FITS文件在天文数据处理中有广泛应用。天文观测中的FITS文件一般由文件头和数据组成文件头含有数据观测的时间系统、起止时刻、观测设备、观测对象信息等大量内容数据则是观测设备记录下的基本观测量。在X射线观测中往往记录的是光子到达探测器的时刻、位置和能量等信息，这些信息在文件头（Header）中都有详细的描述。

Heasoft软件包是X射线数据处理中的一款流行软件，能够完成多种处理功能。现仅简要介绍一些基本功能的实现过程。在完成具体功能前，首先要进行Heasoft软件包的初始化，由于主要针对NuSTAR（The Nuclear Spectroscopic Telescope Array）任务进行数据分析，因此参考作者的使用每次运行Heasoft软件包前都要运行以下语句：

其中：第1句初始化CALDB校准数据库；第2句针对使用XSELECT功能的初始化；第3句对Heasoft软件包的初始化。

Heasoft软件包中的FV（Fits Viewer）命令能够通过图像或列表的形式对FITS文件进行查看、修改和选择。如要读取脉冲星PSR B1509-58在观测ID 400240 01002下的数据时，在Linux系统终端输入：

此时会打开fv界面显示相关数据如图1所示。从该界面可以看到该文件由3部分组成，分别是EVENTS、GTI和REG部分。EVENTS部分包含19 067个X射线光子到达时间，每个光子特有属性共有14列的说明，除时间外还包括能量、坐标、通道等信息；GTI（Good Time Interval）表示可用光子数据的时间范围；REG为光子数据的位置范围利用图像模式的X射线探测器获得的数据需要先根据观测源的位置，选定数据位置范围才能获得较高信噪比的观测源数据。

图1 利用FV读取B1509-58脉冲星的X射线观测数据Fig. 1 Using FV to read the X-ray observation data of B1509-58 pulsar

2 X射线数据初步处理

经过空间观测，能够得到X射线光子到达时间的原始数据，但由于像素坏点、坐标系转换、虚光、像素位置误差等因素的影响必须要对数据进行校准。美国加州理工大学的科学数据运行中心为NuSTAR任务[1]编制了用于数据校准的软件NuSTARDAS（NuSTAR Data Analysis Software），利用该软件可以完成NuSTAR数据的处理和分析等操作。

NuSTAR任务配备了两个相同的X射线光子焦平面收集组件（Focal Plane Module，FPM），分别为FPMA和FPMB，有效探测能级为3.0～79.0 keV。在NuSTARDAS软件中重要的处理命令为nupipeline，主要完成坐标系转换、姿态校准、坏像素剔除、能量校准、天球位置计算、根据校准文件库的数据剔除、数据可用时间段计算、由时间分辨率引起的无效时间的校准、曝光时间计算等处理，调用格式如下：

其中：indir为观测数据的位置；steminputs为观测数据的ID；outdir为经处理后的数据存放位置。

在完成数据校准后便得到了可用于科学分析的数据，由于空间背景的存在使得信噪比水平较低。因此，需要进行脉冲星源的选取，需要去掉背景辐射光子才能进行脉冲星的时间分析和谱分析。首先，需要用到DS9软件选取星源区域，再用Xselect命令选取指定区域的光子序列，获得的指定区域的X射线脉冲星光子序列是光子到达航天器的时间（Local Time）。由于航天器高速飞行产生的多普勒效应会使折叠出的脉冲轮廓并不精确，为了进行精确的时间分析必须将光子到达航天器的时间转换到光子到达太阳系质心的时间，并修正星载时钟的时差。Heasoft软件包集成了该时间转换命令barycorr调用格式为

其中：infile为需要进行太阳系质心时间转换的光子序列文件；outfile为时间转换的输出文件；orbitfiles为卫星轨道文件；clockfile为星载时钟的钟差文件；卫星任务组会定期发布，RA为脉冲星的赤经；DEC为脉冲星的赤纬。

barycorr软件中进行的时间转换为

其中：tb为光子到达太阳系质心的时间即脉冲星导航中利用的光子到达时间；tobs为光子到达航天器的观测时间；（clock）为将航天器的本地时间系统转换为地球质心力学时，包括星载时钟的修正，该项与具体的X射线观测望远镜有关。

在RXTE任务中，只要进行适当的时钟修正，星载记录时间就可以直接表达为地球质心力学时。（dispersion）为色散改正，对于X射线来讲色，散改正量为0，（geometric）为几何学时间延迟量，（"Einstein"）为由相对论效应引起的爱因斯坦延迟，（"Shapiro"）为沙皮诺延迟是由太阳系内引力场的光线弯曲引起的。

经上述一系列的数据处理后就可以进行时间分析和谱分析。时间分析主要进行周期检测、周期的一阶和二阶导数拟合、光子计数及计数率统计、基于历元折叠的轮廓分析、周期跃变、双星系统轨道参数分析等。工作谱分析主要进行星源能量谱的模型拟合、背景谱分析、物理辐射模型建立与验证等工作，往往需要结合具体的X射线望远镜的相关参数和文件进行分析。

3 X射线时间分析

3.1 周期检测

X射线时间分析的重要内容是进行星源的周期检测，包括脉冲星自身的旋转周期检测和双星系统中轨道周期的检测。这能够帮助我们了解脉冲星的年龄等物理参数，进而了解演化进程。在X射线脉冲星导航数据库的构建中，周期检测是基础，只有首先确定脉冲星星源的信号周期性以及轮廓特征才能得到导航脉冲星数据库的精确参数。

离散傅里叶周期检测适用于在周期未知情况下的检测，但计算量较大。因此，在已知脉冲星周期范围时往往采用周期检测或周期检测，这2种方法都在一定的周期范围内进行搜索极大值点，即与脉冲星周期相比离散傅里叶变换方法精度更高。周期检测方法较周期检测方法计算量大但有着更高的搜寻精度。

在周期范围了解的情况下，可以采用X射线天文学中广泛使用的检测标准（卡方检测）进行周期检测检验。按照某一周期折叠的脉冲轮廓是否具有足够好的分布，即：脉冲轮廓各个Bin中的光子数目距平均值的离散程度。在指定范围内获得离散程度最高的轮廓，其对应的周期值即为脉冲星的自转周期。

用于脉冲星轮廓检测的卡方计算标准定义为

图2 卡方周期检测Fig. 2 Period detection using Chi-square test

双星系统的轨道周期较长，每个Bin的时间间隔也相对较长。由于光子探测无效时间的存在，如果简单地利用各个Bin中的光子数目进行检测往往会造成很大的误差。因此，在双星系统轨道周期检测中采用光子计数率（单位时间内的光子数目）来进行卡方分布统计。定义式为

在多个搜索周期下，进行折叠将得到一系列的统计量S。由于光子到达时间为随机量，在背景噪声的影响下，卡方检测的最大值不一定是最优的结果。但大量S值能够表现出一定的统计特性。利用含有周期项的统计函数进行拟合，将使结果更具有可信性。Leahy提出的Epoch Folding方法[2]能够对卡方分布的S量进行拟合，使结果更为精确。在光子较少的情况下Epoch Folding方法同样具有很高的精度，Epoch Folding方法的拟合函数为

利用Epoch Folding方法进行4U 0142+61脉冲星的周期搜索，该脉冲星位于仙后座中，是非常规X射线脉冲星的一员，距地球约1.3万光年、自转周期约为8.69 s、周期一阶导数为0.2 × 10–11s s–1、两级磁场强度为

在2～10 keV能段，在非常规X射线脉冲星中具有较高的光度，约为1 × 1035erg/s。从红外波段到硬X射线波段均有辐射。在观测中也发现了该脉冲星的脉冲轮廓有缓慢的变化，并且存在X射线暴现象。NuSTAR任务在2014年3月对该脉冲星进行了观测。Epoch Folding方法在进行周期拟合后得到自转周期为P= 8.689 157 5 ± 3.799 ×10－6s，同其他文献中得到的自转周期一致，如图3所示。

图3 0142+61自转周期检测Fig. 3 0142+61 rotation cycle period detection

根据脉冲星时间相位模型每个光子到达时间都对应一个脉冲相位即

其中：N为光子总数；n为谐波分量的个数。

3.2 脉冲比例

在完成周期检测后，需要进行脉冲比例（Pulsed Fraction）的分析，用以评价折叠出的脉冲轮廓的强度及好坏。该分析结果在X射线时间分析中有重要作用。比较不同轮廓间的差别来分析脉冲辐射相对于背景辐射的强度。若要在脉冲星导航数据库建立中起到参考作用更倾向于选择脉冲比例大的脉冲星。脉冲比例的计算方法一般有以下几种。

图4 方法进行周期检测Fig. 4 Period detection using test

1）峰值脉冲比例（Peak-to-Peak Pulsed Fraction）

峰值脉冲比例以脉冲轮廓中的最小值和最大值进行脉冲轮廓的评估。

图5 正弦波形脉冲轮廓Fig. 5 Sinusoidal pulse profile

则峰值脉冲比例为

2）均方根脉冲比例（Root Mean Square Pulsed Fraction）

均方根脉冲比例，基于傅里叶变换通过傅里叶能谱评估脉冲轮廓的强度。由于该评价方法利用了轮廓中所有的bin值，相比峰值脉冲比例，有着更好的统计性能。

其中：N为脉冲轮廓中相位间隔的数目；pi为每个相位间隔中的光子数；k为傅里叶变换的谐波数目；和分别为ak和bk的不确定度。

3）面积脉冲比例（Area Pulsed Fraction）

从表达式上很容易理解面积脉冲比例计算的是脉冲振幅的面积与整个脉冲轮廓面积的比值反映了脉冲轮廓相对于噪声的强度。

在不同能级下，同一颗脉冲星的脉冲轮廓存在差异，这时便可以利用脉冲比例进行脉冲强度的分析。例如在不同能级下，4U 0142+61具有不同的累积脉冲轮廓，按照3～5、5～8、8～20、20～35、35～50、50～79 keV划分能量段可得到不同能级下的脉冲轮廓，如图6所示，从中可以看出在不同能级上脉冲轮廓的形态有很大的差异。

按照上述的能段划分利用Area方法进行脉冲比例的计算，可以得到不同能段的脉冲比例。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）。可知Area脉冲比例随着能级的增加持续增长，这种趋势同Hartog[4]得出的结论一致。

4 X射线谱分析

XSPEC软件是美国航空航天局（National Aeronautics and Space Administration，NASA）管理和发布的关于X射线谱分析的软件，能够处理多个X射线探测任务的数据、完成数种辐射模型的建立和拟合，并对拟合结果进行统计评估。通过图形显示和数值分析，帮助人们推断天体辐射的物理特性。XSPEC也能够模拟特定X射线探测任务的能谱，在科学任务分析和规划中起到重要作用，在高能X射线谱分析中应用十分普遍。在脉冲星导航中，可以用于之前X射线探测任务的结果，进行对比分析以确认脉冲星导航探测器是否存在问题，也可以在导航的同时兼顾科学分析。

在仿真中需要用到两个重要文件，即RMF文件和ARF文件，描述X射线探测设备的性能。RMF文件是响应矩阵文件（Response Matrix File），描述了X射线探测设备对X射线能量及收集通道的响应。ARF是辅助响应文件（Ancillary Response File），包含了不同能级。X射线光子的有效面积响应及其参数影响因素，主要有X射线探测设备的设置、视场位置、有效面积、量子效率、滤波器转换效能、晕影改正等因素。

在获得两个文件之后便可以建立辐射模型。XSPEC能根据辐射模型和设备的响应参数进行随机数打靶，获得较高置信度的模拟X射线能谱。对该能谱进行拟合，便能给出该探测设备在指定的观测时间内，所能获得的探测目标物理参数的精度。其中，几个比较重要的命令为model命令、fakeit命令和fit命令。

model命令可以根据设定的参数，进行光子能量谱模型的建立，用于分析星源的辐射模型。常用的辐射模型为黑体辐射谱模型（blackbody）和幂律谱模型（powerlaw），黑体谱可由黑体辐射引起，幂律谱可由热轫致辐射、同步加速辐射和逆康普顿散射引起。

XSPEC中的黑体辐射谱模型为

其中：k为玻尔兹曼常数；T为黑体绝对温度，单位为keV；为星源光度，单位为1039erg/s，为到星源的距离，单位为10 kpc；E为辐射能量。

图6 4U 0142+61不同能级的脉冲轮廓Fig. 6 4U 0142+61 pulse profiles with different energy level

XSPEC中的幂律谱模型为

fakeit命令能够根据X射线探测器的RMF和ARF文件，仿真出该X射线探测器能捕捉到的X射线光子能量数据，这可以作为模拟数据进行模型的拟合，具有很高的可信度，用于观测任务论证等场合。调用格式为

其中，nicer.rmf和nicer.rmf为NICER探测器对应的RMF和ARF文件，J1847.fak为仿真出的数据文件，200 000为仿真的观测时长，单位为s、y和blah，为格式输入。

fit命令根据建立的辐射模型，利用观测或模拟生成的数据进行模型拟合，并给出多种统计量以评判拟合结果的可信性，当可信性较高时，便可以接受模型，推断出星源的物理辐射形式。

仿真中用到的RMF文件和ARF文件，由美国宇航局戈达德航天中心提供。XSPEC中模型的参数需要参考相关探测目标的研究文献。但有些探测目标由于没有进行观测或观测时间较短，或X射线望远镜的观测能力有限，无法从文献中查找，此时可以简单地通过黑体辐射热力学和普朗克函数进行估计。

普朗克黑体辐射定律为

如果已知星源的半径、温度和距离，便能通过普朗克黑体辐射定律的积分求出指定能量段的辐射流量，进而得到黑体辐射模型的系数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。