基于GprMax软件的道路路基空洞探地雷达正演模拟

尹光辉， 冯雨宁， 张怀凯， 冯兴乐

(1.长安大学　信息工程学院,西安　710064；2.北京大学　地球与空间科学学院,北京　100871)



基于GprMax软件的道路路基空洞探地雷达正演模拟

尹光辉1， 冯雨宁2， 张怀凯1， 冯兴乐1

(1.长安大学信息工程学院,西安710064；2.北京大学地球与空间科学学院,北京100871)

探地雷达具有高效、无损、分辨率高、使用灵活的特点，在道路检测方面的应用越来越广泛。根据时域有限差分(FDTD)算法，建立理想条件下道路路基充气和充水空洞模型，用GprMax及Matlab编程分别对不同形状、不同激励源频率、不同大小条件下空洞模型进行正演模拟，并对空洞模型正演模拟雷达图像进行分析。模拟结果有助于对实际探测空洞识别和雷达图像的解释，为探地雷达在路基空洞检测应用提供参考。

探地雷达； 空洞； 正演模拟； FDTD算法； GprMax

0　引言

近年来，我国的道路空洞形成的塌陷事故愈演愈烈，给人民的生命财产安全造成了极大危害。与传统钻芯取样道路检测方法相比，由于探地雷达具有高效、无损、分辨率高的特点，且道路结构层及土基和下方潜伏的充气空洞、充水空洞之间的介电常数差异大，因此探地雷达技术适合应用在道路空洞的检测。在实际应用过程中，缺乏探地雷达空洞异常图像评判客观标准，道路空洞与探地雷达异常图像的对应关系缺乏系统地研究，且雷达资料解释具有多解性和主观经验性，导致同一个异常图像不同的人得出不同的解释结果，导致探地雷达技术在道路地下空洞探测中的应用难以获得理想效果[1]。因此，利用计算机开展道路路基空洞正演模拟，研究各种空洞模型的探地雷达探测图像特征和规律，不必花费大量人力、物力和时间进行实地探测，就可以得到接近实地探测的效果。这对于认识探地雷达实测结果，识别空洞具有实际意义[2]。作者通过建立道路路基下充气与充水空洞模型，用GprMax及Matlab编程实现道路路基空洞正演模拟，对空洞模型正演模拟图像形态和特征进行分析，并通过实例验证正演模拟的有效性。这里模拟结果有助于对实际探测空洞识别和探地雷达图像地解释，提高检测结果的可靠性和检测评判标准的统一性，为探地雷达在路基空洞检测应用积累理论经验。

1　时域有限差分法正演的基本理论

探地雷达探测路基空洞中高频电磁波以宽频带短脉冲形式，通过发射天线被定向送入地下，经存在电性差异的地下地层或空洞反射后返回地面，由接收天线所接收。高频电磁波在介质中传播时，其传播路径、电磁场强度与波形随待测介质的电性特征与几何形态而变化[4]。因此，通过采集时域波形、并进行处理和分析，可确定地下分界面或地质体的空间位置及结构。

为了模拟探地雷达探测时的电磁波在地下介质层及空洞中传播，K.S.Yee提出一种计算电磁波数值的时域有限差分法，基于Yee[3]网格空间离散方式，将地下待探测区域划分为由一定数量的空间网格组成的仿真空间。在FDTD离散中，Yee网格空间中电场和磁场各节点的空间排布见图1。每个磁场分量由四个电场分量环绕，每个电场分量由四个磁场分量环绕。然后用有限差分式代替Maxwell方程时域场旋度方程中的微分式，在设置初始的场值和边界条件后，通过逐步递进的方法求得各个时刻的空间电磁波电磁场的分布数值，就获得了电磁波在整个模型仿真空间的传播模拟结果。

图1　FDTD算法空间Yee网格Fig.1　Yee space grid of FDTD algorithm

由于FDTD算法是用有限差分式代替Maxwell方程组，离散后的差分方程组的解必须收敛稳定。在FDTD中，时间步长Δt和空间网格步长Δx、Δy、Δz不是相互独立的，它们的取值满足一定的关系，才能避免数值结果的不稳定[5-8]，解得稳定性条件如式(1)所示。

