基于經(jīng)驗(yàn)小波變換的目標(biāo)加速度估計(jì)算法

陳浩，郭軍海，齊巍

(北京跟蹤與通信技術(shù)研究所，北京100094)

脈沖雷達(dá)測(cè)速通常采用細(xì)譜線跟蹤技術(shù)，導(dǎo)彈等高動(dòng)態(tài)目標(biāo)的加速度和加加速度會(huì)使回波多普勒譜線展寬甚至出現(xiàn)混疊，導(dǎo)致雷達(dá)測(cè)速系統(tǒng)很難正確跟蹤.因此為了提高脈沖雷達(dá)多普勒測(cè)速精度，估計(jì)目標(biāo)的加速度和加加速度并進(jìn)行相位補(bǔ)償至關(guān)重要［1－2］.當(dāng)目標(biāo)作加速運(yùn)動(dòng)時(shí)，回波信號(hào)為相位具有高階項(xiàng)的非平穩(wěn)信號(hào).目標(biāo)加速度和加加速度分別對(duì)應(yīng)于回波信號(hào)的二階相位系數(shù)和三階相位系數(shù).常用非平穩(wěn)信號(hào)參數(shù)估計(jì)的方法有時(shí)頻分析和基于 Hilbert－Huang變換(HHT)［3］參數(shù)估計(jì)方法.

Wigner－Hough［4］變換、Radon－Wigner［5］變換、Radon－Ambiguity 變換［6－7］、分級(jí)傅里葉變換［8－9］等時(shí)頻分析方法需要進(jìn)行一維搜索或二維搜索，計(jì)算量較大.經(jīng)驗(yàn)?zāi)Ｊ椒纸?EMD)［10－11］是一種適用于非平穩(wěn)信號(hào)處理的新方法.信號(hào)通過(guò)EMD可分解為不同的本征模態(tài)函數(shù)(IMF).IMF是Huang［10］提出的一種調(diào)幅調(diào)頻(AM－FM)信號(hào)，其目的是為了對(duì)IMF作Hilbert變換得到有意義的瞬時(shí)頻率(IF).對(duì)分解得到的各級(jí)IMF作Hilbert變換即可得到對(duì)應(yīng)的瞬時(shí)頻率，再利用直接型主成分提取方法(DPCE)［12］或能量型主成分提取方法(EPCE)［13］就能從得到的各組瞬時(shí)頻率中得到調(diào)頻信號(hào)的頻率主成分.在EMD的基礎(chǔ)上，Wu等人提出了一種叫作總體經(jīng)驗(yàn)?zāi)Ｊ椒纸?EEMD)［14］的噪聲輔助數(shù)據(jù)分析方法.EEMD能解決EMD的模態(tài)混疊(mode mixing)問(wèn)題，同時(shí)對(duì)分解的IMF求總體平均，使分解更具有魯棒性，但由于要進(jìn)行多次EMD分解，也會(huì)有計(jì)算量大的問(wèn)題.結(jié)合EMD的自適應(yīng)性和小波分析的理論框架，Gilles提出了一種稱(chēng)為經(jīng)驗(yàn)小波變換(EWT)［15］的自適應(yīng)信號(hào)處理方法.其核心思想是通過(guò)對(duì)信號(hào)的Fourier譜進(jìn)行自適應(yīng)劃分，建立合適的小波濾波器組來(lái)提取信號(hào)不同的AM－FM成分.對(duì)不同的信號(hào)成分作Hilbert變化可以得到有意義的瞬時(shí)頻率，利用類(lèi)似于EMD的EPCE方法就能提取到調(diào)頻信號(hào)的瞬時(shí)頻率.EWT方法是在小波框架下建立的方法，所以其計(jì)算量遠(yuǎn)小于EMD方法，而且具有較強(qiáng)的魯棒性.

1 經(jīng)驗(yàn)小波變換

