GPRMAX2D知识点汇总.doc

GprMax是爱丁堡大学的Antonis Giannopoulos于1996年推出来的一种基于时域有限差分（FDTD）算法和理想匹配层（PML）边界吸收条件的探地雷达正演数值模拟软件，用于探地雷达成像研究。其中，GprMax2D是二维正演，GprMax3D为三维正演。
　　该软件可以在Windows、Linux和MacOS三个平台上使用。本文主要针对Windows平台进行说明。
一、软件获得

　　该软件为免费软件，可以去GprMax官网下载。也可点此直接下载。

二、准备工作

　　软件无需安装，下载后用解压工具解压，找到Windows文件夹，直接双击GprMax2D.exe即可运行。
　　但是为了方便批量模拟，建议把Windows文件夹下的文件（cygwin1.dll、GprMax2D.exe和GprMax3D.exe）全部复制到系统盘系统搜索路径下，例如：C:\Windows。如果只进行二维正演，只复制cygwin1.dll和GprMax2D.exe即可。
　　当然也可放在任意路径下，只要使用时包含所在路径就OK了。笔者习惯放在系统要目录下，即C:\。
　　另外，最好建立一个输入文件和模拟结果存放的专用文件夹，并且把tools文件夹下的文件（gprmax.m、gprmax2g.m、gprmax3g.m、gprmaxde.m和gprmaxso.m）全部复制到该专用文件夹下，注意文件夹名最好使用英文。如果只进行二维正演，只复制gprmax.m和gprmax2g.m即可。笔者习惯使用D:\GPR。

　　正演结果需要用MATLAB进行绘图，因此需要安装有MATLAB软件。软件下载地址和安装方法此处不再说明，以后使用将假设读者已经成功安装MATLAB软件。

三、软件使用

　　一般进行数值模拟时通常都是若干个对比模型进行模拟，因此本文只介绍批量模拟的方法。如此一来，笔者将认为读者的软件和笔者一样放在C:\，输入文件和gprmax.m、gprmax2g.m文件已经放在D:\GPR。
　　欲进行正演模拟需要先建立输入文件，输入文件的建立将在后面介绍，此处假设在D:\GPR已经有若干输入文件，文件名分别为fname1.in，fname2.in，fname3.in。
　　打开记事本，写入以下两行直线间的内容，另存为*.bat文件，例如：GPR.bat。其中括号里的内容为本条语句的说明。
------------------------------------------------------------------------------------------
path C:\;%path% (把软件所在路径包含进系统搜索路径)
cd/d D:\GPR (设置专用文件夹路径为活动路径)
gprmax2d D:\GPR\fname1.in （按顺序进行批量正演模拟）
gprmax2d D:\GPR\fname2.in
gprmax2d D:\GPR\fname3.in
------------------------------------------------------------------------------------------
　　完成后保存并关闭，然后双击运行*.bat文件即可进行批量正演模拟

四、建立输入文件

　　输入文件是纯文本文件，可以用任何文字处理程序编辑，而且内容全部为英文。打开记事本，写入输入文件内容，另存为*.in即可。

　　输入文件必需包括模拟模型的所有必要信息，每一个有效命令行都必须以符号（#）开始，否则该命令行无效。所有命令的一般语法为：
　　　#命令名: 参数1 参数2 参数3 ……

　　为了方便介绍GprMax2D的命令及参数，做了以下约定：
　　　·f 代表实数，如1.5
　　　·i 代表整数，如15
　　　·c 代表单个字母，如y
　　　·str 代表字条串，如air
　　　·file 代表文件名，如test.in
　　　·所有的长度单位均为m
　　　·所有的时间单位均为s
　　　·所有的频率单位均为Hz
　　　·模型空间的坐标原点（0,0）在左下角
　　输入文件中除了#analysis:和#end_analysis:命令外，先后顺序一般是没有要求的。

　　GprMax2D 2.0版一共有32个命令：
------------------------------------------------------------------------------------------
#title: str
输出文件的标题，包含在输出文件中。

#domain: f1 f2
指定模型空间大小，单位是米。
·f1 x方向的大小
·f2 y方向的大小

#dx_dy: f1 f2
指定网格步长，单位是米。
·f1 x方向的步长，即Δx
·f2 y方向的步长，即Δy

模型划分的网格数目是由#domain:和#dx_dy:两个命令共同决定的。网格划分越小模型越精确，但模拟的计算量也越在，需要权衡确定。

#time_step_stability_factor: f1
时间步长的稳定系数，取值范围（0，1］，实际使用时间为f1×Δt。

#time_window: f1 or i1
时窗大小，单位是秒。
·f1 为时窗，如f1=16e-9，时窗为16ns
·i1 为时间步长数，时窗为i1×Δt。

#messages: c1
是否开启屏幕信息。y(es)开启n(o)关闭，默认n。

#number_of_media: i1
介质数目，i1应大于10。系统分配了10种介质（自由空间即空气free_space，金属pec和8种用户自定义介质）的内存空间，当模型中使用介质数量超过10种时需要使用该命令。

