GprMax中文说明指导书.doc

GprMax中GprMax2D使用方法
V1.3
试验环境：
操作系统：Windows7
软件版本：MATLAB7.1&GprMaxV2.0

参考资料：
[1]UserGuideV2.pdf
一、GprMax2D软件使用
1.1 直接运行
..\GprMaxV2.0\Windows文件夹下GprMax2D.exe文件，会出现以下窗口(也能够在命令提醒符窗口输入命令运行) ：

1.2 输入文件名
注意：要输入文件全路径；*.in文件只要出现任何语法错误或路径错误，软件全部会自动关闭，不会有任何错误提醒。
出现以下画面(以自带例子文件bre1.in为例，*.in命令参考前面文章或[1])：

运行完成会发觉..\GprMaxV2.0\Windows文件夹下多了两个文件bre1.out、bre1.geo文件，复制到tools文件夹。

二、数据成像
tools文件夹下有五个m文件：gprmax.m，gprmax2g.m, gprmax3g.m, gprmaxde.m , gprmaxso.m。
这里只讲gprmax2g.m，gprmax.m这两个文件使用方法，其它三个以后有空再写。gprmax3g.m是处理GprMax3D几何数据；gprmaxde.m用来计算Debye公式(参考[1])介电常数；gprmaxso.m用于计算激励函数。
2.1 gprmax2g.m使用方法
gprmax2g函数用于读取GprMax2D软件仿真探地雷达模型生成二进制几何数据。
gprmax2g函数原型：
[mesh,header,media] = gprmax2g( 'filename' )
filename是.geo文件名；
media: 存放介质类型，media.type；
header: 存放模型几何参数；
header.title: 模型名称；
header.dx: 模型在X轴每次偏移大小(单位:m)；
header.dy: 模型在Y轴每次偏移大小(单位:m)；
header.dt: 最大许可时间步长(单位:秒)；
header.nx: 模型在X轴偏移次数；
header.ny: 模型在Y轴偏移次数；
例子：如输入文件*.in中定义：
#domain: 2.5 0.65
#dx_dy: 0.0025 0.0025
那么：header.dx=0.0025; header.dy=0.0025;
header.dt = 1/(c*sqrt(1/header.dx^2+1/header.dy^2)); (其中c=，为光速，公式参考[1]);
header.nx=2.5/0.0025=1000; header.ny=0.65/0.0025=260;
mesh: 存放模型数据，M x N数组，其中M为Y轴方向Yee单元数目，N为X轴方向Yee单元数目；M=header.nx，N=header.ny；
gprmax2g.m使用例子：
filegeo = 'bre1.geo';
[meshdata,header,media]=gprmax2g(filegeo);
figure(1);
[MM,NN]=size(meshdata);
imagesc((1:NN)*header.dx,(1:MM)*header.dy,meshdata)
axis('equal');
xlabel('x(m)');
ylabel('y(m)');

