彭正紅, 楊德貴,*, 王 行, 王 浩, 朱政亮

(1. 中南大學(xué)航空航天學(xué)院, 湖南 長沙 410083; 2. 中南大學(xué)自動化學(xué)院, 湖南 長沙 410083;3. 浙江大學(xué)工業(yè)控制技術(shù)國家重點實驗室, 浙江 杭州 310027;4. 廈門大學(xué)水聲通信與海洋信息技術(shù)教育部重點實驗室, 福建 廈門 361005)

0 引 言

現(xiàn)代戰(zhàn)場環(huán)境日益復(fù)雜,誘餌與彈頭的結(jié)構(gòu)與特性越來越相似,如何有效識別真假目標成為研究現(xiàn)代電子戰(zhàn)的重中之重。隨著超寬帶雷達技術(shù)、雷達信號處理技術(shù)等的快速發(fā)展,現(xiàn)代雷達能夠提供更加精確的目標信息,其中就包括因目標微動引起的微多普勒(micro-Doppler,m-D)特征,這些特征的提取問題也就變得十分重要。通常,m-D特征提取研究工作分為微動信號建模、m-D曲線分離及微動特征提取3個方面。

雷達目標微動及m-D的概念是由Chen首先提出的,并針對這一問題開展了大量研究[1-2],為后來者的研究提供了堅實基礎(chǔ)。文獻[3]針對進動目標,推導(dǎo)了其m-D頻率理論公式,并對錐體彈頭m-D特性進行了仿真分析。文獻[4]就有翼彈頭和無翼彈頭進動對雷達回波的影響,進行了理論分析和微波暗室實驗。文獻[5]針對目標自旋、進動及章動3種微動形式下瞬時頻率的差異提取出4種特征,采用支持向量機(support vertor machine,SVM)分類器對3種微動形式進行了有效分類。

對于微動特征提取的問題,目前應(yīng)用較為廣泛的是模態(tài)分解方法[6-9]。利用經(jīng)驗?zāi)B(tài)分解(empirical mode decomposition, EMD)[10]、變分模態(tài)分解(variational mode decomposition,VMD)、改進型模態(tài)分解方法[11-13]等分解出m-D曲線中的固定模態(tài),經(jīng)過相應(yīng)計算獲得目標的微動特征。然而,由于目標上不同區(qū)域散射點微動形式存在差異,產(chǎn)生的m-D曲線會出現(xiàn)復(fù)雜的交疊現(xiàn)象,如果不能對m-D曲線進行精確分離,就會給后續(xù)微動特征和參數(shù)估計造成困難。

為了解決m-D曲線的交疊問題,學(xué)者們提出了許多方法。目前較為常見的方法是根據(jù)能量峰值進行提取分離[14],但是對曲線交疊處的分離情況較差。也有利用Viterbi算法[15]、最近鄰原則[16]逐條提取m-D曲線的方法,但同樣存在曲線交疊處分離效果差的情況。為解決交疊處m-D曲線分離效果差這一問題,文獻[17-18]提出了一種曲線路徑重組加固有線性調(diào)頻成分分解方法,分離出m-D曲線,但該方法需要提前知曉m-D曲線的數(shù)量,并且針對的是一維信號數(shù)據(jù),具有一定的局限性。

針對目標m-D曲線交疊處分離效果差這一問題,本文提出一種基于曲線趨勢估計的m-D曲線分離和特征提取算法。首先建立目標等效散射點章動模型以及回波模型,獲取目標m-D曲線;然后通過形態(tài)學(xué)處理獲得穩(wěn)定精細的二值化曲線,利用曲線光滑性及插值理論對曲線趨勢進行精確估計,并對其進行有效分離;最后在獲得各散射點的m-D曲線的基礎(chǔ)上,利用VMD算法分解出各個曲線的所有模態(tài)并與EMD分解結(jié)果進行對比,估算出各個散射點的各種微動頻率。

1 復(fù)雜運動目標微動回波模型

1.1 復(fù)雜運動目標散射點微動模型構(gòu)建

