0 引 言

探地雷达(Ground penetrating radar，GPR)是利用高频宽带电磁波来探测地下结构和特性或者确定不可视物体的结构和特性的一种地球物理方法(曾昭发等，2010)，它具有效率高、分辨率高、无损性、抗干扰能力强、探测结果直观等特点，在军事、交通、环境、考古、水利水电工程等众多领域(Arias et al., 2007；冯德山和戴前伟，2008；Lubowiecka et al., 2009)中得到了广泛的应用.

探地雷达的正演模拟技术是探地雷达技术的基础，它一直是探地雷达科学研究的热点内容.目前，许多学者已经运用很多种不同的方法开展了探地雷达正演模拟研究，主要包括：时域有限差分法(Finite Difference Time Domain Method，FDTD)(Shaari et al., 2010； Zhu et al., 2013)、有限单元法(冯德山等，2012)、时域多分辨算法(冯德山等，2007)、交替方向隐式时域有限差分法(Feng and Dai， 2011)、错格时域伪谱法(李展辉等，2009)、辛分块龙格库方法(Fang and Lin， 2012)和无单元法(冯德山等，2013).其中，时域有限差分法以其存储空间小、计算效率高、适用性广泛和结果简单直观等特点，现已成为最为重要的探地雷达正演模拟方法之一.时域有限差分法的雏形发源于Yee的论文“Numerical Solution of Initial Boundary Value Problems Involving Maxwell’s Equation in Isotripic Media”，经过近二十年的发展，才逐渐走向成熟，其中Taflove和GoodMan D等人做出了重要贡献(Taflove，1992；GoodMan，1994).近年来，国内方广有等做了在满足德拜关系的频散介质中的探地雷达2.5维FDTD正演模拟，分析了影响探地雷达探测深度的因素(方广有等，1998).戴前伟等由地质雷达正演原理着手，做了加载超吸收边界条件的探地雷达时域有限差分法正演模拟(戴前伟等，2004).刘四新、曾昭发在采用Sullivan理论的基础之上，做了在三维频散介质中的探地雷达FDTD模拟(刘四新等，2006；刘四新和曾昭发，2007).刘磊、昌彦君等基于单轴各向异性理想匹配层(UPML)，推导出色散介质的二维FDTD迭代公式，详细分析了地下土壤分别为均匀介质与色散介质时，地下管线的电磁波响应特征(刘磊等，2009).冯德山等提出基于UPML边界条件的ADI-FDTD三维探地雷达数值模拟算法，运用MATLAB平台开发了ADI-FDTD三维模拟程序，并对三维复杂GPR模型进行正演模拟(冯德山和谢源，2011).本文作者运用基于UPML边界条件的FDTD二维GPR数值模拟程序开展正演模拟研究，并将该正演结果应用到水电站引水隧洞衬砌检测中，解决了工程实测数据解释中遇到的问题，为实测资料的合理解释提供了科学的理论依据.

1 时域有限差分法正演原理

根据电磁波理论可知，广谱、高频宽带电磁波在真空和媒质中的传播规律是满足麦克斯韦(Maxwell)方程组的，Maxwell方程组描述了当场随时间变化时互相耦合的磁场与电场的变化关系.在无源区域，Maxwell方程组的两个旋度方程可以表示为：

在直角坐标系中，令 E=E _x a _x+ E _y a _y+ E _z a _z，H=H _x a _x+ H _y a _y+ H _z a _z(a _x，a _y和a _z分别为x，y和z三个坐标的单位矢量)，则Maxwell旋度方程可展开为

应用K.S.Yee式网格模型，通过差分法把旋度方程中的微分格式近似地转变成差分格式，推导出二维TM波时域有限差分法探地雷达正演模拟方程为

