劉勁，韓雪俠，寧曉琳，陳曉，康志偉

1. 武漢科技大學(xué) 信息科學(xué)與工程學(xué)院，武漢 430081

2. 北京航空航天大學(xué) 儀器科學(xué)與光電工程學(xué)院，北京 100191

3. 上海衛(wèi)星工程研究所，上海 200240

4. 湖南大學(xué) 信息科學(xué)與工程學(xué)院，長(zhǎng)沙 410082

脈沖到達(dá)時(shí)間(Time-of-Arrival, TOA)是X射線脈沖星導(dǎo)航的基本觀測(cè)量[1-2]。航天器上的X射線探測(cè)器收集X射線脈沖星輻射的光子[3]，經(jīng)歷元折疊形成累積脈沖輪廓[4-5]，再結(jié)合標(biāo)準(zhǔn)脈沖輪廓、脈沖星計(jì)時(shí)模型等信息估計(jì)脈沖到達(dá)時(shí)間。但航天器運(yùn)動(dòng)和脈沖星周期躍變[6]等導(dǎo)致接收到的脈沖信號(hào)周期發(fā)生變化，使得按固有周期累積的脈沖輪廓發(fā)生畸變，進(jìn)而導(dǎo)致脈沖TOA發(fā)生漂移。脈沖周期估計(jì)對(duì)脈沖TOA高精度估計(jì)具有重要意義，目前已成為一個(gè)研究熱點(diǎn)。

脈沖星周期估計(jì)方法可分為2類：一類利用脈沖星TOA漂移估計(jì)周期誤差；另一類根據(jù)脈沖輪廓畸變反演周期誤差，是目前的主流方法。

基于脈沖星TOA漂移的周期誤差估計(jì)的基本思路如下：建立周期誤差與脈沖星TOA漂移之間的關(guān)系模型。在此基礎(chǔ)上，利用最小二乘法估計(jì)TOA與周期誤差[7-8]；利用天文信息補(bǔ)償[9]或優(yōu)化卡爾曼濾波器，如：EKF-CMBSEE(Extended Kalman Filter-Correlated Measurement Bias and State Estimation Error)[10]、跟蹤濾波器[11]，以達(dá)到抑制周期誤差影響的目的。

基于輪廓畸變的脈沖星周期估計(jì)方法的基本思路如下：嘗試按不同周期折疊脈沖星信號(hào)，得到一系列畸變脈沖累積輪廓，然后找出畸變度最小的脈沖累積輪廓，其對(duì)應(yīng)的周期就是固有周期。如：時(shí)域χ2法、快速折疊法[12]、傅立葉頻譜法[13]、循環(huán)平穩(wěn)信號(hào)相干統(tǒng)計(jì)量法[14]、最大相關(guān)方差搜索法[15]、基于Lomb的χ2算法[16]、脈沖輪廓熵算法[17]和輪廓特征法[18]等。上述方法無一例外都需要嘗試按不同周期折疊脈沖星信號(hào)。但是，脈沖星輻射信號(hào)數(shù)據(jù)量龐大，導(dǎo)致計(jì)算量大等問題，不適合星載計(jì)算機(jī)在軌運(yùn)行。

壓縮感知(Compressed Sensing，CS)[19-20]是一種新興的信號(hào)處理方法，對(duì)稀疏信號(hào)的處理能力很強(qiáng)，恰巧脈沖星信號(hào)是典型的稀疏信號(hào)。近年來，基于CS的脈沖星信號(hào)處理成為一個(gè)研究熱點(diǎn)[21-23]?；诳焖俑道锶~變換-壓縮感知(FFT-CS)的脈沖星周期快速估計(jì)方法[24]利用CS壓縮脈沖星輻射信號(hào)，再直接根據(jù)脈沖星輪廓畸變度估計(jì)脈沖固有周期。該方法避免了多次脈沖星信號(hào)折疊，大幅減少了計(jì)算量，使得脈沖星周期在軌估計(jì)成為可能。

