高效和无损的探地雷达(ground penetrating radar,GPR)技术已广泛应用于公路及铁路工程结构体厚度及缺陷的试验检测中^[1-2]，然而介质中雷达波速标定带来的误差常常严重影响试验检测的精度，使得高效无损的探地雷达技术未充分发挥其高分辨率优势，在公路及铁路工程的质量评定过程中仅作为辅助的参考手段，为提高试验检测效率和精度，规范波速校正的试验检测过程十分必要。相对介电常数和电导率是影响雷达波速的主要介电参数，探索它们之间的相关关系是解决该问题的关键，数值试验是细致研究介电参数特征和雷达波传播规律的有效途径，目前常用的数值实验工具主要为有限元法(finite element method,FEM)和有限差分法(finite difference method,FDM)，基于麦克斯韦方程组的时域有限差分法以其节约计算时间和存储空间、直接时域范围计算、计算程序通用性强等优点^[1-2]，在模拟雷达波传播中得到了广泛的应用，Teixeira等^[3]在完全匹配层吸收边界条件下利用FD法对各向异性的良导介质进行了三维GPR数值模拟，Levent Gürel等^[4]在复杂介质中进行了三维GPR数值模拟，何兵寿等^[5]对矿井地电模型进行了GPR数值模拟，李静等^[6]研究了探地雷达的高阶FD数值模拟，底青云等^[7]推导了含衰减项的探地雷达有限元方程，实现了复杂介质的探地雷达有限元数值模拟，冯德山等^[8-9]研究了基于时域多分辨法的三维GPR数值模拟，刘新荣等^[10]研究了隧道衬砌空洞三维探地雷达正演模拟，为雷达图谱的解释提供了依据，王绪本等^[11]对岩溶洞穴进行了GPR物理模拟，宋华等^[12]采用探地雷达对海堤模型进行了物理模拟试验，Giannopoulos^[13]编写了基于Yee氏网格的Gprmax探地雷达数值模拟工具，张彬等^[14]利采用旋转交错网格FD方法模拟了雷达波在衰减夹层介质中的传播情况，提高了模拟效率的同时改善了对低频隐失波的有效吸收，田钢等^[15]则将FD算法应用于存在反射干扰特征的探地雷达数值模拟，为电磁反射抗干扰研究提供了理论依据。

文中采用二阶精度的FD方法，对不同电导率和不同相对介电常数的钢筋混凝土双介电参数进行了数值实验，着重研究了钢筋混凝土介质的相对介电常数、电导率和波速三者之间的相关关系，总结了双介电参数对雷达波速校正的经验性结果，为减少波速校正带来的误差和规范公路及铁路工程的探地雷达试验检测过程提供参考。

1、 FD法及TM_z差分格式

时域有限差分法于1966年由Yee首次提出^[16]，通过将电场分量E、磁场分量H在时间和空间上进行交替离散化，即每一个电磁E(或H)分量周围由4个电磁场H(E)分量环绕，如图 1所示，并在每个离散点上使用差商来代替微商，将求解包含时间变量的2个麦克斯韦旋度方程和2个散度方程转化为求解6个有限差分方程组，并通过时间轴的推进逐步求取空间上的电磁场值，而实际上，麦克斯韦方程组中2个散度方程可以由2个旋度方程导出，故在数值试验的计算中以旋度方程为基础。

在无源场区域，麦克斯韦方程组中的2个旋度方程表示为如下的形式^[17]：

$\nabla \times E\text{=}-\mu \frac{\partial H}{\partial t}-{{\sigma }_{m}}H,$

(1)

$\nabla \times H\text{=}\varepsilon \frac{\partial E}{\partial t}+{{\sigma }_{e}}E,$

(2)

其中:ε为相对介电常数;σ_m为等效磁阻率,Ω/m，在TM_z极化模式下，只含E_z、H_x和H_y分量的独立方程组为

$\frac{\partial {{E}_{z}}}{\partial t}=\frac{1}{\varepsilon }\left( \frac{\partial {{H}_{y}}}{\partial x}-\frac{\partial {{H}_{x}}}{\partial y}-{{\sigma }_{e}}{{E}_{z}} \right),$

(3)

$\frac{\partial {{H}_{x}}}{\partial t}=-\frac{1}{\mu }\left( \frac{\partial {{E}_{z}}}{\partial y}+{{\sigma }_{m}}{{H}_{x}} \right),$

(4)

$\frac{\partial {{H}_{y}}}{\partial t}=\frac{1}{\mu }\left( \frac{\partial {{E}_{z}}}{\partial x}-{{\sigma }_{m}}{{H}_{y}} \right)\circ $

(5)

采用二阶空间精度的中心差商代替微商，得到二维时域差分方程，即为TM_Z模式下的FDTD更新方程组：