(1)

用差分法对麦克斯韦方程进行数值计算将引起波的色散，即在FDTD网格中波的传播速度将随波长而改变。这种色散将导致非物理原因引起的脉冲波形畸变，使得时域数值计算产生误差。因此，空间网格的几何长度要满足数值色散条件：

(2)

2　GprMax正演模拟

GprMax是基于FDTD算法的探地雷达正演模拟软件，可以用于各向同性均匀介质中电磁波的传播规律以及电磁波与目标体的相互作用，通过编写正演模拟模型的输入文件程序，保存后缀为*.in的文件，输入程序文件路径进行仿真。计算机进行数值计算后得到关于目标体的探地雷达正演模拟后缀为*.out和*.geo的二进制文件。通过编写Matlab程序读入二进制文件数据，处理后缀*.geo文件获得模型结构图，处理后缀*.out文件生成探地雷达正演模拟图像(图2)。

图2　GprMax正演模拟流程图Fig.2　Flow of GprMax forward simulation

2.1道路结构模型参数

为了建立地下待探测区域的物理模型，需要设置相对介电常数、电导率参数。根据实际道路结构，把道路模型定为三层：第一层是道路面层，由厚度为10 cm的沥青材料构成，沥青材料相对介电常数=4，电导率σ=0.005 S/m；第二层为道路基层，由厚度为15 cm混凝土材料构成，相对介电常数εr=9，电导率σ=0.05 S/m；第三层为土基层，由厚度为20 cm压实处理后的土壤构成并向下延伸到2 m，相对介电常数=12，电导率σ=0.1 S/m。水的相对介电常数εr=81，电导率σ=0.03 S/m。

3　道路路基空洞模型正演模拟

3.1不同形状充气及充水空洞正演模拟

根据道路结构模型，分别构造位于土基层的矩形、圆形空洞(图3)，空洞中心点深度距离地表为0.6 m。矩形空洞左下顶点坐标(1,1.2)，右上顶点坐标(1.5,1.6); 圆形空洞圆心坐标(1.25,1.4)，直径为0.4 m;结构激励源为300 MHz的ricker子波，其他参数如表1所示。正演模拟结果如图4所示。

图3　不同形状类型空洞结构图Fig.3　Structure of cavity in different shapes(a)矩形空洞; (b)圆形空洞

表1　FDTD正演模拟参数

对比分析不同形状空洞正演模拟图像，矩形、圆形充气空洞正演模拟结果均呈现双曲线特征，曲线开口向下。双曲线顶部反射振幅最强，两端振幅较弱。地下空洞在地面的水平中心投影位置可以由双曲线顶点确定。矩形空洞由于顶面水平方向具有较长延伸，顶点反射较为明显，双曲线有水平弧度。总的来说圆形、矩形正演模拟图像形成的双曲线较为近似，很难通过双曲线形状估计地下空洞的大致形状。

提取矩形充气及充水空洞中心测线电磁波幅值数据得到电磁波波形图(图5、图6)。发射天线电磁波从空气射入面层、基层、路基时，介电常数依次增大，在每层分界面电磁波反射系数都为负值，因此在0～1 500测点数电磁波发生四次反向，充水型矩形空洞介电常数为81，远大于路基层介电常数，电磁波反射系数为负值，中心测线处电磁波进入充水型空洞发生反向，充气型空洞反射系数为正值，电磁波不发生反向。对比相同形状但是内部填充介质不同的充水空洞，充水类型空洞正演模拟结果形成两条双曲线波形。以圆形充水空洞为例，电磁波在水中的传播速度0.033 m/ns，传播时长约为25 ns，两者相乘除以2即为空洞的直径0.412 5 m与模型参数设置的0.4 m接近。因此，实际探地雷达探测过程中，可以根据双曲线的条数确定空洞内部填充介质是空气还是水，充水型空洞可以根据两条双曲线的间隔时间，计算空洞的垂直距离。

