魯逸杰，宮志華，張 群，王劍欽，李開明

(1.空軍工程大學(xué)信息與導(dǎo)航學(xué)院，陜西 西安 710077;2.中國人民解放軍63850部隊，吉林 白城 137001;3.中國人民解放軍95333部隊，湖南 長沙 410114)

0 引言

彈道目標(biāo)在飛行時，為了保持姿態(tài)穩(wěn)定，需繞自身對稱軸作自旋運動，同時由于彈箭分離以及誘餌釋放時形成的橫向作用力，目標(biāo)還將繞空間某定向軸進行進動。美國海軍研究實驗室(United States Naval Research Laboratory, NRL)的Victor C. Chen最早將進動定義為微動的一種[1]。目標(biāo)的微動會對雷達回波進行頻率、相位等調(diào)制，產(chǎn)生相應(yīng)的微多普勒效應(yīng)。微多普勒效應(yīng)反映了目標(biāo)微動特點，是目標(biāo)區(qū)分的有效獨特效應(yīng)。近年來，基于微多普勒特性進行彈道目標(biāo)分類識別已成為一個重要研究方向，從雷達回波中準(zhǔn)確提取目標(biāo)微多普勒頻率及其微動特征成為目標(biāo)分類識別的關(guān)鍵。

目前，已有文獻對等效散射中心模型下的微多普勒特性進行了研究。文獻[2]基于雙視角下的距離像序列推導(dǎo)出目標(biāo)參數(shù)的解析解。文獻[3]利用多散射中心徑向距離歷程的關(guān)聯(lián)結(jié)果獲得了目標(biāo)的三維運動和三維成像結(jié)果，基于三維運動的歐氏重構(gòu)采用序列二次規(guī)劃優(yōu)化和非線性最小二乘循環(huán)迭代的方法估計目標(biāo)的進動參數(shù)。文獻[4]基于中段目標(biāo)的電磁散射特性，建立了此類目標(biāo)在不同姿態(tài)下的ISAR成像結(jié)果與目標(biāo)二維幾何特征的映射關(guān)系，提出一種基于ISAR像的中段目標(biāo)二維幾何特征反演方法。文獻[5]通過分析散射點一維距離像隨時間變化的規(guī)律，提出了一種基于時間-距離像的彈道目標(biāo)進動特征提取方法。文獻[6]提出了一種基于快速傅里葉變換和自適應(yīng)濾波并結(jié)合 Hilbert變換來估計雷達多普勒信號頻率的方法。文獻[7]基于時頻分布，對低信噪比下雷達回波信號瞬時頻率進行了估計。

需要指出的是，寬帶信號在獲得目標(biāo)微動特征方面的研究較多，而窄帶雷達回波多普勒譜中包含的散射點頻率變化同樣可以反映目標(biāo)的微動特征。在低信噪比條件下，微動特征的有效提取估計目前仍是一個有待解決的問題。本文針對此問題，提出了基于變分模態(tài)分解的進動目標(biāo)微多普勒特征提取方法。

1 變分模態(tài)分解

變分模態(tài)分解(VMD)算法是2014年由Dragonmiretskiy等人提出的一種自適應(yīng)的信號處理方法[2]。該方法在獲取分解分量的過程中通過迭代搜尋變分模型最優(yōu)解來確定每個分量的頻率和中心帶寬，從而能夠自適應(yīng)地實現(xiàn)信號的頻域分割以及各分量的有效分離，VMD的分解過程主要分為以下三步：

步驟一 構(gòu)建目標(biāo)函數(shù)。它分為三步進行：1)對每個模態(tài)通過希爾伯特變換計算其解析信號來獲得單邊譜；2)將各個模態(tài)的解析信號與各自指數(shù)校準(zhǔn)的中心頻率混頻，將頻譜搬移至基帶；3)計算解調(diào)信號的范數(shù)，并估計各模態(tài)信號的帶寬，構(gòu)建出目標(biāo)函數(shù)如下：

(1)

式(1)中，{uk}:={u1,…,uk}為所有模態(tài)集合，{wk}:={w1,…,wk} 為各模態(tài)對應(yīng)的中心頻率。δ(t)為沖激函數(shù)，?t表示對t求偏導(dǎo)數(shù)。

步驟二 引出二次懲罰因子α及拉格朗日乘子λ(t)，把約束性的變分問題轉(zhuǎn)化成非約束性的變分問題，其中的二次懲罰因子在高斯噪聲存在情況下可保證信號重構(gòu)精度，拉格朗日乘子則使得約束條件保持嚴(yán)格性，擴展拉格朗日的表達式如下：

(2)

步驟三 利用交替方向乘子法來解決上述變分問題，求解上式的迭代過程如下：

(3)

再將上式變換至頻域，得：

(4)

該二次優(yōu)化問題的解為：

(5)

(6)

將上式同樣的轉(zhuǎn)化至頻域可得：

(7)

該二次優(yōu)化問題的解為：

(8)

(9)

(10)

(11)

3) 根據(jù)式(19)更新λ;

(12)

4) 重復(fù)2)和步驟3),直至滿足迭代中止條件：

(13)

