基于匈牙利算法的多分量微動信號提取方法

王 歡 李 峰 宋思盛 姜 洋

(1.西安電子工程研究所 西安 710100;2.空軍工程大學 西安 710077)

0 引言

在雷達探測彈道目標場景下，彈道導(dǎo)彈在飛行過程中除了平動飛行以外，同時會繞自身對稱軸做自旋運動并有可能繞某一定向軸作錐旋運動。這些精細的運動會對雷達回波產(chǎn)生除平動多普勒以外的附加頻率調(diào)制。這一現(xiàn)象被稱作微多普勒效應(yīng)[1]。微多普勒為目標識別[2]和雷達成像[3]的研究提供了有效信息。在過去的十年中，窄帶雷達彈道目標微動特征提取已得到了深入的研究[4]。

在實際應(yīng)用中，雷達回波采樣的缺損和微動信號在時頻域交疊已經(jīng)成為阻礙微動特征有效提取的兩個突出問題。在回波受到強干擾期間，采樣的數(shù)據(jù)不得不被丟棄以免影響后續(xù)的數(shù)據(jù)處理，當雷達出現(xiàn)故障和異常值時也會影響回波信號的有效采集。盡管基于稀疏重構(gòu)的方案為機動目標和穩(wěn)定運動目標提供了缺損數(shù)據(jù)恢復(fù)的解決途徑，比如貝葉斯壓縮感知[5]和自適應(yīng)稀疏分數(shù)模糊函數(shù)[6]，然而由于微動信號的瞬時頻率會隨時間快速變化，上述方法不能有效地應(yīng)用于微動目標。文獻[7]利用一種非參數(shù)聯(lián)合時頻字典從損壞的微動信號中重構(gòu)出聚焦良好的時頻分布，但是沒有提供針對具有多散射點微動目標信號的提取方法。多散射點微動目標的信號提取實質(zhì)上可以當作多分量微動信號分解問題處理，一般情況下，基于多分量微動信號分解的參數(shù)化方法都會提前設(shè)定微動的瞬時多普勒頻率多項式曲線[8]或者正弦曲線[9]模型，但是彈道目標微動信號有可能呈現(xiàn)出更為復(fù)雜的調(diào)制特性(例如滑動散射點的回波信號)，需要設(shè)置多個參數(shù)來擬合調(diào)制特性，從而增加了算法的復(fù)雜度。而大多數(shù)非參數(shù)化的方法，例如變分模態(tài)分解[10]、結(jié)合短時傅里葉變換的同步擠壓變換[11]、稀疏短時傅里葉變換[12]，都無法處理微動信號在時頻域出現(xiàn)交疊的情況，然而彈道目標的多分量微動信號在時頻域常常出現(xiàn)交疊的情況。文獻[13]提出一種嶺跡重排的方法提取時頻域的瞬時頻率曲線，并結(jié)合本征線性調(diào)頻分量分解(ICCD)算法[14]有效分解了在時頻域分解重疊的微動信號，但是在時頻分布因采樣缺損出現(xiàn)散焦和缺損的情況下，嶺跡重排算法無法有效運行；短時變分模態(tài)分解(STVMD)[15]算法可以在瞬時頻率不規(guī)則調(diào)制和低信噪比情況下正確分解具有交疊狀態(tài)的微動信號，但是當微動信號在時頻域完全重疊時，分解效果不佳。

為了解決上述問題，本文提出一種采樣缺失條件下的多分量微動信號分解方法，首先將采樣缺失微動信號的時頻分布重構(gòu)問題轉(zhuǎn)換為一個稀疏正則化問題，并采用軟閾值迭代算法(ISTA)求解；其次提取重構(gòu)時頻分布的極大值點，將極大值點關(guān)聯(lián)為瞬時頻率軌跡的問題構(gòu)造為指派問題；最后采用卡爾曼濾波和匈牙利算法求解，得到每個散射點的瞬時頻率軌跡，實現(xiàn)了雷達采樣缺失條件下對多分量微動信號的分解。

1 信號回波模型

在高頻電磁條件下，金屬目標可以由點散射模型表示[7]，回波信號缺損條件下的微動目標窄帶雷達回波信號表示為

(1)

2 微動信號分解

