金家偉,阮懷林,孫 兵
(國防科技大學(xué)電子對抗學(xué)院,安徽 合肥 230031)
彈道目標(biāo)在飛行時,為了保持飛行的穩(wěn)定性,不僅會繞自身對稱軸作自旋運動,還會由于受到橫向作用力繞空間定向軸進(jìn)行進(jìn)動[1]。微動會對雷達(dá)回波產(chǎn)生頻率調(diào)制,從而產(chǎn)生微多普勒效應(yīng)[2-3]。微多普勒頻率具有時變性,頻譜分析法無法展現(xiàn)頻率與時間的關(guān)系,因此,時頻分析法被廣泛應(yīng)用于分析時變的微多普勒頻率,因為時頻分析法可以將一維信號映射到二維時頻平面來同時顯示信號的時頻信息[4]。文獻(xiàn)[5]提出一種基于參數(shù)化時頻分析的進(jìn)動錐裙目標(biāo)微多普勒曲線提取方法。文獻(xiàn)[6]結(jié)合Gabor時頻分布和變分模態(tài)分解,來估計目標(biāo)的自旋頻率和錐旋頻率;文獻(xiàn)[7]利用短時分?jǐn)?shù)階傅里葉變換分離肢體和軀干的微多普勒信號。然而,由于時頻分析法是基于基函數(shù)比較固定傅里葉變換理論,導(dǎo)致其缺乏自適應(yīng)性或者自適應(yīng)性差,在提取微多普勒信息時會存在一些不可避免的不足,難以同時保證時頻分辨率以及避免交叉項[8]。
經(jīng)驗?zāi)B(tài)分解(empirical mode decomposition,EMD)沒有固定的基函數(shù),可以將目標(biāo)回波信號分解為不同尺度的本征模態(tài)函數(shù)(intrinsic mode functions,IMF),這些IMF一般具有明顯的物理意義,表現(xiàn)了信號所包含的真實物理過程。在這過程中不存在人為干預(yù),具有良好的自適應(yīng)性[9]。但是在實際應(yīng)用中,EMD算法并不完美,同樣存在一些問題,如模態(tài)混合、端點效應(yīng)等。為了解決模態(tài)混合問題,文獻(xiàn)[10]通過加入不同有限振幅白噪聲來輔助數(shù)據(jù)分析,提出了總體經(jīng)驗?zāi)B(tài)分解(ensemble empirical mode decomposition,EEMD)算法。
EEMD雖然大大減少了模態(tài)混疊現(xiàn)象,但其分解出的模態(tài)之和不能完美重構(gòu)出目標(biāo)信號,為了解決該問題,文獻(xiàn)[11]在每一階段的分解中都加入一個特定的白噪聲,再通過計算殘差得到該階段的模態(tài),提出了自適應(yīng)噪聲完備總體經(jīng)驗?zāi)B(tài)分解(complete ensemble empirical mode decomposition with adaptive noise,CEEMDAN)算法,該算法實現(xiàn)了可以忽略的重建誤差,并解決了加入的噪聲不同時產(chǎn)生不同模態(tài)數(shù)量的問題。EMD及其衍生算法EEMD和CEEMDAN近些年來已廣泛應(yīng)用于微多普勒特征分析的研究中。文獻(xiàn)[12]使用希爾伯特-黃變換對振動目標(biāo)進(jìn)行雷達(dá)微多普勒特征分析;文獻(xiàn)[13]提出了一種基于經(jīng)驗?zāi)B(tài)分解的圓錐目標(biāo)雷達(dá)微動特征提取方法;文獻(xiàn)[14]通過對多目標(biāo)回波信號進(jìn)行CEEMDAN分解,結(jié)合小波閾值去噪方法,對各本征模態(tài)函數(shù)(IMF)進(jìn)行分層處理并累加,分離出了彈頭和碎片回波。
然而,CEEMDAN得到的IMF中仍含有一定的殘留噪聲,并且其中的信號信息出現(xiàn)得較晚,這些不足會極大地影響微多普勒特征提取的準(zhǔn)確性和實效性。文獻(xiàn)[15]通過將估計模態(tài)替換為估計局部均值,并將真實模態(tài)定義為當(dāng)前殘差與其局部均值的平均值之間的差值,提出了改進(jìn)的自適應(yīng)噪聲完備總體經(jīng)驗?zāi)B(tài)分解(improved complete ensemble empirical mode decomposition with adaptive noise,ICEEMDAN)算法。本文針對CEEMDAN在微多普勒特征提取中存在的問題,提出基于ICEEMDAN的微多普勒特征提取方法。
定義Ek(·)為通過EMD獲取第k個IMF的算子,k=1,…,K;M(·)為獲取目標(biāo)信號局部均值的算子;s(t)為目標(biāo)的雷達(dá)回波信號,則:
E1(s(t))=s(t)-M(s(t))
(1)
ICEEMDAN算法通過將模態(tài)估計替換為局部均值估計,并將其從原始信號中減去,從而減少了模態(tài)中存在的噪聲量。其算法的具體實現(xiàn)步驟如下:
步驟1 當(dāng)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 當(dāng)k=3,…,K時,計算出第k個殘差為:
rk(t)=〈M(rk-1(t)+βk-1Ek(w(i)(t)))〉
(6)
步驟5 計算得出第k個IMF分量為:
(7)
步驟6 對于下一個k,轉(zhuǎn)到步驟4。
重復(fù)步驟4到步驟6,直到獲得的殘差不能被進(jìn)一步分解(即殘差滿足IMF的條件或者其極值小于等于2)。此時目標(biāo)回波信號s(t)可以表示為:
(8)
空間進(jìn)動錐體目標(biāo)模型如圖1所示,坐標(biāo)系(U,V,W)為雷達(dá)坐標(biāo)系,雷達(dá)靜止于原點Q。O為目標(biāo)質(zhì)心,以O(shè)為原點,目標(biāo)對稱軸Oz為z軸建立目標(biāo)本體坐標(biāo)系O-xyz。以O(shè)為原點建立參考坐標(biāo)系O-XYZ,以初始時刻與目標(biāo)對稱軸Oz、進(jìn)動軸OZ共面且垂直于OZ的方向為Y軸,X軸根據(jù)右手準(zhǔn)則確定。目標(biāo)在平動的同時,以角速度ωs繞對稱軸z軸做自旋運動,同時以角速度ωc繞軸OZ做錐旋運動(ωs和ωc均采用參考坐標(biāo)系中的表達(dá)式),自旋軸和錐旋軸之間的夾角為進(jìn)動角θ。
圖1 進(jìn)動目標(biāo)模型Fig.1 Precession target model
在光學(xué)區(qū),雷達(dá)目標(biāo)的整體散射特性通常可以等效為若干個散射中心的疊加。不失一般性,假設(shè)目標(biāo)的微動是周期性進(jìn)動,目標(biāo)等效為K個散射中心,各散射中心各向同性,雷達(dá)發(fā)射的電磁波為連續(xù)單頻波,載頻為f0。雷達(dá)向由i個散射點組成的目標(biāo)發(fā)射電磁波,不考慮目標(biāo)平動帶來的影響,返回的基帶信號可表示為:
(9)
式(9)中,λ是雷達(dá)載波波長,Ri是第i個散射中心到雷達(dá)原點的距離,σi(t)是第i個散射中心的散射強(qiáng)度,取決于目標(biāo)的運動姿態(tài)。則第i個散射中心的瞬時微多普勒頻率為:
(10)
(11)
式(11)中,p表示柯西主值,則其解析信號Zj(t)可被構(gòu)造為:
(12)
得到幅值函數(shù)和相位函數(shù)為:
(13)
從而可以求出其瞬時頻率為:
ωj(t)=2πfj(t)=dθj(t)/dt
(14)
而由式(13)可看出,振幅和相位都是時間的函數(shù),如果把振幅表示在時間—頻率平面上,第j個IMF分量的希爾伯特時頻譜為:
(15)
式(15)中,Re(·)是表示實部。目標(biāo)回波信號s(t)的希爾伯特時頻譜為:
(16)
式(16)可記為:
(17)
式(17)中,H(ω,t)精確地描述了回波信號的幅值隨時間、頻率的變化規(guī)律,通過對H(ω,t)的分析可進(jìn)一步提取出目標(biāo)的微多普勒特征。
分別以ICEEMDAN算法和CEEMDAN算法對進(jìn)動錐體目標(biāo)回波信號進(jìn)行分解,然后對得到的IMF進(jìn)行希爾伯特譜分析,最后對比兩種算法的提取效果,從而驗證ICEEMDAN算法在微多普勒特征提取方面的性能提升。
仿真條件:假設(shè)雷達(dá)工作在10 GHz,且雷達(dá)坐標(biāo)系中本體坐標(biāo)系原點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)上有兩個散射點:第一散射點A位于錐頂,在本體坐標(biāo)系中的坐標(biāo)為(0,0,1) m;第二個散射點B位于錐底的尾翼,在本體坐標(biāo)系中的坐標(biāo)是(0.5,0,-0.5) m。雷達(dá)照射時間為8 s,旋轉(zhuǎn)頻率為fs=1 Hz,圓錐運動頻率為fc=0.5 Hz。
當(dāng)不存在噪聲時,基于CEEMDAN和ICEEMDAN的微多普勒特征提取結(jié)果分別如圖2和圖3所示,其中的頻率都是取模后的頻率,也就是正頻率。圖2(a)和圖3(a)分別是兩種算法獲取的回波信號希爾伯特譜,展現(xiàn)了回波信號時頻分布的細(xì)節(jié);圖2(b)和圖3(b)分別是兩種算法獲取的IMF前三個分量的希爾伯特譜。
圖2 基于CEEMDAN的估計結(jié)果(無噪聲)Fig.2 Estimation result based on CEEMDAN (no noise)
圖3 基于ICEEMDAN的估計結(jié)果(無噪聲)Fig.3 Estimation result based on ICEEMDAN (no noise)
由圖2和圖3可知,兩種算法都能夠準(zhǔn)確提取出目標(biāo)的微多普勒周期為2 s,這體現(xiàn)了兩種算法都具備提取微多普勒特征的功能。但從圖2中可以看到,在沒有噪聲的情況下,CEEMDAN的IMF中仍然殘留著部分噪聲。這些噪聲是在在CEEMDAN算法的實現(xiàn)過程中,主動加入信號進(jìn)行輔助分解的,但該算法卻未能在分解過程中將所主動加入的噪聲從IMF中清理掉,這勢必會影響到微多普勒特征提取的準(zhǔn)確性。而從圖3可以看到,ICEEMDAN的IMF中不存在這些殘留噪聲,表明ICEEMDAN克服了CEEMDAN的這一不足。
當(dāng)雷達(dá)回波信號中混雜著噪聲時,同樣的仿真條件下,假設(shè)信噪比SNR=5 dB,圖4(a)為基于CEEMDAN的回波信號希爾伯特譜,圖4(b)是基于CEEMDAN得到的IMF前三個分量的希爾伯特譜。
圖4 基于CEEMDAN的估計結(jié)果(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的估計結(jié)果(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的信號信息出現(xiàn)得比CEEMDAN早。
進(jìn)一步降低信噪比,圖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分量的希爾伯特譜提取出微多普勒周期。這表明,在低信噪比環(huán)境中,ICEEMDAN具備良好的性能。
本文提出基于ICEEMDAN的進(jìn)動目標(biāo)微多普勒特征提取方法,該方法以ICEEMDAN將進(jìn)動目標(biāo)回波信號分解為IMF,然后對IMF進(jìn)行希爾伯特譜分析得到時頻譜,用以估計目標(biāo)的微動周期。仿真實驗表明,該方法提取特征信息更快,實用性好,其IMF中的殘余噪聲較少,具備更好的抗噪性能。
本 刊 聲 明
中國知網(wǎng)發(fā)起設(shè)立的“學(xué)術(shù)不端文獻(xiàn)檢測中心”,其功能是以《中國學(xué)術(shù)文獻(xiàn)網(wǎng)絡(luò)出版總庫》和大量國際學(xué)術(shù)文獻(xiàn)為全文比對資源,輔助檢查抄襲、一稿多投、不當(dāng)署名、偽造、篡改等學(xué)術(shù)不端行為。我刊作為中國知網(wǎng)的合作單位,有義務(wù)為凈化學(xué)術(shù)空氣,制止學(xué)術(shù)不端行為作出貢獻(xiàn),請各位讀者、作者大力支持,與我們共同努力,從根本上鏟除學(xué)術(shù)腐敗的土壤,樹立全民求真、求實的科學(xué)態(tài)度。
本刊編輯部