CN109186612A

CN109186612A - 基于压缩感知和克拉美劳界的快速脉冲星周期估计方法

Info

Publication number: CN109186612A
Application number: CN201811036881.7A
Authority: CN
Inventors: 刘劲; 刘兰兰; 李阳阳; 何宇飞; 康志伟
Original assignee: Wuhan University of Science and Engineering WUSE
Current assignee: Wuhan University of Science and Engineering WUSE
Priority date: 2018-09-06
Filing date: 2018-09-06
Publication date: 2019-01-11
Anticipated expiration: 2038-09-06
Also published as: CN109186612B

Abstract

一种基于压缩感知和克拉美劳界的快速脉冲星周期估计方法，包括计算脉冲星周期估计误差的克拉美劳界，进行冗余字典的构建，利用傅立叶变换的低频部分设置观测矩阵，获取量测矢量；实现平滑匹配，包括基于感知矢量和量测矢量得到匹配矢量，构成匹配度矩阵，利用平滑模板对匹配度矩阵进行均值滤波，基于平滑后的匹配度矩阵求取畸变展宽度粗估计，利用超解析度估计求出更精确的畸变展宽度估计，计算得到相应的频偏和脉冲星周期误差。该方法精度高，接近于克拉美劳界。由于避免了脉冲星辐射光子信号的反复折叠，本发明计算速度快。

Description

基于压缩感知和克拉美劳界的快速脉冲星周期估计方法

技术领域

本发明属于航天器自主导航领域，特别涉及一种基于压缩感知和CRB的快速X射线脉冲星周期估计方法。

背景技术

X射线脉冲星导航是一种新兴的航天器自主导航方法。X射线脉冲星导航是利用脉冲星辐射信号进行导航。X射线脉冲星对外辐射稳定的、周期性的脉冲信号。将这些信号按照周期进行累积，可获得稳定的脉冲轮廓。将该轮廓与标准脉冲轮廓进行比较，就能获得脉冲星导航的量测量——脉冲到达时间(TOA，time-of-arrival)。

精确的脉冲星周期是获取高精度脉冲累积轮廓的前提。若脉冲星周期存在较大误差，脉冲星累积轮廓会出现明显“畸变”，进而影响脉冲到达时间估计精度。事实上，受周期跃变等因素的影响，脉冲星周期并不稳定。近年来，脉冲星周期估计已经称为一个研究热点。李建勋等人在《物理学报》上分别于2010年2012年发表了《基于循环平稳信号相干统计量的脉冲星周期估计新方法》和《一种脉冲星周期的时域估计新方法》，周庆勇等人于2013年在《物理学报》上发表了《非等间隔计时数据的X射线脉冲星周期快速搜素算法》，张新源等人于2015年在《宇航学报》上发表了《脉冲星导航轮廓折叠失真与周期估计算法》；以上这些方法都需要进行多次周期试探折叠。但是，X射线脉冲星辐射数据量极大，10⁷量级以上，这导致计算量极大，难以在轨实现。宋佳凝等人于2017年在《光子学报》上发表了《基于TOA信息的脉冲星信号周期估计方法》。该方法将一组观测脉冲光子到达时间序列进行等时间间隔分段，在时域上对每段时间间隔序列进行相位估计，获得一系列TOA，并以TOA与时间图像的斜率作为依据，估计脉冲星周期。这种方法本质上是短时TOA估计，根据一组TOA解算脉冲星周期。在X射线探测器面积小的情况下，该方法在短时间内无法获得有效的脉冲星累积轮廓。受载重限制，航天器上X射线探测器面积不大。

综上，若要实现快速高精度的脉冲星周期估计方法，必须避免脉冲星光子信号的折叠。

发明内容

为了克服现有技术的缺陷，本发明提出一种基于压缩感知和CRB的快速X射线脉冲星周期估计方法，通过检测脉冲星累积轮廓的畸变，从而避免脉冲星辐射光子信号的多次折叠。

本发明技术方案提供一种基于压缩感知和克拉美劳界的快速脉冲星周期估计方法，其特征在于，包括以下步骤：

步骤1，计算脉冲星周期估计误差的克拉美劳界，包括确定脉冲星频率f的估计的方差的下限，确定脉冲星周期T的估计的方差的下限；

步骤2，进行冗余字典的构建，实现如下，

设脉冲星畸变轮廓的冗余字典表示为，

其中，m为畸变展宽度，m＝{0，1,2,…M-1}，M为最大畸变展宽度；

为字典中的元素，是N×1的列矢量，N为脉冲星轮廓间隔数，表达式为

其中，表示循环互相关，G_m表示扩展矢量，H为脉冲星标准离散轮廓；

采用超分辨率估计方式，将扩展函数宽度以间隔的整数倍设置，设G_m(n)表示扩展矢量中第n个元素，n为变量，n＝{0，1,2,…N-1}，取值如下，

元素对应的频偏Δf_m为，

其中，T_obs为脉冲星观测时间；

步骤3，利用傅立叶变换的低频部分设置观测矩阵Φ，包括从大小为N×N的傅立叶变换矩阵W_N中取前L行和后L-1行构成观测矩阵Φ，其中，L为预设的参数