2.1 重構(gòu)時頻表征

s(t)是一個典型的非平穩(wěn)信號。時頻分析是一種在時頻域表征非平穩(wěn)信號能量的有效方法，非平穩(wěn)信號可以在時頻域上表征出頻域中無法表示的特征。短時傅里葉變換[16]具備線性和沒有交叉項等諸多優(yōu)勢，因此本文采用短時傅里葉變換分析回波信號。短時傅里葉變換可以表示為

(2)

其中TF0(t,f)代表時頻表征；f代表瞬時頻率；gσ代表長度為σ的窗函數(shù)，用來保證加窗信號足夠平穩(wěn)，窗函數(shù)有多種形式，本文使用高斯窗，它可以很容易地用其他形式重寫，比如矩形窗和漢明窗等。為了方便描述重構(gòu)算法，結(jié)合式(1)，將式(2)轉(zhuǎn)換為矩陣的形式[12]

(3)

其中H=diag{W,…,W}；WN×N是傅里葉變換矩陣；N=Mms,其中ms=1是滑動窗口的步數(shù)；Q=[P1…Pm…PM]T是滑動的窗函數(shù)，其中Pm=diag{pm(1),…pm(n),…pm(N)},n∈[1,N]是對角矩陣；B=diag{b(1),…b(n),…b(N)},b(n)∈{0,1}是表示信號缺損的對角矩陣；沒有缺損的理想回波信號為s0=[s0(1),…s0(n),…s0(N)]T，(·)T代表轉(zhuǎn)置。根據(jù)式(3)，采樣缺損時的信號為

s=BGf+n

(4)

其中G是HQ的偽逆矩陣；n=[η(1),…η(n),…η(N)]T是加性噪聲。為了從幅度起伏和采樣缺損的回波s中得到理想的時頻表征f，令s=[s(1),…s(n),…s(N)]T，式(4)可以被描述成一個優(yōu)化問題

(5)

(6)

其中‖·‖1表示1范數(shù)運算符。L1正則化模型可以被認為是一個套索回歸問題，可以采用ISTA[17]求解。ISTA是一種近端梯度下降法，鄰近算子為

(7)

(8)

soft(z,λ1)=sign(z)max{0,|z|-λ1}

(9)

由圖1所示，通過軟閾值迭代法，得到了缺失采樣條件下的稀疏、去噪且銳化的時頻譜分布。

圖1 缺失采樣信號和重構(gòu)的時頻分布

2.2 瞬時頻率軌跡的提取

將重構(gòu)的時頻譜分布f取模并重排為由N個時刻頻譜組成的時頻圖TF(t,f)。TF(t,f)中極大值點的位置坐標為qt,g(t)=(t,fg(t))，其中t=t0,…,tN-1，fg(t)代表第t時刻的第g(t)個極值點在頻譜上的坐標，t時刻的頻譜總共有G(t)個極值點。極大值點的坐標qt,g(t)由式(10)得到[18]

(10)

其中α代表閾值參數(shù)。為了避免交叉點附近極值點對后續(xù)算法的影響，采用式(11)消除交叉點的極值點

|fg(t)-fg_other(t)|≥df

(11)

其中fg_other(t)代表第t個頻譜除了g(t)以外的其他極大值點在頻譜上的位置，將不滿足式(11)的極大值點剔除，得到qt,g(t)，其在時頻平面的位置由圖2所示，在本文中df=N/20。

圖2 提取的極大值點在時頻平面的位置

受到Simple online and realtime tracking (Sort)算法[19]的啟發(fā)，引入卡爾曼濾波[20]和匈牙利算法[21]相結(jié)合的多目標跟蹤算法將時頻平面的極大值點關(guān)聯(lián)成每個散射點的瞬時頻率軌跡。選取G(tm)=P時刻的極大值點qtm,k作為卡爾曼濾波的初始測量值，瞬時頻率的位置和速度由線性狀態(tài)空間描述。為了簡化分析，選用勻速(CV)模型，第k個散射點瞬時頻率軌跡在tm+1時刻的狀態(tài)向量為

(12)

(13)

(14)