但是，傅立葉變換將信號(hào)能量分散在數(shù)量龐大的傅立葉級(jí)數(shù)中。FFT-CS采用低頻傅立葉變換，所需的傅立葉級(jí)數(shù)仍然上千，這導(dǎo)致測(cè)量矩陣尺寸極大。為提高計(jì)算速度，必須大幅減小尺寸。經(jīng)驗(yàn)?zāi)B(tài)分解(Empirical Mode Decomposition, EMD)[25-26]根據(jù)原始信號(hào)的固有特征，自適應(yīng)地將其分解成有限個(gè)固有模態(tài)函數(shù)(Intrinsic Mode Function，IMF)。IMF包含了原始信號(hào)不同時(shí)間尺度的局部特征信號(hào)。而脈沖輪廓畸變度這一微弱局部特征可體現(xiàn)在某些IMF中。與傅立葉變換相比，IMF數(shù)量大幅減少，僅為101～102量級(jí)。因此，本文將IMF構(gòu)成的測(cè)量矩陣替換FFT-CS中的低頻傅立葉變換矩陣，提出了一種基于EMD-CS的脈沖星周期超快速估計(jì)方法，旨在保持估計(jì)精度的同時(shí)，減小計(jì)算量，加快運(yùn)行速度。

1 最優(yōu)IMF測(cè)量矩陣

為減小測(cè)量矩陣尺寸，本文利用EMD分解得到的IMF構(gòu)造測(cè)量矩陣來替代低頻傅立葉矩陣，并在此基礎(chǔ)上提出了一種EMD-CS方法。EMD-CS方法的基本流程與FFT-CS方法類似，二者唯一不同在于測(cè)量矩陣。本節(jié)重點(diǎn)研究測(cè)量矩陣。

在CS中，測(cè)量矢量的獲取是通過測(cè)量矩陣來實(shí)現(xiàn)的。通常，測(cè)量矩陣尺寸越小，計(jì)算量越少。因此尋找一種測(cè)量矩陣，用更少的測(cè)量次數(shù)達(dá)到預(yù)期的結(jié)果變得尤為重要。

將這一系列IMF序列組合而成一個(gè)IMF原始測(cè)量矩陣Φ0(MIMF×N)，N為一個(gè)脈沖周期內(nèi)的脈沖間隔數(shù)；MIMF為不同畸變輪廓分解的IMF序列數(shù)總和,可表示為

(1)

Φ0可以表示為

(2)

(3)

(4)

(5)

MEIMF的值為

(6)

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

迭代剔除法在地面進(jìn)行，篩選得到的最優(yōu)IMF測(cè)量矩陣存入器載計(jì)算機(jī)，無需在軌篩選。

2 EMD-CS

本節(jié)將介紹EMD-CS方法估計(jì)脈沖星周期的整個(gè)算法流程，如圖1所示分成4個(gè)模塊：① 最 優(yōu)IMF測(cè)量矩陣的設(shè)計(jì)；② 畸變字典的構(gòu)造；③ 脈沖頻率初始偏置；④ 最大元素的超分辨率稀疏恢復(fù)。其中，最優(yōu)IMF測(cè)量矩陣已經(jīng)在第1節(jié)中描述。下面將介紹其他3個(gè)模塊：

圖1 所提算法流程圖

1) 畸變字典的構(gòu)造

文獻(xiàn)[24]已證明，畸變脈沖輪廓可視為標(biāo)準(zhǔn)脈沖輪廓與門函數(shù)的循環(huán)互相關(guān)。因此，構(gòu)造了MD個(gè)不同畸變度的脈沖輪廓hm，表達(dá)式為

hm=Gm?h

(7)

式中：Gm(N×1維)為擴(kuò)散矢量；m為畸變脈沖輪廓序號(hào)；h(N×1維)為標(biāo)準(zhǔn)脈沖輪廓；?表示循環(huán)互相關(guān)。Gm滿足：