式(13)中，ε為判別精度，ε>0。

5) 輸出結(jié)果，得到K個模態(tài)。

通過以上算法，可以獲得K個模態(tài)與其對應(yīng)的時域波形，即可得到不同頻率的信號，由VMD分解得到的信號分量IMF(Intrinsic Mode Function)一般都具有明確的物理意義。在理想的條件下，我們可以通過對時頻圖形上分離出的微多普勒(micro-Dopple, m-D)曲線進行VMD分解并從中得到一些微動特征參數(shù)。

2 基于VMD的進動目標(biāo)微多普勒頻率提取方法

進動錐體目標(biāo)的等效散射中心主要包括錐頂散射點和錐底散射點，如圖1。

圖1 雷達接收機與低空目標(biāo)的幾何關(guān)系示意圖Fig.1 The geometry of radar and target

假設(shè)目標(biāo)以角速度ωs繞自身對稱軸自旋的同時，以角速ωc繞空間某定向軸作錐旋運動。以雷達位置O為坐標(biāo)原點構(gòu)建雷達坐標(biāo)系(X,Y,Z)，以目標(biāo)質(zhì)心O′為坐標(biāo)原點構(gòu)建與雷達坐標(biāo)系平行，且只隨目標(biāo)的平動而運動的參考坐標(biāo)系(U,V,W)，以目標(biāo)質(zhì)心O′為原點構(gòu)建隨目標(biāo)的運動而運動的目標(biāo)本地坐標(biāo)系(x,y,z)，目標(biāo)平動速度為v。由于散射中心A位于錐頂，其自旋線速度為0，因此錐體的進動導(dǎo)致A所對應(yīng)的微多普勒頻率隨時間呈現(xiàn)正弦變化[9]。散射中心B除繞進動軸進動外還繞自旋軸進行自旋，所表現(xiàn)的運動為自旋和錐旋的合成運動。

從內(nèi)容與寫作技法言，全詞詠茶，而不著一個茶字，卻連用若干典故，老練精熟，此筆法在黃庭堅詩中十分常見，將此筆法移之入詞，屬順理成章之事。

(14)

式(14)中，t為時間，σ為散射系數(shù)，fc為載頻，c為光速。對Φ(t)關(guān)于時間t求導(dǎo)，得到對應(yīng)的多普勒頻率：

(15)

(16)

(17)

式(17)中，

(18)

由式(17)可以看出，進動產(chǎn)生的微多普勒頻率變化規(guī)律表現(xiàn)為多正弦分量疊加，其主要頻率分量為自旋頻率Ωs，錐旋頻率Ωc，及其代數(shù)和差|Ωs+Ωc|，|Ωs-Ωc|。

隨后，根據(jù)上述進動目標(biāo)散射中心的回波模型，采用時頻分析的方法對回波進行處理。Gabor 變換是一個時間-頻率分析方法，它采用高斯時間窗口在時域的滑動截取該時段的頻譜信息，簡單說就是窗口可移動，但是形狀不改變。Gabor 變換是高斯窗口短時Fourier 變換，Gabor 變換可以達到時-頻局部化的目的，它能夠在整體上提供信號的全局頻譜信息又能提供局部時刻的頻譜變化信息。

通常多散射中心目標(biāo)的時頻曲線會出線交叉，本文首先需要將散射中心A，B的微動曲線的分離開來，再對散射中心B的微動曲線進行分解。文獻[10]介紹了一種骨架提取和曲線分離算法，通過對曲線進行行高斯平滑濾波、模值二值化等處理來實現(xiàn)。本文采用此算法對時頻圖像上的微動曲線進行曲線分離和骨架提取，得到了兩個散射中心微動曲線。

高斯平滑濾波具體處理過程是：構(gòu)建以(x,y)為期望，σ為標(biāo)準(zhǔn)差的二維高斯函數(shù)：

(19)

對二維高斯函數(shù)進行離散化，得到該二維高斯函數(shù)離散射后的函數(shù)值。隨后將各個位置的坐標(biāo)帶入到離散射化后高斯函數(shù)中，得到高斯平滑模板的系數(shù)值。設(shè)窗口模板的大小為：(2k+1)×(2k+1)，模板中各個元素值的計算公式如下：

(20)

式(20)中，i,j分別為時頻平面坐標(biāo)值，經(jīng)過對時頻平面數(shù)據(jù)的二維卷積，實現(xiàn)高斯平滑濾波。

時頻平面數(shù)據(jù)的二值化是對時頻平面復(fù)數(shù)據(jù)取模，對于給定的閾值T，通過調(diào)整合適的閾值T，用T將圖像的數(shù)據(jù)分成兩部分：大于T的像素群和小于T的像素群。將大于T的像素群的像素值設(shè)定為白色(或者黑色)，小于T的像素群的像素值設(shè)定為黑色(或者白色)。最終得到了兩個散射中心微動曲線。

