基于魯棒Capon波束形成的探地雷達成像算法

周 琳 黃春琳 粟 毅

(國防科技大學電子科學與工程學院 長沙 410073)

1 引言

探地雷達(Ground Penetrating Radar, GPR)成像技術可以對地下目標進行非破壞性的高分辨率檢測和識別，已在軍事和民用領域得到廣泛應用[1,2]。目前常用的 GPR成像方法主要有反向投影(Back Projection, BP)算法[3,4]、距離遷移(Range Migration,RM) 算法[5,6]、衍射層析(Diffraction Tomography,DT)算法[7,8]等，現有算法均屬于獨立于數據的非自適應信號處理方法。盡管使用這些算法可以實現比較高的分辨率，但與依賴于數據的自適應信號處理方法相比，后者具有更高的旁瓣和雜波抑制能力。

Capon波束形成方法(Capon Beamforming, CB)作為典型的自適應信號處理方法，已廣泛應用于陣列信號處理、SAR成像等領域[9]。但是在某些實際應用中，信號方向矢量往往不能準確估計，此時 CB的性能將急劇下降，甚至比非自適應方法的性能還差。為此許多文獻對改進CB的魯棒性開展了研究，試圖解決由方向矢量不確定性導致的估計不準問題，例如基于線性約束的協方差矩陣漸變方法[10]，基于二次約束的對角加載方法[11]，基于噪聲協方差矩陣先驗知識的子空間自適應波數形成方法[12]等，但這些方法并未直接針對方向矢量的不確定性。

魯棒 Capon波束形成(Robust Capon Beamforming, RCB) 方法在CB的基礎上考慮了基于球形或橢球形不確定集的方向矢量，克服了原有CB進行估計時魯棒性差的缺點，在穿墻成像、醫學成像等領域已得到成功應用[13-16]。本文將上述自適應方法應用于GPR成像處理中，提出了基于RCB理論的GPR成像算法，并得到了很好的旁瓣和雜波抑制效果。全文按照如下結構組織，第2節介紹了RCB基礎理論，第3節提出了基于RCB的GPR成像方法，第4節給出了實測數據實驗結果，最后是結束語。

2 RCB基礎理論

對于一個包含M個陣元的線性陣列，假設其接收數據的協方差矩陣為R, RCB方法的關鍵是獲得陣列加權矢量w，該加權矢量可通過在如下線性約束條件下尋找最小值得到[13]：

其中as為陣列輸出方向矢量，上標H在本文中均表示共軛轉置。由于實際情況中as不能精確已知，關于as僅有的先驗知識是其在附近一個很小的鄰域內波動，即，其中e為一個很小的正數。令s2為后向散射信號能量，則上式可轉化為

由于式(2)中s2與as均為未知，故二者之間存在尺度模糊，為消除這種影響，應使用如下范數約束[14]：

對于任意時刻t，式(2)的解為

這樣式(2)可退化為如下二次最優問題：

式(5)中為避免as=0，需設定限制條件[14]：

為求解式(5)，使用 Lagrange乘數法，令目標函數為

其中l≥0為實Lagrange乘數，且滿足R-1+lI正定。這樣對于固定的l,f(as,l)最小化的條件是[14]

并且有

聯立式(8)與式(9)可得到關于 Lagrange乘數l的方程：

為求解式(10)，首先可將協方差矩陣R進行特征值分解：

其中U為R的特征矢量矩陣，其每一列都為R的特征矢量，Γ為對角化特征值矩陣，其主對角線上的各元素對應于R的M個特征值：g1≥g2≥…≥gM。令z=[z1,z2,… ,zM]T=UH，則式(10)可轉化為

則后向散射信號能量可通過將式(13)代入式(4)求得，而陣列加權矢量w可表示為

3 基于RCB的GPR成像算法

