简介EM3DVP 是一套基于 Matlab 的三维地电磁建模与反演可视化软件包面向从事电磁法研究的地球物理科研人员重点解决三维反演流程中模型构建、数据准备与结果展示的繁琐问题。包内共 243 个文件以 239 个 m 脚本为主体涵盖数据读写、模型生成、反演参数保存及响应绘图等核心功能另附说明文档、许可证和示意图压缩包仅 741KB结构精炼便于快速部署。目前已有 231 人浏览学习适合具备一定 Matlab 操作基础、需要开展结构化网格电磁建模与反演任务的用户。借助该工具用户可以方便地导入 EDI 等电磁数据在图形界面中设置模型参数、生成反演输入文件并对模拟响应和切片图进行直观解读从而显著提升三维地电磁数据处理与结果分析效率。1. EM3DVP 在三维地电磁建模反演里解决什么问题在三维地电磁勘探项目里正演的物理问题已经被有限差分和有限元程序解决了很多年真正卡住人的往往是可视化模型是怎么建起来的异常体埋在哪个深度反演出的电阻率体是不是被边界或地形污染了。EM3DVP 就是专门把网格剖分、属性赋值、正演响应、迭代反演和结果成图接进同一套可视界面里的软件包。它把三维地电磁建模反演的中间状态都变成可以看、可以查、可以改的对象而不是一堆只在文本日志里出现的数值。适合三类人做 MT 或 CSEM 解释的工程师要批量测试初始模型的研究生以及需要向前方项目组交付可复核电阻率体的勘查团队。EM3DVP 的价值不在新算法而在把所有环节放进同一个可视化闭环里让你在提交结果前就看出问题在哪。2. EM3DVP 的模型构造从地质解释到计算网格2.1 为什么 EM3DVP 把模型当体素而不是面模型地电磁正演在求解 Maxwell 方程时电阻率在每个网格单元内部是常量。因此 EM3DVP 里的核心对象不是三角网也不是地质层面模型而是一个三维正交网格上的属性数组。很多从专业地质建模软件转过来的用户不适应这一点觉得断层、尖灭没法表达。常见做法是先把地质解释结果栅格化再对每个网格单元赋电阻率。这样做的第一个好处是网格与求解器天然一致不再需要二次离散第二个好处是可视化时能按深度逐层切异常体边界与观测网格一一对应。缺点也明显复杂构造会被阶梯状近似所以建模前必须想清楚网格尺度。我一般把中心测区网格边长设为最小测点距的 0.51 倍外围网格按 1.21.5 倍比例向外扩张直到模型外边界超过目标频率对应趋肤深度的 23 倍。这里要注意 z 轴惯例EM3DVP 通常把地表以下作为正方向空气层单独放在最上层不要把空气层和地下单元混在一起赋值。2.2 用 Python 生成 EM3DVP 属性体一个可跑的示例下面这段脚本用 numpy 构造一个 50×50×30 的三维电阻率体在背景中埋入一个低阻立方体。它的运行结果可以直接作为 EM3DVP 可视化导入的数据源也是后续正演和反演的起算模型。import numpy as np # 网格单元尺寸东向、北向、垂向单位米 dx, dy, dz 40.0, 40.0, 25.0 nx, ny, nz 50, 50, 30 # 背景电阻率 1000 欧姆米单位为欧姆米 model np.full((nz, ny, nx), 1000.0, dtypenp.float32) # 生成三维索引网格 zs, ys, xs np.mgrid[0:nz, 0:ny, 0:nx] # 在水平中心、深度第 13~15 层放一个 30 欧姆米低阻体 anomaly ( (np.abs(xs - 25) 5) (np.abs(ys - 25) 5) (np.abs(zs - 14) 2) ) model[anomaly] 30.0 # 保存为 EM3DVP 可读的 npz 文件 np.savez( em3dvp_model.npz, resistivitymodel, origin(-1000.0, -1000.0, 0.0), cell_size(dx, dy, dz) )逻辑说明代码用布尔数组anomaly圈定目标区域比逐个单元赋值更快也更容易和地质图件对照。如果要放倾斜体可以先构造水平体的布尔掩模再用scipy.ndimage.rotate旋转但旋转后必须做最近邻重采样否则会产生半格点。origin记录模型西南角坐标cell_size决定每个单元的地球物理尺寸这两项在可视化定位测点时非常重要。若 EM3DVP 的导入界面只接受 netCDF 或文本格式把model数组按深度逐层写成文本行即可关键是保持层序从地表向下递增。2.3 空气层和外边界可视化里看不见却决定成败的部分EM3DVP 里建立模型不只是导入地下电阻率还要定义空气层。空气层的作用是让地表电流回路有足够的空间而不是真实模拟大气电性。常见设置是电阻率取 1e12 Ω·m 的近似绝缘体模型顶部至少加 1020 层空气网格。为了省内存删掉空气层会让地表响应产生明显畸变反演结果在浅层出现虚假导体。外边界距离通常用趋肤深度做参考δ ≈ 503 · sqrt(ρ / f)其中 ρ 是围岩电阻率单位 Ω·mf 是频率单位 Hz。外边界建议取 23 倍 δ。以电阻率 1000 Ω·m、频率 0.3 Hz 为例δ 接近 29 km外边界至少 60 km。不同场景的参考值如下频率围岩电阻率趋肤深度 δ外边界建议10 Hz100 Ω·m1.59 km≥ 4 km1 Hz1000 Ω·m15.9 km≥ 35 km0.1 Hz1000 Ω·m50.3 km≥ 110 km0.01 Hz1000 Ω·m159 km≥ 320 km这个外边界不是地下的深度而是测区四周的延展范围。网格从中心向四周以对数方式拉长屏幕上看起来非常空洞但正是这一段没人关心的空白区域决定了 TE 模式下电流回路是否完整。提示在 EM3DVP 中把空气层设为透明仅保留地表网格显示能明显降低场景渲染负担同时不会影响模型检查。3. EM3DVP 正反演参数设计频率、误差下限与迭代控制3.1 先确定频率与极化模式再谈反演三维地电磁模拟通常分两类MT 平面波模拟和 CSEM 可控源模拟。EM3DVP 对前者按频率逐点计算对后者还要定义源位置、源长度和收发距。频率选择不一定要密按对数间隔每十倍频取 35 个频率即可。频率太少深度分辨率不够频率太多计算时长成倍增加。极化模式方面MT 正演要计算两个相互垂直的极化方向用于形成完整的阻抗张量CSEM 则要计算多个源的响应并把电场、磁场分量分别做归一化。可视化软件包通常会把四种曲线画在一起方便你看相位和振幅是否连续。如果后续反演只用标量数据记得在预处理阶段把两条极化曲线合并成平均阻抗避免一会儿用 TE、一会儿用 TM 造成数据不协调。3.2 反演参数表先跑默认值再改三个关键项把一次反演的主要参数放在一张表里前 4 项每一个都要结合数据质量判断参数含义常见初值什么时候调error_floor数据误差下限百分比2%数据噪声明显时升到 5%reg_factor正则化权重0.03迭代发散时增大但不要超过 0.3smooth_x/y/z光滑约束尺度0.5 / 0.5 / 0.2横向延续性好时降低 z 向光滑max_iteration最大迭代步数30前 5 步 RMS 平缓时可减到 20n_cpu并行线程数8和内存带宽匹配不要只追核数说明error_floor是百分比不是欧姆米绝对值。它决定每个数据点的权重设太低会让反演拼命拟合噪声图像出现孤立斑块设太高则把浅层细节全部抹平。reg_factor影响模型光滑度与数据拟合的平衡前几次迭代内需要人工观察等 RMS 进入平稳段后再恢复默认。3.3 可复现的反演命令骨架下面用三段式命令把建模、正演、反演串起来具体子命令以你安装的 EM3DVP 版本帮助为准# 生成工程绑定观测文件和初始模型 em3dvp create project --survey obs.edi --model em3dvp_model.npz --out m3 # 正演一次性计算 1.0、0.3、0.1 Hz 三个频点 em3dvp forward m3/project.h5 \ --frequencies 1.0 0.3 0.1 \ --mode tmte \ --log forward.log # 反演从初始模型出发带误差下限和正则化权重 em3dvp inverse m3/project.h5 \ --obs obs.edi \ --init em3dvp_model.npz \ --error-floor 0.02 \ --reg 0.03 \ --max-iter 30 \ --log inverse.log命令说明create project把 EDI 文件和模型文件绑定在同一个工程中EDI 里包含测点坐标、频率和阻抗张量没有 EDI 时也可以从 CSV 构造观测文件。--mode tmte同时计算 TM 和 TE这是 MT 的标准配置可控源场景应该改成源列表。反演的--init与正演模型可以不同实际项目里经常用均匀半空间或一维反演结果作为初始模型。--error-floor 0.02表示 2%、不是 0.02 欧姆米这个误区最容易让新手觉得参数没生效。3.4 迭代监控别只盯 RMS还要看模型变化量反演日志里通常同时输出数据拟合残差 RMS 和模型变化量 model change。只看 RMS 小于 1 就停止容易把数值振铃当成异常体收进最终模型。更稳的做法是看 model change如果相邻迭代模型的相对变化超过 20%说明模型主体还在被改写如果低于 5% 且 RMS 不再下降可以结束迭代。下面是一个从日志里提取指标的小脚本import re def check_log(path, change_tol0.05): last_change None with open(path) as f: for line in f: m re.search( riter(\d) rms([\d.]) change([\d.]), line ) if not m: continue it, rms, chg int(m[1]), float(m[2]), float(m[3]) print(fiter{it:3d} rms{rms:.3f} change{chg:.3f}) last_change chg return last_change is not None and last_change change_tol print(converged:, check_log(inverse.log))正则表达式按iter... rms... change...的格式提取三列。change是相邻两次迭代的模型相对变化小于 0.05 时函数返回True提示可以收手。总迭代次数多不代表模型可靠在 EM3DVP 里旋转三维视图往往比孤立的 RMS 曲线更能发现深层不合理的柱状体。4. 用 EM3DVP 查看反演结果切片、断面与数据导出4.1 深切片脚本把电阻率体转成图件反演完成后得到的不只是几条剖面而是整个三维电阻率体。在 EM3DVP 中拖动深度滑块看平面图是最快的检查方式但交付图件需要用脚本按指定深度输出import numpy as np import matplotlib.pyplot as plt data np.load(m3/inv_model.npz) rho data[resistivity] # 形状 (nz, ny, nx) dz data[cell_size][2] # 垂向单元厚度米 for k in [3, 8, 12, 20]: depth k * dz fig, ax plt.subplots(figsize(6, 6)) im ax.imshow( np.log10(rho[k]), originlower, cmapSpectral_r ) ax.set_title(f{depth:.0f} m below surface) fig.colorbar(im, labellog10(resistivity / Ohm*m)) fig.tight_layout() fig.savefig(finv_slice_{k:02d}.png, dpi200)originlower保证北方向朝上和 EM3DVP 的平面视图保持一致否则切片图会上下颠倒。色标取log10而不是原始电阻率因为地电磁响应的视电阻率常常跨越三个数量级线性色标会把 30 Ω·m 的异常并进背景。文件名带上深度层号后续用 GIS 做图册时才不会搞错顺序。4.2 在 EM3DVP 里叠加测点位置检查异常体边界光切深度还不够实际项目里要判断异常体边界需要把测点位置叠加到切片上。一般在 EM3DVP 的图层管理器里添加 survey stations 图层点符号设为空心圆圈再与电阻率切片叠加。位于测点稀疏区的异常通常是插值假象反演对该区域没有数据约束。要特别关注测线之间的画饼效应。可视化软件包为了显示平滑会在网格之间做插值但导出数据时仍保留真实反演单元值。如果交付图件上出现非常圆滑的高阻透镜先检查它是否落在相邻测线的空白区。可以用网格加密前后两版模型做差差异最大的地方往往就是数据约束最弱的区域。4.3 把局部坐标模型转成 GeoTIFF交付给钻井或工程方的时候对方通常只要几张带地理参考的 GeoTIFF而不是一个 h5 或 npz 文件。常见做法是用 GDAL 写出仿射变换from osgeo import gdal import numpy as np rho np.load(m3/inv_model.npz)[resistivity] origin_x, origin_y -1000.0, -1000.0 res 40.0 driver gdal.GetDriverByName(GTiff) ds driver.Create( res_model_depth_500m.tif, rho.shape[2], rho.shape[1], 1, gdal.GDT_Float32 ) ds.SetGeoTransform((origin_x, res, 0, origin_y, 0, -res)) ds.GetRasterBand(1).WriteArray(rho[4]) ds.FlushCache()说明SetGeoTransform的六个参数依次是左上角东向坐标、x 分辨率、x 旋转系数、左上角北向坐标、y 旋转系数、y 分辨率。反演体通常以局部直角坐标存储转成全球坐标前先确认单位是米还是经纬度UTM 坐标可以直接写经纬度则需要先投影到平面。4.4 可视化验证三件套响应拟合、L 曲线、伪数据验证项看什么合格线响应拟合观测阻抗 vs 预测阻抗落在误差带内曲线趋势一致L 曲线模型光滑度 vs 数据拟合残差取拐点附近的正则化参数伪数据图正演预测响应与真实观测的差异无相干条纹或系统偏移伪数据图是最容易被跳过的一步。EM3DVP 可以正演当前模型得到伪响应把伪响应和真实观测做对比。如果伪数据图里出现明显条带说明模型有些结构只是为拟合噪声而生深切片再好看也不能直接交付。5. EM3DVP 跑批时最该盯的 4 个问题网格、空气层、初始模型与数据权重5.1 网格加密区必须覆盖测点而不是异常体常见误解是异常体在哪网格就在哪加密。但可控源和 MT 反演都是通过测点处的灵敏度更新模型测点稀疏区的网格再细也没有足够数据支撑只会增加内存和并行调度负担。我一般把测点周边 23 倍测点距的范围设为细网格区异常体位置只影响属性赋值不影响网格密度。检查方法在 EM3DVP 中打开网格剖分线让所有测点落在细网格内部并且细网格边界到最近的测点至少保留一个单元。否则沿测网边缘会出现条带状高阻或低阻异常这类假象最常被解释成构造边界。5.2 空气层电阻率不是越大越好空气层取 1e81e12 Ω·m 都可以正常工作。如果把空气层电阻率设为 1e20线性求解器容易出现数值病态迭代次数反而增加。地表附近网格纵横比也不要超过 5:1否则误差集中在表层浅层切片会出现黑白交替的棋盘格。在 EM3DVP 界面里空气层往往默认显示为透明蓝色体。很多用户拖动视角时误以为它是高阻异常体直接删掉或者让它参与反演更新。正确的做法是把它设成一个独立分组并在反演参数中锁死该分组保证正演时保留、反演时不更新。5.3 初始模型选均匀半空间还是分层模型均匀半空间最容易收敛但会把浅层静态位移或局部噪声反演成假低阻体。如果想压低这种源效应可以先用一维反演拟合每个测点的视电阻率曲线再把插值后的分层模型作为 EM3DVP 的初始模型。注意分层模型的层界面要尽量与网格层中心对齐。如果层界面落在网格层内部反演求解器会额外引入平滑梯度模型会认为这里有真实电阻率变化延长迭代收敛时间。5.4 数据权重里的坏道比正则化更容易引发面条型异常反演对数据权重非常敏感。同一个频点如果观测误差被低估 5 倍它的权重是正常点的 25 倍模型会优先拟合这几个点结果就是在深切片上出现垂直细条或倒三角异常。所以跑批前先做数据质量打分对每个测点、每个频率统计重复观测的相对误差和相位相关性。我一般会先设error_floor0.02再把明显脱群点的误差统一抬高到不小于 0.05。和直接删掉坏道相比抬升误差墙更稳因为删除数据会改变观测矩阵的权重平衡。批量替换可以这样写import numpy as np def clip_error(csv_path, floor0.02, cap0.08): # 假设 CSV 列x, y, freq, error_pct, imp_r, imp_i arr np.loadtxt(csv_path, delimiter,) arr[:, 3] np.clip(arr[:, 3], floor, cap) out csv_path.replace(.csv, .clip.csv) np.savetxt(out, arr, delimiter,) return outnp.clip把误差百分比限制在 2%8% 之间低于下限的数据点被抬高避免过度拟合高于上限的离群点被压低防止它主导模型。这里修的是 CSV 里的误差列如果你直接在 EM3DVP 界面里编辑记得最后导出新的观测文件否则反演还是会读旧权重。6. 把 EM3DVP 接进自动化反演流程的三个技巧6.1 用轮询日志实现收敛即停手动盯着inverse.log再决定是否停止不是不行但批量扫描多个初始模型时效率太低。可以写一个轮询脚本每 60 秒读一次日志当 model change 连续两次小于 0.03 时终止进程并保存当前代模型while true; do python check_log.py inverse.log pkill -f em3dvp inverse break sleep 60 done这里check_log.py复用上一章的函数返回True时pkill结束任务。要注意pkill不能匹配到别的项目里的同名进程给命令加上--out参数区分工程更安全。6.2 批量扫描初始模型的次序初始模型试验不要随机跑。我一般按一维分层模型、均匀半空间、测点平均电阻率三个方向各跑一组每组仅改动--init路径其余参数全部固定。这样日后的差异只能来自初始模型而不是正则化权重。对比时可以写一个 Python 脚本读取每组结果的 RMS 和历史曲线再把 EM3DVP 导出的第一层切片拼成缩略图网格用视觉判断哪一组最符合已知地质规律。6.3 用透明异常体做交付图交付图并非把每个深度的电阻率切片全部摆出来就合格。EM3DVP 通常支持按电阻率等值面生成透明体我一般取 100 Ω·m 的低阻等值面设置透明度为 0.35再叠加测点位置和地表高程。这样钻井方一眼就能看到异常体的走向和埋深不需要查看多张切片。导出时将视图旋转到视倾角 25°、方位角 135°出图时按实际尺寸标尺和北向箭头最后把地物坐标重叠到透明体边界上整个交付包就完整了。本文还有配套的精品资源点击获取