基于IALM的穿墻雷達成像算法

張珊珊，趙建華

西安工業(yè)大學 電子信息工程學院，西安710021

穿墻雷達成像（Though-the-Wall Radar Imaging，TWRI）是一種利用低頻電磁波的穿透性對復雜場景或障礙物后隱藏目標進行有效探測、定位和成像的技術(shù)，在應(yīng)急救援、醫(yī)療監(jiān)護、城市安保和反恐斗爭等領(lǐng)域都發(fā)揮著巨大作用[1]。

穿墻成像中常采用超寬帶（Ultra Wide Band，UWB）信號作為發(fā)射信號來提高距離向分辨率和探測精度[2]，但傳統(tǒng)基于香農(nóng)-奈奎斯特采樣定理（Shannon/Nyquist Sampling Theory）的成像算法均存在著采樣率高、數(shù)據(jù)量大、數(shù)據(jù)傳輸、存儲及后續(xù)數(shù)據(jù)快速處理較為困難等不足[3]。此外，雷達成像近似視為對探測場景進行的最小二乘估計，使得基于匹配濾波思想的傳統(tǒng)成像算法會出現(xiàn)較高旁瓣，使得成像結(jié)果模糊，無法區(qū)分鄰近目標，成像質(zhì)量不高[4]。

近年來，國內(nèi)外學者通過研究與壓縮感知（Compressive Sensing，CS）理論結(jié)合的雷達成像算法，以此來改善傳統(tǒng)算法的不足。美國Rice大學的學者Baraniuk在研究雷達成像時引入CS理論，并在后續(xù)的仿真模擬中，利用正交匹配追蹤（Orthogonal Matching Pursuit，OMP）算法成功重建成像場景，驗證了壓縮感知雷達成像的可行性，這也是研究者第一次通過理論證明了壓縮感知成像算法優(yōu)于其他傳統(tǒng)成像算法[5]。美國的Gurbuz第一次通過仿真模擬實驗驗證了壓縮感知技術(shù)在探地雷達成像的可行性，利用地下目標空間的稀疏性，運用壓縮感知思想，從有限次觀測數(shù)據(jù)中重構(gòu)出原始目標像，驗證了該方法相比于傳統(tǒng)BP方法在成像質(zhì)量、魯棒性和測量帶寬要求上更有優(yōu)勢[6]。我國的余慧敏等學者則側(cè)重于壓縮感知雷達三維成像領(lǐng)域的研究，結(jié)合高斯脈沖體制的探地雷達，采用稀疏約束的正則化算法對隨機孔徑位置壓縮測量的回波數(shù)據(jù)進行恢復[7]。2010年，壓縮感知技術(shù)被Huang等人應(yīng)用到穿墻雷達成像中[8]，提出了利用壓縮感知算法隨機采樣返回到各個天線位置的所有回波數(shù)據(jù)，再通過重構(gòu)算法獲得圖像。維拉諾瓦大學的Yoon等在研究穿墻雷達成像時提出了隨機選擇若干個天線位置上相同的頻率，再采樣進行壓縮感知，重建各個天線位置的數(shù)據(jù)信息，最后進行各個部分的拼接成像[9]。

基于CS理論的雷達成像系統(tǒng)可采用低速的A/D轉(zhuǎn)換器進行采樣，極大程度地降低了雷達接收端對采樣速率的要求，簡化了硬件設(shè)計，同時基于CS理論的成像算法通過引入稀疏的正則化約束，直接建立成像場景與測量數(shù)據(jù)之間的線性關(guān)系，使得基于CS理論的成像算法往往成像效果更好，但在重構(gòu)過程中須嚴格遵循非線性的稀疏約束，對信噪比要求較高，計算過程較復雜，且在TWRI實際應(yīng)用中，墻體反射回波能量遠遠大于目標回波能量，導致目標信號常被雜波信號所覆蓋，這使得目標的檢測變得十分困難。因此，設(shè)計一種同時兼顧抑制墻雜波和目標成像的算法顯得尤為重要。

本文針對傳統(tǒng)CS成像技術(shù)對噪聲敏感、計算過程復雜的問題，提出一種CS框架下的基于非精確增廣拉格朗日乘子（Inexact Augmented Lagrange Multiplier，IALM）的穿墻成像算法。將TWRI問題視為正則化最小二乘優(yōu)化問題，利用成像場景的稀疏性，分析回波信號數(shù)據(jù)矩陣，引入IALM進行迭代求解，先建立步進頻穿墻雷達信號模型，構(gòu)造目標回波信號稀疏表示字典；再利用隨機測量矩陣實現(xiàn)降采樣，最后將目標成像和墻雜波抑制視為一個含核范數(shù)和l1范數(shù)的復合優(yōu)化問題，減少迭代次數(shù)，保證成像實時性。

