CN111324972B - 一种基于gpu并行的探地雷达电磁波数值模拟计算方法 - Google Patents

Abstract

本发明提供了一种基于GPU并行的探地雷达电磁波数值模拟计算方法，将电场和磁场在相同的空间网格节点和相同的时间步长上离散，通过GPU将并行计算技术和曲面共形技术应用到含有CPML吸收边界的辛计算方法中，形成了基于GPU并行加速的含有CPML吸收边界的共形辛计算方法，在提高数值模拟精度的同时，极大地提高了探地雷达电磁波数值模拟计算效率，大大节约了模拟计算的时间，能精确高效的解释实测的探地雷达数据，有助探地雷达在无损检测中应用。

Description

一种基于GPU并行的探地雷达电磁波数值模拟计算方法

技术领域

本发明涉及地球物理领域的探地雷达电磁波数值模拟技术领域，尤其涉及一种基于GPU并行的探地雷达电磁波数值模拟计算方法。

背景技术

探地雷达是利用电磁波的穿透能力来确定地下结构中物质分布的一种无损检测方法。利用探地雷达正演模拟技术，对地下目标的电磁响应进行数值模拟和分析，得到地下目标在地表的反射波形，了解地下结构中电磁波的传播规律。同时，可以加深对探地雷达实测剖面的认识，从而提高探地雷达探测数据的处理和解释精度。此外，正演模拟仍然是基于探地雷达测量信号确定地下目标结构参数反演的基础。准确高效的正演模拟对提高反演精度和速度具有重要意义。

目前，探地雷达正演模拟方法主要有有限元法、射线追踪法、伪谱时域法、快速多极子法、时域有限差分法、交替方向隐式时域有限差分法、辛欧拉法等。尽管这些方法能够精确地模拟探地雷达电磁波在地下结构中的传播，但这些计算方法在效率和精度上都存在一些不足。例如，有限元法适用于求解复杂边界条件问题或复杂介质定解问题，但其计算程序复杂，求解时间长，容易产生伪解现象。时域有限差分法(FDTD)是目前应用最为广泛的时域电磁场数值模拟技术，具有实现简单、应用范围广、适于并行计算等优点。然而，该方法的计算效率受到Courant-Friedrichs-Lewy(CFL)稳定性条件的限制，对于大尺寸或精细结构目标的计算效率不高。对于二维问题，辛欧拉计算方法只能用两个方程来描述电磁场。辛计算方法的计算效率和计算机占用率均优于FDTD方法，但仍然是一种差分计算方法，因而受到CFL稳定性条件的限制。因此，开发一种精确高效的探地雷达正演模拟计算方法迫切需要的，这对于探地雷达理论研究，探测方式的研究具有重要意义。

发明内容

针对现有技术的不足，本发明开发了一种基于GPU并行的探地雷达电磁波数值模拟计算方法。

一种基于GPU并行的探地雷达电磁波数值模拟计算方法，具体步骤如下：

步骤1、将待求解的地下结构模型区域划分为大小相同的网格，电场和磁场定义在网格正中心；在计算机的CPU上分配存储空间，用于储存二维平面地下结构离散区域的长度和宽度、离散区域的网格步长和时间步长、磁场的磁导率、电场的电导率、介电常数参数，并对所述参数初始化；同时在GPU上也分配存储空间用于存储：二维平面地下结构离散区域的长度和宽度、离散区域的网格步长和时间步长、磁场的磁导率、电场的电导率、介电常数参数，并对这些参数初始化；其中在CPU上定义地下结构模型区域的离散划分的各网格磁场和电场的结构体数组，用于在GPU上存储通过磁场和电场公式(1)并行计算得到的各时间步各网格的磁场值和电场值，所述的时间步，就是把整体的时间分为n段，每一段就是一个时间步，在第一个时间步内，计算各网格的磁场和电场，然后保存数据，在第二个时间步内再全部计算各网格的磁场和电场，在计算第二个时间步时需要第一个时间步的磁场和电场的值，依次一直计算到第n个时间步；这些并行计算得到的场值通过传输函数cudaMemcpy传输到CPU上定义的各磁场和电场的结构体数组，其中所述磁场和电场的初始值都为0；

步骤2、利用下面公式计算地下结构模型区域离散各网格的磁场和电场的值，二维辛计算方法的磁场A和电场U的计算公式如式(1)所示：