經(jīng)驗(yàn)小波本質(zhì)上是根據(jù)信號(hào)頻譜特性選擇的一組帶通濾波器.為了確定帶通濾波器的頻率范圍，對(duì)信號(hào)的Fourier譜進(jìn)行分割.假設(shè)將Fourier支撐［0，π］分割成 N 個(gè)連續(xù)的部分 Λn=［wn－1，wn］，n=1，2，…，N(w0=0，wN= π)，wn選擇為信號(hào)Fourier譜相鄰兩個(gè)極大值點(diǎn)之間的中點(diǎn)，那么，如圖1所示.圖中的過(guò)渡區(qū)域以wn為中心，寬度為2τn.確定分割區(qū)間Λn后，根據(jù)Meyer小波確定的經(jīng)驗(yàn)尺度函數(shù)和小波函數(shù)分別如式(1)和式(2)所示.

圖1 Foureier坐標(biāo)系的分割Fig.1 Partitioning of the Fourier axis

式中，τn= γwn;γ ＜ minn［(ωn－1－ ωn)∕(ωn+1+ωn)］;β(x)=x4(35－84x+70x2－20x3).

原始信號(hào)可被重構(gòu)為

式中:* 為卷積符號(hào);Wεx(0，t)為逼近系數(shù);Wεx(n，t)為x(t)的經(jīng)驗(yàn)小波變換，經(jīng)驗(yàn)?zāi)Ｊ絰k(t)可定義為

2 基于EWT速度與加速度估計(jì)方法

2.1 問(wèn)題分析

設(shè)脈沖雷達(dá)發(fā)射的單載頻脈沖串信號(hào)為

式中，Tr為脈沖持續(xù)時(shí)間;fc為載波頻率;φ0為初始相位.對(duì)應(yīng)的回波信號(hào)為

式中，φd=2π∫fd(t)dt;fd(t)=2v(t)/λ為多普勒頻率，λ為信號(hào)波長(zhǎng).通過(guò)正交解調(diào)得到的輸出信號(hào)為

假設(shè)積累N個(gè)脈沖，每個(gè)脈沖采樣一個(gè)點(diǎn)來(lái)求多普勒速度:

若在t1～tN內(nèi)目標(biāo)作加速運(yùn)動(dòng)，速度近似為v(t)≈v0+at+ja/2t2，則

式中，f0=2v0/λ;k=2a/λ;ka=ja/λ，a 為加速度，ja為加加速度，此時(shí)回波的相位具有三階相位系數(shù).

2.2 頻率主成分提取

令含噪回波信號(hào)的模型為

式中，A為幅度;f0為初始頻率;k為二次相位系數(shù);ka為三次相位稀疏;w(t)為引入FM信號(hào)中的高斯白噪聲.文中仿真信號(hào)設(shè)置A=1，f0=100 Hz，k=400，ka=300，采樣率 fs=2500 Hz.

對(duì)X(t)作EWT變換得到不同的經(jīng)驗(yàn)?zāi)Ｊ絏k(t)，即 X(t)=Xk(t).對(duì) Xk(t)作 Hilbert變換，可得

各經(jīng)驗(yàn)?zāi)Ｊ降慕馕龊瘮?shù)為

式中，ai(t)為瞬時(shí)幅度;θi(t)為瞬時(shí)相位;且

則原信號(hào)可表示為

式中fj(t)為第j個(gè)經(jīng)驗(yàn)?zāi)Ｊ降乃矔r(shí)頻率.得到M組瞬時(shí)頻率fj(t)后，可利用基于能量型主成分提取方法來(lái)提取信號(hào)的頻率主成分.對(duì)于同一時(shí)刻τ，選出各IMF在時(shí)刻τ對(duì)應(yīng)的瞬時(shí)幅度最大值ak對(duì)應(yīng)的經(jīng)驗(yàn)?zāi)Ｊ降乃矔r(shí)頻率作為頻率主成分，即

則信號(hào)的瞬時(shí)頻率為f(t)=fk(t)(t).

2.3 基于EWT的參數(shù)估計(jì)性能分析

在利用EWT變換進(jìn)行FM信號(hào)參數(shù)估計(jì)時(shí)，實(shí)際得到的是一組信號(hào)的瞬時(shí)頻率點(diǎn)f(t):

