贾运泽,刘 劲,王奕迪,潘 超

(1.武汉科技大学信息科学与工程学院,武汉 430081;2.国防科技大学空天科学学院,长沙 410003;3.湖北经济大学信息与通信工程学院,武汉 430205)

0 引言

脉冲星是一颗高稳定度的中子星,而脉冲星导航[1-3]是一项不依赖于地面站的自主导航技术。同时,脉冲星能够为近地空间、太阳系乃至星际旅行的航天器提供定轨、守时等功能,满足航天任务在不同轨道下的持续高精度导航需求。脉冲星导航的基本原理为[4]:在脉冲星观测周期内,航天器上安装的X射线探测器会持续收集到脉冲星发出的X射线信号,并将这些信号在太阳系质心(solar system barycenter,SSB)处折叠成一个固定的脉冲星周期[5-6]。这个过程被称为历元折叠(epoch folding,EF),经过历元折叠后会形成一个累积脉冲轮廓,将之与脉冲标准轮廓进行比较可以得到一个相位差,进而可以将相位差转换为脉冲到达时间(time-of-arrival,TOA),将脉冲TOA 作为测量量即可求得航天器任意时刻的速度与位置信息。由此可以看出,脉冲星TOA 是脉冲星导航中的重要参数,但是因为航天器的运动加上脉冲传播过程中的噪声干扰[7-8],会导致接收到的脉冲星信号周期发生变化,使积聚在固有周期内的脉冲星累积轮廓发生畸变,从而引起脉冲星TOA 漂移[9]。因此,脉冲星周期的观测对高精度导航有着重要的意义,当前如何快速且高精度地获取脉冲星周期已经成为了研究的热点。

脉冲星周期估计的方法可分为两类:一类是利用脉冲星TOA 的漂移来估计周期误差,另一类是利用脉冲星轮廓畸变来反演脉冲星周期,后者为当前的主流。其基本原理为:尝试按不同周期折叠脉冲星信号,得到一系列脉冲畸变轮廓,再找出畸变度最小的脉冲累积轮廓,其对应的周期就是固有周期。目前最典型的就是时域χ2统计法[10]和傅立叶频谱法[11]。除此之外,还有文献[12]提出的一种基于频率细分和连续伦姆周期图(CLP)法,文献[13]提出的基于TOA 信息的脉冲星周期估计法等等。以上方法都有一个共同的问题,那就是多次EF 带来的高计算量,而对于深空计算设备来说,高计算量带来的高负载是必须要避免的问题之一。

压缩感知(compressive sensing,CS)[14-15]是一种新兴的信号处理方法,它可以通过极少的观测值重构脉冲轮廓,从而解决脉冲星数据量过大这一问题。CS对于稀疏信号有较好的恢复能力,而脉冲信号就是一个稀疏信号,因此,将CS应用到脉冲信号处理上已经成为一个热点话题。在CS过程中,测量矩阵的构造和字典的设计是关键,如孙海峰等[16]将哈达玛矩阵作为测量矩阵重构脉冲星信号,吴春艳等[17]提出了基于粗估计与精估计两级字典的时延估计法。刘劲等[18]将CS应用在脉冲星周期估计中,提出一种基于FFT-CS的脉冲星周期快速估计方法,并将其与蝶形算法结合,减小了计算量并提高了精度。除此之外,CS 在其他方面也得到了应用,如TOA 估计[19]、航天器定位[20]等。

为了进一步减少计算量并提高计算精度,本文提出了一种基于DCT-CS的脉冲星周期快速估计方法。通过构造不同程度的畸变度字典来直接估计脉冲星周期。脉冲星信号的能量主要集中在低频部分,将信号通过DCT变换为测量矩阵,然后经过超分辨率稀疏恢复直接估计出X射线脉冲星周期。

1 轮廓畸变对DCT的影响

本文提出的基于DCT-CS 的脉冲星周期超快速估计方法需要对脉冲星畸变轮廓进行DCT 计算来得到脉冲星周期。

在脉冲星导航的实际应用中,由于航天器运动会造成多普勒效应,导致观测信号周期改变。此时如果仍然使用提前预估的脉冲星周期而不是脉冲星的固有周期进行历元折叠,就会导致累积脉冲轮廓产生轮廓畸变[21]。