式中，σ表示的电导率，ε表示介电常数，μ表示磁导率，Δt为模拟计算的时间步长，Δx为二维地下结构模型区域离散网格在x方向的长度，Δy为二维地下结构模型区域离散网格在y方向的长度，i为二维平面地下结构模型模拟区域离散网格在y方向的编号，j为二维平面地下结构模型模拟区域离散网格在x方向的编号，

表示网格(i，j)处磁场A第n个时间步长时的值，

表示网格(i，j)处电场U第n个时间步长时的值，

表示网格(i，j)处磁场A第n+1个时间步长时的值，

表示网格(i，j)处电场U第n+1个时间步长时的值；

步骤3、在模拟区域网格上加载探地雷达脉冲，然后使用GPU上设置的线程并行计算公式(1)中第一式的磁场A，得到各网格处的磁场A的值；

步骤4、利用步骤3计算得到的磁场A的值，在GPU上利用设置的并行线程计算公式(1)第二式中的电场U，得到各网格处的电场U的值；最后进行电场U的边界处理，采用坐标伸缩完全匹配层吸收边界条件，二维辛计算方法的磁场U的坐标伸缩完全匹配层吸收边界条件的计算公式如式(2)所示：

式中κ_x为x方向的有效延伸因子，κ_y为y方向的有效延伸因子，其中CA和CB为公式系数，具体为

和

为中间变量函数可由如下公式计算

式中的a_x为x方向第一公式系数，b_x为x方向第一公式系数，a_y为y方向第一公式系数，b_y为y方向第二公式系数，系数具体求解式如下：

式中，式中κ_i为场值在i方向的有效延伸因子，κ_max为κ_i的最大值，本发明中κ_max取5；α_i为场值在i方向的坐标延伸因子的自由度，α_max为α_i的最大值，本发明α_max取0.008；σ_i为坐标伸缩完全匹配层吸收边界区域中在i方向的电导率，σ_max为σ_i的最大值，其计算公式为

式中的m＝4，δ为离散网格的长度；d为坐标伸缩完全匹配层吸收层厚度，d^m是d的m次方；

步骤5、将计算网格的电场和磁场的值存储到GPU全局内存，并重复步骤3和步骤4，计算二维地下结构模型区域离散每个网格不同时间步的电场A和磁场U，从时间步T＝0迭代直到最大时间步数T＝Tmax，Tmax为模拟计算的最大时间步数；

步骤6、将步骤5计算得到的各网格的电场和磁场的值从GPU设备上拷贝到CPU上，对电场和磁场的值进行存储，将电场的值导入到MATLAB中，通过MATLAB画图分析探地雷达电磁波在地下结构中的传播规律。

进一步地，所述的基于GPU并行的探地雷达电磁波数值模拟计算方法，其中：所述步骤(1)中，电场和磁场离散定义在相同的空间网格节点上，同时电场和磁场离散定义在相同的时间步长上，由电磁波传播理论可知，探地雷达遵循的麦克斯韦方程组的表达式为：

式中，H和E分别代表磁场和电场向量，我们引入矢量磁位H＝▽×A，并令E＝-U，则上述的麦克斯韦方程组为

在电磁波传播理论中,麦克斯韦方程组二维TM波的电场波表示为

式中，U_z和A_z为场组分U和A沿z方向的分量，▽²为拉普拉斯算子，应用一阶辛算法求解公式(12)得到公式(1)。

进一步地，所述的基于GPU并行的探地雷达电磁波数值模拟计算方法，其中：在二维地下结构模型划分的网格包含两种不同材质时，计算各网格的介电常数、电导率、磁导率的公式如下:

式中σ^eff表示的等效电导率，ε^eff表示等效介电常数，μ^eff表示等效磁导率，Δx和Δy是模拟区域的划分的网格的宽度和高度，S_xy1和S_xy2分别是二维地下结构模型区域中不同材料介质1和2在网格中所占的面积，ε₁、σ₁、μ₁分别为材料介质1的介电常数、电导率、磁导率，ε₂、σ₂、μ₂分别为材料介质2的介电常数、电导率、磁导率；当二维地下结构模型划分的网格包含两种不同材质时，我们利用公式(1)计算各网格电场和磁场时，所述的电导率、介电常数、磁导率换算成等效电导率、等效介电常数、等效磁导率。

