X射线脉冲星导航的最优观测周期确定

苏剑宇，方海燕，，高敬敬，赵良

1．西安电子科技大学 空间科学与技术学院，西安 710126 2．北京临近空间飞行器系统工程研究所 空间物理重点实验室，北京 100101

X射线脉冲星导航（X-ray Pulsar-based Navigation，XPNAV）是一种新型的自主导航技术，可为航天器提供位置、时间等导航信息［1-4］，实现航天器高精度自主导航。目前，国内外相继开展了X射线脉冲星导航试验，如中国空间实验室天宫二号（Tiangong-2， TG-2）的γ暴偏振探测科学实验［5］，中国首颗X射线脉冲星导航试验卫星（XPNAV-1）在轨开展的X射线脉冲星的探测与脉冲星导航体制的验证［6-7］，以及国内首颗空间X射线天文硬X射线调制望远镜卫星（Hard X-ray Modulation Telescope，HXMT） Insight-HXMT的脉冲星定轨精度验证实验［8］。如美国NICER（Neutron star Interior Composition Ex⁃plorer）项 目 的SEXTANT （Station Explorer X-ray Timing And Navigation Technology）搭载国际空间站（International Space Station， ISS）开展的定轨精度验证工作［9］。可见，X射线脉冲星观测与导航验证工作正进入蓬勃发展的时期。XPNAV的基本原理可描述为［10］：在脉冲星观测周期内，航天器上安装的X射线探测器会记录到一串光子到达时间（Time of Arrivals，TOA），利用观测的光子TOA，通过一定的算法提取出在这段观测时间内的某一时刻处航天器所接收的脉冲星信号的相位以及多普勒频率，由于航天器在任意时刻观测的脉冲信号相位和频率都可以用该时刻的位置、速度以及太阳系质心（Solar System Barycenter， SSB）处的脉冲星信号模型准确表示，因此以估计得到的航天器在某一时刻的观测信号相位和多普勒频率作为导航观测量，并将其表示为有关航天器在该时刻位置和速度的关系式作为导航观测方程，获得任意时刻的观测相位即可得到相应时刻的位置信息。高精度的观测脉冲相位估计是整个导航的难点，为此，国内外学者也相继提出许多脉冲相位估计算法［11-18］。

在上述导航算法中，为获得高精度的位置与速度估计，需要实现观测相位和频率的高精度估计。文献［19］提出在星际转移轨道段可将航天器的运动轨迹简化为匀加速模型，基于此建模了航天器轨道状态误差与其所引起的脉冲TOA偏差之间的关系，并对导航观测方程进行了修正以消除TOA偏差的影响。但线性化过程或产生高阶截断误差，且该截断误差随脉冲星观测周期的增大而增大。所以通常为了将轨道位置误差进行线性化，选择压缩脉冲星观测周期或将观测周期分段。如文献［20］提出将观测时段划分为短时段，在每个较短时段内利用针对恒定多普勒频移模型的相位估计算法估计相位，再借助相位跟踪的方法提高单段相位估计精度并获得动态的多普勒频移估计。文献［21］提出了由轮廓形变来求解速度误差引起的多普勒频移的新思路，定义了轮廓熵以衡量累积轮廓的不同形变程度，推导了存在速度误差时，频率偏差与轮廓熵之间的关系，利用最优化方法通过使熵最小（轮廓最为清晰）从而求解出多普勒频移。该方法本质上仍需要假定在观测时段内的观测频率恒定，即航天器的速度恒定，所以观测周期不能过大。Anderson和Pines［22］研究了基于锁相环的脉冲相位跟踪方法，并给出了针对小流量脉冲星的相位跟踪方法。文献［23］将较短观测时段内的预报位置误差建模为随时间近似线性变化，且将这种线性变化的位置误差对光子到达时间修正的影响建模为恒定的相位偏差加上恒定的多普勒频移，在此基础上，提出了一种基于相位和多普勒频移联合观测的X射线脉冲星导航方法，并进行了实验验证与分析。以上算法在进行观测相位与多普勒频率估计时，未对脉冲星观测周期的选取进行深入研究。而最优观测周期的确定对脉冲星导航计算及观测计划制定等具有重要的意义，本文对该问题进行研究与讨论，以期为观测方案的制定提供参考。

