基于脉冲星辐射信号先验信息的脉冲星自转频率估计方法

1.本发明属于信号处理技术领域，涉及一种x射线脉冲星自转频率估计方法，具体涉及一种基于脉冲星辐射信号先验信息的脉冲星自转频率估计方法，可用于x射线脉冲星导航。


背景技术：

2.脉冲星是一种高速旋转的中子星，致密的结构使其具有极强的引力和磁场，以及超高温和超高压；脉冲星的自转稳定性极强，其辐射波束周期性地扫过太阳系。x射线脉冲星导航xpnav作为一种新型导航方式，就是利用脉冲星稳定的自转周期进行导航，它可为航天器提供位置和时间等导航信息。脉冲星的辐射信号是开展脉冲星导航的基本数据资料。脉冲星辐射的x射线信号能量高，且其所需的探测设备能够实现小型化，可以搭载到航天器上，因此成为脉冲星导航的数据基础。当接收到x射线脉冲星的辐射信号时，探测器记录下每个x射线脉冲星光子到达探测器的时刻，得到x射线脉冲星光子到达时间toa序列，这是一个随机分布的非等间隔时间序列，该序列就是x射线脉冲星辐射信号的基本观测信息。获得x射线脉冲星光子toa序列后，利用脉冲星自转频率估计方法可从中得到脉冲星自转频率估计值，然后利用自转频率估计值对x射线脉冲星光子toa序列进行轮廓折叠以获得航天器处的脉冲星辐射信号的轮廓波形，将之与脉冲星辐射信号的标准轮廓波形进行比较就可获得脉冲星信号的时延估计值，这是xpnav的基本观测量，继而获得航天器位置和时间等导航信息。脉冲星自转频率的估计精度直接影响后续导航算法的计算精度，而且通过较短时间的观测获得较高精度的脉冲星自转频率估计值是xpnav能够顺利开展的基础。
3.目前对脉冲星自转频率的估计方法主要立足于beran提出的检验信号周期性的理论，该理论通过对x射线脉冲星光子toa序列采用自转频率估计值进行周期折叠，若该频率估计值接近真实频率，则折叠轮廓接近脉冲星辐射信号的标准轮廓，反之则折叠轮廓趋于一条直线，由折叠轮廓和标准轮廓的关联性可以设置代价函数进行评价。基于该理论的最常用方法就是pearsonχ2检验、z
m2
检验和基于fft的功率谱积分方法。
4.pearsonχ2检验是一种均匀性检验，属于时域方法，该方法先要将航天器处的x射线脉冲星光子toa序列转换为太阳系质心ssb处的光子toa序列，然后以脉冲星自转频率估计值对转换后的序列进行周期折叠，若频率估计值与真实频率相差较大，则折叠轮廓趋于一条直线，此时直流量增加、峰值量减小，χ2检验统计量较小，反之则折叠轮廓接近标准轮廓，直流量减小、峰值量增加，χ2检验统计量较大，当χ2检验统计量达到最大时，其对应的频率值就为最优频率估计值；该方法虽简洁高效，但易受高频噪声干扰，当观测数据信噪比较低时，χ2检验统计量会产生毛刺现象，导致统计量最大值的误判。
5.z
m2
检验是一种频域方法，该方法根据脉冲星辐射信号的有用信息主要集中在其频域的低频成分这一特点，直接计算x射线脉冲星光子toa序列的功率谱，对其功率谱的前m次谐波进行积分求和并以该和值作为z
m2
检验统计量，当自转频率估计值最优时，对应的z
m2
统计量最大；该方法避免了高频噪声对频率估计的影响，使其频率估计精度高于pearsonχ2检
验，但其计算量巨大，且计算机无法对x射线脉冲星光子toa序列这样的非等间隔数据进行fft运算。
6.根据以上两种基本方法的优势与缺陷，人们研究出改进的脉冲星自转频率估计方法，例如孙雄于2018年在他的硕士毕业论文《x射线脉冲星频率搜索及计时模型研究》中公开了一种基于fft的功率谱积分频率搜索方法。该方法是对z
m2
检验的改进，其先以自转频率估计值对x射线脉冲星光子toa序列进行周期折叠，然后对折叠轮廓进行fft变换并进一步计算得到其功率谱，之后设置频率截断点，对功率谱中一次谐波分量至截断点之间的成分进行积分求和并将该和值作为统计量，当自转频率估计值为最优时对应的统计量达到最大；该方法将非等间隔的x射线脉冲星光子toa序列转化为等间隔的折叠轮廓，使用fft辅助功率谱计算，降低了z
m2
检验的计算量，并基于k-t法确定了脉冲星辐射信号频谱的频率截断点的位置，为频率截断点的选取提供理论依据；但在60s以内的较短观测时间、脉冲星辐射信号数据累积量较少或有效观测数据量较少的应用条件下，观测数据信噪比的降低，导致统计量的计算受到负面影响，进而导致脉冲星自转频率估计值的精度降低。