(8)

從式(7)和式(8)可看出，Gm相當(dāng)于MD個(gè)寬度不同的矩形窗。脈沖星輪廓畸變度定義為畸變寬度與一個(gè)脈沖周期內(nèi)的周期間隔數(shù)N之比，即m/N。隨著m增大，矩形窗變寬，脈沖累積輪廓的畸變度變大。

脈沖星導(dǎo)航常用于深空探測(cè)巡航段，此時(shí)的軌道模型近似于線性的[10]，所以可把脈沖星周期誤差視為非時(shí)變的。式(8)所建立的模型適用于這種情況。

將第m個(gè)畸變脈沖輪廓循環(huán)移位構(gòu)成第m個(gè)畸變脈沖輪廓子字典φm(N×P維)，表達(dá)式為

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

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

(9)

式中:P為最大相位偏移量；φmi(t)為子字典φm中的第i個(gè)原子；hm(t)為第m個(gè)畸變脈沖輪廓?？蓪D個(gè)子字典φm合成一個(gè)三維矩陣字典Ψ(N×P×MD維)。

2) 脈沖頻率初始偏置

累積脈沖輪廓是根據(jù)脈沖固有周期T0將脈沖星輻射光子進(jìn)行折疊而形成的，這一過程稱為歷元折疊。但是，多普勒效應(yīng)以及脈沖星周期躍變會(huì)使航天器接收到的脈沖周期存在誤差ΔT，與頻率偏移之間的關(guān)系可表示為

(10)

式中：Δf為頻率偏移；f0為固有頻率。

文獻(xiàn)[24]已經(jīng)證明，不能直接根據(jù)脈沖星輪廓畸變辨別Δf的正負(fù)號(hào)。為了解決此問題，引入脈沖星周期偏移Toffset，且滿足Toffset?ΔT。這樣，-Toffset+ΔT和-Toffset-ΔT的符號(hào)都是負(fù)號(hào)，但幅度不同。所以，脈沖輪廓按照周期T0-Toffset+ΔT進(jìn)行歷元折疊形成累積脈沖輪廓。

3) 基于最大值的超分辨率稀疏恢復(fù)

在壓縮感知中，信號(hào)稀疏恢復(fù)是通過感知矩陣與測(cè)量矢量相匹配，根據(jù)最大匹配值實(shí)現(xiàn)的。感知矩陣T是最優(yōu)IMF測(cè)量矩陣與字典的乘積：

T=ΦEΨ

(11)

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

測(cè)量矢量y(MEIMF×1維)的獲得方式為

(12)

感知矩陣與測(cè)量矢量的乘積為匹配矩陣S，S可表示為

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

i=1,2，…，P;

j=1,2,…,MD

(13)

(14)

為了使估計(jì)更加準(zhǔn)確，運(yùn)用基于最大值的超分辨率稀疏恢復(fù)：

(15)

(16)

(17)

3 復(fù)雜度分析

本節(jié)分析EMD-CS方法的計(jì)算復(fù)雜度。在該方法中，脈沖星輪廓累積、測(cè)量矢量獲取和匹配在航天器上運(yùn)行。其中，脈沖星累積過程可在脈沖星信號(hào)的觀測(cè)時(shí)間內(nèi)進(jìn)行。而測(cè)量矢量獲取和匹配是在觀測(cè)結(jié)束后才進(jìn)行的?？紤]到航天器高速飛行的影響，需盡快完成這2個(gè)過程，才能保證算法的實(shí)時(shí)性。因此，本文分析二者計(jì)算復(fù)雜度即可。

測(cè)量矢量的獲取由式(12)確定。最優(yōu)IMF測(cè)量矩陣尺寸為MEIMF×N，脈沖累積輪廓為長(zhǎng)度為N。這需要MEIMF×(2N-1)MAC，MAC(Multiply/Accumulate Computation)表示加法乘法的計(jì)算次數(shù)。

