基于抗噪快速压缩感知的脉冲星到达时间估计方法及系统与流程

本发明属于航天器自主导航领域，特别涉及一种基于压缩感知的脉冲星到达时间估计方法及系统。

背景技术：

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

脉冲星到达时间是脉冲星导航系统能正常工作的基础。目前，脉冲星到达时间估计是脉冲星导航领域中的研究热点。近年来，已有学者开始将压缩感知应用于脉冲星信号处理。如：苏哲等人于2011年在《中国科学:物理学力学天文学》上发表了《基于压缩感知的脉冲星轮廓构建算法》；沈利荣等人于2016年在《Optik》上发表了《A robust compressed sensing based method for X-ray Pulsar Profile construction，鲁棒压缩感知方法在脉冲星信号重构中的应用》。这两种方法主要关注于脉冲星信号重构。其实，脉冲星信号去噪，重构的最终目的都是获得高精度脉冲星到达时间。黎胜亮等人于2014年在《Optik》上发表了《Fleet algorithm for X-ray pulsar profile construction and TOA solution based on compressed sensing，基于压缩感知的脉冲星轮廓与TOA快速算法》；余航等人于2015年在《aerospace science and technology》上发表了《A sparse representation-based optimization algorithm for measuring the time delay of pulsar integrated pulse profile，采用稀疏表示的脉冲星累积脉冲轮廓时间延迟测量方法》。这两篇文章将压缩感知应用于脉冲星到达时间估计，取得了较好的效果。但是，在测量矢量的选择上都是随机的，某些测量矢量可能无法发挥作用或作用很小。这样，大量无用的观测矢量被选用，而有用的观测矢量则被遗漏。

相比于传统压缩感知在其他领域的应用，当前脉冲星到达时间估计方法中的压缩感知有以下几个特点：

(1)字典中的各个元素不是正交的，而是冗余的。字典中的元素是拥有不同相位的同一脉冲轮廓。

(2)脉冲星轮廓信号被视为1阶稀疏信号。

为便于在航天器上实现脉冲星到达时间估计方法，需满足以下几个条件：

(1)幅度和相位是脉冲星轮廓的两个参数。由于脉冲星辐射强度不稳定，获得的脉冲星轮廓幅度也是未知的。脉冲星TOA估计重点是获得相位。因此，基于压缩感知的脉冲星TOA估计算法不应受幅度干扰。

(2)为了便于硬件实现，观测矩阵系数需为{1,-1}。

(3)为了减小计算量，应尽量减少观测矢量数目。

但是目前尚未有合适的技术方案出现。

技术实现要素：

本发明提出了一种基于抗噪快速压缩感知的脉冲星到达时间估计技术方案，旨在为导航系统快速地提供更加准确的脉冲到达时间估计。

为达到上述发明目的，本发明的技术方案是：

一种基于抗噪快速压缩感知的脉冲星到达时间估计方法，包括以下步骤，

步骤1，冗余字典的构建，实现如下，

设脉冲星累积轮廓的冗余字典表示为，

其中，n为序号，n＝{0,1,2,…N-1}，N为脉冲相位子间隔个数，冗余字典ψ为N×N的矩阵；

为字典中的第n个列矢量，是N×1的列矢量，表达式如下，

其中，s为脉冲星标准轮廓，i为变量，i＝{0,1,2,…N-1}，为的第i个元素，mod()表示取模；

步骤2，设定观测矩阵，实现如下，

设H为N阶哈达玛矩阵，

H＝[h0；h1；h2；…hn；…hN-1]

其中，hn是哈达玛矩阵中的行矢量，为1*N，n＝{0,1,2,…N-1}；

定义某个行矢量hn的观测能量en为该行矢量与字典各列矢量的乘积变化范围，en表达式如下，

en＝max(hn·Ψ)-min(hn·Ψ)

其中，max(hn·Ψ)表示行矢量hn与字典各列矢量的乘积的最大值，min(hn·Ψ)则表示行矢量hn与字典各列矢量的乘积的最小值；

根据预设的观测能量门限T，若某个行矢量hn的观测能量en大于T，该行矢量被选作观测矩阵的一行，否则舍弃；

设观测矩阵表示为Φ，为m×N的矩阵，m为哈达玛矩阵H中被选中的行数；

步骤3，观测矢量的获取，实现如下，

y＝Φ·x

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

步骤4，匹配追踪，实现如下，

先求出观测矢量y与各字典列矢量的相关性，记为n＝{0,1,2,…N-1}，上标T表示转置；

与观测值相关性最大的字典矢量对应序号nτ的表达式如下，

脉冲星到达时间估计值表达式为

而且，步骤3中，脉冲星累积轮廓x中第i个相位子间隔的表达式如下，

x(i)＝Poisson(A·t·λb/N+A·t·λs·s(mod(i-τ·N/p,N)))

其中，Poisson为泊松分布随机信号发生器，λb为背景噪声光子流量密度，λs为脉冲星辐射光子流量密度，A为X射线探测器有效面积，t为脉冲星辐射信号观测时间，τ为时间延迟，p为脉冲星辐射信号周期。