雷達目標的m-D效應(yīng)是由目標的微動引起的,通常目標微動可分解為自旋、錐旋以及章動。為簡化分析,假設(shè)彈頭的形狀為圓錐形,誘餌的形狀為圓錐形和桿狀,建立坐標系如圖1所示。圖1中,O-XYZ是雷達坐標系;O′-XYZ是參考坐標系,原點O′位于雷達坐標系方位角為ν、俯仰角為α、距離為R0的位置,與雷達坐標系平行;O′-xyz是目標本體坐標系,z軸為錐體目標對稱軸,也是其自旋軸l;u軸是目標的進動軸,進動角為θ;w軸是章動軸,是l與u的叉乘。

圖1 目標微動模型Fig.1 Target micro-motion model

假設(shè)l軸在本體坐標系中的單位矢量為l=(lx,ly,lz)T,u軸在參考坐標系中的單位矢量為u=(uX,uY,uZ),雷達視線在參考坐標系中的矢量表示為r=(rX,rY,rZ),目標上一點在本體坐標系中的坐標為p0=(x0,y0,z0),目標自旋角速度為ω1,進動角速度為ω2,章動角速度為ω3。根據(jù)文獻[19]中所提方法,可以得出目標在運動過程中,各點在參考坐標系中的實時坐標為

(1)

(2)

則目標上各點到雷達的實時距離為

R(tm)=R0+P(tm)n

(3)

式中:n=[-cosαcosν,-cosαsinν,-sinα]

1.2 雷達回波模型構(gòu)建及回波生成

采用大時寬帶寬積的線性調(diào)頻信號:

(4)

(5)

對回波信號作解線頻調(diào)處理,設(shè)參考信號為

(6)

(7)

按照文獻[20]中的方法補償?shù)羰Ｓ嘁曨l相位項和斜置項,得到目標混頻信號為

(8)

對快時間進行傅里葉變換,得

exp[j2πfc(τ0-τi)]

(9)

2 基于趨勢估計的m-D曲線分離

對于包含多種分量的微動信號而言,其m-D曲線都有可能存在交疊情況。假設(shè)某一微動信號x1(t)包含一個零頻分量和兩個相位相反的正弦分量,計算信號的m-D,結(jié)果如圖2所示。觀察圖2可以發(fā)現(xiàn),3種分量的m-D曲線在某一點(或區(qū)間)處存在交疊情況。

圖2 信號x1(t)的m-D圖Fig.2 m-D diagram of signal x1(t)

現(xiàn)有m-D曲線分離方法諸如最近鄰原則法、Viterbi算法、峰值法等,在進行分離時,均會產(chǎn)生非完全分離、分離效果差的現(xiàn)象,如圖3所示。為解決這一問題,本文提出了一種基于曲線趨勢估計的分離算法。首先通過骨架提取獲得穩(wěn)定的二值化曲線,再利用曲線光滑性和插值理論對曲線進行分離,具體算法如下所示。

圖3 傳統(tǒng)算法分離效果Fig.3 Separation effects of traditional algorithms

2.1 骨架提取

首先根據(jù)目標回波信號獲取其m-D曲線,然后運用形態(tài)學(xué)圖像處理中的骨架提取算法對m-D曲線進行“變細”操作。骨架提取算法的基本思想就是“腐蝕”運算和“開”運算,可表示為

(10)

式中:Θ表示“腐蝕”運算;°表示“開”運算;|Sif(f,tm)|ΘkB表示連續(xù)k次對|Sif(f,tm)|進行 “腐蝕”運算;K=max{k|(|Sif(f,tm)|ΘkB)≠?}表示使|Sif(f,tm)|不是空集的最大“腐蝕”次數(shù)。

為防止m-D曲線毛刺對骨架提取算法的干擾,先將獲得的m-D曲線進行二值化和高斯平滑處理,再進行骨架提取。

2.2 曲線趨勢估計及分離

