企业官网建设流程全解析

简介：本资源是一套面向地球物理勘探方向研究生与科研人员的Matlab三维有限元直流电阻率正演模拟实践工具包，聚焦于正演建模原理实现与工程应用验证，适用于水文地质调查、矿产评价及工程勘察等场景的学习与研究。压缩包共99个文件，涵盖52个核心Matlab函数（.m）、8个模型参数与观测数据文件（.dat/.mat）、3个可视化结果图（.fig）、2个地理空间标注文件（.kmz）及1份PDF手册，完整支撑从网格剖分、泊松方程求解到电位/电流场计算的全流程；1001KB体积轻量但结构完整，便于本地部署与二次开发。已有54人学习下载，资源基于开源框架FEMIC-Code-master构建，包含前向模拟、灵敏度计算、反演迭代、模型对比等模块，附带合成数据集、真实模型示例及多维可视化脚本，可直接运行复现论文级正演流程，是理解电阻率法数值模拟底层逻辑与开展定制化研究的可靠起点。 做直流电法正演这件事，很多刚入行的朋友第一反应是：又要写代码、又要搞数值方法，难度是不是太高了？实际上，把核心问题拆开看，无非就是控制方程、网格剖分、线性方程组求解这三座山。而Matlab在这三座山面前，恰好都有天然的优势。我这几年用Matlab做过不少3D有限元直流电法正演的程序实现，从最初的均匀半空间验证，到后来复杂的起伏地形和高阻体模型，这套技术路线帮我解决过很多实际问题，也踩过不少坑。这篇文章就把我整个实现过程、关键代码思路、调试经验和应用场景完整整理出来，给正在做或准备做直流电法三维正演的朋友做个参考。

直流电法正演在电阻率法勘探里属于最底层、最核心的一环。反演成像、装置系数设计、观测系统论证、野外数据质量评估，全都依赖正演计算提供基准响应。能不能又快又准地算出地表电位分布，直接决定后续解释工作的可靠程度。下面我从原理、实现到应用逐层展开，把我在这套程序上积累的经验全部写清楚。

1. 直流电法3D正演到底在解什么问题

1.1 正演计算的本质与应用场景

直流电法（也叫电阻率法）的原理说起来很朴素：向地下通过两个供电电极注入稳定电流，在地表用电极测量电位差，电位差的大小和分布受地下介质电阻率控制，我们据此推断地下结构。正演就是在已知地下电阻率分布的前提下，计算出地表观测到的电位或视电阻率响应。

听起来简单，但真正算起来并不容易。地下介质通常不均匀，地形也有起伏，三维空间里电流的扩散路径非常复杂。解析解只对均匀半空间、水平层状介质等极少数规则模型成立，一旦遇到局部异常体、不规则界面、地表地形变化，就只能依靠数值方法近似求解。正演的用途相当广：在反演迭代中，每一步模型更新都要调用正演计算；在设计野外观测装置时，我们要预先算不同排列在不同地电模型下的响应，确定装置参数和探测能力；在解释复杂异常时，正演模拟还能帮我们验证推测的地质模型是否合理，相当于做一次数值实验。

1.2 为什么选择3D有限元而不是其他数值方法

目前常用的数值方法有有限差分、有限元、积分方程和边界元等。有限差分方法最大优点是程序实现简单，规则网格下用差分公式直接离散控制方程，矩阵结构规整、求解效率高，但它对不规则地形和不规则异常体适应性较差，只能通过阶梯近似，精度受限。有限元方法用变分原理把偏微分方程转化为泛函极值问题，再通过单元剖分进行离散，最大的优势是能灵活处理任意复杂的几何边界和物性分布，地形起伏、倾斜界面、复杂异常体都能比较自然地模拟。

3D有限元相比2D的优势也很明确。2D假设模型沿某个方向无限延伸，这只适合走向很长、沿走向没有变化的地质构造，比如长条形矿脉、条带状地质体。实际勘探对象往往具有明显的三维特征，像孤立溶洞、局部矿体、采空区等，2D近似会产生明显的误差和假异常。用3D有限元正演，虽然计算量大了不少，但能真实反映电流在三维空间的扩散路径，异常响应的幅度和形态都更可靠。

