王 歡 李 峰 宋思盛 姜 洋

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

0 引言

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

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

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

1 信號(hào)回波模型

在高頻電磁條件下，金屬目標(biāo)可以由點(diǎn)散射模型表示[7]，回波信號(hào)缺損條件下的微動(dòng)目標(biāo)窄帶雷達(dá)回波信號(hào)表示為

(1)

2 微動(dòng)信號(hào)分解

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

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

(2)

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

(3)

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

s=BGf+n

(4)

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

(5)

(6)

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

(7)

(8)

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

(9)

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

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

2.2 瞬時(shí)頻率軌跡的提取

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

(10)

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

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

(11)

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

圖2 提取的極大值點(diǎn)在時(shí)頻平面的位置

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

(12)

(13)

(14)

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

(15)

圖3 各散射點(diǎn)的瞬時(shí)頻率軌跡

3 仿真實(shí)驗(yàn)

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

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

圖4 基于滑動(dòng)散射模型的微動(dòng)信號(hào)分解

最后，采用參數(shù)δ來(lái)評(píng)價(jià)兩種方法的估計(jì)精度

(25)

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

表1 不同方法分解微動(dòng)信號(hào)的精度

4 結(jié)束語(yǔ)

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