張榮政，王 勇

（哈爾濱工業(yè)大學(xué)電子工程技術(shù)研究所，哈爾濱 150001）

引 言

隨著雷達(dá)技術(shù)的發(fā)展，雷達(dá)的信號(hào)質(zhì)量和分辨率逐漸提升，傳統(tǒng)逆合成孔徑雷達(dá)（Inverse synthesis aperture radar，ISAR）成像技術(shù)已無(wú)法滿(mǎn)足清晰化、精細(xì)化的成像需求。轉(zhuǎn)臺(tái)成像模型通常假設(shè)所有的散射點(diǎn)都隨目標(biāo)整體一起運(yùn)動(dòng)，但是ISAR成像的常見(jiàn)目標(biāo)上經(jīng)常有些獨(dú)立運(yùn)動(dòng)的部件，例如飛機(jī)的螺旋槳葉片。這些部件相對(duì)于目標(biāo)整體高速旋轉(zhuǎn)，造成了回波多普勒譜的展寬，形成了所謂的微多普勒效應(yīng)［1?4］。若直接采用傳統(tǒng)的ISAR成像算法，由旋轉(zhuǎn)部件造成的多普勒譜展寬會(huì)在成像結(jié)果上形成條帶干擾，從而影響成像質(zhì)量。

針對(duì)這一問(wèn)題，有一類(lèi)解決方法是在成像前預(yù)先從回波中分離由剛體散射點(diǎn)形成的信號(hào)和由高速旋轉(zhuǎn)散射點(diǎn)形成的信號(hào)。文獻(xiàn)［5］提出一種基于低調(diào)頻率匹配濾波的信號(hào)分離方法，用搜索的方式挑選出有效的剛體信號(hào)。文獻(xiàn)［6］提出了一種基于多重測(cè)量矢量（Multiple measurement vectors，MMV）和壓縮感知（Compressive sensing，CS）的含旋轉(zhuǎn)部件目標(biāo)ISAR成像方法，根據(jù)剛體與旋轉(zhuǎn)部件信號(hào)的MMV模型稀疏性的不同分離信號(hào)。文獻(xiàn)［7］使用Radon?Wigner變換方法對(duì)信號(hào)的Wigner?Ville分布進(jìn)行直線線積分，以此識(shí)別散射點(diǎn)類(lèi)型，再利用CLEAN技術(shù)完成了信號(hào)分離。文獻(xiàn)［8］利用復(fù)數(shù)局部均值分解將復(fù)雜信號(hào)分解成若干個(gè)單分量信號(hào)，從而分離出微多普勒信號(hào)。

另一類(lèi)方法是將微多普勒信號(hào)視作干擾直接濾除，再進(jìn)行ISAR成像。文獻(xiàn)［9］提出了一種基于短時(shí)傅里葉變換（Short?time Fourier transform，STFT）時(shí)頻分布和L?統(tǒng)計(jì)量的ISAR成像去除微多普勒效應(yīng)方法，該方法利用剛體和旋轉(zhuǎn)散射點(diǎn)時(shí)頻特征的不同，從時(shí)頻圖上刪除代表旋轉(zhuǎn)部件那部分元素，再將處理后的時(shí)頻圖還原成信號(hào)并完成ISAR成像。文獻(xiàn)［10］利用時(shí)間遞歸迭代自適應(yīng)方法（Time re?cursive iterative adaptive approach，TRIAA）［11?12］替換原方法的STFT時(shí)頻分析，從而提高了時(shí)頻分析的分辨率，更有效地去除了微多普勒效應(yīng)。文獻(xiàn)［13?15］中，壓縮感知技術(shù)被用于信號(hào)還原過(guò)程，從而在原方法基礎(chǔ)上提升了成像質(zhì)量。

本文基于上述研究工作，首先利用TRIAA提高時(shí)頻分析的頻率分辨率，從而更有效去除干擾，再使用壓縮感知稀疏重構(gòu)技術(shù)還原去除干擾后的信號(hào)，最終提出了一種基于TRIAA和壓縮感知的ISAR成像去除微多普勒效應(yīng)方法，在原方法基礎(chǔ)上提高了ISAR成像質(zhì)量。

1 含旋轉(zhuǎn)部件的ISAR信號(hào)模型

