基于CEMD的旋翼微動目標(biāo)雜波抑制方法

夏賽強(qiáng)，向虎，陳文峰，楊軍，陳一暢

空軍預(yù)警學(xué)院，武漢 430019

當(dāng)目標(biāo)中存在旋轉(zhuǎn)、振動、進(jìn)動和章動等微動形式時，雷達(dá)回波將產(chǎn)生除目標(biāo)剛體部分引起的多普勒頻率外的附加頻率調(diào)制，稱為微多普勒現(xiàn)象(Micro-Doppler Effect，MDE)。自2000年Chen[1]提出微多普勒概念以來，對微動目標(biāo)微多普勒現(xiàn)象的研究迅速發(fā)展起來，其中，旋翼旋轉(zhuǎn)作為典型的微動形式，受到國內(nèi)外學(xué)者的廣泛研究[2-4]。實(shí)際中，旋翼微動目標(biāo)的背景通常含有雜波，因此需要對雜波進(jìn)行抑制。目前針對旋翼微動目標(biāo)的雜波抑制方法大多在時域完成，如文獻(xiàn)[5]直接在多普勒域中利用頻域?yàn)V波的方法對雜波進(jìn)行了抑制，但是由于該方法采用的是理想濾波器，信號由頻域逆變換到時域時會出現(xiàn)拖尾效應(yīng)。文獻(xiàn)[6-7]利用延遲線技術(shù)建立動目標(biāo)顯示(Moving Target Indicator,MTI)濾波器，在時域?qū)π盘栠M(jìn)行對消，達(dá)到雜波抑制的效果。文獻(xiàn)[8-9]利用雜波相鄰距離、方位單元內(nèi)的相關(guān)性來進(jìn)行雜波對消，即求取時域回波的協(xié)方差矩陣，然后利用特征分解和主元分析法獲得地雜波的估計值，最后用原始數(shù)據(jù)減去估計值，達(dá)到雜波抑制的效果。

然而，上述時域雜波抑制方法存在回波信號能量損失大的缺點(diǎn)，針對這一問題，本文將回波信號進(jìn)行時頻分析，研究利用復(fù)數(shù)據(jù)經(jīng)驗(yàn)?zāi)B(tài)分解(CEMD)進(jìn)行雜波抑制，從能量的角度選擇微動信號的本征模態(tài)函數(shù)(IMF)分量對微動信號進(jìn)行重構(gòu)來抑制旋翼微動目標(biāo)回波中的地雜波，并與MTI雜波抑制方法進(jìn)行了對比。理論分析和仿真結(jié)果均表明該方法既能夠有效保留微動信號，又可以有效抑制雜波。

1 雜波背景下旋翼微動目標(biāo)信號建模分析

基于窄帶雷達(dá)，在觀測時間較短的情況下，旋翼目標(biāo)不會出現(xiàn)距離單元走動，若觀測時間較長，旋翼目標(biāo)出現(xiàn)距離單元走動現(xiàn)象，此時對旋翼目標(biāo)進(jìn)行運(yùn)動補(bǔ)償后，旋翼目標(biāo)的運(yùn)動狀態(tài)可等效為懸停狀態(tài)[10]。為簡化分析，假設(shè)目標(biāo)與雷達(dá)處于同一個平面(若不處同一平面，只需要將目標(biāo)投影到雷達(dá)所處平面上)，旋翼目標(biāo)模型如圖1所示，雷達(dá)到旋翼目標(biāo)旋轉(zhuǎn)中心C的距離為RC，假設(shè)旋翼微動目標(biāo)為散射點(diǎn)模型，旋翼某一個葉片上的散射點(diǎn)P到旋轉(zhuǎn)中心的距離為r(0≤r≤l,l為旋翼葉片長度)，點(diǎn)P到雷達(dá)的距離為RP，同時點(diǎn)P以角速度ω繞旋轉(zhuǎn)中心C旋轉(zhuǎn)，初始時刻點(diǎn)P的初始旋轉(zhuǎn)角為θ。

設(shè)雷達(dá)發(fā)射信號為窄帶線性調(diào)頻信號[2]

(1)