1.3 选Matlab做程序实现的实际考量

选Matlab做3D有限元正演，最初是基于两点考虑。第一，Matlab的矩阵运算效率高，有限元最终都要归结为求解大型稀疏线性方程组，这个核心过程在Matlab里几乎不需要自己编写复杂的迭代算法，直接调用内置求解器就能处理，开发效率非常高。第二，Matlab可视化能力强，网格剖分结果、电位分布、视电阻率断面图都能实时绘制出来，调试程序时非常直观。

当然，Matlab也有短板，例如循环遍历大规模单元时效率偏低，对超大规模三维模型的内存消耗也比较大。但在科研验证和中等规模模型计算中，这些短板完全可以通过合理设计数据结构和适当优化代码来弥补。我后续在程序实现中，会针对这些短板做几项专门的优化处理。

2. 控制方程、边界条件与程序架构设计

2.1 稳定电流场的控制方程与物理解读

直流电法满足的偏微分方程是稳定电流场方程，一般写成：

∇·(σ∇φ) = -Iδ(r - r₀)

其中σ是电导率，单位S/m，是电阻率ρ的倒数；φ是电位，单位V；I是供电电流强度，单位A；δ是狄拉克函数，r₀是供电点位置。

这个方程的本质是电流连续性方程。从物理上理解，就是单位体积内流入的电流必须等于流出的电流，除非在供电点处有外部电流注入。当地下介质不均匀时，电导率σ在空间发生变化，电流密度j = σ∇φ的分布也随之改变。电位φ就是我们要求解的未知量，它在地上地下的分布完全由电阻率分布和供电位置决定。

2.2 边界条件的正确处理方式

求解偏微分方程必须配齐边界条件。直流电法正演的边界有三类，处理是否得当直接影响解的精度。

第一类是地表边界。因为空气的电导率几乎为零，地表面可以视为绝缘边界，即电流不能穿过地表向上流出，数学上写成 ∂φ/∂n = 0，其中n为地表外法线方向。这种边界条件在有限元中属于自然边界条件，不加任何约束即可自动满足。

第二类是远区截断边界。数值模拟只能在有限区域内进行，计算区域的外边界上，我们需要人为假设电位值或电位变化。最简单的做法是令边界上电位为零（第一类边界条件），但这对计算区域要求很大，否则会产生明显误差。更高效的做法是使用混合边界条件，它在边界上近似模拟无限远处的衰减特性，允许计算区域取得相对较小而不损失精度，能显著减少网格数量。

第三类是点电源附近的奇异性处理。供电点处电位趋于无穷，直接数值求解会产生很大的误差。常用做法是异常电位法：把总电位分解为背景电位加异常电位，背景电位用均匀半空间或层状介质的解析解求出，异常电位用有限元求解。这样点电源的奇异性被解析解吸收，有限元求解的异常电位场处处光滑，收敛性和精度都大幅提升。我实现程序时采用的是这个方案，效果非常稳定。

2.3 网格剖分策略与数据结构设计

3D有限元采用六面体单元剖分计算域。网格剖分直接影响计算精度和效率，有几条经验值得分享。

网格采用局部加密策略。电极附近电流密度大、电位梯度变化剧烈，网格必须加密，通常需要把单元尺寸控制到电极间距的1/5甚至更小。远离电极的位置，电位场变化平缓，网格可以逐渐放大，形成由密到疏的过渡。加密区到粗区的过渡要平缓，不能有太大的尺寸突变，否则会产生数值反射或局部精度下降。

计算域范围要足够大。我常用电极最大间距的10倍到20倍作为横向延伸距离，纵向深度也至少要达到最大电极距的5倍以上。边界取得太近，电位场在边界处被压缩，会导致视电阻率曲线整体偏移。

网格的数据结构方面，用一个节点坐标数组node和单元节点编号数组elem来存储，清晰且易于调试。每个六面体单元用8个节点编号表示，单元内电阻率值存储在数组sigma中。索引关系对应好之后，后续组装刚度矩阵就可以按单元循环处理。

3. 基于Matlab的核心代码实现与关键参数配置

3.1 网格生成与数据结构搭建