ISAR系統(tǒng)通常采用發(fā)射線性調(diào)頻信號(hào)再用Dechirp處理的方式提高發(fā)射信號(hào)的帶寬，從而增加距離分辨率。設(shè)載頻為fc，脈沖重復(fù)時(shí)間為T(mén)，雷達(dá)信號(hào)的發(fā)射時(shí)刻就是tm=mT，m=0，1，2，…，于是發(fā)射信號(hào)就可以表示為

式中：rectT p(u)表示長(zhǎng)度為T(mén) p的矩形窗，t?表示脈沖內(nèi)的快時(shí)間，k表示調(diào)頻率，c表示電磁波傳播速度。假設(shè)tm時(shí)刻第i個(gè)散射點(diǎn)距離雷達(dá)的相對(duì)距離是Ri(tm)，回波就可以表示為

式中A i表示每個(gè)散射點(diǎn)的回波強(qiáng)度。經(jīng)過(guò)Dechirp處理以及脈沖壓縮，雷達(dá)信號(hào)轉(zhuǎn)化為

式中：λ表示發(fā)射信號(hào)的波長(zhǎng)，ΔRi(tm)表示每個(gè)散射點(diǎn)與成像參考點(diǎn)的瞬時(shí)距離。假設(shè)運(yùn)動(dòng)補(bǔ)償已經(jīng)完成且散射點(diǎn)跨域距離單元的問(wèn)題可以忽略，距離單元內(nèi)的每個(gè)散射點(diǎn)對(duì)應(yīng)的信號(hào)可以表示為

接下來(lái)，分析含有高速旋轉(zhuǎn)部件的距離單元內(nèi)信號(hào)的具體形式。圖1給出了含旋轉(zhuǎn)散射點(diǎn)的ISAR幾何模型，這里平動(dòng)補(bǔ)償已經(jīng)完成，目標(biāo)整體繞O點(diǎn)轉(zhuǎn)動(dòng)構(gòu)成轉(zhuǎn)臺(tái)模型。某一距離單元內(nèi)有剛體散射點(diǎn)A以及高速旋轉(zhuǎn)散射點(diǎn)B，A隨目標(biāo)整體以角速度ω0轉(zhuǎn)動(dòng)，B除了隨目標(biāo)整體的運(yùn)動(dòng)外，還繞B以角速度ωs高速轉(zhuǎn)動(dòng)。OA和AB在初始時(shí)刻與OX軸的夾角分別為θA和θB。

圖1 含有旋轉(zhuǎn)散射點(diǎn)的ISAR幾何模型Fig.1 ISAR imaging model with rotat?ing parts

所以散射點(diǎn)A與參考點(diǎn)的距離為

式中ΔR0為轉(zhuǎn)臺(tái)原點(diǎn)O到參考點(diǎn)的距離。由于成像時(shí)間內(nèi)，目標(biāo) 整 體 的 旋 轉(zhuǎn) 角 度ω0tm非 常 小，因 此sin(ω0tm)≈ω0tm，cos(ω0tm)≈1。代入式（5），A在距離單元內(nèi)的信號(hào)可表示為

sA(tm)顯然是一個(gè)單頻信號(hào)，下面用同樣的方法可以得到散射點(diǎn)B到參考點(diǎn)的距離為

同理，可得B在距離單元內(nèi)的信號(hào)可表示為

接下來(lái)考察信號(hào)sB(tm)的頻率，可以得到

于是，發(fā)現(xiàn)旋轉(zhuǎn)散射點(diǎn)B的多普勒頻率以正弦形式變化，所以旋轉(zhuǎn)散射點(diǎn)對(duì)應(yīng)的信號(hào)為正弦調(diào)頻信號(hào)，這就是所謂的微多普勒效應(yīng)。而如式（6）所示，隨目標(biāo)整體運(yùn)動(dòng)的散射點(diǎn)A對(duì)應(yīng)單頻信號(hào)，因此將信號(hào)轉(zhuǎn)換到頻譜就能將A的信號(hào)壓縮為一個(gè)獨(dú)立的點(diǎn)，這就是距離多普勒方法（RD）成像的基本原理。但是假如距離單元內(nèi)同時(shí)存在高速旋轉(zhuǎn)散射點(diǎn)B，其對(duì)應(yīng)的正弦調(diào)頻信號(hào)會(huì)在成像結(jié)果上形成條帶干擾。