技术实现要素：

7.本发明的目的在于针对上述现有技术的不足，提出了一种基于脉冲星辐射信号先验信息的脉冲星自转频率估计方法，用于解决现有技术中在短时观测条件下脉冲星自转频率估计精度较低的技术问题。
8.为实现上述目的，本发明采取的技术方案包括如下步骤：
9.(1)初始化参数：
10.初始化在观测时间[ta,tb]内航天器所探测的i个x射线脉冲星光子到达自身的toa序列为其中ti为第i个x射线脉冲星光子到达航天器的toa时刻，i＞10；
[0011]
(2)获取x射线脉冲星辐射信号的折叠轮廓：
[0012]
以h为步长对脉冲星自转频率估计范围[f
l
,fh]进行遍历，得到p个预设频率值通过用每个预设频率值f
p
对x射线脉冲星光子toa序列进行bin块数为m的周期折叠，得到对应的x射线脉冲星辐射信号的折叠轮廓的集合其中f
l
和fh分别为脉冲星自转频率估计范围的下限和上限，f
p
＝f
l
+(p-1)h，p＝1,2,
…
,n+1，p＝n+1，y
p
为f
p
对应的包含m个bin块的x射线脉冲星辐射信号的折叠轮廓，m＝2s,s∈n
+
，n
+
为正整数集；
[0013]
(3)获取x射线脉冲星辐射信号标准轮廓与折叠轮廓的时域互相关函数的频谱：
[0014]
以x射线脉冲星辐射信号的标准轮廓x为先验信息，对x、每个预设频率值f
p
所对应的x射线脉冲星辐射信号的折叠轮廓y
p
分别进行m点快速傅里叶变换，得到x的频谱f
x
、折叠轮廓的频谱的集合将f
x
和每个f
yp
的共轭结果相乘，得到x射线脉冲星辐射信号标准轮廓x与折叠轮廓的时域互相关函数的频谱的集合其中
[0015]
(4)获取每个预设频率对应的代价函数值：
[0016]
根据设定的x射线脉冲星辐射信号标准轮廓的频谱的截断点q，对每个时域互相关函数的频谱f
rp
中0点至q之间的频谱分量进行快速傅里叶逆变换，得到频谱截断的时域互相关函数的集合并将每个的最大值作为每个预设频率值f
p
对应的代价函数值，其中q∈n
+
,q∈[15,45]，τ∈n
+
,τ∈[1,m]；
[0017]
(5)获取x射线脉冲星自转频率的估计结果：
[0018]
对以p个预设频率值为自变量、以p个预设频率值所对应的p个代价函数值为因变量所组成的代价函数进行单峰值拟合，并将代价函数最大值所对应的频率值作为脉冲星自转频率的估计值
[0019][0020]
本发明与现有技术相比，具有如下优点：
[0021]
本发明由于选用了脉冲星辐射信号的标准轮廓作为先验信息，以周期折叠轮廓和标准轮廓的时域互相关函数的最大值作为代价函数值，不同于现有技术中仅利用周期折叠轮廓计算统计量的方式，弥补了在短时观测条件下观测数据信噪比低对统计量的计算所造成的负面影响，从而提高了在短时观测条件下脉冲星自转频率的估计精度。
附图说明
[0022]
图1为本发明的实现流程图。
[0023]
图2为本发明中对x射线脉冲星光子toa序列进行周期折叠的原理图。
[0024]
图3为本发明实施例中crab脉冲星辐射信号标准轮廓的频谱图。
[0025]
图4为本发明实施例中脉冲星自转频率估计值对应的代价函数图。
[0026]
图5为本发明实施例中脉冲星自转频率估计结果对应的周期折叠轮廓图。
[0027]
图6为本发明和现有技术对脉冲星实测数据的自转频率估计精度对比图。
[0028]
图7为根据本发明和现有技术对脉冲星实测数据的频率估计值对实测数据进行周期折叠所得的折叠轮廓与标准轮廓的相关系数对比图。
具体实施方式
[0029]
以下结合附图和具体实施例，对本发明作进一步详细描述。
[0030]
参照图1，本发明包括如下步骤：
[0031]
(1)初始化参数：
[0032]
初始化在观测时间[ta,tb]内航天器所探测的i个x射线脉冲星光子到达自身的toa序列为其中ti为第i个x射线脉冲星光子到达航天器的toa时刻，i＞10。
[0033]
本实例中，以crab脉冲星为例，设定crab脉冲星的标准频率为29.79113213hz，据
此通过计算机仿真获得一段观测时间为60s的脉冲星光子toa序列，对该仿真光子toa序列采用本发明进行脉冲星自转频率估计。
[0034]
(2)获取x射线脉冲星辐射信号的折叠轮廓：
[0035]
以h为步长对脉冲星自转频率估计范围[f
l
,fh]进行遍历，得到p个预设频率值通过用每个预设频率值f
p
对x射线脉冲星光子toa序列进行bin块数为m的周期折叠，得到对应的x射线脉冲星辐射信号的折叠轮廓的集合其中f
l
和fh分别为脉冲星自转频率估计范围的下限和上限，f
p
＝f
l
+(p-1)h，p＝1,2,
…
,n+1，p＝n+1，y
p
为f
p
对应的包含m个bin块的x射线脉冲星辐射信号的折叠轮廓，m＝2s,s∈n
+
，n
+
为正整数集。
[0036]
用预设频率值f
p
对x射线脉冲星光子toa序列进行周期折叠，周期折叠的原理图如图2所示，实现步骤为：
[0037]
(2a)将x射线脉冲星光子toa序列按照预设周期值t
p
分为k个子序列，其中t
p
＝1/f
p
，k＝(t
b-ta)/t
p
，如图2(a)所示；
[0038]
(2b)对第k个子序列进行m等分，如图2(b)所示，统计第k个子序列中第m个bin块内的x射线脉冲星光子数ck(tm)，其中m∈[1,m]，k∈[1,k]；
[0039]
(2c)将k个子序列中第m个bin块内的x射线脉冲星光子数ck(tm)相加，如图2(c)所示，得到一个对应于预设频率值f
p
的x射线脉冲星辐射信号的折叠轮廓y
p
，如图2(d)所示：
[0040][0041]
本实例中设定脉冲星自转频率的估计范围为[29.50,29.95]，遍历步长h为0.00001，由此获得一系列预设频率，bin块数m为1024。
[0042]
(3)获取x射线脉冲星辐射信号标准轮廓与折叠轮廓的时域互相关函数的频谱：
[0043]
根据信号处理知识，对于两个信号s1和s2，它们的互相关运算和卷积运算存在关系式：
[0044]rs1s2
(τ)＝s1(τ)*s2(-τ)
[0045]
对等式两边进行傅里叶变换，则有：
[0046][0047]
其中fr(jω)、f
s1
(jω)和分别为r
s1s2
(τ)、s1(τ)和s2(-τ)的频谱。而在时域直接进行互相关运算的时间复杂度高于频域，因此可以先在频域进行互相关运算，再将频域的运算结果经傅立叶逆变换得到时域的运算结果。
[0048]
以x射线脉冲星辐射信号的标准轮廓x为先验信息，对x、每个预设频率值f
p
所对应的x射线脉冲星辐射信号的折叠轮廓y
p
分别进行m点快速傅里叶变换，得到x的频谱f
x
、折叠轮廓的频谱的集合将f
x
和每个f
yp
的共轭结果相乘，得到x射线脉冲星辐
射信号标准轮廓x与折叠轮廓的时域互相关函数的频谱的集合其中
[0049]
本实例中选用crab脉冲星辐射信号的标准轮廓作为先验信息。
[0050]
(4)获取每个预设频率对应的代价函数值：
[0051]
由信号处理知识可知，信号的频谱分布主要集中在低频，而噪声的频谱分布在整个频域。脉冲星信号的频谱是幂率谱，如图3所示，其有用信息集中在低频部分，可以采取频率截断的方式截掉高频噪声成分。
[0052]
根据设定的x射线脉冲星辐射信号标准轮廓的频谱的截断点q，对每个时域互相关函数的频谱f
rp
中0点至q之间的频谱分量进行快速傅里叶逆变换，得到频谱截断的时域互相关函数的集合并将每个的最大值作为每个预设频率值f
p
对应的代价函数值，其中q∈n
+
,q∈[15,45]，τ∈n
+
,τ∈[1,m]。
[0053]
本实例中，设定频率截断位置q为35。
[0054]
(5)获取x射线脉冲星自转频率的估计结果：
[0055]
对以p个预设频率值为自变量、以p个预设频率值所对应的p个代价函数值为因变量所组成的代价函数进行单峰值拟合，并将代价函数最大值所对应的频率值作为脉冲星自转频率的估计值
[0056][0057]
本实例中所获得的代价函数如图4所示，横坐标为设定的脉冲星自转频率估计范围，纵坐标为每个预设频率值f
p
对应的代价函数值，此处对代价函数值完成了归一化。从图中可以看出代价函数在准确频率处具有尖锐的峰值，说明本发明对于预设频率范围内的频率变化较为敏感。最终经过代价函数拟合，得到脉冲星自转频率的最优估计值为29.79113194hz，与仿真所用的标准频率对比，该估计值具有较高的准确度。用该最优估计频率对仿真光子toa序列进行周期折叠，所得结果如图5所示，其尖锐清晰的折叠轮廓表明本发明频率估计的精度较高。
[0058]
以下结合卫星实测数据进行实验，对本发明的技术效果进行说明。
[0059]
1.实验条件与内容
[0060]
本实验选用罗西rxte卫星对crab脉冲星的观测数据，数据包编号为90802_02_09_00，实验所用的计算机参数为：cpu为i7-8700型、运行内存为16gb，计算机操作系统为windows 10版本，软件为matlab 2019版。
[0061]
实验内容一为：采用rxte卫星观测的crab脉冲星实测数据，对0.1s-1s时间范围内以0.1s为间隔、1s-10s时间范围内以1s为间隔和10s-60s时间范围内以10s为间隔的24个观测时长，每个观测时长截取50段观测数据，分别采用基于fft的功率谱积分方法和本发明进行脉冲星自转频率估计，以crab脉冲星计时模型外推的观测时段中点时刻频率作为参考值，计算不同观测时长下的频率估计结果的均方根误差，所得结果如图6所示。
[0062]
实验内容二为：采用rxte卫星观测的crab脉冲星实测数据，对0.1s-1s时间范围内以0.1s为间隔的10个观测时长，每个观测时长截取50段观测数据，分别采用基于fft的功率谱积分方法和本发明的频率估计值对观测数据进行周期折叠，计算折叠轮廓与标准轮廓的相关系数，所得结果如图7所示。
[0063]
2.实验结果分析
[0064]
参照图6，横坐标为卫星实测数据的观测时间长度，设定了从0.1s至60s之间的24个观测时长节点，纵坐标为每个观测时长下频率估计结果的均方根误差。从图中可以看到，随着观测时间的延长，两种频率估计方法的均方根误差都在逐步降低，尤其是在对观测时长为3s以内的卫星实测数据的频率估计中，本发明具有更低的均方根误差，说明在极短观测时间的条件下，本发明对脉冲星自转频率具有更高的估计精度。
[0065]
参照图7，横坐标为卫星实测数据的观测时间长度，设定了从0.1s至1s之间的10个观测时长节点，纵坐标为每个观测时长下周期折叠轮廓与标准轮廓的相关系数。从图中可以看到，在0.4s-1s观测时间内，以本发明的频率估计值对卫星实测数据进行周期折叠得到的轮廓与标准轮廓的相关系数更高，同样体现本发明在极短观测时间的条件下具有更高频率估计精度的特点。
[0066]
以上描述仅是本发明的一个具体实例，并未构成对本发明的任何限制，显然对于本领域的专业人员来说，在了解了本发明内容和原理后，都可能在不背离本发明原理、结构的情况下，进行任何形式和细节上的各种修改和改变，但是这些基于本发明思想的修正和改变仍在本发明的权利要求保护范围之内。