gprMax

定义

gpr是Ground Penetrating Radar（探地雷达）的缩写，gprMax是一款模拟电磁波传播的开源软件，它使用有限差分时域 （FDTD） 方法求解 3D 麦克斯韦方程组。gprMax 专为模拟探地雷达 （GPR） 而设计，但也可用于模拟许多其他应用的电磁波传播。很遗憾，目前为止它还没有图形用户界面（GUI），它的建模关键在于in文件的编写。

安装步骤

详细安装步骤参考：gprMax最详安装步骤及常见问题解答 - 知乎

安装步骤简述：

• 下载Anaconda，并配置好环境变量

  因为gprMax是基于python语言编写的，所以需要安装python及所需python包

• 安装git

  目的是使用git克隆github上的gprMax仓库，下载gprMax

• 下载gprMax包

  git clone https://github.com/gprMax/gprMax.git

  最好在D盘新建一个文件夹my_grpmax，然后右键，选择在终端打开, 然后再输入上述命令，然后一个文件夹grpMax就会被下载到新建的文件夹下。然后执行如下操作。

  cd grpMax
  conda env create -f conda_env.yml 

  指令效果是从一个环境文件 (conda_env.yml) 创建一个新的 Conda 环境。这个命令会读取 conda_env.yml 文件中的配置信息，包括依赖包及其版本，并安装这些依赖到新创建的环境中。

• 安装一个支持OpenMP的C编译器

• 安装和建立gprMax

  conda activate gprMax //激活grpMax
  python setup.py build 
  python setup.py install//安装grpMax

• 测试并验证

  python -m gprMax user_models/cylinder_Ascan_2D.in #运行并模拟一个模型
  python -m tools.plot_Ascan user_models/cylinder_Ascan_2D.out #绘制模型的A扫描

.in文件

由上述测试可知，我们必须向 gprMax 提供一个输入文件（.in文件）,该文件应包含运行 GPR 模型所需的所有信息，包括仿真场景、材料属性（介电常数、磁导率和电导率等）、几何结构（如目标体）、源（如天线）以及输出选项等，这些信息通过一条一条的指令来描述，指令和关联的参数应占据输入文件的一行，并且每行只允许一个指令。当指令需要多个参数时，应使用空格字符分隔这些参数。它是一个 ASCII 文本文件，可以使用任何文本编辑器或文字处理程序进行编写，比如，记事本，typora。

更多内容参考:输入文件命令 — gprMax 文档

基本命令

大多数命令是可选的，但有一些基本命令是构建任何模型所必需的。

• #domain

  允许您指定仿真域尺寸的大小。该命令的语法为：

  #domain: f1 f2 f3

  f1,f2,f3分别是模型在 X、Y 和 Z 方向上的大小。例如，要指定 500 x 500 x 1000 毫米的模型，请使用：#domain: 0.5 0.5 1.0

  如果建立的是二维模型，那么在第三维用一个单位网格表示，如：我建立以x,y方向的二维模型，每个网格单元边长为0.002，z就设置为0.002,（时域有限差分方法是要将对象离散成单位网格）。

  域的大小应该包含足够的体积，加上10（默认值）个单元格的 PML （完全吸收边界条件）边界

  PML 与任何对象间大约要相隔10个单元网格，域一般要足够大，才能满足建模条件。对于 2D 模型，在z方向上也应该指定一个单元格。

• #dx_dy_dz

  #dx_dy_dz: f1 f2 f3

  定义网格单元大小，也叫最小空间分辨率, 至于网格单元大小如何确定，有一定的规律，经验法则要求探地雷达的空间分辨率应为模型中最小波长的十分之一。为了确定最小波长，需要知道模型中存在的最高频率与最低速度,参考公式：

  λ=v/f

  那该如何确定最高频率和最低速度呢？

  最高频率不是雷克子波的中心频率，而是检查雷克子波的完整频谱，找到存在的最高频率。

  最低速度通常是指电磁波在仿真域中传播的最慢速度，这取决于介质的相对介电常数 ϵr和磁导率 μr。对于大多数非磁性材料，磁导率接近于 1，因此主要考虑介电常数的影响。

  电磁波在介质中的传播速度 v 可以通过以下公式计算：

  

  c 是光速，约为 3×10^8 米/秒。ϵr 是相对介电常数,μr 是相对磁导率。