最后，運用VMD算法對兩個散射中心的微動曲線進行分解，對于錐底散射點B，在假設(shè)平動已補償?shù)那疤嵯?，其m-D曲線為多個正弦分量疊加[9]，主要頻率成分為自旋頻率Ωs，錐旋頻率Ωc，及其代數(shù)和差|Ωs+Ωc|，|Ωs-Ωc|。對散射點B作VMD分解得到的IMF1為頻率成分|Ωs+Ωc|的分量，IMF2分量即為角頻率成分為Ωc與Ωs中較大的分量，IMF3分量為角頻率Ωc與Ωs中較小的分量，IMF4為成分為角頻率|Ωs-Ωc|的分量。與目前已有的用EMD算法來分解m-D曲線的做法相比，VMD算法可以克服EMD算法存在的模態(tài)混疊問題，并且具有更良好的抗噪聲性能。

3 仿真對比與分析

為了驗證VMD方法的有效性和抗噪性，分別進動目標(biāo)進行無噪和有噪兩種情況進行仿真對比與分析。

1)無噪聲情況下進動目標(biāo)微動特征分析能力仿真。

設(shè)目標(biāo)地面半徑為0.5 m，高為0.5 m，雷達發(fā)射信號載頻fc=10 GHz，目標(biāo)與雷達初始距離為(300 km,100 km,500 km)，本地坐標(biāo)系與參考坐標(biāo)系的初始Euler角為(π/3 rad,π/4 rad,π/5 rad)。目標(biāo)平動速度為0，自旋角速度ωs=(19.233 2,-11.107 2,22.214 4)Trad/s，Ωs=10π rad/s,自旋周期Ts=0.2s；錐旋角速度ωc=(8.8813,-10.8773,12.5601)Trad/s，Ωs=6π rad/s。錐旋周期Tc=0.33s，目標(biāo)上有兩個強散射點，即錐頂散射點A(0 m,0 m,0.5 m)和錐底散射點B(0.5 m,0 m,-0.5 m)。脈沖重復(fù)率為3 000 Hz，回波信號時長為2 s，圖2給出了時頻平面上的m-D曲線。對m-D曲線進行高斯平滑濾波和二值化處理的結(jié)果如圖3所示。對二值化處理的結(jié)果進行骨架提取的結(jié)果如圖4所示。對圖4中的兩條m-D曲線進行曲線分離，可以準(zhǔn)確提取到2條m-D曲線。

圖2 時頻分析結(jié)果Fig.2 Time frequency analysis result

圖3 二值化處理結(jié)果Fig.3 Binarization processing result

圖4 骨架提取結(jié)果Fig.4 Skeleton extraction result

分離出的錐底散射點B的m-D曲線，對其分別做VMD分解和EMD分解獲得的IMF如圖5所示。

從圖5可以看出，VMD算法可以清晰地分離出m-D曲線所包含的四種分量。IMF1為8 Hz，IMF2為5 Hz，IMF3為3 Hz，IMF4為2 Hz。而EMD算法的分解結(jié)果明顯得出現(xiàn)了模態(tài)混疊的現(xiàn)象。雖然IMF1仍能得出正確的高頻分量，但是由于頻率分量比較接近，使得EMD無法準(zhǔn)確分辨不同的低頻成分。而VMD算法作為一種比較新的非線性信號處理算法，利用變分原理使各個模態(tài)不斷調(diào)整中心頻率，比EMD算法具有更好的自適應(yīng)性，解決了EMD方法的缺陷，分解效果更準(zhǔn)確、有效。

圖5 VMD與EMD分解結(jié)果Fig.5 VMD and EMD decomposition results

2)有噪聲情況下進動目標(biāo)微動特征分析能力仿真。

在噪聲條件下，進一步驗證VMD算法和EMD算法的表現(xiàn)。圖6給出了在信噪比為-10 dB的情況下。VMD算法和EMD算法的分解結(jié)果?？梢钥闯?，VMD算法仍然可以提取出除了|Ωs-Ωc|以外的其他分量，而EMD算法不僅僅出現(xiàn)了嚴(yán)重的模態(tài)混疊，還出現(xiàn)了無意義的虛假分量。這是因為EMD方法具有二進制濾波特性，對于復(fù)雜信號噪聲和脈沖的干擾會影響包絡(luò)線計算，不斷傳播擴大的誤差會產(chǎn)生模態(tài)混疊從而導(dǎo)致出現(xiàn)錯誤的IMF分量，失去了其具體物理意義[11]。而VMD算法本質(zhì)是一組自適應(yīng)維納濾波器，具有較好的噪聲魯棒性。因此在噪聲條件下，其分解結(jié)果要比EMD算法更準(zhǔn)確有效。

圖6 -10 dB噪聲下VMD與EMD分解結(jié)果Fig.6 VMD and EMD decomposition results under -10 dB noise

4 結(jié)論

本文提出了基于VMD的微多普勒頻率提取方法。該方法先對Gabor 變換后時頻分布結(jié)果的進行骨架提取，通過高斯平滑濾波法并結(jié)合VMD將回波信號分解成若干獨立曲線，并由分解得到的曲線結(jié)合其微多普勒曲線特征來估計目標(biāo)的自旋頻率和錐旋頻率。仿真結(jié)果表明，該方法能夠在低信噪比條件下將目標(biāo)的微動曲線分解成若干獨立曲線，并獲得其微動頻率。