GPU探地雷达复杂介质DGTD正演方法

摘要

本发明公开了一种GPU探地雷达复杂介质DGTD正演方法，包括建立DGTD正演模型并初始化；将正演模型三角剖分并设置基函数的阶次和节点；设置激励源和接收点的位置；计算正演所需要的参数并开启GPU并行计算；开始时间积分；计算通量并更新电磁场的场值；在UPML区域更新电磁场辅助场；增加时间步并重复上述步骤直至完成整个时间步的模拟；数据输出；重复上述步骤直至所有激发完成。本发明方法具有高精度、能与非结构化网格结合、可用来模拟复杂介质正演的优点，提高了雷达正演效率，能够更好的模拟复杂地质，而且兼顾了计算效率和计算精度。

权利要求

1.一种GPU探地雷达复杂介质DGTD正演方法，包括如下步骤：S1.建立DGTD正演模型，并对模型进行初始化；

S2.将步骤S1建立的正演模型进行三角剖分，并设置基函数的阶次和节点；具体为采用如下步骤设置基函数的阶次和节点：A.基函数设置为正交的Legendre多项式；

B.二维情况下，采用如下算式作为标准三角形定义的N阶基函数：式中ξ为标准三角形I内节点的横坐标；η为标准三角形I内节点的纵坐标；(i,j)≥0，且i+j≤N； 为n阶Jacobi多项式，且当α1＝β1＝0时 为Legendre多项式；N为基函数阶次；a1和b1为Legendre多项式的变量，且 b1＝η；

C.采用如下算式作为扭曲函数ω(r)，并利用等间距节点和LGL节点来产生新的插值节点：

LGL

式中NP为节点的个数，且在二维情况下 ri 为第i个LGL节点；

m m

ri为第m个三角单元的第i个等距节点， 为基于ri的Lagrange多项式；r为新插值节点的位置；

S3.设置激励源和接收点的位置；

S4.计算正演所需要的参数，并开启GPU并行计算；

S5.开始时间积分；

S6.计算通量，并更新电磁场的场值；

S7.在UPML区域更新电磁场辅助场；

S8.增加时间步，并重复步骤S5～S7，直至完成整个时间步的模拟；

S9.数据输出；

S10.重复步骤S3～S9，直至所有激发完成，从而完成DGTD正演，得到最终的雷达剖面。

2.根据权利要求1所述的GPU探地雷达复杂介质DGTD正演方法，其特征在于步骤S3所述的设置激励源和接收点的位置，具体为采用如下步骤设置激励源和接收点的位置：在二维TM模式下，采用如下算式表示电流源项积分：式中Ωm为加载电流源的所属单元；Jz为加载的线电流源，且在z方向上加载线电流源时Jz＝f(t)δ(x‑x0,y‑y0)；f(t)为延时的雷克子波的波函数；上标T表示向量的转置， 为基函数在此单元节点上的值所组成的向量；(x0,y0)为源点的坐标；δ(x,y)为狄拉克函数。

3.根据权利要求2所述的GPU探地雷达复杂介质DGTD正演方法，其特征在于步骤S4所述的计算正演所需要的参数，并开启GPU并行计算，具体为在进行激励源激发，并求取场值之前，先根据已知的数据，计算得到DGTD正演所需要的参数；同时，将所有预先计算并得到的DGTD正演参数存入到GPU中，并进行后续场值的计算，从而对正演模型的模拟进行加速，节约计算时间。

4.根据权利要求3所述的GPU探地雷达复杂介质DGTD正演方法，其特征在于步骤S5所述的开始时间积分，具体为采用五级四阶低储存Runge‑Kutta 方法进行时间积分。

5.根据权利要求4所述的GPU探地雷达复杂介质DGTD正演方法，其特征在于步骤S6所述的计算通量，并更新电磁场的场值，具体为采用如下步骤计算通量并进行更新：a.通量采用迎风通量，并采用如下算式计算迎风通量：式中 为数值通量；上标‑表示单元内部信息；上标+表示单元外部信息；

为节点处的场值；为边界上的法线向量；α为保证整体稳定性的系数且取值范围为[0,1]，同时当α＝0时数值通量 转换为迎风通量；

b.采用如下算式进行DGTD更新式中ε为介电常数；σ为电导率；M为单元的质量矩阵且 Ez为z轴方向的电场值；Sx为x轴方向的单位刚性矩阵且 Sy为y轴方向的单位刚性矩阵且 Hx为x轴方向的磁场值；Hy为y轴方向的磁场值；F为边界单元矩阵且 κh、κe、ve和vh均为通量系数，且为Ωm边界处的外法向向量在x方向上的投影； 为Ωm边界处的外法向向量在y方向上的投影； 为基函数；上标+表示单元外部信息；εr为介质的相对介电常数；μr为相对磁导率；μ为磁导率。

6.根据权利要求5所述的GPU探地雷达复杂介质DGTD正演方法，其特征在于步骤S7所述的在UPML区域更新电磁场辅助场，具体为采用如下步骤进行更新：(1)采用如下算式表示时间域上引入了辅助场P、Q的UPML方程：式中上标‑1表示逆运算；κ、a、b、c和d均为对角张量，且κ＝diag(κx,κy,κz)，其中κi为改善PML对表面