式中e(t)為頻率噪聲.EWT變換等價(jià)于一組帶通濾波器，EWT變換得到的不同IMF的瞬時(shí)頻率分屬不同的頻率區(qū)間.通過(guò)頻率主成分提取方法可將不同頻率區(qū)間內(nèi)的瞬時(shí)頻率點(diǎn)提取出來(lái)，頻率噪聲僅與信號(hào)的信噪比相關(guān).利用最小二乘方法可以對(duì)參數(shù)α=(f0，k，ka)進(jìn)行估計(jì)，因此參數(shù)的估計(jì)精度僅與回波信號(hào)的信噪比相關(guān).由最小二乘估計(jì)的性質(zhì)可知，估計(jì)量α^為參數(shù)α的線性無(wú)偏估計(jì).根據(jù)式(13)確定待估參數(shù)的 C－R下界［16］為

EWT方法是一種小波分析算法，利用Mallat算法［17］可實(shí)現(xiàn)快速計(jì)算.而傳統(tǒng)的FRFT方法是二維搜索算法，其計(jì)算量遠(yuǎn)大于EWT算法，導(dǎo)致數(shù)據(jù)處理速度較慢.因此，EWT算法比FRFT算法更適用于實(shí)時(shí)數(shù)據(jù)處理.

3 仿真與分析

3.1 理論數(shù)據(jù)仿真

理論仿真分為2部分.第1部分僅考慮信號(hào)具有二次相位項(xiàng).

式中w(t)為加性高斯白噪聲.

圖2(a)為信號(hào)經(jīng)EWT變換后得到的各個(gè)經(jīng)驗(yàn)?zāi)Ｊ剑瑘D2(b)為信噪比等于10 dB時(shí)用EWT方法提取的頻率主成分.從圖2(b)可以看出利用EWT方法提取的頻率主成分受噪聲干擾較小，絕大部分瞬時(shí)頻率點(diǎn)都分布在瞬時(shí)頻率直線兩端.圖3、圖4分別為5dB＜信噪比＜15dB時(shí)，分別用EWT，EEMD－PCA和分?jǐn)?shù)階傅里葉變換(FRFT)方法Monte Carlo仿真50次得到的一次、二次相位系數(shù)估計(jì)均方根誤差(RMSE)圖，并將各參數(shù)估計(jì)誤差與C－R下界作比較.從圖3、圖4中可以看出，基于EWT的參數(shù)估計(jì)方法估計(jì)精度要遠(yuǎn)遠(yuǎn)高于傳統(tǒng)的FRFT方法和EEMD－PCA方法.傳統(tǒng)的FRFT方法在估計(jì)線性調(diào)頻(LFM)信號(hào)初始頻率和高階相位系數(shù)時(shí)，其估計(jì)精度主要受采樣點(diǎn)數(shù)N和旋轉(zhuǎn)角度搜索間隔決定，受信噪比影響不大，因此圖中FRFT方法估計(jì)的參數(shù)誤差隨信噪比變化不大.EWT方法估計(jì)的參數(shù)誤差最為接近C－R下界，且隨著信噪比的增大，EWT逐漸逼近于C－R下界.同等硬件條件下EWT方法運(yùn)行一次的計(jì)算時(shí)長(zhǎng)為0.0079 s，F(xiàn)RFT方法的計(jì)算時(shí)長(zhǎng)為2.711 s，EEMD方法的計(jì)算時(shí)長(zhǎng)為6.328 s.EWT算法由于采用了Mallat小波快速算法，計(jì)算速度要遠(yuǎn)遠(yuǎn)快于FRFT方法和EEMD方法.

圖2 利用經(jīng)驗(yàn)小波變換(EWT)算法得到的經(jīng)驗(yàn)?zāi)Ｊ胶退矔r(shí)頻率Fig.2 Estimated intrinsic mode functions(IMF)and instantaneous frequency(IF)using empirical wavelet transform(EWT)method

圖3 一次相位系數(shù)f0估計(jì)誤差Fig.3 Estimating error of coefficient f0

圖4 二次相位系數(shù)k估計(jì)誤差Fig.4 Estimating error of coefficient k

考慮三階相位項(xiàng)時(shí)，假設(shè)ka=300，則