1 基于非精確增廣拉格朗日乘子的穿墻雷達成像

1.1 CS基本原理

設(shè)x是一個長度為N、稀疏度為K的離散信號，則x可用一組基ΨT線性組合表示：

其中，α為x在Ψ上的稀疏系數(shù)矩陣，Ψ為x的稀疏變換基或字典矩陣。

構(gòu)造一個合適的采樣矩陣，將高維信號x投影到低維空間上進行降維采樣，則通過降維采樣矩陣Φ=[φ1,φ2,…,φm,…,φM]后的采樣信號可表示為：

其中，y是采樣后的M×1維的觀測信號，Φ是M×N維的測量矩陣，x是N×1維的離散型原始信號矩陣，M為采樣的維度，M?N。則可推出：

其中，A是M×N維傳感矩陣或感知矩陣。重構(gòu)準確性和穩(wěn)定性決定于A是否滿足約束等容特性（Restricted Isomety Property，RIP）[10]條件。

信號重構(gòu)可通過對式（3）求逆運算得系數(shù)α，再恢復原始信號x。由于Φ的列數(shù)遠大于行數(shù)，使得未知數(shù)個數(shù)遠大于求解方程個數(shù)，求解式（3）等同于求l0范數(shù)約束的優(yōu)化問題，解的個數(shù)不唯一。若放寬約束條件，將非凸的l0范數(shù)松弛為凸的l1范數(shù)，即轉(zhuǎn)化為具有唯一解的凸優(yōu)化問題[11]：

其中，‖‖1表示l1范數(shù)。

1.2 信號字典矩陣的構(gòu)造

系統(tǒng)利用合成孔徑技術(shù)（SAR）的思想，利用單站雷達的運動產(chǎn)生與目標之間的相對位移，把小的實天線孔徑合成為一個等效的大虛擬孔徑，提高方位向分辨率。如圖1所示為穿墻雷達系統(tǒng)探測示意圖。

圖1 穿墻雷達系統(tǒng)探測示意圖

依據(jù)探測模型，雷達系統(tǒng)采用的是超寬帶步進頻信號體制，設(shè)采樣測量點有M個，且每個測量點會發(fā)射N個頻率步進的雷達脈沖，則第n(n=0,1,2,…,N-1)個發(fā)射信號頻率為：

其中，f0為初始頻率，Δf為頻率步進量。

天線沿平行于墻面的水平界面上勻速移動，則第m(m=0,1,2,…,M-1)個位置第n個頻點的目標回波可表示為：

其中，P表示探測區(qū)域內(nèi)散射目總數(shù)，σp表示目標點p的復散射系數(shù)，τm,p表示第m個采樣點接收的第p個目標的雙程時延。

雙程時延可由基于幾何光學的穿墻信號模型[12]求得。如圖2所示，以天線所在的陣列為橫軸，垂直于墻面左側(cè)所在的直線為縱軸建立直角坐標系，利用快速墻體補償法求墻體表面折射點進而確定回波時延。

圖2 基于幾何光學的穿墻信號模型

由ΔBA1A2與ΔBxm A3相似，ΔAB1B2與ΔAOB3相似，可得：

電磁波在穿墻傳播時速度會衰減為：

其中，v表示電磁波在墻體內(nèi)部的傳播速度，εw表示墻體相對介電常數(shù)。

由穿透墻體時的折射點可確定傳播路徑，進而求解單程反射時延τ：

當天線沿水平方向位移時，天線到墻體的距離保持不變，因此，墻體回波與傳播時延為一固定值，墻體反射回波與天線的測量位置無關(guān)，可認為在空域上基本保持不變[13]，則墻體反射回波表示為：

其中，σw表示墻反射系數(shù)，τw表示天線到墻體的傳播時延。

則雷達在第m個位置第n個工作頻點的天線接收信號可表示為：

在TWRI的實際應(yīng)用中，成像區(qū)域內(nèi)目標個數(shù)通常為單個或少數(shù)的，因此，穿墻雷達目標信號符合稀疏特性。將成像區(qū)域近似看作由一個個網(wǎng)格點組成的二維平面[14]。再把成像場景沿著橫軸方向和縱軸方向平均劃分為Q個網(wǎng)格點，每一個網(wǎng)格點代表一個像素點，再縱向拉伸，按列排成一個Q×1維矩陣，Q為網(wǎng)格總數(shù)，如圖3所示為成像區(qū)域網(wǎng)格劃分示意圖。

