脉冲探地雷达后向投影成像方法的研究

来五星 唐文静 孙少山 钟升

摘 要： 对脉冲探地雷达数据进行成像处理有利于地下目标的定位和探测。详细阐述了探地雷达后向投影成像方法的实现原理与过程，并针对该方法存在计算效率低的问题，提出基于子孔径划分的改进后向投影成像方法，最后用频率?波数成像、有限差分成像、标准后向投影成像以及改进的后向投影成像方法分别对仿真和实测数据进行成像处理并比较各成像方法的性能，验证了改进算法的有效性。

关键词： 探地雷达； 成像方法； 后向投影； 子孔径划分； 合成孔径； 成像效率

中图分类号： TN958?34 文献标识码： A 文章编号： 1004?373X（2018）15?0061?04

Research on back?projection imaging method of impulse ground penetrating radar

LAI Wuxing， TANG Wenjing， SUN Shaoshan， ZHONG Sheng

（School of Mechanical Science and Engineering， Huazhong University of Science and Technology， Wuhan 430074， China）

Abstract： The imaging processing of impulse ground penetrating radar （GPR） data is helpful for the location and detection of underground targets. The implementation principle and procedure of the standard back?projection imaging method of GPR are described in detail， and an improved back?projection imaging algorithm based on sub?aperture division is proposed to improve the computational efficiency of the standard imaging method. The frequency and wave number imaging， finite difference imaging， standard back?projection imaging and improved back?projection imaging methods are used to process the simulation data and actual measured data， their performances are compared， and the effectiveness of the improved algorithm is verified.

Keywords： ground penetrating radar； imaging method； back?projection； sub?aperture division； synthetic aperture； imaging efficiency

0 引 言

探地雷达（GPR）也称地质雷达，是对埋地目标进行定位与探测的有效工具，它集探测速度快、分辨率高、穿透力强、无损检测等优点于一身，被广泛应用于市政管线检测、道路与桥梁灾害检查、考古发掘探测、地雷与未爆弹药的排查、矿物资源勘探等诸多方面[1?2]。脉冲探地雷达的原理是发射天线向地下辐射脉宽为纳秒级的高频电磁波，通过对接收到的回波信号的幅度、频率以及相位等特征进行分析处理，得到被测物体的位置和大小信息。探地雷达成像就是对接收到的回波信号进行合成孔径处理，从而准确地反映地下目标的位置。由于探地雷达回波数据与地震波数据具有相似性，地震数据处理使用的偏移技术可以应用在探地雷达的合成孔径成像处理中，偏移可以将从地面接收到的来自地下的绕射波或反射波归位到实际位置[3?4]，常用的偏移成像技术包括频率?波数（[ω?k]）成像和后向投影（BP）成像。

[ω?k]成像基于波动方程以及“爆炸反射模型”的波场反向外推 [5?6]。该方法要求地下介质均匀，即电磁波在地下介质中的传播速度要保持恒定，因此适用于浅地层目标探测或非层状的均匀介质。[ω?k]算法的优势是采用FFT，成像效率高。因此，[ω?k]算法在实际工程中广泛使用。BP成像算法对地下分层介质的适应能力强、成像精度高，在探地雷达成像领域的应用也较多，但该方法计算量大、实时性差。本文基于子孔径的改进BP成像算法在保证成像质量的前提下，可大幅减少BP算法的计算量，提高成像效率。

1 标准BP成像算法

后向投影（Back?Projection，BP）起源于医学领域的CT技术，是McCorkle根据CT成像的投影切片理论推导出来的一种合成孔径雷达时域成像方法[7]。其原理是对成像区域内的任意一个待成像点先计算出该点距离各天线孔径位置的时延，然后在各道接收回波上搜索对应该时延点的幅值，最后将各道回波中相应位置处的幅值进行时域相干叠加，从而完成成像过程[8?9]。因此，BP成像算法的核心思想就是“时延?求和”。

建立如图1所示的探地雷达BP成像算法的场景模型。方位向沿[x]轴向右，距离向沿[z]轴向下，[k]个探地雷达收发天线与[x]轴平行，且与地面相距[h]，不同时刻的合成孔径天线位置用倒三角来表示，其中图示深色倒三角形表示当前发射和接收的第[i]个天线位置，其坐标记为[xi，-h]，浅色倒三角形为其余天线位置。成像场景被直线[z=0]分为两部分，直线以上为空气，直线以下假定为各向同性的均匀地下介质。由于空气与地下介质的介电常数不相等，电磁波在空气与地面分界处会发生折射，因此电磁波从空气到地下介质的实际传播路径为一条经过折射点[（xr，0）]的折线，而不再是经过点[（xl，0）]的直线。

令[θi]和[θr]分别为入射角和折射角，介质相对介电常数为[εr]，成像场景中待成像点[P]的坐标为[xP，zP]，由Snell折射定律，有如下关系式：