1 观测相位估计误差

最优脉冲星观测周期是为得到高精度的观测相位估计，所以将观测相位估计误差最小时对应的脉冲星观测周期称为最优脉冲星观测周期。航天器在(t0，tf)的脉冲星观测周期内，任意时刻的观测相位［10］为

根据式（1），对ϕsc(t)的估计误差取决于对 δr的估计误差。任意时刻的位置误差δr为

式中：T为脉冲星观测周期长度，T=tf−t0；δr0为航天器初始时刻t0的位置误差。则δr引起的相位误差δϕ为

式中：δϕ0为对应t0时刻的初始预测相位。将式（3）展开到2阶，表示为

根据上述分析，对δϕ的估计将同时存在方差与偏差，采用均方误差准则，并用mse(δϕ)表示δϕ的估计均方误差，则均方误差mse(δϕ)为

式中：var(δϕ)为δϕ的估计方差；b2为δϕ的估计偏差，。方差是由θ的估计方差引起的，所以求出估计方差，即可得到使mse(δϕ)最小的T的表达式。

根据文献［12］，利用探测器探测到的光子到达时间序列，使用最大似然估计可实现对脉冲相位延迟的估计，利用最大似然估计法对θ的估计可表示为

式中：λs、λb分别为脉冲星源的流量与背景流量；h(ϕ)为相位ϕ∈[0，1)的脉冲星的标准轮廓，满足=1，h(n+ϕ)=h(ϕ)（n为整数）［12］。Emadzadeh证明了利用式（6）的最大似然估计对脉冲相位延迟的估计方差即为克拉美罗界（Cramer-Rao Lower Bound， CRLB）。所以下面推导对θ估计的CRLB。

2 最优脉冲星观测周期确定

2.1 θ的CRLB推导

任意时刻的相位误差δϕ为

将θ的CRLB表示［13］为

式中：I(θ)为Fisher信息矩阵，根据Emadzadeh和Speyer［13］的结论，I(θ)的第i行第j列元素Iij的表达式为

其中：

由Iij的表达式可得Iij=Iji，所以I为对称矩阵，且

为确定I(θ)矩阵，需要求I11、I12、I13、I23、I335个元素。Emadzadeh给出了I11、I12、I13的表达式，附录A给出了I23、I33的推导过程，结果如式（A19）所示。

得到I11、I12、I13、I23、I33的表达式，即可求得的CRLB(θ)。附录B给出了CRLB(θ)的推导过程。直接求解化解CRLB(θ)难以化解，这里首先化解

式中：a为I(θ)的行列式，将I11、I12、I13、I23、I33代入式（12），化解得到

所以CRLB(δϕd0)=。同理得到

所以

综上可得

式（17）即为对参数θ估计的CRLB。

2.2 最优观测周期确定

得 到θ=的CRLB后，可 确 定δϕ的估计误差。根据式（17）、式（4）中δϕ的估计方差var(δϕ)可表示为

最优观测周期是指对δϕ的估计均方误差最小时的观测周期，即所选观测周期应使式（19）中mse(δϕ)最 小。对 式（19）求 极 值，得 到 使mse(δϕ)最小的观测周期的表达式为

至此得到了最优脉冲星观测周期。需指出，式（20）所定义的最优观测周期是在式（1）所示的航天器处任意时刻观测相位模型下推导，以该观测相位估计误差最小为准则得到。

3 仿真实验验证

为了验证式（20）的正确性，对同一颗脉冲星在不同观测周期下进行仿真实验，并给出仿真的mse(δϕ)随观测周期变化的曲线，通过对曲线进行拟合，得到最小mse(δϕ)对应的观测周期，并与利用式（20）所预测的最优观测周期进行对比。

3.1 θ估计方差的仿真

脉冲星观测信号的仿真是后续实验的基础，所以首先需要仿真航天器在脉冲星观测周期内的脉冲星信号。实验中采用的轨道根数如表1所示，所使用的脉冲星参数如表2所示。光子TOA仿真方法参考文献［24］。

