基于STFT和EMD的多普勒信號(hào)分離算法

楊彥利，馬 德，權(quán)建峰

(1.天津工業(yè)大學(xué)光電檢測技術(shù)與系統(tǒng)重點(diǎn)實(shí)驗(yàn)室，天津 300387；2.西安機(jī)電信息技術(shù)研究所，陜西 西安 710065)

0 引言

多普勒效應(yīng)揭示了反射波的頻率與物體運(yùn)動(dòng)的關(guān)系[1]。多普勒信號(hào)攜帶了目標(biāo)信息，通過檢測多普勒信號(hào)可以實(shí)現(xiàn)對(duì)目標(biāo)的探測。多普勒雷達(dá)、多普勒引信等就是利用了多普勒信號(hào)來探測和感知目標(biāo)。多普勒信號(hào)的可靠檢測與有效識(shí)別已成為該體制雷達(dá)和引信的核心技術(shù)。

當(dāng)波源離目標(biāo)很遠(yuǎn)時(shí)，可將目標(biāo)看成點(diǎn)目標(biāo)，此時(shí)傅里葉變換等方法可實(shí)現(xiàn)對(duì)多普勒頻率的估計(jì)；而當(dāng)波源與目標(biāo)逐漸靠近時(shí)，需用多點(diǎn)目標(biāo)模型來描述多普勒回波信號(hào)[2]，尤其是在近場區(qū)，需要考慮目標(biāo)的體效應(yīng)[3]。文獻(xiàn)[4]通過對(duì)實(shí)測信號(hào)的分析，指出目標(biāo)效應(yīng)和噪聲干擾對(duì)彈目交會(huì)信號(hào)有較大影響。由于目標(biāo)和波源的相對(duì)運(yùn)動(dòng)，多普勒信號(hào)的幅度和頻率具有時(shí)變性。在足夠低的信噪比下，只依靠幅度特征已不能有效地對(duì)目標(biāo)信號(hào)進(jìn)行檢測和提取[5]。文獻(xiàn)[6]研究了經(jīng)驗(yàn)?zāi)B(tài)分解(EMD)[7]對(duì)多普勒信號(hào)的提取技術(shù)。文獻(xiàn)[8]提出了一種基于快速傅里葉變換(FFT)和自適應(yīng)濾波并結(jié)合Hilbert變換來估計(jì)雷達(dá)多普勒信號(hào)頻率的方法?；跁r(shí)頻分布，文獻(xiàn)[9]對(duì)低信噪比下雷達(dá)回波信號(hào)瞬時(shí)頻率的估計(jì)進(jìn)行了有益的嘗試。在時(shí)頻分布上，利用峰值檢測與追蹤有助于信號(hào)瞬時(shí)頻率的估計(jì)[10-11]。利用局部極值點(diǎn)，通過設(shè)計(jì)合適的濾波器可以將復(fù)合信號(hào)分解成多個(gè)窄帶的子信號(hào)[12]，文獻(xiàn)[13]則基于Wigner-Ville分布提出了一種復(fù)合信號(hào)的盲成分分離算法。

近場區(qū)多普勒回波信號(hào)受體目標(biāo)效應(yīng)的影響而表現(xiàn)為多個(gè)多普勒成分的疊加，準(zhǔn)確分離和提取這些成分是進(jìn)行多普勒頻率估計(jì)的關(guān)鍵。近場區(qū)多普勒信號(hào)的識(shí)別和瞬時(shí)頻率的估計(jì)目前仍是一個(gè)有待解決的問題。本文針對(duì)此問題，提出了基于短時(shí)傅里葉變換(STFT)和EMD的多普勒信號(hào)分離算法。

1 多普勒回波信號(hào)模型

在遠(yuǎn)場區(qū)可以用點(diǎn)目標(biāo)模型來描述多普勒回波信號(hào)，可將其表示為[3]

ud=adcos(2πfdt+φd)

(1)

式(1)中，ad為多普勒信號(hào)幅值，fd為多普勒信號(hào)頻率，φd為多普勒信號(hào)的初始相位。ad與接收機(jī)的探測靈敏度SA、發(fā)射波波長λ、天線的方向性系數(shù)D、天線的方向圖F(φ)、目標(biāo)的雷達(dá)截面積σ和波源與目標(biāo)的距離R等因素有關(guān)，用公式表示為[3]：

(2)