雜波背景下旋翼葉片上散射點(diǎn)的回波信號可表示為

(2)

式中：K為葉片數(shù)目；i為旋翼上的第i個葉片；N為目標(biāo)單個葉片上的散射點(diǎn)數(shù)目；j為葉片上的第j個散射點(diǎn)；σij為散射點(diǎn)的散射系數(shù)；c(t)為雜波；n(t)為噪聲；c為光速；Rij(tm)為旋翼上散射點(diǎn)到雷達(dá)的距離，即[2]

Rij(tm)=RC+rijcos(ωtm+θij)

(3)

雜波c(t)是具有高斯型雜波譜的地雜波，高斯型雜波歸一化功率譜密度函數(shù)[9]可表示為

(4)

討論具有高斯型功率譜的雜波譜寬，定義高斯型功率譜的3 dB寬度為雜波的譜寬Sw，即

(5)

從表1中雜波速度起伏展寬值σv的典型值和式(5)來看，地雜波的譜寬較窄，若雷達(dá)發(fā)射波長λ=0.3 m，那么Sw只有幾赫茲到幾十赫茲，并且在時頻域上，雜波能量和頻譜主要集中在零頻附近。

式(2)脈沖壓縮后的結(jié)果可表示為

(6)

表1 σv的典型值Table 1 Typical values of σv

且本文定義的信噪比(SNR)和信雜比(SCR)均為脈壓后的信噪比和信雜比。

由多普勒項(xiàng)可以看出，旋翼目標(biāo)微多普勒相位可表示為

(7)

由式(7)可知，散射點(diǎn)的微多普勒呈現(xiàn)余弦函數(shù)形式，而余弦函數(shù)的周期與旋翼旋轉(zhuǎn)周期是一致的。且不同葉片上散射點(diǎn)引起的微多普勒曲線有著不同的初相，而同一葉片散射點(diǎn)初相是一致的，差異在于多普勒頻率峰值大小的不同。對應(yīng)式(7)，可知旋翼目標(biāo)微多普勒頻率為

(8)

從式(8)中可以看出，當(dāng)sin(ωtm+θij)=1且rij=l時，旋翼目標(biāo)在該散射點(diǎn)的瞬時多普勒頻率最大，此時

(9)

假設(shè)葉片轉(zhuǎn)動速度ω=6π rad/s，葉片長度l=6 m，雷達(dá)發(fā)射信號波長λ=0.3 m，此時旋翼目標(biāo)瞬時微多普勒最大頻率fm-Dmax=754 Hz，并且其變化規(guī)律受到了正弦的調(diào)制。而根據(jù)表1中雜波速度起伏值知道地雜波的典型值均小于1 m/s，根據(jù)式(5)計算出雜波在時頻域的展寬最大值為11.1 Hz，并且分布在零頻附近，因此微動目標(biāo)信號與雜波在時頻域上存在展寬以及分布位置上的區(qū)別，因此可在時頻域上將雜波和微動信號區(qū)分開。

2 基于CEMD的雜波抑制方法原理

基于CEMD的旋翼微動目標(biāo)雜波抑制方法，首先,需要對脈壓后的回波進(jìn)行CEMD處理，得到不同的IMF分量，通過定義的選擇方式，選擇出高頻IMF分量在時域?qū)ξ有盘栠M(jìn)行重構(gòu);然后,對重構(gòu)的信號進(jìn)行時頻分析和過門限處理[12];最終,得到雜波抑制后的微動曲線。處理基本過程如圖2所示。

由于雜波和微動目標(biāo)信號在時頻域上存在明顯的展寬以及分布位置的區(qū)別，對于這種情況，可以利用CEMD方法將雜波和微動目標(biāo)信號分離開來。