波的吸收的参数，i取值为x，y或z；σi表示PML层内i方向的电导率，i取值为x，y或z；ε0为真空中的介电常数；

(2)在二维TM模式下，只有Hx、Hy、Ez三个分量，辅助电场P中Px＝Py＝0；辅助磁场Q中Qz＝

0；

(3)采用DGTD算法来求解步骤(1)中的算式，得到二维TM模式下UPML区域的半离散强格式更新方程如下所示：

式中 Pz为辅助

电场P的z轴分量；

Qx为辅助磁场Q的x轴分量；

Qy为辅助磁场Q的y轴分量。

说明书

GPU探地雷达复杂介质DGTD正演方法

技术领域

[0001] 本发明属于地球物理探测领域，具体涉及一种GPU探地雷达复杂介质DGTD正演方法。

背景技术

[0002] 探地雷达(ground penetrating radar，GPR)是一种广泛使用的地球物理勘探技术，它利用高频电磁(electromagnetic，EM)来成像浅层地下介质的特性。凭借高分辨率、高
效率、直观结果和无损检测等各种优势，GPR已用于各种市政、地质、民用、地球物理和环境
应用中的近地表勘探。

[0003] 但是，由于介质分布不均、自然干扰和其他人为干扰，各种反射、衍射以及其他干扰杂波并存，并且在雷达剖面中彼此混淆，因此很难单独区分和提取。为了提高勘探精度并
促进GPR数据的解释，涉及不规则地质体或复杂介质的精细GPR正演模拟对于理解EM波的传
播，尤其是掌握EM回波的详细特征至关重要。

[0004] FDTD算法作为GPR模拟中应用最常用的一种方法，它以易于编程、低成本和节省时间等优点而引起了极大的关注，被广泛用于分析非均匀介质中的电磁散射特性。无条件稳
定的交变方向隐式时域有限差分(alternating direction implicit finite‑difference 
time‑domain，ADI‑FDTD)算法也被开发了出来，专门用于GPR仿真，有效地改善了传统FDTD
算法的时间色散限制。为了降低时间限制，Zhang等将旋转交错网格(rotated staggered 
grid，RSG)FDTD方法和未拆分卷积完美匹配层(unsplit convolutional perfect matched 
layer，UCPML)相结合，实现了GPR仿真。但是，FDTD算法无法与非结构化网格很好地结合。电
磁波在具有不规则几何边界的地质体之间的传播本质上是一个非常复杂的过程，区分出正
确与错误的反射波仍然是一个难题。FETD方法将求解空间划分为一系列单元，利用插值函
数来代替解函数。FETD算法能与非结构网格结合，可以很好地拟合不规则地质体的边界。

[0005] DGTD方法继承了FETD可以与非结构化网格结合、易于实现h‑p型自适应的优点。不同于FETD算法，DGTD勿需求解大型刚度矩阵、易于实行并行计算。DG最早是由Reed和Hill提
出的。Hu等将DGTD算法应用到波的传播中，并开展了由空间离散化引起的误差分析。作为继
任者，许多其他学者开始将DGTD算法应用于电磁学。Bernacki等已经实现了基于非结构化
网格剖分的并行计算，大大减少了计算时间。将DGTD算法应用于GPR仿真时，截断边界的处
理是无限空间仿真的必要前提和关键问题。许多学者研究了DGTD算法的吸收边界条件
(absorption boundary conditions，ABC)。Kabakian将最简单的ABC吸收边界用于DGTD方
法；Xiao和Liu将PML吸收边界引入到DGTD算法中。此外，广泛使用的单轴PML(uniaxial 
PML，UPML)也得到了改进，并应用于DGTD。在GPR模拟方面，Lu等推导了Debye色散材料在
UPML区域的Maxwell方程；Angulo等在简单的3D模型GPR仿真中使用DGTD算法，证明了该方
法的高度收敛性；Yang等使用具有正交四面体多项式基的高阶DGTD方法进行GPR仿真。

[0006] 但是，现有的DGTD方法的正演效率不高，而且无法同时兼顾计算的效率和计算的精度，在一定程度上制约了DGTD方法的应用。

发明内容

[0007] 本发明的目的在于提供一种能够提高雷达正演效率、模拟复杂地质且兼顾计算效率和计算精度的GPU探地雷达复杂介质DGTD正演方法。

[0008] 本发明提供的这种GPU探地雷达复杂介质DGTD正演方法，包括如下步骤：

[0009] S1.建立DGTD正演模型，并对模型进行初始化；

[0010] S2.将步骤S1建立的正演模型进行三角剖分，并设置基函数的阶次和节点；

[0011] S3.设置激励源和接收点的位置；

[0012] S4.计算正演所需要的参数，并开启GPU并行计算；

[0013] S5.开始时间积分；

[0014] S6.计算通量，并更新电磁场的场值；

[0015] S7.在UPML区域更新电磁场辅助场；

[0016] S8.增加时间步，并重复步骤S5～S7，直至完成整个时间步的模拟；

[0017] S9.数据输出；

[0018] S10.重复步骤S3～S9，直至所有激发完成，从而完成DGTD正演，得到最终的雷达剖面。

