CN104990548A - 基于历元折合的动态脉冲星信号的处理方法 - Google Patents

Abstract

本发明提供了基于历元折合的动态脉冲星信号的处理方法，包括步骤：根据脉冲星的方向矢量、航天器相对于太阳的预估位置，去除光子到达时间序列的轨道效应，产生修正光子到达时间序列；确定修正光子到达时间序列的周期；根据周期对改进光子到达时间序列进行历元折合，获得经验轮廓；根据经验轮廓，计算得到航天器探测初相。本方法充分考虑了航天器轨道运动对光子到达时间序列的影响，误差小，可操作性强，更具有实用性。

Description

基于历元折合的动态脉冲星信号的处理方法

技术领域

本发明涉及导航方法领域，特别地，涉及一种动态脉冲星信号的处理方法，适用于不同飞行任务的X射线脉冲星导航任务。

背景技术

X射线脉冲星导航是一种新兴的天文导航方法，具有精度高、适用范围广的特点。X射线脉冲星是一种高速自旋的中子星，长周期稳定度极佳，且能辐射X射线的脉冲信号。对脉冲信号进一步处理可以获得航天器相对于脉冲星的位置，其可作为航天器导航的基准。

虽然脉冲星的信号被认为是脉冲型，可由于脉冲星的流量及其微弱，航天器仅能接收到一序列光子到达时间(Time of Arrival，TOA)。如何从光子TOA中提取出脉冲到达时间是脉冲星导航的一项关键技术。该技术可通过两类方法实现，即历元折合法(Epoch Folding，EF)和直接利用光子到达时间法(Direct Use of photon TOAs，DUPT)。若航天器静止或者以已知的速度做匀速直线运动，EF和DUPT均可获得较好的效果。若匀速直线运动的速度未知，EF失效。若航天器未进行匀速直线运动，现行的EF和DUPT均失效。在实际情况中，航天器均沿着轨道运行，航天器的速度实时变化，无法简单近似为匀速直线运动。因此提出可应用于该情况的动态脉冲星信号处理方法，对于促进脉冲星导航的发展十分必要。

针对航天器轨道运动的情况，现有的方法大多是从修正导航模型出发。文献1(褚永辉,王大轶,熊凯,黄翔宇.X射线脉冲星导航测量延时补偿方法研究[J].2012,33(11):1617-1622.)通过添加附加轨道影响项来修正脉冲到达时间。然而，该方法并未考虑航天器轨道运动导致脉冲星信号的周期变化，讨论的问题较为简单，且并未涉及具体的脉冲到达时间或脉冲初相的计算。文献2(刘劲,曾宪武,房建成,宁晓琳.基于星光多普勒的脉冲星脉冲到达时间补偿[J].2014,42(1):129-132.)改进了文献1提出的方法，然而该方法同样并未考虑轨道运动导致的脉冲星信号周期变化，仅将轨道影响看做对脉冲到达时间的累积影响。此外，文献1和文献2所提出的方法均通过引入其他导航设施的辅助信息来实现目的，并未从根本上解决单独脉冲星导航的动态信号求解问题。

发明内容

本发明要解决的技术问题为：为了从信号处理的角度，完善解决深空探测器的动态脉冲星信号处理问题，本专利提出一种基于历元叠加的动态脉冲星信号处理方法。

为解决上述技术问题，本发明采用以下技术方案：

一种基于历元叠加的动态脉冲星信号处理方法，步骤为：

(1)根据公式(I)去除光子到达时间序列的轨道效应，产生修正光子到达时间序列

$t_i-t_0+\frac{1}{c}n\·\left(\tilde{r}\right(t_i)-\tilde{r}(t_0\left)\right)=\stackrel{\‾}{t_i}-\stackrel{\‾}{t_0}---\left(I\right)$

其中，c为光速，n为脉冲星的方向矢量， $\tilde{r}\left(t_i\right)$ 为航天器在t_i时刻相对于太阳的预估位置，t₀为观测时段的起点。

(2)确定修正光子到达时间序列的周期。

可采用χ²检验等方法(葛明玉.脉冲星的X射线辐射特性研究[D].中国科学院研究生院,2012)来搜索修正光子到达时间序列的周期。

(3)根据搜索得到的的周期P_t，对进行历元折合，获得经验轮廓p。

根据文献(Emadzadeh,A.A.,Speyer,J.L..Navigation in Space by X-ray Pulsars[M].SpringPress,2011.)给出的方法，对进行历元折合，从而获得经验轮廓。