进一步地，所述的基于GPU并行的探地雷达电磁波数值模拟计算方法，其中：GPU上并行计算的线程的数量对应二维地下结构模型区域离散网格的数量，每个网格由唯一的一个线程计算其场值，在处理二维地下结构模型区域时，需GPU中两个并行计算的内核函数，内核函数由模型网格计算所需的线程组成.第一个内核函数计算各网格的磁场A，第二个内核函数更新各网格的电场U。

进一步地，所述的基于GPU并行的探地雷达电磁波数值模拟计算方法，其中：探地雷达的脉冲信号为中心频率为1GHz单位振幅Ricker子波。

本发明的一种基于GPU并行的探地雷达电磁波数值模拟计算方法，结合二维辛计算方法、曲面共形技术和GPU加速技术，建立二维探地雷达电磁波在地下结构中的传播模，实现了复杂地下结构电磁响应的精细数值模拟。适用于二维复杂地下结构，采用曲面共形网格技术，将GPU并行计算技术应用到含有CPML吸收边界的辛计算方法中，形成了基于GPU并行加速的含有CPML吸收边界的辛计算方法，在保持共形辛计算方法与吸收边界的优良特性的同时，极大地提高了计算效率，大大节约了计算时间，达到了探地雷达探实测数据的处理和解释精度。本发明技术方案的实施提高了探地雷达探实测数据的处理和解释精度，为下一步探地雷达三维反演成像准确、高效地确定复杂地下结构中物质的分布提供了一个高效、准确的正演模型。

附图说明

图1为本技术方案步骤(1)所述电场和磁场的空间分布示意图；

图2为本技术方案二维并行共形辛计算方法的线程安排示意图；

图3为本技术方案基于CPU并行共形辛计算方法的具体处理流程图；

图4中心频率为1GHz的单位波幅Ricker子波波形图；

图5并行共形辛算法得到的探地雷达剖面图；

图6串行非共形辛算法得到的探地雷达剖面图。

具体实施方式

一种基于GPU并行的探地雷达电磁波数值模拟计算方法，具体步骤如下：

步骤1、将待求解的地下结构模型区域划分为大小相同的网格，电场和磁场定义在网格正中心，如图1所示；在计算机的CPU上分配存储空间，用于储存二维平面地下结构离散区域的长度和宽度、离散区域的网格步长和时间步长、磁场的磁导率、电场的电导率、介电常数参数，并对所述参数初始化；同时在GPU上也分配存储空间用于存储：二维平面地下结构离散区域的长度和宽度、离散区域的网格步长和时间步长、磁场的磁导率、电场的电导率、介电常数参数，并对这些参数初始化；其中在CPU上定义地下结构模型区域的离散划分的各网格磁场和电场的结构体数组，用于在GPU上存储通过磁场和电场公式(1)并行计算得到的各时间步各网格的磁场值和电场值，所述的时间步，就是把整体的时间分为n段，每一段就是一个时间步，在第一个时间步内，计算各网格的磁场和电场，然后保存数据，在第二个时间步内再全部计算各网格的磁场和电场，在计算第二个时间步时需要第一个时间步的磁场和电场的值，依次一直计算到第n个时间步；这些并行计算得到的场值通过传输函数cudaMemcpy传输到CPU上定义的各磁场和电场的结构体数组，其中所述磁场和电场的初始值都为0；

步骤2、利用下面公式计算地下结构模型区域离散各网格的磁场和电场的值，二维辛计算方法的磁场A和电场U的计算公式如式(1)所示：

式中，σ表示的电导率，ε表示介电常数，μ表示磁导率，Δt为模拟计算的时间步长，Δx为二维地下结构模型区域离散网格在x方向的长度，Δy为二维地下结构模型区域离散网格在y方向的长度，i为二维平面地下结构模型模拟区域离散网格在y方向的编号，j为二维平面地下结构模型模拟区域离散网格在x方向的编号，

表示网格(i，j)处磁场A第n个时间步长时的值，

表示网格(i，j)处电场U第n个时间步长时的值，

表示网格(i，j)处磁场A第n+1个时间步长时的值，

表示网格(i，j)处电场U第n+1个时间步长时的值；

