基于IALM的穿墻雷達(dá)成像算法

張珊珊，趙建華

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

穿墻雷達(dá)成像（Though-the-Wall Radar Imaging，TWRI）是一種利用低頻電磁波的穿透性對復(fù)雜場景或障礙物后隱藏目標(biāo)進(jìn)行有效探測、定位和成像的技術(shù)，在應(yīng)急救援、醫(yī)療監(jiān)護(hù)、城市安保和反恐斗爭等領(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]。此外，雷達(dá)成像近似視為對探測場景進(jìn)行的最小二乘估計，使得基于匹配濾波思想的傳統(tǒng)成像算法會出現(xiàn)較高旁瓣，使得成像結(jié)果模糊，無法區(qū)分鄰近目標(biāo)，成像質(zhì)量不高[4]。

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

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

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

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

1.1 CS基本原理

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

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

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

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

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

信號重構(gòu)可通過對式（3）求逆運(yùn)算得系數(shù)α，再恢復(fù)原始信號x。由于Φ的列數(shù)遠(yuǎn)大于行數(shù)，使得未知數(shù)個數(shù)遠(yuǎn)大于求解方程個數(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）的思想，利用單站雷達(dá)的運(yùn)動產(chǎn)生與目標(biāo)之間的相對位移，把小的實(shí)天線孔徑合成為一個等效的大虛擬孔徑，提高方位向分辨率。如圖1所示為穿墻雷達(dá)系統(tǒng)探測示意圖。

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

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

其中，f0為初始頻率，Δf為頻率步進(jìn)量。

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

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

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

圖2 基于幾何光學(xué)的穿墻信號模型

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

將式（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可通過奇異值軟閾值法進(jìn)行求解，s可通過l1范數(shù)最小化技術(shù)進(jìn)行計算。λw、λs分別為低秩矩陣和稀疏矩陣的正則化參數(shù)。但該方法涉及奇異值分解運(yùn)算較為復(fù)雜，導(dǎo)致閾值迭代速度較慢。

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

1.4.1 增廣拉格朗日乘子法

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

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

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

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

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

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

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

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

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

在實(shí)際迭代過程中，不需要求得A和E的精確值，而只需迭代一次就可恢復(fù)原始矩陣的局部最優(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?()表示矩陣化運(yùn)算，表示把NM×1維向量轉(zhuǎn)化為N×M維矩陣；Φ?為Φ的偽逆矩陣；Y0=0，S0=0，μ0=0，ρ>1，k=0。

（2）更新Zw：

（3）更新S：

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

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

2 仿真結(jié)果分析

為驗(yàn)證本文所提成像算法的有效性和準(zhǔn)確性，利用圖1所建立的穿墻雷達(dá)探測模型進(jìn)行算法的驗(yàn)證。本文所有仿真實(shí)驗(yàn)均在Windows 10系統(tǒng)下進(jìn)行，成像部分的仿真均在MATLAB R2014a環(huán)境下進(jìn)行。仿真模型設(shè)置如下：雷達(dá)起始頻率fl設(shè)置為1 GHz，步進(jìn)間隔為10 MHz，終止頻率fh為3 GHz，則步進(jìn)頻信號共有201個采樣頻點(diǎn)。合成孔徑長度為2 m，每隔5 cm設(shè)置一個天線，共有41個天線測量位置。設(shè)墻體厚度為0.2 m的混凝土材質(zhì)，相對介電常數(shù)為6，且距離天線0.15 m，墻體散射系數(shù)設(shè)為1，目標(biāo)散射系數(shù)設(shè)為0.1，并加入10 dB噪聲。設(shè)理想目標(biāo)位于水平向1 m、距離向2 m處，邊長為0.05 m。設(shè)置成像區(qū)域?yàn)? m×3 m，將成像區(qū)域平均分為60×60個網(wǎng)格。從201×41=8 241個頻域樣本數(shù)據(jù)中隨機(jī)抽取出1 648個數(shù)據(jù)（保證降采樣）進(jìn)行成像。迭代過程中，設(shè)置正則化參數(shù)λ為0.2，懲罰系數(shù)μ為1，ρ為2，1E?6為迭代運(yùn)算的判斷閾值。重構(gòu)結(jié)果均采用歸一化到[0，1]區(qū)間的圖像強(qiáng)度像素值的絕對值來表示。如圖5所示為仿真成像結(jié)果對比，其中圖5（a）為目標(biāo)真值圖，圖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)目標(biāo)像，驗(yàn)證了該理論在穿墻雷達(dá)成像中具有的利用降維采樣數(shù)據(jù)恢復(fù)原始回波信號的獨(dú)特優(yōu)勢。其中，圖5（b）與圖5（c）中清晰可見墻體的位置，而圖5（d）中墻體回波較弱，不易直接從圖中觀察到，對比可以看出，利用直接CS成像的效果最差，成像結(jié)果中有大量噪聲，目標(biāo)受周圍雜波影響較大，目標(biāo)像分辨率較低，而直接ISTA成像效果最好，但圖中目標(biāo)像的模糊散焦現(xiàn)象較嚴(yán)重。而抑制雜波處理后提高了目標(biāo)成像質(zhì)量，其中，分別對比圖5（b）與圖5（e）、圖5（c）與圖5（f）可以看出，抑制雜波后目標(biāo)成像質(zhì)量更高，說明了很好地消除墻體對目標(biāo)回波的影響，重構(gòu)場景中雜波較少，目標(biāo)像較清晰，圖像質(zhì)量得到了很大的提高，對比圖5（d）與圖5（g）可看出，目標(biāo)成像質(zhì)量基本沒有太大的改變，抑制雜波處理的效果不大，說明了抑制雜波對直接ISTA成像提高成像質(zhì)量影響較小。

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

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

TCR定義為：

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

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

從表1對比看出，各種算法在抑制雜波后均可提高成像質(zhì)量，但相應(yīng)地也增加了運(yùn)行時間。其中，直接CS成像由于調(diào)用CVX工具箱，使得運(yùn)行速度較慢。直接CS成像在未進(jìn)行墻雜波抑制處理時的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成像在運(yùn)行時間上相差不大，但I(xiàn)STA類算法在此場景下恢復(fù)目標(biāo)信號的信雜比更高。

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

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

表1 成像重建性能對比

從仿真成像結(jié)果和重建性能對比中可看出本文所提算法在對目標(biāo)成像過程中能夠有效地抑制墻雜波，TCR較好，說明信號重建誤差較小，重構(gòu)精度較高，成像質(zhì)量更高。這是因?yàn)楸疚乃惴ǖ暮诵乃枷胧腔夭ㄐ盘柧仃嚪纸鉃楹繕?biāo)散射信息的稀疏矩陣及含墻雜波的低秩矩陣和之和，在算法迭代過程中，很好實(shí)現(xiàn)了雜波抑制和目標(biāo)成像相融合。

3 結(jié)束語

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