简介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 文件里到底存了什么EDIElectromagnetic Data Interchange是 MT 领域使用最广的文本交换格式一个测点一个文件全部是 ASCII 文本。文件按 section 组织每个 section 以SECTNAME开头。和三维反演前处理直接相关的有四个 sectionINFO 记录测点号、采集日期、仪器类型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.0000e00 0 3.1623e00 0 station01 523184.6 3620188.2 5 5 5 5 2.10e01 3.20e-01 -4.00e00 -1.22e01 3.10e00 8.80e-01 -2.05e01 2.10e00 1.95e01 7.10e-01 ... 1.80e01 -1.20e00 ...按行拆开看前三行是三个周期“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/21/3cell 过大会把浅部异常平滑掉空气层层数510 层厚度按 ×1.2×1.5 递增顶层厚度需远超空气趋肤深度水平 padding 层数612 层让边界反射影响降到 1% 以下padding 增长因子1.31.6更大会浪费 cell更小则网格不够大最大深度最低频趋肤深度的 35 倍保证深部响应完全衰减趋肤深度公式是 δ 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_psection4.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 导致权重无穷大反演迭代步长震荡。排查顺序建议从第三条开始因为它是纯数据问题修复成本最低再查投影最后怀疑网格设置。本文还有配套的精品资源点击获取
EM3DVP：从EDI到ModEM的大地电磁三维反演前处理工具
简介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 文件里到底存了什么EDIElectromagnetic Data Interchange是 MT 领域使用最广的文本交换格式一个测点一个文件全部是 ASCII 文本。文件按 section 组织每个 section 以SECTNAME开头。和三维反演前处理直接相关的有四个 sectionINFO 记录测点号、采集日期、仪器类型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.0000e00 0 3.1623e00 0 station01 523184.6 3620188.2 5 5 5 5 2.10e01 3.20e-01 -4.00e00 -1.22e01 3.10e00 8.80e-01 -2.05e01 2.10e00 1.95e01 7.10e-01 ... 1.80e01 -1.20e00 ...按行拆开看前三行是三个周期“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/21/3cell 过大会把浅部异常平滑掉空气层层数510 层厚度按 ×1.2×1.5 递增顶层厚度需远超空气趋肤深度水平 padding 层数612 层让边界反射影响降到 1% 以下padding 增长因子1.31.6更大会浪费 cell更小则网格不够大最大深度最低频趋肤深度的 35 倍保证深部响应完全衰减趋肤深度公式是 δ 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_psection4.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 导致权重无穷大反演迭代步长震荡。排查顺序建议从第三条开始因为它是纯数据问题修复成本最低再查投影最后怀疑网格设置。本文还有配套的精品资源点击获取