GPR建模

Georadar（地下探测的无线电技术设备，GPR，探地雷达）目前使用非常广泛-从绘制兔子洞和研究蜥蜴到寻找地雷 ，仍然是一种相当昂贵的娱乐。


GPR显示屏（来自英国电视节目“时间的指挥”的帧）

但是可以在不购买/租用“铁”设备的情况下评估其能力并研究地雷达的电磁场与环境相互作用的各个方面。 在GNU GPL v3许可下分发的gprMax软件包（ gpr -GPR的缩写， Max -James Clerk Maxwell的名字的首字母，奠定了电动力学的基础）已经可以为我们提供帮助。
该项目的作者始于1996年，是诺森比亚大学的Craig Warren和爱丁堡大学的Antonis Giannopoulos。 该软件包最初是用C开发的，然后完全用Python 3 / Cython组合重写。


安装该软件包需要已安装的编译器，这些编译器支持OpenMP（Microsoft Visual C ++ 2015 Build Tools（建议使用此版本！）对于Windows / gcc for Linux），NumPy库和Cython编译器。 从GitHub上的存储库下载并解压缩项目的源代码后，转到根文件夹并执行以下命令：

python setup.py build python setup.py install 

作为“快速开始”，我们考虑使用一个简单的二维示例来处理包装：脉冲雷达的发射天线T（脉冲GPR）发射电磁脉冲，其一部分能量以直接波（DW-直接波）的形式直接到达接收天线R，并且有些穿透通过沙，它从导电圆柱体的表面反射，并以反射波（RW-反射波）的形式到达接收天线：


输入文件格式
在项目的根文件夹中创建models文件夹，在其中放置文本文件hello.in，其中包含用于执行模拟的命令（以下给出的命令与项目的当前（第三）版本相对应）。

每个团队具有以下形式：

 #: _1 _2 _3 ... 

一行上只能写入一个命令，并且包含该命令的行的第一个字符必须为＃。

命令可以带有注释：

 ## 

命令的顺序对于用于构造对象的命令很重要-这些命令以它们在输入文件中出现的顺序执行。

脉冲形状

地质雷达发射的电磁脉冲持续几分之一纳秒，此外，最常用的是三种脉冲形状：

1. 一个正弦波周期（正弦）
2. 高斯动量（高斯）
3. 墨西哥帽（ricker）与高斯函数的二阶导数成比例，这种脉冲的曲线形状类似于草帽（这种脉冲形式在1953年被美国地球物理学家诺曼·里克尔（Norman Ricker）用于研究地震信号）
   

在我们的示例中，我们选择一个具有中心频率的高斯脉冲（脉冲类型-高斯） $$$f_c = 1 \\ GHz = 1 \ cdot 10 ^ 9 Hz$ 由命令定义：

 #waveform: gaussian 1 1e9 pulse 

（1-条件脉冲幅度，脉冲-脉冲标签）

在这种情况下，模拟中使用的动量由以下表达式描述：

$$