$\begin{align} & \tilde{E}_{z}^{n+1}\left( i,j \right)=CA\left( i,j \right)\cdot \tilde{E}_{z}^{n}\left( i,j \right)+CD\cdot CB\left( i,j \right)\cdot \left[ H_{y}^{n+\frac{1}{2}}\left( i+\frac{1}{2},j \right)- \right. \\ & \left. H_{y}^{n+\frac{1}{2}}\left( i-\frac{1}{2},j \right)+H_{x}^{n+\frac{1}{2}}\left( i,j-\frac{1}{2} \right)-H_{x}^{n+\frac{1}{2}}\left( i,j+\frac{1}{2} \right) \right], \\ \end{align}$

(6)

$H_{x}^{n+\frac{1}{2}}\left( i,j+\frac{1}{2} \right)=H_{x}^{n+\frac{1}{2}}\left( i,j+\frac{1}{2} \right)+CD\cdot \left[ \tilde{E}_{z}^{n}\left( i,j \right)-\tilde{E}_{z}^{n}\left( i,j+1 \right) \right],$

(7)

$H_{x}^{n+\frac{1}{2}}\left( i,j+\frac{1}{2} \right)=H_{x}^{n+\frac{1}{2}}\left( i,j+\frac{1}{2} \right)+CD\cdot \left[ \tilde{E}_{z}^{n}\left( i+1,j \right)-\tilde{E}_{z}^{n}\left( i,j \right) \right],$

(8)

其中系数项$CA\left( i,j \right)=\frac{1-\frac{\sigma \left( i,j \right)\text{d}t}{2\varepsilon \left( i,j \right)}}{1+\frac{\sigma \left( i,j \right)\text{d}t}{2\varepsilon \left( i,j \right)}},CB\left( i,j \right)=\frac{{{\varepsilon }_{0}}}{\varepsilon \left( i,j \right)+\frac{\sigma \left( i,j \right)\text{d}t}{2}},$$CD=\frac{\text{d}t}{\text{d}t}\cdot \frac{1}{\sqrt{{{\varepsilon }_{0}}{{\mu }_{0}}}},\left( i,j \right)$表示计算区域中x轴坐标第i，y轴坐标第j个网格。

2、 介电参数数值实验的算法

从更新方程可以看出，电场分量Ez在(n+1)Δt时刻获取更新，而磁场分量H_x和H_y则均在(n+0.5)Δt时刻获取更新，在各个边界处理上，采用了完全匹配层 (PML) 吸收边界条件^[18-19]，详细的介电参数数值实验流程如图 2所示。

3、 典型混凝土模型的数值实验 (3.1) 钢筋混凝土路面模型的单介电参数数值实验

不同钢筋埋深的混凝土介质模型如图 3所示，模型中设置钢筋的横纵向间距分别为10 cm和5 cm，最上部的钢筋深度为10 cm，最底部的钢筋深度为35 cm，模拟时采用天线中心频率为800 MHz，设置混凝土的相对介电常数ε₁为9.0，电导率σ₁为0.000 1 S/m，钢筋为导电属性的理想电导体(PEC)介质，直径均为0.02 m，网格的横向空间迭代步长为0.003 m，纵向空间步长为0.001 5 m，时间步长均为0.005 ns，各钢筋占据8个网格，吸收边界区域设置为5个网格，共完成了40道数据的模拟采集。

数值实验的模拟剖面结果如图 4所示，图中可见绕射波同相轴呈双曲形状，双曲线的曲率随钢筋异常体埋深的增大而减小，同样，绕射波的能量随钢筋异常体埋深的增加而减小，最底部钢筋的双曲线形态难以分辨，其中因多次反射波和绕射波的产生，深部的绕射信号难以分辨，该模拟剖面真实地反映了钢筋混凝土模型中不同钢筋埋深的基本特征。

图 5为该模拟数据的第20道单道波的电场分量Ez，图 6为Ez电场分量在不同传播时刻的波场快照，4个时刻的能量显示已作均一化处理，从模拟剖面、单道波信号和波场快照结果可以看出，钢筋介质的电导率严重影响了雷达波能量继续往深部传播，成为了制约雷达波波速的主要因素。

从不同时刻的Ez分量波场快照可以看出，雷达波在相对介电常数为9.0的背景介质中传播时，因混凝土介质中钢筋埋深各异，雷达拟球面波在钢筋位置产生了绕射异常，并同时以拟球面波的形式向外扩散。从第20道信号的单道波振幅中看出，雷达波的能量衰减严重，且随深度增加，雷达波在钢筋混凝土介质内部的传播速度呈现递减趋势，根据均匀介质中雷达波速理论，在5.0 ns时刻，雷达拟球面波的波前面应传播至模型的下边界，从图 6所示的数值实验结果来看，在6 ns时刻，雷达波的波前面仅传播至第210个网格位置，即传播至210×0.001 5 m=0.315 m位置，波前面的扩散速度受到了电导介质(钢筋)的严重影响。

