基于椭圆包络线算法的探地雷达成像

邵泉杰，欧阳缮

(桂林电子科技大学 信息与通信学院，广西 桂林541004)

探地雷达(ground penetrating radar，简称GPR)是一种地下目标无损检测技术，具有探测速度快、分辨率高、操作灵活方便、探测范围广及探测费用低等优点，被广泛应用于公路检测［1］、桥梁的无损检测［2］、地质工程探测和考古探测等［3］。GPR的基本原理为:发射天线在地面以宽频带窄脉冲的形式向地下发射高频电磁波，入射波在介电性差异的2种介质分界面(地层界面或目的体)产生反射波，接收天线接收记录反射波的波形、振幅及到达时刻，并以雷达图象方式显示探测结果;根据测量的双程走时和波速计算出目标体深度，连续测量剖面各点的反射波，形成雷达图像。对探地雷达回波数据进行处理，最终目的是解释地下异常的情况，其中目标成像是最直观有效的数据处理方式。已有的成像算法研究主要包括反向投影(back projection，简称BP)算法［4］、时间反转(time reversal，简称TR)成像算法［5］、压缩感知(compressive sensing，简称CS)算法［6-7］、反演(inversion)成像［8］等。其中BP算法最为经典，该算法对成像域进行像素点划分，通过计算回波时延进行能量累加成像，实现简单，但运算量大。TR算法充分利用电磁波传播的物理过程，适合复杂多散射环境，具有较高的成像分辨率。反演成像一般基于Maxwell方程，能计算出目标的位置、性质和尺寸等，但其过程需迭代等运算，复杂度高。上述方法大都计算量大，且不适合实时处理。CS算法虽可降低数据存储和运算复杂度，以较少数据对信号进行重构，但对测量矩阵和重构算法要求严格。

椭圆包络线算法［9］是一种目标边界成像算法，该算法基于收发分置天线提出，利用天线与目标边界所构成的几何相切关系成像，实时性好，计算效率高。穿墙雷达只需考虑墙体的影响，对墙体进行补偿后，电磁波相当于在单一的空气介质中传播。但探地雷达对目标体进行探测过程中，由于地下介质的复杂度高，地下情况未知，可能存在多层介质，且多层介质的形状不像墙体那样水平，补偿的方法明显不适用。为此，以一种近似构造椭圆的方法将包络线算法扩展到探地雷达地下目标成像，使其适用于多层介质的地下环境。近似后的椭圆包络线算法对多层介质地下目标成像深度有微小的误差，通过分析误差影响因素，可发现算法的最佳适用环境。

1 椭圆包络线目标成像算法

1.1 穿墙雷达中的椭圆包络线算法

包络线算法用于穿墙雷达时，墙体是对成像产生误差的主要外界因素。由于墙体比较均匀，厚度一定，可采用补偿方法消除墙体影响。基于收发分置的椭圆包络线算法用于穿墙雷达的几何原理如图1所示，图中墙体已进行补偿。收发天线沿着侧线即x轴自左向右移动，目标物体位于空气介质中，假设在一次天线位置，发射天线坐标为F1(XT，0)，接收天线坐标为F2(XR，0)，电磁波照射目标体的C(x，y)点，天线与目标体之间的回波往返距离为Y=Y1+Y2。

图1 穿墙雷达几何原理图Fig.1 Geometry relationship of through-wall radar

几何学上椭圆的定义为:平面上到两点距离之和为常数的点的轨迹。为此，在固定天线位置处，收发天线坐标视为2个定点，电磁波在空气中传播的距离为常数，便可得天线、回波路径以及目标边界之间的椭圆关系。将发送天线与接收天线坐标看成椭圆的2个焦点坐标F1(XT，0)和F2(XR，0)，则半焦距为。天线与目标点之间的回波路径长度Y=c t=2a。其中:c为电磁波在空气中的传播速度;t为电磁波从发射天线到达目标边界点后返回接收天线所用的时间。由此计算椭圆长轴a=v t/2，根据椭圆的性质，可得短轴。如此，回波路径与天线坐标可建立几何关系，得到椭圆方程表达式:

收发天线均以某一固定间隔距离ΔX沿横轴自左向右同时移动，发送天线坐标集合XT={X1，X2，…，XN}，对应的接收天线坐标集合XR={X1－ΔX，X2－ΔX，…，XN－ΔX}。每个天线位置对应一个椭圆，并与目标边界的一个点相切。众多椭圆包络簇实现回波双曲线的能量聚焦，基本无误差地实现对目标上边界的还原。

1.2 探地雷达中的椭圆包络线算法

