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

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

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]是一個較好的選擇。

5 結 論

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

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

2) 估計精度高。在相同測量矩陣尺寸下，EMD-CS的估計精度優于FFT-CS。若要二者達到相同精度，FFT-CS的測量矩陣比EMD-CS大一個數量級以上。

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

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

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

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