#nips_number: i1
指定分配的空间。存储重要信息模型的数组大小可以在内部计算，但通常比实际需要的大，可以使用较小的数字节省空间。如果分配过小，系统会发出错误并提醒分配更多空间

#media_file: file1
要包含的介质文件的文件名（如果有必要还包括文件路径）。介质文件格式为：
　　　# f1 f2 f3 f4 f5 f6 str1
参数的具体含义参考#medium:命令。

#geometry_file: file1
存储模型文件（*.geo）。

#medium: f1 f2 f3 f4 f5 f6 str1
在输入文件中定义介质属性。
·f1 静态时介质的相对介电常数
·f2 在理论上的无限频率的相对介电常数
·f3 介质的弛豫时间（秒）
·f4 静态时介质的电导率（西门子/米）
·f5 相对磁导率的介质
·f6 介质的磁导率
·str1 介质的名称
如果不想使用debye介质，可设置f3为0.0，此时系统将只使用f1和f4来描述介质介电性能。
如果介质为非磁性，则设置f5为1.0，f6为0.0。
如果使用f3的值，该值应始终高于模型中使用的时间步长Δt。
------------------------------------------------------------------------------------------
该组命令不建议使用，除非是熟悉ABC的人用来定制和优化吸收边界条件。

#abc_type: pml
恢复默认的PML边界。

#abc_order: i1
ABC次序。取值范围1，2，3，默认为3。次序越低，性能越差。

#abc_stability_factors: f1 f2 f3
指定ABC的稳定系数。

#abc_optimization_angles: f1 f2 f3
优化ABC为特定的入射角。

#abc_mixing_parameters: f1 f2 f3 f4 f5 f6
改变使用ABC的时间和偏导数的离散化。默认使用box。

#pml_layers: i1
指定PML使用的网格数量。默认为8。数量越多性能越佳。
------------------------------------------------------------------------------------------
#box: f1 f2 f3 f4 str1
使用定义好的介质定义一个矩形模型。
·f1 矩形左下角x坐标，单位米
·f2 矩形左下角y坐标，单位米
·f3 矩形右上角x坐标，单位米
·f4 矩形右上角y坐标，单位米
·str1 使用的介质的名称

#cylinder: f1 f2 f3 str1
使用定义好的介质定义一个圆形模型。
·f1 圆心x坐标，单位米
·f2 圆心y坐标，单位米
·f3 圆半径R，单位米
·str1 使用的介质的名称

#x_segment: f1 f2 f3 f4 f5 str1
使用定义好的介质定义一个圆沿x轴裁剪的模型。
·f1 圆心x坐标，单位米
·f2 圆心y坐标，单位米
·f3 沿x轴开始的位置
·f4 沿x轴结束的位置
·f5 圆半径R，单位米
·str1 使用的介质的名称

#y_segment: f1 f2 f3 f4 f5 str1
使用定义好的介质定义一个圆沿y轴裁剪的模型。
·f1 圆心x坐标，单位米
·f2 圆心y坐标，单位米
·f3 沿y轴开始的位置
·f4 沿y轴结束的位置
·f5 圆半径R，单位米
·str1 使用的介质的名称

#triangle: f1 f2 f3 f4 f5 f6 str1
·f1 第一顶点x坐标，单位米
·f2 第一顶点y坐标，单位米
·f3 第二顶点x坐标，单位米
·f4 第二顶点y坐标，单位米
·f5 第三顶点x坐标，单位米
·f6 第三顶点y坐标，单位米
·str1 使用的介质的名称

模型的建立是有先后顺序的，后建立的模型会把先建立的模型覆盖住。
------------------------------------------------------------------------------------------
该组命令是数值模拟的核心。一个输入文件中可以出现多组数值模拟命令，但是数值模拟的控制命令必须放在#analysis:和#end_analysis:之间。

#analysis: i1 file1 c1
数值模拟开始的命令。
·i1 总步数
·file1 模拟结果输出文件名（*.out）
·c1 类型，取值范围a（文本文件）或b（二进制文件）

#end_analysis:
数值模拟结束的命令，无参数。

#tx: f1 f2 str1 f3 f4
发射天线的属性。
·f1 发射天线位置的x坐标，单位米
·f2 发射天线位置的y坐标，单位米
·str1 定义过的发射天线的名称
·f3 发射天线的延迟时间
·f4 发射天线的清除时间

#rx: f1 f2
接收天线的属性。
·f1 接收天线位置的x坐标，单位米
·f2 接收天线位置的y坐标，单位米

#rx_box: f1 f2 f3 f4 f5 f6
接收区域的属性。
·f1 接收区域左下角的x坐标，单位米
·f2 接收区域左下角的y坐标，单位米
·f3 接收区域右上角的x坐标，单位米
·f4 接收区域右上角的y坐标，单位米
·f5 接收区域在x轴的步长，即Δx，单位米
·f6 接收区域在y轴的步长，即Δy，单位米