設(shè)sq表示加權(quán)指標函數(shù)，定義為：

圖3 成像區(qū)域網(wǎng)格劃分示意圖

當某一個網(wǎng)格點與目標位置重疊時，則該網(wǎng)格點對應(yīng)的散射強度值被定義為目標復散射系數(shù)σp，否則為0。由于目標僅占據(jù)某一部分的網(wǎng)格，則由目標的復散射系數(shù)組成的二維反射系數(shù)矩陣同樣具有稀疏性，將其按列排列堆疊成一個新的一維列向量s，則系統(tǒng)在第m個觀測位置處的目標回波信號矩陣可表示為：

其中，Ψm=[ψm(n,q)]是一個由e-j2πfnτmq通過列堆疊而成的N×Q維矩陣，τmq為第m個測量位置到第q個網(wǎng)格的雙程反射時延，可由式（9）求得。

將M個測量位置的所有目標回波信號表示為一個MN×Q的矩陣：

同理，處理各個位置的墻雜波信號，系統(tǒng)接收信號組成的字典矩陣表示為：

其中，Zw、Zt和V分別表示墻雜波矩陣、目標回波矩陣和噪聲矩陣。

通常探測區(qū)域中有用的目標僅占據(jù)較少的空間位置，目標信號具有稀疏性，故Zt為稀疏矩陣，只有少量非零元素；而墻體位置相對固定，墻雜波信號在每個天線位置的強度近似相等，具有較強的相關(guān)性，則可近似認為Zw為低秩矩陣，只有很少的非零奇異值。

1.3 測量矩陣的構(gòu)造

本文選擇一種簡單且易于操作的壓縮測量方法，即隨機選取天線測量位置，同時隨機選取脈沖工作頻率進行降采樣，則從MN×MN的單位陣中隨機選取K行來組成相應(yīng)的測量矩陣。壓縮感知采樣示意圖如圖4所示。

降采樣后的數(shù)據(jù)為：

其中，Φ是測量矩陣，v(Z)表示將矩陣按列堆疊成列向量的運算。

將式（14）、（15）代入式（16）可得：

圖4 壓縮感知采樣示意圖

即轉(zhuǎn)化為求一個低秩矩陣和稀疏矩陣的優(yōu)化問題[15]：

其中，‖‖*表示核范數(shù)，‖‖1表示l1范數(shù)，ε為重構(gòu)允許的最大誤差，λ為低秩矩陣和稀疏矩陣約束平衡的正則化參數(shù)。

針對這樣一個含核范數(shù)和l1范數(shù)的帶約束條件的凸優(yōu)化問題，常采用拉格朗日乘子法構(gòu)造無約束條件的優(yōu)化問題：

其中，Zw可通過奇異值軟閾值法進行求解，s可通過l1范數(shù)最小化技術(shù)進行計算。λw、λs分別為低秩矩陣和稀疏矩陣的正則化參數(shù)。但該方法涉及奇異值分解運算較為復雜，導致閾值迭代速度較慢。

1.4 基于IALM的成像重構(gòu)

1.4.1 增廣拉格朗日乘子法

將一個原始數(shù)據(jù)矩陣D拆分成一個低秩矩陣A和稀疏噪聲矩陣E之和，矩陣恢復模型可通過拉格朗日算法重構(gòu)，可表示為：

其中，λ表示正則化參數(shù)，rank(A)是矩陣的秩函數(shù)，一般?。╩、n為矩陣的維度數(shù)），是高度非線性非凸優(yōu)化的組合優(yōu)化問題，沒有有效解決方法。通?？赏ㄟ^對目標函數(shù)進行凸松弛，將核范數(shù)逼近矩陣的秩范數(shù)，矩陣的l1范數(shù)逼近矩陣的l0范數(shù)，從而轉(zhuǎn)換成凸優(yōu)化問題進行求解：

其中，‖‖*表示核范數(shù)，‖‖1表示l1范數(shù)，λ是調(diào)節(jié)低秩矩陣和稀疏矩陣的正則化參數(shù)。

對于這樣一個優(yōu)化模型，一般可用主成分追蹤法、加速近似梯度法、增廣拉格朗日法[16]及交替方向乘子法等來求解。

可得到增廣拉格朗日函數(shù)形式：

其中，Y表示拉格朗日乘子，表示標準內(nèi)積，μ代表懲罰系數(shù)，‖‖F(xiàn)表示Frobenius范數(shù)。

