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

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

0 引言

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

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

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

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

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

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

(1)

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

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

(2)

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

(3)

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

(4)

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

(5)

(6)

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

(7)

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

(8)

(9)

(10)

(11)

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

(12)

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

(13)

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

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

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

2 基于VMD的進動目標微多普勒頻率提取方法

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

圖1 雷達接收機與低空目標的幾何關系示意圖Fig.1 The geometry of radar and target

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

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

(14)

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

(15)

(16)

(17)

式(17)中，

(18)

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

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

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

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

(19)

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

(20)

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

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

最后，運用VMD算法對兩個散射中心的微動曲線進行分解，對于錐底散射點B，在假設平動已補償?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方法的有效性和抗噪性，分別進動目標進行無噪和有噪兩種情況進行仿真對比與分析。

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

設目標地面半徑為0.5 m，高為0.5 m，雷達發(fā)射信號載頻fc=10 GHz，目標與雷達初始距離為(300 km,100 km,500 km)，本地坐標系與參考坐標系的初始Euler角為(π/3 rad,π/4 rad,π/5 rad)。目標平動速度為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，目標上有兩個強散射點，即錐頂散射點A(0 m,0 m,0.5 m)和錐底散射點B(0.5 m,0 m,-0.5 m)。脈沖重復率為3 000 Hz，回波信號時長為2 s，圖2給出了時頻平面上的m-D曲線。對m-D曲線進行高斯平滑濾波和二值化處理的結果如圖3所示。對二值化處理的結果進行骨架提取的結果如圖4所示。對圖4中的兩條m-D曲線進行曲線分離，可以準確提取到2條m-D曲線。

圖2 時頻分析結果Fig.2 Time frequency analysis result

圖3 二值化處理結果Fig.3 Binarization processing result

圖4 骨架提取結果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算法的分解結果明顯得出現(xiàn)了模態(tài)混疊的現(xiàn)象。雖然IMF1仍能得出正確的高頻分量，但是由于頻率分量比較接近，使得EMD無法準確分辨不同的低頻成分。而VMD算法作為一種比較新的非線性信號處理算法，利用變分原理使各個模態(tài)不斷調(diào)整中心頻率，比EMD算法具有更好的自適應性，解決了EMD方法的缺陷，分解效果更準確、有效。

圖5 VMD與EMD分解結果Fig.5 VMD and EMD decomposition results

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

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

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

4 結論

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