匹配由式(13)確定。測(cè)量矢量與一個(gè)原子匹配，所需計(jì)算量為2MEIMF-1 MAC。字典中原子個(gè)數(shù)為MD×P。測(cè)量矢量與字典匹配所需計(jì)算量是(2MEIMF-1)×MD×PMAC。

總計(jì)算量為MEIMF×(2N-1)+(2MEIMF-1)×MD×P≈2MEIMF×(N+MD×P) MAC。

以上可以看出，計(jì)算量與MEIMF呈正比關(guān)系。若能大幅減小MEIMF，計(jì)算量也將大幅減少。此外，與FFT-CS相比，EMD-CS僅有實(shí)數(shù)部分，計(jì)算量將更小。

4 仿真實(shí)驗(yàn)

將基于EMD-CS的脈沖星周期估計(jì)方法與基于FFT-CS的周期估計(jì)進(jìn)行仿真比較。首先，篩選出最優(yōu)的測(cè)量矩陣，并分析其有限等距性質(zhì)；其次，將其與基于FFT-CS的脈沖星周期估計(jì)方法作比較，體現(xiàn)本文方法的優(yōu)越性；再次，分析脈沖星和背景噪聲的光子流量密度、X射線探測(cè)器面積以及觀測(cè)時(shí)間對(duì)脈沖星周期估計(jì)精度的影響。

4.1 仿真條件

脈沖星數(shù)據(jù)來自于歐洲脈沖星網(wǎng)(EPN)，仿真實(shí)驗(yàn)在處理器為英特爾Core i5-7500@3.40 GHz四核、內(nèi)存為8 GB的計(jì)算機(jī)上運(yùn)行。最大畸變度MD=21，即構(gòu)造21個(gè)脈沖星畸變輪廓，最大相位偏移量P=21。仿真條件如表1所示。

表1 PSR B0531+21仿真條件

4.2 最優(yōu)測(cè)量矩陣

表2給出了剔除某些IMF的脈沖星周期估計(jì)誤差，表的首行和首列表示剔除的IMF序號(hào)。從表2可以看出，若是剔除第1個(gè)IMF，脈沖星周期估計(jì)誤差大，所以必須保留第1個(gè)IMF；而同時(shí)剔除殘差和第5個(gè)IMF時(shí)，誤差最小(70.177 8 ps)。

表2 第1次剔除IMF的誤差

本文保留第1個(gè)IMF，剔除殘差和第5個(gè)IMF后，在剩下的IMF中繼續(xù)實(shí)施剔除操作，仿真結(jié)果如表3所示。從表3中可以看出，剔除第6和第7個(gè)IMF，誤差最小(62.224 8 ps)。表4給出了剔除第5、6、7個(gè)IMF以及殘差之后的周期估計(jì)。實(shí)驗(yàn)結(jié)果表明，再剔除第3和第8個(gè)IMF，估計(jì)誤差更小(57.063 6 ps)。表5給出了在此基礎(chǔ)上，4次剔除之后的結(jié)果。此次剔除，使得誤差增大。因此，本文將IMF剔除第3、5、6、7、8個(gè)IMF以及殘差之后的IMF序列組合成最優(yōu)IMF測(cè)量矩陣ΦE，該矩陣尺寸為85×33 400。相對(duì)于原始測(cè)量矩陣Φ0，該矩陣的尺寸大大減少，而精度也得到了提高。

表3 第2次剔除IMF的誤差

表4 第3次剔除IMF的誤差

表5 第4次剔除IMF的誤差

為了體現(xiàn)迭代剔除法的優(yōu)越性，將迭代剔除后的最優(yōu)測(cè)量矩陣與隨機(jī)抽取85個(gè)IMF組成的測(cè)量矩陣進(jìn)行比較。圖2給出了觀測(cè)時(shí)間為100～10 000 s時(shí)，2種測(cè)量矩陣的脈沖星周期估計(jì)誤差對(duì)比。與隨機(jī)抽取相比，本文方法得出的測(cè)量矩陣可使脈沖星周期估計(jì)誤差更小。