离散余弦变换(discrete cosine transform,DCT)是从离散傅里叶变换(discrete Fourier transform,DFT)推导而来,是输入函数为偶函数情况下的一种特殊的DFT,相当于一个长度为其2倍的实偶序列的DFT。相比于DFT,DCT 主要有以下三个优点:1)DCT 中不存在虚部,计算主要为实数变换,具有更好的计算效率。2)在不引入间断的同时对信号施加了周期性。在DFT 中,当时间信号被截断并假定其周期性时,会在时域中引入不连续性,而DCT 即使假设信号存在对称性,也不会在信号中引入不连续性和相关操作。3)DCT 比DFT 具有更好的能量聚集度,能够将能量聚集在低频部分后进行处理。一维DCT 的表达式如式(1)所示。

式中,f(x)为输入信号;F(k)为输出信号;k为输出信号点位数;Y为信号长度。

图1为标准脉冲轮廓和畸变轮廓进行DCT 的结果对比。图1(a)为标准轮廓进行DCT 的结果,图1(b)为使用最大轮廓畸变度M=21的累积轮廓进行DCT 的结果,其中轮廓畸变度为畸变宽度与一个脉冲周期内的周期间隔数N之比,表示轮廓畸变的程度。DCT 系数为脉冲轮廓信号经DCT 后得到的数值,其基本特性为直流和低频系数数值大,高频系数数值小,幅值表示信号强度。由图1可以看出,随着轮廓的累积,信号强度会逐渐增大,同时信号能量会向周边分散。因此,相比于直接对标准轮廓进行DCT 计算,对累积轮廓进行DCT 能更好地利用脉冲信号的低频部分。

图1 标准轮廓与畸变轮廓DCT结果Fig.1 DCT results of standard contour and distorted contour

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

提出了一种基于DCT-CS 的脉冲星周期快速估计算法。在传统脉冲星周期估计方法使用的历元折叠方法中,脉冲信号将会以不同的周期进行折叠并生成一系列的畸变轮廓,通过求解最小失真度来求解脉冲星周期。而本算法不需要进行尝试性历元折叠,该算法的核心思想为通过检测脉冲星轮廓的畸变度来直接估计脉冲星周期误差。

本算法包括以下四个部分:1)构建畸变轮廓字典。2)设计低频离散测量矩阵。3)获取频率偏移累积轮廓。4)使用最大元素超分辨率稀疏恢复。算法示意图如图2所示。其中,测量矩阵和字典的构建是提前在地面完成的;脉冲星轮廓的积累和超分辨率稀疏恢复在航天器上进行。脉冲星轮廓的积累与脉冲星信号的采集是同步的。在收集完X射线脉冲星光子后,才会开始计算匹配矩阵和使用最大元素超分辨率稀疏恢复。该算法的目标是在脉冲星信号采集后快速获得脉冲星周期。因此,减少算法第四部分的计算负荷是主要任务。

图2 基于DCT-CS的脉冲星周期超快速估计算法Fig.2 Ultra fast estimation algorithm for pulsar period based on DCT-CS

在本算法中也用到了历元折叠,为了与传统脉冲星周期估计方法中的历元折叠区分,将本算法中的历元折叠命名为基本历元折叠,将传统算法中的历元折叠命名为尝试性历元折叠。基本历元折叠相比于尝试性历元折叠有以下区别:1)计算量不同,基本历元折叠在整个脉冲星周期估计算法中只进行了一次,而传统脉冲星周期估计方法中需要进行多次尝试性历元折叠。2)计算时间不同,基本历元折叠可以在脉冲星信号的采集过程中同步进行,而尝试性历元折叠需要在脉冲星信号的采集之后进行,这大大节省了脉冲星周期估计的计算时间。

下面将对算法流程的四个部分进行逐个介绍。

(1)构造畸变轮廓字典

不同频率偏移引起的累积脉冲星轮廓的畸变程度不同,所以构造畸变轮廓字典是必要的。

设h(φ)为归一化的标准脉冲轮廓,其满足以下条件

式中,φ为脉冲相位;h(φ)为一个周期为1的函数,则有以下公式

畸变脉冲轮廓可视为标准脉冲轮廓与门函数的循环互相关,因此,构造了M个不同畸变度的脉冲轮廓φm,表达式为

