X射线脉冲星短时观测高信噪比轮廓获取方法与流程

本发明属于信号处理技术领域，特别涉及是一种短时观测获取高信噪比脉冲轮廓的方法，可用于对天文数据的处理。

背景技术：

x射线脉冲星导航是未来最具潜力的深空自主导航，现阶段正在向工程应用发展，美国国家航空航天局nasa在2017年11月开展的在轨实验中，选取4颗毫秒脉冲星作为信标，对每个脉冲星信标大约观测5～15min，然后自主旋转至下一个脉冲星信标，为期2天的实验表明，该全自动导航系统实现了世界上首次脉冲星三维定位验证，且实现了精度为16km以内的预定目标，最高精度约4.8km。高信噪比脉冲轮廓是实现脉冲星导航应用的关键之一。目前在脉冲星观测数据中获取脉冲轮廓的方法有以下两种：

第一种方法是周期折叠，该方法是将每个光子到达时间按时间周期折叠到一个脉冲周期内等分的bin块中，统计每个bin块中的光子数，以获得累积脉冲轮廓。

第二种方法是历元折叠，该方法是先计算每个光子到达时间相对于参考历元的相位，并将其归算到一个相位周期内的各个bin块中，再统计每个bin块中的光子数，以获得累积脉冲轮廓。

以上这两种方法，由于航天器收到的光子包含脉冲星信号和噪声，需要长时间的观测脉冲星来获取高信噪比脉冲轮廓，而长观测时间将影响x射线脉冲星导航xpnav的实时性、准确性和增加计算量，且增加了导航难度，很难满足实际工程要求。

技术实现要素：

本发明的目的在于针对上述现有技术的不足，提出一种x射线脉冲星短时观测高信噪比轮廓获取方法，以减少脉冲观测时间，并获得高精度的脉冲轮廓。

为实现上述目的，本发明的技术方案包括如下步骤：

(1)对光子到达时间序列进行折叠，获取时域波形：

1a)选取一段观测时间t的光子到达时间ti，i＝1,2,3...，根据脉冲星星历的自转频率参数，通过脉冲星计时模型计算光子到达时间ti所对应的相位φ(ti)；

1b)选用(2¹¹，∞)个bin块对光子到达时间ti进行周期折叠或对ti对应的相位φ(ti)进行历元折叠，获得时域波形；

(2)对时域波形进行傅里叶变化获得频谱；

(3)对频谱进行信号与噪声的分离：

3a)选用代价函数，计算频谱的截断点；

3b)根据信号的频域主要集中在低频，而噪声分布在整个频域且幅值较小的不同频域特性，利用截断点对频谱进行截断，选取0点至截断点之间的频谱，使得信号与噪声分离；

(4)恢复高信噪比时域观测轮廓：

4a)对截取的频谱去掉直流分量，即去掉零次频率分量的值；

4b)对去掉零次频率分量的频谱进行傅里叶逆变换得到不包含直流分量的时域复信号，对不包含直流分量的时域复信号取实部得到不包含直流分量的时域观测轮廓；

4c)对零次频率分量的值进行傅里叶逆变换得到直流分量；

4d)不包含直流分量的时域观测轮廓加上直流分量得到高信噪比时域观测轮廓。

本发明与现有技术相比具有以下优点：

1、本发明根据信号的频域主要集中在低频，而噪声分布在整个频域且幅值较小的不同频域特性，通过选取部分频谱对信号和噪声进行有效分离，实现了短时观测获取高精度观测轮廓；

2、本发明由于选取了(2¹¹，∞)个特殊的bin块数，并通过频域截断点有效分离其信号和噪声，在获取相同信噪比轮廓的情况下，缩短了计算时间，提高了计算效率。

附图说明

图1为本发明的实现总流程图；

图2本发明中周期折叠原理图。

具体实施方式

以下参照附图对本发明的实施例作进一步详细描述。

