在轨X射线脉冲星计时模型构建方法与流程

本发明涉及天文导航领域，尤其涉及一种在轨x射线脉冲星计时模型构建方法。

背景技术：

x射线脉冲星导航是一种新型自主导航方式，可以在单一仪器上自主完成航天器位置、速度和时间等信息测量，能够为近地、深空以及行星际飞行航天器提供有效的导航与授时解决方案，被认为是未来深空自主导航领域最具潜力的导航与授时方式之一。

构建x射线计时模型是脉冲星导航系统的基础性工作。对于脉冲星自主导航而言，随着传输距离的增加以及星际遮挡等原因，无法依靠地面测控系统实时、连续地传输射电计时模型参数，此外，当导航脉冲星存在周期跃变等自转周期不规则现象时，需要利用x射线观测数据，重新评估原模型参数，并建立新的脉冲星计时模型。因此，利用x射线观测数据独立自主地建立计时模型对航天器自主导航具有明显的支撑作用。

构建x射线计时模型是获得脉冲星时的关键技术。脉冲星辐射的信号的重复周期均匀性好、频率间隔精度高、并且能连续稳定运行，符合公认的钟的原理与特征，在建立脉冲星时过程中，需要描述脉冲星的自转频率、脉冲到达时间与原子时的数学关系，因此，精确的脉冲星计时模型是分析和建立脉冲星时的关键环节。

构建x射线计时模型还是脉冲星辐射机理研究的重要补充手段。对于有些类型的脉冲星，比如psrj0205+6449，先测到x射线接着才测到射电波段信号，随着x射线巡天的深入开展，未来会发现更多此类脉冲星，在未知其射电计时模型时，建立其x射线计时模型研究其辐射机理就显得至关重要。同时，建立x射线计时模型还可以与射电计时模型相互验证，为射电计时模型中存在的色散等误差源提供参考。

2016年10月8日，我国的“脉冲星试验卫星”成功发射，目前也已将构建x射线计时模型列为此次观测试验的一项科学和技术目标。

由于脉冲星在射电和x射线电磁波段均有辐射，所以脉冲星计时模型又可以分为射电计时模型和x射线计时模型，但与射电计时数据处理相比，x射线观测主要存在以下两个主要问题：

(1)x射线观测平台的动态性高。射电观测平台是地面的射电望远镜，而x射线观测平台则是星载x射线探测器。由于脉冲星观测可以认为是一种相对观测，地球自转周期24小时，而近地轨道卫星的自转周期仅100分钟，射电观测站的角加速度要远远低于卫星平台，两种平台的动态性会给信号处理带来较大差异，从而影响时间校正方法。

(2)x射线信号强度弱，射电信号通过大型天线也能转换成电信号，而x射线是随机分立的光子信息，这种弱信号的特性给x射线数据处理方法带来了很大的不同，包括到达时间校正方法、频率搜索方法、脉冲toa估计方法、以及在初始相位不精确条件下，计时模型的代价函数设计等，这些方法均需要重新改正与设计。

现有的射电计时模型处理方法较为成熟，但由于x射线数据形式、平台动态、信号强度等与射电观测有很大的差异，尤其是以上两个问题的存在，需要改进和设计x射线计时模型构建方法。

技术实现要素：

本发明的目的在于提供一种利用x射线观测数据计算脉冲星计时模型的方法框架。

本发明提供一种在轨x射线脉冲星计时模型构建方法，步骤一，x射线脉冲星光子到达时间预处理：将近地轨道航天器探测的光子到达时间的固有时转换到太阳系质心坐标系下的坐标时，对于航天器处测量的x射线光子到达时间假设用τobs表示，那么相应的ssb处的光子到达时间为

tssb＝τobs+△i+△c+△p+△r+△e+△s

公式中，△i为仪器响应延迟，用于修正x射线探测器的硬件响应时间延迟；

△c是时钟校正，该项校正由于记录x射线到达时间的原子钟长期稳定度降低而引入的时间误差；

△p是周年视差校正，校正地球绕太阳周年运动所产生的视差；

△r是roemer延迟，在真空中传播的同一光脉冲从航天器位置到ssb处的时间差；

△e为太阳系einstein延迟，该项是校正由于地球的运动，以及地球的引力势在卫星位置和ssb处的不同而导致的时延；

△sshapiro延迟，该项是校正光在经过太阳系内大质量天体时所引入的时间延迟；