• #time_window

  允许您指定所需的总模拟时间。该命令的语法为：

  #time_window: f1

  或者

  #time_window: i1

常规指令

• #title

  模型的标题

  #title: str1

  了易于辨识，最好命名为材料加模型形态，且用英文。

• #messages

  允许您控制运行 gprMax 时屏幕上显示的信息量。该命令的语法为:

  #messages: c1

• #output_dir

  允许您控制输出文件的路径。该命令的语法为：

  #output_dir: str1

  可以是绝对路径，也可以是相对于输入文件的相对路径，默认存放在输入文件所在目录。

• #num_threads

  允许您控制在运行模型时使用多少个 OpenMP 线程（通常是可用的物理 CPU 内核数），默认使用计算机上所有可用的物理 CPU 内核。gprMax 计算量最大的部分，即 FDTD 求解器循环。该命令的语法为：

  #num_threads: i1

材料指令

• #material

  用来定义一种材料。空气已经作为内置材料存在于 gprMax 中，可以使用 free_space 标识符来访问它。

  该命令的语法如下：

  #material: f1 f2 f3 f4 str1

  • f1: 相对介电常数（无量纲）。
  • f2: 电导率（单位：S/m，西门子每米）。
  • f3: 相对磁导率（无量纲）。
  • f4: 磁损耗（单位：Ω/m，欧姆每米）。
  • str1: 材料的标识符（字符串）。

• #add_dispersion_debye

  允许您基于多极 Debye 公式为已定义的材料添加色散特性,该命令的语法如下：

  #add_dispersion_debye: i1 f1 f2 ... str1

  • i1: Debye 极点的数量（整数）。
  • f1, f2, …: 每个极点对应的参数。
  • str1: 材料的标识符（字符串）

• #add_dispersion_drude

  允许您基于多极 Drude 公式为已定义的材料添加色散特性,Drude 色散通常用于金属等材料,该命令的语法如下：

  #add_dispersion_drude: i1 f1 f2 ... str1

  • i1: Drude 极点的数量（整数）。
  • f1, f2, …: 每个极点对应的参数。
  • str1: 材料的标识符（字符串）。

• #add_dispersion_lorentz

  允许您基于多极洛伦兹公式为已定义的材料添加色散特性。Lorentz 色散适用于更复杂的介质，在不同频率下表现出不同的行为,

  该命令的语法如下：

  #add_dispersion_lorentz: i1 f1 f2 f3 ... str1

  • i1: Lorentz 极点的数量（整数）。
  • f1, f2, f3, …: 每个极点对应的参数。
  • str1: 材料的标识符（字符串）。

• #soil_peplinski

  允许您使用 Peplinski 提出的土壤混合模型来创建具有真实电介质和几何特性的土壤。此模型适用于 0.3 GHz 至 1.3 GHz 的频率范围。该命令的语法如下：

  #soil_peplinski: f1 f2 f3 f4 f5 f6 str1

  • f1: 土壤中的沙子比例（无量纲）。
  • f2: 土壤中的粘土比例（无量纲）。
  • f3: 土壤的体积密度（单位：g/cm³）。
  • f4: 沙粒的密度（单位：g/cm³）。
  • f5, f6: 土壤的体积含水量范围（无量纲）。
  • str1: 土壤的标识符（字符串）。

构造几何对象（目标体）

创建一个确定目标体需要确定它的位置，大小，材质。

