基于EMD-CS的脉冲星周期超快速估计

2020-09-10刘劲韩雪侠宁晓琳陈晓康志伟

航空学报 2020年8期

关键词：脉冲星畸变矢量

刘劲，韩雪侠，宁晓琳，陈晓，康志伟

1. 武汉科技大学信息科学与工程学院，武汉 430081

2. 北京航空航天大学仪器科学与光电工程学院，北京 100191

3. 上海卫星工程研究所，上海 200240

4. 湖南大学信息科学与工程学院，长沙 410082

脉冲到达时间(Time-of-Arrival, TOA)是X射线脉冲星导航的基本观测量[1-2]。航天器上的X射线探测器收集X射线脉冲星辐射的光子[3]，经历元折叠形成累积脉冲轮廓[4-5]，再结合标准脉冲轮廓、脉冲星计时模型等信息估计脉冲到达时间。但航天器运动和脉冲星周期跃变[6]等导致接收到的脉冲信号周期发生变化，使得按固有周期累积的脉冲轮廓发生畸变，进而导致脉冲TOA发生漂移。脉冲周期估计对脉冲TOA高精度估计具有重要意义，目前已成为一个研究热点。

脉冲星周期估计方法可分为2类：一类利用脉冲星TOA漂移估计周期误差；另一类根据脉冲轮廓畸变反演周期误差，是目前的主流方法。

基于脉冲星TOA漂移的周期误差估计的基本思路如下：建立周期误差与脉冲星TOA漂移之间的关系模型。在此基础上，利用最小二乘法估计TOA与周期误差[7-8]；利用天文信息补偿[9]或优化卡尔曼滤波器，如：EKF-CMBSEE(Extended Kalman Filter-Correlated Measurement Bias and State Estimation Error)[10]、跟踪滤波器[11]，以达到抑制周期误差影响的目的。

基于轮廓畸变的脉冲星周期估计方法的基本思路如下：尝试按不同周期折叠脉冲星信号，得到一系列畸变脉冲累积轮廓，然后找出畸变度最小的脉冲累积轮廓，其对应的周期就是固有周期。如：时域χ2法、快速折叠法[12]、傅立叶频谱法[13]、循环平稳信号相干统计量法[14]、最大相关方差搜索法[15]、基于Lomb的χ2算法[16]、脉冲轮廓熵算法[17]和轮廓特征法[18]等。上述方法无一例外都需要尝试按不同周期折叠脉冲星信号。但是，脉冲星辐射信号数据量庞大，导致计算量大等问题，不适合星载计算机在轨运行。

压缩感知(Compressed Sensing，CS)[19-20]是一种新兴的信号处理方法，对稀疏信号的处理能力很强，恰巧脉冲星信号是典型的稀疏信号。近年来，基于CS的脉冲星信号处理成为一个研究热点[21-23]。基于快速傅里叶变换-压缩感知(FFT-CS)的脉冲星周期快速估计方法[24]利用CS压缩脉冲星辐射信号，再直接根据脉冲星轮廓畸变度估计脉冲固有周期。该方法避免了多次脉冲星信号折叠，大幅减少了计算量，使得脉冲星周期在轨估计成为可能。

但是，傅立叶变换将信号能量分散在数量庞大的傅立叶级数中。FFT-CS采用低频傅立叶变换，所需的傅立叶级数仍然上千，这导致测量矩阵尺寸极大。为提高计算速度，必须大幅减小尺寸。经验模态分解(Empirical Mode Decomposition, EMD)[25-26]根据原始信号的固有特征，自适应地将其分解成有限个固有模态函数(Intrinsic Mode Function，IMF)。IMF包含了原始信号不同时间尺度的局部特征信号。而脉冲轮廓畸变度这一微弱局部特征可体现在某些IMF中。与傅立叶变换相比，IMF数量大幅减少，仅为101～102量级。因此，本文将IMF构成的测量矩阵替换FFT-CS中的低频傅立叶变换矩阵，提出了一种基于EMD-CS的脉冲星周期超快速估计方法，旨在保持估计精度的同时，减小计算量，加快运行速度。

