基于三角形网格频域有限元的Debye频散介质探地雷达数值模拟

Numerical Simulation of Debye Dispersive Ground Penetrating Radar Based on Triangular Grid Frequency Domain Finite Element

摘要: 通过数值模拟研究探地雷达(GPR)高频电磁波在频散介质规律，对指导实测雷达资料的处理与解释具有重要意义。频率域数值模拟通过计算一系列频率的响应特征并利用傅里叶变换转换到时间域，具有计算稳定，精度高的优点，而在GPR数值模拟中得到广泛应用。本文发展了基于三角形网格频率有限元法的Debye频散介质GPR数值模拟算法。从Debye频散介质模型满足的频率域电磁波动方程出发，在复拉伸坐标系下通过合理构造合理辅助微分方程，推导了完全匹配层(PML)边界条件的二阶电磁波方程。在此基础上，利用Galerkin法推导频率域三角形网格的有限元方程。数值模拟的结果表明：本文构建的PML边界条件可较好地吸收边界截断边界处的超强反射波；相较于非频散介质，GPR高频电磁波Debye频散介质中传播时能量衰减更快，子波持续时间更长，传播速度和分辨率降低，分析结果有利于GPR数据解释和反演。

Abstract: Numerical investigation into the propagation characteristics of high-frequency electromagnetic waves generated by GPR in dispersive media is of critical importance for guiding the processing and interpretation of field-measured radar data. Numerical simulation in the frequency domain, which computes electromagnetic response characteristics across discrete frequencies and transforms them into the time domain via Fourier transform, offers enhanced computational stability and high accuracy, thereby enabling its widespread application in GPR modeling. This study develops a GPR numerical simulation algorithm based on the frequency-domain finite element method with triangular meshes for Debye dispersive media. Starting from the frequency-domain electromagnetic wave equations governing the Debye dispersive medium model, the second-order electromagnetic wave equations with PML boundary conditions were systematically derived in the complex-stretched coordinate system through the appropriate construction of auxiliary differential equations. Building upon this foundation, the finite element equations for the frequency-domain triangular meshes were rigorously derived using the Galerkin method. The numerical simulation results demonstrate that the PML boundary conditions constructed in this study effectively absorb strong spurious reflections at truncated boundaries. Compared to non-dispersive media, high-frequency GPR electromagnetic waves propagating in Debye dispersive media exhibit faster energy attenuation, prolonged wavelet duration, and reduced propagation velocity and resolution. These analytical findings provide valuable insights for GPR data interpretation and inversion.

1. 引言

探地雷达(GPR)作为一种高效的无损探测技术，已广泛应用于地质勘探、工程检测和考古研究等领域。然而，实际探测环境中，地下介质往往表现出频散特性，即介电常数随频率变化，这会导致电磁波传播速度的改变和波形的畸变，进而影响GPR探测的精度和分辨率[1]。

在频率域中，电位移矢量可以表示为电场强度和复介电常数的乘积，因此在频率域模拟GPR电磁波具有原理简单和易实现的优点，1900年，Pratt和Worthington [2]对频率域正演进行了深入研究。他们采用常规二阶差分法(5点有限差分格式)对声波方程进行离散化处理，构建了阻抗矩阵。得出结论，每个波长至少需要13个网格点，才能将相速度误差控制在1%以内。这项研究为频率域正演奠定了基础，但由于数值频散严重和计算量较大等问题，频率域正演并未得到广泛应用。Shin [3]提出了旋转9点差分格式，使得每个波长只需4个网格点即可达到1%的精度要求。这一改进极大地促进了频率域正演的发展。Shin [4]进一步将9点差分格式扩展到25点差分格式：将单个波长内的网格点减少到2.5个，进一步提高了模拟精度。在单频波场的求解中，LU分解约占用88%的运行时长，主要计算量集中在对阻抗矩阵的分解算法上面，2007年，吴国忱等率先开始关于频率域正演的模拟研究，并将并将25点最优化加权平均差分算子成功运用到弹性波和各向异性介质中的高精度有限差分模拟研究中[4]；随后，任浩然等对吸收边界条件和大型稀疏矩阵存储技术进行了研究，有效地减少了计算所需的存储量，推动了频率域正演的发展[5]。

