基于改進CEEMDAN的微多普勒特征提取

金家偉，阮懷林，孫 兵

(國防科技大學電子對抗學院，安徽 合肥 230031)

0 引言

彈道目標在飛行時，為了保持飛行的穩定性，不僅會繞自身對稱軸作自旋運動，還會由于受到橫向作用力繞空間定向軸進行進動[1]。微動會對雷達回波產生頻率調制，從而產生微多普勒效應[2-3]。微多普勒頻率具有時變性，頻譜分析法無法展現頻率與時間的關系，因此，時頻分析法被廣泛應用于分析時變的微多普勒頻率，因為時頻分析法可以將一維信號映射到二維時頻平面來同時顯示信號的時頻信息[4]。文獻[5]提出一種基于參數化時頻分析的進動錐裙目標微多普勒曲線提取方法。文獻[6]結合Gabor時頻分布和變分模態分解，來估計目標的自旋頻率和錐旋頻率；文獻[7]利用短時分數階傅里葉變換分離肢體和軀干的微多普勒信號。然而，由于時頻分析法是基于基函數比較固定傅里葉變換理論，導致其缺乏自適應性或者自適應性差，在提取微多普勒信息時會存在一些不可避免的不足，難以同時保證時頻分辨率以及避免交叉項[8]。

經驗模態分解(empirical mode decomposition，EMD)沒有固定的基函數，可以將目標回波信號分解為不同尺度的本征模態函數(intrinsic mode functions,IMF)，這些IMF一般具有明顯的物理意義，表現了信號所包含的真實物理過程。在這過程中不存在人為干預，具有良好的自適應性[9]。但是在實際應用中，EMD算法并不完美，同樣存在一些問題，如模態混合、端點效應等。為了解決模態混合問題，文獻[10]通過加入不同有限振幅白噪聲來輔助數據分析，提出了總體經驗模態分解(ensemble empirical mode decomposition，EEMD)算法。

EEMD雖然大大減少了模態混疊現象，但其分解出的模態之和不能完美重構出目標信號，為了解決該問題，文獻[11]在每一階段的分解中都加入一個特定的白噪聲，再通過計算殘差得到該階段的模態，提出了自適應噪聲完備總體經驗模態分解(complete ensemble empirical mode decomposition with adaptive noise，CEEMDAN)算法，該算法實現了可以忽略的重建誤差，并解決了加入的噪聲不同時產生不同模態數量的問題。EMD及其衍生算法EEMD和CEEMDAN近些年來已廣泛應用于微多普勒特征分析的研究中。文獻[12]使用希爾伯特-黃變換對振動目標進行雷達微多普勒特征分析；文獻[13]提出了一種基于經驗模態分解的圓錐目標雷達微動特征提取方法；文獻[14]通過對多目標回波信號進行CEEMDAN分解，結合小波閾值去噪方法，對各本征模態函數(IMF)進行分層處理并累加，分離出了彈頭和碎片回波。

然而，CEEMDAN得到的IMF中仍含有一定的殘留噪聲，并且其中的信號信息出現得較晚，這些不足會極大地影響微多普勒特征提取的準確性和實效性。文獻[15]通過將估計模態替換為估計局部均值，并將真實模態定義為當前殘差與其局部均值的平均值之間的差值，提出了改進的自適應噪聲完備總體經驗模態分解(improved complete ensemble empirical mode decomposition with adaptive noise，ICEEMDAN)算法。本文針對CEEMDAN在微多普勒特征提取中存在的問題，提出基于ICEEMDAN的微多普勒特征提取方法。

1 改進的自適應噪聲完備總體經驗模態分解

定義Ek(·)為通過EMD獲取第k個IMF的算子，k=1,…,K；M(·)為獲取目標信號局部均值的算子；s(t)為目標的雷達回波信號，則：

E1(s(t))=s(t)-M(s(t))

(1)

ICEEMDAN算法通過將模態估計替換為局部均值估計，并將其從原始信號中減去，從而減少了模態中存在的噪聲量。其算法的具體實現步驟如下：

步驟1 當i=1,2,…,I時，通過EMD分別計算s(i)(t)=s(t)+β0E1(w(i)(t))的I個局部均值，得到第一個殘差為：

r1(t)=〈M(x(i)(t))〉

