文獻標識碼： A
DOI：10.16157/j.issn.0258-7998.2015.11.038
中文引用格式： 晉良念，錢玉彬，申文婷，等. 基于改進OMP的超寬帶穿墻雷達稀疏成像方法[J].電子技術應用，2015，41(11)：135-139.
英文引用格式： Jin Liangnian，Qian Yubin，Shen Wenting，et al. Sparse imaging for ultra-wideband through-the-wall radar based on modified OMP algorithm[J].Application of Electronic Technique，2015，41(11)：135-139.
0 引言
穿墻雷達成像（Through the Wall Radar Imageing，TWRI）使用電磁波穿透墻體和內部封閉的建筑物進行監測，以確定建筑結構布局，區分建筑物內部的活動情況，檢測、識別和跟蹤運動目標。TWRI技術在搜救、墻后目標檢測和城區環境的監視等方面備受關注[1-4]。TWRI系統通常利用超寬帶信號和大孔徑天線陣列提供高的距離向和方位向分辨率，能夠在較短的采樣時間內利用較小的存儲空間獲得高質量圖像，與之對應的后向投影一類算法需要大的空間和時間（或頻率）采樣數據才能獲得高質量成像。
近年來，壓縮感知理論的興起，國內外學者將該理論很好地應用到TWRI中，提出了一些高分辨稀疏成像方法[5-6]。它們首先將成像區域劃分為有限個網格，當目標恰好落在網格點上，稀疏成像效果很好。一旦目標偏離了網格點，基矩陣不匹配，從而引起成像模糊，甚至散焦。為此，文獻[5]通過建立off-gird稀疏表示模型，提出了一種基于L21范數懲罰項的FISTA算法，該算法雖能較準確地恢復目標散射系數，但是算法的性能受人工參數設置的影響，在實際應用中很難正確地選擇。文獻[6]針對穿墻雷達偏離網格目標成像中提出了一種貝葉斯推理方法，克服了FISTA算法性能受人工參數設置的影響。但是，該方法在成像區域劃分網格數過大時處理速度會很慢，所以需要尋求新的方法提高處理速度。
本文提出一種改進的正交匹配追蹤方法（GBTOMP算法）實現偏離網格目標超寬帶穿墻雷達稀疏成像。該方法以正交匹配追蹤（Orthogonal Matching Pursuit，OMP）算法為基礎，先通過最速上升梯度優化方法搜索到目標的真實位置，逐漸減小網格偏移量，修正模型中的感知矩陣，再由修正的感知矩陣通過最小二乘方法計算目標散射系數，但在此計算過程中，會引入冗余下標。本文利用貝葉斯假設檢驗去冗余，這樣既可以減小重構過程的運算量，還可以增強重構的精度和抗噪性能。仿真和實驗結果驗證了本文方法的有效性和可行性。
1 稀疏表示模型
一般的穿墻傳播模型可以等效為如圖1所示的兩層傳播模型[7]。假設均勻介質墻體的厚度為d，相對介電常數為?著。將成像區域在方位向和距離向上離散成Nx Ny個網格點，對應的像素值用?滓k，l，k=1，…，Nx，l=1，…，Ny表示，當網格點有目標時k，l不為0，反之k，l等于0。為了方便表示，將k，l按照順序重新排列成Nx Ny維的目標散射系數列向量，令所有元素的網格位置用矩陣表示，其中的元素用p=[xp，yp]T表示，對應的像素值用p表示，則第n個天線回波信號的離散數據模型為：
式中，yn=[yn(t1)，…，yn(tM)]T，An為M×Nx Ny的矩陣，第m行的元素為[s(tn，m-n，1)，…，s(tn，m)]，這里n，p為第n個天線到第p個網格點之間的傳播時延。
根據圖1所示的傳播模型，n，p可以近似表示為：
式中，xn、yn分別表示天線的坐標位置，n，p為傳播的入射角，當成像區域距離天線較遠時可以近似為：
將N個天線到成像區域中所有像素點的傳播時延按照式（2）和式（3）來計算，并把所有天線的接收數據進行堆疊，由此所構造的離散數據模型為：
然而，實際情況中目標很可能偏離預設網格，如圖1所示，這種偏移會影響A′的精確計算。定義成像空間所有網格的偏移量矩陣為，A′關于的函數為A′。另外，對于穿墻成像場景來說滿足稀疏特性，因此可以利用較少的測量值恢復。令表示從y中隨機選擇選擇Q1個天線單元和每個天線單元隨機選擇Q2個數據得到的Q1Q2個數據矢量，并定義這種隨機選擇形式構成的測量矩陣用表示[8]，由此得到大小為Q1Q2×1的測量數據（或稀疏成像）模型（為了書寫方便將A為：
根據式（5）的稀疏成像模型可以看出，感知矩陣A有關，因此求解式（4）反演需要事先準確估計?專參數。當且僅當A滿足約束等距性質時，可以通過求解如下優化問題將高維信號從低維的測量值中準確重構出來。
由于最小0范數問題是一個非確定性多項式難問題，目前的重構算法都是基于0范數的變換或逼近處理，其中貪婪算法以其運算速度快，實現簡單的優勢被用于穿墻雷達成像中。
2 GBTOMP的稀疏成像方法
2.1 GBTOMP算法
本文提出的GBTOMP算法是在正交匹配追蹤算法[9]的迭代過程中采用梯度優化最速上升算法校正網格偏移量，然后采用貝葉斯假設檢驗[10]剔除冗余下標。該算法整體上可以分為兩大步：第一，根據殘差與感知矩陣的內積最大值找到最匹配原子位置（預設網格位置）。由于目標可能會偏離預設的網格位置，因此在該位置需要進行梯度優化的最速上升，以搜索到目標的真實位置，從而以目標真實位置建立新的感知矩陣，并由該矩陣利用最小二乘法求得散射系數向量的近似值，同時通過交替迭代更新下標集求解來得到目標散射系數。第二，OMP算法在有噪聲和沒有噪聲的情況下都會存在冗余下標問題。為解決此問題，GBTOMP算法最后采用貝葉斯假設檢驗理論給出判決準則對第一步輸出的預選下標集中原子進行篩選，得到等于或者大于信號真實下標集估計后，再用最小二乘算法進行重構，以改善算法的性能，增強其抗噪能力。GBTOMP算法流程如下。
初始化：構造測量向量，網格偏移矩陣，感知矩陣A。
(1)采用OMP算法找到匹配原子對應的空間網格位置，利用最速上升梯度優化方法對偏離網格進行估計，再利用最小二乘法對目標散射系數進行粗估計，輸出預選下標集P和粗重構目標散射系數X。
(2)從預選下標集P出發，通過貝葉斯假設檢驗準則設置一個適當門限Thj對P中的元素進行篩選，將大于門限的元素都保留在最終的下標集F中。
輸出：根據最終下標集F，利用最小二乘法得到目標散射系數最終的估計值。
2.2 網格偏移量的估計
在GBTOMP算法的迭代過程中需要估計網格偏移量，其代價函數的構造是關鍵。OMP算法是以測量數據和感知矩陣A構造代價函數的。為了能夠使得估計準確，并使處理過程更加穩定，本文使用OMP算法相關函數的平方開方項來表示[11]：
首先初始化迭代更新次數k=1，殘差向量r。計算感知矩陣A所有列和殘差向量r的內積值，找到內積值最大時所對應的列下標pk，接著進入搜索目標真實位置過程，滿足搜索條件|Ji-Ji-1|<a結束，這里的i代表搜索過程中的迭代次數，a代表很小的一個正數。在此過程中，從?仔中取出第pk列賦給?仔0，以該點坐標(xp，yp)計算出與隨機選擇Q1個天線組成新集合的第q個索引號天線位置之間的雙程傳播時延為：
為求解代價函數J關于所選擇坐標(xp，yp)的梯度，需要事先求解?子q，p關于(xp，yp)的偏導，分兩種情況進行：
(1)當xq=xp時，即q，p=0，q，p簡化為：
(2)當xq≠xp時，即q，p≠0，q，p關于xp和yp的偏導數為：
假設輻射脈沖波形為高斯脈沖s(t)，則第q個天線和所選目標坐標(xp，yp)對應的感知矩陣Aq、B是大小為Q2×1的矩陣，它們的第q2行分別為：
式中是當前(xp，yp)對應的感知矩陣，僅僅是A的一個子集，這樣處理是為了提高運行速度，減小存儲空間。B1x和B1y表示對矩陣分別關于xp和yp的偏導數，它們的具體表達式與s(t)有關。以梯度搜索方法更新所選空間的坐標為：
式中是搜索步長，本文中以網格尺寸為參考的很小的正數。由此得到更新后的重新計算代價函數Ji，搜索迭代次數i=i+1，一旦代價函數滿足條件|Ji-Ji-1|<a就終止搜索過程。終止后，將所選的網格空間位置?仔1所對應的感知矩陣保存在A2矩陣中，同時將所選的感知矩陣列下標pk添加到預選下標集P中，由A2通過最小二乘法計算出第k次粗估計的散射系數xpre，并更新殘差向量。
繼續搜索下標并校正網格偏移量，如此循環迭代，k=k+1，直到殘差向量r滿足條件，粗估計結束。將k次循環求得的散射系數xpre對應其下標集P賦給Nx Ny維列向量X。
X{P}=xpre(18)
2.3 貝葉斯假設檢驗去冗余
為剔除OMP算法輸出的預選集P中的冗余下標，根據貝葉斯假設檢驗模型[10]構造似然函數，其似然比檢驗公式為：
根據判決準則可知，當(zj-mj)2>Thj時，下標j對應的散射系數?滓j的值為非零，更新此時的j為最終下標，反之，下標j對應的散射系數?滓j為零，則從預選下標集中剔除此時的j值。對預選下標集中元素進行篩選，最終剩下的真實下標集為F，接下來利用最小二乘法進行目標散射系數重構完成整個GBTOMP算法。
3 仿真與實驗結果分析
3.1 仿真結果分析
利用GPRMAX產生仿真數據[12]，發射信號采用高斯脈沖，脈沖中心頻率為1 GHz。墻體厚度為0.20 m，相對介電常數為6，電導率為0.03 S/m。
場景1：天線陣列由16個收發共置單元組成，隨機選擇其中的10個天線單元，陣元間距為0.1 m。成像區域是[-1 m，1 m]×[0.3 m，1.3 m]，距離向和方位向的網格尺寸均為0.05 m。假設成像區域中有三個點目標，其坐標分別為(-0.414 m，0.514 m)，(0.025 m，0.725 m)，(0.325 m，0.965 m)。圖2為OMP方法和本文GBTOMP方法成像結果。從圖2(a)可以看出，OMP方法成像出現嚴重散焦，同時帶來虛假像，從圖2(b)可以看到，估計的目標位置和預設目標位置完全吻合，同時目標像聚焦程度高，很明顯本文方法性能優于傳統OMP方法。
場景2：天線陣列由45個收發共置單元組成，隨機選擇其中的10個天線單元，陣元間距為0.04 m，與墻體之間的距離為0.2 m。圖3給出了2維GPRMAX模型，圖中左邊的長方形目標中心位于(0.7 m，3.25 m)，長為0.4 m，寬為0.3 m，右邊的圓形目標體的中心位置位于(1.25 m，3.05 m)，半徑為0.2 m。GPRMAX網格單元尺寸為0.005 m，采樣時窗為30 ns。假設回波信號由有目標時減去無目標時的信號進行模擬。
圖4給出了兩個擴展目標的成像結果，成像區域的方位向和距離向的網格尺寸設為0.03 m。圖中的矩形和圓形虛線框代表兩目標的真實位置。從圖4可以看出，傳統OMP方法的成像效果較差，雖然能夠大致識別目標位置，但是引入了很多的虛假目標像，而本文方法GBTOMP成像效果明顯改善，目標像清晰。
3.2 實驗結果分析
使用美國GSSI公司的探地雷達SIR-20搭建穿墻實驗場景，如圖5所示。實驗墻體為實心磚墻，墻體厚度0.20 m，相對介電常數6.4。選用1 GHz的喇叭天線，架高1.2 m，貼著墻體。水平移動喇叭天線共掃描21個測點，數據處理時隨機選擇10個測點，測點間距為0.05 m，合成孔徑長度為1 m，要求在每個天線孔徑測試點處測量2次，包括有人和無人的場景。人體目標高度為1.80 m，體型寬度約0.40 m，站在墻后1 m處。SIR-20系統數據采集的參數設置：每道采樣1 024點，每秒60道，時間窗15 ns。
首先將SIR-20系統采集的數據波形做平均處理，然后經過去噪和雜波相消、自動增益控制等信號處理得到較好的目標回波，所有成像都是以天線陣列中心為坐標原點。圖6給出實驗數據成像結果比較，圖中的矩形虛線框代表目標的真實位置。從圖中可以看出，傳統OMP方法目標成像散焦現象較嚴重，目標區域不能確定，且周圍還有較多的虛假像。而本文GBTOMP方法成像則不同，雜波相對較少，且目標位置清晰可見，由此看來通過校正目標偏離網格能夠很好解決傳統稀疏成像所出現的散焦現象。
4 結論
本文提出的GBTOMP稀疏成像方法，交替迭代得到目標真實空間位置和目標散射系數，目標位置的校正保證了基矩陣的匹配，實現目標的高質量稀疏成像。仿真和實測數據處理結果表明：與OMP成像相比，本文GBTOMP方法成像效果在目標位置準確度還是目標聚焦程度有明顯的提高。但是，由于以點稀疏模型對擴展目標進行成像，存在圖像不夠平滑，也不連續等問題，下一步工作將對擴展目標的稀疏成像方法展開研究。
參考文獻
[1] TIVIVE F H C，BOUZERDOUM A，AMIN M G.A subspaceprojection approach for wall clutter mitigation in through--wall radar imaging[J].IEEE Transactions on Geoscienceand Remote Sensing，2015，53(4)：2108-2122.
[2] LI G，BURKHOLDER R J.Hybrid matching pursuit for distributed through-wall radar imaging[J].IEEE Transactionson Antennas and Propagation，2015，63(4)：1701-1711.
[3] 蔣留兵，韋洪浪，許騰飛，等.EEMD生命探測雷達人體數量識別技術[J].電子技術應用，2014，40(5)：122-125.
[4] LIU J，CUI G，JIA Y，et al.Sidewall detection using mul-tipath in through-wall radar moving target tracking[J].IEEEGeoscience and Remote Sensing Letters，2015，12(6)： 1372-1376.
[5] XIA S，LIU F.Off-grid compressive sensing through-the-wall radar imaging[C].SPIE Defense+Security.International Society for Optics and Photonics，2014：90771F-90771F-8.
[6] 晉良念，錢玉彬，劉慶華，等.超寬帶穿墻雷達偏離網格目標稀疏成像方法[J].儀器儀表學報，2015，36(4).
[7] JIN T，CHEN B，ZHOU Z.Imaging-domain estimation of wall parameters for autofocusing of through-the-wall SAR imagery[J].IEEE Transactions on Geoscience and Remote Sensing，2013，51(3)：1836-1843.
[8] LAGUNAS E，AMIN M G，AHMAD F，et al.Joint wall mitigation and compressive sensing for indoor image reconstruction[J].IEEE Transactions on Geoscience and Remote Sensing，2013，51(2)：891-906.
[9] 李少東，裴文炯，楊軍，等.貝葉斯模型下的OMP重構算法及應用[J].系統工程與電子技術，2015，37(2).
[10] 甘偉，許錄平，蘇哲，等.基于貝葉斯假設檢驗的壓縮感知重構[J].電子與信息學報，2011，33(11)：2640-2646.
[11] GURBUZ A C，TEKE O，ARIKAN O.Sparse ground-penetrating radar imaging method for off-the-grid target problem[J].Journal of Electronic Imaging,2013，22(2)：021007-021007.
[12] XIE J L，XU J L.Ground penetrating radar-based experimental simulation and signal interpretation on roadway roof separation detection[J].Arabian Journal of Geosciences，2014：1-8.