2 時(shí)間遞歸迭代自適應(yīng)以及壓縮感知的微多普勒效應(yīng)去除方法

根據(jù)第1節(jié)的介紹可知旋轉(zhuǎn)部件會(huì)引起微多普勒效應(yīng)，處理微多普勒效應(yīng)的第一步是時(shí)頻分析，文獻(xiàn)［9］采用了STFT時(shí)頻分析。文獻(xiàn)［10］利用TRIAA得到了更清晰的時(shí)頻分布，又利用剛體散射點(diǎn)與旋轉(zhuǎn)體散射點(diǎn)信號(hào)時(shí)頻特征的區(qū)別，剔除了含有旋轉(zhuǎn)體散射點(diǎn)信號(hào)的元素。本文采用壓縮感知方法將處理后的時(shí)頻分布還原成信號(hào)頻譜，直接得到了ISAR結(jié)果。

2.1 時(shí)間遞歸迭代自適應(yīng)方法

對(duì)于一個(gè)距離單元內(nèi)的離散信號(hào)x(n)，n=1，2，…，N，可以用滑窗分割的方式將其分段為

這里窗長(zhǎng)為L(zhǎng)，因此一共分成了M=N-L+1段。接下來(lái)根據(jù)窗長(zhǎng)設(shè)置頻率導(dǎo)向矢量

這里ωk=2πk/K，k=0，1，…，K-1，K表示頻率采樣點(diǎn)數(shù)，一般K>L。α(ωk)表示sL(m)在頻率ωk處的復(fù)振幅，其協(xié)方差矩陣R就可以寫(xiě)成

式中|α(ωk)|2表示信號(hào)在頻率ωk的功率。接下來(lái)，采用加權(quán)最小二乘估計(jì)來(lái)估計(jì)復(fù)振幅α(ωk)。首先構(gòu)造優(yōu)化問(wèn)題

式中Q k=R-|α(ωk)|2f(ωk)f*(ωk)表示廣義噪聲協(xié)方差矩陣，其物理意義在于假設(shè)信號(hào)在頻率ωk處的復(fù)振幅是α(ωk)，把信號(hào)中除頻率ωk外成分都視為噪聲，因此當(dāng)估計(jì)復(fù)振幅時(shí)加權(quán)矩陣應(yīng)取Q-1k。當(dāng)加權(quán)的剩余平方和最小時(shí)，就得到了復(fù)振幅的估計(jì)值。具體細(xì)節(jié)可參考文獻(xiàn)［16?17］。代入最小二乘估計(jì)解的一般形式，該優(yōu)化問(wèn)題的最優(yōu)解為

應(yīng)用矩陣求逆引理，最終式（14）簡(jiǎn)化為

通過(guò)以上推導(dǎo)可以看出，首先設(shè)|α(ωk)|2=1，k=1，2，…，K，代入式（12）和式（15）得到初步結(jié)果再將結(jié)果代回式（12）和式（15）并反復(fù)迭代多次，即可得到x L(m)的頻譜對(duì)滑窗分割形成的每一段數(shù)據(jù)進(jìn)行同樣處理，就能觀察到x(n)的頻譜隨時(shí)間變化的情況，這就是TRI?AA的基本原理?？紤]到相鄰窗之間，如x L(m)和x L(m+1)之間的頻譜差別不大，因此就能將前一段的結(jié)果作為下一段的迭代初始值，從而極大減少迭代次數(shù)。

TRIAA的具體步驟如下：

步驟1將數(shù)據(jù)s(n)用滑窗的方式分段，每段數(shù)據(jù)的具體形式參考式（10），于是得到下列矩陣

步驟2設(shè)|α(ωk)|2=1，k=1，2，…，K，將x L(1)代入式（12）和式（15），得到初步結(jié)果再將和x L(1)代 回 式（12）和 式（15）并 反 復(fù) 迭 代 多 次，即 可 得 到x L(1)的 頻 譜

步驟3設(shè)將x L(2)代入式（12）和式（15），得到x L(2)的頻譜并以此類(lèi)推，設(shè)將x L(m)代入式（12）和式（15），即可得到x L(m)的頻譜一直計(jì)算到x L(M)。