圖2 迭代剔除與隨機(jī)抽取的對(duì)比

4.3 有限等距性質(zhì)

本節(jié)研究IMF測(cè)量矩陣是否滿足有限等距性質(zhì)(Restricted Isometry Property, RIP)。每個(gè)測(cè)量矢量與其他測(cè)量矢量之間的最大距離與最小距離如圖3所示。從圖中可以看出，單個(gè)測(cè)量矢量的最小距離在0.15～0.27之間，最大值在0.32～0.38之間，波形起伏小。這表明算法性能幾乎與累積輪廓的相位和畸變度無關(guān)。因此，等距常量為0.85，測(cè)量矩陣滿足1階RIP。

圖3 測(cè)量矢量之間的距離

接著，又研究了利用單個(gè)輪廓EMD分解，構(gòu)造測(cè)量矩陣。從圖4中可以看出，最大距離與最小距離相差約1～2個(gè)數(shù)量級(jí)，差距大。雖然測(cè)量矩陣也滿足1階RIP，但是，最小距離太小，易受噪聲干擾。

圖4 單個(gè)輪廓對(duì)應(yīng)的測(cè)量矢量距離最值

4.4 基于FFT-CS的脈沖星周期估計(jì)

本文將基于EMD-CS的脈沖周期估計(jì)方法與FFT-CS方法進(jìn)行對(duì)比，實(shí)驗(yàn)結(jié)果如表6所示。

表6 EMD-CS與FFT-CS的對(duì)比

從表6可看出，隨著測(cè)量矩陣行數(shù)的增加，F(xiàn)FT-CS的估計(jì)誤差減小。若測(cè)量矩陣行數(shù)為85，F(xiàn)FT-CS的估計(jì)誤差為279 ps，遠(yuǎn)大于EMD-CS的57.9 ps。若要達(dá)到58 ps的精度，F(xiàn)FT-CS需要的矩陣尺寸為2 499×33 400，遠(yuǎn)大于EMD-CS的85×33 400。

從運(yùn)行時(shí)間上看，當(dāng)測(cè)量行數(shù)都為85的時(shí)候，F(xiàn)FT-CS方法運(yùn)行時(shí)間是0.003 0 s，EMD-CS方法運(yùn)行時(shí)間是0.002 1 s，若要達(dá)到EMD-CS方法中測(cè)量行數(shù)是85時(shí)候的精度，F(xiàn)FT-CS方法所需運(yùn)行時(shí)間是0.015 5 s。

因此，EMD-CS只需較小的測(cè)量矩陣就可以獲得較高的估計(jì)精度，且運(yùn)行速度更快，明顯優(yōu)于FFT-CS方法。此外，當(dāng)測(cè)量矩陣行數(shù)較大(上千)運(yùn)行，F(xiàn)FT-CS的運(yùn)行時(shí)間與矩陣行數(shù)幾乎成正比，符合理論分析結(jié)論。

值得一提的是，器載計(jì)算機(jī)的計(jì)算能力有限，因此，這2種算法在軌運(yùn)行時(shí)間必將大幅增長(zhǎng)。所以，縮短運(yùn)行時(shí)間是有一定實(shí)際意義的。

4.5 光子流量密度

脈沖星輻射光子流量密度和背景噪聲流量密度是影響脈沖星周期估計(jì)精度的重要因素。本節(jié)研究?jī)烧邔?duì)周期估計(jì)精度的影響，結(jié)果如圖5所示。從圖5中可以看出，脈沖星周期估計(jì)誤差隨背景噪聲的增大而增大；而隨輻射光子流量密度的增大而減小。

圖5 不同流量下的周期估計(jì)誤差

4.6 探測(cè)器面積與觀測(cè)時(shí)間

