低信噪比下快速分裂Bregman 迭代ISAR 超分辨率成像*

朱燕麗，周 鵬?，張振華，張 晰，張 杰，王 影

（1 中國石油大學(xué)（華東）海洋與空間信息學(xué)院 青島 266580 2 北京遙測技術(shù)研究所 北京 100076 3 自然資源部第一海洋研究所 青島 266061）

引言

逆合成孔徑雷達(dá)ISAR 利用非合作運(yùn)動目標(biāo)相對于雷達(dá)的轉(zhuǎn)角進(jìn)行高分辨成像。它能夠全天候、全天時(shí)地工作，在軍事和民用領(lǐng)域具有重要的應(yīng)用價(jià)值[1]。在ISAR 系統(tǒng)中，距離向的高分辨是通過大時(shí)間帶寬積的波形實(shí)現(xiàn)的，方位向的高分辨一般通過長相干積累時(shí)間內(nèi)雷達(dá)觀測目標(biāo)的轉(zhuǎn)角而獲得。然而，在現(xiàn)實(shí)世界中，觀察環(huán)境可能很復(fù)雜?？赡軙袕?qiáng)噪聲源污染雷達(dá)回波，導(dǎo)致信噪比SNR 低。也可能很難獲得連續(xù)較長的積累時(shí)間，導(dǎo)致方位分辨率較低，同時(shí)造成脈沖積累后的信噪比相對較低。當(dāng)ISAR系統(tǒng)的方位分辨率較低時(shí)，將嚴(yán)重降低ISAR 圖像質(zhì)量，對目標(biāo)的分類與識別精度造成非常不利的影響，因此研究ISAR 系統(tǒng)的超分辨率成像算法是非常必要的。

目前，研究ISAR 超分辨率成像的思路主要基于壓縮感知CS（Compressive Sensing）理論[2]。基于CS的成像算法以其結(jié)構(gòu)簡單、圖像質(zhì)量高、計(jì)算復(fù)雜度低的優(yōu)勢，被認(rèn)為是一種優(yōu)秀的數(shù)據(jù)處理算法。如果信號是“稀疏的”，它就可以從很少的觀測數(shù)據(jù)中被精確地重構(gòu)出來。CS的框架主要分為三個部分：測量矩陣、稀疏表示、重構(gòu)算法。其中重構(gòu)算法的選取對重構(gòu)精度和重構(gòu)速度都有很大的影響。常用的重構(gòu)算法包括凸優(yōu)化算法、貪婪算法和統(tǒng)計(jì)學(xué)習(xí)算法等幾大類。目前常用的線性LBI 算法[3]、OMP算法[4]、SBL 算法[5]分別是上述三類算法的典型代表。然而，上述算法在低信噪比的情況下都存在著重構(gòu)性能不佳、計(jì)算復(fù)雜度高等不足。

近年來，快速分裂Bregman 迭代算法FSBI（Fast Split Bregman Iteration）因其抗噪性能強(qiáng)的優(yōu)點(diǎn)，被廣泛用于解決大規(guī)模優(yōu)化問題，在壓縮感知領(lǐng)域取得了良好的效果[6]。受此啟發(fā)，為實(shí)現(xiàn)低信噪比下的ISAR 快速超分辨成像，本文提出了一種基于快速分裂Bregman 迭代的ISAR 超分辨成像算法。首先，在正則化框架下引入L1范數(shù)表示目標(biāo)的稀疏先驗(yàn)信息，將方位分辨率提高問題轉(zhuǎn)化為正則化問題；其次，為實(shí)現(xiàn)快速重構(gòu)，將托普利茲矩陣的低位移秩特征[7]和Gohberg-Semencul 表示[8]相結(jié)合，實(shí)現(xiàn)了矩陣的快速求逆。通過上述兩方面的措施，既能利用分裂Bregman 迭代算法在低信噪比條件下的重構(gòu)能力[9]，又能保證快速成像。最后，通過仿真和實(shí)測數(shù)據(jù)處理驗(yàn)證了本文提出算法的性能。

1 壓縮感知ISAR 成像模型的建立

