简介面向地球物理勘探、应用数学及计算力学领域的科研人员与研究生这份PDF资源复现了基于有限元法的三维电阻率测深数值模拟论文可用于复杂地电结构建模与视电阻率计算。内容涵盖点源电场边值问题、六面体单元剖分、三线性插值、单元刚度矩阵推导及高斯积分实现并附带完整可运行的Python代码与逐步解释从形函数、雅可比矩阵到刚度矩阵组装与稀疏方程组求解步骤一应俱全帮助读者掌握有限元离散到数值计算的完整流程。压缩包仅含1个PDF文件大小456KB便携易用适合在野外勘探数据解释、地电异常体响应分析等场景中作为理论参考与仿真工具。目前已有75人学习下载对于需要复现三层水平地层、低阻矿体、多矿体等典型算例并拓展至实际应用的研究者而言具有较强的实践指导价值。1. 从论文公式到可复现代码三维电阻率测深正演到底卡在哪很多做物探方向的同学第一次翻到“基于有限元法的三维电阻率测深数值模拟”这类论文时心里是先喜后慌喜的是公式看起来不太多慌的是看完摘要发现还是不知道点电流源怎么变成刚度矩阵右端项、视电阻率怎么从几百个节点电位里换算出来。这种论文往往还附带一句“含详细代码及解释”结果附录里要么是十几页 Fortran 片段要么代码缺了最关键的边界处理。这也是华为杯数学建模大赛与研究生数学建模里经常被选中的一类真题——三维正演代码要能复现本质上考验的是把偏微分方程离散、网格生成、求解器和视电阻率公式串成一条完整流水线的工程能力。这篇笔记不转述任何一篇具体论文的原文而是按“有限元法求电位 → 复杂地电结构建模 → 视电阻率计算 → 避坑 → 验证”这条主线把一套常见的 Python 三维电阻率测深正演方案讲清楚。读者里如果有做工程勘探、地热、矿产调查的工程师或者正在复现同类论文的学生可以直接照着把代码跑通再根据自己模型去改网格和电性参数。2. 控制方程与有限元离散让泊松方程变成可解的稀疏方程组三维电阻率测深的正演问题本质上是在已知电阻率分布和供电电流的前提下求解地表的电位分布。这个问题的数学核心是直流电法满足的偏微分方程而有限元法要做的就是把连续方程离散成节点电位未知数的线性方程组。2.1 直流电场控制方程与三类边界条件无源区域内的电位 φ 满足∇·(σ∇φ) 0其中 σ 是电导率等于电阻率 ρ 的倒数。供电电极处电流注入控制方程右端要加上点源项写成 −∇·(σ∇φ) I·δ(r−r⁺) − I·δ(r−r⁻)其中 r⁺ 和 r⁻ 分别是正极和负极位置。边界条件是这个问题的重头戏。地表与空气接触电流不能穿过地表因此地表满足诺伊曼条件 ∂φ/∂n 0。计算区域的侧向边界和底边界不可能取得无限远常见做法是给混合边界条件 Robin∂φ/∂n βφ 0β 一般取 1/rr 是边界到供电源的距离。这个条件比单纯把边界设为零电位或绝缘边界更接近物理实际能减少截断边界带来的误差。在地电分界面上电位连续电流密度的法向分量连续有限元法的弱形式天然满足这个条件不需要额外处理。这就是有限元相对有限差分的一个优势对于地形起伏、断层、倾斜低阻体这类复杂界面单元的物理属性可以随单元位置任意变化不必像有限差分那样把网格边界对齐到地电界面。2.2 加权残量法推导弱形式并引入单元插值对 ∇·(σ∇φ) 0 乘以任意试探函数 w并在计算域 Ω 上积分再用格林公式分部积分得到∫Ω σ∇φ·∇w dV − ∮Γ w·σ(∂φ/∂n) dS 0把 Robin 边界条件代入边界积分项后可以得到包含边界修正项的弱形式。注意点电流源项没有被直接积分因为在离散时它表现为右端向量中对应节点上的集中电流值。这是有限元处理点源的标准手段也被称为“源项分配到节点”。接下来把计算域剖分为六面体单元。每个单元内电位 φ 用节点值和形函数 N 展开φ(x) Σ N_i(x)·φ_i。把形函数代入弱形式就得到单元刚度矩阵K_e[i][j] ∫e σ(x)·∇N_i·∇N_j dV这里的积分用高斯数值积分完成。三维六面体单元一般取 2×2×2 个高斯点对线性八节点单元来说这个积分精度已经足够。每个单元按各自电阻率赋值因此异常体边界不需要和网格边界重合只要异常体覆盖的单元被标记出来就行。2.3 八节点六面体单元刚度矩阵的 Python 实现单元刚度矩阵是有限元里最容易被写错的部分尤其雅可比矩阵的行列式符号和形函数导数的节点顺序。下面这段代码用标准等参单元计算 8×8 的局部刚度矩阵每个高斯点处计算物理空间梯度再累加积分结果。import numpy as np def hex8_local_matrix(nodes_xyz, sigma): nodes_xyz: 8x3 数组, 单元8个节点的三维坐标 sigma: 该单元的电导率 ref_node: 局部坐标下8个节点的参考坐标 ref_node np.array([ [-1, -1, -1], [1, -1, -1], [1, 1, -1], [-1, 1, -1], [-1, -1, 1], [1, -1, 1], [1, 1, 1], [-1, 1, 1] ], dtypefloat) # 2x2x2 高斯点 gp np.array([ [-1, -1, -1], [1, -1, -1], [1, 1, -1], [-1, 1, -1], [-1, -1, 1], [1, -1, 1], [1, 1, 1], [-1, 1, 1] ]) / np.sqrt(3.0) Ke np.zeros((8, 8)) for q in gp: # dN: 3x8, 形函数对局部坐标的导数 dN np.zeros((3, 8)) for i, r in enumerate(ref_node): dN[0, i] 0.125 * r[0] * (1 r[1]*q[1]) * (1 r[2]*q[2]) dN[1, i] 0.125 * r[1] * (1 r[0]*q[0]) * (1 r[2]*q[2]) dN[2, i] 0.125 * r[2] * (1 r[0]*q[0]) * (1 r[1]*q[1]) J dN nodes_xyz # 3x3 雅可比矩阵 detJ np.linalg.det(J) invJ np.linalg.inv(J) grad_phi invJ dN # 物理空间梯度, 3x8 Ke sigma * detJ * (grad_phi.T grad_phi) return Ke这段代码的关键在于每个高斯点上先通过形函数局部导数 dN 与节点坐标相乘得到雅可比矩阵 J再求逆变换到物理空间。detJ 必须为正如果输入的节点坐标顺序是逆序或单元严重畸变detJ 会变成负值组装出来的矩阵会有一排负特征值求解直接发散。另一个容易被忽略的细节是 sigma 按单元而不是按节点赋值这样才能表达复杂地电结构。2.4 全局组装把单元矩阵放到稀疏总矩阵里单元矩阵算完后要按全局节点编号放进大矩阵。三维网格的节点编号使用 i jnx knx*ny 这种方式最方便其中 i、j、k 是节点在 x、y、z 方向的下标。下面的组装代码用 coordinate (COO) 格式收集非零元素最后转成 CSR 格式给稀疏求解器用。from scipy.sparse import coo_matrix def assemble_global(nodes, elements, sigma_e): nn len(nodes) rows, cols, vals [], [], [] for eid, ele in enumerate(elements): ke hex8_local_matrix(nodes[ele], sigma_e[eid]) for a in range(8): for b in range(8): rows.append(ele[a]) cols.append(ele[b]) vals.append(ke[a, b]) K coo_matrix((vals, (rows, cols)), shape(nn, nn)).tocsr() return K组装完成后线性方程组 K·Φ F 的左端就是大型稀疏矩阵。实测经验是网格节点数在 20 万以内时直接稀疏 LU 分解求解器还能承受再往上就要换迭代求解器加预处理了。组装阶段尽量用向量化循环替代逐单元双重循环不然纯 Python 三重循环会把绝大部分时间耗在组装而不是求解上。3. 复杂地电结构建模网格策略与单元属性赋值三维电阻率测深的正演精度一半取决于求解器另一半取决于网格怎么剖、电性参数怎么填。很多刚上手的人拿着论文里的模型图直接套均匀网格结果电极之间距离只有几个网格步长源点附近的电位解完全失真。网格不是越细越好关键是让电极附近和浅层的单元足够小深部和远离测线的区域逐步放大。3.1 网格生成策略测线方向加密垂向指数拉伸合理的网格需要照顾两个精度要求一是供电极和测量极所在位置的单元不能太大否则点源近似误差迅速放大二是深部单元的尺寸可以放大因为电位随距离衰减对浅层网格的依赖远大于深层。水平方向采用“测线附近细、两侧粗”的过渡网格垂向从地表开始按指数因子增长。下面给出一个常用的轴向节点生成函数把测线覆盖范围和电极间距作为输入自动生成非均匀节点坐标。def build_axis_nodes(center, half_width, dense_width, dx_min, ratio1.4): nodes [] # 从测线中心向两侧生成 x center nodes.append(x) while x center half_width: dx dx_min if abs(x - center) dense_width: dx dx_min * ratio ** (abs(x - center) / dense_width) x min(x dx, center half_width) nodes.append(x) left [center - x for x in nodes[1:]] left.reverse() return left nodes # 垂向节点: 地表 z0, 向下指数增长到 zmax def build_depth_nodes(zmax, dz00.5, ratio1.25): nodes [0.0] z 0.0 while z zmax: dz min(dz0 * ratio ** len(nodes), zmax - z) z max(dz, 0.01) nodes.append(z) return np.array(nodes)水平方向的 dx_min 一般取最小电极距的 1/5 到 1/10。电极只在 100 米范围内布置时dx_min 取 2 米就够如果最小电极距只有 5 米dx_min 至少要到 0.5 米。垂向第一层 dz0 取最小电极距的 1/5 左右然后按 1.21.4 的倍率往下拉。地电模型里如果有大型低阻体异常体边界附近的网格密度同样要增加否则视电阻率剖面会在边界处出现明显的振荡。3.2 把复杂地电模型映射到单元层状背景加异常体标记复杂地电结构在代码里的本质就是给每个单元一个电阻率值。层状背景很容易实现先按深度索引给所有单元填背景值再把位于某个深度范围内的单元改成另一层电阻率。不规则异常体要用位置条件去圈定比如球体、倾斜板、断层破碎带本质上是对单元中心坐标做集合判断。下面这段代码演示了如何在背景半空间上叠一个倾斜低阻板和一个低阻球体。所有单元先按层位填背景电阻率再用布尔索引覆盖异常体区域。def assign_rho_to_cells(xc, yc, zc, cells, rho_bg100.0): xc, yc, zc: 单元中心点坐标数组, 长度等于单元数 cells: 单元节点编号列表, 用于后续组装 rho: 返回长度等于单元数的电阻率数组 rho np.full(len(xc), rho_bg, dtypefloat) # 层状: 0~30 米为冲积层, 电阻率 30 ohm-m layer1 zc 30.0 rho[layer1] 30.0 # 倾斜低阻板: 用平面不等式圈定 # 板中心 (x0, y0, z0), 倾角 dip, 厚度 t x0, y0, z0 40.0, 50.0, 45.0 dip np.deg2rad(20.0) normal np.array([np.sin(dip), 0.0, np.cos(dip)]) dist ((xc - x0) * normal[0] (yc - y0) * normal[1] (zc - z0) * normal[2]) slab (np.abs(dist) 8.0) (np.abs(xc - x0) 25.0) (np.abs(zc - 35.0) 30.0) rho[slab] 8.0 # 低阻球体: 半径 12 米, 位于 (30, 60, 70) sphere (xc - 30.0)**2 (yc - 60.0)**2 (zc - 70.0)**2 12.0**2 rho[sphere] 5.0 return rho这里有个容易翻车的点异常体的判断用的是单元中心坐标而不是单元体积覆盖。当异常体边界斜切一个六面体单元时整个单元只有一种电性边界实际上被阶梯化。网格足够细时这种误差可以接受网格很粗时边界阶梯化会造成视电阻率剖面在异常体边界附近出现一个假高值或假低值。如果需要更精确的边界表达就得把网格剖分对齐到地电界面或者提高局部网格密度。3.3 网格与电极的联合生成顺序建模的顺序建议是先定电极位置再生成网格最后给单元赋电阻率。电极坐标要先确定是因为网格节点数组里必须包含电极位置的节点。更稳的做法是直接把电极坐标逐一插入节点列表再对整张表去重排序。如果网格已经生成完电极落到了某个单元内部就要用插值近似电位精度会下降。这个问题在避坑章节会专门展开。正确顺序如下elec_x np.array([10.0, 20.0, 30.0, 40.0, 50.0]) elec_y np.zeros_like(elec_x) elec_z np.zeros_like(elec_x) # 先把电极坐标并进节点数组 all_x np.unique(np.r_[x_nodes, elec_x]) all_y np.unique(np.r_[y_nodes, elec_y]) all_z np.unique(np.r_[z_nodes, elec_z]) # 再建立节点坐标到编号的映射 # 最后用 map 在后续电极电位提取时直接查节点号 node_map {} for idx, (xx, yy, zz) in enumerate(zip(all_x, all_y, all_z)): node_map[(xx, yy, zz)] idx这一步骤虽然简单却是后面减少坑的关键。电极必须精确落在网格节点上这是三维电阻率正演代码能否得到平滑视电阻率剖面的一条分界线。很多论文复现里学生把网格设成均匀粗网格电极随便取了最近节点最后剖面出现锯齿还以为是求解器问题。4. 源项处理与视电阻率计算从电位解到一条可对比的测深曲线矩阵方程 K·Φ F 建好以后真正决定正演结果是否可用的是右端项 F 的构造和视电阻率的换算公式。这一章把两件事分别拆开先把供电电极的点源电流正确放进节点再用温纳装置或三极装置的装置系数把节点电位转换成视电阻率。4.1 点电流源的节点分配点电流源在有限元方程里是 δ 函数直接积分没有意义标准做法是把电流按电极所在节点加入到右端项 F。假设 A 极注入电流 IB 极注入电流 −I那么 F 中对应 A 极节点加 I对应 B 极节点减 I。如果电流值取 1 安培最后求得的电位差再除以电流就能得到视电阻率计算上更干净。电极必须落在网格节点上这点在建模阶段已经保证。但还有一个常见问题当供电极和测量极距离很近而网格不够细时两个节点之间的电位差会严重偏低原因是离散解无法精确表达点源附近的奇异电位。解决手段只能是局部加密网格而不是靠后处理修正。4.2 温纳装置与三极装置的视电阻率公式求得电位分布后沿测线按观测装置提取电位。温纳装置四个电极等间距排列顺序为 A–M–N–B间距均为 a。先解 A 供电、B 为回流电极时的电位 Φ_AB测量 M 与 N 之间的电位差 ΔV Φ_AB(M) − Φ_AB(N)视电阻率按下式计算ρa 2πa · ΔV / I三极装置 AMN 的 B 极在“无穷远”计算区域上把 B 放在远离测线的边界上。装置系数为 K 2π · AM · AN / MN。当 AM a、MN a 时K 4πa。这个系数不要记错很多人第一次算三极视电阻率时少乘了 2。下面这个表格整理了两种装置常用的电极距与装置系数对应关系方便查对。装置电极布置顺序极距关系装置系数 K温纳A–M–N–BAMMNNBa2πa三极A–M–NB无穷远AMaMNa4πa偶极-偶极A–B–M–NABMNaBMnaπa·n(n1)(n2)三维正演里测深数据一般是固定装置、逐步增加极距。极距越大探测深度越大但计算区域也要跟着扩大否则边界反射会影响长极距的视电阻率。4.3 求解稀疏方程组并提取电极电位线性方程组使用 scipy.sparse.linalg.spsolve 做直接求解。节点数在 10 万以内时效率完全够节点数更大时建议用 ILU 预处理的共轭梯度法。以下是完整的求解与视电阻率换算流程。from scipy.sparse.linalg import spsolve def solve_potential(K, nodes, source_list, fix_node): nn K.shape[0] F np.zeros(nn) for node_id, cur in source_list: F[node_id] cur # 远端节点固定电位为 0, 消除诺伊曼边界带来的奇异 K K.tolil() K[fix_node, :] 0.0 K[:, fix_node] 0.0 K[fix_node, fix_node] 1.0 K K.tocsr() F[fix_node] 0.0 phi spsolve(K, F) return phi def wenner_apparent_rho(phi, node_m, node_n, a, I1.0): delta_v phi[node_m] - phi[node_n] return 2.0 * np.pi * a * delta_v / I固定远端节点电位为零这一步非常关键。如果组装矩阵只有诺伊曼边界整个矩阵是奇异的spsolve 会直接报错或者给出满屏 NaN。把一个远离测区的节点钳到零电位不影响任何电位差计算但把矩阵拉回了非奇异区间这在工程上是标准做法。4.4 多极距测深曲线的一次性计算实际测深要在一个测点上跑多个极距比如温纳装置从 a5 米扫到 a100 米。每个极距都要重新生成网格、重新组装矩阵吗如果网格只加密一次大极距的测量电极可能落在较粗的网格区域精度不足每个极距都重建网格又太慢。常见折中方案是对最大极距设计一套网格小极距的电极都落在这套网格的加密区内所有极距共用同一套矩阵。这样做的效率收益很大矩阵组装好一次之后每个极距只需要改右端项和电极节点索引。计算多个极距的测深数据时实际上是一次组装、多次回代求解能省下大概率百分之六十以上的计算时间。# 同一套网格, 不同极距只改源项 for a in [5.0, 10.0, 20.0, 40.0, 80.0]: nodes_a electrode_node_map[(x0 - a, y0, 0)] nodes_m electrode_node_map[(x0 - a/3, y0, 0)] nodes_n electrode_node_map[(x0 a/3, y0, 0)] nodes_b electrode_node_map[(x0 a, y0, 0)] phi_ab solve_potential(K, nodes, [(nodes_a, 1.0), (nodes_b, -1.0)], far_node) rho_a wenner_apparent_rho(phi_ab, nodes_m, nodes_n, a)这里要注意每个极距都要重新计算供电电极位置的电位不能把小极距的电位移用到大极距上因为源位置变了全区域电位分布都变了。这个逻辑理清以后你就能把一条完整的温纳测深曲线串起来了。5. 三维电阻率测深正演的避坑清单五个容易翻车的细节这一章把前面提到的和没提到的常见问题集中整理成五条踩坑记录每一条都是“现象 → 原因 → 解决”的结构。这些问题我基本都在不同类型的三维电阻率模拟代码里见过有些是别人问我的有些是我自己调试时踩过的。5.1 电极不在网格节点上视电阻率剖面出现锯齿现象同一测点、同一模型只把网格稍微平移几个厘米视电阻率就出现可见跳变异常体边界附近尤其严重。原因电极从网格内部强行取最近节点电位源位置和测量位置都引入了位置误差。点电流源附近的电位变化非常剧烈最近节点和真实电极位置之间哪怕差半个单元边长电位误差也足够让视电阻率偏移几个百分点。解决在建模阶段把电极坐标先并入节点数组再去重、排序生成节点编号。后处理阶段想补救只能用电极周围一圈节点的电位做插值但精度远不如直接在节点上建电极。我的习惯是永远先布电极再剖网格不做任何依赖后处理的假设。5.2 计算域边界截断太近长极距视电阻率系统性偏低现象均匀半空间模型下小极距的视电阻率接近背景值大极距时视电阻率逐渐偏离且在边界方向不对称。原因侧向边界被强行设为绝缘或零电位电流在边界处被约束实际等效为一个比真实模型更复杂的边界反射。极距越大边界距离相对越小反射影响越大。解决横向计算域范围至少取最大极距的 5 到 10 倍垂向深度至少取最大极距的 2 到 3 倍。文献里的经验值通常写 5 倍如果追求数据光滑我一般取 8 倍。用 Robin 混合边界可以把域缩小到 3 倍左右但实现时要小心边界项的正负号加错符号会在边界产生电流注入比不加还糟。5.3 远端节点固定电位导致矩阵行/列操作破坏稀疏结构现象网格节点数接近十万量级时固定电极参考点的那几行代码运行特别慢内存占用也比预期高很多。原因K.tolil() 把 CSR 矩阵转成 LIL 格式的代价不低而且对整行、整列清零后原来非零的结构被打乱再转回 CSR 时需要重新压缩。解决如果网格规模不大这种操作可以接受大规模网格建议在组装时就直接跳过参考节点的非对角元素只保留对角元。还可以在组装阶段记录参考节点在 COO 列表里的所有索引组装完成后统一去除。这是血泪经验我自己的第一个版本因为这步操作慢了将近两倍。5.4 低阻异常体参数设置太小求解时出现发散或负电阻率现象模型里设置了电阻率 0.01 Ω·m 的超低阻体求解结果出现巨大负值视电阻率曲线在异常体位置突然下降到负数千。原因电导率极大导致局部矩阵元素量级悬殊直接求解器在数值上失去精度产生虚假振荡。负电阻率在物理上不合理只出现在数值不稳健的状态。解决实际地电模型的电阻率很少低于 1 Ω·m模拟时把低阻体的电阻率控制在合理范围比如 5~20 Ω·m。如果必须模拟极端低阻考虑用对数电阻率参数化或对单元电导率做上限截断。另外还要检查低阻体边界是否与电极单元重叠供电极和低阻体在同一个单元内时局部解会更敏感。5.5 组装矩阵阶段速度极慢纯 Python 循环扛不住现象网格规模只有 50×50×50 个单元组装却花了几个小时的运算时间而实际求解只要几秒。原因每个单元内部 8×8 双重循环可能产生 640 次列表 append单元数超过十万后Python 的列表操作和对象开销成了瓶颈。解决组装阶段尽量按单元批量计算局部矩阵用向量化的方式同时处理一批单元。另一个更实用的技巧是用单元中心电阻率直接插值到高斯点而不是每个高斯点都去做条件判断。如果仍然嫌慢把单元刚度矩阵的计算放到 Numba 的 jit 函数里通常能获得几十倍的速度提升。6. 从均匀半空间解析解开始给三维正演做体检除非你天赋异禀能一次写对所有细节否则每个三维电阻率正演代码都应该先通过均匀半空间模型的解析解验证。地表点电流源在均匀半空间中的电位解析解是 φ ρI/(2πr)r 是观测点到源的水平距离。这个公式在直流电法里是基础得不能再基础的东西但它就是验证程序的“标准砝码”。具体做法是在均匀模型上跑一个温纳测深极距从最小扫到最大计算每个极距的视电阻率误差误差应小于 0.1% 才说明网格和求解器基本可靠。验证代码可以写得很短# 均匀半空间: rho100 ohm-m, 极距 a10 rho_bg 100.0 K_coef 2 * np.pi * 10.0 phi_m mu * I / (2 * np.pi * 10.0) phi_n mu * I / (2 * np.pi * 20.0) # 解析值: rho_a K_coef * (phi_m - phi_n) / I 100.0如果验证曲线在最大极距处偏差超过 2%优先检查计算域边界距离和网格密度而不是去动解析解的公式。这是我每次写完新模块都会做的一步不通过这步就往下加地形或各向异性后面的问题会纠缠成一团乱麻。通过解析验证之后再往模型里加断层、低阻球体、起伏地形。这时还有一个值得做的验证对比不同网格密度下同一模型的视电阻率剖面如果两次加密网格后的曲线差异已经很小说明网格收敛了结果可以对外采信。最后说一个进阶方向各向异性介质的模拟并不难只要把单元电导率 σ 从标量换成 3×3 张量局部刚度矩阵里把 σ 换成张量即可。但各向异性会引入更多参数边界条件对网格方向的依赖更强调试难度会上一个台阶。建议先把各向同性模型的整套流程吃透再去碰各向异性和地形改正。我的习惯是每改一块网格就重新跑一次半空间验证把误差控制住再继续往下走希望帮到你。本文还有配套的精品资源点击获取