骨架提取得到最“細”的二值化m-D曲線后,利用“清除”思想將各散射點的m-D曲線分離出來。設(shè)獲得的m-D曲線數(shù)據(jù)矩陣為NNr×Na,逐列提取NNr×Na中數(shù)據(jù)。理想情況下,估計后續(xù)點的值可通過判斷前一段區(qū)間內(nèi)的導(dǎo)數(shù)值(即曲線光滑性)來進行,區(qū)間長度i0可根據(jù)實際分離曲線的周期進行修整。然而,由于交疊區(qū)間的存在及骨架提取得到的m-D曲線交疊區(qū)間并非單位長度,導(dǎo)致上述方法無法進行有效分離,因此本文將曲線的交疊區(qū)域和非交疊區(qū)域分開估計,具體步驟如下。

步驟 1計算需分離的曲線數(shù)量

(11)

式中:p(Numi)是各列極值點數(shù)量出現(xiàn)的概率。

步驟 2數(shù)據(jù)優(yōu)化處理,初步確定交疊區(qū)間

當numi>curve_number,剔除多余極值點;

當numi

當numi=curve_number,判斷相鄰兩條曲線之間的距離,若小于閾值ΔmD(整數(shù)),視為交叉區(qū)間,N(:,i)=0,否則不作處理,即:

N(rm(m),i)-N(rm(n),i)<ΔmD,m=1,2,…,numi-1;
n=m+1

(12)

上述確定的交叉區(qū)間的左右端點分別存儲在行向量qd和hd中。

步驟 3截取數(shù)據(jù),優(yōu)化交疊區(qū)間。截取中間93%以上的數(shù)據(jù),截取后第一列記為ii;合并相鄰較近及單位長度的交疊區(qū)間,即當qd(i)-hd(i-1)≤3或hd(i)-qd(i)=1時,刪除qd(i)、hd(i)或qd(i)、hd(i-1)。交疊區(qū)間的優(yōu)化,可有效減少因骨架提取缺陷所產(chǎn)生的偽交疊區(qū)間以及合并被分割的交疊區(qū)間,能夠提高算法運行速度,減少計算量。

步驟 4曲線估計。依次提取每列數(shù)據(jù),當在非交疊區(qū)間內(nèi),轉(zhuǎn)至步驟5;當在交疊區(qū)間內(nèi),轉(zhuǎn)至步驟6;當分離最后一條曲線時,轉(zhuǎn)至步驟7。

步驟 5非交疊區(qū)間曲線估計

(1) 當i-ii+1=2時

(13)

(2) 當3≤i-ii+1

(14)

(3) 當i0≤i-ii+1≤Na時

(15)

步驟 6交疊區(qū)間曲線估計

(16)

式中:{A|B}表示在滿足B的條件下進行A操作。

(17)

交疊區(qū)間內(nèi)的值為

(18)

步驟 7最后一條曲線估計。若qd(jj)

步驟 8清除當前曲線,返回步驟4。

具體流程圖如圖4所示。

圖4 曲線估計與分離過程Fig.4 Curve estimation and separation process

相比于傳統(tǒng)m-D曲線分離方法,本文所提方法利用骨架提取方法及交疊區(qū)間附近的曲線趨勢,避免了諸如峰值法在交疊處的非首條曲線能量急劇減少導(dǎo)致分離錯誤的問題、最近鄰原則法在交疊嚴重處極易錯誤估計的問題以及Viterbi算法在交疊處不具備有效區(qū)分度的問題等。

為更加直觀地說明算法的有效性,對信號x1(t)及信號x2(t)進行曲線分離。其中,信號x1(t)組成分量如第2節(jié)開頭所述,信號x2(t)由一個線性分量加兩種不同頻率的正弦分量組成,分離結(jié)果如圖5所示?？梢园l(fā)現(xiàn),本文算法對兩種信號的分離效果均很好,有效解決了曲線交疊問題。

圖5 兩種示例信號的分離效果Fig.5 Separation effect of two sample signals

3 基于VMD算法的微動分量提取