1 最优IMF测量矩阵

为减小测量矩阵尺寸，本文利用EMD分解得到的IMF构造测量矩阵来替代低频傅立叶矩阵，并在此基础上提出了一种EMD-CS方法。EMD-CS方法的基本流程与FFT-CS方法类似，二者唯一不同在于测量矩阵。本节重点研究测量矩阵。

在CS中，测量矢量的获取是通过测量矩阵来实现的。通常，测量矩阵尺寸越小，计算量越少。因此寻找一种测量矩阵，用更少的测量次数达到预期的结果变得尤为重要。

将这一系列IMF序列组合而成一个IMF原始测量矩阵Φ0(MIMF×N)，N为一个脉冲周期内的脉冲间隔数；MIMF为不同畸变轮廓分解的IMF序列数总和,可表示为

(1)

Φ0可以表示为

(2)

(3)

(4)

(5)

MEIMF的值为

(6)

由式(1)和式(6) 可看出，ΦE的尺寸小于Φ0。

迭代剔除法在地面进行，筛选得到的最优IMF测量矩阵存入器载计算机，无需在轨筛选。

2 EMD-CS

本节将介绍EMD-CS方法估计脉冲星周期的整个算法流程，如图1所示分成4个模块：① 最优IMF测量矩阵的设计；② 畸变字典的构造；③ 脉冲频率初始偏置；④ 最大元素的超分辨率稀疏恢复。其中，最优IMF测量矩阵已经在第1节中描述。下面将介绍其他3个模块：

图1 所提算法流程图

1) 畸变字典的构造

文献[24]已证明，畸变脉冲轮廓可视为标准脉冲轮廓与门函数的循环互相关。因此，构造了MD个不同畸变度的脉冲轮廓hm，表达式为

hm=Gm⊗h

(7)

式中：Gm(N×1维)为扩散矢量；m为畸变脉冲轮廓序号；h(N×1维)为标准脉冲轮廓；⊗表示循环互相关。Gm满足：

(8)

从式(7)和式(8)可看出，Gm相当于MD个宽度不同的矩形窗。脉冲星轮廓畸变度定义为畸变宽度与一个脉冲周期内的周期间隔数N之比，即m/N。随着m增大，矩形窗变宽，脉冲累积轮廓的畸变度变大。

脉冲星导航常用于深空探测巡航段，此时的轨道模型近似于线性的[10]，所以可把脉冲星周期误差视为非时变的。式(8)所建立的模型适用于这种情况。

将第m个畸变脉冲轮廓循环移位构成第m个畸变脉冲轮廓子字典φm(N×P维)，表达式为

φm={φmi(t)|φmi(t)=hm(t±i)}

i=0,1,…,(P-1)/2

(9)

式中:P为最大相位偏移量；φmi(t)为子字典φm中的第i个原子；hm(t)为第m个畸变脉冲轮廓。可将MD个子字典φm合成一个三维矩阵字典Ψ(N×P×MD维)。

2) 脉冲频率初始偏置

累积脉冲轮廓是根据脉冲固有周期T0将脉冲星辐射光子进行折叠而形成的，这一过程称为历元折叠。但是，多普勒效应以及脉冲星周期跃变会使航天器接收到的脉冲周期存在误差ΔT，与频率偏移之间的关系可表示为

(10)

式中：Δf为频率偏移；f0为固有频率。

文献[24]已经证明，不能直接根据脉冲星轮廓畸变辨别Δf的正负号。为了解决此问题，引入脉冲星周期偏移Toffset，且满足Toffset≫ΔT。这样，-Toffset+ΔT和-Toffset-ΔT的符号都是负号，但幅度不同。所以，脉冲轮廓按照周期T0-Toffset+ΔT进行历元折叠形成累积脉冲轮廓。

3) 基于最大值的超分辨率稀疏恢复

