GprMax 3.1.5 建模的in文件编写详解（1）

gprMax是一款优秀的基于时域有限差分方法（FDTD）的电磁波数值模拟软件，目前为止，它还没有图形用户界面（GUI），它的建模关键在于in文件的编写。基于此，在浏览gprMax官方网站link后，加上自己的理解，这篇帖子希望对gprMax爱好者有所帮助。

文章目录

GprMax 3.1.5 建模的in文件编写详解（1）
- 一、一个简单的in文件范例
- 二、逐行解释上面命令的意思
- - 1、标题
  - 2、建模域
  - 3、网格大小
  - 4、时窗
  - 5、模型属性
  - 6、子波
  - 7、源
  - 8、接收点
  - 9、模型
  - 10、显示模型
  - 11、波场快照
- 三、gprMax运行结果

一、一个简单的in文件范例

这个例子相当于gprMax中的“hello,world!"，基本包含了in文件写入的要素，它说明了如何建立介电半空间中金属圆柱模型。
具体的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#waveform: ricker 1 1.5e9 my_ricker
#hertzian_dipole: z 0.100 0.170 0 my_ricker
#rx: 0.140 0.170 0#box: 0 0 0 0.240 0.170 0.002 half_space
#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

二、逐行解释上面命令的意思

1、标题

标题的书写格式如下：

#title:A-scan from a metal  cylinder buried in a dielectric half-space

标题的标识符为：#title:，在冒号后面可以定义你自己建模的名称，为了易于辨识，最好命名为材料加模型形态，且用英文。

2、建模域

域的书写格式如下：

#domain:0.240 0.210 0.002

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

域的大小应该包含足够的体积，加上10个单元格的 PML （完全吸收边界条件）边界（默认值是10个单元），PML 与任何对象间大约要相隔10个单元网格，域一般要足够大，才能满足建模条件。对于 2D 模型，在z方向上也应该指定一个单元格。

3、网格大小

网格大小的书写格式如下：

#dx_dy_dz:0.002 0.002 0.002

其中，0.002 0.002 0.002分别为x方向、y方向、z方向的网格单元大小。至于网格单元大小如何确定，有一定的规律，经验法则要求探地雷达的空间分辨率应为模型中最小波长的十分之一。为了确定最小波长，需要知道模型中存在的最高频率与最低速度，最高频率不是雷克子波的中心频率，而是检查雷克子波的完整频谱，找到存在的最高频率。假设水平-40dB的中心频率，最高频率是它的2~3倍，在这种情况下，模型中可能出现的最高频率可能在 4 GHz 左右。半空间（速度最低）4 GHz 的波长是：
λ=3×108÷(4×109×60.5)≈31mm\lambda = 3\times{10^8}\div(4\times{10^9}\times6^{0.5} )\approx31 mm λ=3×108÷(4×109×60.5)≈31mm
空间分辨率大约为31÷10≈3mm31\div10\approx3mm31÷10≈3mm
然而，我们设置的圆柱的直径是20 mm ，20÷3≈720\div3\approx720÷3≈7
这将使圆柱剖分成7个单元，因此，更好的选择是 2 mm，这将使钢筋的直径为10个单元。

4、时窗

时窗的书写格式如下：

#time_window:3e-9