步骤4，获取量测矢量如下，

y＝Φ·x

其中，y为量测矢量，为L×1矢量；x为脉冲累积轮廓，为1×N矢量；

步骤5，实现平滑匹配，包括基于感知矢量a_m和量测矢量y得到匹配矢量z_m，构成匹配度矩阵S，利用平滑模板对匹配度矩阵进行均值滤波，基于平滑后的匹配度矩阵S′求取畸变展宽度粗估计利用超解析度估计求出更精确的畸变展宽度估计计算得到相应的频偏Δf和脉冲星周期误差ΔT。

而且，步骤1的实现如下，

确定脉冲星频率f的估计的方差的下限，

确定脉冲星周期T的估计的方差的下限，

其中，Δf为脉冲星频率误差，h()为标准脉冲轮廓，为相位。

而且，步骤5中，基于感知矢量a_m和量测矢量y得到匹配矢量z_m，实现如下，感知矢量a_m可表示为：

匹配矢量z_m为量测矢量与感知矩阵的乘积，表示为，

z_m(k)＝y(k)×a_m(k)，k＝{0,1,2,…L-1}

其中，k为变量，z_m(k)为匹配矢量z_m的第k个元素，y(k)为量测矢量y的第k个元素，a_m(k)为感知矢量a_m的第k个元素；

补零，使得z_m长度为N，

z_m(k)＝0，k＝{L,L+1,…,N-1}。

而且，步骤5中，构成匹配度矩阵S的实现方式为，该矩阵第m列S_m表示为，

其中，表示傅立叶拟变换；||表示取幅度。

而且，步骤5中，基于平滑后的匹配度矩阵S′，畸变展宽度粗估计由下式求出，

其中，表示矩阵S′最大值所在的行，S′(n,m)为平滑后的匹配度矩阵的第n行第m列的元素；

在此基础上，利用超解析度估计，求出更精确的畸变展宽度估计

其中，为平滑后的匹配度矩阵的第行第列的元素，为平滑后的匹配度矩阵的第行第列的元素，为平滑后的匹配度矩阵的第行第列的元素。

而且，步骤5中，计算得到相应的频偏Δf和脉冲星周期误差ΔT如下，频偏Δf为，

脉冲星周期误差为，

ΔT＝Δf·T²+T_B。

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

(1)脉冲星周期估计精度高。脉冲星周期估计误差接近于克拉美劳界，即接近于理论极限。因此，脉冲星周期估计精度高。

(2)估计算法的速度快。该方法避免了脉冲星光子信号的折叠，而是检测脉冲星累积轮廓的畸变。因此，本发明的计算量小，可满足在轨计算实时性的要求。

(3)本发明所需X射线探测器面积小。由于本发明无需解算短时TOA，即使X射线探测器面积较小，在一段脉冲星观测时间内也能累积出一个脉冲轮廓。

附图说明

图1是本发明实施例的脉冲星到达时间估计流程图；

图2是本发明实施例的PSR B0531+21标准脉冲轮廓示意图；

图3是本发明实施例的脉冲星周期克拉美劳界示意图。

具体实施方式

以下结合附图和实施例具体说明本发明的技术方案。

实施例以X射线脉冲星PSR B0531+21为例。脉冲轮廓数据可从欧洲脉冲星网站(The European Pulsar Network)下载。将该数据归一化并插值，即可得到标准脉冲轮廓。图2给出了PSR 0531+21的标准脉冲轮廓。

在本实施例中，X射线探测器面积A为400cm²，时间分辨率为10^-6s。脉冲星周期T为0.033s，脉冲星频率为f＝1/T。脉冲星轮廓间隔数N为33000。脉冲星观测时间T_obs为1000s，在此期间内脉冲星周期数K约为30303个。X射线脉冲星光子流量密度λ_s为1.54ph/cm²/s，X射线背景流量密度λ_b为0.005ph/cm²/s。

步骤1：计算脉冲星周期估计误差的CRB(Cramer-Rao bound，克拉美劳界)，包括确定脉冲星频率f的估计的方差的下限，确定脉冲星周期T的估计的方差的下限。

脉冲星频率f的估计的方差需满足式(1)。

其中，Δf为脉冲星频率误差，h( )为标准脉冲轮廓，为相位。

脉冲星周期T的估计的方差则需满足式(3)。

图3给出了不同频偏下的脉冲星周期估计误差的CRB。从图3中可以看出，频偏对CRB影响不大。

步骤2：以克拉美劳界作为参考，构建脉冲星畸变冗余字典。

脉冲星PSR B0531+21畸变轮廓的冗余字典可表示为：

其中，m为畸变展宽度，m＝{0，1,2,…M-1}，M为最大畸变展宽度，在本实施例中，M＝11。

冗余字典ψ为N×M的矩阵。

为字典中的元素，是N×1的列矢量，

其表达式为：

其中，表示循环互相关，G_m表示扩展矢量。H为脉冲星标准离散轮廓。本实施例中，X射线脉冲星周期误差CRB在10^-11s量级，观测时间为1000s，扩展效果约为0.3个间隔。