• #box

  用于创建一个矩形盒状物体，指定位置的同时也能确定大小。

  #box: xmin ymin zmin xmax ymax zmax material_id

  • xmin, ymin, zmin: 盒子一角的坐标 (米)。
  • xmax, ymax, zmax: 盒子对角线另一角的坐标 (米)。
  • material_id: 分配给该盒子的材料标识符。

• #cylinder

  用于创建一个圆柱体，指定位置和大小。

  #cylinder: x1 y1 z1 x2 y2 z2 r material_id

  • x1, y1, z1: 圆柱底面中心点的坐标 (米)。
  • x2, y2, z2: 圆柱顶面中心点的坐标 (米)。
  • r: 圆柱半径 (米)。
  • material_id: 分配给该圆柱的材料标识符。

• #sphere

  用于创建一个球体。

  #sphere: x y z r material_id

  • x, y, z: 球心位置的坐标 (米)。
  • r: 球体半径 (米)。
  • material_id: 分配给该球体的材料标识符。

• #fractal_box

  用于创建分形填充的盒子，通常用于模拟复杂或随机分布的目标体，如土壤中的颗粒。

  #fractal_box: xmin ymin zmin xmax ymax zmax material_id fractal_dim n

  • xmin, ymin, zmin: 分形盒子一角的坐标 (米)。

  • xmax, ymax, zmax: 分形盒子对角线另一角的坐标 (米)。

  • material_id: 分配给该分形盒子的材料标识符。

  • fractal_dim: 分形维度，决定了填充模式的复杂度。

  • n: 迭代次数，影响生成结构的精细程度。

• geometry_view

  利用这个命令可以将我们建立模型的几何信息输出到文件中，这些文件使用开源可视化工具包（.VTK）格式，可以在许多免费阅读器（如ParaView,下载对应的zip文件然后解压即可）中查看。这个命令可以用来创建模型的几个3D视图，这对于检查模型是否已经按照预期构建是很有用的。

  • x_min, y_min, z_min: 几何视图体积的左下角坐标（米），定义了视图的起点。

  • x_max, y_max, z_max: 几何视图体积的右上角坐标（米），定义了视图的终点。

  • dx, dy, dz: 网格在每个空间上的大小。通常，这些值与模型的空间离散化相同，但可以根据需要调整为更粗略（Courser）。

  • filename: 输出文件的名称。几何视图将存储在与输入文件相同的目录中。

  • mode：可以是 N（正常）或 F（精细），指定如何输出几何信息

    N (Normal): 按单元格输出几何信息，适用于大多数情况。

    F (Fine): 按单元格边缘输出几何信息，适合查看占用小体积的几何体的详细部分。注意，这种模式会生成较大的文件。

源命令和输出命令

• #waveform

  用于定义一种可被应用于源的波型，需要自定义波形类型、振幅和中心频率

  #waveform: waveform_type amplitude frequency identifier

  • waveform_type: 波形类型（如高斯脉冲 gaussian）。
  • amplitude: 波形的振幅。
  • frequency: 中心频率（单位：Hz）。确定网格大小的时候需要用到 。
  • identifier: 波形的唯一标识符。

• #excitation_file

  允许您从文件中读取激励信号，而不是使用内置波形，不过gprMax其实内置了挺多的波形的。

  #excitation_file: filename identifier

  • filename: 包含时间序列数据的文件名。
  • identifier: 激励信号的唯一标识符。

• #hertzian_dipole

  用于创建赫兹偶极子源，这是最常见的电场源之一。

  #hertzian_dipole: polarization x y z identifier [start_delay stop_delay]

  • olarization: 源的极化方向（x, y, 或 z）。极化方向决定了它发射的电磁波的电场分量如何定向
  • x, y, z: 源在模型中的坐标（米）。
  • identifier: 应与该源一起使用的波形标识符。
  • start_delay, stop_delay: 可选参数，分别是启动延迟时间和停止延迟时间（秒）。

