景素雅，晉良念,2，劉慶華

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

0 引言

近年來，穿墻雷達成像在軍事和救援安全中具有廣泛的應用，建筑物透視雷達技術得到了來自科研界和工業(yè)界的廣泛關注和研究，并在這十多年間取得了快速的發(fā)展。現(xiàn)有穿墻雷達成像算法有后向投影(Back Projection，BP)成像算法，延時求和(Delay and Sum，DAS)和壓縮感知稀疏成像算法。BP 算法[1]由于成像原理簡單，易于實現(xiàn)，廣泛應用于穿墻雷達成像領域。DAS[2]由于其快速、直觀的實現(xiàn)，成為雷達應用中最常用的圖像形成技術之一，其通過在特定像素點處將由所有通道引起的后向散射貢獻累加，來估計該像素點的反射系數(shù)。然而BP算法和DAS算法都存在成像旁瓣多、分辨率低等問題。

由于在穿墻探測場景中目標相對于整個成像空間具有稀疏性，因此壓縮感知被用于其成像。常用到的稀疏成像算法有貪婪類算法、凸優(yōu)化以及貝葉斯算法。文獻[3]提出了一種在貪婪類正交匹配追蹤(Orthogonal Matching Pursuit，OMP)基礎上改進的動態(tài)閾值弱正交匹配追蹤算法，該方法成像的性能有所改善。然而貪婪類方法對弱目標的成像不太友好，BP成像中的柵旁瓣問題同樣會出現(xiàn)在此類算法中[4]。文獻[5]通過將前一雷達位置的測量值作為邊信息加入到當前位置雷達處的場景重構中，有效地實現(xiàn)了隱藏在墻后的靜態(tài)人體目標檢測。然而該方法只是針對點目標的成像?？紤]到實際場景中擴展目標更為常見，文獻[6]采用可分離逼近結構稀疏恢復算法，充分利用擴展目標信號的結構稀疏先驗對其進行稀疏成像重構。然而該方法面臨著需要預先設置懲罰參數(shù)的問題。為了在保留擴展目標邊緣特性的同時, 避免預先人工設置參數(shù)，文獻[7]進一步提出全變分(Total Variation，TV)約束下的擴展目標稀疏成像方法，利用了目標的先驗信息，通過“Integrate-out”方法得到擴展目標的最大后驗(Maximum a Posteriori，MAP)估計，然后交替循環(huán)迭代至最優(yōu)解。然而，在穿墻雷達稀疏成像中，由于字典矩陣的存在，成像過程中面臨字典存儲所占容量大、多次字典求逆和乘法運算導致計算復雜度高等問題。

針對以上問題，本文針對TV約束的MAP稀疏成像方法，考慮到字典矩陣的構建由發(fā)射信號和時延共同決定，通過對發(fā)射波序列進行變換，將迭代式中包含字典矩陣的所有項用線性卷積和哈希(Hash)表進行代替，從而將迭代循環(huán)中繁雜的計算過程進行簡化。這種方法可以有效地減少運算時間，并在空間復雜度方面得到了改善，同時在成像質量性能上也得到了提高。

1 信號模型

使用L個收發(fā)共置天線單元合成線性陣列，對SOI(scene-of-interest)進行探測，陣列平行于墻體。將SOI離散化為N個像素點，堆疊L個天線的數(shù)據(jù)，得到雷達回波信號模型為

y=Ax+e

(1)

Al(m,i)=s(mTs-τli)

(2)

采用包含TV約束項的MAP成像模型來重構x，既保留了擴展目標邊緣特性，又避免了LASSO問題中懲罰參數(shù)需預先人工設定參數(shù)問題[7]。假設噪聲矢量e服從e～N(0,βI)，βI為噪聲的協(xié)方差矩陣，則TV約束下關于x的MAP的目標函數(shù)為[7]

Nlog(‖Dx‖1)+MLlogβ

(3)

2 算法描述

首先對式(3)的前三項分別進行優(yōu)化最小化(Majorize-Minimize，MM)替代[9]，然后對替代后的目標函數(shù)J(x,β)求解關于x和參數(shù)β的偏導數(shù)，并令其倒數(shù)為0，得到β和x的第i個元素xi的迭代式[10]，即