步骤3、在模拟区域网格上加载探地雷达脉冲，然后使用GPU上设置的线程并行计算公式(1)中第一式的磁场A，得到各网格处的磁场A的值；

步骤4、利用步骤3计算得到的磁场A的值，在GPU上利用设置的并行线程计算公式(1)第二式中的电场U，得到各网格处的电场U的值；最后进行电场U的边界处理，采用坐标伸缩完全匹配层吸收边界条件，二维辛计算方法的磁场U的坐标伸缩完全匹配层吸收边界条件的计算公式如式(2)所示：

式中κ_x为x方向的有效延伸因子，κ_y为y方向的有效延伸因子，其中CA和CB为公式系数，具体为

和

为中间变量函数可由如下公式计算

式中的a_x为x方向第一公式系数，b_x为x方向第一公式系数，a_y为y方向第一公式系数，b_y为y方向第二公式系数，系数具体求解式如下：

式中，式中κ_i为场值在i方向的有效延伸因子，κ_max为κ_i的最大值，本发明中κ_max取5；α_i为场值在i方向的坐标延伸因子的自由度，α_max为α_i的最大值，本发明α_max取0.008；σ_i为坐标伸缩完全匹配层吸收边界区域中在i方向的电导率，σ_max为σ_i的最大值，其计算公式为

式中的m＝4，δ为离散网格的长度；d为坐标伸缩完全匹配层吸收层厚度，d^m是d的m次方；

步骤5、将计算网格的电场和磁场的值存储到GPU全局内存，并重复步骤3和步骤4，计算二维地下结构模型区域离散每个网格不同时间步的电场A和磁场U，从时间步T＝0迭代直到最大时间步数T＝Tmax，Tmax为模拟计算的最大时间步数；

步骤6、将步骤5计算得到的各网格的电场和磁场的值从GPU设备上拷贝到CPU上，对电场和磁场的值进行存储，将电场的值导入到MATLAB中，通过MATLAB画图分析探地雷达电磁波在地下结构中的传播规律。

图2为本技术方案二维并行共形辛计算方法的线程安排示意图；图3为本技术方案基于CPU并行共形辛计算方法的具体处理流程图；

进一步地，所述步骤(1)中，电场和磁场离散定义在相同的空间网格节点上，同时电场和磁场离散定义在相同的时间步长上，由电磁波传播理论可知，探地雷达遵循的麦克斯韦方程组的表达式为：

式中，H和E分别代表磁场和电场向量，我们引入矢量磁位H＝▽×A，并令E＝-U，则上述的麦克斯韦方程组为

在电磁波传播理论中,麦克斯韦方程组二维TM波的电场波表示为

式中，U_z和A_z为场组分U和A沿z方向的分量，▽²为拉普拉斯算子，应用一阶辛算法求解公式(12)得到公式(1)。

其中：在二维地下结构模型划分的网格包含两种不同材质时，计算各网格的介电常数、电导率、磁导率的公式如下:

式中σ^eff表示的等效电导率，ε^eff表示等效介电常数，μ^eff表示等效磁导率，Δx和Δy是模拟区域的划分的网格的宽度和高度，S_xy1和S_xy2分别是二维地下结构模型区域中不同材料介质1和2在网格中所占的面积，ε₁、σ₁、μ₁分别为材料介质1的介电常数、电导率、磁导率，ε₂、σ₂、μ₂分别为材料介质2的介电常数、电导率、磁导率；当二维地下结构模型划分的网格包含两种不同材质时，我们利用公式(1)计算各网格电场和磁场时，所述的电导率、介电常数、磁导率换算成等效电导率、等效介电常数、等效磁导率。

另外，GPU上并行计算的线程的数量对应二维地下结构模型区域离散网格的数量，每个网格由唯一的一个线程计算其场值，在处理二维地下结构模型区域时，需GPU中两个并行计算的内核函数，内核函数由模型网格计算所需的线程组成.第一个内核函数计算各网格的磁场A，第二个内核函数更新各网格的电场U。

此外，所述探地雷达的脉冲信号为中心频率为1GHz单位振幅Ricker子波。

下面将本发明应用于含有裂缝、空洞等病害的三层路面结构模型中来进一步对本发明详细说明，