表1 轨道参数Table 1 Orbit parameters

表2 脉冲星参数Table 2 Pulsar parameters

利用上述脉冲星各项参数与轨道进行光子序列仿真，航天器初始的偏差与噪声方差设置如表3所示。

表3 初始误差Table 3 Initial error

利用预测位置偏差对光子TOA进行时间校正，得到修正的光子TOA，然后代入式（6）的算法求̇。采用利用格点搜索法［25］求式（6），可知求解过程为三维的搜索过程，的3σ搜索 范 围 分 别 为与。

重复仿真光子TOA，利用每一次仿真的光子TOA进行一次搜索，得到θ的第i次估计值，进 行800次 重 复 试 验，即i=1，2，…，800。然 后 利 用800次 结 果，可 求 得δϕ0、δf0、δḟ0的 均 方 误 差 mse(δϕ0)、mse(δf0)、mse(δḟ0)。根据式（18）、式（19），可得mse(δϕ)的仿真结果为

3.2 不同观测周期下的mse(δϕ)

改变观测周期，重复3.1节中计算mse(δϕ)的步骤，可得到mse(δϕ)随观测周期变化的曲线。脉冲星PSR B1821-24仿真结果如图1、图2所示。

图1给出了对θ估计的均方根误差（Root Mean Square Error， RMSE）结果，根据2.1节的推导，CRLB即为对θ估计的RMSE的理论值，根据图1的仿真结果，的估计均方根误差随观测周期的增加首先逐渐靠近其CRLB，证明推导的θ的CRLB的正确性。但当观测周期继续增大，的估计均方根误差逐渐偏离了CRLB，这是因为对式（3）进行了截断，且截断误差近似以代替，该误差会随观测周期的增加而增大。

图1 对θ估计的均方根误差（PSR B1821-24）Fig. 1 Root mean square error of estimation of θ（PSR B1821-24）

图2（a）为mse(δϕ)随观测周期变化的曲线，其横纵坐标都进行了对数运算。根据仿真结果，mse(δϕ)随观测周期先减小后增大，且存在使mse(δϕ)最小的点，在该点附近拟合，得到图2（b），根据拟合曲线可明显看出极小值点的存在，并可得到mse(δϕ)最小时对应的观测周期为620 s，利用式（20）得到仿真条件下的最优观测周期Topt=673 s，与实验结果相差53 s，可见预测最优观测周期在实际最优观测周期附近，证明了该公式的正确性。

图2 mse(δϕ)随观测周期变化曲线（PSR B1821-24）Fig. 2 Curve of mse(δϕ) with observation period （PSR B1821-24）

同理，图3、图4给出了脉冲星PSR B0531+21的仿真结果。图3为θ估计的均方根误差结果，可看出，δϕ0、δf0、的估计均方根误差随观测周期的增加首先逐渐靠近其CRLB，证明推导的θ的CRLB的正确性。当观测周期继续增大，由于截断误差的影响，δϕ0、δf0、的估计均方根误差逐渐偏离了CRLB。

图3 对θ估计的均方根误差（PSR B0531+21）Fig. 3 Root mean square error of estimation of θ（PSR B1821-24）

图4（a）为仿真得到的mse(δϕ)随观测周期变化的曲线，其横纵坐标都进行了对数运算。根据仿真结果，mse(δϕ)随观测周期先减小后增大，且存在使mse(δϕ)最小的点；图4（b）为在最小值附近拟合的结果，曲线的最小值即为仿真得到的最优观测周期。

图4 mse(δϕ)随观测周期变化曲线（PSR B0531+21）Fig. 4 Curve of mse(δϕ) with observation period （PSR B0531+21）

根据拟合曲线可明显看出极小值点的存在，并可得到mse(δϕ)最小时对应的观测周期为435 s，利用式（20）得到仿真条件下的最优观测周期391 s，与实验结果相差44 s，可见预测最优观测周期在实际最优观测周期附近，证明了该公式的正确性。

