简介：EM3DVP是一套基于Matlab开发的三维地电磁建模与反演可视化工具包，主要面向地质电磁法研究人员与工程师，用于简化三维反演代码的输入模型、数据与参数文件准备，并提供结果模型和电磁响应的绘制界面。压缩包共收录243个文件，以239个m脚本文件为主体，涵盖模型创建、数据读写、响应绘制等模块，另含1个说明文档、1个许可证文件、1个Markdown说明和1张示意图，整体仅741KB，轻量易部署。目前已有231人学习浏览，适合具备一定Matlab操作经验、需要对结构化网格电磁数据进行建模和反演可视化分析的科研用户。借助该工具可快速导入地质数据、设置参数、运行反演并直观查看结果，还能通过绘图接口理解模型与数据间的对应关系，在资源勘探、工程地质调查等场景中提升三维电磁数据处理效率。
1. EM3DVP 是给三维地电磁反演省时间的，但先要读懂它的边界
EM3DVP 是给做大地电磁（MT）三维地电磁建模和反演的人省时间的：手里几十个测点的 EDI 文件，要先转成反演代码能读的格式，再逐行校对阻抗张量；网格参数改一个方向，模型文件就要全量重排。这个 Matlab 脚本包用 create_model_gui.m 和 create_resultviewer_gui.m 两个图形界面，把“EDI 读入→结构化网格建模→ModEM 数据文件输出→反演结果绘制”整条链路串起来。先说边界：它只支持结构化网格，地形起伏、非规则测网要谨慎评估。适合熟悉 Matlab、知道 ModEM 格式但没时间手搓脚本的科研人员，也能当三维反演入门教具。
2. EDI 读入与 ModEM 数据生成：read_edi.m 与 save_modemdata.m 的链路
2.1 EDI 文件里到底存了什么
EDI（Electromagnetic Data Interchange）是 MT 领域使用最广的文本交换格式，一个测点一个文件，全部是 ASCII 文本。文件按 section 组织，每个 section 以>SECTNAME开头。和三维反演前处理直接相关的有四个 section：INFO 记录测点号、采集日期、仪器类型；HMEAS 与 VMEAS 定义磁道与电道的方位角、倾角和电极距，本质上是在描述测量坐标系；ZMEAS 声明阻抗张量 ZXX、ZXY、ZYX、ZYY 的编排顺序；ZXDATA 之后才是按周期排列的阻抗实部、虚部和误差。
初次接触 EDI 的人最容易栽在 ZMEAS 上。不同采集系统对四个阻抗分量的排列不完全一致，有的按“ZXX 实、ZXX 虚、ZXY 实、ZXY 虚……”的顺序，有的则把 ZXY 与 ZYX 位置互换。EM3DVP 的 read_edi.m 实质上是把 ZXDATA 段解析成 Matlab 结构体数组，并做两个隐式约定：一是把误差值里的 0 替换成统一下限，避免单点权重无穷大；二是把周期统一换算成秒并强制升序排列。这两件事看起来琐碎，却决定了后续反演会不会在第一步就发散。
| EDI section | 主要内容 | EM3DVP 中的用途 |
|---|---|---|
| INFO | 测点 ID、日期、仪器 | 站点命名与日志 |
| HMEAS / VMEAS | 磁道/电道方位角、电极距 | 判断测点坐标系 |
| ZMEAS | 四个阻抗分量的排列顺序 | 决定数据列解析顺序 |
| ZXDATA | 周期、阻抗实虚部、误差 | 直接进入反演数据文件 |
2.2 read_edi.m 的解析逻辑与调用方式
read_edi.m 的写法不复杂，核心就是一个按行扫描的文本解析器：用 startsWith 定位 section 标记，进入 ZXDATA 后逐行 sscanf。下面是与它等价的简化骨架，保留了最关键的分支逻辑。
% 简化版 EDI 阻抗解析骨架（对应 read_edi.m 的核心路径） function d = read_edi_basic(ediFile) fid = fopen(ediFile, 'r'); lines = textscan(fid, '%s', 'Delimiter', '\n', 'Whitespace', ''); fclose(fid); lines = lines{1}; seg = find(startsWith(lines, '>ZXDATA'), 1); % 定位阻抗数据段 k = 0; for i = seg + 1 : numel(lines) s = strtrim(lines{i}); if startsWith(s, '>'), break; end % 遇到下一 section 停止 vals = sscanf(s, '%f'); if numel(vals) < 3, continue; end k = k + 1; d.period(k) = vals(1); % 单位：秒 d.Zxx(k) = vals(2) + 1i*vals(3); % 实部 + i*虚部 if numel(vals) >= 9 d.Zxy(k) = vals(4) + 1i*vals(5); d.Zyx(k) = vals(6) + 1i*vals(7); d.Zyy(k) = vals(8) + 1i*vals(9); end end end这段代码说明三件事。第一，ZXDATA 行前 9 个数是“周期、ZXX 实、ZXX 虚、ZXY 实、ZXY 虚、ZYX 实、ZYX 虚、ZYY 实、ZYY 虚”的标准顺序，实际使用前必须拿 ZMEAS 段核对一遍。第二，sscanf 返回的 vals 是按空白切分后的浮点数数组，文件里如果混入制表符或全角空格，numel(vals) 会异常，建议解析前先统一做空白字符规整。第三，read_edi.m 在 EM3DVP 里按“一个 EDI 一个结构体”的方式工作，调用上直接d = read_edi('site01.edi')，得到的是带 period、Zxx、Zxy、Zyx、Zyy、err 字段的结构体，后续绘图和保存都吃这个结构。
2.3 save_modemdata.m 的输出约定与坐标转换
ModEM 的数据文件与 EDI 完全是两套约定。文件头先声明周期个数，然后逐周期给误差模式，接着每个站点一行：站点名、东坐标、北坐标、四个分量的误差百分比，后面按“实部 虚部 实部 虚部”间隔排列阻抗值。用 save_modemdata.m 生成的片段大致是这样：
3 3.1623e-01 0 1.0000e+00 0 3.1623e+00 0 station01 523184.6 3620188.2 5 5 5 5 +2.10e+01 +3.20e-01 -4.00e+00 -1.22e+01 +3.10e+00 +8.80e-01 -2.05e+01 +2.10e+00 +1.95e+01 +7.10e-01 ... +1.80e+01 -1.20e+00 ...按行拆开看：前三行是三个周期，“0”表示该周期直接用站点行尾部的百分比误差，不做统一覆盖；第 4 行前两个数是 UTM 东、北坐标，单位米，后面四个 5 表示 ZXX、ZXY、ZYX、ZYY 的误差均为 5%。与 EDI 不同，ModEM 只认平面坐标，所以 save_modemdata.m 最关键的一步是坐标转换：把 EDI 里的经纬度按测区中央经线投影到 UTM 或本地直角坐标系。这一步出错，反演结果会和实际测点错位几公里甚至几十公里，而且从曲线上很难看出来。
注意：save_modemdata.m 内部的投影参数是按测区中央经线计算的，测网跨度超过 6° 时建议拆成两个文件分别处理。
我的习惯是保存前先打印坐标极值，确认东坐标在几十公里到几百公里的量级，而不是 118.3 这样的度数：如果看起来像经度，说明投影未生效，回头检查中央经线设置。
3. 结构化网格建模：create_model_gui 的网格设计与参数取舍
3.1 为什么 EM3DVP 只做结构化网格
结构化网格的意思是地下被离散成规则正交的长方体单元，每个 cell 有 6 个固定邻居，模型协方差算子、有限差分算子和 Jacobian 组装都能用稀疏矩阵高效实现。非结构化网格当然能更好地贴合地形和任意测点，但网格生成、插值与反演方程的装配复杂度高出一个量级，对 GUI 工具来说交互成本也大得多。EM3DVP 的定位是给 ModEM 这类结构化网格反演代码准备输入，锁定 rectilinear 网格是合理取舍。
但结构化网格有一个隐藏约束：测点必须落在网格内部，最好落在核心区而不是 padding 区。如果测网横跨数度经度，直接用经纬度建网格，单元在东西方向会被拉长，等效电阻率出现畸变。惯用做法是先把测点投影到 UTM，再在平面坐标下建网格，最后把模型写回经纬度用于成图。create_model_gui.m 内部也按这个顺序处理，用户第一步输入的就是投影后的测点坐标文件。
3.2 网格参数：空气层、趋肤深度与 padding 增长因子
create_model_gui.m 里最值得花心思的不是 GUI 布局，而是网格参数表怎么填。三维 MT 反演的网格一般分三部分：空气层（z 为负）、核心区（测点覆盖范围）、padding 区（向四周和深部扩展的过渡网格）。各参数的实际效果如下：
| 参数 | 常规取值 | 选择依据 |
|---|---|---|
| 核心区最小 cell | 测点间距的 1/2～1/3 | cell 过大会把浅部异常平滑掉 |
| 空气层层数 | 5～10 层，厚度按 ×1.2～×1.5 递增 | 顶层厚度需远超空气趋肤深度 |
| 水平 padding 层数 | 6～12 层 | 让边界反射影响降到 1% 以下 |
| padding 增长因子 | 1.3～1.6 | 更大会浪费 cell，更小则网格不够大 |
| 最大深度 | 最低频趋肤深度的 3～5 倍 | 保证深部响应完全衰减 |
趋肤深度公式是 δ = 503·sqrt(ρ·T)，单位米。假设围岩电阻率 100 Ω·m，最低周期 1000 s，δ ≈ 503×sqrt(10^5) ≈ 159 km，最大深度至少要到 500 km。网格最大深度不够的典型症状是反演在最低频段始终拟合不上去，误差棒压不下去，因为模型边界上的等效半空间已经不成立。
3.3 从 GUI 拖拽到模型文件落盘
用户在 create_model_gui.m 里输入参数后点击生成，界面内部做的工作可以概括为三步：沿三个方向生成节点向量，给每个 cell 赋初始电阻率，把网格和电阻率按 ModEM 模型文件格式写盘。深度方向的生成逻辑通常长这样：
% 构造深度方向节点：空气层取负值，地下按递增因子加密 zAir = -fliplr(cumsum([20, 20*1.5.^(1:9)])); % 10 层空气，向高空疏散 zEar = cumsum([0, 20*1.15.^(0:49)]); % 地下 50 层，首层 20 m z = [zAir, zEar]; rho = 100 * ones(nx-1, ny-1, nz-1); % 初始半空间 100 Ω·m rho(:,:, 1:10) = 1e8; % 前 10 层为空气层这里的关键是空气层阻值。不能取 0，有限差分求解器在电阻率为 0 的单元上会直接 NaN；也不建议小于 1e6，否则空气层参与电流分配，浅部视电阻率曲线会出现不该有的下降。ModEM 模型文件的写法是：第一行 Nx Ny Nz，随后是 x、y、z 三个方向的节点坐标向量，最后逐层写电阻率，楼层顺序（从空气层往下还是从最深层往上）各版本反演代码有差异，写盘前先和手册核对一版。
提示：填网格时如果测点坐标是经纬度，先投影到 UTM 再填，模型文件里的坐标单位是米，不是度。
保存模型的同时，GUI 会调用 save_data.m 把观测数据一并写出，保证网格与数据文件的测点坐标严格对应。这一步对应关系是三维反演里最多发的低级错误来源。
4. 反演响应与结果的判读路径：plot_resp、plot_sounding 与 plot_psection
4.1 plot_resp.m：先看拟合，再谈地质
反演迭代收敛不等于结果可信。create_resultviewer_gui.m 打开结果后，第一件事是用 plot_resp.m 对比观测与预测的阻抗分量，并叠加误差棒。判读顺序我习惯固定为：先看 ZXY、ZYX 两个主模式，它们在三维反演中占据主要拟合权重；再看低频段是否系统性偏移——如果是全站点的系统偏移，多半是静态位移或网格边界问题，而不是深部构造。RMS 失配的定义是 RMS = sqrt( (1/N)·Σ((obs-resp)/err)² )，EM3DVP 的结果查看界面按站点给出这个值，作为初筛，RMS 小于 2 可接受，大于 3 就要回头查数据或网格。
% 用 plot_resp 的思路核对单个测点主模式拟合 figure('Color', 'w'); loglog(d.period, abs(d.Zxy), 'o', 'MarkerSize', 5); hold on loglog(r.period, abs(r.Zxy), '-', 'LineWidth', 1.5); set(gca, 'XScale', 'log', 'YScale', 'log', 'FontSize', 11); legend('观测', '响应', 'Location', 'northwest'); xlabel('周期 (s)'); ylabel('|Zxy|');这里必须用双对数坐标：MT 阻抗在宽频带上跨越两三个数量级，线性坐标下低频段的偏差几乎看不见，双对数才能同时暴露高频浅部和低频深部的失配结构。plot_resp.m 在 EM3DVP 里还支持多站点批量平铺成小多图，适合快速扫描整条测线哪些站点有问题。
4.2 plot_sounding.m：单点视电阻率与相位曲线
plot_sounding.m 针对单个测点，把阻抗换算成视电阻率和相位两条曲线。换算关系是 ρa = |Z|²/(μω)，其中 μ = 4π×10⁻⁷ H/m，ω = 2π/T；相位是 φ = atan2(Im(Z), Re(Z))。对视电阻率曲线，我一般关心三件事：曲线整体水平对应围岩电阻率；中高频段的形态变化对应浅部结构；末端是否发散或翻转，发散往往说明该周期信噪比不够，反演时应考虑降低权重。
相位曲线的作用和视电阻率互补。视电阻率对低阻层的响应是“先下降后抬升”，相位则直接给出低频段的斜率信息。如果一个测点视电阻率曲线很光滑但相位剧烈抖动，常见原因不是地质而是仪器相位标定问题，这类点在反演前就要标记出来。plot_sounding.m 允许在曲线旁标注站点名与 RMS，导出 PNG 后可以直接贴进报告。
| 曲线特征 | 可能原因 | 处理建议 |
|---|---|---|
| 高频段整体平行上移 | 静态位移或近地表电性不均 | 相位不受影响，考虑空间滤波 |
| 低频段末端发散 | 信噪比不足 | 降低该周期权重 |
| 视电阻率光滑但相位抖动 | 仪器相位标定问题 | 反演前标记或剔除 |
4.3 plot_psection.m：伪断面不是深度剖面
伪断面是沿测线方向，横轴为测点位置、纵轴为周期的等值线图，颜色表示视电阻率或阻抗幅值。最常被误解的一点是伪断面的纵轴不是深度。周期与探测深度的关系依赖电阻率假设，均匀半空间下趋肤深度 δ = 503·sqrt(ρ·T)，同样的周期在低阻区和高阻区对应的深度相差甚远。所以伪断面适合看横向电性分带和静态位移特征，不适合直接当成剖面解释。
% 伪断面绘制骨架：横轴站点，纵轴 log10(周期)，颜色 log10(视电阻率) [X, T] = meshgrid(station_x, d.period); rhoa = abs(d.Zxy).^2 ./ (4*pi*1e-7 .* (2*pi ./ T)); % ρa = |Z|²/(μω) contourf(X, log10(T), log10(rhoa), 25, 'LineColor', 'none'); colorbar; xlabel('测点位置 (km)'); ylabel('log10(周期, s)');代码里的纵轴取 log10(周期) 是实用考量：MT 周期通常从 0.001 s 到 1000 s 跨越六个数量级，不取对数则所有低频信息挤在图底部。如果伪断面出现沿测线方向的竖直条带，而且宽度与测点间距相当，优先怀疑静态位移而非真实地质体，配合相位曲线能进一步确认。
5. 把 EM3DVP 接进反演闭环的三个实战技巧
5.1 批量读 EDI 与公共周期检查
EM3DVP 单测点用起来容易，真正的效率在批量。用 dir 拿到全部 EDI 文件名，循环调用 read_edi.m，再交给 save_data.m 合并输出。合并前先检查不同站点的周期集合是否一致：不一致时宁可对低频段做公共周期插值，也不要让数据文件出现缺测周期，否则 ModEM 会直接报错或把缺测位置当零权重计算。
files = dir('*.edi'); for i = 1:numel(files) sites(i) = read_edi(files(i).name); end save_data(sites, 'all_sites.data');5.2 save_edi.m：把反演响应导回标准格式
反演出结果后，plot_resp.m 看到的是内部结构体；如果需要和商业软件或他人代码对比，用 save_edi.m 把预测响应写成标准 EDI 文件，直接放进任何支持 EDI 的成图工具。这个文件与观测 EDI 结构完全一样，只是阻抗值换成了预测值，对比时把两者叠绘即可量化拟合。
5.3 三个高频踩坑点
实验里最常遇到的三类问题供排查：
- 坐标投影不一致导致测点与网格错位，表现为所有测点 RMS 同时偏高；
- 空气层电阻率设置过小（低于 1e6 Ω·m），浅部曲线异常下掉；
- 误差下限未设置，个别周期误差为 0 导致权重无穷大，反演迭代步长震荡。
排查顺序建议从第三条开始，因为它是纯数据问题，修复成本最低；再查投影；最后怀疑网格设置。
本文还有配套的精品资源，点击获取