用本文的方法提取得到的瞬時(shí)頻率如圖5所示，對(duì)獲取的頻率點(diǎn)用抗差最小二乘擬合可估計(jì)各階相位系數(shù).對(duì)具有三階相位項(xiàng)非平穩(wěn)信號(hào)，傳統(tǒng)的FRFT方法不再適用.因此，將基于EWT方法的參數(shù)估計(jì)的估計(jì)誤差與EEMD－PCA方法和各待估參數(shù)的C－R下界作比較.一二階相位系數(shù)估計(jì)誤差與3.1節(jié)類(lèi)似，三次相位項(xiàng)的估計(jì)誤差如圖6所示.信號(hào)的瞬時(shí)頻率為二次曲線，從圖5可以看出，提取的瞬時(shí)頻率集中點(diǎn)分布在二次曲線兩側(cè).從圖6中可以看出，利用本文提出的方法估計(jì)得到的信號(hào)的三次相位系數(shù)估計(jì)精度最接近C－R下界，且隨著信噪比的提高，逐漸逼近于C－R下界.

圖5 具有三階相位系數(shù)信號(hào)的瞬時(shí)頻率Fig.5 Estimated instantaneous frequency(IF)of signal with the third order phase coefficients

圖6 三階相位系數(shù)估計(jì)誤差Fig.6 Estimating error of the third order phase coefficient

3.2 實(shí)測(cè)數(shù)據(jù)仿真

利用某C波段窄帶脈沖雷達(dá)測(cè)量得到的飛行器主動(dòng)段數(shù)據(jù)進(jìn)行仿真.為了滿足目標(biāo)作加速運(yùn)動(dòng)或加加速度運(yùn)動(dòng)的假設(shè)，僅積累50個(gè)回波的I/Q數(shù)據(jù)作為滑動(dòng)窗口，每個(gè)回波采一個(gè)點(diǎn)，窗口每次滑動(dòng)一個(gè)點(diǎn).利用本文的加速度估計(jì)算法得到的加速度作相位補(bǔ)償估計(jì)的速度誤差如圖7所示.利用真實(shí)速度作21點(diǎn)中心平滑得到的加速度均方根誤差為0.03 m/s2，將本文得到的加速度與速度中心平滑得到加速度誤差如圖8所示.從圖中可以看出，利用本文的方法估計(jì)的加速度進(jìn)行相位補(bǔ)償效果很好，使估計(jì)的速度誤差小于0.05 m/s，最大加速度誤差小于0.4 m/s2.

圖7 加速度補(bǔ)償后的速度誤差Fig.7 Velocity error after acceleration compensation

圖8 實(shí)測(cè)數(shù)據(jù)加速度誤差Fig.8 Estimated acceleration error of measured data

4 結(jié)論

本文在經(jīng)驗(yàn)小波變換的基礎(chǔ)上，提出EWT方法對(duì)目標(biāo)徑向加速度進(jìn)行估計(jì)，仿真和理論數(shù)據(jù)對(duì)該算法進(jìn)行驗(yàn)證表明:

1)仿真表明該算法在不同的信噪比條件下均能以較高的精度估計(jì)信號(hào)的參數(shù)，估計(jì)精度高于傳統(tǒng)FRFT算法和EEMD算法，且估計(jì)誤差逼近于C－R下界;

2)計(jì)算速度要遠(yuǎn)遠(yuǎn)快于傳統(tǒng)算法;

3)脈沖雷達(dá)實(shí)測(cè)I/Q數(shù)據(jù)表明，該算法估計(jì)的加速度誤差小于0.4 m/s2，加速度的補(bǔ)償后估計(jì)的速度誤差小于0.05 m/s.

References)

［1］ 袁斌，陳曾平，徐世友，等.基于距離單元篩選快速最小熵的含旋轉(zhuǎn)部件目標(biāo)相位補(bǔ)償方法［J］.電子與信息學(xué)報(bào)，2013，35(5):1128－1134.Yuan B，Chen Z P，Xu S Y，et al.Phase compensation for targets with rotating parts based on range bins selection in fast minimum entropy［J］.Journal of Electronics ＆ Information Technology，2013，35(5):1128－1134(in Chinese).