在压缩感知中，信号稀疏恢复是通过感知矩阵与测量矢量相匹配，根据最大匹配值实现的。感知矩阵T是最优IMF测量矩阵与字典的乘积：

T=ΦEΨ

(11)

式中：T的尺寸为MEIMF×P×MD。

测量矢量y(MEIMF×1维)的获得方式为

(12)

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

S(i,j)=T(:,i,j)Ty

i=1,2，…，P;

j=1,2,…,MD

(13)

(14)

为了使估计更加准确，运用基于最大值的超分辨率稀疏恢复：

(15)

(16)

(17)

3 复杂度分析

本节分析EMD-CS方法的计算复杂度。在该方法中，脉冲星轮廓累积、测量矢量获取和匹配在航天器上运行。其中，脉冲星累积过程可在脉冲星信号的观测时间内进行。而测量矢量获取和匹配是在观测结束后才进行的。考虑到航天器高速飞行的影响，需尽快完成这2个过程，才能保证算法的实时性。因此，本文分析二者计算复杂度即可。

测量矢量的获取由式(12)确定。最优IMF测量矩阵尺寸为MEIMF×N，脉冲累积轮廓为长度为N。这需要MEIMF×(2N-1)MAC，MAC(Multiply/Accumulate Computation)表示加法乘法的计算次数。

匹配由式(13)确定。测量矢量与一个原子匹配，所需计算量为2MEIMF-1 MAC。字典中原子个数为MD×P。测量矢量与字典匹配所需计算量是(2MEIMF-1)×MD×PMAC。

总计算量为MEIMF×(2N-1)+(2MEIMF-1)×MD×P≈2MEIMF×(N+MD×P) MAC。

以上可以看出，计算量与MEIMF呈正比关系。若能大幅减小MEIMF，计算量也将大幅减少。此外，与FFT-CS相比，EMD-CS仅有实数部分，计算量将更小。

4 仿真实验

将基于EMD-CS的脉冲星周期估计方法与基于FFT-CS的周期估计进行仿真比较。首先，筛选出最优的测量矩阵，并分析其有限等距性质；其次，将其与基于FFT-CS的脉冲星周期估计方法作比较，体现本文方法的优越性；再次，分析脉冲星和背景噪声的光子流量密度、X射线探测器面积以及观测时间对脉冲星周期估计精度的影响。

4.1 仿真条件

脉冲星数据来自于欧洲脉冲星网(EPN)，仿真实验在处理器为英特尔Core i5-7500@3.40 GHz四核、内存为8 GB的计算机上运行。最大畸变度MD=21，即构造21个脉冲星畸变轮廓，最大相位偏移量P=21。仿真条件如表1所示。

表1 PSR B0531+21仿真条件

4.2 最优测量矩阵

表2给出了剔除某些IMF的脉冲星周期估计误差，表的首行和首列表示剔除的IMF序号。从表2可以看出，若是剔除第1个IMF，脉冲星周期估计误差大，所以必须保留第1个IMF；而同时剔除残差和第5个IMF时，误差最小(70.177 8 ps)。

表2 第1次剔除IMF的误差

本文保留第1个IMF，剔除残差和第5个IMF后，在剩下的IMF中继续实施剔除操作，仿真结果如表3所示。从表3中可以看出，剔除第6和第7个IMF，误差最小(62.224 8 ps)。表4给出了剔除第5、6、7个IMF以及残差之后的周期估计。实验结果表明，再剔除第3和第8个IMF，估计误差更小(57.063 6 ps)。表5给出了在此基础上，4次剔除之后的结果。此次剔除，使得误差增大。因此，本文将IMF剔除第3、5、6、7、8个IMF以及残差之后的IMF序列组合成最优IMF测量矩阵ΦE，该矩阵尺寸为85×33 400。相对于原始测量矩阵Φ0，该矩阵的尺寸大大减少，而精度也得到了提高。

表3 第2次剔除IMF的误差

表4 第3次剔除IMF的误差

表5 第4次剔除IMF的误差