如果我们想要看到圆柱反射回来的电磁波，时间窗口必须合适，使电磁波从源发射，通过半空间传播到圆柱，然后反射回接收器，即。
t=0.180÷((3×(108)×(6−0.5))≈1.5nst = 0.180\div((3\times(10^{8})\times(6^{-0.5}))\approx1.5 nst=0.180÷((3×(108)×(6−0.5))≈1.5ns
这是所需的最小时间，由于源波形的脉宽为 1.2 ns，为了让整个源波形反射回接收机，我们设置时间窗口为 3 ns （时间窗口宁可设大，不可设小）。

5、模型属性

#material:6 0 1 0 half_space

我们的示例中有三种不同的材料：空气、介电半空间、金属圆柱体，空气已经作为内置材料存在于 gprMax 中，可以使用 free_space 标识符来访问它。我们假设金属圆柱体是一个完美的电导体，它也是gprMax的内置材料，可以使用pec标识符访问。所以唯一需要定义的材料就是介电半空间，设它是一种非磁性材料，相对介电常数为6，电导率为0，相对磁化率为1，磁损为0，使用标识符 half_space。

6、子波

使用#waveform:命令指定源发射的波形，这里选择Ricker波形，指定振幅为1，中心频率为1.5 GHz，my_ricker为任意标识符。gprMax中还有其他的子波，如gaussian、gaussiandot、gaussiandotnorm、gaussiandotdotnorm、gaussianprime、gaussiandoubleprime、sine、contsine等。

#waveform:ricker 1 1.5e9 my_ricker

7、源

#hertzian_dipole:z 0.100 0.170 0 my_ricker

gprMax所用的源在建模之前就需要确定，通常基于GPR系统或正在建模的场景。低频用于穿透深度大的情况，高频用于较低的穿透深度和更高的分辨率。这里我们用#hertzan_dipole:命令创建赫兹磁偶极源，指定z方向激化（如果创建b扫描，则垂直于测量方向），定位在平板表面，并使用已经创建的Ricker波形。

8、接收点

#rx:0.140 0.170 0

利用#rx:命令创建接收点的位置，0.140、0.170、0分别是接收点的x、y、z坐标。在接收点的位置，模型迭代次数的电场和磁场分量值将被保存到文件中，其后还可以加两个参数str1和str2，str1是接收方的标识符，str2是这个接收器的输出列表，它可以从Ex、Ey、Ez、Hx、Hy、Hz、Ix、Iy、Iz中任意选择。如果没有设置这两个参数时，默认输出电场和磁场的六个分量。

在gprMax 3.1.5 中，源或接收点是可以移动的，使用命令#rx_steps:f1、f2、f3，f1、f2、f3是移动接收点的增量（x,y,x)，这个命令不可以移动具有复杂几何形状的模型。

9、模型

#box:0 0 0 0.240 0.170 0.002 half_space
#cylinder:0.120 0.080 0 0.120 0.080 0.002 0.010 pec

我们用#box:命令来设置背景介质，#box:是引入正交六面体的命令，前三个参数为六面体的左下顶点的坐标，之后三个参数是六面体的右下坐标，三维坐标系遵守笛卡尔坐标系和右手定则。half_space是标识符，与我们前面定义的材料标识符对应。用#cylinder:命令引入圆柱体，圆柱体的轴线方向是任意的（非正交），前三个参数是圆柱体一个面的中心坐标（x,y,z），接下来的三个参数是圆柱体另一个面的中心坐标（x,y,z），第7个参数是圆柱体的半径，pec是材质标识符，是gprMax的内置材料，金属圆柱体被建模成一个完美的电导体。

10、显示模型

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

利用#geometry_view:命令将我们建立模型的几何信息输出到文件中，这些文件使用开源可视化工具包（.VTK）格式，可以在许多免费阅读器（如Paraview）中查看。这个命令可以用来创建模型的几个3D视图，这对于检查模型是否已经按照预期构建是很有用的。该命令的第1、2、3个参数是坐下（x,y,z)坐标，单位为米；第4、5、6个参数是有右上（x,y,z）坐标；第7、8、9个参数是网格大小，通常情况下，这与模型中网格大小是相同的；第10个参数是文件的文件名，这个文件保存路径与in文件的路径相同；第11个参数可以是n或f，n是“正常”，f是“细”，指定输出文件的大小和图形的好坏。当你只想检查模型的几何模型的形状时，可以使用可选命令参数-geometry-only运行gprMax，这样不会运行正演模拟。

11、波场快照

#snapshot:0 0 0 0.240 0.170 0.002 0.002 0.002 0.002 1e-9 snap1

这个命令用于生成波场快照文件（VTK格式），波场快照就是在指定时刻的电磁场信息，用动图显示一段时间的电磁场信息，这个文件可以在许多免费软件（如Paraview）中查看。该命令的第1、2、3个参数是坐下（x,y,z)坐标，单位为米；第4、5、6个参数是有右上（x,y,z）坐标；第7、8、9个参数是网格大小，通常情况下，这与模型中网格大小是相同的；第10个参数是以秒为单位的时间（浮点数）或迭代次数（整数），表示将进行快照的时间点；第11个参数是存储快照文件的名称，快照文件会自动的储存在一个目录中，该目录的名称为输入文件名称，后面附有“-snaps”，对于多个模型运行，每个模型运行都有自己的目录。可以利用python脚本生成一系列快照。例如，要以0.1 ns的间隔创建30个从0.1 ns开始到3 ns的快照，可以在in文件中使用以下代码片段：

#python:
for i in range(1, 31):print(f."#snapshot:x1,y1,z1,x2,y2,z2,dx,dy,dz,{},snapshot{}".(i/10*1e-9, i))
#end_python:

三、gprMax运行结果

在Pycharm或VScode中运行gprMax的方法我已经在上篇帖子中讲过，在这里我们在cmd中运行命令。
打开cmd（gprMax需要添加到环境变量中），输入命令：

$ activate gprMax

进入到gprmax，然后再输入：

$ python -m gprMax user_models/cylinder_Ascan_2D.in

user_models/cylinder_Ascan_2D.in是指你的In文件路径，可以直接用鼠标将in文件拖到cmd中，按回车键运行。

运行完后再输入命令：

$ python -m tools.plot_Ascan user_models/cylinder_Ascan_2D.out

user_models/cylinder_Ascan_2D.out是指你的out文件路径，可以直接用鼠标将out文件拖到cmd中，按回车键运行。

我们运行的结果如下：

后续我还会发布一些高级操作，希望这篇入门级介绍对大家有所帮助，可以灌注我后面的帖子。