探地雷达正演模拟是计算电磁场的辐射、散射等开放域问题，由于受到计算机资源的限制，探地雷达正演模拟不可能模拟无穷大的空间，所以只能用有限的空间模拟半无穷大的空间.因此，在模型边界的网格截断处一定要设置适当的边界条件，否则模型边界上一定会产生很强的边界反射.边界条件的好坏将直接影响正演模拟的精确度.好的边界条件要求反射小、稳定且不增加很大的计算量.具有宽频带吸收特性的PML边界条件是目前最适合FDTD法的吸收边界条件，它主要包括：Berenger PML、UPML和CPML等.UPML具有节约计算机存储量，高精度，高效率，无需进行电场和磁场分裂，迭代公式简单等优点，本文作者选用UPML边界条件结合FDTD法进行GPR正演模拟.

根据电磁场与电磁波理论可知：二维GPR波的UPML有4个平面区和4个棱角区，无源时，Maxwell方程可以表示为

如图 1所示，模型1是一个2.0 m×0.5 m的矩形区域，模型上部是空气层，厚度为0.05 m；下部为衬砌混凝土，其相对介电常数为8，电导率σ=0.005 s/m，混凝土层中存在2根不同型号的钢筋，电导率σ=0.01 s/m，位于坐标(0.75 m，0.3 m)处为Φ12钢筋，位于坐标(1.35 m，0.3 m)处为Φ18钢筋，Ricker子波主频率为900 MHz.

图 1

Fig. 1

图 1 模型1示意图 Fig. 1 Sketch map of model 1

如图 2为模型1正演剖面图，图中可见：在钢筋所在位置异常表现明显，存在2个典型的双曲线绕射波，双曲线弧顶为钢筋所在位置，2个双曲线绕射波不存在明显差异，在剖面图中未能区分出Φ12钢筋与Φ18钢筋，这是由于钢筋型号差别不大.

图 2

Fig. 2

图 2 模型1正演剖面图 Fig. 2 Simulate profile map of model 1

3.2 模型2

如图 3所示，模型2是一个2.0 m×0.5 m的矩形区域，模型上部是空气层，厚度为0.05 m；下部为衬砌混凝土，其相对介电常数为8，电导率σ=0.005 s/m，混凝土层中存在2个空洞，空洞相对介电常数为1，电导率σ=0.01 s/m，矩形空洞大小为0.2 m×0.1 m，三角形空洞底为0.2 m，高为0.1 m，Ricker子波主频率为900 MHz.

图 3

Fig. 3

图 3 模型2示意图 Fig. 3 Sketch map of model 2

如图 4为模型2正演剖面图，图中可见：在空洞所在位置异常表现明显，可以清晰地分辨出空洞异常体的上、下界面反射波.空洞下界面反射波能量比上界面弱，上界面在剖面图中依然是平界面，上界面的两个棱角位置出现双曲线绕射波.矩形空洞上界面棱角的双曲线绕射波能量比三角形空洞上界面棱角的双曲线绕射波能量强，矩形空洞上界面棱角的双曲线绕射波弧形两翼比三角形空洞上界面棱角的双曲线绕射波的弧形两翼宽.

图 4

Fig. 4

图 4 模型2正演剖面图 Fig. 4 Simulate profile map of model 2

3.3 模型3

如图 5所示，模型3是一个2.0 m×0.5 m的矩形区域，模型上部是空气层，厚度为0.05 m；中部为衬砌混凝土层，其相对介电常数为8，电导率σ=0.005 s/m，混凝土层中在0.35 m处平均分布10根Φ18型号钢筋，钢筋间距为0.2 m. 模型下部为泥岩，其相对介电常数为20，电导率σ=0.01 s/m，泥岩层上部存在2个空洞，空洞相对介电常数为1，电导率σ=0.01 s/m，矩形空洞大小为0.2 m×0.1 m，三角形空洞底为0.2 m，高为0.1 m，Ricker子波主频率为900 MHz.

图 5

Fig. 5

图 5 模型3示意图 Fig. 5 Sketch map of model 3