步骤二，高精度的脉冲星自转频率计算：假设第i个光子到达时间分别用ti表示，对于假定的自转频率f，那么每个光子相位可以表示为

φi＝2πmod(fti,1)

mod为取余运算符，那么一种高精度的方法的代价函数表示为

其中的和分别为

最优的参数m的快速计算方法为

式中：为允许的积分平均误差，参数m可以通过给定的直接由该式计算；

步骤三，计时模型的初值区间估计：在建模x射线脉冲星计时模型时，需要计算并估计频率参数的初值区间；

步骤四，观测脉冲轮廓重构：在获得步骤二中的各数据包的自转频率后，利用历元折叠的方法获取观测脉冲轮廓；

步骤五，高精度的脉冲到达时间估计：根据观测脉冲轮廓模型估计初始相位的似然函数，计算脉冲峰值的达到时间。

与相关技术相比，本发明提供的在轨x射线脉冲星计时模型构建方法避免了传统的卡方量估计方法精度低带来的问题，同时解决了解决传统小波变换对脉冲轮廓消噪处理导致的信号急剧变化部分产生的震荡现象，同时保留脉冲轮廓的线性相位，提高计算精度。

附图说明

图1为本发明提供的在轨x射线脉冲星计时模型构建方法流程图；

图2本发明利用方法处理卫星实测数据得到的代价函数图；

图3为对7组卫星数据搜索的频率拟合的结果示意图；

图4为初始相位误差为0.2时，传统频率搜索代价函数得到的结果示意图；

图5为初始相位误差为0.2时，本发明提出的频率搜索代价函数得到的结果示意图。

具体实施方式

以下将参考附图并结合实施例来详细说明本发明。

如图1所示，所述在轨x射线脉冲星计时模型构建方法包括如下步骤：

1、x射线脉冲星光子到达时间预处理：将近地轨道航天器探测的光子到达时间的固有时转换到太阳系质心坐标系下的坐标时，对于航天器处测量的x射线光子到达时间假设用τobs表示，那么相应的ssb处的光子到达时间为

tssb＝τobs+△i+△c+△p+△r+△e+△s

公式中，△i为仪器响应延迟，用于修正x射线探测器的硬件响应时间延迟；

△c是时钟校正，该项校正由于记录x射线到达时间的原子钟长期稳定度降低而引入的时间误差；

△p是周年视差校正，校正地球绕太阳周年运动所产生的视差；

△r是roemer延迟，在真空中传播的同一光脉冲从航天器位置到ssb处的时间差；

△e为太阳系einstein延迟，该项是校正由于地球的运动，以及地球的引力势在卫星位置和ssb处的不同而导致的时延；

△sshapiro延迟，该项是校正光在经过太阳系内大质量天体时所引入的时间延迟。

该步骤与射电计时观测数据处理的差别：1)射电需要进行大气色散校正、星际色散校正，而x射线计时观测数据处理并不需要；2)射电仅需要对某个脉冲到达时间的固有时校正，而x射线需要对所有光子到达时间的固有时进行校正。

2、高精度的脉冲星自转频率计算：假设第i个光子到达时间分别用ti表示，对于假定的自转频率f，那么每个光子相位可以表示为

φi＝2πmod(fti,1)

mod为取余运算符，那么一种高精度的方法的代价函数表示为

其中的和分别为

最优的参数m的快速计算方法为

式中：为允许的积分平均误差，参数m可以通过给定的直接由该式计算。

利用方法得到的代价函数如图2所示，取代价函数最大值对应的频率即为最优频率，不过，为了提高频率估计的精度通常在最大频率附近采用插值拟合的方法，计算最优频率。

该步骤与现有频率计算方法的差别在于，这里提出一种的高精度频率估计方法，而且其中的参数m可以通过给定的允许的积分平均误差快速计算，避免了传统的卡方量估计方法精度低带来的问题。

3、计时模型的初值区间估计：在进行x射线脉冲星计时模型建模时，需要利用频率参数的初值区间。这里给出了一种基于最小二乘法的频率参数初值估计方法。

(1)将步骤二估计的m个数据包的脉冲星自转频率参数记为{(△t1,f1),(△t2,f2),…,(△tm,fm)}，其中△ti＝ti-t0，即观测时间与参考时间t0的时间间隔，那么测量数据可以表示为y＝[f1,f2,…,fm]^t，要估计的参数向量表示为β＝[f^(l)/l！,f^(l-1)/(l-1)！,…,f⁽⁰⁾/0！]^t，根据最小二乘法，有以下关系成立