面层裂缝、基层空洞是道路工程中常见的病害问题，考虑面层裂缝、基层空洞情况，建立了三层路面结构模型。该模型第一层为厚度30cm的沥青面层，在路面模型正中间有一条宽度为1cm的裂缝，第二层为厚度50cm的水泥基层，在路面左右分别有一个30cm×40cm矩形空洞和一个半径为5cm的圆形空洞，第三层为厚度100cm的土底基层。假定所有材料的相对磁导率μ都为1，各层的介电参数如表1所示：

表1各介质的介电参数

探地雷达脉冲信号选取为接近实际工程的中心频率为1GHz单位振幅Ricker子波。依据路面结构模型的大小，空间间隔增量和时间步长分别采用Δx＝Δy＝0.005m和dt＝0.01ns。在圆形空洞边界处，离散网格包括两种材料时，我们采用共形技术，应用等效电导率、等效介电常数、等效磁导率代替电导率、介电常数、磁导率。

我们设置了80个探地雷达脉冲信号的发射点和接收点，发射点和接收点之间的距离为10cm，且它们离沥青面层的距离为0.5cm，沿雷达前进的方向每隔1cm发射和接收一次，图5和图6是使用基于GPU并行共形辛方法和串行非共形辛方法模拟该路面模型时获得的探地雷达剖面图。其中基于GPU并行共形辛方法需要0.1771h，串行非共形辛方法需要2.7357h。模拟研究结果表明，基于GPU并行共形辛方法可以大大减少计算时间，相比串行非共形辛方法节约93.5％的计算时间。从图5和6可以看出，三层路面模型的层界面清晰可见，裂缝处出现明显断界面曲线，在矩形和圆形空洞处出现明显双曲线绕射波，与实测的路面结构的探地雷达剖面特征吻合。结果表明，本发明可以高效的模拟探地雷达地磁波在不同路面结构中的传播，实现探地雷达电磁波在路面结构无损检测结果的准确分析和解释。为下一步探地雷达三维反演成像准确、高效地确定复杂地下结构中物质的分布提供了一个高效、准确的正演模型。本发明可以高效精确的解释实测的探地雷达数据，有助探地雷达在地下结构无损检测中的应用。

Claims (5)

1.一种基于GPU并行的探地雷达电磁波数值模拟计算方法，其特征在于，包括如下步骤：

步骤1、将待求解的地下结构模型区域划分为大小相同的网格，电场和磁场定义在网格正中心；在计算机的CPU上分配存储空间，用于储存二维平面地下结构离散区域的长度和宽度、离散区域的网格步长和时间步长、磁场的磁导率、电场的电导率、介电常数参数，并对所述参数初始化；同时在GPU上也分配存储空间用于存储：二维平面地下结构离散区域的长度和宽度、离散区域的网格步长和时间步长、磁场的磁导率、电场的电导率、介电常数参数，并对这些参数初始化；其中在CPU上定义地下结构模型区域的离散划分的各网格磁场和电场的结构体数组，用于在GPU上存储通过磁场和电场公式(1)并行计算得到的各时间步各网格的磁场值和电场值，所述的时间步，就是把整体的时间分为n段，每一段就是一个时间步，在第一个时间步内，计算各网格的磁场和电场，然后保存数据，在第二个时间步内再全部计算各网格的磁场和电场，在计算第二个时间步时需要第一个时间步的磁场和电场的值，依次一直计算到第n个时间步；这些并行计算得到的场值通过传输函数cudaMemcpy传输到CPU上定义的各磁场和电场的结构体数组，其中所述磁场和电场的初始值都为0；

步骤2、利用下面公式计算地下结构模型区域离散各网格的磁场和电场的值，二维辛计算方法的磁场A和电场U的计算公式如式(1)所示：

式中，σ表示的电导率，ε表示介电常数，μ表示磁导率，Δt为模拟计算的时间步长，Δx为二维地下结构模型区域离散网格在x方向的长度，Δy为二维地下结构模型区域离散网格在y方向的长度，i为二维平面地下结构模型模拟区域离散网格在y方向的编号，j为二维平面地下结构模型模拟区域离散网格在x方向的编号，

表示网格(i，j)处磁场A第n个时间步长时的值，

表示网格(i，j)处电场U第n个时间步长时的值，

表示网格(i，j)处磁场A第n+1个时间步长时的值，

表示网格(i，j)处电场U第n+1个时间步长时的值；