参照图1，本实例的实现步骤如下：

步骤1.选取光子到达时间ti，并计算其对应的相位φ(ti)。

1.1)选取一段观测时间t的光子到达时间ti：

罗西rxte卫星在轨16年，提供了大量的x射线脉冲星观测数据，这些观测数据存档在美国国家航空航天局nasa飞行中心高能物理研究中心的数据档案库中。

本实例选用蟹云crab脉冲星90802_02_02_00数据包，该数据包的观测时长为906秒，截取其前20秒的观测时长，作为光子到达时间ti；i＝1,2,3；

1.2)计算光子到达时间ti所对应的相位φ(ti)

根据脉冲星星历的自转频率参数，通过脉冲星计时模型计算光子到达时间ti所对应的相位φ(ti)，公式如下：

其中：φ(t0)是参考时间t0对应的初始相位，f^(l)分别是自转频率的l阶偏导数，l取值为2，f是脉冲星自转频率，取值为29.7925564233089hz；f⁽¹⁾是频率一阶导数，取值为-3.73397e-10；f⁽²⁾是频率二阶导数，取值为8.00e-21。

步骤2.对光子到达时间ti进行折叠。

折叠包括周期折叠和历元折叠，其中：

周期折叠，是将每个光子到达时间按时间周期折叠到一个脉冲周期内等分的bin块中，统计每个bin块中的光子数，以获得累积脉冲轮廓；

历元折叠，是先计算每个光子到达时间相对于参考历元的相位，并将其归算到一个相位周期内的各个bin块中，再统计每个bin块中的光子数，以获得累积脉冲轮廓。

参考图2，其折叠原理如下：

将光子到达时间或者对应的相位按照周期分成p份，按bin块数将每个周期进行等分，计算每个bin块内的光子数，将每个bin块内的光子数累加到第一个周期，进行归一化处理得时域波形；

所述对光子到达时间ti进行周期折叠，步骤如下：

2.1)将观测时长t为20秒的光子到达时间按脉冲周期ts分成p份，如图2(a)所示，ts值为0.0335654311027秒；

2.2)按bin块数m为2²³对每个周期进行等分，如图2(b)所示，计算第j个周期内第i个bin块内的光子数cj(ti)，i取值为[1，m]；j取值为[1,p]；

2.3)将第j个周期内第i个bin块内的光子数cj(ti)累加到第一个周期，得到未归一化的时域观测信号，如图2(c)所示；

2.4)将未归一化的时域观测信号进行归一化，得到时域观测信号波形如图2(d)所示，的表示式如下：

式中，m＝2²³为bin块数，tb＝ts/m，ts为脉冲星周期，值为0.0335654311027秒；cj(ti)为第j个周期内第i个bin块内的光子数。

所述对光子到达时间ti对应的相位φ(ti)进行历元折叠，其步骤如下：

第一步，将光子到达时间ti对应的相位φ(ti)按周期分成q份，其中，φ(ti)的周期值为1；

第二步，按bin块数m为2²³将每个相位周期进行等分，计算第j个整数相位第i个bin块内的光子数cj,i，i取值为[1，m]；j取值为[1,q]；

第三步，将第j个整数相位第i个bin块内的光子数cj,i累加到[0,1)之间，得到未归一化的时域观测信号；

第四步，将未归一化的时域观测信号进行归一化，得到时域观测信号波形表示式如下：

式中：cj,i为第j个整数相位第i个bin块内的光子数。

本实例采用但不限于对光子到达时间ti对应的相位φ(ti)进行历元折叠。

步骤3.对时域波形进行傅里叶变换。

根据傅里叶变换公式，对时域波形进行傅里叶变化获得频谱，其公式如下：

式中：m＝2²³为bin块数；为时域观测信号波形，为时域观测信号波形傅里叶变换后的频谱。

步骤4.计算频谱的截断点，并利用截断点分离信号与噪声。