非平穩(wěn)信號由多個單分量信號組成,研究非平穩(wěn)信號的組成單分量十分有必要。VMD不同于傳統(tǒng)模態(tài)分解方法的遞歸模式,其將信號的分解過程轉(zhuǎn)化到了變分模式中,主要過程是將一個非平穩(wěn)信號分解為指定個數(shù)的相互獨立的模態(tài)uk中,每個模態(tài)對應(yīng)一個中心頻率wk。除此之外,為使分解出的各分量的帶寬之和最小,需要通過希爾伯特變換計算出每個模態(tài)的單邊頻譜,再根據(jù)各自的中心頻率得到各模態(tài)的基帶信號,最后根據(jù)高斯平滑統(tǒng)計各個模態(tài)帶寬之和:

(19)

式中:{uk}={u1,u2,…,uk}為所有模態(tài)的集合;{wk}={w1,w2,…,wk}為所有模態(tài)中心頻率的集合;?t表示對t求偏導(dǎo)。

通過引入二次懲罰項和拉格朗日乘子將問題轉(zhuǎn)化為無約束問題,式(19)將轉(zhuǎn)化為增廣拉格朗日形式:

(20)

式中:N是期望分解出的模態(tài)數(shù)量;β是懲罰因子;λ是拉格朗日乘子。在頻域中使用交替方向乘子法,得出u,w,λ在頻域的更新公式:

(21)

(22)

(23)

合理設(shè)置參數(shù)N,β,λ,最大程度上對非平穩(wěn)信號進行分解,得出本征模態(tài)函數(shù)(intrinsic mode function,IMF)。

4 仿真分析

設(shè)共有4個目標,各自屬性如表1所示,當目標為錐形結(jié)構(gòu)時,以錐鼻以及一對對稱的錐尾邊緣散射點構(gòu)建目標散射點模型;當目標為桿狀結(jié)構(gòu)時,以兩端散射點構(gòu)建目標模型。

表1 目標基本屬性

不考慮目標的平動,雷達發(fā)射信號為脈沖體制的線性調(diào)頻信號,脈沖寬度為10 μs,脈沖重復(fù)頻率為200 Hz,載頻為9 GHz,帶寬為2 GHz。各目標散射點散射強度為[0,1]中的隨機值。

通過第2節(jié)和第3節(jié)所述方法,4個目標對應(yīng)的m-D曲線均能實現(xiàn)較好的分離,且能正確提取出目標的微動特征,能有效驗證本文所提算法。由于篇幅的限制,現(xiàn)選擇目標3和目標4進行結(jié)果展示及詳細分析。圖6是目標3在f-tm平面上的m -D曲線。對圖6中曲線進行二值化,再利用1×1像素的高斯空間掩膜對其進行高斯平滑處理,通過骨架提取算法得到圖7所示圖像。設(shè)置曲線分離參數(shù)i0=15,ε=9,ΔmD=4,得到分離結(jié)果,并與峰值法、最近鄰原則法及Viterbi算法進行比較,結(jié)果如圖8所示。

圖6 原始m-D曲線(目標3)Fig.6 Original m-D curve (target 3)

圖7 二值化及骨架提取得到的m-D曲線(目標3)Fig.7 m-D curve obtained by binarization and skeleton extraction (target 3)

圖8 傳統(tǒng)算法與本文算法的分離效果(目標3)Fig.8 Separation effect by traditional algorithms and proposed algorithm (target 3)

觀察圖8可以發(fā)現(xiàn),峰值法對能量依賴較大,當待分離曲線的能量差別較小時,會出現(xiàn)圖8(a)所示的結(jié)果;最近鄰原則法在曲線交疊處基本無效;Viterbi算法在兩交疊點之間曲線路徑相差較小的情況下易出現(xiàn)錯誤分離現(xiàn)象;而本文算法在分離較為復(fù)雜的交疊曲線時能夠取得較好的效果。

利用本文算法分離出的m-D曲線,分別采用EMD和VMD算法,合理設(shè)置參數(shù),得到該目標m-D曲線的分解結(jié)果如圖9所示。

圖9 目標3 m-D曲線EMD及VMD分解結(jié)果Fig.9 M-D curve EMD and VMD decomposition results of target 3