(4)

和

(5)

2.1 β(t+1)的計算

m=1,2,…，M

(6)

記

gl[m]?Alx(t)[m]

(7)

s(mTs)?u[m]

(8)

利用沖激函數(shù)的抽樣特性，得到

(9)

(10)

(11)

(12)

令z[M+1-m]?s(mTs)，則有

(13)

2.3 Hi的計算

由式(2)可知，

(14)

令f[M+1-m]?s2(mTs)，則有

(15)

同樣，Hi的運算轉換為nl(nl[m]表示Al的第m行中非零元素的數(shù)量)和f的卷積計算，也避免了A的構建、存儲和計算。而且，通過對這里的s2(mTs)進行替換，使所提方法同樣對任何超寬帶發(fā)射波形都可以將Hi轉換為卷積的計算，而不限于對稱的發(fā)射波。

值得一提的是，如果要快速得到Hi，那么nl的計算是至關重要的。因為發(fā)射信號是超寬帶信號，也就是時域為短時信號，所以對于長度為M的發(fā)射序列u[m](即s(mTs))而言，存在一個整數(shù)Q，當滿足mTs≤QTs時，使得u[m]不為零?？紤]矩陣Al的第i列為Al=[s(Ts-τli),s(2Ts-τli),…,s(MTs-τli)]T，則元素Al(m,i)非零當滿足

1≤m-mli≤Q

(16)

且離散時延mli滿足

mmin≤mli≤mmax

(17)

所以有

max(mmin,m-Q)≤mli≤min(mmax,m-1)

(18)

即

(19)

式中，φn={i=1,2,…,N|n=mli}表示第l個天線關于N個像素網(wǎng)格點的時延，它是滿足{mli=n|mli=[ml1,ml2,…,mlN]}的mli的個數(shù)，可通過Hash表[11]求得。

3 復雜度分析

從前面的分析可知，若文獻[10]采用100%回波數(shù)據(jù)重構，在迭代更新β時，它和本文的不同點體現(xiàn)在計算Ax上。文獻[10]更新Ax所需乘法MLN次，所需加法為MLN-ML次，而采用本文方法需調用L次卷積函數(shù)和L次accummary函數(shù)(求pl[n]時用到)。文獻[10]更新Hi的過程中涉及到求解字典矩陣A中每行非零元素個數(shù)，只能對字典A逐行搜索求解，則本文方法利用發(fā)射波的短時特性，只需調用L次accummary函數(shù)計算nj,j=1,2,…,ML。在已知nj的情況下，文獻[10]更新全部的H=[H1,H2,…,HN]所需乘法為2MLN次，所需加法為MLN-N次，而采用本文方法需調用L次卷積函數(shù)，然后對不同的像素點i，取式(15)中L個卷積序列所對應的第γ個元素進行累加，即可得到所有的Hi。在迭代更新G時，更新全部的G=[G1,G2,…,GN]所需乘法為2MLN次，所需加法為2MLN-N次，而采用本文方法需調用L次accummary函數(shù)(求pl[n]時用到)，2L次卷積函數(shù)以及ML次加法運算，然后對不同的i，取式(13)中L個卷積序列所對應的第γ個元素進行累加，即可得到所有的Gi，降低了計算復雜度。

4 仿真及實驗結果分析

4.1 仿真結果

仿真場景模型如圖1(a)和圖2(a)所示，距離天線1 m處長為3 m的墻體由均勻介質材料構成，其墻厚為0.2 m、相對介電常數(shù)為4.5、電導率為0.01 S/m。由GprMax電磁仿真軟件產(chǎn)生仿真數(shù)據(jù)，仿真中采用發(fā)射波中心頻率為1 GHz的窄高斯脈沖信號。仿真場景中，矩形目標體長均為0.5 m、寬為0.2 m，目標距墻體距離為1 m。

