简介三维电阻率测深是地球物理勘探的重要方法之一。这份资料面向地球物理勘探、应用数学与计算力学领域的科研人员和研究生提供基于有限元法的三维电阻率测深数值模拟完整复现方案。资源为单一PDF文档大小仅456KB内容包含论文复现的详细推导与可运行Python代码及逐步解释。目前已有75人学习下载。文档系统梳理了点源电场边值与变分问题、六面体单元剖分、三线性插值、单元刚度矩阵组装、高斯积分与雅可比矩阵等关键环节并附带三层水平地层、低阻矿体、多矿体等典型地电模型算例可帮助读者直观理解节点电位求解与地表视电阻率计算的完整流程为复杂地质条件下的电阻率法勘探数据解释提供可靠的数值模拟工具。1. 三维电阻率测深数值模拟到底难在哪有限元法为什么是绕不开的选择“三维电阻率测深数值模拟”听起来只是换了一个维度真正做起来却会让不少人翻车地下不是一个均匀半空间起伏界面、低阻异常体、有限延深断层会把电位场扭曲得面目全非解析公式只够用来做半空间校验。用有限元法做三维数值模拟等于把连续的地下电性结构切成有限个单元把“求电位”变成“解线性方程组”再按电极装置系数换算成视电阻率。这篇笔记就按我复现这类论文的流程写从泊松方程离散、复杂地电结构建模、视电阻率计算到最常踩的坑和校验技巧。照着跑通一个小规模算例你就能往自己的模型上迁移。2. 从泊松方程到有限元离散控制方程、边界条件与刚度矩阵组装2.1 控制方程与边界条件点电源在三维地电模型中怎么描述直流电阻率法的工作方式可以概括为“供电、测电位、反推电阻率”。在三维空间中稳定电流场的电位φ满足∇·(σ∇φ) -I₀δ(r - r_s)其中σ为电导率S/mI₀为供电电流Ar_s为点电源位置。对于块状电阻率模型σ在单元内部为常数在界面上发生跳变而电位本身连续、电流法向分量连续。求解这个方程之后任意测量电极处的电位值就能通过插值得到进一步算视电阻率。实际建模时边界条件需要分开说。地面z0是空气与地下的分界面空气电导率几乎为0电流不能穿出地面所以满足诺依曼条件∂φ/∂n0。模拟区域的左右、前后、底面都是人为截断边界如果距离源足够远电位趋近于0可以用狄利克雷条件φ0。不过截断距离不够远时这种“硬边界”会反射电流测深曲线尾支偏大或偏小所以论文里更常见的是混合边界条件把边界上的电位与法向导数线性组合起来。我在复现时先把边界放到足够远至少几倍最大电极距验证算例再用解析解校验混合边界留到后面做。有限元求解的不是微分方程本身而是它的等效变分形式。稳定电流场对应的泛函是J(φ) ∫_Ω (1/2 σ|∇φ|² - I₀δ(r-r_s)φ) dV极大极小化这个泛函得到的φ就是满足控制方程和自然边界条件的解。这样做的优势在于电阻率突变界面不需要特殊处理边界条件也会在单元积分和边界积分中自动体现。2.2 六面体单元与形函数把连续电位变成节点未知量有限元的第一步是把连续区域划分成单元。这里常用的是规则六面体单元原因有三个一是网格坐标和节点编号可以用三重循环生成便于实现二是复杂地电结构通过给不同单元赋不同的电阻率值来表达不需要动网格三是在测深方向上电极附近网格可以适当加密远处逐渐拉伸减少未知量。线性六面体单元的每个单元有8个节点节点编号顺序决定了形函数和几何映射。在标准局部坐标系(ξ,η,ζ)下8个节点的局部坐标都是±1形函数为N_i (1ξξ_i)(1ηη_i)(1ζζ_i)/8这个式子对每个方向独立所以三维单元刚度矩阵可以拆成三个方向的一维组合也可以通过数值积分一次性算。局部坐标到全局坐标的映射由节点坐标线性插值完成。单元刚度矩阵的元素为K_{e,ij} ∫_e σ_e ∇N_i·∇N_j dV其中σ_e是单元电导率。这里要注意如果σ_e在单元内不变可以提到积分号外但形函数梯度仍随局部坐标变化。三维情况下用手推解析式很繁琐通常用高斯积分。每个方向取2个高斯点三维共8个积分点积分点坐标为±1/√3权重为1。2.3 组装刚度矩阵的代码把单元贡献装进稀疏矩阵下面是我在复现时用的最小实现。第一步定义网格与模型import numpy as np import scipy.sparse as sp import scipy.sparse.linalg as spla nx, ny, nz 21, 21, 16 # 网格节点数 x np.linspace(-80, 80, nx) y np.linspace(-80, 80, ny) z np.linspace(0, 30, nz) # z 向下为正 dx np.diff(x) dy np.diff(y) dz np.diff(z) nnode nx * ny * nz nelem (nx - 1) * (ny - 1) * (nz - 1) def node_id(i, j, k): return i j * nx k * nx * ny def fine_nodes(ex, ey, ez): 返回单元 (ex,ey,ez) 的8个全局节点编号顺序与局部坐标一致 i, j, k ex, ey, ez return [ node_id(i, j, k), node_id(i1, j, k), node_id(i1, j1, k), node_id(i, j1, k), node_id(i, j, k1), node_id(i1, j, k1), node_id(i1, j1, k1), node_id(i, j1, k1), ]网格范围x、y都是-80到80z从0到30。对一个小规模测深算例最大电极距在30m以内边界距离超过3倍基本能压低边界效应。该网格共21×21×167056个节点、6000个单元未知量不大能说明完整流程。node_id按“x方向最快、y方向次之、z方向最慢”编号这决定了组装时非零元素的分布对稀疏矩阵的迭代求解性能有影响。第二步计算单元刚度矩阵def hex8_stiffness(sigma_e, dx, dy, dz): 规则六面体单元刚度矩阵sigma_e 为单元电导率 gp 1.0 / np.sqrt(3.0) gauss [ (-gp, -gp, -gp, 1.0), ( gp, -gp, -gp, 1.0), ( gp, gp, -gp, 1.0), (-gp, gp, -gp, 1.0), (-gp, -gp, gp, 1.0), ( gp, -gp, gp, 1.0), ( gp, gp, gp, 1.0), (-gp, gp, gp, 1.0), ] Ke np.zeros((8, 8)) vol dx * dy * dz for (xi, eta, ze, w) in gauss: grad np.zeros((8, 3)) for i in range(8): # 局部坐标顺序: x方向 - y方向 - z方向 sx (i % 2) * 2 - 1 sy ((i // 2) % 2) * 2 - 1 sz (i // 4) * 2 - 1 dN_dxi sx / 8.0 * (1 sy * eta) * (1 sz * ze) dN_deta sy / 8.0 * (1 sx * xi) * (1 sz * ze) dN_dze sz / 8.0 * (1 sx * xi) * (1 sy * eta) grad[i] [dN_dxi / dx, dN_deta / dy, dN_dze / dz] Ke w * sigma_e * (grad grad.T) * vol return Ke规则单元的雅可比行列式为dx·dy·dz所以形函数对全局坐标的导数就是局部导数除以对应边长。高斯点取±1/√3每个方向两个点三维组合后有8个积分点权重都是1。这里的关键参数是dx、dy、dz如果某个方向网格尺寸变化过大比如近电极处0.5m远处直接跳到10m单元刚度矩阵数值差异会拉大全局矩阵条件数变差。第三步组装全局矩阵并施加狄利克雷边界A_mat sp.lil_matrix((nnode, nnode)) rhs np.zeros(nnode) for e in range(nelem): ex e % (nx - 1) ey (e // (nx - 1)) % (ny - 1) ez e // ((nx - 1) * (ny - 1)) sigma_e 1.0 / rho_model[ex, ey, ez] Ke hex8_stiffness(sigma_e, dx[ex], dy[ey], dz[ez]) nodes fine_nodes(ex, ey, ez) for a in range(8): ia nodes[a] for b in range(8): ib nodes[b] A_mat[ia, ib] Ke[a, b] # 狄利克雷边界把最外层节点电位置零 boundary [] for k in range(nz): for j in range(ny): boundary.extend([node_id(0, j, k), node_id(nx - 1, j, k)]) for k in range(nz): for i in range(nx): boundary.extend([node_id(i, 0, k), node_id(i, ny - 1, k)]) for i in range(nx): for j in range(ny): boundary.append(node_id(i, j, nz - 1)) for bid in boundary: A_mat[bid, :] 0 A_mat[:, bid] 0 A_mat[bid, bid] 1.0 rhs[bid] 0.0 A_mat A_mat.tocsc()这里用LIL矩阵逐个累加适合第一次写通逻辑。网格稍大后建议改成COO格式先记录所有非零坐标和值再用sp.coo_matrix((data,(row,col)))一步组装速度会快一个量级。边界处理完还要用tocsc()转换因为后面直接求解或迭代都更高效。3. 复杂地电结构建模电阻率模型到网格单元的映射与源项处理3.1 网格剖分策略起伏界面、异常体和测深范围怎么取舍复杂地电结构建模的难点不在于“画一个多边形”而在于把真实地质体的尺度、位置和电阻率值映射到有限元单元上。规则网格下异常体边界只能用阶梯状近似这是绕不开的代价。我在复现论文时通常先把区域分成背景层和异常体背景层用水平分层或起伏界面描述异常体用三维掩膜数组描述。网格剖分要兼顾三个尺度电极距、异常体尺寸和围岩体积。对测深来说浅部电极距小、电位梯度大所以地表附近网格要密深部电位变化平缓网格可以逐渐拉伸。但网格拉伸比不宜超过1.5倍否则单元形状畸变会让刚度矩阵病态。x和y方向至少覆盖最大电极距的3到5倍z方向要延伸到最深探测深度以下。若模型里存在低阻异常体它对电流有“吸引”作用边界距离要取更大否则低阻体会把电流引向边界尾支视电阻率明显偏低。还有一个取舍是要不要为起伏界面做贴体网格论文里有时用坐标变换生成贴体网格但工程落地时更常见的是用加密的规则网格去逼近界面。只要界面附近单元尺寸足够小视电阻率曲线对界面位置的变化就足够敏感。我的习惯是先做两组网格一组加密一倍对比结果不超过3%再继续算测深。3.2 从电阻率模型生成单元电导率一个最小代码假设背景为100Ω·m表层5m内为20Ω·m的覆盖层在坐标(-15,15)×(-15,15)×(10,18)范围内有一个5Ω·m的低阻异常体。用规则网格表达这个模型只需一步布尔运算# 单元中心坐标 xc (x[:-1] x[1:]) / 2.0 yc (y[:-1] y[1:]) / 2.0 zc (z[:-1] z[1:]) / 2.0 Xc, Yc, Zc np.meshgrid(xc, yc, zc, indexingij) # 背景 100 Ω·m, 浅层覆盖 20 Ω·m rho_model np.ones((nx - 1, ny - 1, nz - 1)) * 100.0 rho_model[:, :, :5] 20.0 # 低阻异常体 mask ( (Xc -15) (Xc 15) (Yc -15) (Yc 15) (Zc 10) (Zc 18) ) rho_model[mask] 5.0这段代码中meshgrid使用indexingij是为了让数组维度顺序与节点编号一致第一个维度对应x第二个对应y第三个对应z。rho_model的形状是(nx-1, ny-1, nz-1)与fine_nodes里的单元循环顺序完全对应避免组装时单元电导率取错。对于起伏界面可以用深度数组topo(xc, yc)描述每个单元中心处界面的z坐标然后比较Zc与topo[Xc, Yc]即可给不同层位赋不同电阻率。这里要注意的是topo数组需要用插值生成并且界面处不要出现锯齿过大的突变否则梯度计算在界面附近容易失真。3.3 源项分配点电源的奇异性和最近节点加载点电源在理论上是δ函数电位在源点处趋于无穷直接数值求解会把全部电流压到一个节点上导致周围节点的电位被污染。论文里常见做法是奇异电位法把电位拆成背景均匀半空间解析解和剩余电位只对剩余电位做有限元求解。这个方案实现起来比较复杂。我在快速验证时先用一个工程近似——“源点体积分配”把供电电流按形函数权重分配到包含源点的单元节点上def distribute_source(sx, sy, sz, current): 把点电源 current 分配到源点所在单元的8个节点 ex min(np.searchsorted(x, sx) - 1, nx - 2) ey min(np.searchsorted(y, sy) - 1, ny - 2) ez min(np.searchsorted(z, sz) - 1, nz - 2) ex max(ex, 0); ey max(ey, 0); ez max(ez, 0) # 源点在单元内的局部坐标 xi (sx - x[ex]) / dx[ex] * 2.0 - 1.0 eta (sy - y[ey]) / dy[ey] * 2.0 - 1.0 zeta (sz - z[ez]) / dz[ez] * 2.0 - 1.0 nodes fine_nodes(ex, ey, ez) for i in range(8): sx_i (i % 2) * 2 - 1 sy_i ((i // 2) % 2) * 2 - 1 sz_i (i // 4) * 2 - 1 N (1 xi * sx_i) * (1 eta * sy_i) * (1 zeta * sz_i) / 8.0 rhs[nodes[i]] current * Nnp.searchsorted返回第一个大于等于sx的位置减1后就是源点左侧的网格索引。对于正好落在节点或网格边界上的源点部分形函数权重为0分配会自动处理。注意这种体积分配法只能缓解奇异性不能完全消除。做高精度论文复现时最终还是要换奇异电位法或者把源点附近网格加密到足够细让视电阻率结果对源点邻域不再敏感。当供电是双极源时A点加IB点加-I。不要忘了B点的负号否则电位差会差出一倍。测量电极M、N通常也在节点之间取值时要对电位做三线性插值不能用最近节点代替否则浅部小极距测深会看到明显台阶。4. 视电阻率计算电极装置、求解流程与测深曲线输出4.1 电极装置与视电阻率公式从电位差到电阻率求解得到节点电位后视电阻率不是直接读某个电极的电位而是通过装置系数换算。地表均匀半空间中点电源的电位解析解是φ I₀ρ / (2πr)这里用2π而不是4π因为电流只在地下半空间扩散。对于任意四极装置供电电极A、B测量电极M、N视电阻率公式为ρ_a K · (φ_M - φ_N) / I₀其中K 2π / (1/AM - 1/AN - 1/BM 1/BN)温纳装置的四个电极等间距排列AMMNNBa代入后K2πa。三极装置中B极放到“无穷远”理论上K4πa当A-Ma且M-Na时但实际模拟时不可能真放无穷远我会把B放在距测点超过20倍a的位置然后仍用通用公式算K这样即使B不够远装置系数也能把几何影响考虑进去。装置系数的物理意义是“把电位差折算成均匀半空间视电阻率”的几何因子。只要电极坐标确定K就能精确计算不需要管地下模型。所以代码里我不会写死温纳或三极的公式而是先算各电极间距离再算Kdef electrode_factor(A, B, M, N): A,B为供电电极坐标M,N为测量电极坐标返回装置系数K AM np.linalg.norm(np.array(A) - np.array(M)) AN np.linalg.norm(np.array(A) - np.array(N)) BM np.linalg.norm(np.array(B) - np.array(M)) BN np.linalg.norm(np.array(B) - np.array(N)) return 2.0 * np.pi / (1.0/AM - 1.0/AN - 1.0/BM 1.0/BN)这里所有电极默认为地表坐标z0。如果后续要在坑道或井中做三维电阻率测深电极坐标会包含z分量但通用公式依然成立只需要把距离改成三维空间距离。4.2 求解主流程加载模型、组装、求解、计算视电阻率我把求解流程封装成一个函数输入是模型数组和电极坐标输出是各测量电极的电位。这样测深曲线就是循环改变电极坐标后反复调用同一个函数def solve_potential(rho_model, src_list): src_list: [(x,y,z,current), ...] # 组装复用前面的代码 A_mat sp.lil_matrix((nnode, nnode)) rhs np.zeros(nnode) for e in range(nelem): ex e % (nx - 1) ey (e // (nx - 1)) % (ny - 1) ez e // ((nx - 1) * (ny - 1)) sigma_e 1.0 / rho_model[ex, ey, ez] Ke hex8_stiffness(sigma_e, dx[ex], dy[ey], dz[ez]) nodes fine_nodes(ex, ey, ez) for a in range(8): ia nodes[a] for b in range(8): ib nodes[b] A_mat[ia, ib] Ke[a, b] # 边界处理 for bid in boundary_nodes: A_mat[bid, :] 0 A_mat[:, bid] 0 A_mat[bid, bid] 1.0 # 源项 for sx, sy, sz, cur in src_list: distribute_source(sx, sy, sz, cur) A_mat A_mat.tocsc() phi spla.spsolve(A_mat, rhs) return phi注意一个细节施加狄利克雷边界时我先把整个矩阵行清空再把对角线设1但没有清rhs。如果某个边界节点上恰好被分配了源项源项会被边界方程吞掉。因此边界处理要在源项分配之前完成或者在对角线置1后把rhs对应位置置0。上面代码的顺序是先边界后源项是正确的。求解器选择上spla.spsolve对几万未知量的问题是秒出的。当网格加密到几十万节点时直接求解会吃内存此时改用cg迭代phi, info spla.cg(A_mat, rhs, rtol1e-8, maxiter1000)cg要求矩阵对称正定。电阻率有限元矩阵在狄利克雷边界下是对称正定的所以可以放心用。maxiter不要小于500否则大网格下还没收敛就退出测深曲线会抖。4.3 测深曲线与拟断面图结果怎么组织和观察以温纳测深为例A、B、M、N沿x轴排列电极距a从2m递增到30m每个a对应一条测深数据def wenner_measure(rho_model, center_x0, a_listnp.arange(2, 31, 2)): rho_a [] for a in a_list: A (center_x - 1.5*a, 0, 0) B (center_x 1.5*a, 0, 0) M (center_x - 0.5*a, 0, 0) N (center_x 0.5*a, 0, 0) phi solve_potential(rho_model, [(A[0], A[1], A[2], 1.0), (B[0], B[1], B[2], -1.0)]) phi_m sample_potential(phi, M[0], M[1], M[2]) phi_n sample_potential(phi, N[0], N[1], N[2]) K electrode_factor(A, B, M, N) rho_a.append(K * (phi_m - phi_n)) return np.array(a_list), np.array(rho_a)sample_potential用三线性插值。如果要画拟断面图就固定一套电极距沿x方向滑动测点把结果填成二维数组再用plt.contourf绘图。拟断面图能直观看到低阻异常体引起的“V形”视电阻率低值区这是三维测深数值模拟最常见的产出图件。需要注意视电阻率值对网格边界和源点分配方式非常敏感。同一模型、同一电极装置用不同网格剖分得到的测深曲线首支可能差5%以上。所以论文复现中必须先把代码放在均匀半空间模型上做校验确认装置系数和边界处理都正确后再开始加复杂结构。5. 三维电阻率测深模拟的常见坑收敛、边界、网格与源项排查5.1 现象电位出现负值或锯齿振荡原因源点附近网格太粗点电源奇异性没有被体积分配吸收造成源点周围电位梯度异常另一个常见原因是单元刚度矩阵计算错误尤其是局部节点顺序和形函数导数顺序不对导致矩阵不对称。解决先把模型设为均匀半空间单独查看距离源点几个单元处的电位是否与解析解一致。如果曲线呈锯齿状检查fine_nodes和hex8_stiffness里的节点顺序是否一一对应。我的经验是写一个单元级测试任意取一个单元计算所有节点电位都相等时的Ke行和理论上每行之和应为0若不为0说明形函数导数积分有误。5.2 现象网格加密后结果反而偏离解析解原因边界和源点处理不一致。很多人加密网格时只加密了电极附近x、y外边界距离没有跟着扩大或者网格过度拉伸。加密后单元体积变小源点体积分配的范围也变小如果网格加密方向与源点分配不匹配结果自然漂移。解决加密网格时必须同时检查边界到最近供电点的最小距离至少保持在最大电极距的3倍以上。另外加密前后对比均匀半空间模型的视电阻率误差应小于2%。如果加密后结果偏离更大要优先怀疑边界条件而不是有限元本身。5.3 现象测深曲线首支或尾支畸变原因首支小电极距畸变通常是源点奇异性影响因为测量电极离供电点太近源点体积分配造成的误差直接落到测量电位上尾支畸变通常是截断边界太近或网格深度不够。解决首支畸变可在源点附近加密网格把最近电极距对应的单元尺寸控制在1/5以下尾支畸变则把z方向延伸加深并加宽x、y范围。还有一个技巧是用对称四极装置时先分别计算A和B单独供电的电位再加权叠加这样能看清是哪一侧的源造成了尾支异常。不要一上来就调装置系数那往往不是问题根源。5.4 现象供电点和接收点太近时的数值爆表原因供电点在A、B接收点在M、N当A和M重合或距离小于一个网格间距时电位差接近两个巨大数值相减浮点误差被放大。此外源点分配把电流加到有限几个节点上离源越近的节点电位越大插值误差也越大。解决在实际测深设计中很少让供电测量电极完全重合但数值模拟中为了扫小电极距经常出现。我一般把最小电极距设为源点周围网格尺寸的2倍以上若必须模拟微小距则采用奇异电位法并单独对源点周围网格加密。另一个兜底办法是改用更高精度的插值三线性插值而不是最近节点来取测量电极电位。5.5 现象CG迭代不收敛或收敛极慢原因网格尺度跨越太大比如近电极处0.2m远处20m矩阵条件数可能到1e8以上或者边界节点没有置1导致矩阵奇异。解决先用np.linalg.cond检查均匀模型的矩阵条件数超过1e10就要考虑网格拉伸比。规则网格下我建议各方向最大拉伸比控制在1.3以内宁肯多加一些单元也不要让条件数失控。迭代时使用rtol1e-10配合maxiter2000实际运行时观察每步残差若前几百步残差下降缓慢多半是网格问题。6. 验证与提速用解析解校验、用对称性压缩计算量6.1 均匀半空间解析解校验最小验证用例无论模型多复杂第一步都该跑均匀半空间。取电阻率100Ω·m电极距a10m的温纳装置理论视电阻率必须等于100Ω·m。我习惯把供电电流设为1A这样装置系数乘以电位差就是视电阻率省去量级转换。如果算出来是100±3%以内说明边界距离、源项分配和装置系数都基本正确。超过5%就先别往下做。6.2 利用对称性压缩计算量三维电阻率模拟中有大量对称性可用。如果地电结构关于x轴对称且测深点位于对称面上那么供电点A和B对称时电位也对称。此时只需组装一次矩阵对两个右侧向量求解甚至可以把源项合并求解减少一半计算。模型完全对称时我还会把y方向网格减半只模拟半空间并把y0对称面设为诺依曼边界自由度直接少一半。网格未知量从几万降到几千迭代求解会快很多。但要注意对称模型必须同时满足几何和电性对称。一旦加入倾向、各向异性或倾斜异常体对称就破坏了不能硬用。6.3 我现在的习惯和踩过的坑以前我为了追求省事把源点直接加到最近节点结果均匀半空间校验就差了10%。后面改成体积分配再把地表附近网格加密误差才降到1%。现在我的固定流程是画网格、建均匀模型、解析解校验、加背景层、加异常体、跑测深曲线、换网格再跑一遍对比。每次改动只动一个变量防止多个坑叠在一起。三维电阻率测深数值模拟的论文复现难点不在“会写公式”而在“知道每个参数动了以后结果会往哪偏”。把边界、源点、网格三件事管住后面的复杂地电结构建模基本就是体力活。希望这套流程能帮你少走一点弯路。希望帮到你。本文还有配套的精品资源点击获取