3.2不同激励源中心频率正演模拟

在不同中心频率正演模拟中，道路及空洞模型与图2的圆形空洞模型一致，改变ricker子激励源中心频率，依次为200 MHz、300 MHz、400 MHz、500 MHz、600 MHz，研究不同中心频率激励源条件下雷达波成像规律，不同频率下正演模拟图像如图7所示。

对比不同频率正演模拟图像，200 MHz图像分辨率最低，随着激励源中心频率的增大，雷达波正演模拟图像分辨率逐渐提升。其中以400 MHz和600 MHz的雷达波正演模拟图像分辨率最好，但不是频率越高正演模拟图像分辨率越好，500 MHz激励源正演模拟图像中出现许多杂波。

3.3不同大小圆形空洞正演模拟

不同大小圆形空洞顶面埋深一致，半径分别为0.2 m、0.4 m、0.6 m，ricker激励源频率为400 MHz，其他参数与前面一致。

图8中的正演模拟结果表明，随着圆形空洞半径的增大，正演模拟图像生成的双曲线雷达波图像长度增大，曲率半径增大。但是双曲线延长线水平位置坐标与实际空洞位置并不符合，不能通过双曲线长度获得空洞水平距离。

图4　不同填充介质正演模拟图Fig.4　Simulation of cavity in different medium(a)矩形充气空洞; (b)矩形充水空洞; (c)圆形充气空洞; (d)圆形充水空洞

图5　矩形充气空洞中心测线波形图Fig.5　Center line waveform of rectangle aeration cavity

图6　矩形充水空洞中心测线波形图Fig.6　Center line waveform of rectangle watery hole

图7　不同频率ricker激励源正演模拟图Fig.7　Simulation of different frequencies ricker excitation(a)频率200 MHz；(b)频率300 MHz；(c)频率400 MHz；(d)频率500 MHz；(e)频率600 MHz

图8　不同半径圆形充气空洞正演模拟图Fig.8　Simulation of aeration circular cavity aeration in different radius(a)半径0.2 m；(b)半径0.4 m；(c)半径0.6 m

4　探测实例正演模拟验证

利用天毅达无线探地雷达对某处沥青路面进行检测，探测结果如图9所示。在雷达增益为110倍时，实时伪彩图显示清晰双曲线，与正演模拟结果呈现一致性，初步推断为空洞。经过探地雷达数据处理，显示地下埋深约为0.5 m，实际验证为一竖直排水口。根据探地雷达实测资料以及参数设置对此排水口进行正演模拟，构造0.6 m×1.3 m×0.85 m三维空间，空间网格大小为0.01 m×0.01 m×0.01 m，路面介质材料为沥青，介电常数为5，电导率为0.01 S/m。构造埋深约为0.5 m的矩形空洞，矩形空洞三维坐标(0.2, 0.55, 0.1)(0.4, 0.75, 0.3)，收发天线坐标(0.3, 0.125, 0.55)(0.3, 0.375, 0.55)，天线沿Y轴步进距离为0.04 m，测线道数21条，ricker激励源发射频率为900 MHz，时窗15 ns。正演模拟图像如图10所示。对比正演模拟图像与实测图像，路面结构层分层时间都在3 ns，双曲线顶点时间都在6 ns，图像吻合度较好，验证了基于时域有限差分法的GprMax结合Matlab进行探地雷达探测空洞正演模拟结果的有效性。

图9　探地雷达实测图像Fig.9　Image of GPR detection

图10　实例正演模拟图像Fig.10　Simulation of instance cavity

5　结论

基于时域有限差分法，通过GprMax及Matlab编程实现对道路路基空洞模型的正演模拟，模拟结果表明：

1)不同形状空洞正演模拟结果双曲线类似，难以通过双曲线反推空洞形状。可根据双曲线顶点计算空洞顶面埋深，根据两条双曲线相差走时计算充水型空洞垂直距离，根据双曲线条数推断空洞充水还是充气类型。

2)不同激励源中心频率，影响探地雷达图像分辨率，一定频率范围内中心频率越高，探地雷达图像分辨率越高。