式中,⊗循环相关;m是畸变脉冲轮廓序号;^h为反向标准脉冲轮廓,表达式如式(5)所示

Gm(N×1)为传播矢量,表达式如式(6)所示

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

将M个φm合成畸变轮廓合成一个畸变字典ψ(N×M),其可以表示为

式中,m={1,2,3…M-1};M为脉冲星最大轮廓畸变度,文中为21。

(2)设计低频离散测量矩阵

因为脉冲星累积轮廓的相位是未知的,所以脉冲星周期估计方法必须不受相位影响。因为DCT系数与相位无关,所以本文选用DCT 来求得测量矩阵,由于脉冲星轮廓的能量集中在低频部分,所以取用低频部分提取测量矩阵。

文中的低频测量矩阵由字典中的一个元素φm经DCT 得出,则低频离散测量矩阵Φ的表达式为

假设对应字典因素φm的频率偏移为Δfm,可以表示为

式中,N为脉冲轮廓间隔数。

可以构造感知矩阵Θ,即测量矩阵与轮廓字典的乘积,可表示为

(3)获取频率偏移累积轮廓

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

式中,Δf是频率偏移;f0是固有频率。

如果不考虑相位,不能仅通过累积轮廓来区分Δf和-Δf,为了解决该问题,在此处引入脉冲星周期偏移量Toff≫ΔT,-Toff+ΔT和-Toff-ΔT的符号是相同的,但是其振幅不同,由此可以区分出频率的正负符号。在基本历元折叠中,采用T0-Toff+ΔT为周期进行轮廓累积并进行如下操作:

收集一个观测周期内的所有光子信号,并将其按照[0,T0-Toff+ΔT]的脉冲星周期进行折叠,这就意味着其在一个观测周期内的计算模量为

然后,将周期持续时间分为等长的几块,计算每块中的光子量为

最后,对得到的光子进行归一化处理,得到累积脉冲轮廓。

(4)使用最大元素超分辨率稀疏恢复

设观测矢量y的表达式为

匹配矢量的第m项zm,(m={0,1,2,…,M-1})可表示为观测矢量和感知矢量的乘积

式中,am为感知矩阵Θ的第m列。

设匹配矩阵为S,其第m列Sm可表示为

式中,L为感知矩阵测量次数,Θ(Μ,1∶L)代表感知矩阵第M行的1到L列。

然后得到S最大值的行下标和列下标

最后,脉冲星周期可估计为

3 计算复杂度分析

本章中将计算整个算法的计算复杂度,算法流程如图2所示。其中,脉冲星累积过程可在脉冲星信号的观测时间内进行,因此此处只计算提取低频部分、匹配字典和超分辨率匹配估计的计算量,即乘-加操作次数(multiply-accumulate operations,MAC)。

1)提取低频部分:在此部分对字典中的单个元素φm进行了DCT来提取信号低频部分,所以计算量集中在DCT 部分,对于一个长度为Y的信号,一维DCT的计算量为Y2。设脉冲星累积轮廓间隔数为N,因此在此部分中计算量为N2MAC。

2)匹配字典:在此部分中进行了一次DCT 计算,如公式(10)所示,所以计算量依然为N2MAC。

3)超分辨率匹配估计:该部分的计算量主要由匹配矢量和匹配矩阵决定,其中匹配矢量如式(16)所示,共需要(M×L)MAC。匹配矩阵包含式(17)和式(18)两方面,设匹配矩阵大小为L,故其中式(15)包含L次计算,式(18)包含L次计算,要计算出max(Sm)需要计算一次式(17)和M次式(18),共需(2L+M×L)MAC。

则算法整体需2N2+2L+2M×L次计算。在预先设定N=33 000,L=5 000,M=21的情况下,共需约2.2×109MAC。