由于 RCB在理論上具有良好的旁瓣和雜波抑制能力，為此本文將其應用到GPR成像中，提出了基于RCB的GPR成像算法。首先建立2維成像場景如圖1所示，整個場景被直線z=0分為兩部分。z＜ 0 的區域為空氣，其介電常數為e1=e0，其中e0是自由空間中的介電常數。z＞0區域為各向同性均勻的土壤，其介電常數為e2=er e0，其中er是土壤的相對介電常數。為簡化模型，假設兩個區域內的磁導率都等于自由空間中的磁導率，即m1=m2=m0，其中m1,m2,m0分別為空氣中、土壤中和自由空間中的磁導率，并假定空氣和土壤均無耗。考慮GPR中典型的B-scan合成孔徑掃描方式，假設合成孔徑數目為M，當前工作的天線位置序號為k，用黑色的三角形表示，坐標為(xk, -h)，其余M-1個孔徑位置由灰色三角形表示。

第k個孔徑處接收到的信號可表示為

其中t= 1Δt, 2 Δt, … ,NΔt, Δt為時域采樣間隔。對圖1中給定的某成像點A，首先要計算出從A到各天線位置處的雙程時延。由于空氣-土壤交界面的影響，該路徑不是一條直線而是一條折線。發射信號從第k個天線位置出發，經過折射點(xr,0)到達點A，再沿原路徑返回天線處。入射角和折射角分別用qi和qt表示。

圖1 成像場景說明

對于任意一個成像點A，根據Snell定律，有如下關系：

由圖1中坐標關系，式(16)可轉化為

式(17)是一個關于xr的四次方程，通過求解xr，可得到雙程時延表達式如下所示：

其中c是自由空間中電磁波的波速，tA,k即對應于點A和第k個天線之間的雙程時延。利用式(18)所得時延，可將原始回波數據進行時延校正如下：

則在某時刻t的M個陣元的校正后接收信號可表示為

為了估計as及w，首先需要估計接收數據的協方差矩陣R。各種基于CB的自適應方法的一個主要差別就是關于協方差矩陣的不同假設及構造方法[15]。一種直觀的構造方法即利用y(t)直接生成協方差矩陣R：

但是因為合成孔徑掃描方式未能對場景提供足夠多的視，按照式(21)構造的協方差矩陣會降低算法抑制旁瓣和雜波的能力，為此需將接收數據劃分為多個子集以增加可用視數[9]。利用文獻[9]提出的劃分子集方法，本文在時域上將M×1維的校正后信號y(t)劃分為L個相互重疊的N×1維的信號yi(t) ,i= 1,2,… ,L，圖 2是該劃分方法的示意圖，由圖2可知L=M-N+ 1。采用降低信號的維度從而生成多視的方法來構造協方差矩陣會一定程度上降低成像分辨率，為減小這種效應，N不應過小，文獻[15]中關于N的經典取值為N=0.8M。使用較大的N會增加算法中矩陣不可逆的可能性，但另一方面會使成像分辨率沒有大幅降低的同時保持算法的高信雜比和高峰值旁瓣比特性。本文中協方差矩陣定義為

圖2 將校正信號劃分為子集示意圖

下面便可按照上節介紹的RCB理論進行GPR成像，當求得加權矢量的估計值后，則t時刻點A處的后向散射信號s(t)為

點A處的后向散射能量可通過式(24)計算[16]：

將后向散射能量p(A)作為點A處的成像值，通過遍歷待成像區域的所有點，便完成成像過程。綜上所述，本文提出的基于RCB的GPR成像算法的具體步驟如下：

步驟 1 對成像區域中任意一點，按照式(18)計算其距離M個合成孔徑位置的雙程時延，再按照式(19)得到校正后信號y(t)；

步驟 2 利用y(t)按照式(22)構造協方差矩陣R；

步驟 3 將R按照式(11)進行特征值分解，從而求得U,Γ及z；

步驟 4 利用Γ和z求解式(12)中的 Lagrange系數l0；

步驟 5 利用l0及式(8)，式(13)，式(14)求解,和；

步驟 7 按照上述步驟，遍歷整個成像區域。

4 實驗結果及分析