經(jīng)過脈沖壓縮和平動補(bǔ)償后的某一距離單元的回波可以表示為：

式中，K是散射點(diǎn)的數(shù)量，Ak為第k個散射點(diǎn)的回波幅度，tm為慢時(shí)間，MT為方位向時(shí)間長度，M為脈沖個數(shù)，fk為第k個散射點(diǎn)相對于參考點(diǎn)的多普勒頻率，n(tm)表示加性噪聲。

為實(shí)現(xiàn)方位向的超分辨率成像，分別定義多普勒采樣間隔和時(shí)間采樣間隔為Δfd和Δt，脈沖重復(fù)頻率為fr，M=TM/Δt。離散化慢時(shí)間序列表示為tm=[1:M]T·Δt，Δt=1/fr，多普勒域?yàn)閇-fr/2,fr/2]，各多普勒單元的散射點(diǎn)的強(qiáng)度用向量σ表示為σ=，σq（q=0,1,···,Q-1）表示位于第q個多普勒單元內(nèi)的散射系數(shù)，=qΔfd-fr/2。只有當(dāng)多普勒單元內(nèi)有散射點(diǎn)時(shí)，σq才不為0。由于包含散射點(diǎn)的多普勒單元的數(shù)量很少，σ可以視為稀疏的。

式中，s為方位向信號向量，n為噪聲向量。

在壓縮感知理論中，對稀疏信號s的觀測不是直接測量信號s本身，而是通過非相關(guān)測量將信號s投影到一組低維的觀測向量，即有：

式中，y為觀測向量，Φ為測量矩陣，M′為觀測向量中元素的個數(shù)。由于M′＜M，直接求解式（3）是一個欠定問題，所以無法直接從y中解出信號s。將式（2）代入式（3）中，有：

式中，=A ΦΨ為感知矩陣，E=Φn為噪聲矩陣?？梢酝ㄟ^求解式（4）得到稀疏系數(shù)向量σ，從而最終重構(gòu)出稀疏信號s。

可以將式（4）轉(zhuǎn)化為如下優(yōu)化問題：

式（5）可以利用傳統(tǒng)的線性規(guī)劃方法進(jìn)行求解，但是噪聲的擾動可能導(dǎo)致估算值遠(yuǎn)離真實(shí)值。因此將式（5）轉(zhuǎn)化為以下正則化形式：

式中，μ是正則化參數(shù)，f(σ) 表示正則化項(xiàng)。

通常，對正則項(xiàng)的選擇取決于目標(biāo)先驗(yàn)信息?？梢砸牒侠淼恼齽t項(xiàng)來提高分辨率。在ISAR 成像場景中，通常只關(guān)注強(qiáng)散射點(diǎn)目標(biāo)。在稀疏假設(shè)下，將式（6）轉(zhuǎn)化為以下形式：

通過正則化項(xiàng)‖σ1‖ 控制的稀疏性，通過-‖Aσ y控制的誤差，同時(shí)保證了解的稀疏性和準(zhǔn)確性，而且能較好地克服噪聲的影響。

由于L1范數(shù)和L2范數(shù)都是凸函數(shù)，可以通過凸優(yōu)化方法對式（7）求解。然而L1范數(shù)是不可微的，所以求解它具有挑戰(zhàn)性。

2 快速分裂Bregman 迭代ISAR 超分辨率成像算法

在本節(jié)中，利用托普利茲矩陣的低位移秩特征和Gohberg-Semencul 定理對分裂Bregman 迭代SBI（Split Bregman Iteration）算法進(jìn)行加速，以實(shí)現(xiàn)L1正則化問題的快速求解。

2.1 分裂Bregman 迭代算法

為了求解式（7），首先引入輔助變量g來代替σ，式（7）變?yōu)椋?/p>

將上面的約束問題進(jìn)一步轉(zhuǎn)化為下面的非約束優(yōu)化問題：

式中，γ是大于0的常數(shù)。表面上，式（8）和式（9）看起來是多余的，然而，它們能極大地降低計(jì)算復(fù)雜度并獲得更有效的迭代策略。