探地雷达针对地下目标体进行探测，假设目标体为圆形，收发分置天线位于地表，电磁波进入地下环境介质中传播，设地下介电常数为μ1，目标体的介电常数为μ2，目标体的上下界面均为曲面，电磁波穿透目标体，上下表面均会产生回波。上表面的回波为单一介质情况，与墙体补偿后的穿墙雷达成像相似，唯一区别为介质中的波速。下表面的回波经过2种介质，即地下环境介质及目标体介质，而且2种介质的分界面均为曲面，这种情况下对下表面进行成像，即多层介质成像。构造椭圆时需要考虑2种介质的情况，其几何关系如图2所示。对于多层介质情况，椭圆焦点坐标仍然为发射天线和接收天线的坐标。由于接收到的回波数据为电磁波的时延，长轴的计算采用不同介质内电磁波运行距离分段相加近似计算的方法，

其中:Y1+Y4为电磁波在介质环境中往返路程;Y2+Y3为目标体内往返路程;t1为目标上表面回波往返时间;t2为目标下表面的回波往返时间;v1为电磁波在环境介质中的传播速度;v2为电磁波在目标介质中的传播速度。椭圆长轴a=Y/2，通过b=计算短轴，构造椭圆包络线，从而实现下界面的边界重建。对于实际中的天线离地较远情况，相当于增加一层空气介质，分段相加的方法仍适用。

图2 地下多层介质几何关系Fig.2 Geometric relationship of underground multi-medias

所有天线位置构造的包络线簇呈现目标体的边界形状。为了更直观地观察，提取反映目标区域的边界，这需确定边界的起始点和终止点以及2个端点之间的值。假设天线移动过程中总计产生M个椭圆，天线离目标体越远，电磁波照射到的目标体的位置越边缘，故以第1个和第2个椭圆的交点横坐标为起始点xmin，第M－1个和第M个椭圆的交点为终止点xmax。目标边界有凹和凸2种情况，与椭圆的对应关系分别为内切和外切，图2所示的目标体下界面为凹形，构造的椭圆与下界面外切，按照文献［9］确定形状与椭圆的相切关系。在划定成像区域范围［xmin，xmax］且y≥0中，椭圆的最小个数为N，椭圆与坐标轴围成的椭圆区域集合为∂E，即∂E=。∂E1与∂EN交点为P，目标边界为∂B。设∂Eα为∂E边界点的集合，∂Eβ为∂E内部点的集合，满足并集关系∂E=∂Eα∪∂Eβ。根据P所在的集合建立边界提取的数学描述，其表达式为:

其中:∂Eout和∂Ein分别为椭圆与目标外切和内切时的成像边界集合;ai、bi、ci和xi分别为椭圆∂Ei的半长轴、半短轴、半焦距和发射天线横坐标。

1.3 算法性能及可行性分析

穿透多层介质时，以一种分段相加近似计算的关系构造椭圆长轴。由Snell折射定律知，电磁波在介电特性不同的2种介质中，入射角和折射角的大小与2种介质的介电特性有关。例如，图2中以Y1+Y2+Y3+Y4近似计算长轴，而实际长轴计算应当为2点间直线距离F1C+CF2，由三角形关系知，三角形任意二边之和大于第三边，所以穿透多层介质时计算的深度略微有所增加。

误差产生的原因主要是多层介质介电特性的差异，反映到距离上，即为电磁波在不同介质中的运行距离。与电磁波运行距离最直接相关的因素是目标体埋深。验证目标体埋深对误差影响，设置不同深度的模型，如图3所示。电磁波从介质1入射，经过交界面进入介质2，介质2内存在目标体，假设介质1和介质2的介电常数不变，目标体深度不同。实际中的长轴为x1，近似构造的长轴为a1+a2，由三角形性质可得如下关系:

在介电常数一定情况下，角度θ是固定的，误差Δ只与天线在介质1的高度以及目标体的埋深有关。深度a1越小，天线离地距离a2越小，误差越小。

图3 不同深度模型Fig.3 Different depth model

2 算法仿真与实测数据分析

2.1 仿真结果与分析

采用GprMaxV2.0进行仿真，验证椭圆包络线算法对目标边界成像的有效性，并探讨探地雷达成像的性能以及深度对多层介质成像的影响。仿真场景为:地下区域长2 m、宽1.2 m，环境的介电常数为3，地下埋一个圆柱体，圆心位于地下0.5 m，半径分别设置为0.2、0.1 m，介电常数为10;天线采用收发分置间距0.1 m，每次移动0.03 m，总计移动60步，频率900 MHz，2种介质的介电特性已知，求得环境介质和目标体内电磁波传播速度，对圆柱上下界面进行包络线成像，仿真结果如图4所示。