考虑到超分辨率估计方法可将精度提高一个数量级，本实施例中扩展函数宽度以间隔的整数倍来设计，设G_m(n)表示扩展矢量中第n个元素，n为变量，n＝{0，1,2,…N-1}，具体取值如下：

字典元素对应的频偏Δf_m为：

步骤3：观测矩阵的设计：以傅立叶基的低频部分作为观测矩阵。

脉冲星累积轮廓的相位存在误差，而傅立叶变换的幅度信息与相位无关。因此，本发明采用傅立叶变换基设计观测矩阵。由于低频段信噪比较高，本发明仅利用傅立叶变换的低频部分。设W_N为大小为N×N的傅立叶变换矩阵。

其中，e为自然底数。

在本实施例中，W_N为33000×33000傅里叶变换矩阵。

观测矩阵Φ取W_N的前L行和后L-1行。其中，L为预设的参数。具体实施时，可根据仿真实验，取最高精度的周期估计对应的L值。本实施例中，L为1300。

步骤4：量测矢量的获取。

由于本方法只能确定频偏的绝对值，而无法确定其正负，本发明缩短累积周期。本实施例中，脉冲星周期T为0.033s，周期误差为3.4*10-¹¹s，周期偏置T_B为3.3*10-¹⁰s。即累积脉冲星辐射光子信号按照T-T_B累积。

量测矢量的计算式如下：

y＝Φ·x (10)

其中，y为量测矢量，为L×1矢量。x为脉冲累积轮廓，为1×N矢量。

步骤5：平滑匹配的实现。

由于信号稀疏阶数为1，本发明匹配算法无需跟踪。为抑制噪声，利用平滑模板对匹配度矩阵进行均值滤波。具体如下：

感知矢量a_m可表示为：

匹配矢量z_m为量测矢量与感知矩阵的乘积，可表示为：

z_m(k)＝y(k)×a_m(k)，k＝{0，1,2,…L-1} (12)

其中，k为变量，z_m(k)为匹配矢量z_m的第k个元素，y(k)为量测矢量y的第k个元素，a_m(k)为感知矢量a_m的第k个元素。

补零，使得z_m长度为N。

z_m(k)＝0，k＝{L,L+1,…,N-1} (13)

匹配度矩阵为S，该矩阵第m列S_m表示为：

其中，表示傅立叶拟变换。||表示取幅度。

利用平滑模板M_s对匹配度矩阵进行均值滤波。M_s的表达式如下：

平滑后的匹配度矩阵为S′。

畸变展宽度粗估计可由下式求出:

其中，表示矩阵S′最大值所在的行，S′(n,m)为平滑后的匹配度矩阵的第n行第m列的元素。

在此基础上，利用超解析度估计方法，可求出更精确的畸变展宽度估计

频偏Δf为：

脉冲星周期误差为：

ΔT＝Δf·T²+T_B (19)

本发明实施例对其进行了500次蒙特卡洛仿真，其误差为2.1×10^-11s，CRB为6.5×10^-12s。

具体实施时，可以采用软件的方式实现自动运行以上流程。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims

1.一种基于压缩感知和克拉美劳界的快速脉冲星周期估计方法，其特征在于，包括以下步骤：

步骤2，进行冗余字典的构建，实现如下，

设脉冲星畸变轮廓的冗余字典表示为，

元素对应的频偏Δf_m为，

其中，T_obs为脉冲星观测时间；

步骤4，获取量测矢量如下，

y＝Φ·x

2.根据权利要求1所述基于压缩感知和克拉美劳界的快速脉冲星周期估计方法，其特征在于：

步骤1的实现如下，

确定脉冲星频率f的估计的方差的下限，

确定脉冲星周期T的估计的方差的下限，

其中，Δf为脉冲星频率误差，h()为标准脉冲轮廓，为相位。

3.根据权利要求1或2所述基于压缩感知和克拉美劳界的快速脉冲星周期估计方法，其特征在于：步骤5中，基于感知矢量a_m和量测矢量y得到匹配矢量z_m，实现如下，

感知矢量a_m可表示为：

匹配矢量z_m为量测矢量与感知矩阵的乘积，表示为，

z_m(k)＝y(k)×a_m(k)，k＝{0,1,2,…L-1}

补零，使得z_m长度为N，

z_m(k)＝0，k＝{L,L+1,…,N-1}。

4.根据权利要求3所述基于压缩感知和克拉美劳界的快速脉冲星周期估计方法，其特征在于：步骤5中，构成匹配度矩阵S的实现方式为，该矩阵第m列S_m表示为，

其中，表示傅立叶拟变换；||表示取幅度。

5.根据权利要求4所述基于压缩感知和克拉美劳界的快速脉冲星周期估计方法，其特征在于：

步骤5中，基于平滑后的匹配度矩阵S′，畸变展宽度粗估计由下式求出，

6.根据权利要求4所述基于压缩感知和克拉美劳界的快速脉冲星周期估计方法，其特征在于：步骤5中，计算得到相应的频偏Δf和脉冲星周期误差ΔT如下，

频偏Δf为，

脉冲星周期误差为，

ΔT＝Δf·T²+T_B。