由于式(2)中多個(gè)參數(shù)隨時(shí)間變化，導(dǎo)致多普勒信號(hào)的幅值具有時(shí)變性。對(duì)于無線電引信而言，由于發(fā)射機(jī)和接收機(jī)處于同一個(gè)載體，故fd可表示為：

fd=±2vr/λ

(3)

式(3)中，vr表示發(fā)射機(jī)和接收機(jī)之間的接近速度。在遠(yuǎn)場區(qū)，多普勒頻率基本不變[14]。

當(dāng)彈目距離逐漸靠近時(shí)，需要用多點(diǎn)目標(biāo)模型來描述多普勒信號(hào)，即

(4)

不同散射點(diǎn)反射的多普勒信號(hào)在幅值和頻率上會(huì)有差異，導(dǎo)致整個(gè)多普勒回波信號(hào)在頻譜上具有一定的帶寬。

在近場區(qū)，體目標(biāo)效應(yīng)會(huì)影響多普勒回波信號(hào)。目標(biāo)的不同形狀，甚至是目標(biāo)不同部位會(huì)具有不同的雷達(dá)截面積，再考慮到目標(biāo)的不同部位導(dǎo)致距離R不同，這些因素會(huì)造成反射的多普勒信號(hào)幅值出現(xiàn)變化。文獻(xiàn)[14]也指出近距離時(shí)由于體目標(biāo)效應(yīng)顯著，多普勒信號(hào)幅值會(huì)發(fā)生畦變。圖1所示是某實(shí)測信號(hào)的STFT時(shí)頻分布，將STFT變換后的幅值沿時(shí)間軸累加求和，得到在不同時(shí)刻回波信號(hào)的幅值變化曲線，如圖2所示。由該圖可知，在體目標(biāo)效應(yīng)的影響下，多普勒幅值總體上是隨著彈目距離接近而增大、反之則減小，但具有很大的波動(dòng)性。李向軍等人[15]建立了近炸引信終點(diǎn)彈道多普勒信號(hào)仿真平臺(tái)，從其仿真的信號(hào)以及對(duì)比的實(shí)測信號(hào)也可以發(fā)現(xiàn)：從時(shí)域上看，在接近目標(biāo)過程中多普勒信號(hào)的幅值總體上是增加的，但有明顯的波動(dòng)。

圖1 某實(shí)測多普勒信號(hào)的時(shí)頻分布Fig.1 Time-frequency distribution of a measured Doppler signal

圖2 STFT變換后幅值沿時(shí)間軸的累加和Fig.2 Cumulative sum of amplitudes along the time axis after STFT transformation

對(duì)于目標(biāo)的不同散射點(diǎn)，彈目的接近速度會(huì)有所差異，而由式(3)描述的多普勒頻率可知，在同一時(shí)刻的回波信號(hào)中存在多個(gè)多普勒成分。近場區(qū)，式(4)也可以用來描述多普勒信號(hào)，不過需要很多點(diǎn)來描述回波信號(hào)。通過圖1也可看到近場區(qū)多普勒回波信號(hào)具有一定的帶寬，這為識(shí)別多普勒信號(hào)和估計(jì)多普勒頻率帶來一些挑戰(zhàn)。

2 多普勒信號(hào)的分離

根據(jù)多普勒信號(hào)的特點(diǎn)，設(shè)計(jì)了一個(gè)基于STFT和EMD的多普勒信號(hào)分離算法，流程圖如圖3所示。該算法包括如下幾個(gè)步驟：首先對(duì)采集的多普勒回波信號(hào)進(jìn)行STFT分析；然后將STFT處理后得到的系數(shù)沿頻率軸進(jìn)行投影；再根據(jù)極值法進(jìn)行濾波處理得到子信號(hào)；將子信號(hào)進(jìn)行EMD處理得到窄帶的多普勒信號(hào)。

眾所周知，給定一個(gè)時(shí)域信號(hào)x(t)，其STFT變換可表示為

(5)

式(5)中，h(τ)是窗函數(shù)。計(jì)算機(jī)處理過程中采用離散STFT變換，得到系數(shù)X(tn,fn)。

圖3 基于STFT和EMD的多普勒信號(hào)分離算法流程圖Fig.3 Flow chart of Doppler signal separation algorithm based on STFT and EMD