正演模拟中的边界条件研究已相对成熟。自Berenger (1994) [6]提出完全匹配层(PML)吸收边界条件以来，PML因其卓越的吸收性能成为电磁波模拟中最广泛使用的边界条件(Liu等，2019)。Chew和Weedon (1994) [7]通过引入复伸展坐标系，将PML边界公式化，为其数值模拟奠定了理论基础。近年来，PML边界条件在声波和弹性波模拟中也得到了广泛应用。为进一步改善PML边界在某些介质中的不稳定性，Kuzuoğlu和Mittra (1996) [8]提出了复频移PML (C-PML)边界条件。在此基础上，学者们研究并实现了非分裂形式的C-PML边界，显著增强了其简捷性和实用性。

本文研究的基于三角形网格频率有限元法的Debye频散介质GPR数值模拟算法相对于时间域的数值模拟来说，有着自身的优势[9]。首先，频率域相对于时间域的计算效率更高，每个频率点的阻抗矩阵只需要计算一次，加入后的计算可以极大地提升计算效率；其次，频率域的模拟是计算区域所有网格点上的全部频率点的解，通过多次傅里叶变换得到模拟结果，误差将分配到所有参与计算的每一个网格点上，不会造成累积误差，而时间域的结果有可能因为某一时间的结果计算错误而产生累积误差，所以在频率域的数值模拟中，可以考虑长时间的GPR电磁波数值模拟。

2. 三角形网格频率域有限元Debye频散介质GPR数值模拟方法

2.1. Debye频散介质中电磁波动方程

根据电磁波传播理论，假设目标体走向为y轴，二维地电条件下的频率域GPR高频电磁波在地下介质中满足的电磁波方程为[10]

$\left(\frac{{\partial }^2{E}_y}{\partial x^2}+\frac{{\partial }^2{E}_y}{\partial z^2}\right)+k^2{E}_y=0$ (1)

式中，$E$ 为电场强度(V/m)，$k=\sqrt{{\omega }^2\mu \epsilon -i\omega \mu \sigma }$ 为电磁波波数，$\epsilon ,\mu ,\sigma$ 为介质的介电常数(F/m)、磁导率(H/m)、和电导率(S/m)。

考虑电阻率的影响，介质相对介电常数$\epsilon$ 随频率变化的Debye频散模型可表示为[11]

$\epsilon \left(\omega \right)={\epsilon }^{\infty }+\frac{\sigma }{i\omega {\epsilon }_0}+\frac{{\epsilon }^0-{\epsilon }^{\infty }}{1+i\omega \tau }$ (2)

为吸收模型截断边界处的超强反射波，我们采用各向异性完全匹配层(Uniaxial Perfectly Matched Layer, UPML)作用于模型边界处。假设复拉伸坐标下的坐标伸缩因子为$S_x=S_x\left(x\right)$ 、$S_z=S_z\left(z\right)$ 。根据UPML理论，式(1)施加UPML边界条件后的方程可表示为：

$\frac{\partial }{{\partial }_x}\left(\frac{S_z}{S_x}\frac{\partial E_y}{{\partial }_x}\right)+\frac{\partial }{\partial z}\left(\frac{S_x}{S_z}\frac{\partial E_y}{\partial z}\right)+S_x{S}_z{k}^2{E}_y=0$ (3)

其中$S_x=1-i{\sigma }_x/{\epsilon }_0\omega$ ，$S_y=1-i{\sigma }_y/{\epsilon }_0\omega$ ；${\epsilon }_0$ 为真空的介电常数；${\sigma }_x$ 和${\sigma }_z$ 的具体表达式如下所示：

${\sigma }_{\xi }={\sigma }_{\xi }{}_{,\max }{\left(\frac{l}{d}\right)}^m,\xi =x,z$ (4)

