基于改進(jìn)CEEMDAN的微多普勒特征提取

金家偉，阮懷林，孫 兵

(國(guó)防科技大學(xué)電子對(duì)抗學(xué)院，安徽 合肥 230031)

0 引言

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

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

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

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

1 改進(jìn)的自適應(yīng)噪聲完備總體經(jīng)驗(yàn)?zāi)B(tài)分解

定義Ek(·)為通過(guò)EMD獲取第k個(gè)IMF的算子，k=1,…,K；M(·)為獲取目標(biāo)信號(hào)局部均值的算子；s(t)為目標(biāo)的雷達(dá)回波信號(hào)，則：

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

(1)

ICEEMDAN算法通過(guò)將模態(tài)估計(jì)替換為局部均值估計(jì)，并將其從原始信號(hào)中減去，從而減少了模態(tài)中存在的噪聲量。其算法的具體實(shí)現(xiàn)步驟如下：

步驟1 當(dāng)i=1,2,…,I時(shí)，通過(guò)EMD分別計(jì)算s(i)(t)=s(t)+β0E1(w(i)(t))的I個(gè)局部均值，得到第一個(gè)殘差為：

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

(2)

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

(3)

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

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

(4)

并得出第二個(gè)IMF分量為：

(5)

步驟4 當(dāng)k=3,…,K時(shí)，計(jì)算出第k個(gè)殘差為：

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

(6)

步驟5 計(jì)算得出第k個(gè)IMF分量為：

(7)

步驟6 對(duì)于下一個(gè)k，轉(zhuǎn)到步驟4。

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

(8)

2 基于ICEEMDAN的微動(dòng)周期提取方法

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

圖1 進(jìn)動(dòng)目標(biāo)模型Fig.1 Precession target model

在光學(xué)區(qū)，雷達(dá)目標(biāo)的整體散射特性通常可以等效為若干個(gè)散射中心的疊加。不失一般性，假設(shè)目標(biāo)的微動(dòng)是周期性進(jìn)動(dòng)，目標(biāo)等效為K個(gè)散射中心，各散射中心各向同性，雷達(dá)發(fā)射的電磁波為連續(xù)單頻波，載頻為f0。雷達(dá)向由i個(gè)散射點(diǎn)組成的目標(biāo)發(fā)射電磁波，不考慮目標(biāo)平動(dòng)帶來(lái)的影響，返回的基帶信號(hào)可表示為：

(9)

式(9)中，λ是雷達(dá)載波波長(zhǎng)，Ri是第i個(gè)散射中心到雷達(dá)原點(diǎn)的距離，σi(t)是第i個(gè)散射中心的散射強(qiáng)度，取決于目標(biāo)的運(yùn)動(dòng)姿態(tài)。則第i個(gè)散射中心的瞬時(shí)微多普勒頻率為：

(10)

(11)

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

(12)

得到幅值函數(shù)和相位函數(shù)為：

(13)

從而可以求出其瞬時(shí)頻率為：

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

(14)

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

(15)

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

(16)

式(16)可記為：

(17)

式(17)中，H(ω,t)精確地描述了回波信號(hào)的幅值隨時(shí)間、頻率的變化規(guī)律，通過(guò)對(duì)H(ω,t)的分析可進(jìn)一步提取出目標(biāo)的微多普勒特征。

3 仿真對(duì)比與分析

分別以ICEEMDAN算法和CEEMDAN算法對(duì)進(jìn)動(dòng)錐體目標(biāo)回波信號(hào)進(jìn)行分解，然后對(duì)得到的IMF進(jìn)行希爾伯特譜分析，最后對(duì)比兩種算法的提取效果，從而驗(yàn)證ICEEMDAN算法在微多普勒特征提取方面的性能提升。

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

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

圖2 基于CEEMDAN的估計(jì)結(jié)果(無(wú)噪聲)Fig.2 Estimation result based on CEEMDAN (no noise)

圖3 基于ICEEMDAN的估計(jì)結(jié)果(無(wú)噪聲)Fig.3 Estimation result based on ICEEMDAN (no noise)

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

本 刊 聲 明

中國(guó)知網(wǎng)發(fā)起設(shè)立的“學(xué)術(shù)不端文獻(xiàn)檢測(cè)中心”,其功能是以《中國(guó)學(xué)術(shù)文獻(xiàn)網(wǎng)絡(luò)出版總庫(kù)》和大量國(guó)際學(xué)術(shù)文獻(xiàn)為全文比對(duì)資源,輔助檢查抄襲、一稿多投、不當(dāng)署名、偽造、篡改等學(xué)術(shù)不端行為。我刊作為中國(guó)知網(wǎng)的合作單位,有義務(wù)為凈化學(xué)術(shù)空氣,制止學(xué)術(shù)不端行為作出貢獻(xiàn),請(qǐng)各位讀者、作者大力支持,與我們共同努力,從根本上鏟除學(xué)術(shù)腐敗的土壤,樹(shù)立全民求真、求實(shí)的科學(xué)態(tài)度。

本刊編輯部