(4)根据第3步得到的经验轮廓，计算得到t₀时刻的航天器探测初相φ₀。

获得经验轮廓之后，可利用互相关法或非线性最小二乘法计算得到t₀时刻的航天器探测初相φ₀。

与现有技术相比，本发明的优点在于：

(1)本发明提出的动态脉冲星信号处理方法可去除航天器轨道运动对光子到达时间序列的影响。由文献[]可知，脉冲星动态信号求解的难点在于，航天器接收到的光子到达时间信息与航天器轨道耦合。由于航天器进行轨道运动，航天器在轨接收到的信号频率呈非线性变化。利用本发明扣除轨道效应后，可将原来的具有非线性变化频率的信号变成具有未知常频率的信号。进而，可通过采用周期搜索的方法，估计得到信号的频率，大幅降低了脉冲星信号动态求解的难度过程如下：

首先，根据下式(I-I)计算t_i时刻脉冲星的预估相位φ_pre(t_i)：

${\mathrm{\φ}}_{pre}\left(t_i\right)=\frac{f_s}{c}n\·\left(\tilde{r}\right(t_i)-\tilde{r}(t_0\left)\right)---(I-I)$

其中，f_s为脉冲星的自转频率，c为光速，n为脉冲星的方向矢量，为航天器在t_i时刻相对于太阳的预估位置，t₀为观测时段的起点；

然后，根据公式(I)去除光子到达时间序列的轨道效应，产生修正光子到达时间序列

$t_i-t_0+\frac{1}{c}n\·\left(\tilde{r}\right(t_i)-\tilde{r}(t_0\left)\right)=\stackrel{\‾}{t_i}-\stackrel{\‾}{t_0}---\left(I\right)$

其中，c为光速，n为脉冲星的方向矢量，为航天器在t_i时刻相对于太阳的预估位置，t₀为观测时段的起点。

公式(I-I)和公式(I)可去除原光子到达时间序列的轨道效应。与现有的脉冲星信号处理方法相比，本发明的方法适用于单独的脉冲星导航，不需要引入外部测量信息，更具有实用性。

从仿真结果看，随着观测时间的增加，本发明方法求得的航天器探测初相的估计误差从0.002逐渐降低至0.0006左右。但现有技术中未扣除轨道效应方法的估计误差则随着时间，从0.002逐渐上升至0.11以上。这说明本申请方法得到的脉冲星信号估计效果可随着时间逐渐提高，更加可靠。

(2)本发明提出的方法可操作性强，易推广和使用，为其推广提供了条件。

附图说明

图1是本发明具体应用实例的流程示意图。

图2给出了本发明方法初相估计误差结果示意图。

具体实施方公式

以下将结合附图和具体实施例对本发明作进一步详细说明。

在具体应用实例中，所采用的航天器为一种深空探测器。

如图1所示，本发明的具体步骤如下：

(1)根据公式(I)去除光子到达时间序列的轨道效应，产生修正光子到达时间序列其过程如下：

首先，根据下式(I-I)计算t_i时刻脉冲星的预估相位φ_pre(t_i)：

${\mathrm{\φ}}_{pre}\left(t_i\right)=\frac{f_s}{c}n\·\left(\tilde{r}\right(t_i)-\tilde{r}(t_0\left)\right)---(I-I)$

其中，f_s为脉冲星的自转频率，c为光速，n为脉冲星的方向矢量，为航天器在t_i时刻相对于太阳的预估位置，t₀为观测时段的起点；

然后，当频率系数K取1时，根据公式(II-I)计算t_i时刻的航天器探测相位φ_det(t_i)：

${\mathrm{\φ}}_{\det }\left(t_i\right)={\mathrm{\φ}}_0+\underset{k=1}{\overset{K}{\Σ}}{f}_{dy}^{\left(k\right)}{(t_i-t_0)}^k+{\mathrm{\φ}}_{pre}\left(t_i\right)---(II-I)$

其中，φ_det(t_i)为航天器接收到的光子相位，φ₀为光子序列的初相，为第k项频率系数；φ_pre(t_i)为预估相位；

再次，将公式(I-I)代入公式(II-I)，得到：

${\mathrm{\φ}}_{d\mathrm{\α}}\left(t_i\right)={\mathrm{\φ}}_0+(f_{dy}^{\left(1\right)}+f_s)(t_i-t_0)+\frac{f_s}{c}n\·\left(\tilde{r}\right(t_i)-\tilde{r}(t_0\left)\right)---\left(1\right)$