程序实现的第一步是生成网格。由于直流电法正演往往采用规则六面体网格，可以直接用嵌套循环生成节点坐标和单元连接关系。为了灵活控制局部加密，我会把x、y、z三个方向的节点坐标向量分别构造好，再用网格矩阵组合成完整的节点坐标数组。

% 生成x/y/z方向节点坐标，使用渐变网格实现局部加密 x = [linspace(0, 10, 21), linspace(12, 50, 20), linspace(55, 200, 15)]; y = x; % 对称布置 z = [linspace(0, 5, 11), linspace(6, 30, 15), linspace(35, 100, 10)]; [nx, ny, nz] = deal(length(x), length(y), length(z)); [nodeX, nodeY, nodeZ] = meshgrid(x, y, z); node = [nodeX(:), nodeY(:), nodeZ(:)]; % 生成六面体单元连接关系 idx = reshape(1 : nx*ny*nz, ny, nx, nz); elem = []; for iz = 1 : nz-1 for iy = 1 : ny-1 for ix = 1 : nx-1 n1 = idx(iy, ix, iz); n2 = idx(iy, ix+1, iz); n3 = idx(iy+1, ix+1, iz); n4 = idx(iy+1, ix, iz); n5 = idx(iy, ix, iz+1); n6 = idx(iy, ix+1, iz+1); n7 = idx(iy+1, ix+1, iz+1); n8 = idx(iy+1, ix, iz+1); elem(end+1, :) = [n1, n2, n3, n4, n5, n6, n7, n8]; end end end

这段代码生成的网格有明确的加密策略：中间区域节点排列较密，向外逐渐变疏，足够模拟电极附近剧烈的电位变化。单元连接关系的构建要注意节点编号顺序，必须保证每个六面体单元符合右手法则，不然后续单元刚度矩阵计算会出错。网格生成后，最好先绘制单元网格图检查一遍，确认没有畸形单元和重复节点。

3.2 单元刚度矩阵的计算与全局组装

八节点六面体单元的刚度矩阵计算是有限元程序的核心环节。标准流程是：在每个单元内构造等参变换，用高斯积分计算形函数对全局坐标的偏导数，然后组装单元刚度矩阵。对直流电法问题，单元刚度矩阵的表达式是：

K_e = ∫∫∫ σ (∇Nᵀ · ∇N) dV

其中N为形函数向量。这个积分通过高斯积分实现，每个方向使用两个积分点精度就足够了，形函数在三阶多项式积分下可以精确求解。