傳統(tǒng)的經(jīng)驗(yàn)?zāi)B(tài)分解(EMD)可以將復(fù)雜的多分量實(shí)信號自適應(yīng)分解為若干個IMF分量的和，對每個IMF分量進(jìn)行Hilbert變換，可求出瞬時頻率和幅值，從而得到信號完整的時頻分布[13]。但是對于雷達(dá)回波而言，回波都是復(fù)信號形式，EMD方法無法對其進(jìn)行處理，文獻(xiàn)[14]針對EMD方法無法處理復(fù)信號的缺點(diǎn)，將EMD方法由實(shí)信號拓展到復(fù)信號，提出了CEMD方法。該方法采用提取零均值的旋轉(zhuǎn)分量來計算信號局部細(xì)節(jié)，并將信號分解為一系列從快變到慢變的旋轉(zhuǎn)矢量。實(shí)質(zhì)上，CEMD的分解過程就是一系列的濾波過程，將信號中存在的不同頻率或趨勢逐級分解開來，產(chǎn)生一系列具有不同頻率尺度的IMF，并且含高頻分量的IMF優(yōu)先被分離出來，也就是說，CEMD分解的過程就是一種從高頻到低頻的自適應(yīng)濾波過程。CEMD分解的示意圖如圖3所示，能有效地將包含了多個頻率分量的信號分解成頻率不同的單頻信號和一個余量[15]。

由分析可知，同時包含了微動目標(biāo)回波和雜波的雷達(dá)回波信號，如果直接用加窗的方法來分開低頻和高頻，將會對微動信號產(chǎn)生影響。但是從時頻域分析來看，其微動目標(biāo)回波受到旋轉(zhuǎn)運(yùn)動的調(diào)制，引起的微多普勒頻率是從負(fù)頻率到正頻率變化的寬譜，而地雜波通常只在零頻附近，存在明顯的區(qū)別，因此通過CEMD方法可以有效地將雜波和微動目標(biāo)回波區(qū)分開來，然后在時域和時頻域進(jìn)行微動目標(biāo)信號的重構(gòu)，從而達(dá)到雜波抑制的效果。CEMD方法具體流程[16]如圖4所示。

圖4中：x(t)為輸入信號,記為Si,k(t)；imfi,k-1為模態(tài)函數(shù)；Pφi,k(t)為在一系列方向φi,k(φi,k=kπ/K, 0≤k≤K)上對imfi,k-1進(jìn)行投影，并取實(shí)部；Sφi,k(t)為將實(shí)信號Pφi,k(t)的IMF在φi,k上投影得到的復(fù)IMF。

通常，在CEMD分解得到了多個IMF分量之后，通過計算各IMF分量與回波信號的相關(guān)系數(shù)來判斷哪個IMF分量是主雜波，并認(rèn)為主雜波之后的IMF分量均為雜波分量[15]。基于該方法運(yùn)算量大的缺點(diǎn)，本文將CEMD分離后的各個模態(tài)函數(shù)IMF(i)取出，分別求取各個模態(tài)函數(shù)內(nèi)信號的能量，并求出各個模態(tài)函數(shù)占全部模態(tài)函數(shù)總能量的比重T

(10)

由第1節(jié)分析可知，雜波只存在于低頻部分，而目標(biāo)信號不僅僅局限于低頻部分，還存在大量的高頻分量。進(jìn)行CEMD處理后，根據(jù)IMF(i)從高頻到低頻的排列方式可對信號進(jìn)行頻率區(qū)分，在求出各分量能量比重后，其比重T突變較快的位置都是可能的雜波和目標(biāo)在頻率上的界限，當(dāng)比重T在高頻區(qū)域急劇下降或者在低頻區(qū)域急劇上升的位置即可判定為雜波和目標(biāo)的界限，在本文方法中，這個界限的經(jīng)驗(yàn)值通常取IMF(5)。將這個分界IMF作為選擇標(biāo)準(zhǔn)來區(qū)分微動信號IMF分量和雜波IMF分量，然后將占據(jù)能量比重較高的低頻雜波IMF分量去除，僅保留具有一定能量的高頻微動信號IMF分量對微動目標(biāo)信號進(jìn)行重構(gòu)，即對符合標(biāo)準(zhǔn)的IMF分量求和。重構(gòu)得到的微動時域信號為

(11)

式中：IMF(j)為選擇后的高頻模態(tài)分量。對重構(gòu)的微動時域信號進(jìn)行Gabor變換可以重構(gòu)出微動目標(biāo)信號的時頻圖，Gabor變換屬于加窗短時傅立葉變換，是短時Fourier變換中當(dāng)窗函數(shù)取為高斯函數(shù)時的一種特殊情況，由于高斯函數(shù)滿足測不準(zhǔn)原理，因此，在所有可能的窗口函數(shù)中，高斯窗口函數(shù)能得到最好的時頻效果。其表達(dá)式[17]為