• #magnetic_dipole

  用于创建磁偶极子源，适用于产生磁场。

  #magnetic_dipole: polarization x y z identifier [start_delay stop_delay]

  • polarization: 源的极化方向（x, y, 或 z）。
  • x, y, z: 源在模型中的坐标（米）。
  • identifier: 应与该源一起使用的波形标识符。
  • start_delay, stop_delay: 可选参数，分别是启动延迟时间和停止延迟时间（秒）。

• #voltage_source

  允许您在电场位置引入电压源。如果内阻为零，则它是硬源；否则是电阻电压源。

  #voltage_source: polarization x y z resistance identifier [start_delay stop_delay]

  • polarization: 源的极化方向（x, y, 或 z）。
  • x, y, z: 源在模型中的坐标（米）。
  • resistance: 内阻（欧姆），设置为 0 表示硬源。
  • identifier: 应与该源一起使用的波形标识符。
  • start_delay, stop_delay: 可选参数，分别是启动延迟时间和停止延迟时间（秒）。

• #transmission_line

  允许您引入一维传输线模型，通常用于激励天线。

  #transmission_line: polarization x y z characteristic_impedance identifier [start_delay stop_delay]

  • polarization: 传输线的极化方向（x, y, 或 z）。
  • x, y, z: 传输线在模型中的坐标（米）。
  • characteristic_impedance: 特征阻抗（欧姆），可以是大于零且小于自由空间阻抗（376.73 欧姆）的任何值。
  • identifier: 应与该源一起使用的波形标识符。
  • start_delay, stop_delay: 可选参数，分别是启动延迟时间和停止延迟时间（秒）。

• #rx

  用于定义接收点，保存电场和磁场分量的时间历史记录。

  #rx: x y z [identifier output_list]

  • x, y, z: 接收器在模型中的坐标（米）。
  • identifier: 接收器的唯一标识符（可选）。
  • output_list: 要输出的电场或磁场分量列表（如 Ex Ey Ez Hx Hy Hz），默认输出所有分量。

• #rx_array

  提供一种简单的方法来定义多个接收点。

  #rx_array: x_min y_min z_min x_max y_max z_max dx dy dz

  • x_min, y_min, z_min: 输出区域的左下角坐标（米）。
  • x_max, y_max, z_max: 输出区域的右上角坐标（米）。
  • dx, dy, dz: 每个方向上的增量（米），可以设置为零以防止在特定方向上出现任何输出点。

• #snapshot

  这个命令用于生成波场快照文件（VTK格式），波场快照就是在指定时刻的电磁场信息，用动图显示一段时间的电磁场信息，这个文件可以在许多免费软件（如Paraview）中查看。

  #snapshot: x_min y_min z_min x_max y_max z_max dx dy dz time_or_iteration file_name

  • x_min, y_min, z_min: 快照体积的左下角坐标（米）。

  • x_max, y_max, z_max: 快照体积的右上角坐标（米）。

  • dx, dy, dz: 网格大小，通常情况下，这与模型中网格大小是相同的。

  • time_or_iteration: 快照的时间点（秒）或迭代编号。

  • file_name: 存储快照的文件名。快照文件会自动的储存在一个目录中，该目录的名称为输入文件名称-snaps

PML 命令

在电磁波传播、声波传播等领域的数值仿真中，仿真域总是有限的，而实际物理现象往往是无限或非常大的。当波到达仿真域边界时，如果没有适当的处理，波会被反射回仿真域内部，导致不真实的干扰,PML 被设计为一个虚拟层，放置在仿真域的边界周围。它能够有效地吸收进入其中的波，使这些波不会被反射回仿真域内，从而模拟了无限空间的效果。