图4 包络线成像的仿真结果Fig.4 Simulation results of envelope imaging

由仿真结果可知，上表面只涉及单一介质，在已知介质波速情况下，不管深度还是形状，其还原均十分准确。由于分段相加近似计算的长轴偏长，在半径为0.2 m情况下，圆的下表面成像产生了接近0.02 m的深度误差。但相同仿真条件下，半径为0.1 m的圆下表面的误差基本可忽略。上述结果验证了深度越小，近似计算的包络线算法误差越小。可见，对于浅层目标成像，近似计算的包络线算法具有较好的性能。

2.2 实测数据成像与分析

实测数据实验场景描述:沙坑中埋有一个圆弧形的铁片，长度为0.18 m，每次移动0.02 m，总计移动41步。使用的实验设备为脉冲源、喇叭天线、数据采集示波器、低噪声放大器等。脉冲源为二阶高斯脉冲，6 db带宽0.32～1.56 GHz;发射接收天线型号ATB1G-12-003，均为极化天线;低噪声放大器型号ADL5523;示波器型号安捷伦ium 80000，采样频率40 GHz，数据经示波器显示并保存。

波速通过已知目标深度换算的方法［10］精确求取。提取回波延时，根据天线移动情况，应用多层介质情况下的包络线成像方法，成像结果如图5所示。由于物体埋深很浅，原始数据采集较为准确，与物体上表面的实线进行对比，成像后误差较小，上表面的形状还原精确。

图5 包络线成像结果Fig.5 Envelope imaging result

与CBP成像算法［11］比较，原始数据经CBP成像后的结果如图6所示。从图6可见，改进的BP算法只能粗略地对能量进行聚焦，虽然物体的大致位置信息等基本正确，但根据成像结果很难分辨出物体的形状，图中目标物体位置出现2个类似目标体的情况，这会导致地下情况的错误解读。BP成像算法需要对众多像素点进行能量叠加，导致计算复杂度较高。在奔腾E5400处理器2 GB内存情况下，Matlab进行BP成像所需时间为49.54 s，包络线成像所需时间为2.58 s。可见，包络线算法的处理速度明显优于BP算法，更适合实时处理。

图6 CBP成像结果Fig.6 CBP imaging result

3 结束语

将包络线成像算法应用于探地雷达地下目标成像，根据多层介质内电磁波照射物体往返路径分段相加的近似方法构造椭圆的长轴，通过理论分析及实测数据验证浅层目标成像，近似相加方法的误差在可允许范围之内。与BP成像算法进行对比，包络线算法对边界的构造准确，运算速度快，更适合实时处理。

［1］卢成明，秦臻，朱海龙，等.探地雷达检测公路结构层隐含裂缝实用方法研究［J］.地球物理学报，2007，50(5):1558-1568.

［2］Wang Z W，Zhou M，Slabaugh G G，et al.Automatic detection of bridge deck condition from ground penetrating radar images［J］.IEEE Transactions on Automation Science and Engineering，2011，8(3):633-640.

［3］曾昭发，刘四新，王者江.探地雷达方法原理及应用［M］.北京:科学出版社，2006:213-223.

［4］Zhou L，Huang C，Su Y.A fast back-projection algorithm based on cross correlation for GPR imaging［J］.Geoscience and Remote Sensing Letters，2012，9(2):228-232.

［5］郑文军，赵志钦，张薇，等.超宽带探地雷达中TRM-SAR成像技术研究［J］.电子科技大学学报，2011，40(3):363-366.

［6］Tuncer M A C，Gurbuz A C.Ground reflection removal in compressive sensing ground penetrating radars［J］.Geoscience and Remote Sensing Letters，2012，9(1):23-27.

［7］Gurbuz A C，Mcclellan J H，JR Scott W R.Compressive sensing of underground structures using GPR［J］.Digital Signal Processing，2012，22(1):66-73.

［8］丁亮，韩波，刘润泽，等.基于探地雷达的混凝土无损检测反演成像方法［J］.地球物理学报，2012，55(1):317-326.

［9］李育晖，欧阳缮，晋良念，等.超宽带穿墙雷达椭圆包络线目标边界成像算法［J］.电子与信息学报，2014，36(7):1532-1537.

［10］邓小燕，王通.探地雷达探测中对媒质相对介电常数的测定［J］.物探与化探，2009，33(1):43-45.

［11］周琳，粟毅.基于互相关的探地雷达反向投影成像算法［J］.电子与信息学报，2011，33(11):2714-2719.