第1個實驗使用GPRMax軟件生成的仿真數據進行分析。實驗場景如圖 3(a)所示，兩根鋼筋棍埋設于相對介電常數為4的各向均勻土壤中，鋼筋均沿y軸放置，場景剖面具體的幾何位置關系見圖3(b)。深度為10 cm的鋼筋直徑為0.5 cm，深度為25 cm的鋼筋直徑為1 cm。GPR天線距地面10 cm,B-scan路徑沿x軸方向，每個合成孔徑掃描位置在圖中用黑色圓點表示。發射信號為高斯微分脈沖，中心頻率為1 GHz。

圖3 仿真數據成像實驗的場景說明

圖 4(a)是使用本文所提 RCB算法進行成像的歸一化結果。考慮到 BP算法是傳統的非自適應算法中的典型代表，且可以充分考慮介質的分層現象，其在 GPR成像理論和實際中都已得到深入研究和廣泛應用，故本文選取BP算法的結果與RCB算法進行比較，如圖4(b)所示。

由圖4可以看出，使用RCB算法可得到視覺上非常直觀的旁瓣和雜波抑制效果。為了定量描述雜波抑制的效果，本文使用積分旁瓣比(Integrated Side Lobe Ratio, ISLR)和峰值旁瓣比(Peak Side Lobe Ratio, PSLR)對成像質量進行評價。ISLR定義為目標所有旁瓣能量和主瓣能量的比值：

其中Etotal和Emain分別是圖像中總的能量和目標主瓣能量。PSLR定義為目標第一副瓣峰值和主瓣峰值的幅度之比：

其中Iside和Ipeak分別是圖像目標處旁瓣的最大幅值和主瓣的最大幅值。

表1給出了兩種算法成像結果中各目標的ISLR和PSLR值，可見與BP算法相比，RCB算法各個指標均有所下降，這說明RCB算法的成像結果中能量更加集中在目標區域，即其更好地抑制了非目標區域的干擾。在圖4中通過觀察兩種成像算法結果中目標峰值處的主瓣-3 dB寬度，對于左上方目標，在x方向上，RCB算法為BP算法的0.81倍，在z方向上，RCB算法為BP算法的0.80倍；對于右下方目標，在x方向上，RCB算法為BP算法的0.89倍，在z方向上，RCB算法為BP算法的0.67倍。主瓣寬度的縮小進一步說明 RCB算法在分辨率上也較BP算法有一定程度的提高。

表1仿真數據成像結果性能對比

第2個實驗使用實測數據進行分析，一個裝滿純凈水的礦泉水瓶被埋于沙坑中，如圖 5(a)所示。沙子的相對介電常數為 2.37。水瓶的長度方向沿y軸放置，中心深度為 17 cm，水瓶主體部分橫截面直徑為7 cm。天線距地面40 cm, B-scan路徑沿x軸方向，橫穿水瓶中心，每個合成孔徑掃描位置在圖5(a)中用黑色圓點表示，水瓶的B-scan橫截面示意圖如圖5(b)所示。本實驗中所使用的發射信號與實驗1中相同。

圖4 仿真數據成像結果對比

圖5 實測數據成像實驗的場景說明

圖6是分別使用RCB算法和BP算法進行成像的結果，表2給出了兩種算法成像結果的ISLR和PSLR，再次說明了RCB算法的干擾抑制性能。

為了更好地觀察到 RCB算法在分辨率上的性能，圖7給出了兩種算法成像結果中峰值點位置在x和z兩個方向的剖面圖。考慮目標峰值處的主瓣-3 dB寬度，在x方向上，RCB算法為BP算法的0.75倍，在z方向上，RCB算法為BP算法的0.71倍，再次說明 RCB算法可以在一定程度上提高成像分辨率。

表2 實測數據成像結果性能對比

5 結束語