[0019] 步骤S2所述的设置基函数的阶次和节点，具体为采用如下步骤设置基函数的阶次和节点：

[0020] A.基函数设置为正交的Legendre多项式；

[0021] B.二维情况下，采用如下算式作为标准三角形定义的N阶基函数：

[0022]

[0023] 式中 为标准三角形I内节点的横坐标；η为标准三角形I内节点的纵坐标；(i,j)≥0，且i+j≤N； 为n阶Jacobi多项式，且当α＝β＝0时， 为Legendre多
项式；N为基函数阶次；a和b为Legendre多项式的变量，且 b＝η。

[0024] C.采用如下算式作为扭曲函数ω(r)，并利用等间距节点和LGL(Legendre‑Guass‑Lobatto)节点来产生新的插值节点：

[0025]

[0026] 式中NP为节点的个数，且在二维情况下 为第i个LGL节点； 为第m个三角单元的第i个等距节点， 为基于 的Lagrange多项式；r为
新插值节点的位置。

[0027] 步骤S3所述的设置激励源和接收点的位置，具体为采用如下步骤设置激励源和接收点的位置：

[0028] 在二维TM模式下，采用如下算式表示电流源项积分：

[0029]

[0030] 式中Ωm为加载电流源的所属单元；Jz为加载的线电流源，且在z方向上加载线电流源时Jz＝f(t)δ(x‑x0,y‑y0)；f(t)为延时的雷克子波的波函数；上标T表示向量的转置； 为
基函数在此单元节点上的值所组成的向量；(x0,y0)为源点的坐标；δ(x,y)为狄拉克函数。

[0031] 步骤S4所述的计算正演所需要的参数，并开启GPU并行计算，具体为在进行激励源激发，并求取场值之前，先根据已知的数据，计算得到DGTD正演所需要的参数；同时，将所有
预先计算并得到的DGTD正演参数存入到GPU中，并进行后续场值的计算，从而对正演模型的
模拟进行加速，节约计算时间。

[0032] 步骤S5所述的开始时间积分，具体为采用五级四阶低储存Runge‑Kutta方法(LSERK)进行时间积分。

[0033] 步骤S6所述的计算通量，并更新电磁场的场值，具体为采用如下步骤计算通量并进行更新：

[0034] a.通量采用迎风通量，并采用如下算式计算迎风通量：

[0035]

[0036] 式中 为数值通量；上标‑表示单元内部信息；上标+表示单元外部信息；为节点处的场值； 为边界上的法线向量；α为保证整体稳定性的系数且取值范围
为[0,1]，同时当α＝0时数值通量 转换为迎风通量；

[0037] b.采用如下算式进行DGTD更新

[0038]

[0039] 式中ε为介电常数；σ为电导率；M为单元的质量矩阵且 Ez为z轴方向的电场值；Sx为x轴方向的单位刚性矩阵且 Sy为y轴方向的单位
刚性矩阵且 Hx为x轴方向的磁场值；Hy为y轴方向的磁场值；F为边界
单元矩阵且 κh、κe、ve和vh均为通量系数，且
为Ωm边界处
的外法向向量在x方向上的投影； 为Ωm边界处的外法向向量在y方向上的投影； 为基函
数；上标+表示单元外部信息。

[0040] 步骤S7所述的在UPML区域更新电磁场辅助场，具体为采用如下步骤进行更新：

[0041] (1)采用如下算式表示时间域上引入了辅助场P、Q的UPML方程：

[0042]

[0043] 式中上标‑1表示逆运算；κ、a、b、c和d均为对角张量，且κ＝diag(κx,κy,κz)，其中κi为改善PML对表面
波的吸收的参数，i取值为x，y或z；σi表示PML层内i方向的电导率，i取值为x，y或z；ε0为真空
中的介电常数；

[0044] (2)在二维TM模式下，只有Hx、Hy、Ez三个分量，辅助电场P中Px＝Py＝0；辅助磁场Q中Qz＝0；

[0045] (3)采用DGTD算法来求解步骤(1)中的算式，得到二维TM模式下UPML区域的半离散强格式更新方程如下所示：

[0046]

[0047] 式中 Pz为辅助电场P的z轴分量；
Qx为辅助磁场Q的x轴分量；
Qy为辅助磁场Q的y轴分量。

[0048] 本发明提供的这种GPU探地雷达复杂介质DGTD正演方法，将传统的DGTD正演方法与非结构化网格相结合，并利用GPU并行策略进行加速，提高计算效率；通过引入半离散的
强格式而不是求解大型的刚度矩阵，并使用Runge‑Kutta时间积分方案，本发明方法可以有
效地克服内存不足的问题，进而解决了不稳定问题；此外，采用UPML进行了扩展以匹配有损
介质，并用作吸收边界条件来模拟开放空间；最后，本发明方法还探究了网格大小和基函数
阶次对所本发明方法的建模精度的详细影响，通过选取合适的网格和阶次，可以在保证精
度的同时，提高计算效率；因此，本发明方法具有高精度、能与非结构化网格结合、可用来模
拟复杂介质正演的优点，提高了雷达正演效率，能够更好的模拟复杂地质，而且兼顾了计算
效率和计算精度。