(2)

步驟2 計算第一個(k=1)IMF分量為：

(3)

步驟3 計算序列r1(t)+β1E2(w(i)(t))局部均值的平均值作為第二個殘差，即：

r2(t)=〈M(r1(t)+β1E2(w(i)(t)))〉

(4)

并得出第二個IMF分量為：

(5)

步驟4 當k=3,…,K時，計算出第k個殘差為：

rk(t)=〈M(rk-1(t)+βk-1Ek(w(i)(t)))〉

(6)

步驟5 計算得出第k個IMF分量為：

(7)

步驟6 對于下一個k，轉到步驟4。

重復步驟4到步驟6，直到獲得的殘差不能被進一步分解(即殘差滿足IMF的條件或者其極值小于等于2)。此時目標回波信號s(t)可以表示為：

(8)

2 基于ICEEMDAN的微動周期提取方法

空間進動錐體目標模型如圖1所示，坐標系(U,V,W)為雷達坐標系，雷達靜止于原點Q。O為目標質心，以O為原點，目標對稱軸Oz為z軸建立目標本體坐標系O-xyz。以O為原點建立參考坐標系O-XYZ，以初始時刻與目標對稱軸Oz、進動軸OZ共面且垂直于OZ的方向為Y軸，X軸根據右手準則確定。目標在平動的同時，以角速度ωs繞對稱軸z軸做自旋運動，同時以角速度ωc繞軸OZ做錐旋運動(ωs和ωc均采用參考坐標系中的表達式)，自旋軸和錐旋軸之間的夾角為進動角θ。

圖1 進動目標模型Fig.1 Precession target model

在光學區，雷達目標的整體散射特性通常可以等效為若干個散射中心的疊加。不失一般性，假設目標的微動是周期性進動，目標等效為K個散射中心，各散射中心各向同性，雷達發射的電磁波為連續單頻波，載頻為f0。雷達向由i個散射點組成的目標發射電磁波，不考慮目標平動帶來的影響，返回的基帶信號可表示為：

(9)

式(9)中，λ是雷達載波波長，Ri是第i個散射中心到雷達原點的距離，σi(t)是第i個散射中心的散射強度，取決于目標的運動姿態。則第i個散射中心的瞬時微多普勒頻率為：

(10)

(11)

式(11)中，p表示柯西主值，則其解析信號Zj(t)可被構造為：

(12)

得到幅值函數和相位函數為：

(13)

從而可以求出其瞬時頻率為：

ωj(t)=2πfj(t)=dθj(t)/dt

(14)

而由式(13)可看出，振幅和相位都是時間的函數，如果把振幅表示在時間—頻率平面上，第j個IMF分量的希爾伯特時頻譜為：

(15)

式(15)中，Re(·)是表示實部。目標回波信號s(t)的希爾伯特時頻譜為：

(16)

式(16)可記為：

(17)

式(17)中，H(ω,t)精確地描述了回波信號的幅值隨時間、頻率的變化規律，通過對H(ω,t)的分析可進一步提取出目標的微多普勒特征。

3 仿真對比與分析

分別以ICEEMDAN算法和CEEMDAN算法對進動錐體目標回波信號進行分解，然后對得到的IMF進行希爾伯特譜分析，最后對比兩種算法的提取效果，從而驗證ICEEMDAN算法在微多普勒特征提取方面的性能提升。

仿真條件：假設雷達工作在10 GHz，且雷達坐標系中本體坐標系原點O的坐標為(400，500，100) km，本地坐標系和參考坐標系之間的初始歐拉角(x-y-z序列)為(30°，60°，45°)。假設目標繞z軸旋轉，目標上有兩個散射點：第一散射點A位于錐頂，在本體坐標系中的坐標為(0，0，1) m；第二個散射點B位于錐底的尾翼，在本體坐標系中的坐標是(0.5，0，-0.5) m。雷達照射時間為8 s，旋轉頻率為fs=1 Hz，圓錐運動頻率為fc=0.5 Hz。

當不存在噪聲時，基于CEEMDAN和ICEEMDAN的微多普勒特征提取結果分別如圖2和圖3所示，其中的頻率都是取模后的頻率，也就是正頻率。圖2(a)和圖3(a)分別是兩種算法獲取的回波信號希爾伯特譜，展現了回波信號時頻分布的細節；圖2(b)和圖3(b)分別是兩種算法獲取的IMF前三個分量的希爾伯特譜。