現有的 GPR成像算法均屬于獨立于數據的非自適應信號處理方法，考慮到依賴于數據的自適應信號處理方法在干擾抑制方面的優良性能，本文將RCB理論應用于GPR成像中，不僅大幅降低了成像結果中旁瓣和雜波的能量，并在一定程度上提高了成像分辨率。但是在成像質量提高的同時，由于RCB算法中涉及大量的矩陣求逆運算，其計算量也急劇增加，為此如何提高RCB算法的計算效率將成為下一步的研究重點。

圖6 實測數據成像結果對比

[1]粟毅, 黃春琳, 雷文太. 探地雷達理論與應用[M]. 北京: 科學出版社, 2006: 1-10.

Su Yi, Huang Chun-lin, and Lei Wen-tai. Theory and Applications of Ground Penetrating Radar [M]. Beijing:Science Press, 2006: 1-10.

[2]Khan U S, Al-Nuaimy W, and Abd El-Samie F E. Detection of landmines and underground utilities from acoustic and GPR images with a cepstral approach [J].Journal of Visual Communication and Image Representation, 2010, 21(7):731-740.

[3]Zhou L and Su Y. A GPR imaging algorithm with artifacts suppression[C]. Proceedings of the 13th Internarional Conference on Ground Penetrating Radar, Lecce, Italy, 2010:331-337.

[4]Chen Lei and Shan Ouyang. A time-domain beamformer for UWB through-wall imaging [C]. IEEE Region 10 Conference,Taipei, China, 2007: 1-4.

[5]Zyada Z, Matsuno T, Hasegawa Y,et al.. Advances in GPR-based landmine automatic detection [J].Journal of the Franklin Institute, 2011, 348(1): 66-78.

[6]Schmelzbach C, Tronicke J, and Dietrich P. Threedimensional hydrostratigraphic models from groundpenetrating radar and direct-push data [J].Journal of Hydrology, 2011, 398(3/4): 235-245.

[7]Dafflon B, Irving J, and Barrash W. Inversion of multiple intersecting high-resolution crosshole GPR profiles for hydrological characterization at the Boise Hydrogeophysical Research Site(BHRS)[J].Journal of Applied Geophysics, 2011,73(4): 305-314.

[8]Hugenschmidt J, Kalogeropoulos A, Soldovieri F,et al..Processing strategies for high-resolution GPR concrete inspections [J].NDT&E International, 2010, 43(4): 334-342.

[9]Benitz G R. High-definition vector imaging [J].Proceedings of the31st Lincoln Laboratory Journal, 1997, 10(2): 147-170.

[10]Zatman M. Comments on theory and applications of covariance matrix tapers for robust adaptive beamforming [J].IEEE Transactions on Signal Processing, 2000, 46(6):1796-1800.

[11]Wu R, Bao Z, and Ma Y. Control of peak sidelobe level in adaptive arrays [J].IEEE Transactions on Antennas and Propagation, 1996, 44(10): 1341-1347.

[12]Lee C C and Lee J H. Robust adaptive array beamforming under steering vector errors [J].IEEE Transactions on Antennas and Propagation, 1997, 45(1): 168-175.

[13]Li J, Stoica P, and Wang Z. On robust Capon beamforming and diagonal loading [J].IEEE Transactions on Signal Processing, 2003, 51(7): 1702-1715.

[14]Xie Y, Guo B, Xu L,et al.. Multi-static adaptive microwave imaging for early breast cancer detection [C]. Proceedings of 39th ASILOMAR Conference on Signals, Systems and Computers, Pacific Grove, Calif, USA, 2005: 285-289.

[15]Ahmad F, Amin M G, and Kassam S A. A beamforming approach to stepped-frequency synthetic aperture throughthe-wall radar imaging [C]. Proceedings of 1st IEEE International Workshop on Computational Advances in Multi-Sensor Adaptive Processing (CAMSAP 2005), Puerto Vallarta, Mexico, 2005: 24-27.

[16]Xie Y, Guo B, Li J,et al.. Novel multistatic adaptive microwave imaging methods for early breast cancer detection[J].EURASIP Journal on Applied Signal Processing, 2006:1-13.