通過迭代的方法將式（23）的最小化公式代入下式進行求解：

通過多次的交替更新迭代方式可直接求得上式的最優(yōu)化解，即首先固定E和Y不變，通過運算式（24）最小化條件下的A，再確定A和Y不變，通過函數(shù)極小化條件下得到E。通過在每一步的迭代中不斷地求矩陣A和E的精確解，最終恢復原始矩陣的方法稱為精確增廣拉格朗日乘子（Exact Augmented Lagrange Multiplier，EALM）法。

1.4.2 基于IALM的成像重構(gòu)

在實際迭代過程中，不需要求得A和E的精確值，而只需迭代一次就可恢復原始矩陣的局部最優(yōu)解，該方法稱為非精確增廣拉格朗日乘子（Inexact Augmented Lagrange Multiplier，IALM）法。

本文提出采用非精確增廣拉格朗日乘子（IALM）法來求解式（18），先固定矩陣S和Y（Y代表拉格朗日乘子），求一個最小化的Zw；接著固定Zw和Y，求一個最小化的S，循環(huán)迭代就可得到收斂的最優(yōu)解。

首先，定義軟閾值算子Λλ(Σ)為：

其中，λ為一常數(shù)。

矩陣奇異值分解（SVD）定義為：

奇異值軟閾值算子（SVT）定義為：

具體迭代步驟為：

（1）參數(shù)初始化：

設(shè)Z0=v*(Φ?y)，其中v?()表示矩陣化運算，表示把NM×1維向量轉(zhuǎn)化為N×M維矩陣；Φ?為Φ的偽逆矩陣；Y0=0，S0=0，μ0=0，ρ>1，k=0。

（2）更新Zw：

（3）更新S：

（4）對參數(shù)μk進行更新：

令k=k+1，當函數(shù)收斂時退出循環(huán)。

2 仿真結(jié)果分析

為驗證本文所提成像算法的有效性和準確性，利用圖1所建立的穿墻雷達探測模型進行算法的驗證。本文所有仿真實驗均在Windows 10系統(tǒng)下進行，成像部分的仿真均在MATLAB R2014a環(huán)境下進行。仿真模型設(shè)置如下：雷達起始頻率fl設(shè)置為1 GHz，步進間隔為10 MHz，終止頻率fh為3 GHz，則步進頻信號共有201個采樣頻點。合成孔徑長度為2 m，每隔5 cm設(shè)置一個天線，共有41個天線測量位置。設(shè)墻體厚度為0.2 m的混凝土材質(zhì)，相對介電常數(shù)為6，且距離天線0.15 m，墻體散射系數(shù)設(shè)為1，目標散射系數(shù)設(shè)為0.1，并加入10 dB噪聲。設(shè)理想目標位于水平向1 m、距離向2 m處，邊長為0.05 m。設(shè)置成像區(qū)域為3 m×3 m，將成像區(qū)域平均分為60×60個網(wǎng)格。從201×41=8 241個頻域樣本數(shù)據(jù)中隨機抽取出1 648個數(shù)據(jù)（保證降采樣）進行成像。迭代過程中，設(shè)置正則化參數(shù)λ為0.2，懲罰系數(shù)μ為1，ρ為2，1E?6為迭代運算的判斷閾值。重構(gòu)結(jié)果均采用歸一化到[0，1]區(qū)間的圖像強度像素值的絕對值來表示。如圖5所示為仿真成像結(jié)果對比，其中圖5（a）為目標真值圖，圖5（b）為直接CS成像，圖5（c）為直接OMP成像，圖5（d）為直接ISTA成像，圖5（e）為抑制雜波[17]后直接CS成像，圖5（f）為抑制雜波后直接OMP成像，圖5（g）為抑制雜波后直接ISTA成像，圖5（h）為ALM成像，圖5（i）為IALM成像結(jié)果。

從圖5可以看出，在壓縮感知框架下利用降采樣數(shù)據(jù)仍能重構(gòu)目標像，驗證了該理論在穿墻雷達成像中具有的利用降維采樣數(shù)據(jù)恢復原始回波信號的獨特優(yōu)勢。其中，圖5（b）與圖5（c）中清晰可見墻體的位置，而圖5（d）中墻體回波較弱，不易直接從圖中觀察到，對比可以看出，利用直接CS成像的效果最差，成像結(jié)果中有大量噪聲，目標受周圍雜波影響較大，目標像分辨率較低，而直接ISTA成像效果最好，但圖中目標像的模糊散焦現(xiàn)象較嚴重。而抑制雜波處理后提高了目標成像質(zhì)量，其中，分別對比圖5（b）與圖5（e）、圖5（c）與圖5（f）可以看出，抑制雜波后目標成像質(zhì)量更高，說明了很好地消除墻體對目標回波的影響，重構(gòu)場景中雜波較少，目標像較清晰，圖像質(zhì)量得到了很大的提高，對比圖5（d）與圖5（g）可看出，目標成像質(zhì)量基本沒有太大的改變，抑制雜波處理的效果不大，說明了抑制雜波對直接ISTA成像提高成像質(zhì)量影響較小。