其中，φ_det(t_i)为航天器接收到的光子相位，φ₀为光子序列的初相，f_s为脉冲星自转频率，为第1项频率系数。

因此，

$t_i-t_0+\frac{1}{c}\frac{f_s}{(f_s+f_{orbit}^{\left(1\right)})}n\·\left(\tilde{r}\right(t_i)-\tilde{r}(t_0\left)\right)\≈t_i-t_0+\frac{1}{c}n\·\left(\tilde{r}\right(t_i)-\tilde{r}(t_0\left)\right)---\left(2\right)$

令 $\stackrel{\‾}{t_i}={\mathrm{\φ}}_{\det }\left(t_i\right){(f_s+f_{orbit}^{\left(1\right)})}^{-1}$ 和 $\stackrel{\‾}{t_0}={\mathrm{\φ}}_0{(f_s+f_{orbit}^{\left(1\right)})}^{-1},$ 可得

$t_i-t_0+\frac{1}{c}n\·\left(\tilde{r}\right(t_i)-\tilde{r}(t_0\left)\right)=\stackrel{\‾}{t_i}-\stackrel{\‾}{t_0}---\left(I\right)$

其中，c为光速，n为脉冲星的方向矢量，为航天器在t_i时刻相对于太阳的预估位置，t₀为观测时段的起点。

(2)确定修正光子到达时间序列的周期。

由的表达式可知，的周期与不同，因此需要周期搜索。本实例采用χ²检验法。按照试验周期P_t将该光子到达时间序列分为M段。然后，将M段长度为试验周期T的光子到达时间序列分为N份，即将一个周期的相位分为N份。各个试验周期的光子到达时间依次累积按相位到第一段上，得到第i个相位的光子数为n_i。建立χ²统计量为

${\mathrm{\χ}}^2=\mathrm{\Σ}{(n_i-\stackrel{\‾}{n})}^2/\stackrel{\‾}{n}---\left(3\right)$

其中，为平均光子数。一直变化试验周期P_t，使χ²达到最大值。χ²最大值对应的P_t即为的周期。

(3)根据搜索得到的的周期P_t，对进行历元折合，获得经验轮廓p(t)。

可采用文献(Emadzadeh,A.A.,Speyer,J.L..Navigation in Space by X-ray Pulsars[M].Spring Press,2011.)给出的方法，对进行历元折合，从而获得经验轮廓。该方法由以下三步构成：1)将收集的光子到达时间序列归化到一个脉冲周期内；2)将脉冲周期划分为Bin个均等的封装段，并将归化后的光子放入相应的封装段；3)统计每个封装段中的光子数量，绘制脉冲轮廓。

(4)根据第3步得到的经验轮廓，计算得到t₀时刻的航天器探测初相φ₀。

本实例采用互相关法来求初相。可通过如下的公式求出两组轮廓间的时延。其中，R(φ₀)为经验轮廓与标准轮廓的互相关函数。互相关函数达到最大值对应的相位即为所求的脉冲星信号初相。

${\widehat{\mathrm{\φ}}}_0={\mathrm{argmax}}_{{\mathrm{\φ}}_0}\left(R\right({\mathrm{\φ}}_0\left)\right)$

$R\left({\mathrm{\φ}}_0\right)=\frac{1}{N}\underset{k=1}{\overset{N}{\Σ}}p\left(kT\right)s(kT-{\mathrm{\φ}}_0)---\left(4\right)$

其中，N为两组轮廓数据的采样点个数，T为相位的采样间隔，s为标准轮廓。

为了表明前述方法的有效性，现以仿真验证。假设一个处于巡航段的深空探测器观测脉冲星PSR B0531-21。脉冲星信号的初相为0.3。航天器的初始导航误差为100km和100m/s。通过1000次蒙特卡洛仿真分析本发明的性能。图2给出了本发明方法和未扣除轨道效应方法的性能对比。由图2可知，随着观测时间的增加，在3000S的观测时间内，本发明方法求得的航天器探测初相的估计误差从0.002逐渐降低至0.0006左右。但现有技术中未扣除轨道效应方法的估计误差则随着时间，从0.002逐渐上升至0.11以上。可以预见，随着观测时间的进一步加长，本发明方法得到的初相估计误差会进一步缩小，而其他方法的误差会进一步扩大。因此，本发明的方法充分考虑了航天器的轨道运动影响，保证了脉冲星信号估计效果可随着时间逐渐提高。