• #pml_cells

  允许您控制模型域六个侧面上使用的 PML 单元数量（即厚度）。请注意，PML 不会增加仿真域的实际大小；它是在定义的仿真域内实现的。

  #pml_cells: i1 [i2 i3 i4 i5 i6]

  • i1: 模型域所有侧面使用的 PML 单元数（可以设为零以完全关闭 PML），或最靠近 x 轴原点一侧的 PML 单元数。
  • i2: 最靠近 y 轴原点一侧的 PML 单元数。
  • i3: 最靠近 z 轴原点一侧的 PML 单元数。
  • i4: 最远离 x 轴原点一侧的 PML 单元数。
  • i5: 最远离 y 轴原点一侧的 PML 单元数。
  • i6: 最远离 z 轴原点一侧的 PML 单元数。

  可以将任意一个参数设为零，以关闭特定边上的 PML。

案例分析

更多内容参考文档：入门级 （2D） 模型 — gprMax documentation

我们在测试使用的user_models/cylinder_Ascan_2D.in文件内容如下：

#title: A-scan from a metal cylinder buried in a dielectric half-space
#domain: 0.240 0.210 0.002 //建立的是二维模型
#dx_dy_dz: 0.002 0.002 0.002 //指定网格的大小
#time_window: 3e-9 //时间窗口

#material: 6 0 1 0 half_space //定义一种材质half_space

#waveform: ricker 1 1.5e9 my_ricker //定义一种波my——ricker
#hertzian_dipole: z 0.100 0.170 0 my_ricker //使用这种波
#rx: 0.140 0.170 0 //定义接收器receiver的位置

//创建长方体box，左下点坐标为(0，0，0)，右上点坐标为(0.24,0.17.0.02),使用的材料是half_space
#box: 0 0 0 0.240 0.170 0.002 half_space 
//创建圆柱体,底部圆心坐标为(0.12,0.08,0)，顶部圆心坐标为(0.12,0.08,0.02),使用的材料是内置的pec
#cylinder: 0.120 0.080 0 0.120 0.080 0.002 0.010 pec

#geometry_view: 0 0 0 0.240 0.210 0.002 0.002 0.002 0.002 cylinder_half_space n

对应的2D图：

金属圆柱体将被建模为完美电导体，它作为 gprMax 中的内置材料存在,可以使用标识符pec表示

因此，唯一需要定义的材料是电介质半空间half_space,它是一种非磁性材料,可定义为

#material: 6 0 1 0 half_space

计算模拟域和最小空间分辨率的大小（单元格大小）的方法在前文指令部分已经讲过了，这里对具体的例子进行计算

通过检查 Ricker 波形的频谱，很明显存在更高的频率，即在距中心频率 -40dB 的电平上，存在 2-3 倍的频率。在这种情况下，模型中存在的最高有效频率可能在 4 GHz 左右。半空间（速度最低）中 4 GHz 的波长为：

这将提供 31/10=3mm 的最小空间分辨率。但是，圆柱体的直径为 20 mm，因此可以解析为 7 个单元。因此，更好的选择是 2 mm，它将钢筋的直径解析为 10 个单元。

域的大小应该包含足够的体积，加上10（默认值）个单元格的 PML （完全吸收边界条件）边界

PML 与任何对象间大约要相隔10个单元网格，域一般要足够大，才能满足建模条件。对于 2D 模型，在z方向上也应该指定一个单元格。

希望看到圆柱体的反射，因此时间窗口必须足够长，以允许电磁波从源头通过半空间传播到圆柱体并被反射回接收器,

其中0.18 = 0.09*2。

这是所需的最短时间，但源波形的宽度为 1.2 ns，为了允许整个源波形反射回接收器，将测试 3 ns 的初始时间窗口。

#time_window: 3e-9 //单位是秒

.out文件

gprMax 会生成一个输出文件，该文件与输入文件同名，但附加了 URL。输出文件使用广泛支持的 HDF5 格式，该格式旨在存储和组织大量数值数据。有许多免费工具可用于读取 HDF5 文件，比如，MATLAB 。

参考资料：