根据三角几何关系，有：

[sin θi=xr-xi（xr-xi）2+h2] （2）

[sin θr=xP-xr（xP-xi）2+z2P] （3）

联立以上三式，可得：

[（xr-xi）2（xr-xi）2+h2?（xP-xr）2+z2P（xP-xr）2=εr] （4）

式（4）为关于折射点[xr]的一元四次方程，直接求解过程繁琐。根据Mast[10]提出的方法，折射点[xr]的近似求解公式如下：

[xr≈x0+（xl-x0）εr] （5）

在折射点位置[xr]求出解后，可以求得各道回波上对应于待成像点[P]的时延为：

[τP，i=2（xr-xi）2+h2c+2（xP-xr）2+z2Pv] （6）

式中：[τP，i ]表示第[i]个合成孔径位置到点[P]的双程时延，[c]为电磁波在自由空间中的波速；[v=cεr]表示电磁波在介质中的波速。

若第[i]道回波数据可表示为[Si（t）]，便可得到[P]点在第[i]道回波上的散射响应幅值：

[xP， i=Si（t=τP， i）] （7）

将[P]点在各道回波数据上的散射响应幅值进行叠加，便完成对[P]点的成像：

[EP=xP， i] （8）

式中[EP]为待成像点[P]的最终成像结果。

通过遍历成像区域中全部待成像点，不断重复上述“时延?求和”操作，即可实现对整个探地雷达图像的BP成像。

2 改进的BP成像算法

假定待成像点[P]的坐标为[x，z]，由标准BP成像算法的过程可知，[P]点成像结果可表述为：

式中：[K]为总的天线数目；[Si（t）]为第[i]个天线孔径位置处的回波数据；[τx，z，i]为[P]点坐标[x，z]到第[i]个天线位置的双程时延。

由前述分析可知，标准BP成像算法存在一定的缺陷：一方面，在整个BP成像过程中，对于每个待成像点都需要遍历全部[K]个天线孔径位置，对大小为[M×N]的回波矩阵，算法的时间复杂度为[Θ（MNK）]，当探测区域较大时，算法运算速度慢、效率低；另一方面，对于成像区域内的某个待成像点，并不是在所有天线位置都会产生散射响应，这些冗余计算进一步增大了BP算法的运算量，因此散射响应小或者没有散射响应的位置在成像时可以不遍历。准确地判断埋地目标回波在接收的回波矩阵中的分布情况后，BP成像只利用目标回波矩阵进行运算，从而实现快速成像。基于此，改进后的BP成像算法将[K]个合成孔径天线按照一定的准则划分为[Nsub]个子孔径，这样式（9）中的求和操作无需遍历全部[K]个合成孔径位置，而是在天线孔径数目为[ Ksub]的天线子孔径内进行相干叠加，得到若干个子孔径图像，最后合并全部子孔径图像即完成整个成像过程。因此，改进的BP算法复杂度为[ΘMNKsub]，运算时间为标准BP算法的[KsubK]。

根据改进的BP成像算法，第[j]个子孔径对应的子孔径图像成像结果可表述为：

式中：[Ejx，z]为第[j]个子孔径图像中坐标点为[x，z]的最终成像结果；[Sj，it]为第[j]个子孔径内第[i]个天线孔径位置的回波数据；[Ksub]为每个子孔径内所需遍历的天线孔径数；[Nsub=KKsub]表示子孔径数。

3 算法性能对比

探地雷达成像算法的性能可以用方位向分辨率、积分旁瓣比（Integrated Side Lobe Ratio，ISLR））、运算时间等指标来衡量[11]。其中，方位向分辨率用归一化扫描能量图来表征，而ISLR是雷达信号处理中常用的一个指标，用来衡量算法对旁瓣和杂波的抑制程度，ISLR值越大则抑制效果越好，成像精度越高[12]。ISLR定义为目标所有旁瓣能量与主瓣能量之比，其表达式如下：

[ISLR=-10lgEtotal-EmainEmain] （11）

式中：[Etotal]为目标总能量；[Emain]为主瓣能量。

下面分别用仿真和实测数据进行实验，并定量地用上述性能指标加以衡量和对比。

3.1 仿真数据实验

GprMax2D为一款基于FDTD算法的探地雷达二维正演仿真软件，常用于探地雷达成像方法的研究[13]。图2为利用GprMax2D软件模拟“钢筋?混泥土”的场景并用不同方法进行成像。

仿真数据成像中各性能参数的大小见表1。

从表1可以看出：方位向分辨率上，标准BP算法和改进的BP算法最高；由ISLR可知，改进的BP算法ISLR值最大，表明其旁瓣和杂波抑制效果最好，成像精度最高，而[ω?k]成像的ISLR值最小，这在图像上表现为寄生能量较BP成像多；从运行时间上看，[ω?k ]算法的时间开销最少，标准BP成像算法的时间最长、效率最低；相对于标准BP成像算法而言，改进后的BP成像算法的運行时间仅为其1/5，不仅在时间消耗上大幅减少，而且在成像质量上亦有所提升，证明了改进算法的有效性。