如图 6为模型3正演剖面图，图中可见：钢筋异常体异常表现明显，连续的双曲线绕射波分布均匀，能很好地反应钢筋的位置与分布情况.虽然由于异常体比较密集，存在反射波、绕射波、多次波等混合叠加，相互干扰等情况，剖面图相对比较复杂，但是仍然能从剖面图中分辨出混凝土与泥岩层的分界面以及空洞异常体的上界面，空洞异常体的上界面反射波是由两个棱点的绕射波和棱点间直线界面的反射波组成，由于空洞上边界边长较短，所以棱点间的直线界面反射波表现得不够明显.从剖面图中能够分辨出三角形空洞异常体的下界面，矩形异常体的下界面则很难分辨出.模型3中的钢筋间距与实际建设中大部分引水隧洞衬砌中设计的钢筋间距相同，正演结果说明在钢筋间距为0.2 m时，探地雷达能够发现衬砌中及衬砌下存在的空洞缺陷.

图 6

Fig. 6

图 6 模型3正演剖面图 Fig. 6 Simulate profile map of model 3

4 工程应用实例与分析

本次隧洞衬砌检测对象为福建省某水电站引水隧洞，引水隧洞剖面形态主要为马蹄形及圆形切割两种形态，洞径10.4 m，钢筋配筋方式为双筋、单筋及构造配筋，钢筋型号主要为Φ12和Φ18，钢筋密度与本文正演模型3相同，钢筋间距为0.2 m.本次隧洞衬砌检测实施方案拟共布置7条测线，拱顶中线布设一条测线，洞顶左右约120°各一条测线，左右洞壁1.5 m高度各布设一条测线，左右洞壁0.5 m高度各布设一条测线.本次隧洞衬砌检测主要检测目标为脱空及不密实检测以及是否存在钢筋缺失及钢筋间距过大.

本次隧洞衬砌检测综合介质特性及探测目标的尺寸与埋深等因素，检测过程中选用900 MHz和400 MHz天线结合，采用点测、自激自收、天线与测线同步平移，每道采样点为1024，采样间隔为0.1 m，利用皮尺加测距轮控制精度.图 7为本次检测合格的典型剖面图，如图所示，剖面图中未见脱空异常，钢筋异常表现非常明显，钢筋排列均匀，与正演模拟结果相同，在剖面图标记的4 m范围内均匀分布20根钢筋，与设计相符.图 8为本次检测发现存在钢筋间距过大的典型剖面图，在剖面图标记的10 m范围内只分布24根钢筋，间距为0.42 m，与设计的0.2 m间距不相符.经过开挖验证，图 8所在实际位置10 m内实际存在24根钢筋，与本次检测结果相同，证实该处实际存在钢筋间距过大的情况.

图 7

Fig. 7

图 7 实测剖面图 1(900 MHz) Fig. 7 Profile map 1 of actual measurement(900 MHz)

图 8

Fig. 8

图 8 实测剖面图 2(900 MHz) Fig. 8 Profile map 2 of actual measurement(900 MHz)

图 9为本次检测发现存在脱空缺陷的典型剖面图，图中所圈定的椭圆形区域为脱空异常区域，区域内反射强烈、无序、信号杂乱，由于实际存在的空洞缺陷几乎不可能是规则形状，所以实际检测的脱空缺陷异常表现与正演所示的脱空缺陷异常表现有所不同，实际空洞由于形状不规则，存在很多棱点，所以存在很多绕射波，且绕射波不完整.但是，它们的基本特征和正演结果相吻合：脱空缺陷位置存在双曲线绕射波，反射波同相轴不连续，多次反射明显.经过业主后期开挖验证，证实圈定区域存在脱空缺陷，与检测结果相符.

图 9

Fig. 9

图 9 实测剖面图 3(400 MHz) Fig. 9 Profile map 3 of actual measurement(400 MHz)

5 结 论

5.1     给出了二维TM波时域有限差分法探地雷达正演模拟方程，介绍了UPML边界条件，推导得出基于UPML边界条件探地雷达二维正演模拟的FDTD表达式，利用正演模拟明确了钢筋及空洞缺陷的异常表现形式，证明了在钢筋间距为0.2 m时探地雷达能够探测钢筋下存在的空洞缺陷.

5.1    GPR正演结果和实测资料表明：钢筋及空洞缺陷在实测雷达剖面图上的异常表现的特征与正演结果相符.GPR正演模拟能够有效的指导探地雷达实测资料解释.