(3.2) 钢筋混凝土模型的双介电参数数值实验

采用同样的方法进行了2组数值实验，针对相对介电常数为9.0的介质，对电导率分别为0.1 S/m、0.01 S/m的钢筋混凝土路面模型进行数值实验，其他参数保持一致，图 7为不同电导率混凝土模型的第20道单道波信号。同时，针对电导率为0.000 1 S/m的介质模型，对相对介电常数分别为30、15、5的钢筋混凝土模型进行相应的数值实验，图 8为不同介电常数混凝土模型的第20道单道波信号。

如图 7中红色曲线所示，当钢筋混凝土介质的电导率为0.1 S/m时，Ez分量的单道波信号迅速衰减，且仅能观测到直达波信号，在第150个样点到第512个样点的深度区间，Ez分量的幅值保持为零，即在该情况下雷达波对混凝土模型中的钢筋异常体无绕射波产生，波形形态发生了严重畸变，电导介质对雷达波的传播速度影响较严重；相反，对于电导率为0.001 S/m的算例，雷达波在钢筋异常体位置的绕射波信号突出，易于识别，对后期的雷达波的衰减影响不大，雷达波的传播速度变化较小；由蓝色曲线可知当电导率为0.01 S/m时，传播后期的Ez分量具有较小的幅值，即电导介质对雷达波的传播速度产生了一定程度的影响。

如图 8所示，在混凝土介质模型的相对介电常数为5时，各个时刻雷达波的响应时间最早，见图中的黑色曲线，表明该情况下雷达波波速趋于常数$\frac{c}{\sqrt{\varepsilon }}$；红色曲线为相对介电常数为30的算例，由图可知，红色曲线的响应出现较严重的滞后现象；而相对介电常数为15的算例则处于二者之间，蓝色曲线存在有一定程度的滞后现象，表明该相对介电常数对雷达波的波速产生了一定程度的影响。

(3.3) 电导率、介电参数和波速的相关关系

根据波场快照的能量图谱求取了对应的综合雷达波波速，并将该结果与统计的介电参数数值实验数据进行关联，三者的离散相关关系如表 1所示。

将上述相关关系的离散点进行统计，分别采用了三次插值和双调和样条插值的方法，将离散的三维曲面进行拟合，2种方法拟合的适合度信息如表 2所示。

从拟合的适合度信息表中看出，2种方法的均方根误差都为1，表明2种方法拟合的相关关系的变量对雷达波速都具有很强的解释能力，即拟合的适合度较高，而双调和样条插值方法的总平方误差小于三次插值的5个数量级，前者拟合的适合度更高。从图 9、图 10可以看出，其相关关系表现为一个分段的曲面，由三次插值和双调和样条插值的分段曲面可知，综合雷达波波速是相对介电常数和电导率共同影响的结果，如图 10所示，在电导率小于0.04 S/m和相对介电常数小于15的曲面内，雷达波速出现一个上升的缓冲区，数据的离散程度较小，即在该曲面范围内，钢筋混凝土介质的雷达波波速的校正可由$\frac{c}{\sqrt{\varepsilon }}$来获取，由相对介电常数所引起的误差较小，而当双介电参数在该曲面范围之外时，由$\frac{c}{\sqrt{\varepsilon }}$来校正雷达波波速所引起的误差将随双介电参数增大而增加。

2种插值方法拟合的介电参数与综合雷达波波速的相关关系曲面分别如图 9、图 10所示。

4、 结 论

1)将双介电参数数值实验的结果应用于探地雷达试验检测中的波速校正，采用FD方法，对不同电导率和不同相对介电常数的钢筋混凝土模型进行了双介电参数数值实验，并从数值模拟剖面、单道波信号、电场分量波场快照3个方面，讨论了钢筋混凝土介质的相对介电常数、电导率和雷达波波速三者之间的关系；

2) 采用三次插值和双调和样条插值方法，拟合了双介电参数与雷达波波速之间的相关关系，得到了双介电参数对雷达波波速的相关关系校正的结论，为相关的探地雷达试验检测提供参考；

3) 双介电参数与雷达波波速的相关关系呈分段曲面分布，在电导率小于0.04 S/m和相对介电常数小于15的曲面内，雷达波波速的数据离散程度较小，采用$\frac{c}{\sqrt{\varepsilon }}$来校正介质中的雷达波波速所引起误差较小；当双介电参数在该曲面范围之外时，由$\frac{c}{\sqrt{\varepsilon }}$来校正的雷达波波速误差较大，不能单纯地依据钻芯法来校正波速，应以相对介电常数和电导率的芯样综合测试结果来进行波速校正。