2.2 gprmax.m使用方法
gprmax函数用于读取GprMax2D和GprMax3D软件仿真探地雷达模型生成二进制波形数据。
gprmax函数原型：
[Header, Fields] = gprmax( 'filename' )
filename是.out文件名；
1) Header：该变量包含以下组员：
Header.title: 模型名称；
Header.iterations: 迭代次数；
Header.dx: 在X轴每次偏移大小；
Header.dy: 在Y轴每次偏移大小；
Header.dt: 最大许可时间步长；
Header.NSteps: 仿真次数；等于*.in文件中#analysis:命令第一个参数；
例子：如输入文件*.in中定义：
#domain: 2.5 0.65
#dx_dy: 0.0025 0.0025
#time_window: 12e-9
……
#analysis: 115 bre1.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:
……
那么，Header.iterations=ceil(Header.removed/Header.dt)；
Header.dx=0.0025；Header.dy=0.0025；
Header.dt= 1/(c*sqrt(1/Header.dx^2+1/Header.dy^2));
Header.NSteps=115；
Header.TxStepX=0.02/Header.dx=8; Header.TxStepY=0.0/Header.dy=0;
Header.RxStepX=0.02/Header.dx=8; Header.RxStepY=0.0/Header.dy=0;
Header.ntx:
Header.nrx:
Header.nrx_box:
Header.tx=0.0875/Header.dx=35；
Header.ty=0.4525/Header.dy=181；
Header.source=’MyLineSource’；
Header.delay=0；(等于#tx:命令第四个参数)
Header.removed=12e-9；(等于#tx:命令第五个参数)
Header.rx=0.1125/Header.dx=45；
Header.ry=0.4525/Header.dy=181；

2)Fields: 该变量包含以下组员：
Fields.t: 每个波形时间轴。数组大小Header.iterations*1；
Fields.ez: Z轴方向磁矢量数据。数组大小Header.iterations*1* Header.NSteps；
Fields.hx: X轴方向电矢量数据。数组大小Header.iterations*1* Header.NSteps；
Fields.hy: Y轴方向电矢量数据。数组大小Header.iterations*1* Header.NSteps；
因为GprMax仿真是基于FDTD算法，以上三者关系为：

具体可参考FDTD算法相关文件。

gprmax.m使用例子：
fileout = 'bre1.out';
[Header,Fields]=gprmax(fileout);
N=1:Header.NSteps; %移动次数
Position=Header.dx*Header.tx+(N-1)*(Header.dx*Header.TxStepX); %天线每次所在位置
Data(:,:)=Fields.ez(:,1,:); %转换数组格式
figure(2);
imagesc(Position,Fields.t*1e9,Data);%画图
colorbar
xlabel('Antena Position (m)');
ylabel('t(ns)');

GprMaxV2.0中GprMax2D输入文件命令
试验环境：
操作系统：Windows 7
软件版本：MATLAB 7.1 & GprMaxV2.0

参考文件：
GprMaxV2.0软件manual文件夹下UserGuideV2.pdf。
GprMaxV2.0下载地址：
说明：翻译得不好，还望大家见谅，因为我也是边看边译。
3.1 GprMax2D命令通常注意事项
为了描述GprMax2D/3D命令及其参数，我们作以下约定：
f表示浮点数(如1.5或15e-1、1.5e1)
i表示整数
c表示字符
str表示字符串
file表示输入文件名

全部空间距离基础单位为米
全部时间基础单位为秒
全部频率参数基础单位是Hz

3.2 GprMax2D 2.0版本共有32条命令：
#title:
#domain:
#dx_dy:
#time_step_stability_factor:
#time_window:
#messages:

#number_of_media:
#nips_number:
#media_file:
#geometry_file:
#medium:

#abc_type:
#abc_order:
#abc_stability_factors:
#abc_optimization_angles:
#abc_mixing_parameters:
#pml_layers:

#box:
#cylinder:
#x_segment:
#y_segment:
#triangle:

#analysis:
#end_analysis:
#tx:
#rx:
#rx_box:
#snapshot:
#tx_steps:
#rx_steps:

#line_source:
#excitation_file:

以下能够运行于GprMax中命令在GprMax2D版本2.0中将不再支持：
#scan:
#csg:
#extra_tx:

另外，部分命令旧参数规则亦发生了改变：
#tx:
#snapshot:
3.3 GprMax2D命令参数
为了愈加好地介绍GprMax2D命令，我们将它们分成四类：
通常命令：包含用于指定大小和模型离散
ABC相关命令: 许可定制和优化吸收边界条件
介质和对象结构命令: 用不一样参数来在模型中引入不一样介质和结构简单几何形状
激励和输出命令：用来放置源代码和模型输出点

运行GprMax2D最低程度命令以下：
•#domain:
•#dx_dy:
•#time_window:
•最少一个#analysis:及和其对应结束命令#end_analysis:
最少一个#tx:和#rx:，或#rx_box:命令
• 为了使#tx: 命令正确运行，同时需要一行新#line_source:命令

3.3.1通常命令
#title:
模型题目
#title: str
Str即是模型题目，必需是单行。

#domain:
模型范围(单位：米)
#domain: f1 f2
f1和f2分别代表x和y轴上量度大小
#domain: 1.0 1.5
表示1.0*1.5大小范围

#dx_dy:
表示x和y轴上偏移量(如△x、△y)
#dx_dy: f1 f2
表示x轴偏移f1和y轴方向上偏移f2

#domain: 1.0 1.5
#dx_dy: 0.1 0.1
则模型单元数为10*15

#domain: 1.0 1.5
#dx_dy: 0.01 0.01
则模型单元数为100*150

最大许可时间步△t和△x、△y约束关系以下：
(3.1)
C为光速，GprMax2D中计算△t使用3.1式等号。
#time step stability factor:
经过该命令，你能够修改GprMax2D所计算△t值，但必需满足3.1式要求。

#time_window:
用于指定所需总模拟时间。语法：
#time_window: f1
或
#time_window: i1

总迭代次数和模拟时间窗口： (3.2)

#number_of_media:
但你需要使用大于10个介质时必需使用该命令，因为GprMax2D只初始化了10个介质使用空间。
#number_of_media: i1
i1大于10。

#media_file:
使用该命令，你能够定义常见介质所描述结构参数文件路径和名称。
#geometry_file:
使用该命令，你能够定义模型几何信息二进制文件。这些信息能够用于创建模型图像和检验创建正确是否。
#geometry_file: model.geo
#messages:
使用该命令，你能够控制软件运行时在屏幕上输出信息。
#messages:c1
#nips_number:
该命令仅当在GprMax2D需要时才用于你输入文件中。
#nips_number:i1

4.1 GPRMAX3D命令通常注释
大多数能够用在GPRMAX3D命令几乎和GPRMAX2D是一样。不过，有部分命令是GPRMAX3D独有。为了简明扼要，这里只叙述和GPRMAX2D不一样命令。
基础空间和临时离散步骤分别为△t和△x，△y，△z。

4.2 命令清淡
在版本2里面，GPRMAX3D一共有42个命令用在3D GPR建模中。她们是：
#title：
#domain:
#dx_dy_dz:
#time_step_stability_factor:
#messages:

#number_of_media:
#nips_number:
#media_file:
#geometry_file:
#medium:

#abc_type:
#abc_order:
#abc_stability_factors:
#abc_optimization_angles:
#abc__mixing_parameters:
#pml_layers:

#box:
#cylinder:
#cylinder_new:
#cylindrical_segment:
#sphere:
#plate:
#edge:
#triangle:
#bowtie:
#thin_wire:

#analysis:
#end_analysis:
#tx:
#rx:
#rx_box
#snapshot:
#tx_steps:
#rx_steps:

#huygens_surface:
#hertzian_dipole:
#voltage_source:
#transmissio_line:
#plane_wave:
#excitation_file:

4.3 GPRMAX3D命令参数
为了愈加好地介绍GprMax2D命令，我们将它们分成四类：
1.通常命令：包含用于指定大小和模型离散
2. ABC相关命令: 许可定制和优化吸收边界条件
3.介质和对象结构命令: 用不一样参数来在模型中引入不一样介质和结构简单几何形状
4.激励和输出命令：用来放置源代码和模型输出点

4.3.1 通常命令
#title：
和GprMax2D使用方法形同
#domain：
用来指定模型尺寸（单位：米）
#domain： f1 f2 f3
f1，f2和f3分别是模型x，y，z方向上大小

#dx_dy_dz:
指定x，y，z方向上增量
#dx_dy_dz: f1 f2 f3
f1是空间步x方向上增量，f2是空间步y方向上增量，f3是空间步z方向上增量。空间离散化控制最大许可时间步△t和为了达成所需要仿真时间窗而提出处理方案。△t和△x，△y，△z之间关系是：
（4.1）
其中c为光速。在GPRMAX3D中，上式取等号。4.1中，一个小△x，△y，△z值造成△t得小值，这个△t小值意味着为了达成所给仿真时间而需要更多迭代次数。不过，需要指出更小△x，△y，△z和△t值会让模型更正确。

#time_step_stability_factor：
使用方法和GPR2D相同。能够修改△t值。
#time_window：
使用方法和GPR2D相同。
#number_of_media:
使用方法和GPR2D相同。
#geometry_file：
使用方法和GPR2D相同。
#messages：
使用方法和GPR2D相同。
#nips_number：
使用方法和GPR2D相同。仅当GPRMAX3D被要求用时，才必需用到她。

4.3.2 ABC 相关命令
在GPRMAX 3D这些影响Higdon ABCs配置和性能命令和GPRMAX 2D使用方法相同。不过，GPRMAX3D用了更有力PML ABC。
#abc_type: pml
#abc_type: higdon 能够设置ABC默认参数。

#pml_layers: i1
i1是pml所占Yee单元数目。默认值为8。所占用Yee单元数越多，PML性能越好，不过所花计算资源也越多。而且，有两点要注意：
a） PML是几何模型一部分。不过，PML层里域不参与计算而且用她们来计算就是错。所以，不要把源和接收器设置在这个区域。PML深度是用YEE 单元来度量而不是距离。
b） 注意：目前PML仅用于非磁介质。所以，假如你模型需要用到磁性参数，你必需用Higdoon ABC而不是PML。
4.3.3 介质和对象结构命令
在GPRMAX3D中，这些命令用来轻易引入模型中不一样介质和介质结构。在处理通常命令以后，GPRMAX3D建立一个初始化为free space（air）任意尺寸模型。

注意，free space和pec已经被定义在GPRMAX3D中,你无须再去定义这两种介质。所以，关键词free_space和pec能够直接用而不需要再定义参数。
其它介质她们参数必需经过下面之一来设置：
l #media_file：命令包含多种多样常常见介质定义（见附录A）
l #medium： 输入文件中命令
#medium：
和GPRMAX2D中使用方法相同。而且，同一个介质文件能够不改变她们结构下既用到GPRMAX2D又用到GPRMAX3D

#thin_wire：
引入细电线模型。细电线被用作介质标示符给#dege：对象结构命令。
#thin_wire： f1 str1
f1是电线半径，而且为了建立适宜细电线物理模型，必需必△l小。Str1是介质ID。Thin wire被认为是pec。
例：#thin_wire： 0.001 MyWire
#box： f1 f2 f3 f4 f5 f6 str1
f1 f2 f3是左下角坐标，f4 f5 f6是右上角坐标。Str1是#medium：定义标示符。

#plate：f1 f2 f3 f4 f5 f6 str1
f1 f2 f3是金属板左下角坐标。f4 f5 f6是金属板右上角坐标。str1 是介质标示符。

#triangle：
三角行块。
#triangle： f1 f2 f3 f4 f5 f6 f7 f8 f9 str1
f1 f2 f3，f4 f5 f6 ，f7 f8 f9分别是三角形三个顶点坐标。str1是介质标示符。

#bowtie：

蝴蝶天线。由两个三角形块组成
#bowtie：c1 c2 f1 f2 f3 f4 f5 str1
C1是蝴蝶天线方向，x，y或z。c2是剩下方向。F1，f2，f3是天线馈电点x，y，z坐标。F4是天线元长度（完整蝴蝶天线长度二分之一）。f5是展开角。Str1是介质标示符。

#edge：
一个天线。这个天线仅仅是YEE元边缘，当建立电阻或细电线是能够用她。
#edge：f1 f2 f3 f4 f5 f6 str1
F1 f2 f3是边缘起始坐标，f4 f5 f6是边缘结束坐标。Str1是介质标示符。

#cylinder：
有限维3D圆柱模型。
#cylinder：c1 f1 f2 f3 f4 f5 str1
C1是圆柱轴方向，能够是x，y或z。f1和f2分别是圆柱轴高低坐标。F3，f4是其它关键坐标用来表示两个圆柱圆形面中心。
l X方向圆柱（f1，f3，f4）为（x1，y，z），（f2，f3，f4）为（x2，y，z）
l y方向圆柱（f3，f1，f4）为（x，y1，z），（f3，f2，f4）为（x，y2，z）
l z方向圆柱（f3，f4，f1）为（x，y，z1），（f3，f4，f2）为（x，y，z2）
f5是圆柱圆盘半径，str1是介质标示符。

#cylinder_new：
有限维3D圆柱。和#cylinder：不一样是，它轴向能够是任意。
#cylinder_new： f1 f2 f3 f4 f5 f6 f7 str1
F1 f2 f3 是圆柱底面圆心坐标，f4 f5 f6是圆柱顶面圆心坐标。F7是圆柱底面半径。Str1是介质标示符。

#cylindrical_segment：
圆柱一段。
#cylindrical_segment：c1 f1 f2 f3 f4 f5 str1 c2 f6 f7
C1是圆柱轴向，能够是x，y，z。f1,f2是圆柱轴向上部和下部坐标。F3，f4是表示圆柱顶，底部其它两个坐标。F5是顶，底部半径。C2是片段改变方向。F6，f7是片段起始和结束点。
#sphere：
球体。
#sphere： f1 f2 f3 f4 str1
F1 f2 f3是球心坐标 f4是球半径。Str1是介质标示符
4.3.4 激励和输出命令
#excitation_file：

许可用户指定单个包含能够激励模型幅度值清单ASCII文件。这些值最少要和迭代次数相等。
#excitation_file：str1
Str1是ASCII文件名字。
例：#excitation_file： mysource.dat

#hertzian_dipole：
定义最简单激励。
#hertzian_dipole：f1 f2 str1 str2
F1 f2分别是源波形幅度和频率。Str1是波形标示符。Str2是源标示符。
例：#hertzian_dipole：1.0 600e6 ricker MyDipole
和GPRMAX2D里面#line_source等价

#voltage_source：
定义电压源。它引入一个电压器件位置，能够是一个硬源或一个内部集成电阻。
#voltage_source：f1 f2 str1 f3 str2
F1 f2是源波形振幅和频率。Str1是波形标示符。 F3是内部电阻R。
例：#voltage_source：1.0 600e6 gaussian 50.0 MyVolt

#transmission_line：
定义能够刺激天线1D两线传输线参数。
#transmission_line： f1 f2 str1 f3 f4 str2
F1 f2是源波形振幅和频率，str1是激励类型。F3是传输线长度。F4是阻抗特征，str2是源标示符
例：#transmission_line：1.0 600e6 gaussian 0.5 200.0 MyLine

#plane_wave：
描述平面波源
#plane_wave：f1 f2 str1 f3 f4 f5 str2 str3
F1 f2是源波形振幅和频率，str1是激励类型。F3 f4分别是0到pi和0到2pi角度。

#huygens_surface：
必需和平面波同时用。
#huygens_surface：f1 f2 f3 f4 f5 f6 str1
F1 f2 f3是huygens表面左下角坐标，f4 f5 f6是右上角坐标。Str1是要用到#plane_wave标示符
例：#huygens_surface：0.2 0.2 0.2 0.8 0.8 0.8 MyHugens

#tx：c1 f1 f2 f3 str1 f4 f5
C1源极化方向，能够是x，y，z。f1 f2 f3是源坐标。Str1是源ID，f4是激励掩饰，f5是源移去时间。

#rx： f1 f2 f3
F1 f2 f3是输出点坐标

#rx_box：f1 f2 f3 f4 f5 f6 f7 f8 f9
F1 f2 f3是输出体积左下坐标，f4 f5 f6是输出体积右上坐标，f7 f8 f9是定义每一方向输出点步数。最小值为f7为△x，f8为△y，f9为△z。

#tx_steps: f1 f2 f3
F1 f2 f3是x，y，y增量
#rx_steps: f1 f2 f3
同上

#snapshot：
取得模型区域电磁场信息。
#snapshot：i1 f1 f2 f3 f4 f5 f6 f7 f8 f9 f10 filed1 c1
或
#snapshot：i1 f1 f2 f3 f4 f5 f6 f7 f8 f9 i2 filed1 c1
i1是一个成为全局户激励点计数器。它决定哪一个snapshot被用激励源位置。她值在1和步数之间，这个步数定义为#analysis第一个参数
f1 f2 f3是左下角坐标；f4 f5 f6是右上角坐标；f7f8f9是x，y，z方向上样本间隔。
f10或i2是snapshot时间，或迭代次数。
file1是存放snapshot文件名
c1 能够是a或b表示snapshot格式是ASCII还是BINARY

GprMax中GprMax3D使用方法
很多网友问我GprMax 3D使用方法，没有那么多时间去一一叙述GprMax 3D原理。
直接贴一个GprMax 3D例子上来，给大家参考吧。
*注：这个例子也是很久之前参考网上一个例子，起源应该是起源这里：

(1)Second.in文件内容以下：
#medium: 6.0 0.0 0.0 0.01 1.0 0.0 concrete
--------------------------------------
#domain: 0.6 1.2 0.65
#dx_dy_dz: 0.01 0.01 0.01
#time_window: 12e-9
--------------------------------------
#box: 0.0 0.0 0.0 0.6 1.2 0.5 concrete
#box: 0.2 0.5 0.1 0.4 0.7 0.3 free_space
--------------------------------------
#hertzian_dipole: 1.0 900e6 ricker MyDipole
#analysis: 21 second.out b
#tx: x 0.3 0.115 0.55 MyDipole 0.0 12e-9
#rx: 0.3 0.365 0.55
#tx_steps: 0.0 0.04 0.0
#rx_steps: 0.0 0.04 0.0
#end_analysis:
--------------------------------------
#geometry_file: second.geo
#messages: y
(2)GprMax 3D界面和运行方法跟GprMax 2D相同，这里不多说了。
Matlab处理文件：

name = 'second.geo';
[mesh,ID,header,media] = gprmax3g(name);
figure
p1=patch(isosurface(mesh,1));
set(p1,'facecolor',[1,0,0],'edgecolor','none');
p2=patch(isosurface(mesh,2));
set(p2,'facecolor',[0,1,0],'edgecolor','none');
[Header, Fields] = gprmax('second.out');
radargram(:,:) = Fields.ez(:,1,:);
t=1e9*Fields.t;
Y = Header.dy * Header.rx: Header.dy * Header.RxStepY :Header.dy * Header.rx + (Header.dy * Header.RxStepY * (Header.NSteps-1));
figure;
imagesc(Y,t,radargram);
title('Ez');
ylabel('Time (ns)');
xlabel('Scan length (m)');
figure;
[x,y]=meshgrid(Y,t);
surf(x,y,radargram,'EdgeColor','none')
view([1,1,0.5]);
title('Ez' );
ylabel('Time (m)');
xlabel('Scan length (m)');
处理结果图像：（见下页）
(3)问题：
geo文件显示不正确，没时间去找为何了，有空再弄吧。