附图说明

[0049] 图1为本发明方法的方法流程示意图。

[0050] 图2为本发明方法的一般三角形与标准三角形之间的映射示意图。

[0051] 图3为本发明方法的扭曲函数作用示意图。

[0052] 图4为本发明方法的等边三角形上优化后的节点示意图。

[0053] 图5为本发明方法的模型1的示意图。

[0054] 图6为本发明方法的模型1的不同方法和解析解之间的A‑scan对比示意图。

[0055] 图7为本发明方法的模型2的两种算法的网格剖分图。

[0056] 图8为本发明方法的模型2的FDTD、DGTD的波场快照和雷达剖面示意图。

[0057] 图9为本发明方法的模型3的网格剖分图及不同基函数阶次的节点示意图。

[0058] 图10为本发明方法的模型3的不同基函数阶次的雷达剖面示意图。

[0059] 图11为本发明方法的模型3的不同基函数阶次的A‑scan对比示意图。

[0060] 图12为本发明方法的模型3的不同计算条件下的计算时间和误差示意图。

[0061] 图13为本发明方法的二维简单模型500个时间步的GPU和CPU计算时间及加速比示意图。

[0062] 图14为本发明方法的隧道地质超前预报模型及雷达正剖面示意图。

[0063] 图15为本发明方法的隧道地质超前预报的wiggle图示意图。

[0064] 图16为本发明方法的DGTD、FETD在隧道超前预报中的A‑scan对比示意图。

具体实施方式

[0065] 如图1所示为本发明方法的方法流程示意图：本发明提供的这种GPU探地雷达复杂介质DGTD正演方法，包括如下步骤：

[0066] S1.建立DGTD正演模型，并对模型进行初始化；

[0067] S2.将步骤S1建立的正演模型进行三角剖分，并设置基函数的阶次和节点；具体为采用如下步骤设置基函数的阶次和节点：

[0068] 一般基函数是建立在标准三角形(等边直角三角形)上的。而实际上，采用delaunay来进行剖分时，剖分的三角形是一般三角形，并非标准三角形；如图2所示，便是本
发明方法的标准三角形和一般三角形之间的映射关系；

[0069] 标准三角形的定义为：

[0070]

[0071] 假设一般三角形的三个顶点为(v1,v2,v3)，这些顶点按照逆时针编号；定义其重心坐标为(λ1,λ2,λ3)；它们具有的性质为：

[0072] 0≤λi≤1,λ1+λ2+λ3＝1

[0073] 三角形的任意一点都可以表示为：

[0074]

[0075] 因为标准三角形I的三个顶点坐标分别为(‑1，‑1)，(1，‑1)和(‑1，1)；所以，I内的点可以表示为：

[0076]

[0077] 根据上述方程，可以得到：

[0078]

[0079] 映射关系则可以表示为：

[0080]

[0081] 由于高阶的简单单项基函数会使所求场值的条件数变差，精度变低，所以本发明采用正交的Legendre多项式；

[0082] 然后，二维情况下，在标准三角形I(直角等腰三角形)内定义N阶的基函数为：

[0083]

[0084] 式中 为标准三角形I内节点的横坐标；η为标准三角形I内节点的纵坐标；(i,j)≥0，且i+j≤N； 为n阶Jacobi多项式，且当α＝β＝0时 为Legendre多项式；N为
基函数阶次；a和b为Legendre多项式的变量，且 b＝η。

[0085] 同时；在有限单元法中，节点大多数为等距的。但是当选用的基函数阶数较大时，等距节点会使插值不准确；设置一个扭曲函数ω(r)，并利用等间距节点和LGL(Legendre‑
Guass‑Lobatto)节点来产生新的插值节点：

[0086]

[0087] 式中NP为节点的个数，且在二维情况下 为第i个LGL节点； 为第m个三角单元的第i个等距节点， 为基于 的Lagrange多项式；r为新
插值节点的位置。

[0088] 图3所示的便是本发明方法各条边的扭曲函数作用图；图4为本发明方法的等边三角形上优化后的节点分布图，对应基函数的阶数分别为N＝2、3、4、5、6、7；从图4中可以看
出，新设置的节点会向三个顶点聚集，相对于等距节点，更加适合插值；相对于LGL节点，它
的分布更加均匀，求解场值时会更加准确

[0089] S3.设置激励源和接收点的位置；具体为采用如下步骤设置激励源和接收点的位置：

[0090] 在DGTD算法中，因为可以设置高阶基函数，因此其节点并不全在单元边界上，有的节点会在单元内部；如果采用时域有限元算法(FETD)的方法，将激励源加载到某一边界节
点上，将会有一定的误差；所以，本发明将采用形函数的加载方式，将源等价加到这一单元
的所有节点内；

[0091] GPR一般采取延时的雷克子波f(t)作为电流源，如若在z方向上加载线电流源Jz，则：Jz＝f(t)δ(x‑x0,y‑y0)，其中(x0,y0)为源点的坐标；二维TM模式下，电流源项积分又可被
表示为：

[0092]