其中${\sigma }_{\xi }$ 是$x$ 和$z$ 方向上的UPML参数；$d$ 为UPML层厚度；$l$ 表示UPML区域内的计算点至UPML区域内边界的距离；${\sigma }_{\xi ,\max }$ 为最大损耗参数；$m$ 是UPML的指数参数，通常在3和4之间，一般取$m=4$ 。R表示法向入射目标的反射系数，通常取$R={\text{e}}^{-10}$ ，最大损耗参数${\sigma }_{\xi ,\max }$ 表示为

${\sigma }_{\xi }{}_{,\max }=-\frac{\left(m+1\right)\ln R}{2{\eta }_{0}d\sqrt{{\epsilon }_r}}$ (5)

其中，${\eta }_0=\sqrt{{\mu }_0/{\epsilon }_0}$ 为自由空间波阻抗。

令$A=S_z/S_x$ ，$B=S_x/S_z$ ，$C=S_x/S_z$ ，$\mu =E_y$ ，式(3)可写为

$\frac{\partial }{\partial x}\left(A\frac{\partial u}{\partial x}\right)+\frac{\partial }{\partial z}\left(B\frac{\partial u}{\partial z}\right)+C{k}^{2}u=0$ (6)

2.2. 网格剖分与插值基函数

由Galerkin方法，式(6)相应的弱解形式为

${\displaystyle \underset{\Omega }{\int }\phi \left[ - \left(A\frac{\partial \phi }{\partial x}\frac{\partial u}{\partial x}+B\frac{\partial \phi }{\partial y}\frac{\partial u}{\partial y}\right)+C{k}^2\phi u-\phi q\right]\text{d}\Omega }={\displaystyle \underset{\Gamma }{\int }\phi \left(A\frac{\partial u}{\partial x}{n}_x+B\frac{\partial u}{\partial y}{n}_y\right)\text{d}\Gamma }$ (7)

其中，$\phi$ 为权函数，而由于计算边界为完美电导，上式右端为0，因此式的弱解形式可简化为

${\displaystyle \underset{\Omega }{\int }\phi \left[ - \left(A\frac{\partial \phi }{\partial x}\frac{\partial u}{\partial x}+B\frac{\partial \phi }{\partial y}\frac{\partial u}{\partial y}\right)+C{k}^2\phi u-\phi q\right]\text{d}\Omega }=0$ (8)

本文利用如图1所示的三角形单元进行模型的网格剖分，其三个顶点按逆时针方向编号，分别记为1，2，3对应的坐标分别为$\left(x_1,z_1\right),\left(x_2,z_2\right),\left(x_3,z_3\right)$ ，其场值分为$u_1,u_2,u_3$ ，其面积为$\Delta$ 。

Figure 1. Schematic diagram of the triangular grid elements

图1. 三角形网格单元示意图

三角形内某个点的场函数可表示为

$u=N_1{u}_1+N_2{u}_2+N_3{u}_3$ (9)

其中，$N_1,N_2,N_3$ 为形函数矩阵中的元素：