3)不同大小规模空洞形成的双曲线曲率不同、延长度不同，空洞规模越大，双曲线延长度越大，曲率越大。

通过探地雷达实测图像与正演模拟结果对比，验证了正演模拟结果的有效性，可用于指导探地雷达空洞探测实例识别，为建立探地雷达空洞模型数据库提供数据资料。

[1]葛广志.北京市道路地下空洞探地雷达探测正演模拟研究[D].北京:中国地质大学,2010.

GE G Z.Forward simulation of underground cavity in Beijing by GPR[D].Beijing:China University of Geosciences,2010.(In Chinese)

[2]郭成超,王复明.探地雷达电磁波正演模拟研究[J].公路交通科技,2008,25(8):37-41.

GUO C C,WANG F M.Study on forward modeling of GPR electromagnetic wave propagation[J].Journal of Highway and Transportation Research and Development,2008,25(8):37-41.(In Chinese)

[3]彭湘佳.探地雷达在北京市道路下方病害检测中的研究[D].北京:中国地质大学,2011.

PENG X J.Research on the disease below road in Beijing by GPR[D].Beijing:China University of Geosciences,2010.(In Chinese)

[4]冯彦谦,王银,刘四新.基于FDTD的探地雷达数值仿真与成像研究[J].铁道工程学报,2009,26(6):17-20.

FENG Y Q,WANG Y,LIU S X.Research on the numerical simulation and imaging of GPR based on the FDTD[J].Journal of Railway Engineering Society,2009,26(6):17-20.(In Chinese)

[5]周奇才,李炳杰,郑宇轩,等.基于GPRMax2D的探地雷达图像正演模拟[J].工程地球物理学报,2008,5(4):396-399.

ZHOU Q C,LI Y J,ZHENG Y X,et al.Forward simulation of GPR image based on GPRMax2D[J].Chinese Journal of Engineering Geophysics,2008,5(4):396-399.(In Chinese)

[6]谢雄耀,于超,赵永辉,等.山岭隧道空洞灾害雷达探测图像影响因素的正演模拟研究[J].工程地质学报,2010,18(1):81-87.

XIE X Y,YU C, ZHAO Y H,et al.The forward simulation research of cavity radar image influence factors in mountain tunnel[J].Journal of Engineering Geology,2010,18(1):81-87.(In Chinese)

[7]徐波,蒙爱军,吴凤昌.利用数值模拟仿真探地雷达地下障碍物探测[J].路基工程,2010,27(5):170-171.

XU B,MENG A J,WU F C.Underground obstacle detection of GPR by numerical simulation[J].Subgrade Engineering,2010,27(5):170-171.(In Chinese)

[8]A.GIANNOPOULOS.Modeling ground penetrating radar by GprMax[J].Construction and Building Materials,2005,19(10):755-762.

Forward simulation of ground penetration radar based on the GprMax for the roadbed cavity

YIN Guang-hui1, FENG Yu-ning2, ZHANG Huai-kai, FENG Xing-le

(1.School of Information Engineering, Chang 'an University, Xi'an710064, China;(2.School of Earth and Space Sciences, Peking University, Beijing100871,China)

Ground penetrating radar has been widely used in road detection due to the features of efficient, non-destructive, high resolution and flexible. Based on the finite different time domain algorithm (FDTD), the ideal models of roadbed cavity filled with air and water are built, GprMax and Matlab programming are used to realize forward simulation which aimed at different shapes, different center frequencies and different sizes, and then the cavity modeling simulation radar images are analyzed. Simulation results contribute to the detection of cavity and the interpretation of the radar image and references are provided for the ground penetrating radar application in cavity detection.

ground penetration radar; cavity; forward simulation; FDTD algorithm; GprMax

2015-06-15改回日期：2015-09-05

国家自然科学基金(41404095)

尹光辉(1990-)，男，硕士，主要研究方向为探地雷达检测，E-mail:412245505@qq.com。

1001-1749(2016)04-0480-07

P 631.3

A

10.3969/j.issn.1001-1749.2016.04.07