基于TRIAA與壓縮感知的ISAR去微多普勒方法

張榮政，王 勇

（哈爾濱工業大學電子工程技術研究所，哈爾濱 150001）

引 言

隨著雷達技術的發展，雷達的信號質量和分辨率逐漸提升，傳統逆合成孔徑雷達（Inverse synthesis aperture radar，ISAR）成像技術已無法滿足清晰化、精細化的成像需求。轉臺成像模型通常假設所有的散射點都隨目標整體一起運動，但是ISAR成像的常見目標上經常有些獨立運動的部件，例如飛機的螺旋槳葉片。這些部件相對于目標整體高速旋轉，造成了回波多普勒譜的展寬，形成了所謂的微多普勒效應［1?4］。若直接采用傳統的ISAR成像算法，由旋轉部件造成的多普勒譜展寬會在成像結果上形成條帶干擾，從而影響成像質量。

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

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

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

1 含旋轉部件的ISAR信號模型

ISAR系統通常采用發射線性調頻信號再用Dechirp處理的方式提高發射信號的帶寬，從而增加距離分辨率。設載頻為fc，脈沖重復時間為T，雷達信號的發射時刻就是tm=mT，m=0，1，2，…，于是發射信號就可以表示為

式中：rectT p(u)表示長度為T p的矩形窗，t?表示脈沖內的快時間，k表示調頻率，c表示電磁波傳播速度。假設tm時刻第i個散射點距離雷達的相對距離是Ri(tm)，回波就可以表示為

式中A i表示每個散射點的回波強度。經過Dechirp處理以及脈沖壓縮，雷達信號轉化為

式中：λ表示發射信號的波長，ΔRi(tm)表示每個散射點與成像參考點的瞬時距離。假設運動補償已經完成且散射點跨域距離單元的問題可以忽略，距離單元內的每個散射點對應的信號可以表示為

接下來，分析含有高速旋轉部件的距離單元內信號的具體形式。圖1給出了含旋轉散射點的ISAR幾何模型，這里平動補償已經完成，目標整體繞O點轉動構成轉臺模型。某一距離單元內有剛體散射點A以及高速旋轉散射點B，A隨目標整體以角速度ω0轉動，B除了隨目標整體的運動外，還繞B以角速度ωs高速轉動。OA和AB在初始時刻與OX軸的夾角分別為θA和θB。

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

所以散射點A與參考點的距離為

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

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

同理，可得B在距離單元內的信號可表示為

接下來考察信號sB(tm)的頻率，可以得到

于是，發現旋轉散射點B的多普勒頻率以正弦形式變化，所以旋轉散射點對應的信號為正弦調頻信號，這就是所謂的微多普勒效應。而如式（6）所示，隨目標整體運動的散射點A對應單頻信號，因此將信號轉換到頻譜就能將A的信號壓縮為一個獨立的點，這就是距離多普勒方法（RD）成像的基本原理。但是假如距離單元內同時存在高速旋轉散射點B，其對應的正弦調頻信號會在成像結果上形成條帶干擾。

2 時間遞歸迭代自適應以及壓縮感知的微多普勒效應去除方法

根據第1節的介紹可知旋轉部件會引起微多普勒效應，處理微多普勒效應的第一步是時頻分析，文獻［9］采用了STFT時頻分析。文獻［10］利用TRIAA得到了更清晰的時頻分布，又利用剛體散射點與旋轉體散射點信號時頻特征的區別，剔除了含有旋轉體散射點信號的元素。本文采用壓縮感知方法將處理后的時頻分布還原成信號頻譜，直接得到了ISAR結果。

2.1 時間遞歸迭代自適應方法

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

這里窗長為L，因此一共分成了M=N-L+1段。接下來根據窗長設置頻率導向矢量

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

式中|α(ωk)|2表示信號在頻率ωk的功率。接下來，采用加權最小二乘估計來估計復振幅α(ωk)。首先構造優化問題

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

應用矩陣求逆引理，最終式（14）簡化為

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

TRIAA的具體步驟如下：

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

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

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

步驟4將每段數據的頻譜拼接為一個矩陣，結果就是用TRIAA方法獲得的時頻分布

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

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

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

這里

于是，某個時刻m的瞬時頻譜可以表示為

式中

這樣，整個TRIAA過程就可以表示為

這里向量x表示完整的時域信號

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

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

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