另一個重要的概念是Bregman 距離。基于凸函數(shù)J(·)上的點(diǎn)u和w的Bregman 距離定義為[10]：

式中，向量? J(u) 是J(u)的梯度。

式（11）中的最小化問題可以轉(zhuǎn)化為對2 個子問題σk+1和gk+1分別進(jìn)行優(yōu)化，即有：

對子問題σk+1的求解可以采用多種方法。當(dāng)問題規(guī)模不大時(shí)，對式（12）第1 個式子中的σ進(jìn)行求導(dǎo)并令導(dǎo)數(shù)等于0，可以得到：

式（13）中涉及到求逆運(yùn)算，計(jì)算復(fù)雜度較高。為簡化計(jì)算，可以使用Gauss Seidel 迭代法來求解[11]。子問題gk+1可直接利用軟閾值法得到最優(yōu)解，其表達(dá)式如下：

其中：

式中，sign()表示符號函數(shù)。在求解式（14）的過程中，式（12）中的第2 個式子可直接對b進(jìn)行更新。

2.2 快速迭代計(jì)算算法

雖然用SBI 求解L1正則化問題需要三個步驟，但主要的計(jì)算復(fù)雜度通常來自于σ的求解，即式（13）。g和b的求解，只涉及到一些簡單的加減運(yùn)算，與式（13）相比，它們的計(jì)算復(fù)雜性可以忽略不計(jì)。因此，這里主要分析式（13）的計(jì)算復(fù)雜度。經(jīng)分析后可知，其計(jì)算復(fù)雜度是O(N3)。

實(shí)際上，可以通過降低矩陣求逆引起的高計(jì)算復(fù)雜度來降低式（13）的計(jì)算復(fù)雜度。首先，將式（13）改寫為：

其中，F(xiàn)=μATA+λI，dk=μATy+γ（g k-bk）。從矩陣A和I的結(jié)構(gòu)來看，可以看出矩陣F具有托普利茲結(jié)構(gòu)。因此，F(xiàn)-1計(jì)算復(fù)雜度的降低可以通過使用合適的Gohberg-Semencul 表示來解決。并且，式（16）的計(jì)算可以通過使用快速托普利茲向量乘法更有效地實(shí)現(xiàn)。

基于Gohberg-Semencul 表示的加速策略可以實(shí)現(xiàn)系數(shù)矩陣的快速求逆，大大降低了矩陣求逆的計(jì)算復(fù)雜度。

首先，自回歸系數(shù)a=（1,a2,···,aM）T和預(yù)測誤差e可以通過Levinson-Durbin 算法從Yule-Walker AR方程獲得[12]，即有：

基于Gohberg-Semencul 表示，托普利茲矩陣的求逆可以表示為[13]：

其中：

因此，式(16)的計(jì)算可以轉(zhuǎn)化為：

由于U和V的特殊結(jié)構(gòu)，UUHdk和VVHdk的計(jì)算可以通過快速傅里葉變換FFT（Fast Fourier Transform）和矩陣截取來實(shí)現(xiàn)。

為了實(shí)現(xiàn)快速計(jì)算，構(gòu)造如下矩陣：

U可以通過截取U1的前M行獲得，UH可以通過截取U2的第M行到第2M-1 行獲得。因此，U與向量的乘積可視為矩陣U1與向量的乘積結(jié)果的第1 到第M行，即U和向量FFT的結(jié)果的第1 到第M個元素，UH與向量相乘也是同理。同樣地，VVHdk也可以通過兩次FFT 和矩陣截取來計(jì)算得到。

2.3 計(jì)算復(fù)雜度分析

本文FSBI 算法的計(jì)算復(fù)雜度分析如下。假設(shè)迭代次數(shù)為K，首先，利用Levinson-Durbin 算法計(jì)算自回歸系數(shù)a和預(yù)測誤差e，計(jì)算復(fù)雜度為O((N-1)2)。式（25）的解可以通過四個托普利茲向量運(yùn)算來實(shí)現(xiàn)，計(jì)算復(fù)雜度為7φ(2N)+φ(N)+O(4N)，dk的計(jì)算復(fù)雜度為φ(N)+O(4N)。因此，經(jīng)過加速后的算法的計(jì)算復(fù)雜度為O(K((N-1)2+4N))+7Kφ(2N)+φ(N)。而SBI算法的計(jì)算復(fù)雜度為O((K+1)N3+5KN2+3KN+NlogN)+φ(N)。計(jì)算復(fù)雜度從O(N3)降低到O(N2)，計(jì)算復(fù)雜度大大降低。