4.1)计算频谱的截断点：

4.1.1)利用平均积分均方误差m.i.s.e构建出如下代价函数：

式中，sk为蟹云crab脉冲星标准轮廓傅里叶变换，为时域观测信号波形傅里叶变换后的频谱，为脉冲星标准轮廓傅里叶变换的共轭，为时域观测信号波形傅里叶变换的共轭，q为频域截断点，为代价函数值；

4.1.2)根据代价函数值随截断点q单调递减到最小值后一直保持单调递增的性质，当代价函数取最小值所对应的q，即为频谱的最优频率截断点，本步骤取代价函数的极值点为22，即截断点q为22；

4.2)信号与噪声分离

根据信号处理知识可知，信号的频域主要集中在低频，而噪声分布在整个频域且幅值较小的不同频域特性，在低频段信噪比高，噪声对信号影响小，高频段信噪比低，噪声对信号影响大，利用截断点对频谱进行截断，选取零点至截断点之间的频谱，使得信号与噪声分离，其公式如下：

式中：m＝2²³为bin块数，为时域观测信号波形，q为频域截断点，为通过截断点截取的频谱。

步骤5.恢复高信噪比时域观测轮廓。

时域观测轮廓包含直流分量和非直流分量，非直流分量即为不包含直流分量的时域观测轮廓，在恢复高信噪比时域观测轮廓时，需要分别恢复直流分量和不包含直流分量的时域观测轮廓，通过截断点截取的频谱包含直流分量的频谱和非直流分量的频谱，需要对截断点截取的频谱去掉直流分量，对去掉直流分量的频谱进行傅里叶逆变换得到不包含直流分量的时域观测轮廓，对直流分量的频谱进行傅里叶逆变换得到直流分量，其实现如下：

5.1)对通过截断点截取的频谱进行处理：

由于时域的直流分量对应频域零次频率分量的幅值，因此，对截取的频谱去掉直流分量，就是去掉零次频率分量的值，其公式如下：

式中：m＝2²³为bin块数；为时域观测信号波形；q为频域截断点，取值为22；为通过截断点截取的频谱，为对去掉零次频率分量的频谱；

5.2)恢复不包含直流分量的时域观测轮廓：

5.2.1)根据傅里叶逆变换公式，对去掉零次频率分量的频谱进行傅里叶逆变换，得到不包含直流分量的时域复信号，其公式如下：

式中：m＝2²³为bin块数；n为时域采样点数，其值为1024；q为频域截断点，取值为22；为通过截断点截取的频谱，为对去掉零次频率分量的频谱，为傅里叶逆变换后的时域复信号；

5.2.2)对时域复信号取实部得到不包含直流分量的时域观测轮廓，其公式如下：

式中：n为时域采样点数，其值为1024；为傅里叶逆变换后的时域复信号，为不包含直流分量的时域观测轮廓；

5.3)获取直流分量，并恢复高信噪比时域观测轮廓：

5.3.1)获取直流分量，根据傅里叶逆变换公式，频域零次频率分量的幅值对应时域的直流分量，其公式如下：

式中：m＝2²³为bin块数，为零次频率分量的幅值；

5.3.2)恢复高信噪比时域观测轮廓，即根据时域观测轮廓包含直流分量和不包含直流分量的时域观测轮廓，得到高信噪比时域观测轮廓，其公式如下：

式中：为不包含直流分量的时域观测轮廓，为高信噪比时域观测轮廓，其信噪比为17.71db，时域轮廓与标准轮廓的相关系数为0.9876。

以上描述仅是本发明的一个具体实例，并未构成对本发明的任何限制，显然对于本领域的专业人员来说，在了解了本发明内容和原理后，都可能在不背离本发明原理、结构的情况下，进行任何形式和细节上的各种修改和改变，但是这些基于本发明思想的修正和改变仍在本发明的权利要求保护范围之内。