交替迭代最小化稀疏穿墻成像快速算法

(1. 桂林電子科技大學信息與通信學院，廣西桂林 541004；2. 廣西無線寬帶通信與信號處理重點實驗室，廣西桂林 541004)

0 引 言

超寬帶穿墻雷達能夠為墻后目標檢測提供良好的手段，因而在災后救援、反恐巷戰等領域得到了廣泛應用[1-3]，經典的成像算法如后向投影(Back Projection，BP)算法、延時求和算法等需要大量的空間與時間數據采樣，耗時且無法避免成像結果中的雜波干擾。為了減少成像數據量并提高成像分辨率，基于壓縮感知的成像技術被引入到穿墻成像中，通過在成像模型中引入稀疏的正則化約束，可將目標的重構問題轉化為1范數優化問題，采用包括貪婪法[4-5]、貝葉斯方法[6-7]、正則化方法[8]等來恢復稀疏信號。文獻[4]和文獻[5]都是根據目標特性構建相應稀疏特性成像字典矩陣，然后通過正交匹配追蹤(Orthogonal Matching Pursuit, OMP)算法完成目標成像，這類算法最重要的一步是從成像場景的最小二乘估計結果中選取支撐集，顯然這一步是最為耗時，且OMP算法極易受到信號輸入信噪比影響。相比于OMP方法，貝葉斯方法能夠在低信噪比情況下取得較好成像，文獻[6]在BCS的基礎上結合擴展目標像素間的結構性信息，然后利用吉布斯采樣方法得到高分辨的成像結果，但該方法加大了成像模型的復雜度致使算法運算時間較長。而文獻[7]通過構建點目標與擴展目標的組合字典，并利用一階泰勒級數展開解決目標網格點偏移問題，算法成像結果具有較高分辨率，然而構建兩個成像字典勢必帶來更多的成像時間。文獻[8]提出一種基于迭代的拉格朗日乘子的穿墻擴展目標成像方法，通過稀疏增強、區域增強以及圖像域優化三步，獲得了高分辨的擴展目標結果，同時拉格朗日乘子項也避開了對先驗信息的獲取，然而該方法同樣存在算法耗時過高的問題。顯然，上述所列的稀疏成像方法均在一定程度上改善了成像結果的分辨率，但他們都不可避免地提高了成像時間，因此，如何保證成像分辨率較高的前提下縮短成像時間對于穿墻應用十分重要。

有鑒于此，本文提出一種基于交替迭代框架的墻后目標成像快速算法，首先利用貝葉斯框架建構起穿墻信號的稀疏表示模型，并利用最大后驗估計準則得到包含參數與目標散射矢量的目標函數，然后在優化最小化(Majorization-Minimization，MM)框架下構造此時目標函數對應的穩健優化函數，最后對聯合概率密度函數的目標散射系數、噪聲功率以超參數依次進行交替迭代優化。與現有的成像方法對比，本文方法在信噪比較低時能實現墻后目標清晰成像的同時，也大幅提升了成像速度。

1 信號模型

使用收發天線共置的雷達系統沿著平行于前墻的方向對墻后目標進行J點的合成孔徑探測，同時將感興趣場景(SOI)在距離向和方位向上劃分為N=Nx×Nz個像素網格，定義xn表示第n像素點的目標反射系數，n=1,2,…,N，則第j個探測位置接收到的目標回波信號表示為

j=1,2,…,J

(1)

對式(1)的sj(t)進行K點采樣，采樣后的數據可以表示成

ej(kTs)，k=1,2,…,K

(2)

式中，Ts表示采樣周期。將式(2)寫成矢量形式表示為

rj=PjDjx+ej

(3)

式中，rj=[rj(0),rj(1),…,rj(K-1)]T表示第j個測量位置處的回波信號矢量，x=[x1,x2,…,xN]T表示散射系數矢量的值，Dj是由距離衰減系數組成的對角矩陣，表示為Dj=diag(αj1,αj2,…,αjN)，Pj為一個K×N的矩陣，其中元素表示為脈沖信號相應平移，具體為

(4)

考慮到所有J個孔徑接收到的回波數據，可以得到以下線性關系:

r=Ψx+e

(5)

y=Φ(Ψx+e)=Αx+n

(6)

式中，Φ表示一個維度是M×JK的測量矩陣，矩陣Α=ΦΨ表示從Ψ中隨機抽取行和列構成新的字典。

2 優化最小化框架的遞歸迭代成像算法

2.1 算法描述

假設穿墻回波信號y滿足貝葉斯模型

p(y|x,η)=N(Ax,ηI)

(7)

式中，η為高斯白噪聲的方差，且該數值即為噪聲的功率，其取值范圍η∈[0,+)，因而可以合理假設其先驗概率密度函數為

p(η)∝1

(8)