步驟4將每段數(shù)據(jù)的頻譜拼接為一個(gè)矩陣，結(jié)果就是用TRIAA方法獲得的時(shí)頻分布

2.2 基于TRIAA和壓縮感知的微多普勒消除方法

假 設(shè) 待 處 理 的 信 號(hào) 是x(n)，n=1，2，…，N，其 離 散 傅 里 葉 變 換 就 可 以 表 示 為X(k)，k=1，2，…，K。對(duì)其進(jìn)行TRIAA時(shí)頻分析，根據(jù)式（17）和式（18），TRIAA上每個(gè)元素可表示為

根據(jù)式（15），R m為直接得到(ωk)的協(xié)方差矩陣，則有如下形式

這里

于是，某個(gè)時(shí)刻m的瞬時(shí)頻譜可以表示為

式中

這樣，整個(gè)TRIAA過(guò)程就可以表示為

這里向量x表示完整的時(shí)域信號(hào)

TRIAA表示大小為M×K的矩陣，即

而KM×N維的系數(shù)矩陣W可以表示為

這里0K×1表示元素全為零的列向量。由于TRIAA采用滑窗形式，因此W1，W2，…也以同樣的滑窗形式分布。然后，用表示傅里葉逆變換矩陣，代入式（25）就有以下關(guān)系

設(shè)X為x的傅里葉變換，最終獲得如下關(guān)系

對(duì)于信號(hào)x(n)的短時(shí)傅里葉變換TRIAA(m，k)，把頻率為ωk的元素都放在同一個(gè)集合Xk(m)里，有

然后計(jì)算集合X k(m)所有元素的絕對(duì)值，并根據(jù)絕對(duì)值的大小將所有元素從小到大重新排列，得到新集 合ψk(m)∈Xk(m)，滿(mǎn) 足|ψk(0)|≤|ψk(1)|≤|ψk(2)|≤…≤|ψk(M-1)|。時(shí) 頻 圖 上 的 單 頻 信 號(hào) 由 于頻率不隨時(shí)間改變，因此重排后順序幾乎不變。而如正弦調(diào)頻信號(hào)的非平穩(wěn)信號(hào)頻率隨時(shí)間改變，所以其元素會(huì)下沉到絕對(duì)值較大的一側(cè)。

接下來(lái)，設(shè)置一個(gè)門(mén)限Q，剔除ψk(m)絕對(duì)值較大的那部分元素，若M Q=int[M(1-Q/100)]表示剩余的元素個(gè)數(shù)，M-M Q表示從ψk(m)中剔除的元素個(gè)數(shù)。根據(jù)第1節(jié)對(duì)距離單元內(nèi)信號(hào)的分析可知，高速旋轉(zhuǎn)散射點(diǎn)對(duì)應(yīng)正弦調(diào)頻信號(hào)，剛體散射點(diǎn)對(duì)應(yīng)單頻信號(hào)。于是經(jīng)過(guò)此番處理，TRIAA上剩余的元素就只包含了剛體散射點(diǎn)產(chǎn)生單頻信號(hào)。

然而，正弦調(diào)頻信號(hào)和單頻信號(hào)可能會(huì)發(fā)生重疊，這就造成了在時(shí)頻圖上剔除正弦調(diào)頻信號(hào)時(shí)也可能損失一部分單頻信號(hào)?？紤]到ISAR成像時(shí)，目標(biāo)一般滿(mǎn)足散射點(diǎn)模型，每個(gè)距離單元內(nèi)的剛體散射點(diǎn)的個(gè)數(shù)也為有限個(gè)。每個(gè)散射點(diǎn)又只對(duì)應(yīng)一個(gè)單頻信號(hào)，所以濾除正弦調(diào)頻分量后的信號(hào)在頻域上有稀疏性。因此依據(jù)式（30），濾除正弦調(diào)頻分量后的TRIAA可用壓縮感知方式直接還原為濾除正弦調(diào)頻分量后信號(hào)的頻譜X′。這樣也直接完成了ISAR成像的最后一步，既方位向壓縮。

令TRIAACS表示刪除含干擾元素后的TRIAA，再省略AFULL上對(duì)應(yīng)的行，就得到恢復(fù)矩陣ACS，并且滿(mǎn)足如下關(guān)系

