简介这份资源面向地球物理勘探、电法勘探方向的学习者与科研人员提供直流电测深一维正演模拟的MATLAB程序用于在给定地层电阻率模型下计算视电阻率响应是理解正演与后续反演的基础工具。压缩包共4个文件均为m脚本与函数整体约2KB体积轻量便于快速上手。其中包含正演核心模型函数、雅可比矩阵计算模块以及两个示例脚本分别演示不同地层参数设置与调用流程用户可修改参数适配不同地质场景。已有230人学习下载说明其在教学与科研中具有一定参考价值。通过阅读与运行这些代码读者可掌握直流电场偏微分方程的数值求解思路理解地层参数变化对视电阻率曲线的影响规律为反演建模、实测数据解释及矿产资源调查提供理论框架与计算支持适合具备一定MATLAB基础的中高级学习者。1. 直流电测深正演从 sdc1dfwd 这个包名说起如果你手头有一个名为sdc1dfwd.rar的压缩包解压后看到的是 Fortran 或 C 写的一维直流电测深正演代码那它大概率属于上世纪九十年代到本世纪初电法勘探教学与科研中流传的一类经典程序。这类程序的核心任务只有一个给定地下介质的层状电阻率模型和电极距算出对应的视电阻率曲线。听起来简单但真正把它跑通、跑对、跑出能用于反演的数据中间隔着一堆参数约定和数值细节。直流电正演在电法勘探里的地位相当于地震勘探里的射线追踪——它是反演的基础也是理解观测数据的钥匙。测深正演尤其如此因为一维层状模型虽然简化却是绝大多数野外电测深数据解释的起点。视电阻率曲线长什么样、在什么位置出现拐点、尾支渐近值是多少这些判断都依赖正演结果的可靠性。这篇文章面向的是需要复现或改造这类程序的从业者不管你是研究生刚接手师兄的代码还是工程师要验证一套野外数据的合理性下面这些内容应该能帮你少走弯路。2. 一维直流电测深正演的理论骨架与数值实现2.1 从拉普拉斯方程到视电阻率核心公式链直流电测深正演的物理基础是稳定电流场在层状介质中的分布。点电流源位于地表时地下电位满足拉普拉斯方程在圆柱坐标系下通过分离变量得到电位表达式。对于 n 层水平层状介质地表电位可以写成含贝塞尔函数 J₀ 的积分形式V(r) (ρ₁I / 2π) ∫₀^∞ K(λ) J₀(λr) dλ其中 K(λ) 是核函数由各层电阻率和厚度通过递推公式确定。视电阻率则通过测量电位差反推ρₐ 2πr · ΔV / I实际编程时核函数 K(λ) 的递推是数值稳定性的关键。常见做法是从最底层向上递推利用双曲函数或指数函数构造反射系数。这里有一个血泪经验如果层厚与电极距比值跨了几个数量级直接算双曲函数容易溢出必须做归一化处理。注意不同教材对核函数的定义可能差一个 ρ₁ 因子移植代码时先确认视电阻率的归一化方式否则曲线整体平移一个常数倍。2.2 数值积分汉克尔变换的离散化选择核函数确定后剩下的就是算积分。被积函数包含振荡衰减的贝塞尔函数直接数值积分效率极低。工程上普遍采用汉克尔变换的离散化方法把积分变成加权求和ρₐ(r) Σᵢ wᵢ · K(λᵢ/r) / r其中 λᵢ 和 wᵢ 是预先算好的滤波系数。常用的有 47 点、61 点、120 点等不同精度的滤波器。点数越多能覆盖的电极距范围越宽但计算量也越大。我一般会先用 47 点滤波器快速验证模型确认无误后再换高精度滤波器出正式结果。下面是一个用 Python 复现核心计算流程的最小示例假设你已经从原程序中提取了滤波系数import numpy as np # 47点汉克尔变换滤波系数示例实际需从原程序或文献获取完整表 lambda_vals np.array([...]) # 长度47 weights np.array([...]) # 长度47 def kernel_function(rho, h, lam): 计算层状介质核函数 K(λ) rho: 各层电阻率列表长度n h: 各层厚度列表长度n-1最后一层为半空间 lam: 积分变量 n len(rho) # 从底层向上递推 # 底层反射系数为0半空间 R 0.0 for i in range(n-2, -1, -1): # 层间反射系数 r_i (rho[i] - rho[i1]) / (rho[i] rho[i1]) # 双曲正切形式避免溢出 th np.tanh(lam * h[i]) R (r_i R * np.exp(-2 * lam * h[i])) / (1 r_i * R * np.exp(-2 * lam * h[i])) return R def apparent_resistivity(rho, h, r): 计算电极距r处的视电阻率 total 0.0 for lam, w in zip(lambda_vals, weights): total w * kernel_function(rho, h, lam / r) return rho[0] * total # 归一化因子视具体定义调整 # 三层模型示例ρ1100, ρ210, ρ3100h110, h220 rho [100.0, 10.0, 100.0] h [10.0, 20.0] r_values np.logspace(0, 3, 50) # 电极距从1到1000米 rho_a [apparent_resistivity(rho, h, r) for r in r_values]这段代码的逻辑说明kernel_function从最底层半空间开始逐层向上计算反射系数 R。递推中用了np.exp(-2*lam*h)而不是直接算双曲函数是为了在 lam*h 很大时避免数值溢出。apparent_resistivity对每个电极距 r把滤波系数和核函数加权求和。参数说明rho和h的长度关系必须满足len(h) len(rho) - 1最后一层厚度无穷大不参与递推。实际使用时滤波系数的精度直接决定结果可信度建议从原程序的数据文件中读取不要自己手算。2.3 层状模型参数化厚度与电阻率的约定原程序sdc1dfwd这类代码通常对输入格式有固定约定。常见的有两种一种是第一行给层数然后每行给厚度和电阻率最后一层只给电阻率另一种是把厚度和电阻率分成两个数组分别读入。移植或调用时最容易翻车的地方是厚度单位——有的程序用米有的用千米还有的用相对第一层厚度的归一化值。视电阻率曲线对厚度单位极其敏感差一个量级整条曲线形态完全变样。我一般会先用一个已知解析解的两层模型做验证第一层厚度 h1电阻率 ρ1底层电阻率 ρ2。当电极距远小于 h1 时视电阻率趋近 ρ1当电极距远大于 h1 时趋近 ρ2。如果跑出来的曲线尾支不对优先检查厚度单位和滤波系数是否匹配。3. 把 sdc1dfwd 跑起来编译、输入文件与结果验证3.1 老代码的编译环境搭建这类正演程序多数是 Fortran 77 或 Fortran 90 写的少数用 C。在 Windows 上原包可能附带.exe但如果你要改参数或移植到 Linux 服务器就需要自己编译。推荐用gfortran大多数老代码稍作修改就能通过。# 安装gfortranUbuntu/Debian sudo apt-get install gfortran # 编译单个源文件 gfortran -o sdc1dfwd sdc1dfwd.f -O2 -stdlegacy # 如果有多个文件 gfortran -o sdc1dfwd main.f kernel.f filter.f -O2 -stdlegacy参数说明-O2开启常规优化-stdlegacy允许老式语法如固定格式、隐式类型。如果报错涉及PAUSE或DATA语句加-stdlegacy通常能解决。编译成功后先别急着跑复杂模型用一个两层模型试一下确认程序能正常输出。3.2 输入文件格式与参数含义假设原程序的输入文件叫model.in典型格式可能是3 10.0 100.0 20.0 10.0 100.0第一行是层数 n接下来 n 行每行两个数厚度和电阻率最后一层厚度省略或写 0。但不同程序差异很大有的要求第一行给电极距个数和最小最大电极距有的把滤波系数单独放一个文件。拿到一个不熟悉的包第一步永远是找一个输入示例文件对照代码里的READ语句逐行确认。提示如果压缩包里没有示例输入用文本编辑器打开源文件搜索READ关键字按变量声明顺序反推输入格式。这是最笨但最可靠的办法。3.3 用两层模型验证正演结果跑通程序后必须做验证。两层模型是最简单的试金石。设 ρ1100 Ω·mh110 mρ210 Ω·m。理论曲线在 r 10 m 时接近 100在 r 100 m 时接近 10中间过渡段单调下降。如果程序输出曲线在中间出现振荡或尾支不收敛问题通常出在滤波系数点数不够或核函数递推溢出。# 读取程序输出并绘图验证 import matplotlib.pyplot as plt # 假设输出文件两列电极距 视电阻率 data np.loadtxt(output.dat) r data[:, 0] rho_a data[:, 1] plt.loglog(r, rho_a, b-, label正演结果) plt.axhline(y100, colorr, linestyle--, labelρ1100) plt.axhline(y10, colorg, linestyle--, labelρ210) plt.xlabel(电极距 r (m)) plt.ylabel(视电阻率 (Ω·m)) plt.legend() plt.grid(True, whichboth) plt.show()逻辑说明用双对数坐标绘图是电测深的标准做法因为视电阻率曲线在双对数坐标下形态更直观。两条水平参考线分别对应第一层和底层电阻率正演曲线应该从 ρ1 平滑过渡到 ρ2。如果曲线在过渡段出现明显上凸或下凹说明中间层参数或滤波系数有问题。4. 正演程序移植与改造中的避坑清单4.1 滤波系数不匹配导致曲线整体偏移现象正演曲线形态正确但整体比理论值高或低一个固定倍数。原因滤波系数对应的归一化方式与代码中的视电阻率公式不一致。不同文献的汉克尔变换滤波器可能对应不同的积分定义有的把 1/r 因子吸收进权重有的没有。解决先用一个均匀半空间模型验证。均匀半空间下视电阻率应恒等于 ρ1与电极距无关。如果跑出来是一条水平线但数值不对调整归一化因子即可。4.2 层厚过大或过小引发数值溢出现象程序运行中报浮点数溢出或输出 NaN。原因核函数递推中exp(-2*lam*h)在 lam*h 很大时下溢到 0本身不报错但如果代码里用了exp(2*lam*h)就会上溢。另外当层厚与电极距比值超过 10^6 时双曲函数直接计算也会不稳定。解决统一用指数衰减形式并在递推前对 lam 做截断。我一般会在核函数里加一个判断如果lam * h 50直接令exp(-2*lam*h) 0避免无谓计算。4.3 电极距数组与滤波系数点数不匹配现象程序只输出前几个电极距的结果后面的全是零或重复值。原因滤波系数的点数决定了能有效计算的 λ 范围进而限制了电极距范围。47 点滤波器通常能覆盖 3~4 个数量级的电极距如果你要算 5 个数量级必须换更多点的滤波器。解决先确认滤波器点数再根据r_min和r_max估算所需范围。常见做法是电极距范围不超过滤波器有效范围的 80%。4.4 输入文件编码与空格问题现象程序读入参数后计算结果完全离谱或直接崩溃。原因老 Fortran 程序对输入格式极其挑剔固定格式下空格和制表符含义不同Windows 换行符\r\n也可能导致读入错误。解决用dos2unix转换换行符输入文件用纯空格分隔避免制表符。如果程序用列表导向读入READ(*,*)对格式容忍度稍高但仍建议保持整洁。4.5 视电阻率定义中的装置系数混淆现象正演结果与文献对比差一个装置系数。原因直流电测深常用对称四极装置视电阻率公式为ρₐ K · ΔV / I其中装置系数K π · (r² - a²) / (2a)或类似形式取决于电极距定义。有的程序把 r 定义为供电电极距的一半有的定义为整个供电极距。解决确认代码中 r 的物理含义再对照装置系数公式。最可靠的办法是用一个已知解析解的三层模型对比商业软件或文献曲线。5. 从正演到反演用视电阻率曲线做模型敏感性分析正演程序跑通只是第一步真正体现价值的地方在于用它做敏感性分析为反演提供先验约束。一个具体技巧是固定层数系统改变某一层的电阻率或厚度观察视电阻率曲线的变化幅度和位置。变化越大的电极距区间说明该参数在该区间最敏感反演时应该给这些数据更大的权重。# 敏感性分析改变第二层电阻率观察曲线变化 rho_base [100.0, 10.0, 100.0] h_base [10.0, 20.0] r_values np.logspace(0, 3, 50) plt.figure() for rho2 in [5, 10, 20, 50]: rho_test [100.0, rho2, 100.0] rho_a [apparent_resistivity(rho_test, h_base, r) for r in r_values] plt.loglog(r_values, rho_a, labelfρ2{rho2}) plt.xlabel(电极距 r (m)) plt.ylabel(视电阻率 (Ω·m)) plt.legend() plt.grid(True, whichboth) plt.show()这段代码的逻辑是保持第一层和第三层不变只改第二层电阻率看曲线在哪个电极距区间分叉最明显。参数说明rho2从 5 到 50 变化覆盖了低阻到高阻的情况。实际解释时如果野外数据在某个电极距区间分辨率高反演结果对该层电阻率的约束就强反之如果曲线对某参数几乎不敏感反演时就应该固定该参数或加正则化约束。我自己的习惯是每拿到一套新数据先跑一遍正演敏感性分析搞清楚哪些参数是数据能约束的哪些是数据约束不住的。这个习惯帮我避免了很多次反演跑飞的情况。正演程序不是黑匣子把它当成理解数据的手段比单纯拿来生成训练集有价值得多。希望帮到你。本文还有配套的精品资源点击获取