[0093] 式中Ωm为加载电流源的所属单元；Jz为加载的线电流源，且在z方向上加载线电流源时Jz＝f(t)δ(x‑x0,y‑y0)；f(t)为延时的雷克子波的波函数；上标T表示向量的转置， 为
基函数在此单元节点上的值所组成的向量；(x0,y0)为源点的坐标；δ(x,y)为狄拉克函数；

[0094] 因为，Jz本身具有δ函数，拥有选择性质；在源所在的单元内，将δ函数与基函数相乘，则变为源在各个节点上的形函数，如此便将源加载到了源所在单元内的各个节点上，分
配权重即为源在各个节点上的形函数；在其他单元上，δ函数均为0；

[0095] S4.计算正演所需要的参数，并开启GPU并行计算；具体为在进行激励源激发，并求取场值之前，先根据已知的数据，计算得到DGTD正演所需要的参数；同时，将所有预先计算
并得到的DGTD正演参数存入到GPU中，并进行后续场值的计算，从而对正演模型的模拟进行
加速，节约计算时间；

[0096] 具体实施时，开展一次正演，除了激发点与接受点的位置不同外，所有关于网格的信息、正演参数都是相同的；因此在进行激励源激发，求取场值之前，先根据已经剖分好的
网格、设置好的基函数和节点等信息，算出正演所需要的参数，可以节约大量时间，避免重
复计算；此外，将所需要的参数存入到GPU中，进行求取场值的计算，可以对正演模型的模拟
进行加速，节约计算时间；

[0097] S5.开始时间积分；具体为采用五级四阶低储存Runge‑Kutta方法(LSERK)进行时间积分；

[0098] 正演模拟的实质是，求取模型中每个时刻每个节点的场值；因此，时间离散对于正演模拟也非常重要；本发明采取的是五级四阶低储存格式Runge‑Kutta(LESRK)来进行时间
上的模拟，具体公式如下：

[0099]

[0100] 其中， 为第n个时间步的场值， 为第n+1个时间步的场值，Δt为时间步长，(i) (i)
p 、k 为中间变量，αi、βi、γi为常数，表1为本发明方法的LESRK方法中的系数值；

[0101] 表1 LESRK方法中的系数值示意表

[0102]

[0103] 由于二阶Runge‑Kutta不具有高阶收敛性，会在计算中产生误差，而五级四阶低储存的Runge‑Kutta方法不仅稳定性高，还可以降低内存消耗，提高计算效率；

[0104] S6.计算通量，并更新电磁场的场值；具体为采用如下步骤计算通量并进行更新：

[0105] 由电磁波传播理论可知，雷达波在有耗媒介中的传播遵循Maxwell方程组；时间域二维TM波的标量波动方程可以表示为：

[0106]

[0107] 其中其中ε为介电常数(F/m)，μ为磁导率(H/m)，t为时间(s)，σ是电导率(S/m)，Hx，Hy分别代表x、y
2
方向上的电场分量(V/m)，Ez是z方向上的磁场分量(A/m)，Jz是z方向上的电流密度(A/m)；
▽·表示散度运算符。

[0108] 采用DGTD法求解二维GPR波动方程；将整个区域Ω划分为M个三角单元，每个单元内部的电性参数相同，则整体解U用分片式局部解Um逼近；在第m个三角单元Ωm中，将上式乘
以试函数 并在单元内部积分，可以得到

[0109]

[0110] 经过分部积分，上式可以写为

[0111]*

[0112] 其中 指的是Ωm的边界， 是Ωm边界处的外法向向量，(F(Um)) 是数值通量；由于每个单元相邻边界处的场值不相同，所以需要引入数值通量将两个相邻单
元的边界值统一起来；上式为DGTD算法的弱解形式，再进行一次分部积分，则可以得到强格
式方程：

[0113]

[0114] 相对于弱解形式方程，强格式不需要试函数 光滑。在边界上定义迎风数值通量矢量为

[0115]

[0116] 其中E和H为当前单元的电磁场，上标(+)表示相邻单元格的值，(*)代表数值通量，κe、vh、κh、ve为通量系数，具体的表达式为：
其中，Z和Y分别代表局部阻抗和局部导纳，εr为介质
的相对介电常数，μr为相对磁导率；

[0117] 在TM模式下，边界处的迎风通量可以写为：

[0118]

[0119] 在单元Ωm内部，将场值函数用线性插值来近似：T
其中， Ez＝(Ez,1,Ez,2,...,Ez,NP) ，Hx＝(Hx,1,Hx,2,...,
T T
Hx,NP) ，Hy＝(Hy,1,Hy,2,...,Hy,NP) ；NP为此单元内节点的个数；

[0120] 同时，通量采用迎风通量，并采用如下算式计算迎风通量：

[0121]

[0122] 式中 为数值通量；上标‑表示单元内部信息；上标+表示单元外部信息；为节点处的场值；为边界上的法线向量；α为保证整体稳定性的系数且取值范围
为[0,1]，同时当α＝0时数值通量 转换为迎风通量；

[0123] 然后，采用如下算式进行DGTD更新：

[0124]