2.4 正則化參數(shù)的選取

通常采用L-curve曲線法來確定正則化參數(shù)[14]。圖1給出利用L-curve曲線確定正則化參數(shù)的示意圖，其中橫坐標(biāo)表示正則化項(xiàng)殘差的范數(shù)，縱坐標(biāo)表示解的范數(shù)，曲線上的拐點(diǎn)代表了最佳正則化參數(shù)的位置。

在計(jì)算式（12）的過程中，需要確定μ和γ兩個正則化參數(shù)。在第一次迭代中，參數(shù)γ的初值是根據(jù)經(jīng)驗(yàn)確定的，然后通過L-curve 曲線法獲得正則化參數(shù)μ的初值。具體而言，是把作為的函數(shù)并作曲線，曲線上的點(diǎn)與正則化參數(shù)μ一一對應(yīng)。通過曲線的拐點(diǎn)處的值獲得對應(yīng)μ的值。同理，參數(shù)γ可以通過式（12）的第二個式子確定，L-curve 曲線法繪制作為的函數(shù)，通過曲線的拐點(diǎn)處的值獲得對應(yīng)μ的值

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

首先，針對一維信號，采用蒙特卡洛模擬的方法，評估并比較了OMP 算法、SBI（Split Bregman Iteration）算法、LBI 算法、本文的FSBI 算法在不同信噪比條件下的稀疏信號重構(gòu)性能；接下來，針對短積累時(shí)間和數(shù)據(jù)隨機(jī)缺失兩種情況，利用仿真的方法生成了低信噪比條件下勻速轉(zhuǎn)動船只的ISAR 模擬回波。評估并比較了RD 算法、OMP 算法、SBI 算法、本文FSBI 算法的成像結(jié)果與圖像性能指標(biāo)；最后，針對實(shí)測的船只ISAR 回波數(shù)據(jù)進(jìn)行了超分辨率成像處理，評估并比較了OMP 算法、SBI 算法、LBI 算法、FSBI 算法的成像結(jié)果與圖像性能指標(biāo)。對于實(shí)測船只ISAR 回波數(shù)據(jù)的處理，同樣分為了短積累時(shí)間和數(shù)據(jù)隨機(jī)缺失兩種情況。

3.1 一維稀疏信號重構(gòu)的性能分析

假設(shè)原始信號x是一個1 024×1的隨機(jī)稀疏向量，其中有50 個非零值。經(jīng)乘以測量矩陣后降維為300×1的觀測向量，通過添加高斯白噪聲使信噪比從-5 dB 到20 dB 變化。對于每種信噪比水平，通過200次的蒙特卡洛模擬實(shí)驗(yàn)，計(jì)算重構(gòu)值與真實(shí)值之間的均方根相對誤差RMSRE（Root Mean Square Relative Error）。這里的RMSRE 定義為：

式中，i為元素序號，N為向量中的元素個數(shù)，xri為重構(gòu)值，xti為真實(shí)值。

圖2 給出了OMP 算法、SBI 算法、LBI 算法、本文FSBI 算法的RMSRE 隨信噪比水平的變化。需說明的是，由于FSBI 算法和SBI 算法的原理相同，二者的RMSRE 曲線相同。從圖中可以看出，本文FSBI 算法和SBI 算法在各種信噪比水平下的RMSRE值均最小。

表1 給出了四種算法的平均運(yùn)行時(shí)間，其中LBI 算法的運(yùn)行時(shí)間最長，接下來依次是SBI 算法、FSBI 算法和OMP 算法。因此，與現(xiàn)有常用算法相比，本文FSBI 算法可以在可接受的運(yùn)行時(shí)間內(nèi)獲得最佳的稀疏重建性能。

