简介面向地球物理勘探与科研人员这份MATLAB压缩包聚焦一维正演模拟覆盖瞬变电磁法TEM、CSAMT法及电测深三种常用非侵入性探测技术。资源针对地下电性结构建模需求帮助使用者理解电磁波在地层中的传播规律掌握由地电模型正演计算地表响应的方法。包内共31个文件核心为15个m脚本如temfwd.m、FJCST.m等承担正演计算、滤波定义与绘图任务配套11个mat数据文件保存模型与中间结果、4个fig图形文件展示曲线界面外加1个dat数据文件辅助验证整体仅121KB结构紧凑易于研读。已有1185人学习下载受到相关领域学习者的关注。通过运行这些代码读者可以直观对比TEM、CSAMT与电测深的正演响应了解不同方法在浅层与深部探测中的差异并可基于现有框架修改地电参数或扩展反演算法是一份兼顾教学与自学的实用参考资料。1. 一维正演TEM、CSAMT、电测深共用的数学地基拿到这套物探一维正演包时第一反应是看文件名构成的边界。temfwd.m、FDEMfwd.m、FJCST.m、dcgsafwd.m、Z_CR.m刚好覆盖瞬变电磁TEM、可控源音频大地电磁CSAMT和直流电测深三条主线而model_result系列mat文件则是每一条线跑完后的产物。一维正演本身计算量不大难的是把“地层电阻率加厚度”翻译成地表可观测的电磁响应这中间涉及分层递推、汉克尔变换、频时转换与视电阻率定义代码里任何一处索引错位都会让曲线形状完全失真。对搞TEM和CSAMT资料解释的人先把这组matlab函数跑通才能谈到反演初始模型和数据质量评估对做工程电法的同行电测深曲线的一维正演又是解释工作的底线。下面按函数在执行路径中的位置拆开讲并给出一套可复现的最小测试流程。2. TEM瞬变电磁正演拆解temfwd.m与FDEMfwd.m的时域变换路径2.1 为什么TEM正演要先从频域走一趟瞬变电磁的一维正演有两条路线。一条是直接用时域解析式把阶跃关断后的二次场写成误差函数与指数核的组合这种方法只对均匀半空间和少数简单层状模型成立另一条是先计算频率域的表面阻抗或垂向磁场再利用汉克尔变换和余弦变换转到时域这套包采用的是后一条路。原因是分层模型中的反射系数递推在频域是显式的每一层只需要做一次向上传播的阻抗更新代码结构非常规整而时域解析式每增加一层就要重推一次公式。频时转换的关键工具是数字滤波。filters.mat里存放的是一组预先计算好的滤波系数包括零阶贝塞尔函数J0的汉克尔滤波系数、一阶J1以及余弦变换系数。数字滤波的基本思路是把积分核与权重系数在对数坐标系上做褶积用很少的采样点逼近整个积分数值复杂度从逐点积分的O(N×M)降到一次褶积的O(NM)。这种做法在地球物理电磁正演里几乎是标准操作区别只在滤波器长度和采样密度。这套包里filters.mat就是为temfwd.m、dcgsafwd.m、FJCST.m等函数共用的因此运行前务必保证工作目录包含该文件否则会直接报出无法加载变量的错误。文件名称内容类型使用位置filters.matJ0、J1、余弦变换滤波系数temfwd.m、dcgsafwd.m、FJCST.mgy.mat、gy1.mat工程模型的中间结果模型对比与曲线验证zeros_J.mat贝塞尔函数零点或滤波零点表汉克尔变换辅助数据model_result系列.mat分层模型正演结果后处理、绘图与反演对照2.2 temfwd.m的参数约定与最小调用实例temfwd.m是瞬变电磁时域正演的主入口其输入输出约定直接决定了其他函数怎么调用它。常见做法是接收三个参数地电模型、观测系统参数和时间轴采样。地电模型用两行矩阵表示第一行是各层电阻率第二行是各层厚度最后一层厚度写0代表无限半空间。观测系统参数则包括发射线圈半径、收发距和接收线圈有效面积。function [tline, V, rho_a] temfwd(rho, h, cfg) % rho : 1×n 各层电阻率(Ω·m) % h : 1×(n-1) 各层厚度(m)末层视为半空间 % cfg : [a, r, A, t_min, t_max, nt] % a 发射半径(m), r 收发距(m), A 接收面积(m²) % t_min/t_max 时间窗(秒), nt 采样点数 % 返回 : 时间轴、感应电压(V)、晚期视电阻率(Ω·m) % S load(filters.mat, J0, COS); freq logspace(-3, 7, 512); % 频域采样 1e-3~1e7 Hz Hf FDEMfwd(rho, h, cfg(1:2), freq); % 频域垂向磁场 Hhz hankel_J0(Hf, cfg(2), S.J0); % 汉克尔变换到收发距上 V cosine_transform(Hhz, freq, tline, S.COS); V -cfg(3) .* V; % 法拉第感应定律 rho_a late_time_appres(tline, V, cfg(1)); % 晚期视电阻率 end代码第一行到第五行在组装模型参数第六行从filters.mat读取滤波系数第七到第九行完成频率采样与频域响应计算第十行以后的汉克尔变换和余弦变换是关键。需要注意一点FDEMfwd.m只负责算频域的H_z响应它自己不做任何时域转换频率点取多少、取到多高直接决定后期时间道是否稳定一般至少需要覆盖发射电流基频的100倍。调用时最实用的是直接读取包里已有的模型结果文件做参数对齐load(model_result111.mat, rho, h); cfg [50, 0, 100, 1e-5, 1, 400]; % 半径50m, 接收面积100m² tline logspace(-5, 0, 400); [t, V, rho_a] temfwd(rho, h, cfg); loglog(t, abs(V)*1e6, LineWidth, 1.2); xlabel(t / s); ylabel(感应电压 / μV); grid on;这里用loglog绘图是因为TEM衰减曲线跨越多个数量级若用线性坐标会在早期道挤成一团。abs取模的原因是感应电压在早期可能为负值直接从符号上无法反映衰减趋势。2.3 FDEMfwd.m与temfwd.m的协作边界FDEMfwd.m在包里的定位是频域通用底层函数它根据输入的电阻率、厚度和频率向量返回分层介质的表面电磁场分量。temfwd.m、FJCST.m都建立在它之上区别只在于后续是走余弦变换还是直接输出阻抗。这种分层的工程设计有一个明显好处你在检查代码时不需要同时关注频域求解和时域转换两套逻辑出问题可以先单独用均匀半空间模型测试FDEMfwd.m的频响是否匹配解析解再回头查滤波变换部分。function Hz FDEMfwd(rho, h, geo, freq) % rho: 各层电阻率, h: 各层厚度 % geo: [a, r] 发射半径与收发距 % freq: 频率数组(Hz) % 返回: 垂向磁场 Hz复数 % w 2*pi*freq; for k 1:numel(freq) u sqrt(1i*w(k)*4e-7*pi./rho); % 各层波数 Z z_te_recursive(u, rho, h); % 由底向上递推 Hz(k) integral_kernel(Z, geo(1), geo(2)); end end代码里z_te_recursive做的是TE模分层递推从最后一层半空间开始逐步向上累积表面阻抗。这里有一个常见误区很多初学者会把波数写成实数的衰减系数实际上在谐变时间因子作用下波数是复数虚部对应感应电流的相位滞后这正是TEM二次场能包含地层电导率信息的物理来源。FDEMfwd.m不关心时间也不做符号判断把频响交给上游函数后由上游决定用实部还是虚部参与后续变换。3. CSAMT与电测深正演复阻抗计算、频率扫描与直流滤波3.1 Z_CR.m如何计算CSAMT复阻抗CSAMT正演的核心不是视电阻率公式本身而是地表复阻抗ZE_x/H_y的求取。对层状介质TE模与TM模阻抗各自独立递推野外实际观测通常是两者的某种混合但在一维正演中多数程序默认只计算TE模因为TM模受静态位移影响大曲线形态也更复杂。Z_CR.m从命名看就是复阻抗计算的独立模块接受层参数和频率返回复数阻抗。function Z Z_CR(rho, h, f) % 计算层状介质地表复阻抗(Ω) % rho: 层电阻率(Ω·m), h: 层厚度(m), f: 频率(Hz) % w 2*pi*f; Z zeros(size(f)); for k 1:numel(f) Z(k) te_impedance(rho, h, w(k)); end end由复阻抗得到卡尼亚视电阻率ρ_a |Z|²/(ωμ₀)相位φatan(Im(Z)/Re(Z))。注意单位换算MATLAB里角频率w默认是rad/s卡尼亚公式里的μ₀取4π×10⁻⁷ H/m输入频率为Hz时不需要再除以2π。很多CSAMT程序在低频段出现视电阻率负值多半不是算法问题而是递推时把TM模或TE模的平方根选错了支根导致相位符号翻转。3.2 FJCST.m与FCST.m的频率扫描实现FJCST.m和FCST.m是一对同源函数差异大概率在源类型或极化方向处理上。常见做法是对频率数组做循环逐个频点调用Z_CR.m获取阻抗再换算成视电阻率和相位最终形成频率测深曲线。这类函数的关键参数是频点数量和对数间隔。CSAMT的观测频段一般在0.1Hz到10000Hz之间每十倍频程取6到12个频点就够一维正演使用频率过密只增加计算时间曲线细节并不会因此改善。function [f, rho_a, phi] FJCST(rho, h, freq, AB) % 输入: % rho 层电阻率向量, h 层厚度向量 % freq 频率序列(Hz), AB 发射极距(m) % 输出: % f 频率轴, rho_a 视电阻率, phi 相位(度) % Z Z_CR(rho, h, freq); rho_a abs(Z).^2 ./ (2*pi*freq*4e-7*pi); phi atan2(imag(Z), real(Z)) * 180/pi; end这个函数的调用很简单但使用时必须保持rho与h的维度一致。如果从model_result.mat加载的模型里厚度少了一层递推到底层半空间时会越界MATLAB不会报维度错误而是返回NaN曲线直接断掉。排错时先把Z_CR单独拿出来跑一次看低频端是否收敛到一个有限值高频端是否符合表层电阻率的渐近值。3.3 dcgsafwd.m与直流电测深的线性滤波正演直流电测深的正演与TEM、CSAMT有一个本质区别它不涉及感应只考虑传导电流数学上简化为对核函数做一阶贝塞尔函数积分。dcgsafwd.m的典型处理方式是用反射系数递推得到核函数K(λ)然后对K做J1汉克尔滤波乘上装置系数即得视电阻率。function rho_a dcgsafwd(rho, h, AB2, MN2) % 对称四极装置直流电测深正演 % rho: 层电阻率, h: 层厚度 % AB2: 供电极半距数组(m), MN2: 测量极半距数组(m) % S load(filters.mat, J1); K kernel_from_layers(rho, h); % 核函数λ域采样 rho_a j1_filter(K, AB2, MN2, S.J1); % 线性滤波求出视电阻率 rho_a rho_a .* (AB2.^2 ./ MN2); % 装置系数换算 end直流电测深最需要注意的是电极距极差。AB2从几米到几百米变化MN2如果一直保持固定小值深部曲线会被浅层信号淹没程序里MN2一般按AB2的比例同步放大即保持MN/AB比值在0.1到0.3之间。包里的compute.fig和cr1dmod.fig大概是这类装置系数换算的图形界面用鼠标调整AB2数组后直接看曲线变化比命令行里改参数要直观。4. 从model_result到成品图件compute.m、batch.m、plotdat.m的数据串联4.1 model_result系列mat文件里到底存了什么打开包目录能看到的model_result.mat、model_result111.mat、model_result1122.mat、model_result11.mat、model_result1111.mat、model_resul11212t.mat后缀数字的差异多半代表不同层数或不同装置下的正演结果。这些文件通常保存两到三个变量最常见的是rho、h和响应矩阵。响应矩阵的行对应时间道或频点列对应不同的收发距或电极距读取后可以直接绘制测深曲线或视电阻率断面。文件名称可能对应的正演内容典型用途model_result.mat默认模型TEM响应主流程演示model_result11.mat两层模型或小极距响应曲线对比model_result111.mat三层地电模型响应层参数验证model_result1111.mat四层模型或加密频点响应深部结构测试model_result1122.mat变观测系统参数后的响应系统参数敏感性分析model_resul11212t.mat更新模型后的临时保存版本迭代中间过程读取这些mat文件时我用了一个小技巧先whos -file model_result1122.mat列出变量名而不是直接load进工作区因为历史遗留的mat文件里经常混进plot用的临时变量直接load会把当前工作空间搞乱。确认变量名后再按需load(model_result111.mat, rho, h)既安全又省内存。4.2 compute.m与batch.m的批量执行流程compute.m处理单次正演batch.m做批量循环。batch.m的存在说明这套代码早就是按项目批量处理思路设计的一次性把多个模型结果算好存成model_result系列后续画图直接从mat文件读取不需要重算。下面给出一段与batch.m风格一致的批量调用代码% 批量正演不同层参数模型并保存结果 cfg [50, 0, 100, 1e-5, 1, 400]; list {[100 200], [100 20 300], [10 100 500 50]}; % 三种模型 for k 1:numel(list) m list{k}; rho m(1:2:end); h m(2:2:end); [t, V] temfwd(rho, h(end), cfg); save(sprintf(model_%d.mat, k), t, V, rho, h); end这里h(end)用得比较巧妙实际上是只取了最后一个厚度值因为temfwd的h参数约定最后一层是半空间不需要传入最后一层厚度。批量脚本最怕参数硬编码把发射半径、时间窗都写死在循环里一旦观测系统改变要改多处。建议把cfg定义提到脚本顶部用结构体存起来不同装置方案用不同cfg结构体传入比散装变量容易维护。4.3 绘图函数与伪断面图的生成方式plotdat.m负责绘制单条测深曲线topview.m负责把多条曲线拼成视电阻率伪断面图。topview从文件名就能看出是把测点沿剖面排列纵轴用时间或频率横轴用测点号数值用色标表达形成俯视图。伪断面图是野外物探资料解释中最常用的展示形式它能从全局角度看到地层电性横向变化单条曲线看不出这种空间连续性。% 使用topview.m的典型流程 x 50:50:500; % 测点位置 rho_all zeros(numel(t), numel(x)); for k 1:numel(x) rho_all(:,k) temfwd(rho, h, cfg); end topview(x, t, rho_all, log);topview.m内部一般用pcolor或contourf绘制色块图纵轴取对数时间轴。MATLAB的pcolor行列逻辑容易混需要确认第一维是时间道还是测点号方向搞反了断面就会横竖颠倒。绘图函数里加载filters.mat的情况少更常见的是直接读取算好的model_result系列数据因此正演脚本与绘图脚本通常分开运行这也是batch.fig与topview.fig同时存在的意义——前者管计算界面后者管出图。5. 一维正演的验收手段与参数调整速查5.1 用均匀半空间解析解验证代码正确性任何正演代码拿到手先用均匀半空间模型跑一条理论曲线校验不要直接套多层模型。TEM中心回线晚期视电阻率满足渐近公式σ μ₀·M/(20π^{3/2}·t^{5/2}·V_t)其中M为发射磁矩V_t为感应电动势。把包内代码算出的晚期段视电阻率与设定电阻率对比误差在5%以内说明频域递推和时域变换链路没有方向性问题。CSAMT则用远区视电阻率平坦段做校验低频渐近值应趋近最深层的电阻率。% 均匀半空间校验脚本 rho0 100; h0 1000; % 半空间近似 cfg [50, 0, 100, 1e-5, 1, 400]; [t, V] temfwd(rho0, h0, cfg); rho_a 1e7 * pi^(3/2) * t.^(5/2) .* abs(V) ... / (4e-7 * pi * 50^2); loglog(t, rho_a); yline(100, --, 理论值); ylim([1 1000]);这里把公式稍作变形反算出视电阻率便于直接与理论值对比。如果曲线在中晚期段出现明显上翘或下坠优先怀疑filters.mat里的余弦滤波系数长度是否覆盖了当前时间窗而不是怀疑地层递推代码。5.2 常见运行问题与参数修正现象可能原因处理方式load(filters.mat)报错工作目录未切换到代码目录用绝对路径或cd到工程目录结果全为NaNrho与h维度不一致检查厚度矩阵是否比电阻率少一行低频段视电阻率成负值波数平方根支根选错检查递推中复数的sqrt符号约定曲线早期段震荡剧烈时间采样点数过少或频点覆盖不足增加nt到400以上频率上限提到1e7伪断面图横轴与纵轴颠倒pcolor行列参数顺序写反用whos确认响应矩阵维度后再绘图Excel打开mat文件乱码mat格式是二进制专有格式用h5read或MATLAB导出为csv后再查看如果是老版本MATLAB升级到新版本后再跑这套代码优先排查filters.mat的加载方式旧代码常用load filters.mat不带变量名新版本对文件格式更严格。算完一条曲线后记得用plotdat.m把正演响应与野外实测曲线放在同一坐标系对比只靠数值检查发现不了装置系数单位错误目视对比虽然土但最有效。本文还有配套的精品资源点击获取