#snapshot: i1 f1 f2 f3 f4 f5 f6 f7 file1 c1
或
#snapshot: i1 f1 f2 f3 f4 f5 f6 i2 file1 c1
获取模型一个区域内，某一时刻的电磁场的快照。
·i1 第几道的扫描，值介于1和#analysis:命令中的总步数。
·f1 区域左下角的x坐标，单位米
·f2 区域左下角的y坐标，单位米
·f3 区域右上角的x坐标，单位米
·f4 区域右上角的y坐标，单位米
·f5 区域在x轴方向上的采样间隔，单位米
·f6 区域在y轴方向上的采样间隔，单位米
·f7 快照的时刻，单位秒
·i2 快照的迭代次数
·file1 快照存储的文件名（*.out）
·c1 类型，取值范围a（文本文件）或b（二进制文件）

#tx_steps: f1 f2
发射天线的移动步长。
·f1 发射天线沿x轴的移动步长，单位米
·f2 发射天线沿y轴的移动步长，单位米

#rx_steps: f1 f2
接收天线的移动步长。
·f1 接收天线沿x轴的移动步长，单位米
·f2 接收天线沿y轴的移动步长，单位米
------------------------------------------------------------------------------------------
#line_source: f1 f2 str1 str2
·f1 发射天线的振幅
·f2 发射天线的频率
·str1 信号波形，有五种：
·cont_sine 连续正弦波
·sine 脉冲正弦波
·gaussian 高斯波
·ricker 雷克子波
·user 用户自定义的激励函数
·str2 发射天线的名称

#excitation_file: str1
选择用户自定义的激励文件。自定义激励文件是文本文件。
------------------------------------------------------------------------------------------

五、MATLAB绘图

　　笔者使用MATLAB7.0，以该版本为例进行说明，其他版本可参考使用。为方便说明，笔者假设读者已对MATLAB具有最基本的了解和使用能力，同时主要以命令来完成操作。读者若对MATLAB一无所知，请另行学习。

1.设定当前目录和搜索路径。
------------------------------------------------------------------------------------------
cd D:\GPR %设定当前目录
addpath D:\GPR %设定本次搜索路径
------------------------------------------------------------------------------------------

2.模型文件绘图
------------------------------------------------------------------------------------------
geo = 'fname1.geo';
[mesh,header,media] = gprmax2g(geo);
[M,N] = size(mesh); %获取模型大小
imagesc((1:N)*header.dx,(1:M)*header.dy,mesh) %绘图
xlabel('Distance/m'); %设置x轴标签
ylabel('Depth/m'); %设置y轴标签
------------------------------------------------------------------------------------------

3.输出文件绘图
------------------------------------------------------------------------------------------
out = 'fname1.out';
[Header,Fields] = gprmax(out);
NN = 1:Header.NSteps; %移动次数
Position = Header.dx*Header.tx+(NN-1)*Header.dx*Header.TxStepX); %天线每次所在位置
Data(:,:) = Fields.ez(:,1,:); %转换数组格式
imagesc(Position,Fields.t*1e9,Data) %绘图
colorbar; %显示色彩条
xlabel('Distance/m'); %设置x轴标签
ylabel('Time/ns'); %设置y轴标签
------------------------------------------------------------------------------------------
　　绘出的图像请及时保存。
例子：
#medium: 6.0 0.0 0.0 0.01 1.0 0.0 concrete
#medium: 20.0 0.0 0.0 0.1 1.0 0.0 wet_sand
---------------------------------------------------------------------
#domain: 2.5 0.5
#dx_dy: 0.0025 0.0025
#time_window: 12e-9
---------------------------------------------------------------------
#box: 0.0 0.0 2.5 0.45 wet_sand
----- Get a cylinder of free space and then put a slab of concrete
to cut it in half -------------------------------------
#cylinder: 1.05 0.3 0.2 free_space
#box: 0.0 0.3 2.5 0.45 concrete
#cylinder: 0.25 0.375 0.0125 pec
#cylinder: 0.45 0.375 0.0125 pec
#cylinder: 0.65 0.375 0.0125 pec
#cylinder: 0.85 0.375 0.0125 pec
#cylinder: 1.05 0.375 0.0125 pec
#cylinder: 1.25 0.375 0.0125 pec
#cylinder: 1.45 0.375 0.0125 pec
#cylinder: 1.65 0.375 0.0125 pec
#cylinder: 1.85 0.375 0.0125 pec
#cylinder: 2.05 0.375 0.0125 pec
#cylinder: 2.25 0.375 0.0125 pec
#triangle: 0.25 0.3 0.65 0.3 0.45 0.1 free_space
#triangle: 1.45 0.3 1.85 0.3 1.65 0.1 free_space
#box: 2.05 0.1 2.25 0.3 free_space
#cylinder: 2.15125 0.10125 0.1 wet_sand
-------------------------------------------------------------
#line_source: 1.0 900e6 ricker MyLineSource
-------------------------------------------------------------
#analysis: 115 bre2.out b
#tx: 0.0875 0.4525 MyLineSource 0.0 12e-9
#rx: 0.1125 0.4525
#tx_steps: 0.02 0.0
#rx_steps: 0.02 0.0
#end_analysis:
------------------------------------------------------------
#geometry_file: bre2.geo
#title: BRE Model 2
#messages: y