简介：这是一份面向地球物理勘探与科研人员的直流电测深（DCR）正演模拟工具包，解决在地下地质模型未知情况下快速计算视电阻率响应的问题。压缩包共4个文件，均为MATLAB脚本（.m），涵盖核心正演计算函数、雅可比矩阵实现以及两个示例脚本，可帮助用户理解正演模型构建、参数化输入与结果输出的完整流程。资源整体仅2KB，轻量易用，适合地球物理、测深技术初学者入门，也可供有经验的研究者快速验证不同地层参数下的视电阻率变化。已有230人学习下载。通过这套程序，用户能掌握直流电场在地下传播的数学模型，并基于示例脚本修改地层电阻率、厚度等参数开展模拟实验，为后续反演算法研究提供理论与计算基础。
1. 从一条视电阻率曲线说起
直流电测深（DCR）正演，简单说就是：给定地层电阻率分层模型，算出地表电极排列下应该测到的视电阻率曲线。sdc1dfwd.rar 里就是一套 MATLAB 实现的一维直流电测深正演代码——四个 .m 文件，主函数、雅可比、两个示例。很多物探同行拿到这套程序，想跑通又不清楚滤波系数、采样间隔怎么调，结果卡在数值震荡上。这套代码的适用边界很明确：一维水平层状介质、点电源近似、对称四极或温纳装置；它解决的是“给定模型预测响应”的问题，也是后续反演和资料解释的公共底座。适合刚接手正演任务的研究生，也适合想把这套代码改造成反演前向算子的工程师。
2. 直流电测深正演的理论骨架与控制方程
2.1 稳定电流场下的电位控制方程
直流电法在地下建立的电场是稳定场，忽略感应效应。在点电源情形下，电位 φ 满足泊松方程，形式如 ∇·(σ∇φ) = -Iδ(r-r_s)。对于一维水平分层的地电模型，电阻率只随深度变化，即 ρ=ρ(z)，那么在柱坐标系下可以分离变量，把三维问题变成对径向距离 r 和深度 z 的二维问题。进一步做汉克尔变换，电位可以表示为贝塞尔函数的积分：φ(r,z) = (1/2π) ∫₀^∞ K(k,z) J₀(kr) dk。这里 K(k,z) 是核函数，由各层电阻率和厚度递推得到；J₀ 是第一类零阶贝塞尔函数。实际计算中不可能解析求这个积分，而是采用线性数字滤波法做数值逼近。
2.2 线性数字滤波法与核函数递推
线性数字滤波法的思想，是把汉克尔积分近似成离散卷积。常见做法是构造一组滤波系数 a_i 和基准距离 r_0，使积分结果表示为：φ(r) = (1/r) Σ_{i=1}^{N} a_i · T(x_i)，其中 x_i 是与 r 和滤波系数的对数采样间隔有关的坐标。滤波系数的点数 N 通常取 40～200，采样间隔一般设为 ln(10)/10，即每个十倍距离取 10 个点。N 越大，覆盖的波数范围越宽，但计算量也线性上升。对直流测深来说，极距从几米到几百米跨越三个量级，所以滤波系数至少要能覆盖 10⁻³～10 rad/m 的波数段。
核函数 K(k,z) 的计算则是标准的反射系数递推。从底层半空间开始，每一层界面的反射系数由相邻层电阻率和厚度决定。递推的起点是从半空间向上逐层回溯，最后得到地表处的核函数。在 sdc1dfwd.m 中，这个递推被封装成局部函数，输入是 rho 和 h，输出是对应离散波数 λ 的核函数值。这里的波数采样必须与滤波系数的对数网格严格对齐，如果错位，曲线的尾支就会失去“渐近到深层电阻率”的物理特征。
2.3 视电阻率的定义与装置的几何因子
正演最终输出的是视电阻率曲线，而不是电位本身。视电阻率由实测参数换算：ρa = K_g · ΔU / I。ΔU 是测量电极 M、N 间的电位差，I 是供电电流，K_g 是装置系数，只与电极几何位置有关。不同电极装置有不同 K_g 表达式，下表列出最常见的三种：
| 装置类型 | 电极排列示意 | 装置系数 K_g |
|---|---|---|
| 温纳 (Wenner) | A-M-N-B 等间距 a | 2πa |
| 施伦贝谢 (Schlumberger) | A-M-N-B，AB/2≫MN/2 | π·(AM·AN)/MN |
| 偶极-偶极 (Dipole-dipole) | A-B-M-N 等间距 | π·n(n+1)(n+2)a |
示例脚本中通常采用温纳或对称四极装置，具体要看 sexample 里的装置系数赋值。如果程序内部用温纳，K_g=2πa，那么 AB2 向量对应的 a 就是 AB2/2。改装置类型时，K_g 和电位差的取点方式必须同步修改，否则曲线畸变。野外仪器输出的是模拟信号，需要经译码器转为数字信号后才能记录存储；正演代码里直接算连续量，但这不意味着离散化不重要——极距取点方式就是等效的“采样率”。
2.4 模拟信号与数字信号的衔接考量
正演代码里算的是连续数学量，但野外仪器接收的是模拟信号电压，经放大器后要由 ADC 做数字采样。采样率不够、没有抗混叠滤波时，数字信号会失真。这个问题在正演中对应的则是对极距的离散采样——AB2 取点太疏会漏掉浅层响应，太密又会放大数值误差。常见做法是对数等间距取点，每十倍距 6～10 个点，保证曲线形态完整。这就像译码器把模拟量转成数字量时的量化步长一样，步长选错了，后续曲线解释都是白搭。
3. sdc1dfwd.rar 的文件结构与核心函数实现
3.1 文件职责与调用链
压缩包解开后是四个 .m 文件，建议先按表格建立全局认识：
| 文件 | 类型 | 职责 |
|---|---|---|
| sdc1dfwd.m | 主函数 | 输入层状参数和极距，输出视电阻率 |
| sdc1dJacob.m | 雅可比函数 | 计算视电阻率对各层电阻率的偏导数，供反演使用 |
| sexample1.m | 脚本 | 单模型正演并绘图，展示最小调用方式 |
| sexample2.m | 脚本 | 多模型对比正演，展示参数扫描效果 |
调用关系是 sexample -> sdc1dfwd；反演循环里 sdc1dJacob -> sdc1dfwd。sdc1dJacob 是站在正演之上的检测工具，不参与正演本身的计算链路。
3.2 sdc1dfwd.m 主函数的关键实现
主函数的典型接口定义如下：
function [rhoa, AB2_out, MN2_out] = sdc1dfwd(rho, h, AB2, MN2) % rho : 各层电阻率，尺寸 1×N，N=层数 % h : 各层厚度，尺寸 1×(N-1)，最后一层是半空间不需要厚度 % AB2 : 供电电极极距的一半 AB/2，单位 m % MN2 : 测量电极极距的一半 MN/2，单位 m % rhoa : 与 AB2 同样长度的视电阻率向量函数内部大体分为三段。第一段根据 rho 和 h 递推核函数；第二段对每个极距做汉克尔滤波得到地表电位；第三段由 M、N 两点的电位差和装置系数换算视电阻率。汉克尔滤波部分通常写成一个小函数，把滤波系数表预存成常量以避免重复计算。如果用扰动法求雅可比，sdc1dfwd 会被调用几十上百次，所以这段预存的效率直接决定反演速度。
下面给出一个典型汉克尔滤波的 MATLAB 风格实现片段（该片段逻辑是行业通用做法，与压缩包内实现等价）：
function phi = hankel_filter(K, r, filter_coeff) % K : 核函数，在离散波数 lambda 上取值 % r : 目标半径 % filter_coeff : 包含 b 和 c 两个向量，b 是滤波系数，c 是系数间隔 b = filter_coeff.b; % 滤波系数 lambda0 = filter_coeff.lambda0; % 基准波数 dlogr = filter_coeff.dlogr; % 对数间隔，常为 ln(10)/10 phi = 0; for i = 1:length(b) lam = (1/r) * exp((b(i)) * dlogr); % 通过插值获取核函数在该波数的值 K_lam = interp1(filter_coeff.lambda, K, lam, 'spline', 0); phi = phi + b(i) * K_lam; end phi = phi / r; % 汉克尔变换结果 end上述代码中，b(i) 是滤波权重，它的构造来自数字化汉克尔变换的经典系数；lambda0 和 dlogr 决定了波数采样网格。interp1 的插值要选中 spline，因为核函数在波数域是光滑变化的；最后一个参数 0 表示波数超出范围时返回零，防止外推引入幻影值。实际程序中这一步也可以不用插值，而是直接把核函数和滤波系数对齐，这样不涉及插值误差，但要求 lambda 采样与滤波系数严格一致。
3.3 sdc1dJacob.m 的雅可比实现
反演需要用到的雅可比矩阵，在示例代码里是以扰动法实现的。实现思路是“差商近似偏导”，代码如下：
function J = sdc1dJacob(rho, h, AB2, MN2) % 基础正演响应 rhoa0 = sdc1dfwd(rho, h, AB2, MN2); n = length(rho); m = length(AB2); J = zeros(m, n); delta = 0.02; % 扰动比例，经验值 for j = 1:n rho_pert = rho; rho_pert(j) = rho(j) * (1 + delta); rhoa_pert = sdc1dfwd(rho_pert, h, AB2, MN2); J(:,j) = (rhoa_pert - rhoa0) / (rho(j) * delta); end end扰动比例取 0.02 是权衡数值噪声和线性近似误差的结果。比例太小，正演本身截断误差会淹没差商；比例太大，又超出线性区间。更稳妥的自适应选择是设置 eps_j = max(rho(j)*0.01, 1e-6)。此外，如果你的电阻率跨多个数量级（比如 1 到 100000 Ω·m），建议改到对数模型下求导，即直接对 ln ρ 扰动。这样可以保证每个分量的相对扰动一致，避免深部低阻层被完全忽略。
3.4 示例脚本的差异化设计
sexample1.m 通常展示一个两层或三层模型下视电阻率曲线的“S”型形态——浅层高阻、深层低阻时，曲线从左到右先平、后降、再平。sexample2.m 则可能循环改变某一层厚度或电阻率，画出多条曲线对比，用来直观理解参数灵敏度。这两个脚本的共同点是都用 loglog 绘图，因为测深曲线在双对数坐标下才有几何相似性。建议把脚本里的默认参数地图化，改成你工区的背景值，再判断正演响应是否符合地质认识。
4. 实战：跑通正演并自定义地层模型
4.1 运行前的环境准备
这套代码只要是 MATLAB R2010 以后的版本都能直接跑，不需要额外工具箱。把解压后的 sdc1dfwd 文件夹加入 MATLAB 路径：
addpath('D:\Downloads\sdc1dfwd'); % 或者用 cd 进入目录 cd('D:\Downloads\sdc1dfwd');在命令行输入sexample1并按回车。如果顺利，会弹出一个双对数坐标的 Figure，横轴是 AB/2，纵轴是视电阻率，曲线形态应该平滑单调。如果报错，先检查工作路径是否含中文或空格，老版本 MATLAB 对路径里的特殊字符很敏感。
4.2 自定义三层模型的标准改法
我来演示一个河流阶梯沉积地层的例子：表层粉质黏土 5 m，电阻率 80 Ω·m；下卧砂层 20 m，电阻率 200 Ω·m；再下面是泥岩基岩，电阻率 40 Ω·m。电极距从 1 m 到 300 m，每十倍距取 8 个点，MN/2 取 AB/2 的 1/10。脚本如下：
% 自定义模型：三层 rho = [80, 200, 40]; h = [5, 20]; % 对数等间距极距，范围 1~300 m AB2 = logspace(0, log10(300), 25); MN2 = AB2 / 10; % 正演计算 [rhoa, AB2_out, MN2_out] = sdc1dfwd(rho, h, AB2, MN2); % 绘制 figure('Color','w'); loglog(AB2_out, rhoa, 'b-o', 'LineWidth', 1.5, 'MarkerSize', 4); xlabel('AB/2 (m)', 'FontSize', 12); ylabel('视电阻率 \rho_a (\Omega\cdotm)', 'FontSize', 12); title('直流电测深正演曲线', 'FontSize', 12); grid on;这段代码的命令解释：rho 向量长度比 h 多 1，这是硬性规定；AB2 用 logspace 生成对数均匀点，保证短极距和长极距在曲线上的采样密度一致；MN2 = AB2/10 是经验比例，实际装置里如果 MN 过大，电位差近似会失效，导致曲线尾支上翘。
参数表：
| 参数 | 含义 | 单位 | 本例取值 | 备注 |
|---|---|---|---|---|
| rho(1) | 第一层电阻率 | Ω·m | 80 | 表层覆盖 |
| rho(2) | 第二层电阻率 | Ω·m | 200 | 含水砂层 |
| rho(3) | 第三层电阻率 | Ω·m | 40 | 基岩半空间 |
| h(1) | 第一层厚度 | m | 5 | 底深 5m |
| h(2) | 第二层厚度 | m | 20 | 底深 25m |
| AB2 | 供电极距一半 | m | 1~300 | 对数均匀 |
| MN2 | 测量极距一半 | m | 0.1~30 | AB2/10 |
4.3 半空间验证与曲线趋势判据
拿到结果先别急着用，用半空间模型校验程序完整性是最省事的做法。半空间就是只有一层，没有厚度：
rho = [250]; h = []; [rhoa_test] = sdc1dfwd(rho, h, logspace(-1, 2, 10), logspace(-2, 1, 10)); disp(max(abs(rhoa_test - 250)));如果输出小于 1e-6，说明主函数里对空厚度向量有正确分支。如果输出接近 0.1 量级，说明滤波系数在极距和波数匹配上还有问题。另一个物理判据是：在最浅极距下，ρa 应接近第一层电阻率；在最大极距下，ρa 应趋近最深层电阻率。对前面的三层模型，曲线左端应接近 80，右端应接近 40。
4.4 直流电与交流电的类比：为什么不用感应理
直流电测深使用的是稳定电流，它不产生显著的电磁感应，因此正演只需解拉普拉斯方程，而无需像频域电磁法那样考虑趋肤深度。这就像安规电容用在直流电路中时，主要考核的是直流耐压和漏电流，而不是像交流应用那样关注容抗损耗。理解了直流和交流在物理机制上的差异，你就不会拿交流电磁法的反演思路来强套这套直流代码。
4.5 失败时看什么
最常出的问题是极距数组不单调或长度不匹配。sdc1dfwd 内部如果假设 AB2 是递增向量，你传一个乱序的 AB2 会导致电位差计算错乱。另一个问题是 h 的长度不等于 length(rho)-1，直接抛错。所以在调用前加一行断言：
assert(length(rho) == length(h)+1, '层数和厚度不匹配'); assert(all(diff(AB2) > 0), 'AB2 必须严格递增');这就是把丑陋的错误提前暴露出来，而不是让数值运算在静默中给你一条看似合理实则错误的曲线。
5. 进阶：雅可比矩阵在反演衔接中的正确打开方式
5.1 反演方程组与雅可比的角色
一次性正演只是起点，真正价值在反演里体现。阻尼最小二乘迭代格式为 Δm=(JᵀJ+βLᵀL)⁻¹Jᵀ(d-F(m))。J 由 sdc1dJacob 产生，它的每一行对应一个观测极距，每一列对应一个层参数。在写反演循环前，先检查 size(J) 是否等于 [length(AB2), length(rho)]，维度错了后面矩阵乘法全是 NaN。
5.2 对数参数化改进数值稳定性
对电阻率直接求导时，若各层电阻率差两三个数量级，雅可比矩阵条件数会非常差。我常用的做法是把模型参数换成 m=lnρ，对雅可比做一次链式变换：J_log(:,j) = rho(j) * J(:,j)。这样每个分量的尺度一致，正则化矩阵也不用专门地按层缩放。实现时只需在扰动法基础上乘上 rho(j)，代价几乎为零，但收敛稳定性改善明显。
5.3 常见数值问题排查表
| 症状 | 病因 | 对策 |
|---|---|---|
| 曲线高频振荡 | 滤波系数点数不足 | 增加滤波系数 N 到 100+ |
| 长极距段 ρa 发散 | 波数下限太高 | lambda_min 降到 1/(100*maxAB) |
| 半空间验证误差大 | 滤波系数与极距尺度不匹配 | 检查滤波系数基准距离是否与极距单位一致 |
| 雅可比矩阵 NaN | 扰动后电阻率变负 | 改用乘法扰动或对数扰动 |
| 反演不收敛 | J 和 F(m) 维度不一致 | 统一通过 sdc1dfwd 输出 |
5.4 用 persistent 缓存加速反演循环
sdc1dfwd 每次调用都要重建滤波系数表，而反演中正演可能被调用上千次。常见的优化是使用 persistent 变量缓存：
function phi = hankel_filter_cached(K, r) persistent coeff; if isempty(coeff) coeff = load_filter_coefficients(); % 只加载一次 end % 后续使用 coeff end这个改动简单但效果明显。把滤波系数从每次计算转移到首次调用加载，1000 次正演计算大概能省去 60% 的冗余工作。如果你的数据加了噪声、要做 20 次迭代的随机反演，这个加速会直接决定你在下班前能不能跑完。
这一步完成后，雅可比矩阵就能直接喂给 Gauss-Newton 迭代循环，配合实测视电阻率数据逐轮更新层参数了。
本文还有配套的精品资源，点击获取