步骤3、在模拟区域网格上加载探地雷达脉冲，然后使用GPU上设置的线程并行计算公式(1)中第一式的磁场A，得到各网格处的磁场A的值；

步骤4、利用步骤3计算得到的磁场A的值，在GPU上利用设置的并行线程计算公式(1)第二式中的电场U，得到各网格处的电场U的值；最后进行电场U的边界处理，采用坐标伸缩完全匹配层吸收边界条件，二维辛计算方法的磁场U的坐标伸缩完全匹配层吸收边界条件的计算公式如式(2)所示：

式中κ_x为x方向的有效延伸因子，κ_y为y方向的有效延伸因子，其中CA和CB为公式系数，具体为

和

为中间变量函数可由如下公式计算

式中的a_x为x方向第一公式系数，b_x为x方向第一公式系数，a_y为y方向第一公式系数，b_y为y方向第二公式系数，系数具体求解式如下：

式中，式中κ_i为场值在i方向的有效延伸因子，κ_max为κ_i的最大值，κ_max取5；α_i为场值在i方向的坐标延伸因子的自由度，α_max为α_i的最大值，α_max取0.008；σ_i为坐标伸缩完全匹配层吸收边界区域中在i方向的电导率，σ_max为σ_i的最大值，其计算公式为

式中的m＝4，δ为离散网格的长度；d为坐标伸缩完全匹配层吸收层厚度，d^m是d的m次方；

步骤5、将计算网格的电场和磁场的值存储到GPU全局内存，并重复步骤3和步骤4，计算二维地下结构模型区域离散每个网格不同时间步的电场A和磁场U，从时间步T＝0迭代直到最大时间步数T＝Tmax，Tmax为模拟计算的最大时间步数；

步骤6、将步骤5计算得到的各网格的电场和磁场的值从GPU设备上拷贝到CPU上，对电场和磁场的值进行存储，将电场的值导入到MATLAB中，通过MATLAB画图分析探地雷达电磁波在地下结构中的传播规律。

2.如权利要求1所述的基于GPU并行的探地雷达电磁波数值模拟计算方法，其特征在于：所述步骤(1)中，电场和磁场离散定义在相同的空间网格节点上，同时电场和磁场离散定义在相同的时间步长上，由电磁波传播理论可知，探地雷达遵循的麦克斯韦方程组的表达式为：

式中，H和E分别代表磁场和电场向量，我们引入矢量磁位

并令E＝-U，则上述的麦克斯韦方程组为

在电磁波传播理论中,麦克斯韦方程组二维TM波的电场波表示为

式中，U_z和A_z为场组分U和A沿z方向的分量，

为拉普拉斯算子，应用一阶辛算法求解公式(12)得到公式(1)。

3.如权利要求1所述的基于GPU并行的探地雷达电磁波数值模拟计算方法，其特征在于：在二维地下结构模型划分的网格包含两种不同材质时，计算各网格的介电常数、电导率、磁导率的公式如下:

式中σ^eff表示的等效电导率，ε^eff表示等效介电常数，μ^eff表示等效磁导率，Δx和Δy是模拟区域的划分的网格的宽度和高度，S_xy1和S_xy2分别是二维地下结构模型区域中不同材料介质1和2在网格中所占的面积，ε₁、σ₁、μ₁分别为材料介质1的介电常数、电导率、磁导率，ε₂、σ₂、μ₂分别为材料介质2的介电常数、电导率、磁导率；当二维地下结构模型划分的网格包含两种不同材质时，我们利用公式(1)计算各网格电场和磁场时，所述的电导率、介电常数、磁导率换算成等效电导率、等效介电常数、等效磁导率。

4.如权利要求1所述的基于GPU并行的探地雷达电磁波数值模拟计算方法，其特征在于：GPU上并行计算的线程的数量对应二维地下结构模型区域离散网格的数量，每个网格由唯一的一个线程计算其场值，在处理二维地下结构模型区域时，需GPU中两个并行计算的内核函数，内核函数由模型网格计算所需的线程组成.第一个内核函数计算各网格的磁场A，第二个内核函数更新各网格的电场U。

5.如权利要求1所述的基于GPU并行的探地雷达电磁波数值模拟计算方法，其特征在于：探地雷达的脉冲信号为中心频率为1GHz单位振幅Ricker子波。

Images