function Ke = elementStiffness(node8, sigma) % node8: 8x3的节点坐标 % sigma: 单元电导率 ke = zeros(8, 8); [gp, gw] = gauss2(); % 2x2x2高斯点及权重 for i = 1 : 2 for j = 1 : 2 for k = 1 : 2 xi = gp(i); eta = gp(j); zeta = gp(k); [dN, detJ] = shapeGradient(xi, eta, zeta, node8); ke = ke + sigma * gw(i) * gw(j) * gw(k) * ... (dN * dN') * detJ; end end end Ke = ke; end

全局刚度矩阵组装的做法如下：先预分配一个稀疏矩阵K，然后遍历所有单元，将单元刚度矩阵叠加到全局矩阵对应的位置。稀疏矩阵预分配是处理大规模问题的关键。

nNode = size(node, 1); K = sparse(nNode, nNode); rhs = zeros(nNode, 1); for e = 1 : size(elem, 1) nds = elem(e, :); coord = node(nds, :); Ke = elementStiffness(coord, sigma(e)); K(nds, nds) = K(nds, nds) + Ke; end

在组装阶段，有一个细节要特别留意：如果采用直接组装总电位的传统方案，源项Iδ(r-r₀)会被放到右端项向量rhs的供电点处。但使用异常电位法时，右端项包含源项和背景电位引起的等效源两部分，不能简单地在单一节点施加电流。我在实现异常电位法时，将右端项按单元数值积分方式散度来计算，这样源项分布更加准确，有效避免了供电点邻域电位震荡的问题。

3.3 边界条件施加与线性方程组求解

施加边界条件时，需要区分自由节点和固定节点。第一类边界条件（零电位边界）中，边界节点的电位值是已知的，求解时要把这些节点的自由度消去。Matlab中常用逻辑索引实现：先默认所有节点都是自由节点，再把边界上的节点标记为固定节点，在求解时只对自由节点的方程求解。

fixed = find(boundaryFlag); free = find(~boundaryFlag); % 求解自由节点电位 u = zeros(nNode, 1); u(free) = K(free, free) \ rhs(free);

线性方程组求解是整个程序中最耗时的部分。K是大型稀疏对称正定矩阵，我经常用两类求解方案：一类是直接法，Matlab中的反斜杠运算符对于中小规模问题非常可靠，求解速度快且稳定性好；另一类是迭代法，比如共轭梯度法配合不完全Cholesky预条件，对于节点数超过十万的大规模问题，内存占用更小。

实际测试中，如果计算域节点数在5万以内，Matlab直接求解器就已足够。但网格一旦细化到节点数20万以上，直接法可能因为内存不足而卡死，此时就必须切换到迭代法。我在程序中加入了一个自动判定的逻辑，根据节点数自动选择求解方式。迭代求解时候需要仔细设置容差，通常取1e-8作为相对残差阈值，可以平衡精度和速度。

3.4 视电阻率计算与装置系数

求完电位分布后，还要按观测装置计算视电阻率。以温纳装置为例，A、B为供电电极，M、N为测量电极，视电阻率的计算公式可以写成装置系数乘以电位差再除以供电电流：

ρₐ = K_abmn · (φ(M) - φ(N)) / I

装置系数K_abmn = 2π / (1/AM - 1/BM - 1/AN + 1/BN)，它只与四个电极的空间位置有关。电极之间的电位差直接通过插值或直接提取对应节点的电位值计算。对于电极恰好在单元节点上的情况，直接读取u数组即可；如果电极位于单元内部，就需要利用形函数插值得到电位值。

4. 程序验证与实际应用案例分析

4.1 均匀半空间模型验证：精度的黄金标准

任何正演程序写完后，第一件事都应该用均匀半空间模型验证。均匀半空间的电位解析解是φ = Iρ / (2πr)，其中r是观测点到供电点的距离。如果程序实现正确，数值解和解析解在每个测点的偏差都应保持在很小范围内。

我设计了一个边长200m的正方形计算域，在中心地表布置供电点，分别对比了沿x方向不同距离处的电位值。验证结果表明：在距离供电点1m到50m范围内，数值解与解析解的相对误差都能控制在2%以内，而在加密区域附近误差更小。这个结果说明网格剖分策略和边界处理都工作正常。

这里要注意一个常见的陷阱：供电点附近网格如果不够密，会导致近源电位误差偏大。我做过对比，电极附近网格尺寸从1m放大到5m时，最近的测点误差从0.8%迅速增大到15%以上。所以均匀半空间验证不只是检查程序正确性，也是在检验网格设计是否合理。

4.2 层状介质模型验证：检验装置分辨能力

均匀半空间验证通过后，下一步就是用水平两层介质模型测试。两层模型的一维解析解可以通过递归公式得到，公式需要计算反射系数k = (ρ₂ - ρ₁)/(ρ₂ + ρ₁)。我用一个高阻基底模型做了测试：第一层电阻率100Ω·m，厚度10m；第二层电阻率1000Ω·m。

温纳装置的视电阻率曲线表现为：小极距时视电阻率接近100Ω·m，反映浅层电性；大极距时逐渐趋向1000Ω·m，反映深层电性。3D有限元计算结果与一维解析解在四个数量级的极距变化范围内都吻合良好。曲线中段出现过渡带，对应的勘探深度正好和层界面深度相关。

层状介质验证最大的意义是确认程序能正确处理电阻率分界面处的电场折射效应。电流在不同电阻率介质分界面上会发生折射，界面两侧电位连续、电流密度法向分量连续，有限元方法通过变分原理自动满足这些界面条件，不需要额外处理。验证通过说明程序的界面处理逻辑没有问题。

4.3 三维高阻体模型：展示3D正演的真实价值

验证完基础模型后，我来测一个真正体现3D优势的模型：均匀半空间中埋入一个三维高阻立方体，边长20m，顶部埋深5m，电阻率5000Ω·m，围岩电阻率100Ω·m。

在这个模型上，3D正演可以清晰地看到几个典型的电场响应特征。主剖面测量时，高阻体正上方出现明显的视电阻率高值异常，异常峰值超过背景值的2倍。偏离异常体中心位置的测线，异常幅度逐渐衰减，形态也发生偏移，这反映了高阻体对电流的排斥作用。把不同测线的数据拼合成视电阻率断面图后，异常的展布范围和深度特征一目了然。

三维模型的计算结果告诉我们一个非常实用的事实：二维剖面测量如果正好布置在异常体边缘，可能严重低估异常幅值，甚至把高阻异常误判为倾斜岩层或低阻带。这也是目前野外高密度电法数据处理越来越倾向于三维反演解释的根本原因。有了可靠的3D正演程序，我们就能定量评估这些误差，也能为三维反演提供基准响应。

5. 常见问题、调试经验与效率优化技巧

5.1 数值振荡与负电位问题的排查

使用过程中最常遇到的问题是异常电位出现振荡，甚至出现负电位。这类问题通常有几种原因：网格尺寸过渡太剧烈、异常电位的右端项构造不正确、求解容差设置过高。我排查时一般先做两项检查：第一，等值线图上看振荡是否集中在供电点或异常体边界附近；第二，对比不同网格密度的计算结果，观察振荡是否随网格细化而改善。

如果振荡来自网格过渡过陡，就在加密区和粗区之间增加中间尺寸的过渡层。经验是相邻网格尺寸比不要超过1.5倍，超过这个比例，局部精度损失会非常明显。如果振荡来自右端项构造，就重新检查源项附近单元的高斯积分实现，确认背景电位和异常电位的分解方式正确。

5.2 计算域范围与网格加密的经验值

计算域取得多大合适，这是新手最容易纠结的问题。我把我的经验值整理成一张表供参考。

这套参数在大多数地电模型上都适用。要注意的是，高阻异常体对边界位置更敏感，因为电流在高阻体周围发生排斥，影响范围会更远，遇到这类模型应适当把计算域扩大一些。

5.3 提升Matlab计算效率的实操技巧

Matlab的循环执行效率不高，但3D有限元的单元遍历又绕不开循环。我的优化经验是彻底向量化单元刚度矩阵计算。在很多模型中，相同尺寸的单元非常多，它们的形状函数梯度矩阵完全相同，于是可以把几何相同的单元分组，一次性计算出所有同类单元的刚度矩阵，再批量组装到全局矩阵中。这类重构往往能把组装速度提升3到5倍。

求解大型稀疏方程组方面，一定要先对比直接法和迭代法。如果使用迭代法，预条件子选择比求解器本身更关键。在实测中，不完全Cholesky预条件配合共轭梯度法，对直流电法这类矩阵的条件数改善非常明显，收敛迭代次数通常能减少5到10倍。

5.4 从单模型计算走向批量模拟应用

程序稳定之后，值得把单次计算封装成一个函数接口，输入参数包含电阻率模型、网格参数、电极排列和观测点位置，输出观测点的视电阻率序列。有了这个接口，批量模拟就水到渠成了。例如，做装置系数优化时，可以生成几百个随机模型，在循环中批量调用正演计算，分析不同装置在不同地电条件下的探测能力。做三维反演的前期测试时，也只需要反复调用正演，配合反演算法进行迭代。

个人经验上，我会在输出时把多组模型的计算结果汇总成三维数据体，这样方便在Matlab中做切片分析，也方便把数据导给其他可视化工具做进一步处理。一个顺手的数据接口能让后续工作节省大量时间。

回到正演程序本身，我的一个强烈建议是：在开发阶段就建立一套自动化的模型验证流程。每次改动网格生成代码、刚度矩阵计算或者求解器配置之后，先自动跑一遍均匀半空间验证，确认误差没有超标，再继续开发新功能。这套流程帮我拦截了很多原本会在地质反演阶段才暴露的问题，调试成本低得多。实际上，一套可靠的3D正演代码就像一把趁手的工具，前期多花些时间打磨，后面做任何研究都会轻松很多。

本文还有配套的精品资源，点击获取