［2］ 夏猛，楊小牛.基于三次相位補(bǔ)償?shù)倪\(yùn)動(dòng)目標(biāo)參數(shù)估計(jì)［J］.電子科技大學(xué)學(xué)報(bào)，2013，42(4):559－564.Xia M，Yang X N.Parameter estimation for moving target based on three－phase compensation［J］.Journal of University of Electronic Science and Technology of China，2013，42(4):559－564(in Chinese).

［3］ Cao S，Bing P，Lu J，et al.Seismic data time－frequency analysis by the improved Hilbert－Huang transform［J］.Oil Geophysical Prospecting，2013，48(2):246－254.

［4］ Wang Z Z，Liu F，Huang Y，et al.Digitized periodic wignerhough transform and its performance analysis［J］.Telecommunications Engineering，2012，52(9):1452－1458.

［5］ 歐國(guó)建，陳玲瓏，何俞璟.一種多分量LFM信號(hào)參數(shù)估計(jì)的快速仿真算法［J］.重慶郵電大學(xué)學(xué)報(bào):自然科學(xué)版，2013，25(4):459－463.Ou G J，Chen L L，He Y J.A simulation fast algorithm for parameters estimationof the multicomponent LFM signals［J］.Journal of Chongqing University of Posts and Telecommunications，2013，25(4):459－463(in Chinese).

［6］ Yu Y.Detection and parameter estimation of linear frequency modulated signals based on radon transform［J］.Journal of Modern Defence Technology，2013，41(1):136－141.

［7］ Zhu Y W，Zhao Y J，Jia W G.Fast parameter estimation method for LFM signal based on ambiguity function slice and FrFT［J］.Journal of Information Engineering University，2012，13(2):218－223.

［8］ Mohindru P，Khanna R K R，Bhatia S S.Analysis of chirp signal with fractional Fourier transform［J］.Majlesi Journal of Multimedia Processing，2013，2(1):314－322.

［9］ Din? E，Duarte F B，Machado J A T，et al.Application of continuous wavelet transform to the analysis of the modulus of the fractional Fourier transform bands for resolving two component mixture［J］.Signal，Image and Video Processing，2013，21(10):1－7.

［10］ Huang N E.The empirical modedecomposition and the Hilbert spectrum for nonlinear andnon－stationary time series analysis［J］.Proceedings of Royal Society of London，1998，A4(54):903－995.

［11］ Flandrin P，Rilling G，Goncalves P.Empirical mode decomposition as a filter bank［J］.Signal Processing Letters，IEEE，2004，11(2):112－114.

［12］ 崔華.一種新的線性調(diào)頻信號(hào)的瞬時(shí)頻率估計(jì)方法［J］.計(jì)算機(jī)應(yīng)用研究，2008，25(8):2532－2533.Cui H.New method for instantaneous frequency estimations of LFM signals［J］.Application Research of Computers，2008，25(8):2532－2533(in Chinese).

［13］ 王燕，鄒男，付進(jìn)，等.基于局部瞬時(shí)能量密度級(jí)的瞬態(tài)信號(hào)檢測(cè)方法［J］.電子與信息學(xué)報(bào)，2013，35(7):1720－1724.Wang Y，Zou N，F(xiàn)u J，et al.Transient signal detection method based on partial instantaneous energy density level［J］.Journal of Electronics ＆ Information Technology，2013，35(7):1720－1724(in Chinese).

［14］ Wu Z H，Huang N E.Ensemble empirical mode decomposition:a noise assisted data analysis method［J］.Advances in Adaptive Data Analysis，2009，1(1):1－41.

［15］ Gilles J.Empirical wavelet transform［J］.IEEE Transactions in Signal Processing，2013，61(16):3999－4010.

［16］ Steven M K.Fundamentals of statistical signal processing，Volume I:estimation theory［M］.Beijing:Publishing House of E－lectronics Industry，2006:25－39.

［17］ 成禮智.小波的理論與應(yīng)用［M］.北京:科學(xué)出版社，2009:75－88.Cheng L Z.Theories and applications of wavelet［M］.Beijing:Science Press，2009:75－88(in Chinese).