基于Wigner-Ville分布，可以實(shí)現(xiàn)對(duì)信號(hào)不同成分的分離，但Wigner-Ville分布存在交叉項(xiàng)問題。由于STFT變換沒有交叉項(xiàng)，再考慮到在很短一個(gè)采樣周期內(nèi)，多普勒信號(hào)的頻率變化較小，因此，這里采用對(duì)STFT的變換結(jié)果進(jìn)行投影，這樣還可以減少計(jì)算量。對(duì)X(tn,fn)沿頻率軸方向投影，計(jì)算得到

(6)

式(6)中，N表示離散點(diǎn)數(shù)。對(duì)Y(fn)利用極值法濾波，得到窄帶子信號(hào)，極值法濾波的流程圖如圖4所示。為了能夠更有效地估計(jì)瞬時(shí)頻率，對(duì)濾波得到的子信號(hào)進(jìn)行EMD處理，以使子信號(hào)滿足本征模態(tài)函數(shù)(IMF)的要求。

圖4 極值法濾波的流程圖Fig.4 Flow chart of the extreme filtering method

3 測試分析

對(duì)所提出的多普勒信號(hào)分離算法用Matlab編程進(jìn)行測試分析，測試信號(hào)為圖1所示的實(shí)測信號(hào)。將STFT系數(shù)沿頻率軸方向投影，得到如圖5所示的一維函數(shù)。根據(jù)該圖的一維函數(shù)，利用極值法對(duì)原始信號(hào)進(jìn)行濾波，再對(duì)濾波結(jié)果進(jìn)行EMD處理，得到一系列子信號(hào)，如圖6所示，其中圖6(b)是圖6(a)子信號(hào)的頻譜圖。

為了進(jìn)一步驗(yàn)證所提算法的有效性，用EMD對(duì)原始信號(hào)進(jìn)行分解，分解結(jié)果如圖7所示。對(duì)比圖6和圖7可知，本文算法得到的子信號(hào)帶寬更窄，這樣更有利于估計(jì)信號(hào)的瞬時(shí)頻率。

圖5 實(shí)測多普勒信號(hào)的幅頻曲線Fig.5 The amplitude-frequency curve of the measured Doppler signal

圖6 本文方法得到的子信號(hào)Fig.6 Sub-signals obtained by the proposed method

對(duì)于窄帶的子信號(hào)xi(t)，這里使用下式進(jìn)行瞬時(shí)頻率估計(jì)[10]：

(7)

瞬時(shí)頻率的估計(jì)結(jié)果如圖8所示。對(duì)比圖1的時(shí)頻分布，可知圖8的瞬時(shí)頻率估計(jì)結(jié)果反映了回波信號(hào)中不同多普勒成分隨時(shí)間的變化情況。

圖7 EMD分解后的子信號(hào)Fig.7 Sub-signals decomposed by EMD

利用Hilbert變換來估計(jì)瞬時(shí)頻率也是常用的方法，該方法需要構(gòu)造解析信號(hào)，即[6]

(8)

(9)

對(duì)式(9)進(jìn)行微分，即可得到估計(jì)的瞬時(shí)頻率。經(jīng)過Hilbert變換后估計(jì)的瞬時(shí)頻率如圖9所示， 部分?jǐn)?shù)據(jù)如表1所示，由該圖可知，Hilbert估計(jì)的瞬時(shí)頻率存在很多尖峰。對(duì)比圖8和圖9可知，雖然圖8估計(jì)的瞬時(shí)頻率有些波動(dòng)，但是相對(duì)于圖9來說，圖8的結(jié)果更加合理。

圖8 本文估計(jì)的瞬時(shí)頻率Fig.8 The instantaneous frequency estimated by the proposed method

圖9 Hilbert法估計(jì)的瞬時(shí)頻率Fig.9 The instantaneous frequency estimated by the Hilbert method

表1 部分Hilbert法估計(jì)的瞬時(shí)頻率值

4 結(jié)論

本文提出了基于STFT和EMD的多普勒信號(hào)分離算法。該算法根據(jù)STFT時(shí)頻分布，通過極值濾波法并結(jié)合EMD將回波信號(hào)分解成若干個(gè)窄帶的子信號(hào)，實(shí)現(xiàn)了對(duì)多普勒信號(hào)的分離，進(jìn)而利用窄帶的子信號(hào)來估計(jì)多普勒信號(hào)的瞬時(shí)頻率。實(shí)測驗(yàn)證表明，該算法能夠?qū)㈩l率成分相近的多普勒回波信號(hào)分解成若干個(gè)窄帶的子信號(hào)，有助于實(shí)現(xiàn)對(duì)近場區(qū)多普勒信號(hào)的提取和識(shí)別。