其中ygk=0代表第g(tm+1)個測量值與第k測量值不匹配；ygk=1代表g(tm+1)個測量值與第k測量值匹配。指派問題的核心問題是如何構(gòu)造效率矩陣，如式(14)所示的指派問題，其效率矩陣Cm+1可以通過高斯核密度函數(shù)構(gòu)造。

(15)

圖3 各散射點的瞬時頻率軌跡

3 仿真實驗

為了證明本文所提算法的有效性，通過計算機仿真進行驗證，并對仿真結(jié)果進行分析。設(shè)雷達的載頻為10GHz，脈沖重復(fù)頻率為512Hz，觀測時間為1s，在加性高斯白噪聲環(huán)境下，信噪比為10dB，短時傅里葉變換的高斯窗函數(shù)長度σ=53。一般情況下，參數(shù)λ1的值越大，ISTA算法生成的時頻分布就越稀疏；參數(shù)α的值越大，從時頻圖中提取的極大值點越稀疏；參數(shù)L的值越大，ICCD算法時頻濾波器組的寬度越大。根據(jù)經(jīng)驗，當參數(shù)值λ1設(shè)置在0.01～0.04之間，α設(shè)置在0.4～0.6之間，L設(shè)置在5～15之間，λ2設(shè)置在4 ～ 6之間時算法可以有效運行。在本文中參數(shù)λ、α、L、λ2分別設(shè)置為0.01、0.5、5、5。

將本文算法應(yīng)用于基于滑動散射點的微動信號分解中。當彈道目標在空間上進動時，一些散射點的位置隨著目標的運動而滑動，稱為“滑動散射點”[22]。設(shè)彈道目標的錐旋轉(zhuǎn)頻率為0.5Hz，錐底半徑為0.5m，錐旋角為10°,從錐頂?shù)綀A錐質(zhì)量中心的距離為2m，從圓錐質(zhì)量中心到圓錐底部中心的距離為0.5m，雷達入射角度為45°。數(shù)據(jù)隨機缺損40%，基于滑動散射模型的缺損采樣信號的時頻分布如圖4(a)所示，可見時頻分布是破損且散焦的；結(jié)合ICCD算法[22]，本文算法分解缺損采樣條件下的滑動散射點信號分量如圖4(b)、圖4(c)、圖4(d)所示。結(jié)果表明，該方法能夠有效地分解出各滑動散射點的微動信號分量。

圖4 基于滑動散射模型的微動信號分解

最后，采用參數(shù)δ來評價兩種方法的估計精度

(25)

表1列出了兩種方法在采樣缺損30%和40%條件下分解理想散射目標回波和滑動散射目標回波的精度，可以看出：一方面，采樣缺損率越大，兩種方法分解得到的微動信號精度越低，這是因為隨著采樣缺損率的增加，時頻圖的散焦程度更加嚴重，導(dǎo)致了兩種算法的性能下降；另一方面，本文方法比STVMD算法分解信號的精度平均提高了1dB左右，這是因為當某個散射點信號與其他散射點信號在時頻域交疊時，STVMD算法可能會將該散射點的信號能量分解到其他散射點的回波中，導(dǎo)致出現(xiàn)交疊時刻信號丟失的情況。

表1 不同方法分解微動信號的精度

4 結(jié)束語

本文提出了一種微動信號采樣缺損和回波時頻域交疊情況下的多分量微動目標回波分解新方法。首先在構(gòu)建缺損采樣多分量微動信號稀疏正則化模型的基礎(chǔ)上，采用ISTA算法對時頻分布進行重構(gòu)；其次結(jié)合卡爾曼濾波和匈牙利算法將重構(gòu)時頻分布極大值關(guān)聯(lián)為瞬時頻率軌跡；最后仿真結(jié)果表明本文相對于STVMD算法方法分解得到的微動信號精度更高。應(yīng)該指出的是，本文方法作為時頻分布的一種后處理方法，其性能受限于時頻分布的分辨率，因此需要選擇合適的時頻分辨率以保證算法的有效運行；其次，與許多已有的多分量信號分解方法相同，本文方法也需要預(yù)先設(shè)定待分解信號中信號分量的個數(shù)。在下一步工作中，我們將重點尋找更有效的字典矩陣以重構(gòu)時頻分布，并改進現(xiàn)有算法使其在不需要預(yù)先設(shè)定信號分量個數(shù)的前提下有效運行。