(12)

從上述分析可知，通過分析目標(biāo)微動信號與雜波的特性和區(qū)別，利用CEMD方法分離出代表雜波和微動目標(biāo)的各個模態(tài)函數(shù)，再通過微動目標(biāo)的模態(tài)函數(shù)重構(gòu)微動信號，達(dá)到抑制雜波、保留微動信號的目的。

3 仿真分析

為了更加清晰得到本文方法對雜波的抑制效果，這里選擇3個旋翼葉片頂點(diǎn)和1個旋轉(zhuǎn)中心點(diǎn)作為旋翼散射點(diǎn)，模型如圖5所示，其葉片長度為l=6 m，旋翼轉(zhuǎn)速ω=6π rad/s。且本文定義的SCR和SNR均為脈壓后的信雜比和信噪比。具體仿真參數(shù)如表2所示。

參數(shù)PRF/HzTp/μsB/MHzfs/MHzTa/sσv/(m·s-1)數(shù)值4 000100120.5120.32

3.1 回波建模仿真

回波建模仿真如圖6所示，其中圖6(a)為脈壓后SNR=5 dB時的回波仿真結(jié)果；圖6(b)為在圖6(a)中添加雜波，使脈壓后SCR=0 dB時的回波；圖6(c)為圖6(b)的距離多普勒表示。

從圖6可以看出，在窄帶條件下，旋翼目標(biāo)回波信號脈壓后均在相同距離單元內(nèi)，不會產(chǎn)生距離單元走動，在所有距離單元內(nèi)都存在雜波和噪聲，仿真的地雜波在方位向是相關(guān)且是慢起伏的。

3.2 雜波抑制性能仿真與分析

雜波抑制性能仿真如圖7所示，其中圖7(a)為SCR=0 dB時，旋翼目標(biāo)回波信號添加雜波后的脈壓結(jié)果；圖7(b)為對圖7(a)中目標(biāo)所在的距離單元進(jìn)行時頻分析的結(jié)果；圖7(c)為進(jìn)行CEMD分解后各個IMF所占總能量的比重T之間的關(guān)系；圖7(d)為雜波抑制后對所有距離單元在時域重構(gòu)的結(jié)果；圖7(e)為對圖7(d)中微動目標(biāo)所在的距離單元進(jìn)行時頻分析以及過門限處理的結(jié)果；圖7(f)為對圖7(e)進(jìn)行過門限處理后的結(jié)果。

從圖7可以看出，在SCR=0 dB時，添加雜波后對目標(biāo)所在的距離單元進(jìn)行時頻分析可以清晰完整地得到時頻圖。其中微動信號的時頻曲線呈現(xiàn)余弦變化規(guī)律，并且從負(fù)頻率到正頻率占據(jù)了很大的譜寬，其值與葉片轉(zhuǎn)速和長度以及雷達(dá)工作波長有關(guān)，而雜波的時頻曲線只是一條在零頻附近占據(jù)了很小一段譜寬的直線，并且在整個時頻域都充滿了噪聲。得到各個IMF所占總能量的比重后，以IMF(5)為界，可以判斷IMF(1)～I(xiàn)MF(4)是旋翼目標(biāo)信號和噪聲的高頻IMF，IMF(6)～I(xiàn)MF(11)是包含雜波、噪聲和旋翼目標(biāo)信號的低頻部分的低頻IMF。在進(jìn)行CEMD分離后，利用本文方法選擇出的高頻IMF可以很好地在時域重構(gòu)出旋翼目標(biāo)信號，有效地進(jìn)行了雜波抑制。并且重構(gòu)出的旋翼目標(biāo)信號與無雜波情況的旋翼目標(biāo)信號幅度類似，并且對信噪比有一定的提升，這是因?yàn)樵谶M(jìn)行CEMD后，噪聲的低頻部分連同雜波可一起被抑制。對重構(gòu)出的時域旋翼目標(biāo)信號進(jìn)行時頻分析，可以有效地得到雜波抑制后的旋翼目標(biāo)信號時頻圖，但是由于僅取出了高頻IMF，導(dǎo)致了旋翼目標(biāo)信號時頻圖在低頻部分有缺失。使用過門限處理方法和骨架提取方法對旋翼目標(biāo)信號重構(gòu)的時頻圖進(jìn)行處理可以得到很清晰的骨架圖像，能夠較為準(zhǔn)確地讀取微動信息，此時旋翼目標(biāo)瞬時微多普勒最大頻率為759 Hz，與理論值754 Hz比較接近。