另外，假設目標散射系數x中的元素為獨立同分布，并服從未知參數為λ的拉普拉斯分布，所以散射系數矢量的概率密度函數表示為

(9)

本文利用最大后驗估計(MAP)準則來重構目標散射系數矢量，先給出目標散射系數的聯合后驗概率分布函數p(x|y,λ,η)：

p(x|y,λ,η)∝p(y|x,η)p(x|λ)p(η)=

(10)

根據MAP準則，需要通過最大化聯合后驗概率密度函數p(x|y,λ,η)，進行實現對多變量目標函數的優化求解。為方便描述，將式(10)取對數后表示為

(11)

通過最大化目標函數ln(p(x|y,λ,η))，同時忽略ln(p(x|y,λ,η))中的常數項，可以求解未知的散射系數矢量x、噪聲功率η以及超參數λ的MAP估計值，即

(12)

在尋找散射系數矢量x最優解過程中，隨著迭代次數增加，目標函數G(x,λ,η)是一個遞減的趨勢[13]，第l+1次迭代得到的x(l+1)均比前一次的x(l)要小，因此，目標函數G(x,λ,η)存在以下不等式成立：

(13)

(14)

(15)

([Ax]m)2≤qk(x,x(l))

(16)

(18)

結合式(17)和式(18)，G(x(l),λ(l),η(l))在當前散射系數估計x(l)下的最優函數可表示為

(19)

且Q(x,λ(l),η(l)；x(l))滿足

G(x,λ(l),η(l))≤Q(x,λ(l),η(l);x(l))

l=1,2,…,L

(20)

充分借助大數據等計算技術，基于課程關鍵字、學員信息等數據的分析，為學員提供個性化資源推薦。依靠技術記錄學員學習時間、終端系統等，形成用戶行為數據，通過自動匯總分析系統將結果反饋給學員及干部網絡教育決策人員。

n=1,2,…,N

(21)

在第2步中，固定散射系數矢量x和超參數λ，對噪聲方差η進行優化，此時的目標函數表示為

(22)

將G(η,λ(l),x(l))對η求偏導，并令其為0，可得

(23)

最后，固定噪聲方差η和散射系數矢量x去估計超參數λ，此時的目標函數G(λ,x(l),η(l))表示為

(24)

計算G(λ,x(l),η(l))對λ的偏導數，并令其為0，可得

(25)

考慮到l1范數在零點處不可導，此處可以對向量范數進行平滑近似[12]，

(26)

此時的λ(l+1)可以改寫為

(27)

式中，ε表示平滑因子，其值采用逐漸遞減參數的方法，即ε(l)=ε(l-1)/10，直到滿足迭代終止條件為止。

2.2 算法實現

本文算法在估計散射系數矢量x的步驟中，利用優化最小化框架避免了直接利用封閉形式求解帶來的求逆運算，有效提高了算法運算效率。算法流程如表1所示。在構建字典矩陣、測量矩陣、測量值向量和初始參數后，進入式(21)、式(23)和式(27)的循環迭代過程，直到滿足條件后退出循環。

表1 算法流程

3 仿真與實驗結果與分析

3.1 仿真結果分析

利用電磁仿真軟件GprMax仿真墻后目標成像得到仿真數據。前墻體長為3 m，墻體的厚度為0.2 m，相對介電常數為6.4，電導率為0.01 S/m。由31個收發共置的天線單元構成線性陣列均勻分布在橫軸0.9～2.1 m、縱軸0.5 m處，陣元間隔0.04 m，線陣距離墻體0.5 m。發射脈沖寬度為1 ns高斯脈沖信號，GprMax的網格單元為0.01 m，時間步長為23 ps，采樣時間窗為25 ns，仿真中成像區域的方位向與距離向網格大小均設置為0.05 m，為模擬真實情況，仿真中引入高斯白噪聲。

圖1給出了GprMax仿真模型，為模擬不同形態的墻后目標成像場景，設置混合目標并存的仿真場景，圖1(a)為兩個目標仿真示意圖，圖中擴展目標長寬分別為0.3 m和0.2 m，且兩個目標彼此之間以及與前墻體的距離分別為0.6 m和1.3 m，圖1(b)為點目標與擴展目標混合仿真示意圖，擴展目標尺寸和距離與圖1(a)一致，點目標分別位于(1.05,2)和(1.35,2.25)處。

(a) 兩目標

(b) 混合目標圖1 GprMax仿真模型