[0125] 式中ε为介电常数；σ为电导率；M为单元的质量矩阵且 Ez为z轴方向的电场值；Sx为x轴方向的单位刚性矩阵且 Sy为y轴方向的单位
刚性矩阵且 Hx为x轴方向的磁场值；Hy为y轴方向的磁场值；F为边界
单元矩阵且 κh、κe、ve和vh均为通量系数； 为Ωm边界处的外法向向量在x
方向上的投影； 为Ωm边界处的外法向向量在y方向上的投影； 为基函数；上标+表示单
元外部信息；

[0126] S7.在UPML区域更新电磁场辅助场；具体为采用如下步骤进行更新：

[0127] 实际上的电磁场传播区域是无限的；在计算机上进行正演模拟时，只能计算有限的区域，因此需要在边界处设置完全匹配层，以免造成人为的虚假反射；

[0128] (1)根据电磁波场理论，采用如下算式表示时间域上引入了辅助场P、Q的UPML方程：

[0129]

[0130] 式中上标‑1表示逆运算；κ、a、b、c和d均为对角张量，且κ＝diag(κx,κy,κz)，其中κi为改善PML对表面
波的吸收的参数，i取值为x，y或z；σi表示PML层内i方向的电导率，i取值为x，y或z；ε0为真空
中的介电常数；

[0131] (2)在二维TM模式下，只有Hx、Hy、Ez三个分量，辅助电场P中Px＝Py＝0；辅助磁场Q中Qz＝0；

[0132] (3)采用DGTD算法来求解步骤(1)中的算式，得到二维TM模式下UPML区域的半离散强格式更新方程如下所示：

[0133]

[0134] 式中 Pz为辅助电场P的z轴分量； 为辅
助磁场Q的x轴分量； Qy
为辅助磁场Q的y轴分量；

[0135] S8.增加时间步，并重复步骤S5～S7，直至完成整个时间步的模拟；

[0136] 上述的步骤推导了在一个时间步上的DGTD算法更新方程，在步骤S5给出了时间积分的相关公式；重复上述步骤S5～S7，即可求得模拟区域内每个时间步上的每个节点的场
值，以此来模拟波场在时域的进程；

[0137] S9.数据输出；主要输出在接收点处的每个时间步的波场值，这便是雷达剖面中的一道数据；

[0138] S10.重复步骤S3～S9，直至所有激发完成，从而完成DGTD正演，得到最终的雷达剖面。

[0139] 以下结合一个实施例，对本发明方法进行进一步说明：

[0140] 与FETD算法对比：

[0141] 如图5所示为本发明方法的模型1的均匀模型示意图：三角形为激励源的位置，圆形为接受点的位置。

[0142] 该模型的大小为1m×1m，其相对介电常数为6，电导率为6mS/m，相对磁导率为1。采用主频为900MHz的雷克子波作为激励源，其位置处于正中心，坐标为(0m，0m)，接收点的位
置为(0.2m，0.2m)。时间步长为0.01ns，模拟时间为10ns。DGTD使用具有二阶基函数的粗网
格进行仿真(3200个单元和1681个节点)。粗网格FETD的网格剖分与DGTD算法相同。此外，还
采用细网格的FETD方法进行仿真(80000个单元和40401个节点)。将这三个结果与解析解进
行比较，并计算误差。

[0143]

[0144] 其中utest是由DGTD算法或FETD算法模拟的波场值，而uREF是解析解的波场值。

[0145] 所有测试均在Windows 10环境中进行。CPU环境是Inter(R)Xeon(R)Platinum 8168CPU@2.70GHz，GPU环境为NVIDIA Quadro P6000。

[0146] 如图6所示，便是本发明方法的模型1的不同方法和解析解之间的A‑scan对比图。表2给出了本发明方法的模型1的不同算法之间的计算条件和精度。

[0147] 表2本发明方法的模型1的不同算法之间的计算条件和精度示意表

[0148] 算法 时间(s) 单元数 自由度 误差‑6
DGTD 20.58 3200 19200 3.16×10
‑3
粗网格FETD 3.02 3200 1681 1.1×10
‑5
细网格FETD 120.53 80000 40401 2.9×10

[0149] 分析图6和表2可知，DGTD算法模拟结果与解析解相吻合，其误差仅为3.16×10‑6。‑3
但是粗网格的FETD算法误差达到1.1×10 ，波形上有较大的波动。细网格的FETD结果与解
析解基本吻合，但是比DGTD算法模拟花费更多的时间。这表明与FETD相比，DGTD不仅对网格
的依赖程度较低，可以使用高阶基函数来补充网格的“不足”，而且精度更高，计算效率显著
增加。

[0150] 与FDTD算法对比：

[0151] 为了说明DGTD算法相对于FDTD算法能更好地与非结构网格结合，建立一个具有粗糙倾斜界面的模型。

[0152] 如图7所示为本发明方法的模型2的网格剖分图：模型大小为20m×10m，三角形和圆形分别表示激励源和接收点的位置；图7(a)为FDTD算法的剖分网格；图7(b)为DGTD算法
的剖分网格。

