基于EMD-CS的脈沖星周期超快速估計

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

(1)

Φ0可以表示為

(2)

(3)

(4)

(5)

MEIMF的值為

(6)

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

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

2 EMD-CS

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

圖1 所提算法流程圖

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

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

hm=Gm?h

(7)

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

(8)

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

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

將第m個畸變脈沖輪廓循環(huán)移位構(gòu)成第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個畸變脈沖輪廓?？蓪D個子字典φm合成一個三維矩陣字典Ψ(N×P×MD維)。

2) 脈沖頻率初始偏置

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

(10)

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

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

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

在壓縮感知中，信號稀疏恢復(fù)是通過感知矩陣與測量矢量相匹配，根據(jù)最大匹配值實現(xiàn)的。感知矩陣T是最優(yōu)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)

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

(15)

(16)

(17)

3 復(fù)雜度分析

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

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

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

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

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

4 仿真實驗

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

4.1 仿真條件

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

表1 PSR B0531+21仿真條件

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

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

表2 第1次剔除IMF的誤差

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

表3 第2次剔除IMF的誤差

表4 第3次剔除IMF的誤差

表5 第4次剔除IMF的誤差

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

圖2 迭代剔除與隨機抽取的對比

4.3 有限等距性質(zhì)

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

圖3 測量矢量之間的距離

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

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

4.4 基于FFT-CS的脈沖星周期估計

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

表6 EMD-CS與FFT-CS的對比

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

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

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

值得一提的是，器載計算機的計算能力有限，因此，這2種算法在軌運行時間必將大幅增長。所以，縮短運行時間是有一定實際意義的。

4.5 光子流量密度

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

圖5 不同流量下的周期估計誤差

4.6 探測器面積與觀測時間

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

圖6 不同探測器面積和觀測時間下的周期估計誤差

為便于分析，本文在100 cm2、1 000 cm2、10 000 cm2這3種典型的探測器面積下研究估計誤差，如圖6(b)所示?？梢钥闯觯琗射線探測器面積增大2個數(shù)量級，脈沖星周期估計精度提高約一個數(shù)量級。觀測時間延長2個數(shù)量級，脈沖星周期估計精度提高約3個數(shù)量級。因此，相比于增加探測器面積，延長觀測時間更加有利于提高估計精度。

值得一提的是，本文方法并未考慮延長觀測時間所帶來的軌道外推誤差影響。這是因為本文方法適用于深空探測巡航段，此時的軌道模型近似于線性的[10]，軌道外推誤差的影響可忽略不計。若在環(huán)繞段長時間觀測，利用天文測角信息補償[9]是一個較好的選擇。

5 結(jié) 論

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

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

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

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

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

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

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