因此，得到如下最小化問(wèn)題

解決該最小化問(wèn)題，即可得到濾除正弦調(diào)頻分量后信號(hào)的頻譜X′。而由于TRIAACS上已無(wú)含有干擾的元素，所以頻譜X′只包含剛體信號(hào)。

對(duì)于消除的比例Q，有兩種設(shè)定的方法。第一種選擇是設(shè)置一個(gè)固定的消除百分比，只要消除比例Q足夠大，可以消除所有范圍內(nèi)的干擾。當(dāng)微多普勒效應(yīng)所占能量較大時(shí)，所需的Q也較大。當(dāng)計(jì)算時(shí)頻分布的窗口長(zhǎng)度較長(zhǎng)時(shí)，正弦調(diào)頻信號(hào)的頻率在每個(gè)窗口內(nèi)變化很大，導(dǎo)致時(shí)頻聚集性變低，所需的Q也更大。

另一種更復(fù)雜的方法是為每個(gè)距離單元中的信號(hào)設(shè)置自適應(yīng)的消除百分比［9］，其原理是當(dāng)信號(hào)的時(shí)頻分布重新排列后，包含微多普勒效應(yīng)的部分比只包含單頻信號(hào)的部分能量增加更大。然而，該方法要求信號(hào)的微多普勒效應(yīng)引起的頻帶擴(kuò)展足夠大，且信號(hào)的信噪比足夠低。

最終信號(hào)分離方法的具體步驟如下：

步驟1計(jì)算信號(hào)的TRIAA，在計(jì)算m時(shí)刻的頻譜時(shí)，記錄對(duì)應(yīng)的協(xié)方差矩陣R m。

步驟2根據(jù)TRIAA上每個(gè)元素的絕對(duì)值大小沿時(shí)間重排。設(shè)置門(mén)限，剔除絕對(duì)值大的那部分元素，從而得到TRIAACS。

步驟3用記錄下來(lái)的R m，構(gòu)造矩陣W。計(jì)算再根據(jù)TRIAACS構(gòu)造恢復(fù)矩陣ACS。

步驟4解決式（33）中的最小化問(wèn)題，得到去除干擾后信號(hào)的頻譜。

3 實(shí)驗(yàn)與分析

3.1 仿真信號(hào)分離

首先用仿真信號(hào)實(shí)驗(yàn)展示本文的信號(hào)分離方法。仿真一個(gè)多個(gè)單頻分量和多個(gè)正弦調(diào)頻分量組成的信號(hào)

3個(gè)單頻信號(hào)的頻率分別為50、30和10 Hz；兩個(gè)正弦調(diào)頻信號(hào)的調(diào)頻率都是90，重復(fù)頻率分別為1.5和1 Hz，信號(hào)的采樣點(diǎn)數(shù)為256，這里設(shè)置采樣頻率為1 Hz，信噪比為20 dB。信號(hào)分離的過(guò)程如圖2所示。

圖2 基于TRIAA和壓縮感知的信號(hào)分離實(shí)驗(yàn)Fig.2 Signal separation experiment based on TRIAA and compressive sensing

圖2（a）表示信號(hào)由TRIAA得到的時(shí)頻圖，可以在圖中看到3條直線和2條正弦曲線，分別表示單頻信號(hào)和正弦調(diào)頻信號(hào)。圖2（b）表示將時(shí)頻圖按每個(gè)元素絕對(duì)值的大小沿時(shí)間維重新排列后的結(jié)果。圖2（c）黑色的部分就是含有干擾元素的位置，這里將去除門(mén)限設(shè)置在60%。圖2（d）表示剔除干擾信號(hào)后的時(shí)頻圖。

然后就能用正交匹配追蹤恢復(fù)清除干擾后信號(hào)的頻譜。這里對(duì)比本文方法與文獻(xiàn)［9］中基于ST?FT的信號(hào)分離方法?；謴?fù)的頻譜結(jié)果如圖3所示。圖3（a）表示使用文獻(xiàn)［9］中的信號(hào)分離方法的結(jié)果，可見(jiàn)得到的頻譜仍然包含部分干擾。圖3（b）表示利用本文方法重構(gòu)后的信號(hào)?？梢?jiàn)本文方法的頻譜干擾幾乎全部被剔除。另外，3個(gè)單頻信號(hào)的幅度應(yīng)該是相近的，而圖3（a）3個(gè)峰出現(xiàn)明顯區(qū)別，這說(shuō)明有效信號(hào)在處理中發(fā)生了損失。從圖3（b）可見(jiàn)本文方法的信號(hào)損失要小很多，這說(shuō)明壓縮感知方法可以更好地恢復(fù)信號(hào)。