$\{\begin{array}{l}{N}_1=\frac{1}{2\Delta }\left(a_{1}x+b_{1}z+c_1\right)\\ N_2=\frac{1}{2\Delta }\left(a_{2}x+b_{2}z+c_2\right)\\ N_3=\frac{1}{2\Delta }\left(a_{3}x+b_{3}z+c_3\right)\end{array}$ (10)

单元质量矩阵$M$ 和单元刚度矩阵$K$ 的计算式如下：

$\{\begin{array}{l}M={\displaystyle \underset{e}{\iint }N{N}^{\text{T}}\text{d}x\text{d}z}\\ K=\frac{1}{\mu \epsilon }{\displaystyle \underset{e}{\iint }\left[\frac{\partial N}{\partial x}{\left(\frac{\partial N}{\partial x}\right)}^{\text{T}}+\frac{\partial N}{\partial z}{\left(\frac{\partial N}{\partial z}\right)}^{\text{T}}\right]\text{d}x\text{d}z}\end{array}$ (11)

计算区域划分单元后，总体积分可以写成各个单元积分之和，可表示为：

${\displaystyle \sum _{e=1}^{N_e}{\displaystyle \underset{{\Omega }^e}{\iint }\left[ - \left(A\frac{\partial \phi }{\partial x}\frac{\partial u}{\partial x}+B\frac{\partial \phi }{\partial y}\frac{\partial u}{\partial y}\right)+C{k}^2\phi u-\phi q\right]\text{d}\Omega }}=0$ (12)

再将单元场向量、激励源及单元系数矩阵扩展成整体矩阵，有如下矩阵表达形式：

$-Ku+Mu=f$ (13)

从而形成多源系统的线性方程组：

$Au=f$ (14)

式中，$A=M-K$ 是与频率相关的所有源公用的系数矩阵；$f$ 是离散化激励源项；$u$ 是对应节点的离散电场值。使用LU分解法求解线性方程组(14)，LU直接只对系数矩阵求逆一次就能同时获得所有源激励的解。

由于A与频率相关，一个新的频率需要重新对A进行计算。$u$ 是二维网格节点离散电场值，为了获得观测数据$d_{cal}$ ，需要知道接收点对应的位置矩阵P，因此测量数据与正演场值关系可表示为

$d_{cal}=pu$ (15)

3. 数值试验

3.1. 均匀模型

为了测试本文构建的Debye频散介质GPR频率域有限元模拟的PML边界条件下的正确性和有效性，建立一个均匀半空间介质下，介质的类型是否为频散介质。

如图2为一个大小为1 m × 4 m的狭长型均匀半空间模型，模型内灰色部分为0.15 m厚PML边界层，计算区域被三角形网格剖分成50,674个单元，三角形网格的最大边长为0.1 m，为保证时间迭代的稳定性，时间步长取0.01n s，非频散介质的相对介电常数为4，电导率为0.0001，频散介质的复相对介电常数为8+4i，电导率为0.0001，发射源为500 Mhz的雷克子波，位置在左上角(0.175, 0.175)，如图中五角星所示，三个接收点坐标为位置分别为(2.175, 0.175)、(2.175, 0.525)、(0.275, 0.825)，如图中黑色圆点所示。

使用基于PML边界条件下的Debye频散介质与非频散介质GPR频率域有限元法进行模拟，计算后的10 ns、15 ns、25 ns和35 ns波场快照如图3所示。

由图4可知，由于发射源靠近上侧和左侧边界的PML边界层，常规PML边界对于低频反射波和掠射波的吸收不完全，容易产生低频反射波和掠射波，在左侧边界和上侧边界周围可以看到较为明显的低频反射波和掠射波，如图4(a)、图4(b)、图4(e)和图4(f)的白色箭头所示；随着传播距离和入射角的增加，低频反射波的能量也随之增大，由图4(a)~(f)可知，相同条件下，GPR高频电磁波在非频散介质中的传播速度比在频散介质中的传播速度快，并且如图4(a)、图4(b)、图4(e)和图4(h)中显示，PML边界条件下，低频反射波在频散介质中远距离偏移处产生的低频反射波能量更弱。

Figure 2. Narrow and uniform half-space model

图2. 狭长均匀半空间模型

Figure 3. Schematic diagram of triangulation of the long and narrow model

图3. 狭长模型三角剖分示意图

Figure 4. Snapshot diagram of the wave field at different times of the dispersive medium and the non-dispersive medium under the PML boundary conditions

图4. PML边界条件下频散介质和非频散介质不同时刻的波场快照图

图5是接收点1所收到的PML和CFS-PML边界条件下的单道波形与参考解的对比图，图5展示的是Debye频散介质的结果，而图5(b)呈现的是非频散介质的结果。参考解的获取方式为：把模型的规模扩大至原来的3倍，同时将发射源移至模型的中心位置，在此过程中保持接收点与发射源的相对位置恒定不变。由于此时模型的尺寸足够大，因而在观测范围所对应的时间窗内，不会接收到源自边界的反射波。从图5中可以看出，在接收点1处，基于PML边界条件所得到的单道波形，无法与参考解实现良好的拟合，存在着较大的误差。CFS-PML边界条件的单道波形与参考解几乎完全拟合，与参考解误差几乎为0；结果表明CFS-PML边界条件对掠射波与低频反射波效果强于PML边界条件，与图4中得到的结果相互印证。

(a) Debye频散介质 (b) 非频散介质

Figure 5. Comparison of single-channel waveform results received by receiving point 1

图5. 接收点1接收的单道波形结果对比图

在接收点2和接收点3处，PM边界条件下L和CFS-PML边界条件下的单道波形与参考解的对比如图6与图7所示，接收点2与接收点3和发射源的相对位置导致其接受的信号为接近垂直入射波或小角度入射波，在此条件下，Debye频散介质和非频散介质在接收点2与接收点3处所接收的单道波形与参考解拟合程度都比较高，误差几乎不见。由此可见，PML边界条件与CFS-PML边界条件对垂直入射波以及小角度入射波的吸收效果差异不大，CFS-PML边界条件并无明显优势。

(a) Debye频散介质 (b) 非频散介质

Figure 6. Comparison of single-channel waveform results received by receiving point 2

图6. 接收点2接收的单道波形结果对比图

3.2. 均匀介质模型

为研究GPR高频电磁波在Debye频散介质中的传播规律，构建了一个3.0 m × 3.0 m的双层地电模型，发射源采用500 MHz的雷克子波，位于模型左上角(0 m, 0 m)处，时间步长$dt=0.01ns$ ，采样点数$N\text{r}=4000$ ，该模型地层1非频散介质的相对介电常数$\epsilon =3$ ，电导率$\sigma$ = 0.001 s/m；频散介质的复相对

(a) Debye频散介质 (b) 非频散介质

Figure 7. Comparison of single-channel waveform results received by receiving point 3

图7. 接收点3接收的单道波形结果对比图

介电常数${\epsilon }_r^0=5$ ，${\epsilon }_r^{\infty }=3$ ，电导率$\sigma$ = 0.001 s/m，弛豫时间$\tau$ = 0.10 ns。地层2非频散介质的相对介电常数$\epsilon =8$ ，电导率$\sigma$ = 0.03 s/m；频散介质的复相对介电常数${\epsilon }_r^0=10$ ，${\epsilon }_r^{\infty }=8$ ，电导率$\sigma$ = 0.03 s/m分别使用非频散介质与Debye频散介质对模型进行GPR正演模拟(图8)。

Figure 8. Schematic diagram of the double-layer geoelectric model

图8. 双层地电模型示意图

将图9所示双层地电模型示意图使用Comsol软件进行非结构化三角剖分，剖分结果如图9所示。整个模型被80,532个节点剖分成160,046个三角形单元。(图10)

Figure 9. Unstructured triangular meshing results of a double-layer geoelectric model

图9. 双层地电模型非结构化三角形网格剖分结果图

(a) 频散介质13 ns；(b) 频散介质35 ns；(c) 非频散介质13 ns；(d) 非频散介质35 ns。

Figure 10. Snapshot results of the Debye dispersive medium and non-dispersive medium of the double-layer geoelectric model

图10. 双层地电模型Debye频散介质与非频散介质波场快照结果图

为探究GPR电磁波在不同介质中的传播特性，在特定实验设置下进行了相关测试。把接收天线分别置于(0.5 m, 0 m)、(1.0 m, 0 m)、(1.5 m, 0 m)、(2.0 m, 0 m)这几个位置，由此获取了偏移距依次为0.5 m、1.0 m、1.5 m、2.0 m的GPR电磁波时域单道波形对比图，具体呈现于图11。从该图能够清晰地观察到，在偏移距一致的情形下，相较于在非频散介质里传播，电磁波于Debye频散介质中传播时所携带的能量更为微弱，这一现象充分表明GPR电磁波在Debye频散介质传播过程中会出现强烈的衰减情况。进一步分析发现，随着偏移距逐渐增大，GPR电磁波在Debye频散介质中传播时的能量衰减愈发显著。这主要是因为偏移距越大，GPR电磁波在Debye频散介质内所需传播的时长就越长，进而致使Debye频散介质的频散效应愈发凸显，此结论可通过对图11(a)~(d)的对比得出。例如，在收发天线偏移距达到2.0 m时，接收天线接收到的Debye频散介质振幅能量大致仅为非频散介质振幅能量的五分之一，这一情况在图11(d)中清晰可见。

(a) 0.5 m; (b) 1.0 m; (c) 1.5m; (d) 2.0 m。

Figure 11. Single-channel waveform of a receiving point with different offsets

图11. 不同偏移距接收点单道波形图

3.3. 复杂模型

为了在复杂地质结构下验证结论的准确性和有效性，故建立如图12所示的起伏界面模型。图12中，上界面为非频散介质，相对介电常数为8，电导率为0.001；下界面为频散介质A，相对介电常数为14+4i，电导率为0.01；其余部分为异常体，使用频散介质B，相对介电常数为32+4i，电导率为0.03。

Figure 12. Schematic diagram of the undulating interface model

图12. 起伏界面模型示意图

由图13(a)和图13(c)可知，10 ns时，电磁波传播至起伏界面分界处，此时，频散介质刚刚接触起伏界面，非频散介质与频散介质产生明显的波速差异，非频散介质的波速要略大于频散介质的波速；由图13(b)和图13(d)的黑色箭头所位置可见，非频散介质中，电磁波在传播至圆形空洞后产生绕射波，但在频散介质中，绕射波不明显。图14为PML边界条件下频散介质和非频散介质的多偏移距剖面，可以明显看出起伏界面中，频散介与非频散介质的波速差异约为3 ns，在图14(b)黑色箭头处也可看到非频散介质所产生的绕射。

Figure 13. Snapshot of the wave field for 10 ns and 20 ns for dispersive and non-dispersive media under PML boundary conditions. (a), (b) are non-dispersive media, (c) and (d) are dispersive media

图13. PML边界条件下频散介质和非频散介质10 ns和20 ns的波场快照图。(a)、(b)为非频散介质，(c)、(d)为频散介质

4. 结论

为了实现GPR频率域数值模拟，了解GPR高频电磁波在不同介质中的传播规律，本文基于复拉伸坐标系下的Debye频散介质电磁波动方程，结合常规PML偏微分边界条件，推导了Debye介质的二阶矢量微分方程。利用Galerkin方法，得到了时谐场的PML吸收边界的Maxwell方程，并采用频域有限元法对GPR电磁波在频散介质中的传播进行了数值模拟。构造了三个不同的模型进行对比实验，分别进行

Figure 14. Multi-Offset profiles of dispersive and non-dispersive media under PML boundary conditions. (a) is a non-dispersive medium and (b) is a dispersive medium

图14. PML边界条件下频散介质和非频散介质的多偏移距剖面。(a)为非频散介质，(b)为频散介质

了算法的验证、CFS-PML边界条件的吸收效果分析、磁波在Debye频散介质中的传播规律研究、四边形网格剖分与非结构化三角形网格剖分的数值模拟精度对比。针对本文的研究，所得到的结论如下：

在狭长模型的数值模拟中，对比两种不同边界条件的吸收效果，通过分析可得：相比于PML边界条件，CFS-PML边界条件对于大角度入射的电磁波产生的低频反射波吸收效果更好。

在双层地电模型中，对于电磁波在非频散介质与Debye频散介质中传播的规律进行对比分析，结果表明，相较于非频散介质，电磁波在Debye频散介质中的传播能量较低，传播速度更慢。

构建复杂起伏界面模型，通过不同时刻的波场快照以及多偏移GPR剖面图，研究复杂介质中电磁波在非频散介质与Debye频散介质中传播的规律。非结构化三角形网格能够对目标区域进行自适应加密，使计算聚焦于目标体或介电常数变化界面附近，数值模拟精度更高，效果更佳。