本发明还提供一种基于抗噪快速压缩感知的脉冲星到达时间估计系统，包括以下模块，第一模块，用于冗余字典的构建，实现如下，

设脉冲星累积轮廓的冗余字典表示为，

其中，n为序号，n＝{0,1,2,…N-1}，N为脉冲相位子间隔个数，冗余字典ψ为N×N的矩阵；

为字典中的第n个列矢量，是N×1的列矢量，表达式如下，

其中，s为脉冲星标准轮廓，i为变量，i＝{0,1,2,…N-1}，为的第i个元素，mod()表示取模；

第二模块，用于设定观测矩阵，实现如下，

设H为N阶哈达玛矩阵，

H＝[h0；h1；h2；…hn；…hN-1]

其中，hn是哈达玛矩阵中的行矢量，为1*N，n＝{0,1,2,…N-1}；

定义某个行矢量hn的观测能量en为该行矢量与字典各列矢量的乘积变化范围，en表达式如下，

en＝max(hn·Ψ)-min(hn·Ψ)

其中，max(hn·Ψ)表示行矢量hn与字典各列矢量的乘积的最大值，min(hn·Ψ)则表示行矢量hn与字典各列矢量的乘积的最小值；

根据预设的观测能量门限T，若某个行矢量hn的观测能量en大于T，该行矢量被选作观测矩阵的一行，否则舍弃；

设观测矩阵表示为Φ，为m×N的矩阵，m为哈达玛矩阵H中被选中的行数；

第三模块，用于观测矢量的获取，实现如下，

y＝Φ·x

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

第四模块，用于匹配追踪，实现如下，

先求出观测矢量y与各字典列矢量的相关性，记为n＝{0,1,2,…N-1}，上标T表示转置；

与观测值相关性最大的字典矢量对应序号nτ的表达式如下，

脉冲星到达时间估计值表达式为

而且，第三模块中，脉冲星累积轮廓x中第i个相位子间隔的表达式如下，

x(i)＝Poisson(A·t·λb/N+A·t·λs·s(mod(i-τ·N/p,N)))

其中，Poisson为泊松分布随机信号发生器，λb为背景噪声光子流量密度，λs为脉冲星辐射光子流量密度，A为X射线探测器有效面积，t为脉冲星辐射信号观测时间，τ为时间延迟，p为脉冲星辐射信号周期。

本发明研究了一种低信噪比条件下，基于抗噪快速压缩感知的脉冲星到达时间估计技术方案，在压缩感知中，不同相位的脉冲轮廓构成冗余的稀疏字典；以哈达玛矩阵为基础，以观测能量为判据，选择少量哈达玛矢量构造观测矩阵，该观测矩阵行数较少，抗噪能力强；最后，利用匹配追踪估计脉冲到达时间。与传统基于压缩感知的脉冲星到达时间估计方法相比，该技术方案对噪声具有较强的鲁棒性，观测矢量数目较少，并且便于硬件实现，处于国际领先地位，具有重要的市场价值。

附图说明

图1为本发明实施例的流程图；

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

图3为本发明实施例的观测能量示意图；

图4为本发明实施例的相关性示意图。

具体实施方式

以下结合附图和实施例对本发明技术方案进行具体描述。

现有技术中，大量无用的观测矢量被选用，而有用的观测矢量则被遗漏。针对这一问题，本发明定义了观测能量，并依据此选择有用的测量矢量，并构成观测矩阵。

实施例中，以X射线脉冲星PSR B1937+21为例对本发明作进一步说明。脉冲轮廓数据可从欧洲脉冲星网站(The European Pulsar Network)下载。将该数据归一化后即为标准脉冲轮廓。图2给出了PSR B1937+21的标准脉冲轮廓。

实施例提供的流程包括以下步骤：

步骤1：冗余字典的构建。

脉冲星PSR B1937+21累积轮廓的冗余字典可表示为：

其中，n为序号，n＝{0,1,2,…N-1}，N为脉冲相位子间隔个数，在本实施例中，N＝1024。冗余字典ψ为N×N的矩阵。

为字典中的第n个列矢量，是N×1的列矢量，其表达式为：

其中，s为脉冲星标准轮廓，其波形如图2所示，波形横坐标为相位，纵坐标为光子密度。i为变量，i＝{0,1,2,…N-1}，为的第i个元素，即字典ψ中第n列第i行的元素。mod()表示取模。

步骤2：观测矩阵的设定。

为了便于硬件实现，观测矩阵系数需为{1,-1}。哈达玛(Hadamard)矩阵是由+1和-1元素构成的正交方阵。本发明提出以哈达玛矩阵为基础，设计观测矩阵。设H为N阶哈达玛矩阵。在本实施例中，H为1024阶哈达玛矩阵。

H＝[h0；h1；h2；…hn；…hN-1] (3)

其中，hn是哈达玛矩阵中的行矢量，为1*N，n＝{0,1,2,…N-1}。

定义某个行矢量hn的观测能量en为该行矢量与字典各列矢量的乘积变化范围，en表达式如下：

en＝max(hn·Ψ)-min(hn·Ψ) (4)