目標3為桿狀目標,兩端點處散射點的m-D特性基本上一致。桿狀結(jié)構(gòu)目標自旋軸是自身,故微動只需考慮進動和章動,根據(jù)表1及理論算法,應(yīng)當分解出1.0 Hz,1.6 Hz,2.6 Hz,3.6 Hz 4種頻率分量。觀察圖9發(fā)現(xiàn),采用EMD分解,混疊現(xiàn)象嚴重,很難觀察到各單頻模態(tài);而采用VMD分解則可分解出1.6 Hz,2.6 Hz,3.6 Hz 3種分量。分析EMD分解結(jié)果,無法確定目標的進動和章動頻率;而分析VMD分解結(jié)果,最高頻應(yīng)是進動與章動頻率之和,為3.6 Hz,則1.6 Hz及2.6 Hz中必有一個是進動或章動頻率,若是1.6 Hz,另一頻率應(yīng)為2.0 Hz,不可能出現(xiàn)2.6 Hz這一分量,從而分析出進動和章動頻率分別為1.0 Hz和2.6 Hz。

目標4為錐體目標,不同于桿狀目標,散射點模型組成有錐鼻和錐尾邊緣之分。按照本文所提分離算法對目標4的m-D曲線進行分離,并與傳統(tǒng)m-D曲線分離算法進行比較,得到如圖10～圖12所示曲線分離結(jié)果。觀察圖12同樣可以發(fā)現(xiàn)本文算法與傳統(tǒng)算法的區(qū)別,在m-D曲線交疊處,傳統(tǒng)算法極易出現(xiàn)錯誤匹配現(xiàn)象,導(dǎo)致分離結(jié)果出錯;而本文算法則能正確分離出交疊的曲線,方便后續(xù)微動特征的提取。

圖10 原始 m-D曲線(目標4)Fig.10 Original m-D curve (target 4)

圖11 二值化及骨架提取得到的m-D曲線(目標4)Fig.11 m-D curve obtained by binarization and skeleton extraction (target 4)

圖12 傳統(tǒng)算法與本文算法的分離結(jié)果(目標4)Fig.12 Separation effect by traditional algorithms and proposed algorithm (target 4)

利用VMD及EMD算法對分離出的錐鼻及錐尾邊緣散射點的m-D曲線進行分解,得如圖13和圖14所示結(jié)果。

圖13 目標4錐鼻散射點m-D曲線EMD及VMD分解結(jié)果Fig.13 m-D curves EMD and VMD decomposition results of cone-nose scatterer of target 4

結(jié)合表1,若所有散射點均存在自旋、進動和章動3種微動,則目標4每個散射點的m-D曲線理論上可分解出0.2 Hz,1.2 Hz,1.4 Hz,1.6Hz,2.6 Hz,2.8 Hz,3.0 Hz,4.0 Hz,4.2 Hz,5.4 Hz 10種模態(tài)。但錐鼻處散射點在其自旋軸上,可不考慮自旋對其影響,其m-D曲線應(yīng)分解出1.4Hz,2.8 Hz,4.2 Hz 3種模態(tài)。觀察圖13可以發(fā)現(xiàn),錐鼻散射點m-D曲線EMD分解可勉強分解出1.4 Hz,其他模態(tài)均混疊嚴重。相較而言,VMD分解能夠很好地分解出1.4 Hz,2.8 Hz,4.2 Hz 3種模態(tài)。根據(jù)VMD分解結(jié)果,可計算出目標4的進動和章動頻率分別為1.4 Hz和2.8 Hz。再觀察圖15,發(fā)現(xiàn)EMD分解結(jié)果混疊嚴重,VMD算法可分解出1.2 Hz,1.6 Hz,2.6Hz,3.0 Hz,4.0 Hz,4.2 Hz,5.4 Hz 7種清晰模態(tài),結(jié)合之前計算出的目標的進動和章動頻率,可以很快計算出目標的自旋頻率為1.2 Hz。