$$W（t）= e ^ {-2 \ cdot {\ pi} ^ 2 \ cdot {f_c} ^ 2 {（t- {1 \ over {f_c}}）^ 2}$$

环境模型和坐标系

在2D建模中，将研究区域划分为给定大小的像元，并且模型坐标系如下所示-X和Y轴形成计算平面（宽度为 $$$w$ 又高 $$$h$ ），则模型沿Z轴的长度的值等于采样步长 $$$\三角洲$ 。


选择采样步骤时，可以遵循经验法则-步长不应超过模型中研究的最小波长的十分之一（“每个波长10个像元”）：

$$

$$\ Delta \ le 0,1 \ cdot {\ lambda} _ {min}$$

为了确定波长，需要知道发射信号频谱中所考虑的最大频率以及所考虑介质中的波速。

电磁波在具有相对介电常数的介质中的传播速度 $$${\ epsilon} _r$ ，以厘米/纳秒为单位-雷达采用的速度单位由以下表达式确定：

$$

$$v = {30 \ over {\ sqrt {{\ epsilon} _r}}}$$

以厘米为单位的波长由以下表达式确定：

$$

$$\ lambda = {v \ over f}$$

（ $$$f$ -以GHz为单位的频率）。
要显示高斯脉冲的形状及其频谱，可以使用以下命令：

 python -m tools.plot_source_wave gaussian 1 1e9 5e-9 1e-12 -fft 

（高斯-脉冲类型，1-脉冲幅度，1e9-中心频率（1 GHz），5e-9-脉冲显示持续时间（5 ns），1e-12-时间步长（1 ps））


根据高斯脉冲的频谱 $$$f_s = 1 \ GHz$ 确定频率为-40 dB $$$f \大约3 \ GHz$ 。
在此示例中，我们是探测对象所在的介质，我们选择具有相对介电常数的干砂 $$${\ epsilon} _r = 3$ 。
求出电磁波在沙子中的传播速度：

$$

$$v = {30 \ over {\ sqrt {3}}} = 17.3 \ cm / ns$$

定义沙子中的波长：

$$

$$\ lambda = {v \ over f} = {17.3 \ over {3}} = 5.8厘米= 58毫米$$

基于此，我们选择所有轴都相同的步骤（ $$$\ Delta = {\ Delta} _x = {\ Delta} _y = {\ Delta} _z$ ），为了方便起见，等于2 mm = 0.002 m（整数步长适合1 cm）：

 #dx_dy_dz: 0.002 0.002 0.002 

将模拟区域限制为宽度为矩形的区域 $$$w$ 等于80厘米= 0.8 m和高度 $$$h$ 等于60厘米= 0.6 m：

 #domain: 0.60 0.60 0.002 

（对于二维模型，要求表示等于一个台阶的厚度（0.002））

仿真区域的大小和空间步长大小决定了模型单元的数量，并因此决定了计算机内存的需求。

我们用比电导率来描述沙子 $$$\ sigma = 0.01 \\ cm / m$ 和相对介电常数 $$${\ epsilon} _r = 3$ 命令：

 #material: 3 0.01 1 0 sand 

（1对应于相对磁导率 $$${\ mu} _r$ 等于1（无磁性），0（无磁损耗）和沙子 （此材料的标签）。

用大部分模拟区域填充沙子（从y = 0到y = 38 cm = 0.38 m）：

 #box: 0 0 0 0.80 0.38 0.002 sand 

（0 0 0-左下角的坐标，0.80 0.38 0.002-右上角的坐标（0.002-采样步长））
其余的是自由空间（标签为free_space），其空气属性几乎相等（ $$${\ epsilon} _r = 1$ ， $$${\ mu} _r = 1$ ， $$${\ sigma} = 0$ ）
模拟区域的边界表示为吸收边界条件（ABC）。

作为目标，我们从理想导体（完全反射电磁辐射）中选择一个圆柱体，其半径为6 cm = 0.06 m，中心位于坐标x = 25 cm = 0.25 m和y = 10 cm = 0.1 m的点上：

 #cylinder: 0.250 0.1 0 0.250 0.1 0.002 0.06 pec 

（pec是一种完美的导电材料）

天线类

模拟的GPR配备了两个天线-发送和接收。
在我们的案例研究中，我们代表一个带有赫兹偶极子的发射天线，该天线的长度等于采样步长（在三维情况下，您可以从一个广泛的库中选择一个天线），它位于区域中间左侧5厘米（x = 35）的距离内（与被测介质接触） cm = 0.35 m，y = 38 cm = 0.38 m）：

 #hertzian_dipole: z 0.35 0.380 0.0 pulse 

（z是偶极极化轴（对于二维情况（2D TMz模式），只有z有效），pulse是天线辐射的脉冲形状的标记）

接收天线通常位于距离接收器小的恒定距离处，这被称为天线单元的底部（天线相对位置的此选项称为“公共偏移”）。 作为基础，选择10 cm的距离，因此水平坐标为35 + 10 = 45 cm = 0.45 m（区域中间右边5 cm）：

 #rx: 0.45 0.38 0.0 

模拟间隔

确定用于建模的时间窗口的选择，以便从目标反射的信号有时间到达接收天线。

我们根据天线到目标的距离，确定在考虑情况下信号所需的大概时间 $$$h = 18 \\ cm$ ：
$$$t \大约{2 \ cdot h \ over {v}} = {2 \ cdot 18 \ over 17.3} = 2.1 \ \ ns$
假设中心频率为1 GHz的雷达高斯脉冲的顶部相对于时间轴的起点偏移1 ns，我们选择5纳秒的时间窗：

 #time_window: 5e-9 

造型

因此，输入文件的内容如下：


我们开始建模过程：

 python -m gprMax models\hello.in 

为了进行仿真，使用了时域有限差分法（FDTD，时域有限差分法）（该方法的基本算法由Kane Yee提出），之后将模型划分为的单元称为Yee单元。 。 通过求解每个单元的麦克斯韦方程，可以在时域中获得数值解。

在二维情况下（2D TMz模式），仅计算分量 $$$E_z$ 电场和分量 $$$H_x$ 和 $$$H_y$ 磁场。

如果超出了可用内存量，则会发出一条有关内存不足的消息以执行模拟：


如果任何对象不在模拟范围内，则会显示错误消息：


要使用所描述的模型执行仿真，它仅占用了约56 MB的RAM（如果将步长减半-最多1毫米-内存需求将增加到99 MB）。

模拟完成后， hello.out文件将显示在models文件夹中，其中包含HDF5格式的模拟结果，该结果是专门为存储数字数据而设计的。

轨道建设

为了使结果可视化，我们构造迹线：

 python -m tools.plot_Ascan models\hello.out 

在打开的窗口中显示的每个轨迹（A扫描）都会在接收天线的位置显示电磁场成分之一的时间图：


在路径上可以看到直接来自发射天线（DW）的直接波和从目标（RW）反射的波。

时间的水平轴与通过物质中电磁波的速度反射信号的目标深度有关。

但是，如果将柱面命令放在输入文件中box命令的前面，会发生什么？


从圆柱体反射的信号消失了-沙子吸收了圆柱体（这是构造物体顺序的重要性的示例）。

资料建立

但是，更有意义的是雷达图（雷达图）-轮廓图（B扫描），它是沿着给定方向移动地质雷达时建立的许多轨迹的组合-这是沿着研究区域移动带有雷达的小车的过程。

通过将天线移动到水平轴的开头来更改其描述：

 #hertzian_dipole: z 0.10 0.38 0.0 pulse #rx: 0.20 0.38 0.0 

我们将移动天线的步骤设置为等于1 cm = 0.01 m：

 #src_steps: 0.01 0 0 #rx_steps: 0.01 0 0 

因此，输入文件的内容如下：


以批处理模式运行模拟：

 python -m gprMax models\hello.in -n 50 

（50是雷达移动的步数）。

启动后，将依次对50个GPR位置执行仿真：


模拟完成后，模型文件夹中有50个文件hello1.out ... hello50.out。
使用以下命令将这些文件合并到hello_merged.out文件中：

 python -m tools.outputfiles_merge models/hello 

建立个人资料：

 python -m tools.plot_Bscan models\hello_merged.out Ez 

（Ez是建立轮廓的电磁场的分量，该分量直接转换为电压）


如您所见，由直射波引起的水平条显示在上方的轮廓上，而由反射波引起的特征双曲线并显示了目标到GPR移动时距离的变化。

右侧的图例显示了匹配的颜色和场强。 $$$E_z$ （显示在轨道上）：
红色- $$$E_z> 0$
白色- $$$E_z = 0$
蓝色- $$$E_z <0$
通过分析这样的轮廓，我们可以得出关于目标的深度，大小甚至形状的结论，并且神经网络也用于此。


GPR显示屏上的路线和配置文件（来自英国电视节目“时间的指挥”的帧）

但是反射不仅会发生在导电物体上，还会发生在介电常数不同的两层边界处。

在模型的下部创建第二个具有介电常数的沙层 $$${\ epsilon} _r = 9$ ：



可以看出，在直接波（DW）的“迹线”下方，出现了从两个砂层界面反射的波（RW）的线性段。

在此示例的范围之外，gprMax 功能包括三维建模，复杂的表面形状，详细的天线模型，电磁波的散射等。此外，该软件包不仅可以用于模拟GPR，还可以用于研究电磁波在各种环境中的传播，以及使用NVIDIA CUDA技术可以加快计算速度，与使用OpenMP的CPU上的并行计算相比，该技术的速度提高了十倍。 为了增加建模的灵活性，您可以在输入文件中放置Python代码块。

使用gprMax软件包的一些示例：

• 性质对GPR杀伤人员地雷搜寻效果的影响评估
• 检查使用探地雷达进行玉米种子播种深度无损确定的可能性
• 在研究公路地下通道结构强度方面提高GPR现场数据的解释质量
• 反射信号的振幅与钢筋直径之间关系的研究
• 雷达检测树干中空隙的适用性的评估

有用的链接：

gprMax官方网站
GprMax用户指南
YouTube上的gprMax