圖2 基于CEEMDAN的估計結果(無噪聲)Fig.2 Estimation result based on CEEMDAN (no noise)

圖3 基于ICEEMDAN的估計結果(無噪聲)Fig.3 Estimation result based on ICEEMDAN (no noise)

由圖2和圖3可知，兩種算法都能夠準確提取出目標的微多普勒周期為2 s，這體現了兩種算法都具備提取微多普勒特征的功能。但從圖2中可以看到，在沒有噪聲的情況下，CEEMDAN的IMF中仍然殘留著部分噪聲。這些噪聲是在在CEEMDAN算法的實現過程中，主動加入信號進行輔助分解的，但該算法卻未能在分解過程中將所主動加入的噪聲從IMF中清理掉，這勢必會影響到微多普勒特征提取的準確性。而從圖3可以看到，ICEEMDAN的IMF中不存在這些殘留噪聲，表明ICEEMDAN克服了CEEMDAN的這一不足。

當雷達回波信號中混雜著噪聲時，同樣的仿真條件下，假設信噪比SNR=5 dB，圖4(a)為基于CEEMDAN的回波信號希爾伯特譜，圖4(b)是基于CEEMDAN得到的IMF前三個分量的希爾伯特譜。

圖4 基于CEEMDAN的估計結果(SNR=5 dB)Fig.4 Estimation result based on CEEMDAN (SNR=5 dB)

從圖4(a)中可以看出，回波信號希爾伯特譜中的微多普勒頻率伴隨著大量的噪聲頻率，難以從中提取微多普勒特征。但以CEEMDAN將信號分解后，噪聲主要集中于IMF第一個分量中，信號大部分被分解到IMF的第三個分量中，IMF第二個分量也包含了一部分的信號分量，如圖4(b)所示，只能從IMF分量中提取微多普勒特征。

在其他條件不變的情況下，改用ICEEMDAN，圖5(a)為回波信號的希爾伯特譜，圖5(b)為其IMF前三個分量的希爾伯特譜。

圖5 基于ICEEMDAN的估計結果(SNR=5 dB)Fig.5 Estimation result based on ICEEMDAN (SNR=5 dB)

與圖4(a)相似，同樣難以從圖5(a)提取微多普勒特征，此時也同樣只能從圖5(b)中的IMF分量中提取信息。但對比圖5(b)和圖4(b)，不難看出，CEEMDAN算法中的微多普勒信號大部分被分解在第三個分量中，而ICEEMDAN算法中的微多普勒信號則主要存在于第二個分量，即ICEEMDAN的信號信息出現得比CEEMDAN早。

進一步降低信噪比，圖6(a)和圖6(b)分別給出了信噪比為1 dB和0 dB時，基于ICEEMDAN得到的IMF前三個分量的希爾伯特譜。

圖6 ICEEMDAN前三個IMF分量的希爾伯特譜Fig.6 Hilbert spectrum of the first three IMFs components of ICEEMDAN

由圖6(b)可知，即使信噪比低至0 dB，仍能通過其第二個IMF分量的希爾伯特譜提取出微多普勒周期。這表明，在低信噪比環境中，ICEEMDAN具備良好的性能。

4 結論

本文提出基于ICEEMDAN的進動目標微多普勒特征提取方法，該方法以ICEEMDAN將進動目標回波信號分解為IMF，然后對IMF進行希爾伯特譜分析得到時頻譜，用以估計目標的微動周期。仿真實驗表明，該方法提取特征信息更快，實用性好，其IMF中的殘余噪聲較少，具備更好的抗噪性能。

本 刊 聲 明

中國知網發起設立的“學術不端文獻檢測中心”,其功能是以《中國學術文獻網絡出版總庫》和大量國際學術文獻為全文比對資源,輔助檢查抄襲、一稿多投、不當署名、偽造、篡改等學術不端行為。我刊作為中國知網的合作單位,有義務為凈化學術空氣,制止學術不端行為作出貢獻,請各位讀者、作者大力支持,與我們共同努力,從根本上鏟除學術腐敗的土壤,樹立全民求真、求實的科學態度。

本刊編輯部