圖3 信號(hào)分離結(jié)果Fig.3 Signal separation results

3.2 實(shí)測(cè)數(shù)據(jù)ISAR成像實(shí)驗(yàn)

接下來(lái)，將本文方法應(yīng)用于實(shí)測(cè)的雷達(dá)回波。具體流程圖如4所示。

圖4 本文方法ISAR成像流程圖Fig.4 Flow chart of ISAR imaging

為了驗(yàn)證所提方法的有效性，使用本方法對(duì)含有微多普勒干擾的飛機(jī)實(shí)測(cè)數(shù)據(jù)進(jìn)行ISAR成像實(shí)驗(yàn)。這里，雷達(dá)載頻為5 520 MHz，帶寬為400 MHz，脈沖重復(fù)頻率為400 Hz。飛機(jī)兩側(cè)機(jī)翼有一對(duì)螺旋槳，可見(jiàn)兩條明顯的條帶干擾。常規(guī)RD算法的成像結(jié)果如圖5所示。

圖5表示An?26飛機(jī)的RD成像結(jié)果，可以看到圖上有兩條嚴(yán)重的條帶干擾，這是由于兩個(gè)螺旋槳的旋轉(zhuǎn)形成的微多普勒效應(yīng)。然后使用文獻(xiàn)［9］方法以及本文方法來(lái)處理干擾最顯著的第65個(gè)距離單元內(nèi)的信號(hào)。實(shí)驗(yàn)結(jié)果如圖6所示。

圖5 An?26飛機(jī)RD成像Fig.5 RD imaging of An?26

圖6（a）表示對(duì)第65距離單元內(nèi)的信號(hào)使用FFT進(jìn)行方位向壓縮后的結(jié)果，可以發(fā)現(xiàn)微多普勒效應(yīng)形成的頻帶展寬。圖6（b）展示了使用文獻(xiàn)［9］的基于STFT方法消除干擾后的頻譜，這里將去除門(mén)限設(shè)置在50%，可見(jiàn)該方法消除了大部分干擾，但是部分干擾仍然存在。圖6（c）表示使用本文方法后的處理結(jié)果，可見(jiàn)干擾幾乎全部被濾除。接下來(lái)處理全部距離單元，成像結(jié)果如圖7所示。

圖6 第65個(gè)距離單元中的信號(hào)Fig.6 Signal in the 65th range cell

圖7 An?26飛機(jī)去微多普勒ISAR成像Fig.7 An?26 aircraft filtered micro?Doppler ISAR imaging

圖7（a）表示利用文獻(xiàn)［9］中的基于STFT方法的ISAR成像結(jié)果，可見(jiàn)大部分條帶干擾被消除。圖7（b）表示使用本文方法去除微多普勒效應(yīng)的ISAR成像結(jié)果，可見(jiàn)本方法基本完全消除了微多普勒效應(yīng)造成的干擾，且飛機(jī)目標(biāo)的輪廓清晰，效果明顯優(yōu)于文獻(xiàn)［9］中的方法。最后，用成像結(jié)果的圖像熵對(duì)比這兩種方法的效果。圖像熵越小說(shuō)明干擾越少，圖像越清晰?；赟TFT的方法得到的圖像熵是6.365 9，而基于TRIAA和壓縮感知方法的圖像熵是4.778。這表明本文方法的成像效果最好。

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

ISAR成像的目標(biāo)上經(jīng)常存在螺旋槳葉片等高速運(yùn)動(dòng)散射點(diǎn)，這些部件形成的微多普勒效應(yīng)會(huì)在成像結(jié)果上形成條帶干擾，從而影響成像質(zhì)量。本文提出了一種基于TRIAA和壓縮感知技術(shù)的ISAR成像去微多普勒方法。仿真和實(shí)測(cè)數(shù)據(jù)的實(shí)驗(yàn)結(jié)果驗(yàn)證了本文方法的有效性。