設X為x的傅里葉變換，最終獲得如下關系

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

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

接下來，設置一個門限Q，剔除ψk(m)絕對值較大的那部分元素，若M Q=int[M(1-Q/100)]表示剩余的元素個數，M-M Q表示從ψk(m)中剔除的元素個數。根據第1節對距離單元內信號的分析可知，高速旋轉散射點對應正弦調頻信號，剛體散射點對應單頻信號。于是經過此番處理，TRIAA上剩余的元素就只包含了剛體散射點產生單頻信號。

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

令TRIAACS表示刪除含干擾元素后的TRIAA，再省略AFULL上對應的行，就得到恢復矩陣ACS，并且滿足如下關系

因此，得到如下最小化問題

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

對于消除的比例Q，有兩種設定的方法。第一種選擇是設置一個固定的消除百分比，只要消除比例Q足夠大，可以消除所有范圍內的干擾。當微多普勒效應所占能量較大時，所需的Q也較大。當計算時頻分布的窗口長度較長時，正弦調頻信號的頻率在每個窗口內變化很大，導致時頻聚集性變低，所需的Q也更大。

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

最終信號分離方法的具體步驟如下：

步驟1計算信號的TRIAA，在計算m時刻的頻譜時，記錄對應的協方差矩陣R m。

步驟2根據TRIAA上每個元素的絕對值大小沿時間重排。設置門限，剔除絕對值大的那部分元素，從而得到TRIAACS。

步驟3用記錄下來的R m，構造矩陣W。計算再根據TRIAACS構造恢復矩陣ACS。

步驟4解決式（33）中的最小化問題，得到去除干擾后信號的頻譜。

3 實驗與分析

3.1 仿真信號分離

首先用仿真信號實驗展示本文的信號分離方法。仿真一個多個單頻分量和多個正弦調頻分量組成的信號

3個單頻信號的頻率分別為50、30和10 Hz；兩個正弦調頻信號的調頻率都是90，重復頻率分別為1.5和1 Hz，信號的采樣點數為256，這里設置采樣頻率為1 Hz，信噪比為20 dB。信號分離的過程如圖2所示。

圖2 基于TRIAA和壓縮感知的信號分離實驗Fig.2 Signal separation experiment based on TRIAA and compressive sensing

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

然后就能用正交匹配追蹤恢復清除干擾后信號的頻譜。這里對比本文方法與文獻［9］中基于ST?FT的信號分離方法?；謴偷念l譜結果如圖3所示。圖3（a）表示使用文獻［9］中的信號分離方法的結果，可見得到的頻譜仍然包含部分干擾。圖3（b）表示利用本文方法重構后的信號?？梢姳疚姆椒ǖ念l譜干擾幾乎全部被剔除。另外，3個單頻信號的幅度應該是相近的，而圖3（a）3個峰出現明顯區別，這說明有效信號在處理中發生了損失。從圖3（b）可見本文方法的信號損失要小很多，這說明壓縮感知方法可以更好地恢復信號。

圖3 信號分離結果Fig.3 Signal separation results

3.2 實測數據ISAR成像實驗

接下來，將本文方法應用于實測的雷達回波。具體流程圖如4所示。

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

為了驗證所提方法的有效性，使用本方法對含有微多普勒干擾的飛機實測數據進行ISAR成像實驗。這里，雷達載頻為5 520 MHz，帶寬為400 MHz，脈沖重復頻率為400 Hz。飛機兩側機翼有一對螺旋槳，可見兩條明顯的條帶干擾。常規RD算法的成像結果如圖5所示。

圖5表示An?26飛機的RD成像結果，可以看到圖上有兩條嚴重的條帶干擾，這是由于兩個螺旋槳的旋轉形成的微多普勒效應。然后使用文獻［9］方法以及本文方法來處理干擾最顯著的第65個距離單元內的信號。實驗結果如圖6所示。

圖5 An?26飛機RD成像Fig.5 RD imaging of An?26

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

圖6 第65個距離單元中的信號Fig.6 Signal in the 65th range cell

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

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

4 結束語

ISAR成像的目標上經常存在螺旋槳葉片等高速運動散射點，這些部件形成的微多普勒效應會在成像結果上形成條帶干擾，從而影響成像質量。本文提出了一種基于TRIAA和壓縮感知技術的ISAR成像去微多普勒方法。仿真和實測數據的實驗結果驗證了本文方法的有效性。