以上所述仅是本发明的优选实施方公式，本发明的保护范围并不仅局限于上述实施例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应该提出，对于本技术领域的普通技术人员来说，在不脱离本发明原理前提下的改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (3)

1.基于历元折合的动态脉冲星信号的处理方法，包括步骤：

A、根据脉冲星的方向矢量、航天器相对于太阳的预估位置，去除光子到达时间序列的轨道效应，产生修正光子到达时间序列；其中，将修正光子到达时间序列定义为扣除轨道效应后的光子到达时间序列；

B、确定修正光子到达时间序列的周期；

C、根据周期，对修正光子到达时间序列进行历元折合，获得经验轮廓；

D、根据经验轮廓，计算得到航天器探测初相。

2.根据权利要求1所述的处理方法，其特征在于，包括步骤：

A、根据公式I去除光子到达时间序列的轨道效应，产生修正光子到达时间序列

$t_i-t_0+\frac{1}{c}n\·\left(\tilde{r}\right(t_i)-\tilde{r}(t_0\left)\right)=\stackrel{\‾}{t_i}-\stackrel{\‾}{t_0}---I$

其中，c为光速，n为脉冲星的方向矢量，为航天器在t_i时刻相对于太阳的预估位置，t₀为观测时段的起点；

B、确定修正光子到达时间序列的周期：

建立χ²统计量为

${\mathrm{\χ}}^2=\mathrm{\Σ}{(n_i-\stackrel{\‾}{n})}^2/\stackrel{\‾}{n}---II$

其中，为平均光子数，n_i为第i个相位的光子数；χ²最大值对应的P_t即为的周期；C、根据搜索得到的的周期P_t，对进行历元折合，获得经验轮廓p(t)；

D、根据经验轮廓，计算得到t₀时刻的航天器探测初相φ₀；

通过如下的公式Ⅲ求出两组轮廓间的时延；

${\widehat{\mathrm{\φ}}}_0={\mathrm{argmax}}_{{\mathrm{\φ}}_0}\left(R\right({\mathrm{\φ}}_0\left)\right)$

$R\left({\mathrm{\φ}}_0\right)=\frac{1}{N}\underset{k=1}{\overset{N}{\Σ}}p\left(kT\right)s(kT-{\mathrm{\φ}}_0)---III$

其中，N为两组轮廓数据的采样点个数，T为相位的采样间隔，s为标准轮廓。

3.根据权利要求1所述的处理方法，其特征在于，步骤A包括以下步骤：

A1、根据下式I-I计算t_i时刻脉冲星的预估相位φ_pre(t_i)：

${\mathrm{\φ}}_{pre}\left(t_i\right)=\frac{f_s}{c}n\·\left(\tilde{r}\right(t_i)-\tilde{r}(t_0\left)\right)---I-I$

其中，f_s为脉冲星的自转频率，c为光速，n为脉冲星的方向矢量，为航天器在t_i时刻相对于太阳的预估位置，t₀为观测时段的起点；

A2、当频率系数K取1时，根据公式II-I计算t_i时刻的航天器探测相位φ_det(t_i)：

${\mathrm{\φ}}_{\det }\left(t_i\right)={\mathrm{\φ}}_0+\underset{k=1}{\overset{K}{\Σ}}{f}_{dy}^{\left(k\right)}{(t_i-t_0)}^k+{\mathrm{\φ}}_{pre}\left(t_i\right)---II-I$

其中，φ_det(t_i)为航天器接收到的光子相位，φ₀为光子序列的初相，为第k项频率系数；φ_pre(t_i)为预估相位；

A3、将公式I-I代入公式II-I，得到：

${\mathrm{\φ}}_{\det }\left(t_i\right)={\mathrm{\φ}}_0+(f_{\mathrm{dy}}^{\left(1\right)}+f_s)(t_i-t_0)+\frac{f_s}{c}n.(\tilde{r}\left(t_i\right)-\tilde{r}\left(t\right)0)---\left(1\right)$

其中，φ_det(t_i)为航天器接收到的光子相位，φ₀为光子序列的初相，f_s为脉冲星自转频率，为第1项频率系数；

A4、令可得：

$t_i-t_0+\frac{1}{c}n\·\left(\tilde{r}\right(t_i)-\tilde{r}(t_0\left)\right)=\stackrel{\‾}{t_i}-\stackrel{\‾}{t_0}$

其中，c为光速，n为脉冲星的方向矢量，为航天器在t_i时刻相对于太阳的预估位置，t₀为观测时段的起点。