y＝xβ+ε

其中的x用矩阵的形式展开为

(2)待估计的参数向量可以表示为

其中的r、q分别是矩阵x进行qr分解的上、下三角矩阵；

(3)的最小二乘估计协方差矩阵为

(4)对于中的每个参数相应的参数估计区间为

其中，βi＝f⁽ⁱ⁾/i！，tα/2(m-l-1)为学生氏t分布，置信水平为α，cii为矩阵(x^tx)^-1第i行、第i列的元素，方差的无偏估计量

如图3所示，利用最小二乘法确定的频率范围(29.48,29.96)，频率一阶导数范围(-5.543×10-10,1.279×10-9)，频率二阶导数范围(-1.564×10-18,1.731×10-19)。经过该步骤，就可以为计时模型中的各频率参数提供初始估计区间。

4、观测脉冲轮廓重构：在获得步骤二中的各数据包的自转频率后，利用历元折叠的方法获取观测脉冲轮廓。

(1)选择各数据包的中点时刻tmid作为历元折叠的参考时间。

(2)计算光子到达时间相对于中点参考时间的相位，对于第i个数据包，其相应的光子到达相位为φj＝fi(tj-tmid)。

(3)借助历元折叠方法计算观测脉冲轮廓

式中：n0为初始相位φ(t0)引入的相位整数，nb是历元折叠的bin块数量，np为积分时间内包含的整周期，tobs≈np/f，ci,k是第i个整相位第k个bin块内的光子数。

(4)借助具有平移不变性质的小波变换对观测脉冲轮廓进行消噪处理，此处使用平移不变小波消噪，以解决传统小波变换对脉冲轮廓消噪处理导致的信号急剧变化部分产生的震荡现象，同时保留脉冲轮廓的线性相位，以便于后续的脉冲到达时间估计处理。

5、高精度的脉冲到达时间估计：利用每个bin内的光子计数值的异方差泊松模型，估计初始相位φ0的似然函数为：

式中：λ(φ)为标准速率函数，观测脉冲轮廓的联合概率密度函数为：

其相对于φ0的最大值即可得到φ0的最大似然估计值的最大似然函数为：

其中：

脉冲到达时间差d可以表示为

脉冲峰值的到达时间为

tpk＝tmid+d。

该步骤提出一种快速最大似然脉冲到达时间估计方法，计算脉冲峰值的到达时间，该步骤中使用快速最大似然方法，具更高的计算精度。

6、鲁棒的计时模型估计：该方法对初始参考频率参数敏感，当初始参考频率参数区间的误差较大时，由于光子到达时间相位取整过程发生变化，导致估计精度降低，甚至产生错误估计。这里提出了一种鲁棒的频率估计方法，频率搜索的代价函数设计如下：

(1)根据步骤三的计算结果，设置中各参数的搜索区间。

(2)设置中各参数的搜索步长，步长越精细，搜索的精度越高，但相应的搜索时间也会增加。

(3)计算在各假定频率参数下的频率搜索的代价函数χ²(β)。

(4)将多维问题转换为一维的量，绘制代价函数χ²(β)曲线

(5)计算代价函数χ²(β)曲线的最小值，计算相应的频率参数值

该方法克服了传统方法在频率区间误差较大时，计算的计时模型出错的问题，而且该方法还可以在初始相位位置的条件下进行计算，计算方法的鲁棒性得到大大加强。

如图4、图5所示，当初始相位误差为0.2时，利用传统的计时模型构建方法得到的代价函数是发散的，代价函数的整体误差较大，在中心频率处也没有有效地降低，利用传统代价函数得到的参数估计结果为：频率估计误差为﹣3.5×10^-4，频率一阶导数估计误差为4.69×10^-13，频率二阶导数估计误差为1.97×10^-20。而利用本发明的方法，代价函数在中心频率处的代价函数值迅速减小，频率估计误差为﹣5.9067×10^-9，频率一阶导数估计误差为9×10^-15，频率二阶导数估计误差为9.7×10^-21。

从参数估计结果的比较来看，本发明提出的方法能够有效地降低初始相位误差对估计精度的影响，具有较好的鲁棒性。

以上所述仅为本发明的实施例，并非因此限制本发明的专利范围，凡是利用本发明说明书及附图内容所作的等效结构或等效流程变换，或直接或间接运用在其它相关的技术领域，均同理包括在本发明的专利保护范围内。