为进一步验证公式的正确性，使用NICER的实测数据进行了验证。数据包为“ni1013010142_0mpu7_cl.enents.gz”［26］。该 数据包为PSR B0531+21的实测数据，包含了40279473个光子TOA。表2给出了脉冲星参数信息，轨道根数与表1数据相同。利用该实测数据对θ=进行搜索，并改变观测周期，重复上述步骤，可得到mse(δϕ)随观测周期变化的曲线，图5（a）为mse(δϕ)随观测周期变化的曲线，其横纵坐标都进行了对数运算。通过在mse(δϕ)最小的点附近拟合，得到拟合结果如图5（b）所示。

图5 mse(δϕ)随观测周期变化曲线（PSR B0531+21实测数据）Fig. 5 Curve of mse(δϕ) with observation period （Ob⁃servational data of PSR B0531+21）

根据拟合曲线可明显看出极小值点的存在，并可得到mse(δϕ)最小时对应的观测周期为460 s，利用式（20）得到仿最优观测周期Topt=410 s，与实验结果相差50 s，可见预测最优观测周期在实际最优观测周期附近，且将实测数据的验证结果与图4的仿真结果对比，可看出在相同条件下实测数据与仿真数据的结果是接近的，证明了式（20）的可靠性。

使用NICER对PSR B1821-24的观测数据，进一步进行实验验证。数据包名称为“ni0070010102_0mpu7_cl_t1”“ni0070010103_0mpu7_cl_t1”“ni0070010104_0mpu7_cl_t1”［26］，共19971个光子TOA，结果如图6所示。

从图6（a）可看出，mse(δϕ)随观测周期的变化规律与图5相同，根据拟合曲线可明显看出极小值点的存在，并可得到mse(δϕ)最小时对应的观测周期为650 s，利用式（20）得到仿真条件下的最优观测周期714 s，与实验结果相差64 s，即预测最优观测周期在实际最优观测周期附近，证明了式（20）的正确性。且将实测数据的验证结果与图2的仿真结果对比，可看出在相同条件下实测数据与仿真数据的结果接近，证明了所提公式的可靠性。

图6 mse(δϕ)随观测周期变化曲线（PSR B1821-24实测数据）Fig. 6 Curve of mse(δϕ) with observation period （Ob⁃servational data of PSR B1821-24）

4 结论

在X射线脉冲星导航中，观测周期的选取是获得高精度观测相位估计的关键，本文对该问题进行了研究，得到以下主要结论：

1）以航天器处观测相位估计的均方误差最小为准则，推导了脉冲星观测周期与观测相位估计的均方误差的关系，给出了最优脉冲星观测周期的近似计算公式。

2）设计了仿真验证实验，仿真结果表明给出的最优观测周期计算公式对脉冲星PSR B1821-24与PSR B0531+21的预测误差约为64、44 s，而由脉冲星PSR B0531+21的实测数据得到的预测误差为50 s，即由该公式给出的最优观测周期与实验中对观测相位估计的均方误差最小时对应的观测周期接近，证明所给公式的正确性。

附录A：

定义

根据Emadzadeh和Speyer［12］的结论，I11、I12、I13可表示为

式中：

以下将推导I23、I33的表达式。

将观测周期分成N+1段，即

图A1 观测周期分段Fig. A1 Observation period segmentation

当n=0，Pn=0，有

将观测时间分成N+1段，将式（A4）的积分表示为图A1中每个小区间的积分求和，I23可表示为

根据图A1，在第n+1个区间内，脉冲星周期为Pn+1，则积分变量t可表示时刻对应的相位ϕ(t)为

将式（A6）代入式（A5），以ϕ为新的积分变量，并利用h的周期性，可化解得到

由于N≫1，所以式（A7）中占据主导地位，即

同理，I33的结果为

下面对I23、I33的表达式进行简化处理。I23的求和项的多项式展开式中，去除包含P0的项，共有项非零项，定义Pˉ23为

相应的，I23、I33可表示为

式中：q23为形如的非零项。所以

综上，由于N≫1，所以可近似认为=0，即。同理可证。

将I11、I12、I13表 示 为I23、I33的 形 式，则I11、I12、I13、I23、I33表达式为

附录B：

I−1(θ)的求解：根据逆矩阵运算法则，I−1(θ)可表示为

则θ=中的每一个元素的CRLB可表示为