表1 不同算法平均運(yùn)行時(shí)間的比較Table 1 Comparison of the average running time for the four algorithms

3.2 仿真數(shù)據(jù)的實(shí)驗(yàn)結(jié)果

針對短積累時(shí)間和數(shù)據(jù)隨機(jī)缺失兩種情況，利用仿真的方法生成了低信噪比條件下勻速轉(zhuǎn)動船只的ISAR 模擬回波。對OMP 算法、SBI 算法、LBI 算法和FSBI 算法的成像結(jié)果和成像性能指標(biāo)進(jìn)行了評估比較。仿真中使用的主要參數(shù)如表2 所示。仿真過程中使用的船只散射點(diǎn)模型如圖3 所示。針對短時(shí)間積累的情況，通過加入高斯白噪聲，使得距離壓縮后的回波信噪比為-5 dB 和0 dB。針對數(shù)據(jù)隨機(jī)缺失的情況，數(shù)據(jù)缺失比例為25%。

表2 仿真過程中使用的主要參數(shù)Table 2 Main parameters used in simulation

首先，測試了短積累時(shí)間情況下利用幾種不同算法獲得的ISAR 圖像。圖4 和圖5 顯示了該情況下的ISAR 回波和利用幾種不同算法獲得的ISAR 圖像。其中，圖4（a）為模擬得到的ISAR 距離壓縮回波，圖4（b）、（c）、（d）、（e）分別為利用RD 算法、OMP 算法、LBI 算法、本文FSBI 算法獲得的ISAR 成像結(jié)果。因?yàn)镕SBI 算法和SBI 算法的原理相同，獲得的ISAR 圖像相同。它們之間的唯一區(qū)別是計(jì)算復(fù)雜性。因此，圖中結(jié)果僅列出了FSBI 算法的成像結(jié)果。從圖中可以觀察到，OMP、LBI、SBI、FSBI 算法等幾種基于壓縮感知的成像算法獲得了比RD 算法更好的聚焦圖像。對比幾種壓縮感知算法獲得的圖像發(fā)現(xiàn)，OMP 算法產(chǎn)生的圖像虛假點(diǎn)相對較多，LBI 算法產(chǎn)生的圖像較暗淡，而FSBI 算法生成的ISAR 圖像獲得了最清晰的船只結(jié)構(gòu)。為了進(jìn)一步定量比較四種算法成像結(jié)果的質(zhì)量和計(jì)算復(fù)雜度，計(jì)算了圖像熵并記錄了運(yùn)行時(shí)間。表3 列出了四種算法的圖像熵和運(yùn)行時(shí)間。從表中可以看出，在SNR=-5 dB的條件下，F(xiàn)SBI 算法和SBI 算法的圖像熵最小，在SNR=0 dB的條件下，LBI 算法的熵值最小，F(xiàn)SBI 算法與LBI 算法的熵值接近。與此同時(shí)，F(xiàn)SBI 算法的運(yùn)行時(shí)間也相對較短。因而取得了最佳的成像質(zhì)量。

表3 短積累時(shí)間情況下不同算法成像結(jié)果的圖像熵和運(yùn)行時(shí)間比較Table 3 Comparison of image entropy and running time of four algorithms in the case of short accumulation time

接下來，測試了數(shù)據(jù)隨機(jī)缺失情況下利用幾種不同算法獲得的ISAR 圖像，設(shè)置數(shù)據(jù)缺失比例為25%。圖6 顯示了該情況下的ISAR 距離壓縮回波和利用幾種不同算法獲得的ISAR 圖像。其中，圖6（a）為模擬得到的ISAR 距離壓縮回波，在這種情況下，由于采樣網(wǎng)格不連續(xù)，RD 算法失效。圖6（b）、（c）、（d）分別為利用OMP 算法、LBI 算法、FSBI 算法獲得的ISAR 成像結(jié)果。同樣計(jì)算了幾幅圖像的圖像熵和運(yùn)行時(shí)間。表4 列出了對應(yīng)的計(jì)算結(jié)果。對比計(jì)算結(jié)果可發(fā)現(xiàn)，F(xiàn)SBI 算法的圖像熵值最低、運(yùn)行時(shí)間相對較短。綜合而言，F(xiàn)SBI 算法的成像性能最優(yōu)。