以單個目標為例，表1給出了文獻[10]分別采用100%、50%和10%的回波數(shù)據(jù)以及本文方法在單次循環(huán)(即更新一次x和β)的時間(Per-Time)、總運行時間(Total Time)和分配的內存(Total Memory)3個方面的對比。從表1可知，對比文獻[10]采用100%回波數(shù)據(jù)重構時的結果，本文方法的運算時間減少2/3左右，分配內存減少3/5左右；而相比文獻[10]采用50%回波數(shù)據(jù)重構的結果，運算時間減少3/5左右，分配內存減少2/5左右；即使文獻[10]采用10%的回波數(shù)據(jù)來重構，運算時間仍減少達2/5左右，分配內存也小到3/4左右。由此看來，本文方法有效地減少了運行時間和分配的內存。

表1 Per-Time、Total Time和Total Memory的比較

圖1(b)到圖1(f)分別為單個目標情況下采用BP算法和文獻[10]100%、50%、10%的回波數(shù)據(jù)以及本文方法的重構結果。圖2(b)到圖2(d)給出了兩個目標情況下BP算法、文獻[10](采用100%的回波數(shù)據(jù))以及本文方法的重構結果?？梢钥闯?，傳統(tǒng)BP方法成像效果最差，雜波最多。本文方法成像最為清晰，且目標邊緣輪廓明顯。為了更直觀地對重構質量進行對比，表2給出單個目標情況下其在TCR、ENT和MFOD三方面的對比，分別反映了目標雜波比、圖像熵和目標的邊緣輪廓特性。對比發(fā)現(xiàn)，本文方法和文獻[10]方法的TCR都比BP成像大，ENT都比BP成像小，說明這兩種方法的成像對比傳統(tǒng)BP方法雜波少，圖像清晰。本文方法的TCR最大，ENT最小，成像性能最優(yōu)。文獻[10]在采用3種不同的回波數(shù)據(jù)重構時，100%的回波重構時性能相對最好，而此時所提方法仍略優(yōu)于文獻[10]。MFOD反映了區(qū)域的平滑性，其值越小，目標的邊緣輪廓特性越強，結合MFOD和表1中總運行時間對比可知，本文方法在減少運算時間的同時，并沒有損壞擴展目標的邊緣特性。

表2 TCR、ENT和MFOD的比較

(a)仿真場景

(a)仿真場景

4.2 實測結果

實際探測場景如圖3所示，實測數(shù)據(jù)集采用維拉諾瓦(Villanova)大學高級通信中心(CAC)與空軍研究實驗室(AFRL)合作收集得到的穿墻雷達成像實驗數(shù)據(jù)集。數(shù)據(jù)使用型號為ENA 5071B的安捷倫網(wǎng)絡分析儀收集。實驗采用以2.5 GHz為中心的帶寬為1 GHz的步進頻率波形，步進間隔為5 MHz，文獻[12]給出了場景的詳細介紹。本文采用其中的一個場景進行實驗驗證，選取平面陣列中高度位于中心位置的一排天線所得到的數(shù)據(jù)進行二維平面成像。圖4和圖5分別為傳統(tǒng)BP和本文方法成像結果，可以發(fā)現(xiàn)，后者重構結果雜波少。在運行時間上，采用BP重構、文獻[10](采用100%的回波數(shù)據(jù))重構和本文方法重構時間分別為3.339,121.346和38.473 s; 后兩者在運行所需內存上分別為2 231 542 Kbit和914 623 Kbit?？梢园l(fā)現(xiàn)，對比文獻[10]，本文方法有效地降低了運行時間和所需內存。

圖3 實驗場景

圖4 BP成像

圖5 本文方法的重構結果

5 結束語

本文提出一種在MM優(yōu)化的TV約束MAP估計框架下，無需構建感知字典的超寬帶穿墻雷達稀疏成像的快速實現(xiàn)方法。仿真結果和實驗結果表明，通過與BP成像和文獻[10]方法重構的圖像相比，得到的圖像有更少的噪聲，能捕捉SOI內主要的散射點，有效抑制了旁瓣和柵瓣。而且，在提升運算速度和占用內存更少的同時，有效地保留了擴展目標的邊緣特性。該方法為大數(shù)據(jù)量成像系統(tǒng)，如三維成像提供了重要的參考。