而圖5（h）與圖5（i）成像質(zhì)量較好，場景中目標散射點成像較清晰，目標旁瓣小，雜波較少，成像精度更高，表明本文所提出的算法在成像同時抑制了墻雜波，在未提前抑制墻雜波的情況下仍能滿足成像要求，分辨率有了顯著提高，能較清晰地重建目標，有更好的聚焦性能和散射點重構(gòu)效果。但由于墻體的存在，均會使目標位置與實際位置存在一定程度的偏差，造成目標成像位置的偏移。

為了定量評價成像重建性能，引入信雜比（Target to Clutter Ratio，TCR）和運行時間作為衡量標準。其中，TCR值越大，說明目標重構(gòu)誤差越小，成像質(zhì)量越好；時間值越小，說明處理速度更快，更符合實際應(yīng)用中實時性要求。

TCR定義為：

其中，Nt和Nc分別表示目標區(qū)域和雜波區(qū)域的像素值個數(shù)，At和Ac分別表示目標區(qū)域和雜波區(qū)域，Iq表示成像區(qū)域內(nèi)第q個像素值。

反復進行50次實驗求得平均值，統(tǒng)計成像重建性能對比結(jié)果如表1所示。

從表1對比看出，各種算法在抑制雜波后均可提高成像質(zhì)量，但相應(yīng)地也增加了運行時間。其中，直接CS成像由于調(diào)用CVX工具箱，使得運行速度較慢。直接CS成像在未進行墻雜波抑制處理時的TCR最小，說明了成像質(zhì)量最差，抑制雜波后的TCR提高了約22 dB，但對應(yīng)處理速度也變慢，平均成像重構(gòu)時間增加了約2倍，TCR增加了約1倍。直接OMP成像經(jīng)雜波抑制處理后，平均成像TCR提高了約13 dB，重構(gòu)時間也相應(yīng)地增加了約20 s。而雜波抑制處理對直接ISTA成像幾乎沒有太大的改善，TCR在處理前后大致相等。直接OMP與直接ISTA成像在運行時間上相差不大，但ISTA類算法在此場景下恢復目標信號的信雜比更高。

而本文所提算法的TCR在表1幾種處理算法中較高，該算法利用降采樣數(shù)據(jù)重構(gòu)目標像時，運算速度有所提高，成像質(zhì)量也有了很好的改善。其中，ALM成像的TCR和IALM成像的TCR相差不大，但處理時間卻慢了約4 s，說明所提算法在提高成像質(zhì)量方面的優(yōu)勢，且可以更好地滿足實時性的要求。

圖5 仿真成像結(jié)果對比

表1 成像重建性能對比

從仿真成像結(jié)果和重建性能對比中可看出本文所提算法在對目標成像過程中能夠有效地抑制墻雜波，TCR較好，說明信號重建誤差較小，重構(gòu)精度較高，成像質(zhì)量更高。這是因為本文算法的核心思想是回波信號矩陣分解為含目標散射信息的稀疏矩陣及含墻雜波的低秩矩陣和之和，在算法迭代過程中，很好實現(xiàn)了雜波抑制和目標成像相融合。

3 結(jié)束語

本文針對傳統(tǒng)CS成像算法的計算過程復雜、對信噪比要求高等不足提出了一種基于CS框架下的非精確增廣拉格朗日的穿墻雷達成像算法。該算法將穿墻雷達成像問題轉(zhuǎn)化為一個含核范數(shù)和l1范數(shù)的復合優(yōu)化問題，再引入非精確增廣拉格朗日乘子法進行迭代求解，數(shù)值仿真實驗結(jié)果表明，相比傳統(tǒng)成像算法提高了目標信雜比，且可以有效抑制旁瓣，很好地兼顧了重構(gòu)精度和計算效率。而在重建過程中正則化參數(shù)是基于測量數(shù)據(jù)的先驗信息和重建結(jié)果質(zhì)量的經(jīng)驗嘗試判斷確定的，因此對正則化參數(shù)的設(shè)計，將是下一步研究工作的重點。