[0153] 该模型的相对磁导率为1，电导率为0，上半部分的相对介电常数为3，下半部分的相对介电常数为9。采用主频为200MHz的雷克子波在顶界面的中心处来进行激发，在0.2m的
深度处每隔0.2m设置一个接收天线，总共接收100道数据，模拟时间为120ns。FDTD的网格数
为400×200，空间步长为0.05m，时间步长为0.08ns。因为FDTD剖分网格是结构化的，所以分
界处具有小型的锯齿，与实际的界面不符；采用DGTD算法的基函数阶次为N＝3，模拟区域总
单元数为7490，节点数为3847，时间步长为0.04ns。由于DGTD算法的剖分为非结构化网格，
所以分界线贴合实际情况。

[0154] 如图8所示为本发明方法的模型2的FDTD、DGTD算法的波场快照和雷达剖面图：图8(a)、(c)、(e)分别为FDTD算法的48ns、56ns和64ns的波场快照；图8(b)、(d)、(f)分别为DGTD
算法的48ns、56ns和64ns的波场快照；图8(g)和(h)分别为FDTD和DGTD算法的雷达剖面。

[0155] 分析图8可知，在图8(c)和图8(e)中，黑色椭圆所圈起的部分，便是由于FDTD算法的直交网格所引起的细小毛刺波，从而导致整个波表面存在振荡现象。在图8(g)中，可以看
到箭头指向的位置为结构化网格所引起的多次波。在DGTD算法的波场快照和雷达剖面中，
由于非结构化网格，波形平滑且无杂波。两种算法的对比表明，DGTD可以与非结构化网格相
结合，可适用于复杂模型的仿真，其精度更高。

[0156] 本发明方法的优势如下：

[0157] DGTD算法可以采用高阶基函数、大网格来进行模拟，也可采用低阶基函数、小网格来进行模拟。本发明详细阐述了基函数阶次和网格大小对DGTD算法建模精度的影响，并给
出了一个经验公式，在保证计算精度的同时，尽可能地提高计算效率；利用GPU加速策略进
行计算，可以节约计算时间；将DGTD算法应用到地质隧道预报中，证明DGTD算法对复杂模型
的适应性

[0158] 1)调节参数优化提高计算精度的优势

[0159] 基函数阶次对精度的影响

[0160] 为了说明基函数的阶次对DGTD算法精度的影响，建立了图9所示的模型3。图9为本发明方法的模型3的网格剖分图及不同基函数阶次的节点图：图9(a)为具有两个圆形异常
体的模型示意图，模拟区域大小为1m×1m，三角形和圆形分别代表激励源和接收器的位置。
图9(b)、(c)、(d)和(e)是图9(a)中虚拟框内的不同基函数阶次节点图，分别对应于3阶、4
阶、5阶和6阶基函数，圆点为节点的位置。

[0161] 该模型的背景介质的相对介电常数为9，电导率为1mS/m，相对磁导率为1；存在两个异常体，左边的空洞介质相对介电常数为1，电导率为0，圆心位置为(‑0.25m，0.1m)，半径
为0.1m；右边为混凝土介质的空洞，其相对介电常数为6，电导率为0.01mS/m，圆心的位置为
(0.25m，0.1m)，半径为0.1m。整个模拟区域的单元数为2594。激励源为100MHz的雷克子波。
发射天线位于顶界面的中心处，接收天线在‑0.45m深度处，从‑0.49m水平位置开始，每隔
0.01m设置一个，共接收100道数据；模拟时间为20ns。

[0162] 图10为本发明方法的模型3不同基函数阶次的雷达剖面图：图10(a)：基函数阶次N＝3；图10(b)：N＝4；图10(c)：N＝5；图10(d)：N＝6。

[0163] 分析图10可知，随着基函数阶次提高，每个单元内部的节点数增多，雷达剖面的信噪比也越来越高。图10(a)中，因为基函数的阶次低，节点少，所以在波的传播过程中出现了
振荡现象，虽然能看出异常体的存在，但是信噪比低，有较多杂波。相比较于图10(a)，图10
(b)中雷达剖面清晰度有所提高，杂波有所减少，异常体所产生的绕射波形更为明显。图10
(c)与图10(d)中，因为基函数阶次高，节点数多，雷达剖面更为清晰，信噪比更高。

[0164] 图11为本发明方法的模型3的不同基函数阶次的A‑scan对比图，皆为图10剖面图的中心道。分析图11可知，从整体上看，它们的数据保持一致，但是将黑色虚线框中的数据
放大，在40～60ns的时间内，基函数的阶次为3时，起伏较大，有一定的误差；N＝4时，也有一
定的起伏，但相对N＝3时较小；N＝5、N＝6的线条基本没有起伏，符合实际情况，这也反映出
同一网格下，基函数的阶次越高，精度越高。

[0165] 网格大小对精度的影响

[0166] DGTD算法正演计算过程中，计算效率与计算精度由网格的大小及基函数的阶次综合影响决定，计算的宗旨是在满足精度要求的前提下，尽可能地提高计算效率。通过大量实
验，表3给出了本发明方法的网格大小与基函数阶次建议对应的关系。其中，d为网格大小，λ
为激励源在背景介质处的波长，N为基函数阶次，网格大小与阶次的关系为d/N≈λ/15。

[0167] 表3本发明方法的网格大小与基函数阶次建议对应的关系示意表