为了对比体现本算法的优势,下面计算同样情况下使用FFT 算法进行脉冲星周期估计所需的计算量。同样计算提取信号低频部分,匹配字典和超分辨率匹配估计三部分的计算量。在提取信号低频部分中,计算量可近似看作FFT 部分的计算量,其计算量为N×log2NMAC。而在匹配字典部分,需要先进行FFT 再进行IFFT,所以需要进行2N×log2NMAC。最后的超分辨率匹配估计,其中匹配矢量共计算M×L次,由于FFT 需要进行二维计算,所以式(18)的计算量变为2L次,则计算出max(Sm)需要进行一次式(17),M次式(18)运算,共需要计算(3L+M×L)MAC。则采用FFT 进行脉冲星周期估计共需:N×log2N+2N×log2N+M×L+3L+M×L=3N×log2N+3L+3M×L,约为1.8×106MAC。

从上述计算量分析可知,基于DCT-CS的脉冲星周期估计方法在前期的低频测量矩阵获取与历元折叠过程中的计算量较大,其原因在于DCT 的计算量与脉冲轮廓间隔数的平方,即N2成正比。当已知航天器的速度和位置信息后,可减少搜索点的数量,届时可大幅减少计算量。而在其后的脉冲星周期估计算法部分里,基于DCT-CS的脉冲星周期估计方法的计算量相较于FFT-CS较少。其原因在于DCT 没有虚数部分且算法流程更为简单,所以计算量较少。从结果来看,基于DCT-CS的脉冲星周期估计方法从计算量角度与基于FFT-CS的周期估计方法各有优劣,需要通过实验验证分析。

4 仿真结果与讨论

本次仿真实验中使用了欧洲脉冲星网络(European pulsar network,EPN)的数据,此数据库收集了超过一千颗脉冲星的标准脉冲轮廓数据,而且,数据都是采用ASCII字符存储,使用者可以直接免费下载其纯文本文件,非常方便。本实验选取的脉冲星为蟹状星云脉冲星PSR B0531+21,其标准脉冲轮廓图如图3所示,之所以选这颗脉冲星,是因为其具有较大的光子流量密度,便于观测和实验。

图3 脉冲星PSR B0531+21标准轮廓Fig.3 Standard profile of pulsar PSR B0531+21

设时间分辨率为1μs。脉冲轮廓间隔数N等于脉冲星周期除以时间分辨率的商,其值为33 000。脉冲星辐射光子流量密度为1.54 ph/(cm2·s),背景噪声光子流量密度是0.005 ph/(cm2·s)。X射线探测器的有效探测面积为1 m2,观测时长为1 000 s,Toff为2.7×10-10s。测量矩阵大小为L×N,其中L为5 000,N为33 000。字典大小为M×L,其中M为21。实验设备为一台处理器为AMD7840@3.80 GHz,RAM 16 G的电脑。

4.1 观测时长和探测器面积的影响

本节研究观测时长和探测器面积对脉冲星周期估计误差的影响,同时通过对比分析不同观测时长和探测器面积对结果的影响程度。

图4 给出了不同观测时长和探测器面积下的仿真结果,图中横纵坐标轴为对数坐标值。由图4可以看出脉冲星周期估计误差随观测时长和探测器面积增大而减少,且观测时长对脉冲星周期估计误差的影响强于探测器面积。由此可以得出,无论是面积的扩大还是观测时间的延长,都有利于脉冲星周期精度的提高。而且,观测时间的延长比观测面积的增大更能提高周期估计的精度。

图4 时间面积对周期估计误差的影响Fig.4 The influence of time and area on period estimation error

其原因为根据脉冲星周期估计误差的克拉美罗下界(Cramer-Rao lower bound,CRLB)推导得出,脉冲星周期估计误差与探测器面积和观测时长的立方成反比[18],其性质可以表示为

式中,Te为脉冲星周期估计误差;Tobs和A分别为观测时长和探测器面积。由式(23)和式(24)可以看出提高观测时长和增大探测器面积均可以提高脉冲星周期估计的精度,且提高观测时长的影响比增大探测面积更大,而实验结果与该性质相符。

4.2 测量次数的影响

本节研究了测量次数对脉冲星周期估计误差的影响,由于脉冲星轮廓具有较低的信噪比,且脉冲星轮廓的能量集中在低频部分,因此本算法通过DCT 提取矩阵中选择低频部分,而不是传统的随机提取。