该行矢量hn与字典ψ各列矢量乘积记为hn·Ψ，具体包括hn与N个列矢量(n＝{0,1,2,…N-1})分别取乘积，共得到N个元素，max(hn·Ψ)表示这些元素中的最大值，即行矢量hn与字典各列矢量的乘积的最大值；min(hn·Ψ)则表示最小值，即表示行矢量hn与字典各列矢量的乘积的最小值。二者差值体现了变化幅度。若该值不变化，或变化幅度很小，则根据该观测值无法区分字典中的元素，即该值对信号恢复无作用。因此，观测能量越大，hn越能在信号恢复中发挥作用，越能抗噪。哈达玛矩阵中的各个行矢量对应的观测能量如图3所示。

定义观测能量门限T，根据门限提取高观测能量的哈达玛行向量。若观测能量大于该门限，hn可选作观测矩阵的一行，否则舍弃。设观测矩阵可表示为Φ，为m×N的矩阵，m为哈达玛矩阵被选中的行数。一般而言，门限T越小，行数m越多，精度越高，但计算量增大；反之亦然。

在本实施例中，若设定门限T为1.9，则m为6，即选择了哈达玛矩阵的第128，192，256，320，384，448行；若设定门限T为1.5，则m为14；若设定门限T为0.5，则m为62；若设定门限为0.2，则m为116；若设定门T限为0，则m为1024。在本实施例中，经过仿真实验，在计算量与脉冲星到达时间估计精度之间折衷，门限设定为0.2。具体实施时，本领域技术人员可自行根据需要预设观测能量门限T的取值。

步骤3：观测矢量的获取。

观测矢量的计算式如下：

y＝Φ·x (5)

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

由于光子的到达时间服从泊松分布，x中第i个相位子间隔的表达式如下：

x(i)＝Poisson(A·t·λb/N+A·t·λs·s(mod(i-τ·N/p,N))) (6)

其中，Poisson为泊松分布随机信号发生器，λb为背景噪声光子流量密度，λs为脉冲星辐射光子流量密度，A为X射线探测器有效面积，t为脉冲星辐射信号观测时间，τ为时间延迟，p为脉冲星辐射信号周期。

在本实施例中，λb为0.005ph/cm²/s，λs为4.99*10^-5ph/cm²/s，A为384cm²，t为1000s，p为1.5578ms，τ为0.38945ms。

步骤4：匹配追踪的实现。

由于信号稀疏阶数为1，本发明匹配追踪算法无需迭代。可先求出观测矢量y与各字典列矢量的相关性，记为n＝{0,1,2,…N-1}，上标T表示转置。不同字典矢量与观测值的相关性如图4所示。与观测值相关性最大的字典矢量对应序号nτ即为所求。综上，nτ的表达式如下：

由于匹配方法与信号幅度值无关，所以，脉冲星轮廓幅度不会对相位估计产生影响。

脉冲星到达时间估计值表达式如下：

综上，该方法满足了应用于航天器的三个条件，具有较好的性能。

对本发明实施例技术方案进行了500次蒙特卡洛仿真，其误差为1.4669us，对应的定位误差为440m。

具体实施时，本发明所提供方法可基于软件技术实现自动运行流程，也可采用模块化方式实现相应系统。本发明实施例相应提供一种基于抗噪快速压缩感知的脉冲星到达时间估计系统，包括以下模块，

第一模块，用于冗余字典的构建，实现如下，

设脉冲星累积轮廓的冗余字典表示为，

其中，n为序号，n＝{0,1,2,…N-1}，N为脉冲相位子间隔个数，冗余字典ψ为N×N的矩阵；

为字典中的第n个列矢量，是N×1的列矢量，表达式如下，

其中，s为脉冲星标准轮廓，i为变量，i＝{0,1,2,…N-1}，为的第i个元素，mod()表示取模；

第二模块，用于设定观测矩阵，实现如下，

设H为N阶哈达玛矩阵，

H＝[h0；h1；h2；…hn；…hN-1]

其中，hn是哈达玛矩阵中的行矢量，为1*N，n＝{0,1,2,…N-1}；

定义某个行矢量hn的观测能量en为该行矢量与字典各列矢量的乘积变化范围，en表达式如下，

en＝max(hn·Ψ)-min(hn·Ψ)

其中，max(hn·Ψ)表示行矢量hn与字典各列矢量的乘积的最大值，min(hn·Ψ)则表示行矢量hn与字典各列矢量的乘积的最小值；

根据预设的观测能量门限T，若某个行矢量hn的观测能量en大于T，该行矢量被选作观测矩阵的一行，否则舍弃；

设观测矩阵表示为Φ，为m×N的矩阵，m为哈达玛矩阵H中被选中的行数；

第三模块，用于观测矢量的获取，实现如下，

y＝Φ·x

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

第四模块，用于匹配追踪，实现如下，

先求出观测矢量y与各字典列矢量的相关性，记为n＝{0,1,2,…N-1}，上标T表示转置；

与观测值相关性最大的字典矢量对应序号nτ的表达式如下，

脉冲星到达时间估计值表达式为

各模块具体实现可参见相应步骤，本发明不予赘述。

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