由于采用稀疏成像，本文選擇全部31個天線單元，每個接收天線的數據隨機選擇200個采樣點。圖2和圖3分別給出了兩目標體和混合目標體的成像結果，圖中的虛線框表示目標真實位置。圖2(a)和圖3(a)的BP成像方法，可以看出BP算法得到的旁瓣雜波較多，目標之間雜波干擾導致分辨率急劇下降，并且在目標真實位置的后面產生了嚴重的虛假像，難以識別出目標位置；圖2(b)和圖3(b)采用文獻[4]的成像方法，成像結果相對BP算法分辨率提高，但仍然無法消除目標之間旁瓣影響，目標體之間輪廓特性不明顯；而圖2(c)和圖3(c)是采用本文方法成像結果，相比較前兩種成像結果，本文方法得到的擴展目標輪廓連貫，點目標成像聚焦性較高，而且目標與目標之間的旁瓣影響被抑制，圖像較為清晰，同時本文的優化方法提高了計算速度。

(a) BP成像

(b) 文獻[4]方法

(c) 本文方法圖2 兩目標成像結果

(a) BP成像

(b) 文獻[4]方法

(c) 本文方法圖3 混合目標成像結果

為了對比算法的成像性能，圖4給出了兩種成像場景下目標雜波比(Target-to-Clutter Ratio, TCR)與信噪比之間的關系圖，從圖4(a)可以看出，隨著信噪比提升TCR數值也呈現上升趨勢，且在相同信噪比下，本文算法的TCR值都高于文獻[4]算法和BP算法，說明本文得到的目標像聚焦性能優于這兩種算法；考慮到圖4(b)中成像場景是點目標與擴展目標共存，成像環境較為復雜，因而成像結果雜波較多，TCR數值相較于純擴展目標場景要小，但本文算法的目標雜波比值仍略高于兩種對比算法成像結果，進一步說明本文方法在復雜環境中也能得到較好的聚焦性能。

(a) 兩目標成像

(b) 混合目標成像圖4 不同成像場景性能對比

為了看出本文算法在成像效率上的改進，將本文算法與對比算法進行對比，需要說明的是，由于本文算法進行的是稀疏信號恢復進而成像，所以這里僅是對比稀疏算法之間的效率問題，因此，可以根據稀疏成像所需數據量來對比程序運行時間，利用程序運行時間(Program Running Time, PRT)來對表征各算法效率，其值越小效率越高。表2分別給出了兩種成像場景下的數據量與PRT之間的關系，從表中可以看出，在數據量比例相同時，本文方法回復信號所需時間遠遠小于文獻[4]的成像結果，說明本文算法有效改進了成像效率。

表2 兩種成像場景下的數據量與PRT之間的關系

3.2 實驗結果分析

使用美國GSSI公司的探地雷達SIR-20搭建實驗場景，實驗墻體為實心磚墻，厚度為0.2 m，相對介電常數為6.4，天線距離墻體0.2 m，選用1 GHz的喇叭天線架高1.2 m，從-0.6 m到0.6 m以間距0.04 m水平移動掃描距前墻體1.8 m處的柜子，柜子長寬分別是0.9 m和0.3 m，共計測量31道，每道采樣點數為1 024。將SIR-20系統收集的回波數據經取平均、去噪、雜波相消和自動增益控制等信號處理得到墻體回波，使用全部31個天線單元位置，每個接收天線隨機選擇200個采樣點進行成像。

圖5給出了實驗場景的示意圖，圖6給出了BP算法、文獻[4]方法以及本文方法的實驗場景數據的成像結果，虛線框表示擴展目標的真實位置。可以看出，實驗數據處理中BP方法成像圖像旁瓣雜波較多，目標邊緣輪廓無法識別，整體圖像分辨率較低，聚焦性較差；由于實測數據的信噪比不高，文獻[4]算法方法對此較為敏感，直接導致目標成像結果散焦嚴重，幾乎無法分辨目標像。相比之下，本文方法由于在稀疏信號貝葉斯建模時充分利用分層模型的特點，使信號進一步稀疏，因此目標圖像的邊緣輪廓較為明顯，圖像旁瓣得到有效抑制，圖像分辨率較高，同時該算法還結合MM優化框架，以較快速度形成成像結果。

(a) 實驗設備

(b) 實驗場景圖5 成像場景

(a) BP成像

(b) 文獻[4]方法

(c) 本文方法圖6 實測數據成像結果

4 結束語

本文提出一種基于交替迭代框架的墻后目標成像算法，利用稀疏信號貝葉斯模型的最大后驗估計準則得到包含參數與目標散射矢量的目標函數，然后在優化最小化(MM)框架下利用目標函數對應的優化函數對目標散射系數、噪聲功率和超參數進行交替迭代求解，貝葉斯分層模型保證了成像結果的分辨率，而基于優化最小化準則的交替迭代方式有效改善了該方法的成像效率。仿真數據和實驗數據的處理結果表明，擴展目標輪廓邊緣清晰，點目標成像聚焦度明顯，同時成像速度得以較大提高。在下一步的工作中，將繼續對該方法進行完善使其能應用到更復雜的場景中。