觀察目標3和目標4的EMD及VMD分解結(jié)果,可以發(fā)現(xiàn)無論是EMD還是VMD分解,都無法完全分解出所有的變量,原因可能如下:① VMD算法自身存在低頻分量無法分解的問題;② 組成各m-D曲線的模態(tài)較為接近,VMD分解過程中可能將相近的模態(tài)合并;③ 曲線分離過程中,交叉區(qū)域使用線性插值估計中間值與原始值存在偏差,導(dǎo)致提取出的m-D曲線存在誤差,從而遺失部分分量。

為了更好地體現(xiàn)曲線分離過程中產(chǎn)生的誤差,以目標4的分離結(jié)果進行誤差分析,結(jié)果如圖15和圖16所示。

圖15 曲線分離結(jié)果與原始圖比較Fig.15 Comparison of curve separation results with original graph

觀察圖15可以發(fā)現(xiàn),本文提出的曲線分離方法得出的曲線在極值處與原始曲線存在一定的誤差,圖16所示結(jié)果也很好地顯示出了這一現(xiàn)象。但觀察圖16可以發(fā)現(xiàn),本文所提曲線分離方法誤差基本在[-1.5,1]區(qū)間內(nèi),屬于正常范圍,因此本文所提方法能有效實現(xiàn)曲線分離。

圖16 誤差曲線Fig.16 Error curve

為了進一步測試算法穩(wěn)定性和計算時間復(fù)雜度,對比分析各算法在不同信噪比下的均方根誤差及運行時間。如圖17所示,在信噪比高于-15 dB時,本文算法均能有效分離不同目標的m-D曲線。這說明此時算法對信噪比不敏感,其原因是本文算法的曲線分離過程是基于骨架提取進行的,當信噪比高于-15 dB時,本文所提方法均能完整地提取目標m-D曲線骨架;但是當信噪比低于-15 dB后,由于骨架提取出現(xiàn)大量斷點,導(dǎo)致后續(xù)分離過程產(chǎn)生偽交疊區(qū)間,進而無法有效分離m-D曲線。反觀另外幾種分離算法,可以發(fā)現(xiàn),Viterbi算法在信噪比大于0 dB時,均方根誤差急劇減小,但仍比本文算法大;最近鄰原則法前期處理同樣需要提取m-D曲線骨架,對信噪比敏感性也不高,但均方根誤差極大,表明分離錯誤率較高;峰值法均方根誤差浮動較大,對曲線能量的要求極高,不適用于各m-D曲線能量相近的情況。

圖17 信噪比分析Fig.17 Signal to noise ratio analysis

算法計算時間復(fù)雜度方面,對每種算法進行500次蒙特卡羅實驗,平均運行時間如表2所示。仿真平臺為64位Windows 10操作系統(tǒng),內(nèi)存16 GB,AMD處理器,內(nèi)核Ryzen 7 4800H,CPU 2.9 GHz,MATLAB R2018b。對比表2所示各種算法的運行時間可以發(fā)現(xiàn),本文算法快于Viterbi算法,慢于最近鄰原則法及峰值法。然而,結(jié)合圖18信噪比分析結(jié)果,本文算法在穩(wěn)定性和實時性方面的幾種算法中是最優(yōu)的。

表2 算法時間復(fù)雜度

5 結(jié) 論

本文針對章動目標,通過等效散射點模型及回波模型獲得回波信號及m-D曲線。針對m-D曲線交疊和分離問題,提出一種基于曲線趨勢估計的分離算法。經(jīng)實驗驗證,所提方法可有效解決曲線交疊問題,正確分離m-D曲線。然后,利用VMD算法對分解出的每一條m-D曲線進行分解,將分解結(jié)果與EMD分解結(jié)果進行對比,發(fā)現(xiàn)VMD算法分解效果更佳,混疊程度更低,更加有利于獲得目標的自旋、進動及章動頻率。最后,通過信噪比及時間復(fù)雜度分析,發(fā)現(xiàn)算法能在信噪比高于-15 dB時實現(xiàn)m-D曲線的快速穩(wěn)定有效分離。