表4 數(shù)據(jù)隨機(jī)缺失情況下不同算法成像結(jié)果的圖像熵和運(yùn)行時(shí)間比較Table 4 Comparison of image entropy and running time of four algorithms in the case of random loss of data

3.3 實(shí)測數(shù)據(jù)的實(shí)驗(yàn)結(jié)果

本節(jié)利用艦船ISAR 實(shí)測數(shù)據(jù)來進(jìn)一步驗(yàn)證本文FSBI 算法的有效性。雷達(dá)信號的中心頻率為9.25 GHz，帶寬為500 MHz，脈沖重復(fù)頻率為100 Hz。完整的雷達(dá)回波包含256 個脈沖，每個脈沖由256 個采樣值組成。針對短時(shí)間積累的情況，通過加入高斯白噪聲，使得距離壓縮后的回波信噪比為0 dB和5 dB。針對數(shù)據(jù)隨機(jī)缺失的情況，數(shù)據(jù)缺失比例為20%。

圖7 和圖8 顯示了短時(shí)間積累情況下的實(shí)測數(shù)據(jù)的ISAR 回波和利用幾種不同算法獲得的ISAR 圖像。其中，兩圖中的子圖（a）為實(shí)測數(shù)據(jù)的ISAR 回波，兩圖中的子圖（c）、（d）、（e）分別為利用OMP 算法、LBI 算法、FSBI 算法獲得的ISAR 成像結(jié)果。計(jì)算了幾幅圖像的圖像熵并記錄了運(yùn)行時(shí)間。表5 列出了相應(yīng)的圖像熵和運(yùn)行時(shí)間。對比計(jì)算結(jié)果可發(fā)現(xiàn)，F(xiàn)SBI 算法的圖像熵值相對較低，因而成像質(zhì)量最優(yōu)且運(yùn)行時(shí)間相對較短。

表5 短時(shí)間積累情況下不同算法對實(shí)測數(shù)據(jù)成像結(jié)果的指標(biāo)值Table 5 Index values of imaging results of the real data corresponding to four algorithms in the case of short accumulation time

圖9 顯示了數(shù)據(jù)隨機(jī)缺失情況下利用幾種不同算法獲得的ISAR 圖像。數(shù)據(jù)缺失比例為20%。其中，圖9（a）為ISAR 距離壓縮回波，圖9（b）、（c）、（d）分別為利用OMP 方法算法、LBI 算法、FSBI算法獲得的ISAR 成像結(jié)果。同樣計(jì)算了幾幅圖像的圖像熵和運(yùn)行時(shí)間。表6 列出了對應(yīng)的計(jì)算結(jié)果。對比計(jì)算結(jié)果可發(fā)現(xiàn)，所提本文方法算法（FSBI）的圖像熵值最低、運(yùn)行時(shí)間相對較短。因而，綜合而言，F(xiàn)SBI 算法的成像質(zhì)量性能最優(yōu)。

表6 數(shù)據(jù)隨機(jī)缺失情況下不同算法對實(shí)測數(shù)據(jù)成像結(jié)果的指標(biāo)值Table 6 Index values of imaging results of the real data corresponding to four algorithms in the case of random loss of data

4 結(jié)束語

針對低信噪比下的ISAR 超分辨率成像問題，提出了一種基于快速分裂Bregman 迭代的ISAR 超分辨率成像算法。所提出的FSBI 算法的主要優(yōu)勢為：充分利用托普利茲矩陣的低位移秩特征和Gohberg-Semencul 表示實(shí)現(xiàn)了矩陣的快速求逆，從而既能利用SBI 算法在低信噪比條件下的重構(gòu)能力，又實(shí)現(xiàn)了快速成像。通過對仿真和實(shí)測數(shù)據(jù)進(jìn)行處理，驗(yàn)證了所提出的FSBI 算法相比于現(xiàn)有常用算法能夠在低信噪比的情況下具有更佳的聚焦性能，且運(yùn)行時(shí)間相對較短。