简介gprMax 是基于有限差分时域FDTD方法求解麦克斯韦方程的开源电磁波模拟软件主要用于探地雷达GPR正演建模也可用于天线与土壤介质等电磁仿真场景。压缩包面向刚接触 GPR 数值模拟的研究生、科研人员与工程技术人员提供可运行的源码与配套算例帮助完成从建模、仿真到结果可视化的入门闭环。包内共 373 个文件压缩后约 33.1MB。其中 59 个 py 脚本与 11 个 pyx、2 个 pxd 文件承载求解器与 Python 接口30 个 rst 文档和 12 份 pdf 构成使用手册与理论说明64 个 .in 文件是天线及典型场景的建模算例另有 png 图像、npz/npy 数据、vtp/vti 可视化结果和 sh、makefile、yml 等构建配置便于按目录定位源码、算例与文档。已有 3211 人学习下载。读者可借此了解 CPU 并行与 GPUCUDA两套求解器的组织方式参照 GSSI 1500、MALA 1200 等天线算例改造模型并用已有输出文件对照验证结果。1. gprMax 能替你先跑一遍雷达从模糊双曲线到可复现正演现场采回来的雷达剖面上挂着一串双曲线标注之前总得先问一句这到底是管线、空洞还是层界面gprMax 就是拿来回答这类问题的开源 FDTD 软件用有限差分时域方法把探地雷达的正演过程跑一遍。给定介电常数、电导率、几何体和天线位置它按 Yee 网格一步步推进电磁场输出能拿去和实测对比的 A-scan、B-scan。软件本体由 Python 驱动核心循环用 Cython 编译支持 Debye/Lorentz 色散土壤、任意几何建模、多天线收发以及 OpenMP、MPI 和 CUDA 三种加速路线。做 GPR 资料解释、天线频带设计、教学验证的工程师和研究生都能用得上想一键出图的人会失望模型得自己一行行搭。2. 环境与第一条 A-scanconda 安装、GPU 开关和 300 道 B-scan先把“跑起来”这件事拆成三步装对包、写对一个最小模型、搞清命令行参数。三步里最容易翻车的是第一步很多人卡在编译扩展上后面干脆放弃。下面的顺序是按依赖从弱到强排的先把 OpenMP 版跑通再碰 GPU。2.1 安装路线OpenMP 版先跑通再谈 CUDAgprMax 的 Python 包依赖 Cython 编译的扩展模块直接pip install在 Windows 上大概率卡在编译器找不到头文件。常见做法是走 conda-forge它把编译好的二进制和运行时一起装下来。# 建一个独立环境避免和系统里的 numpy、h5py 打架 conda create -n gprmax python3.10 -y conda activate gprmax # conda-forge 的包自带编译好的 Cython 扩展和 OpenMP 运行时 conda install -c conda-forge gprmax -y # 确认模块能导入、命令行能起来 python -c import gprmax; print(gprmax.__file__) python -m gprMax --help这几条命令的逻辑很简单第一条隔离环境第二条装包后两条做冒烟测试。--help里要重点看两个东西——有没有-mpi选项决定你这套包是不是 MPI 构建有没有-gpu决定 CUDA 路径通不通。如果--help正常但一跑模型就报ImportError: cannot import name cython之类基本是环境里混进了 pip 版的旧包删环境重装比修快。GPU 路线要额外的 CUDA 工具链和 pycuda# CUDA 驱动由系统提供这里只装 Python 侧的绑定 conda install -c conda-forge pycuda -y # 初始化设备并打印型号能出结果才说明驱动、运行时、绑定三层对齐 python -c import pycuda.autoinit, pycuda.driver as drv; print(drv.Device(0).name())pycuda.autoinit会在导入时自动创建上下文这一句能跑通-gpu才有意义。版本对齐是这里唯一的玄学CUDA 运行时大版本必须和驱动支持的上限匹配装完跑不通就先nvidia-smi看驱动上限再回退 pycuda 的构建版本。MPI 则是另一条路conda 包不一定带--help里没有-mpi就别试硬加参数只会报未知选项。注意Windows 上 MPI 和编译工具链的配置成本很高我一般直接在 Linux 或 WSL 里做省下来的时间够跑十几组参数扫描。2.2 建第一个 .in砂土里埋一根金属管gprMax 的模型文件是纯文本一行一个指令#开头。下面这个模型是 2D 近似0.6 m × 0.6 m 的砂土半空间中间埋一根直径 6 cm 的金属管理想赫兹偶极子做发射沿 x 方向步进采集。#title: sand pipe demo #domain: 0.6 0.6 0.002 #dx_dy_dz: 0.002 0.002 0.002 #time_window: 8e-9 #pml_cells: 10 10 10 10 0 0 #material: 6 0.005 1 0 sand #material: 1 1e7 1 0 metal #box: 0 0 0 0.6 0.6 0.002 sand #cylinder: 0.20 0.30 0 0.20 0.30 0.002 0.03 metal #waveform: ricker 1 1.2e9 my_ricker #hertzian_dipole: z 0.05 0.10 0 my_ricker #rx: 0.05 0.10 0 #src_steps: 0.002 0 0 #rx_steps: 0.002 0 0 #geometry_view: 0 0 0 0.6 0.6 0.002 0.002 0.002 0.002 pipe_demo n逐行说明关键参数。#domain是物理尺寸单位米z 方向只给 0.002 是为了做 2D 近似网格在 z 上只有一层。#dx_dy_dz是空间步长0.002 m 对应中心频率 1.2 GHz 时每波长约 40 个网格余量足够。#time_window是记录长度8 ns 大约能覆盖 0.6 m 深度内砂土中的双程走时。#pml_cells六个数字对应 x-、x、y-、y、z-、z 六个方向的吸收层厚度z 方向只有一层网格所以填 0 而不是 10。#material四个参数依次是相对介电常数、电导率S/m、相对磁导率、磁损耗。砂土取 εr6、σ0.005 是干燥砂的常见量级金属用 εr1、σ1e7 近似理想导体。#box和#cylinder的后三个参数是最小角和最大角坐标圆柱还要在坐标后跟半径。#waveform定义一个 Ricker 子波两个数字是振幅和中心频率最后是名字。#hertzian_dipole的方向z表示极化沿 z后面是位置和波形名。#src_steps和#rx_steps让收发天线每次 A-scan 后沿 x 移动 2 mmgprMax 会一直步进到天线接近 PML 边界为止。提示#geometry_view最后那个n表示输出 vtk 而不是 HDF5 体素跑之前先看一眼几何比跑完发现管子位置错了再回头省事得多。2.3 命令行参数与并行开关别把 -n 和 -mpi 搞混# 单次运行输出到指定目录 python -m gprMax user_models/pipe_demo.in -n 1 -o output/ # 只生成几何、不推进时间步用来检查建模有没有越界 python -m gprMax user_models/pipe_demo.in -geometry-only # GPU 加速-n 表示重复运行次数-gpu 表示走 CUDA 路径 python -m gprMax user_models/pipe_demo.in -n 1 -gpu -o output/ # 把同一模型的多道 A-scan 合并成一个 B-scan 文件 python -m tools.outputfiles_merge output/pipe_demo-n是重复运行的次数用于同一模型换随机种子做统计不是并行线程数这点最容易搞混。真正的并行有两套OpenMP 版通过环境变量OMP_NUM_THREADS控制线程数MPI 版用-mpi 4指定进程数但前提是这套包本身是 MPI 构建。-gpu走单卡和-mpi互斥同时给会直接报错。-geometry-only是个被低估的开关几何体越界、材料名拼错、PML 和结构重叠都能在这一步暴露出来不用等几十分钟的时间步推进。outputfiles_merge这个工具按文件名前缀匹配所以同一目录下别放多个同名模型的不同参数版本否则合并出来的道序会乱。跑完先看终端输出的网格数、时间步数和预计内存占用网格数乘上 6 个场分量再乘单精度 4 字节就是 GPU 上最保守的显存下限超过显存会直接抛 CUDA 错误而不是降级到 CPU。3. 模型文件怎么写几何叠加、Debye 色散与天线激励最小模型跑通只证明工具链没坏真正决定仿真有没有物理意义的是几何怎么叠、土壤色散怎么给、天线怎么激励。这一章三件事按建模顺序排先把结构摆对再给材料参数最后处理源和天线。3.1 几何体的覆盖顺序与粗糙面gprMax 的几何是“后写覆盖先写”的规则同一个位置后添加的物体材料会替换掉先前的。所以标准写法是先铺大背景再放目标体最后处理表面。#box: 0 0 0 0.6 0.6 0.15 sand #add_surface_roughness: 1.5 0.003 0.05 7 #box: 0 0 0.15 0.6 0.6 0.6 clay #cylinder: 0.30 0.30 0.15 0.30 0.30 0.16 0.05 metal#add_surface_roughness作用在最近一个#box的上表面四个参数是分形维数、均方根高度、相关长度、随机种子。第一段先铺 0.15 m 厚的砂层并给它一个粗糙顶面第二段在粗糙面之上填黏土第三段把金属管埋进界面附近。覆盖顺序决定了最终材料分布如果把背景 box 写在目标体之后目标体会被整块覆盖掉几何视图里只剩一片均匀介质这是新手最常见的“模型跑完什么都没有”的原因。粗糙面的参数不是随便填的均方根高度要显著小于波长否则界面散射会强到把目标回波盖住相关长度决定粗糙度的空间尺度太短会引入网格级锯齿。种子固定住同一组参数才能复现做参数对比时别每次换种子。3.2 土壤色散参数怎么定Debye 极点与电导率干砂当成非色散介质处理误差不大含水土壤就必须要色散模型否则波形拖尾和频散对不上实测。gprMax 用#add_dispersion_debye给材料挂 Debye 极点。#material: 4 0.01 1 0 wet_sand #add_dispersion_debye: 1 8 0.4e-9 wet_sand #material: 15 0.05 1 0 clay #add_dispersion_debye: 1 6 1.2e-9 clay#add_dispersion_debye的三个参数是极点数量、介电常数增量 Δε、弛豫时间 τ。上面第一组表示湿砂在静态介电常数 4 的基础上再叠加 8总静态 εs 约 12弛豫时间 0.4 ns第二组黏土总量级更高、弛豫更慢。工程上常取的量级如下表具体场地还是得以取样或 TDR 实测为准。介质相对介电常数 εr电导率 σ (S/m)典型 τ干砂4 ~ 60.001 ~ 0.01可忽略湿砂10 ~ 200.01 ~ 0.050.2 ~ 0.6 ns粉质黏土15 ~ 250.02 ~ 0.10.5 ~ 1.5 ns饱和黏土25 ~ 350.1 ~ 0.51 ~ 2 ns参数给错的典型表现是目标回波位置对但相位反了或者 B-scan 上出现整体下移的拖尾。前者多半是 εr 偏差导致波速不对后者往往是漏了电导率或弛豫时间设得太长。色散极点数量不是越多越好单极点拟合在 GPR 频段内通常够用多极点会显著增加每步的内存访问GPU 上尤其明显。3.3 从理想偶极子到真实天线馈电与匹配电阻#hertzian_dipole是理想点源做层位验证足够但它的辐射波形和真实天线差得远尤其是振铃和方向图。要复现实测得自己搭天线结构。#edge: 0.10 0.10 0 0.10 0.20 0 pec #edge: 0.12 0.10 0 0.12 0.20 0 pec #waveform: gaussian 1 1.0e9 my_pulse #excitation_file: user_models/tx_voltage.txt my_pulse #add_lumped_element: 0.10 0.10 0 0.10 0.12 0 50 resistor#edge用细线近似金属条带两端坐标加材料名这里是理想导体。#excitation_file从外部文本读激励波形适合把实测馈电电压灌进去。#add_lumped_element在指定线段上加集总元件最后一个参数是阻值50 欧姆是常见匹配值。天线建模的血泪经验是不加匹配电阻馈电点会来回反射时域波形上挂一串等间隔振铃看着像浅层多次波实际是天线自己振。加了电阻振铃压下去但辐射效率也降了这个折中要在频带和目标深度之间权衡。天线模型还有个体积问题细线结构在网格里是零厚度实际金属条带厚度接近网格尺寸时要么加密网格要么在#edge外面再包一层薄#box否则会低估天线等效长度主频偏移个百分之十几很常见。4. 避坑与排查五个让仿真白跑的细节这一章的五条都是我踩过或者帮别人排查过的共同点是模型能跑完、文件也生成了但结果不能用。每条按现象、原因、解决写。4.1 时间窗截断回波被切在记录末端现象是 B-scan 上目标双曲线只画了一半或者深部层的同相轴在记录末尾突然消失。原因十有八九是#time_window不够最大双程走时没有覆盖最深的界面。c 3e8 eps_r 12.0 # 土壤相对介电常数 depth 0.6 # 最深处埋深单位米 v c / eps_r ** 0.5 # 介质中波速 t_max 2 * depth / v # 双程走时 print(f最短时间窗建议 {t_max * 1.5e9:.1f} ns至少留 50% 余量)这段计算里波速按非色散近似给色散介质实际群速会略低所以乘 1.5 的余量不是保守是必要。如果目标上方还有低损耗层天线和界面之间会来回多次反射记录尾部会被多次波占满这时候要么加时间窗要么在接收端加时域窗把晚到的能量压掉。改完时间窗记得同步检查输出文件大小网格不变的情况下记录变长显存和磁盘占用一起涨。4.2 网格色散每波长不足 10 格现象是波形尾部出现高频振荡或者同一个模型加密网格后目标回波位置移动了。原因是空间步长相对最小波长太粗FDTD 的数值色散把波形拖变形了。c 3e8 f 1.5e9 # 最高有效频率不是中心频率 eps_max 20.0 # 模型中最高的相对介电常数 lam_min c / (f * eps_max ** 0.5) print(f最小波长 {lam_min*100:.1f} cm建议 dx {lam_min/10*1000:.1f} mm)注意这里的频率要用频带高端而不是中心频率Ricker 子波在中心频率两倍附近还有可观能量。判断标准是每个最小波长至少 10 个网格做定量对比时给到 15 到 20 更稳。加密网格的代价是网格数按立方增长2D 模型还好3D 模型从 2 mm 加密到 1 mm内存直接乘 8。折中办法是保持粗网格验证几何再用局部细网格做重点区域或者干脆把模型切小。4.3 PML 与几何重叠边界反射伪装成目标现象是 B-scan 两侧出现倾斜的强同相轴随天线移动而移动位置看着像浅层目标。原因是 PML 吸收层和几何体重叠了或者 PML 厚度不够、吸收不干净。解决方法分两步先用#geometry_view把模型切出来确认所有几何体的最大角坐标离 domain 边界至少有 PML 厚度加几个网格的距离再看 PML 层数GPR 频段一般 10 层起步介电常数高的介质要加到 15 层。发射和接收天线如果落在 PML 里问题更严重波还没传播就被吸收掉所以#src_steps的起点要留出足够余量。这一条在 2D 模型里容易被忽略因为 z 方向 PML 设了 0x、y 方向一旦贴边反射照样进来。4.4 直达波淹没目标收发距与天线振铃现象是 B-scan 上一整条水平强轴目标双曲线埋在它下面看不清。原因有两种收发距太小直达波和目标回波在时间上重叠或者天线没匹配振铃拖尾太长。先调收发距把接收器从 0.10 挪到 0.20 甚至更远用空间换时间分离再不行就在后处理里做背景去除用逐道减去平均道把不随位置变化的直达波压掉但这一步会同时削掉水平层界面用之前想清楚目标类型。天线振铃的根治还是匹配电阻和激励波形整形后处理只能治标。做定量对比时直达波的位置本身也是校验项它和理论到时对不上说明源点或时间零点有问题。4.5 GPU 显存与 -mpi 混用现象是 GPU 跑到一半报显存不足或者同一个模型 CPU 能跑、GPU 直接退出。原因通常是显存估算不足或者把-gpu和-mpi一起给了。显存下限粗略按网格数乘场分量再乘系数估3D 模型还要乘上色散极点带来的附加存储。跑到一半才爆的多半是输出文件累积或者多道记录同时在显存里-o指到磁盘而不是内存盘能缓解一部分。-gpu和-mpi互斥这一点报错信息很明确但有人为了“用满机器”硬凑参数最后得到的是进程互相抢显存。单卡跑大模型时把-n拆成多次小批次跑再在合并阶段拼起来比一次性硬扛稳得多。5. 输出文件解析HDF5 结构、B-scan 成像与时间零点校正跑完的.out文件是 HDF5 格式不打开看一眼很容易把它当黑匣子。这一章从文件结构讲到成像最后用解析解校一遍幅值。5.1 .out 里到底存了什么import h5py import numpy as np path output/pipe_demo_merged.out with h5py.File(path, r) as f: print(顶层组:, list(f.keys())) # 通常是 rxs print(接收器:, list(f[/rxs].keys())) # rx1, rx2 ... rx1 f[/rxs/rx1] print(场分量:, list(rx1.keys())) # Ex, Ey, Ez, Hx ... ez rx1[Ez][:] # 形状 (道数, 时间采样) dt f.attrs[dt] # 时间步长单位秒 dx f.attrs[dx_dy_dz] # 空间步长 print(Ez 形状:, ez.shape, dt:, dt, dx:, dx)合并后的文件把多道 A-scan 堆在第一个维度上第二个维度是时间采样。dt存在文件属性里不要自己按#time_window除采样数去猜迭代步长和记录长度不是整除关系。单道文件和合并文件的区别只在第一个维度长度读法一样。如果Ez是空的检查#rx指令里有没有写分量gprMax 默认记录全部电场分量但合并工具只搬它认识的组。读到数组之后先看首道的直达波位置这是后面做时间零点校正的基准。直达波到时约等于收发距除以介质波速对不上就回头查源点和接收点坐标。5.2 后处理脚本背景去除、时间零点与增益import h5py import numpy as np import matplotlib.pyplot as plt with h5py.File(output/pipe_demo_merged.out, r) as f: ez f[/rxs/rx1/Ez][:] dt f.attrs[dt] # 1) 时间零点校正以首道直达波峰值为基准整体平移 first_trace np.abs(ez[0]) shift int(np.argmax(first_trace)) - int(0.1e-9 / dt) ez np.roll(ez, -shift, axis1) # 2) 背景去除逐道减去所有道的平均道压制水平直达波 ez_bg ez - ez.mean(axis0, keepdimsTrue) # 3) 增益按时间指数补偿介质衰减和球面扩散 ns ez.shape[1] t np.arange(ns) * dt gain np.exp(t / (3e-9)) ez_gain ez_bg * gain[None, :] plt.figure(figsize(8, 5)) plt.imshow(ez_gain.T, aspectauto, cmapgray, extent[0, ez.shape[0], t[-1]*1e9, 0], vmin-np.percentile(np.abs(ez_gain), 98), vmaxnp.percentile(np.abs(ez_gain), 98)) plt.xlabel(道号) plt.ylabel(双程走时 (ns)) plt.colorbar(labelEz (V/m)) plt.tight_layout() plt.savefig(output/bscan.png, dpi150)三步处理的顺序不能换。时间零点必须最先做否则平均道里各道的直达波错位减完反而留下残差。背景去除会削掉所有水平同相轴层状界面也会被压如果目标就是水平层跳过这一步改用中值滤波。增益的指数常数按介质衰减估湿地取小一点干砂取大一点给太猛会把噪声一起放大双曲线看着更清楚但幅值信息全丢了。extent里时间轴反向是为了让走时向下增大和商业软件显示习惯一致。5.3 用 Fresnel 反射系数验证幅值import numpy as np eps1, eps2 6.0, 15.0 # 上层砂、下层黏土 r (np.sqrt(eps1) - np.sqrt(eps2)) / (np.sqrt(eps1) np.sqrt(eps2)) print(f界面反射系数 {r:.3f}功率反射 {r**2:.3f})这是垂直入射的 Fresnel 公式用来做量级校验。做法是在同一模型里放一个水平界面取界面回波的峰值除以入射波峰值和公式结果比。差在 10% 以内算正常差一倍以上通常是三个原因天线方向性导致入射角不是垂直、界面附近网格太粗、或者背景去除把界面能量一起削了。这个校验不需要额外工具但它能区分“模型跑通了”和“模型算对了”后者才是能拿去做反演和解释的基础。6. 收敛性验证与参数扫描把仿真当实验台用仿真结果能不能信取决于你有没有做过收敛性测试。我的习惯是同一几何用两套网格跑看波形是否重合。# 粗网格 python -m gprMax user_models/pipe_demo.in -n 1 -o output/coarse/ # 把 #dx_dy_dz 改成一半domain 和几何坐标保持不变 python -m gprMax user_models/pipe_demo.in -n 1 -o output/fine/两套结果的差异看两个指标目标回波峰值的相对误差、双程走时的偏移。峰值差在 5% 以内、走时差在半个采样点以内就可以认为当前网格够用。加密后差异反而变大多半是 PML 层数没跟着加或者色散参数在细网格下暴露了拟合误差。验证项粗网格加密一倍判定标准目标回波峰值基准对比相对误差 5%双程走时基准对比偏移 0.5 个采样点直达波振铃长度基准对比不随网格显著变化边界反射幅值基准对比低于目标回波 20 dB参数扫描用.in文件里的#python:块比手改文件可靠循环变量直接进几何坐标和材料参数跑完一批再用脚本批量读输出。#python: import numpy as np for eps in [6.0, 10.0, 15.0]: print(f#material: {eps} 0.01 1 0 scan_soil) #end_python:#python:块在解析阶段执行print出来的内容会被当成模型指令插入所以它能动态生成材料、几何甚至整个扫描列表。用它的前提是变量名不和外层冲突且不要在块里做耗时计算解析阶段卡住会让人误以为是求解慢。从那以后我每次新建模型都先跑一个只有半空间和 PML 的空场景确认直达波什么时候到、边界有没有回波再往上叠几何和色散。这一步十分钟能省掉后面几小时的返工。希望帮到你。本文还有配套的精品资源点击获取