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

王 歡 李 峰 宋思盛 姜 洋

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

0 引言

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

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

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

1 信號回波模型

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

(1)

2 微動信號分解

2.1 重構時頻表征

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

(2)

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

(3)

其中H=diag{W,…,W}；WN×N是傅里葉變換矩陣；N=Mms,其中ms=1是滑動窗口的步數；Q=[P1…Pm…PM]T是滑動的窗函數，其中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代表轉置。根據式(3)，采樣缺損時的信號為

s=BGf+n

(4)

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

(5)

(6)

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

(7)

(8)

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

(9)

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

圖1 缺失采樣信號和重構的時頻分布

2.2 瞬時頻率軌跡的提取

將重構的時頻譜分布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)

其中α代表閾值參數。為了避免交叉點附近極值點對后續算法的影響，采用式(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]的啟發，引入卡爾曼濾波[20]和匈牙利算法[21]相結合的多目標跟蹤算法將時頻平面的極大值點關聯成每個散射點的瞬時頻率軌跡。選取G(tm)=P時刻的極大值點qtm,k作為卡爾曼濾波的初始測量值，瞬時頻率的位置和速度由線性狀態空間描述。為了簡化分析，選用勻速(CV)模型，第k個散射點瞬時頻率軌跡在tm+1時刻的狀態向量為

(12)

(13)

(14)

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

(15)

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

3 仿真實驗

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

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

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

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

(25)

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

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

4 結束語

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