简介这份资源面向地球物理勘探领域的学习者与工程技术人员提供直流电测深一维正演模拟的完整工具集用于在给定地层电阻率模型下计算对应的视电阻率响应是理解反演流程与解释实测曲线的理论基础。压缩包共4个文件均为MATLAB脚本.m整体约2KB其中包含核心正演模型函数、雅可比矩阵计算模块以及两个不同参数设置的示例脚本分别演示地层参数输入、正演计算与灵敏度关系求解的调用方式。已有230人学习下载适合作为教学演示或算法二次开发的起点。读者可借助示例脚本快速跑通正演流程掌握视电阻率随地层参数变化的规律并在此基础上理解反演中雅可比矩阵的作用为地质勘探与矿产资源调查中的数据处理提供可复用的计算框架。1. 直流电测深正演从 sdc1dfwd.rar 说起视电阻率曲线到底怎么算出来的野外跑直流电测深最怕的不是电极接地电阻大而是解释成果和实际地层对不上。你布好极距测出一串视电阻率画在双对数纸上可这条曲线到底对应几层地、每层多厚、电阻率多少光靠量板法在电脑时代已经不够用了。sdc1dfwd.rar 这类一维直流电正演程序干的就是这件事给定层状地层模型算出理论视电阻率曲线。它解决的是“正演”问题——已知地下结构求地表响应。适合做电法勘探解释的、写反演代码的、以及想搞明白视电阻率曲线为什么长那样的从业者。正演是反演的基础正演都跑不通反演就是盲人摸象。这一章先把直流电测深正演的核心逻辑讲清楚后面再动手复现。2. 直流电测深正演的数学内核从拉普拉斯方程到视电阻率2.1 点电流源在层状大地中的电位分布直流电法正演的理论根基是稳定电流场的拉普拉斯方程。在地面点电流源 A 供电时地下介质满足 ∇·(σ∇U) -Iδ(r-r_A)。对于水平层状介质电阻率只随深度变化σ σ(z)方程在柱坐标系下变成偏微分方程。求解方法通常是分离变量加汉克尔变换。电位 U(r,z) 在每层内满足 ∇²U 0层与层之间满足电位连续和电流密度法向连续。最终地表电位可以写成含贝塞尔函数 J₀(λr) 的积分形式积分核函数由各层电阻率和厚度递推得到。这个递推关系就是正演程序的核心通常用传播矩阵或递推公式实现。实际编程时积分核函数 B(λ) 的递推从最底层开始往上算。最底层是半空间核函数为 1。每往上一层用该层的反射系数和衰减因子更新。反射系数 k_i (ρ_{i1} - ρ_i)/(ρ_{i1} ρ_i)衰减因子 e^{-2λh_i}。递推公式为 B_i (B_{i1} k_i e^{-2λh_i}) / (1 B_{i1} k_i e^{-2λh_i})。这个公式在多数电法教材里都能找到但写代码时容易在指数项符号上翻车。2.2 视电阻率的定义与电极装置系数视电阻率不是真电阻率是假设均匀半空间时按同样装置测得的等效电阻率。对于对称四极装置 AMNB视电阻率 ρ_s K · ΔU / I其中装置系数 K π · (AM·AN) / MN。温纳装置是 AMMNNBK 2πa。施伦贝格装置是 MN 很小K ≈ π·(L² - l²)/(2l)L 是供电极距一半l 是测量极距一半。正演程序通常直接算电位差 ΔU再乘 K 得到 ρ_s。注意不同装置的 K 公式不同程序里如果写死一种装置换装置时一定要改 K 的计算。视电阻率曲线画在双对数坐标下横轴是极距 AB/2纵轴是 ρ_s。对于二层模型曲线首支趋近第一层电阻率 ρ₁尾支趋近第二层 ρ₂中间过渡形态取决于 ρ₂/ρ₁ 和 h₁。三层模型会出现更复杂的极值比如 H 型ρ₁ρ₂ρ₃有极小值K 型ρ₁ρ₂ρ₃有极大值。正演程序的价值就在于快速生成这些理论曲线用于反演拟合或解释对比。2.3 数值积分汉克尔变换的实用处理方法理论公式里的汉克尔变换积分 ∫₀^∞ B(λ) J₀(λr) dλ 没有解析解必须数值计算。常见做法有两种一是直接数值积分把 λ 离散化用高斯求积或辛普森法二是用数字滤波法把积分写成卷积形式用预先算好的滤波系数做快速计算。数字滤波法速度快但滤波系数需要高精度网上流传的系数版本很多用错了曲线尾部会震荡。直接积分法慢但可控适合教学和验证。我一般用直接积分法做验证再用滤波法做批量计算。直接积分时λ 的上限取到 10/r_min 左右下限取 10^{-6}/r_max步长按对数等间隔。积分核 B(λ) 在 λ 很大时趋近 1J₀(λr) 震荡衰减所以截断误差可控。如果曲线在极距很大时出现异常上翘多半是积分上限不够或者核函数递推溢出了。3. 用 Python 复现一维直流电测深正演从模型到曲线3.1 核心递推与积分代码下面这段代码实现三层模型的视电阻率正演装置为温纳直接数值积分。代码不依赖特殊库只用 numpy 和 scipy 的贝塞尔函数。import numpy as np from scipy.special import j0 from scipy.integrate import quad def kernel_b(lam, rho, h): 计算层状介质汉克尔变换核函数 B(λ) rho: 各层电阻率列表从第一层到最底层半空间 h: 各层厚度列表长度比 rho 少 1 n len(rho) b 1.0 # 最底层半空间核函数为 1 for i in range(n-2, -1, -1): k (rho[i1] - rho[i]) / (rho[i1] rho[i]) exp_term np.exp(-2 * lam * h[i]) b (b k * exp_term) / (1 b * k * exp_term) return b def rho_s_wenner(rho, h, ab2): 温纳装置视电阻率正演 ab2: AB/2 极距单位与 h 一致 a ab2 / 1.5 # 温纳装置 AMMNNBa, AB/2 1.5a r a # 测量极距 MN 中点距供电极中点距离温纳取 a def integrand(lam): return kernel_b(lam, rho, h) * j0(lam * r) integral, _ quad(integrand, 0, np.inf, limit200) rho_s 2 * np.pi * a * integral # 温纳 K 2πa return rho_s # 三层 H 型模型ρ1100, ρ210, ρ3100, h15, h210 rho [100, 10, 100] h [5, 10] ab2_list np.logspace(0, 2.5, 30) rho_s_list [rho_s_wenner(rho, h, ab2) for ab2 in ab2_list] for ab2, rs in zip(ab2_list, rho_s_list): print(fAB/2{ab2:.2f} ρs{rs:.2f})这段代码的逻辑说明kernel_b 从最底层往上递推每次用反射系数和衰减因子更新。rho_s_wenner 里温纳装置的 a AB/2 / 1.5因为温纳的 AB 3a所以 AB/2 1.5a。积分用 scipy.integrate.quad 做自适应积分上限设无穷quad 会自动处理震荡。参数说明rho 列表长度必须比 h 多 1因为最底层是半空间没有厚度。h 的单位和 ab2 一致通常用米。如果曲线在 AB/2 很小时偏离 ρ₁检查积分下限是否足够小。3.2 装置系数与极距换算的代码实现不同装置的 K 和极距换算不一样下面把三种常见装置写成函数方便切换。def rho_s_schlumberger(rho, h, ab2, mn): 施伦贝格装置视电阻率正演 ab2: AB/2 mn: 测量极距 MN L ab2 l mn / 2 r np.sqrt(L * l) # 等效极距取几何平均 def integrand(lam): return kernel_b(lam, rho, h) * j0(lam * r) integral, _ quad(integrand, 0, np.inf, limit200) K np.pi * (L**2 - l**2) / (2 * l) rho_s K * integral return rho_s def rho_s_dipole_dipole(rho, h, ab2, mn): 偶极-偶极装置AB 和 MN 分开 ab2: AB/2 mn: MN 长度 # 简化处理取等效极距 r ab2 mn / 2 def integrand(lam): return kernel_b(lam, rho, h) * j0(lam * r) integral, _ quad(integrand, 0, np.inf, limit200) K np.pi * ab2 * (ab2 mn) * (ab2 mn mn) / mn # 近似公式 rho_s K * integral return rho_s逻辑说明施伦贝格装置的 K 公式是 π(L²-l²)/(2l)L 是 AB/2l 是 MN/2。偶极-偶极的 K 公式更复杂这里用近似实际项目要查规范。参数说明mn 是测量极距施伦贝格要求 mn 远小于 ab2否则 K 公式误差大。偶极-偶极的 ab2 和 mn 是独立参数极距组合多正演时注意 r 的取值。如果换装置后曲线形态突变先检查 K 公式和 r 的换算。3.3 用 pygimli 做交叉验证自己写的代码容易有隐藏 bug用成熟库交叉验证是血泪经验。pygimli 是地球物理常用库里面有一维直流电正演。import pygimli as pg import numpy as np # 定义三层模型 model pg.Vector([100, 10, 100]) # 电阻率 thk pg.Vector([5, 10]) # 厚度最后一层半空间 # 创建一维模型 mesh pg.meshtools.createMesh1D(len(model), 1) data pg.DataContainer() data.set(ab2, np.logspace(0, 2.5, 30)) # 正演 f pg.core.DC1dModelling(data, mesh) rhoa f.response(model) print(rhoa)逻辑说明pygimli 的 DC1dModelling 封装了递推和滤波直接给模型返回视电阻率。参数说明model 长度比 thk 多 1和前面一致。ab2 是极距数组。如果 pygimli 结果和自己代码差很多先检查装置类型是否一致pygimli 默认可能是施伦贝格或温纳需要设置。交叉验证能发现积分精度、递推符号、K 系数三类错误。4. 正演程序的避坑与排查那些年我们踩过的雷4.1 曲线尾部震荡或上翘现象AB/2 很大时视电阻率曲线出现周期性震荡或异常上翘不符合物理规律。原因数值积分上限截断太早或者核函数递推在 λ 很大时数值溢出。解决把积分上限提高到 100/r_min 以上或者改用数字滤波法。检查核函数递推里的 exp(-2λh)当 λh 很大时 exp 趋近 0但如果 h 很小而 λ 很大浮点下溢可能产生 NaN。加一个判断当 exp_term 1e-15 时直接置 0。4.2 首支不趋近第一层电阻率现象AB/2 很小时视电阻率不等于 ρ₁偏差很大。原因积分下限太大漏掉了 λ 很小的贡献或者装置系数 K 算错。解决积分下限取到 1e-6/r_max确保覆盖小 λ 区间。检查 K 公式温纳的 K2πa施伦贝格的 Kπ(L²-l²)/(2l)偶极的 K 更复杂。用均匀半空间模型验证如果 ρ₁ρ₂ρ₃100正演结果应该恒等于 100否则 K 或积分有问题。4.3 多层模型递推顺序写反现象三层模型正演结果和理论量板对不上曲线形态完全错误。原因核函数递推从第一层往下算而不是从最底层往上算。递推必须从半空间开始因为半空间的核函数已知为 1。解决循环从 n-2 到 0每次用下一层的 B 值更新当前层。如果写反了反射系数的符号和衰减因子的方向都会错。用二层模型验证ρ₁100ρ₂10h₁5曲线尾支应该趋近 10首支趋近 100。4.4 极距换算和装置类型不匹配现象换装置后曲线整体平移或形态突变。原因AB/2 和 a 的换算关系搞混。温纳的 AB/2 1.5a施伦贝格的 AB/2 L偶极的极距定义更复杂。解决在代码里把装置类型作为参数统一换算到等效极距 r。温纳 ra施伦贝格 rsqrt(L·l)偶极 r 取几何平均。画图前先确认横轴是 AB/2 还是 a两者差 1.5 倍双对数坐标下平移明显。4.5 数值积分速度太慢现象批量计算几百个极距时直接积分法耗时太长。原因quad 自适应积分对每个极距都重新计算核函数重复调用多。解决把核函数 B(λ) 预先算在固定 λ 网格上然后用插值或滤波系数做快速卷积。数字滤波法速度能快几十倍适合反演迭代。如果只是画几条理论曲线直接积分够用。批量反演时建议用 pygimli 或自己实现滤波法。5. 从正演到反演视电阻率曲线的进阶用法正演跑通之后下一步就是反演。反演是给定实测视电阻率曲线反推层状模型参数。常见做法是最小二乘拟合目标函数是正演响应和实测数据的残差平方和用阻尼最小二乘或 Occam 反演迭代。正演速度直接决定反演效率所以滤波法比直接积分更实用。我一般先用直接积分验证几条曲线确认代码无误后换成滤波法做反演。验证正演程序是否正确有一个简单办法用均匀半空间模型视电阻率应该恒等于真电阻率误差小于 0.1%。再用二层模型和教材上的量板对比首支尾支和过渡形态一致。三层模型可以检查极值位置H 型极小值出现在 AB/2 ≈ h₁ h₂ 附近K 型极大值类似。如果极值位置偏移检查厚度和电阻率单位是否一致。进阶技巧用正演程序生成理论曲线库然后对实测曲线做相关分析快速估计层数和大致参数。比如 H 型曲线极小值对应的 AB/2 约等于第一层厚度加第二层厚度极小值大小和 ρ₂/ρ₁ 有关。这些经验规则能帮你在反演前缩小搜索范围。另外正演时注意极距范围要覆盖实测范围外推的曲线不可靠。我自己的习惯是每写一个新正演程序先用均匀半空间和二层模型做单元测试再用三层模型和 pygimli 交叉验证。测试通过后才拿去跑实际数据。这个习惯帮我省了很多后悔药因为正演里的符号错误和积分截断问题在复杂模型下很难肉眼发现。希望帮到你。本文还有配套的精品资源点击获取