图5表明了测量次数与脉冲星周期估计误差的关系,图中纵坐标轴为对数值。可以看出,随着测量次数的增加,脉冲星周期估计误差逐渐减小,直至某一点到达最小值,对于不同的探测器面积,达到最小值的测量次数均为约5 000,大于5 000时误差值趋近于定值。其原因为脉冲星轮廓的能量集中在低频部分,因此小尺寸的测量矩阵就包含了绝大部分能量,后续再增加测量次数能得到的能量很小且会引入噪声。因此本次实验中采用L=5 000作为测量矩阵的测量次数。

图5 测量次数对周期估计误差的影响Fig.5 The influence of observation number on period estimation error

4.3 光通量与背景噪声的影响

光通量和背景噪声是脉冲星周期估计中的重要参数,本节研究了光通量和背景噪声对脉冲星周期估计误差的影响。

图6给出了3 000次蒙特卡罗实验中不同光通量和背景噪声下脉冲星周期估计的误差结果,图中横、纵坐标轴为对数坐标值。由图6可以看出,脉冲星周期估计误差随着脉冲星光子通量的增加而减小,而随着背景噪声的增加而增大,脉冲星周期估计误差随着脉冲星光子通量的增加而减小,而随着噪声级的增加而增大。因此,无论是增大脉冲星光子通量还是减小背景噪声,都有利于提高脉冲星周期估计精度。

图6 光通量和噪声对周期估计误差的影响Fig.6 Impacts of flux and noise on period estimation error

4.4 DCT与FFT的对比

本算法中使用DCT 替换了传统的FFT 来进行脉冲星周期估计,现对比这两种算法的优劣。

表1给出了不同观测时长和探测器面积的情况下,采用FFT 和DCT 进行脉冲星周期估计的误差结果和进行历元折叠需要的计算时间。由表1可以看出,DCT 的计算时间少于FFT 且与参数无关,这是因为本算法的计算时间与光子量无关。而在周期估计误差方面,在观测时长较短和探测器面积较小的情况下,采用FFT 的误差结果会远大于采用DCT 得到的结果,但随着观测时长的增加和探测器面积的增大,两者的估计误差结果会逐渐接近,最终采用FFT 得到的结果会优于采用DCT 方法得到的结果。其原因在于FFT 较DCT 多出来一个虚数部分,计算精度损失低于DCT,具有较高的鲁棒性,在高分辨率的情况下,达到的结果会好于DCT。

表1 两种周期估计方法的对比Tab.1 Comparison of two cycle estimation methods

从理论角度出发,DCT 对比FFT 有以下几个优势:1)DCT 相比于FFT 有更好的能量聚集度,能更好地处理脉冲信号。2)DCT 没有虚部,相比于FFT 拥有更小的计算复杂度。从实际角度出发,考虑到航天器高速飞行过程中产生的多普勒效应的影响,需尽快完成脉冲星周期估计的计算过程,才能保证算法的实时性。根据表1的实验结果表明,基于DCT 的脉冲星周期估计方法需要的计算时间更短。而基于FFT 的周期估计方法只有在高分辨率的情况下的误差估计结果才会优于基于DCT 的脉冲星周期估计方法,限制条件过高。因此本文采用DCT 进行脉冲星周期估计。

5 结论

本文提出了一种基于DCT-CS 的脉冲星周期超快速估计算法,该算法不需要进行尝试性历元折叠,而是通过低频测量矩阵、轮廓畸变字典和超分辨率稀疏恢复算法来直接估计脉冲星周期。

本算法有三个优点:1)计算负荷低。DCT 对低频信号有着优秀处理能力,该算法舍弃了常用的FFT方法,而是通过DCT 来进行低频矩阵的提取和构建测量矩阵。理论分析表明基于DCT的脉冲星周期估计方法计算负荷要低于基于FFT的脉冲星周期估计方法。2)计算精度高。文中通过实验评估了基于DCT的脉冲星周期估计方法的计算精度,实验结果表明其计算精度达到了3.82×10-12s,优于传统脉冲星周期估计方法。3)计算时间短。由于基于DCT的脉冲星周期估计方法的计算量较低,所以相应的计算时间也较短。实验结果表明基于DCT的脉冲星周期估计方法的计算时间仅为9.31 ms,要少于基于FFT的脉冲星周期估计方法的14.9 ms。

综上所述,基于DCT-CS的脉冲星周期快速估计算法可在降低系统计算负担的同时提高周期估计精度。