3.3 雜波抑制性能對比仿真與分析

為了直觀地對本文方法與文獻(xiàn)[6]中所提方法(簡稱方法1)進(jìn)行雜波抑制性能對比，仿真分為兩組進(jìn)行，第1組是在無噪聲情況下進(jìn)行的，此時SCR=-54 dB；第2組是在噪聲背景下進(jìn)行的，此時SCR=-54 dB，SNR=5 dB。圖8和圖9分別為第1組仿真和第2組仿真的結(jié)果，其中圖8(a)和圖9(a)為方法1雜波抑制后的時域?qū)οY(jié)果，圖8(b)和圖9(b)為對目標(biāo)所在的距離單元進(jìn)行時頻分析的結(jié)果，圖8(c)和圖9(c)為本文方法雜波抑制后對全部距離單元在時域進(jìn)行重構(gòu)的結(jié)果，圖8(d)和圖9(d)為本文方法雜波抑制后對微動信號在時頻域進(jìn)行重構(gòu)的結(jié)果。

從仿真結(jié)果可以看到，經(jīng)本文方法CEMD雜波抑制后，只存在較少的雜波剩余，并且對信噪比有一定提升，這是由于本文方法在抑制雜波的同時對低頻部分的噪聲同樣具有抑制效果，但是經(jīng)方法1對雜波進(jìn)行對消后，信號幅度明顯變小且存在明顯雜波剩余，并且利用本文方法在時頻域重構(gòu)出的微動信號時頻曲線基本保持完整，僅在低頻部分出現(xiàn)了少許缺失。而利用方法1進(jìn)行雜波抑制后，目標(biāo)所在距離單元的時頻分析結(jié)果不僅丟失了微動信號時頻曲線低頻部分的信息，還對不同頻率位置的微動曲線有削弱和放大的效果，并且對高頻部分的噪聲放大效果明顯，信噪比惡化嚴(yán)重，這是由于構(gòu)建的方法1中的濾波器幅頻特性導(dǎo)致的。

對比兩種雜波抑制方法，在無噪聲和有噪聲情況下，本文方法均能夠很好地抑制雜波，雜波剩余較少，并且能夠重構(gòu)出較為完整的旋翼微動信號，而利用方法1在進(jìn)行雜波抑制時存在明顯雜波剩余情況，并且丟失了大量低頻附近的信號，只能重構(gòu)出部分旋翼微動信號。因此，本文方法在雜波抑制時性能明顯優(yōu)于方法1。

4 結(jié) 論

1) 雷達(dá)回波信號中的雜波具有明顯不同于旋翼微動目標(biāo)信號的時頻特征，并且時頻分析方法常用來處理旋翼微動目標(biāo)信號，因此可以在時頻域?qū)﹄s波進(jìn)行抑制。

2) 基于CEMD的微動目標(biāo)雜波抑制方法，在對雜波進(jìn)行了抑制的同時，還對噪聲具有一定的抑制效果，能夠有效提高雷達(dá)回波信號中旋翼微動目標(biāo)信號的信雜比和信噪比。

3) 基于CEMD的微動目標(biāo)雜波抑制方法，通過改變IMF分量的選擇標(biāo)準(zhǔn)，在時域和時頻域?qū)π砦幽繕?biāo)信號進(jìn)行重構(gòu)，避免了傳統(tǒng)雜波抑制方法對回波信號能量的大量損失，與現(xiàn)有雜波抑制方法相比具有更好地抑制雜波能力和信號重構(gòu)效果。