为了体现迭代剔除法的优越性，将迭代剔除后的最优测量矩阵与随机抽取85个IMF组成的测量矩阵进行比较。图2给出了观测时间为100～10 000 s时，2种测量矩阵的脉冲星周期估计误差对比。与随机抽取相比，本文方法得出的测量矩阵可使脉冲星周期估计误差更小。

图2 迭代剔除与随机抽取的对比

4.3 有限等距性质

本节研究IMF测量矩阵是否满足有限等距性质(Restricted Isometry Property, RIP)。每个测量矢量与其他测量矢量之间的最大距离与最小距离如图3所示。从图中可以看出，单个测量矢量的最小距离在0.15～0.27之间，最大值在0.32～0.38之间，波形起伏小。这表明算法性能几乎与累积轮廓的相位和畸变度无关。因此，等距常量为0.85，测量矩阵满足1阶RIP。

图3 测量矢量之间的距离

接着，又研究了利用单个轮廓EMD分解，构造测量矩阵。从图4中可以看出，最大距离与最小距离相差约1～2个数量级，差距大。虽然测量矩阵也满足1阶RIP，但是，最小距离太小，易受噪声干扰。

图4 单个轮廓对应的测量矢量距离最值

4.4 基于FFT-CS的脉冲星周期估计

本文将基于EMD-CS的脉冲周期估计方法与FFT-CS方法进行对比，实验结果如表6所示。

表6 EMD-CS与FFT-CS的对比

从表6可看出，随着测量矩阵行数的增加，FFT-CS的估计误差减小。若测量矩阵行数为85，FFT-CS的估计误差为279 ps，远大于EMD-CS的57.9 ps。若要达到58 ps的精度，FFT-CS需要的矩阵尺寸为2 499×33 400，远大于EMD-CS的85×33 400。

从运行时间上看，当测量行数都为85的时候，FFT-CS方法运行时间是0.003 0 s，EMD-CS方法运行时间是0.002 1 s，若要达到EMD-CS方法中测量行数是85时候的精度，FFT-CS方法所需运行时间是0.015 5 s。

因此，EMD-CS只需较小的测量矩阵就可以获得较高的估计精度，且运行速度更快，明显优于FFT-CS方法。此外，当测量矩阵行数较大(上千)运行，FFT-CS的运行时间与矩阵行数几乎成正比，符合理论分析结论。

值得一提的是，器载计算机的计算能力有限，因此，这2种算法在轨运行时间必将大幅增长。所以，缩短运行时间是有一定实际意义的。

4.5 光子流量密度

脉冲星辐射光子流量密度和背景噪声流量密度是影响脉冲星周期估计精度的重要因素。本节研究两者对周期估计精度的影响，结果如图5所示。从图5中可以看出，脉冲星周期估计误差随背景噪声的增大而增大；而随辐射光子流量密度的增大而减小。

图5 不同流量下的周期估计误差

4.6 探测器面积与观测时间

X射线探测器面积和观测时间均为脉冲星周期估计精度的影响因素。本节给出了二者与脉冲周期估计误差的关系，仿真结果如图6所示。由图6(a)可看出，观测时间和探测器面积的增加均能减小脉冲星周期估计误差。因此，增加X射线探测器面积与观测时间，有助于提高精度。

图6 不同探测器面积和观测时间下的周期估计误差

为便于分析，本文在100 cm2、1 000 cm2、10 000 cm2这3种典型的探测器面积下研究估计误差，如图6(b)所示。可以看出，X射线探测器面积增大2个数量级，脉冲星周期估计精度提高约一个数量级。观测时间延长2个数量级，脉冲星周期估计精度提高约3个数量级。因此，相比于增加探测器面积，延长观测时间更加有利于提高估计精度。

值得一提的是，本文方法并未考虑延长观测时间所带来的轨道外推误差影响。这是因为本文方法适用于深空探测巡航段，此时的轨道模型近似于线性的[10]，轨道外推误差的影响可忽略不计。若在环绕段长时间观测，利用天文测角信息补偿[9]是一个较好的选择。