[0168] 网格大小(d) 基函数阶次(N) d/Nλ/14 1 λ/14
λ/7 2 λ/14
λ/5 3 λ/15
λ/4 4 λ/16
λ/3 5 λ/15
2λ/5 6 λ/15
λ/2 7 λ/14
3λ/5 8 3λ/40

[0169] 为了寻求最好的网格大小、基函数阶次与计算精度、正演效率之间的关系，根据表3，设置不同的网格大小，采用不同阶次的基函数来模拟图9中的模型，网格信息和计算效率
如表4所示。为了比较误差，将极细网格的FDTD算法模拟效果作为标准(400×400网格)。误
差函数在模型1中已给出。

[0170] 表4网格信息和计算效率示意表

[0171]

[0172] 图12为本发明方法的模型3的不同条件下的计算时间和误差：柱形代表时间；带点的虚线表示计算的自由度；带有正方形的实线表示误差。

[0173] 分析图12和表4可知，不同阶次下的误差都在4‑6×10‑5这个区间，证明整体模拟精度基本在同一水平。但当基函数为1次或2次时，由于计算自由度大，导致计算时间急剧增
加；当基函数阶次提高时，计算自由度逐渐减小，计算效率显著提高，当基函数大于3阶后，
误差变化不大，计算效率呈缓慢提高。由此可见，用大网格高阶次来进行模拟比小网格低阶
次更好一些。但网格剖分过大，可能会导致对模型拟合不好，损失边角信息，从而导致误差。
因此，开展DGTD算法正演时，应该首先根据精度要求选取合适的网格大小，再根据网格大小
来选取合适的基函数阶次。

[0174] 2)GPU加速策略优势

[0175] 与FETD方法不同，DGTD方法的值在单元边界处是不连续的，将引入更多的自由度。在提高准确性的同时，它将占用更多的计算机内存并降低效率。因此，有必要采取措施来提
高效率。近年来，随着计算机技术的发展，使用GPU来提高计算效率非常流行，这也是提高
DGTD算法计算效率的有力工具。GPU和CPU的体系结构不同。CPU适用于处理少量复杂任务，
而GPU适用于处理大量简单任务。为了比较GPU并行和传统CPU串行的特性，本发明选择了一
个线源位于计算域中心的2D简单模型，对比了两种情况下的DGTD的仿真效率。

[0176] 图13为本发明方法的500个时间步下的GPU与CPU计算时间及加速比图：横坐标为计算自由度的开方，左纵坐标表示计算时间，右纵坐标为加速比。

[0177] 分析图13可知，当计算量比较小，在20000以下时，GPU计算时间要稍大于CPU计算时间，大部分时间用来传输数据；随着数据量的增大，GPU和CPU的计算时间都在增大，但GPU
计算时间为线性增长，而CPU计算时间为指数增加，最后加速比稳定在10左右。

[0178] 3)复杂模型适用性优势

[0179] 图14为本发明方法的隧道地质超前预报模型的网格剖分图及其正演剖面：图14(a)为16m×12.5m的模型，点和三角形分别代表了激励源和接收点的位置；图14(b)为雷达
剖面图。

[0180] 该模型是针对含水围岩裂隙和空洞两类病害而建立的。背景介质的相对磁导率为1。最上层是厚度为0.5m的空气层，相对介电常数为1，电导率为0。上方浅黄色的区域为黏土
覆盖层，其相对介电常数为5，电导率为2mS/m；下方土黄色区域为围岩，其相对介电常数为
9，电导率为0.5mS/m；中间灰色带区域为断层，其相对介电常数为16，电导率为10mS/m；左侧
蓝色区域为围岩破碎带，内部充满水，其相对介电常数为81，电导率为1.8mS/m；右侧白色区
域为空气溶洞，其相对介电常数为1，电导率为0。计算域由10555个单元和5397个节点组成，
计算的自由度为158325。采用主频为100MHz的雷克子波作为激励源，时间步长为0.01ns，模
拟时间为240ns。收发天线距地表0.02m，收发距为0.2m，自左向右同步移动，采样间隔为
0.15m，共记录100道数据。采用基函数阶次N＝4来进行模拟。GPU计算时间为179.89分钟，
CPU计算时间为1007.7分钟，加速比为5.6。

[0181] 分析图14(b)可知，在40到60ns之间有一条明显的波形，这是两层介质的分界面。在左下方，有大量不规则的抛物线状反射信号，这是由图14(a)左侧的细小岩石破裂带引起
的。中间的黑色实线勾勒出了断层的位置，与图14(a)相对应。在右侧上方，由椭圆形圈出的
位置，上方是来空洞的上界面反射，下面是为其下界面反射。图15是该模型的Wiggle图。

[0182] 为了验证复杂模型的DGTD算法的准确性，还使用FETD算法来正演模拟此模型。整个区域中的单元数为43019，节点数为21768。总仿真时间为430.35分钟。图16为图14(a)中
红色虚线处的两种算法的A‑scan比较图。由图可见，两者的整体数据是一致的。点A为直接
波，点B是上曲线界面的反射，点C是上空洞界面的反射，点D是空洞下界面的反射。从波形到
达时间和传播速度可以计算出该记录信号与每个异常体之间的距离基本相同，这说明了
DGTD在GPR复杂模型模拟中的可行性。