3.2 實测数据实验

实测数据成像如图3所示，图3a）为国外某反雷技术研究中心实测的探地雷达原始B?Scan数据，图像大小为512×98，反映的是对埋藏在地下的PMN?2地雷进行探测的结果。

实测数据成像中各算法的性能参数大小见表2。

从表2可以看出：实测数据实验的处理结果与仿真实验的结论不谋而合，即三种成像算法中，[ω?k]偏移成像算法耗时最短，实时性最好；标准BP成像算法耗时最长；改进的BP算法成像精度最高，耗时较标准BP算法明显减少，但仍要大于[ω?k]算法。因此，当系统对成像算法的实时性要求较高时，应首选[ω?k]成像算法；而当注重成像质量时，可采用改进的BP成像算法。

4 结 语

[ω?k ]成像基于波动方程和“爆炸反射模型”，成像质量弱于BP算法，但实时性好。标准BP算法实现原理简单，成像精度高，但运算效率低，实时性较差。仿真和实测数据实验结果表明，本文改进的BP算法相对于标准BP算法，不仅成像质量有所提升，而且成像重建的运算速度加快，成像效率得到明显提升。

Fig. 3 Imaging results of actual measured data

参考文献

[1] TREES H L V. Ground penetrating radar theory and applications [M]. US： Elsevier Science， 2009.

[2] 张春城.浅地层探地雷达中的信号处理技术研究[D].成都：电子科技大学，2005.

ZHANG Chuncheng. Study on signal processing technology in ground penetrating radar of shallow layer [D]. Chengdu： UESTC， 2005.

[3] 张金花.车载探地雷达天线特性分析及其成像处理研究[D].成都：西南交通大学，2014.

ZHANG Jinhua. Analysis of antenna characteristic of vehicle ground penetrating radar and study on imaging processing [D]. Chengdu： SWJTY， 2014.

[4] 蔚建斌，陈自力，江涛.基于偏移技术的探地雷达SAR成像方法[J].信号处理，2010，26（5）：778?782.

WEI J B， CHEN Z L， JIANG T. The SAR imaging method of GPR based on migration [J]. Signal processing， 2010， 26（5）： 778?782.

[5] CHEN G Y， BUI T D. Multiwavelets denoising using neighbo?ring coefficients [J]. IEEE signal processing， 2003， 10（7）： 211?214.

[6] ENGIN E， CIFTCIOGLU B， ?ZCAN M， et al. High resolution ultrawideband wall penetrating radar [J]. Microwave & optical technology letters， 2007， 49（2）： 320?325.

[7] 郑文军.复杂多散射环境下TRM成像技术研究[D].成都：电子科技大学，2010.

ZHENG W J. Study on TRM imaging technology under complex multi?scattering environment [D]. Chengdu： UESTC， 2010.

[8] 雷文太.脉冲GPR高分辨率成像算法研究[D].长沙：国防科技大学，2006.

LEI W T. Study on pulse GPR high?resolution imaging algorithm [D]. Changsha： NUDT， 2006.

[9] 孔令讲.浅地层探地雷达信号处理算法的研究[D].成都：电子科技大学，2003.

KONG L J. Study on signal processing algorithm of ground penetrating radar of shallow layer [D]. Chengdu： UESTC， 2003.

[10] MAST J E， JOHANSSONE M. Three?dimensional ground pe?netrating radar imaging using synthetic aperture time?domain focusing [J]. Proceedings of SPIE： the international society for optical engineering， 1994， 2275： 205?214.

[11] OZDEMIR C， DEMIRCI S， YIGIT E， et al. A review on migration methods in B?scan ground penetrating radar imaging [J]. Mathematical problems in engineering， 2014 （1）： 1?16.

[12] YIGIT E， DEMIRCI S， OZDEMIR C， et al. Short?range ground?based synthetic aperture radar imaging： performance comparison between frequency?wavenumber migration and back?projection algorithms [J]. Journal of applied remote sen?sing， 2013， 7（1）： 382?385.

[13] 周奇才，李炳杰，郑宇轩，等.基于GPRMax2D的探地雷达图像正演模拟[J].工程地球物理学报，2008（4）：396?399.

ZHOU Q C， LI B J， ZHENG Y X， et al. Forward simulation of GPR image based on GPRMax2D [J]. Chinese journal of engineering geophysics， 2008（4）： 396?399.

[14] 李广场.有限差分法探地雷达波动方程偏移[D].南京：河海大学，2004.

LI Guangchang. Wave equation offset of ground penetrating radar with finite difference method [D]. Nanjing： HHU， 2004.