X射線探測(cè)器面積和觀測(cè)時(shí)間均為脈沖星周期估計(jì)精度的影響因素。本節(jié)給出了二者與脈沖周期估計(jì)誤差的關(guān)系，仿真結(jié)果如圖6所示。由圖6(a)可看出，觀測(cè)時(shí)間和探測(cè)器面積的增加均能減小脈沖星周期估計(jì)誤差。因此，增加X射線探測(cè)器面積與觀測(cè)時(shí)間，有助于提高精度。

圖6 不同探測(cè)器面積和觀測(cè)時(shí)間下的周期估計(jì)誤差

為便于分析，本文在100 cm2、1 000 cm2、10 000 cm2這3種典型的探測(cè)器面積下研究估計(jì)誤差，如圖6(b)所示。可以看出，X射線探測(cè)器面積增大2個(gè)數(shù)量級(jí)，脈沖星周期估計(jì)精度提高約一個(gè)數(shù)量級(jí)。觀測(cè)時(shí)間延長(zhǎng)2個(gè)數(shù)量級(jí)，脈沖星周期估計(jì)精度提高約3個(gè)數(shù)量級(jí)。因此，相比于增加探測(cè)器面積，延長(zhǎng)觀測(cè)時(shí)間更加有利于提高估計(jì)精度。

值得一提的是，本文方法并未考慮延長(zhǎng)觀測(cè)時(shí)間所帶來的軌道外推誤差影響。這是因?yàn)楸疚姆椒ㄟm用于深空探測(cè)巡航段，此時(shí)的軌道模型近似于線性的[10]，軌道外推誤差的影響可忽略不計(jì)。若在環(huán)繞段長(zhǎng)時(shí)間觀測(cè)，利用天文測(cè)角信息補(bǔ)償[9]是一個(gè)較好的選擇。

5 結(jié) 論

本文提出了一種基于EMD-CS的脈沖星周期超快速估計(jì)方法。考慮到IMF包含了脈沖輪廓畸變度這一微弱局部特征，該方法利用畸變脈沖輪廓的IMF構(gòu)造測(cè)量矩陣，使得采樣率大幅下降，從而提高計(jì)算速度。該方法具有以下優(yōu)點(diǎn)：

1) 計(jì)算量小，約2 ms。究其原因，一方面，與FFT-CS類似，EMD-CS僅利用一個(gè)累積輪廓進(jìn)行估計(jì)，避免了多次脈沖信號(hào)累積；另一方面，與FFT-CS相比，測(cè)量矩陣尺寸減小了一個(gè)數(shù)量級(jí)以上，降低了采樣和匹配估計(jì)的計(jì)算量。

2) 估計(jì)精度高。在相同測(cè)量矩陣尺寸下，EMD-CS的估計(jì)精度優(yōu)于FFT-CS。若要二者達(dá)到相同精度，F(xiàn)FT-CS的測(cè)量矩陣比EMD-CS大一個(gè)數(shù)量級(jí)以上。

3) 純實(shí)數(shù)運(yùn)算。EMD無虛數(shù)運(yùn)算，IMF為實(shí)數(shù)。與FFT-CS相比，EMD-CS進(jìn)一步降低了計(jì)算復(fù)雜度。

4) 原信號(hào)長(zhǎng)度不受限制。哈達(dá)瑪矩陣的階數(shù)N應(yīng)滿足N、N/12、N/20必須是2的整數(shù)次冪。由于尺寸受限，實(shí)際難與脈沖輪廓長(zhǎng)度匹配，無法完成采樣。而IMF長(zhǎng)度必與原信號(hào)相等。

綜上，基于EMD-CS的脈沖星周期超快速估計(jì)方法計(jì)算量小，估計(jì)精度高，適合于深空探測(cè)在軌計(jì)算。

基于EMD-CS的脈沖星周期超快速估計(jì)方法適用于時(shí)不變的軌道模型。而環(huán)繞段、捕獲段等是時(shí)變的，如何在時(shí)變段實(shí)現(xiàn)脈沖星周期估計(jì)是一個(gè)值得研究的問題。