简介本资源是一套面向地球物理勘探专业学生、科研人员及工程技术人员的电法勘探正演建模实践代码包聚焦MATLAB环境下电阻率法正演模拟的核心实现解决地质参数建模、电场响应预测与数据可视化等关键问题。压缩包共含4个文件3个MATLAB脚本1幅TIFF插值图总大小仅23KB轻量实用firstwork.m为主控脚本统筹流程DataFFT.m实现基于快速傅里叶变换的高效正演计算elecltick.m负责电极布设与电流分布模拟Lagrange_Interpolation_6_150dpi.tif则直观呈现拉格朗日插值生成的二维电阻率分布效果。已有344人学习下载内容紧扣AKQH二维电位法模型涵盖MATLAB基础编程、泊松方程数值求解有限差分、滤波插值等实操环节提供可直接运行、参数可调的完整正演工作流是理解正演原理与开展反演前仿真实验的优质入门材料。1. 电法勘探正演不是“画图填数”而是用物理方程驱动的数值实验当你在野外布设完电极、记录下视电阻率曲线后真正决定解释精度的往往不是后期反演软件的参数调优而是正演模型是否真实复现了地下电性结构对电流场的响应。mt_电法勘探编程_电法勘探正演_这个标题指向的是一套基于麦克斯韦方程组低频近似即准静态近似构建的、可编程实现的直流/频率域电法正演计算流程——它不依赖商业软件界面而是通过离散化控制方程、组装稀疏矩阵、求解大型线性系统输出理论电位分布或视电阻率响应。这类编程实践常见于高校地电实验室、物探院所算法验证、以及国产地球物理软件底层模块开发。适合有数值方法基础如有限差分、有限元、熟悉 Python 或 MATLAB 的地球物理方向工程师对仅会调用 GeoStudio 或 Res2Dinv 的用户而言这是从“操作员”转向“建模者”的关键跃迁。正演不是黑箱输出它是你对地下构造假设的第一次数学检验如果正演结果连均匀半空间都拟合不了反演再快也无意义。2. 用有限差分法在规则网格上实现直流电法正演的最小可运行代码2.1 为什么选有限差分而非有限元——从地质建模需求出发的选型逻辑电法正演方法中有限元法FEM对复杂地形和任意边界适应性强但网格生成与刚度矩阵组装开销大对初学者调试不友好而有限差分法FDM在规则矩形网格上直接离散拉普拉斯方程矩阵结构高度稀疏且规律性强便于手动推导系数、快速验证物理逻辑。对于mt_电法勘探编程_电法勘探正演_这类教学级或原型验证场景FDM 是更务实的选择它能清晰暴露电流守恒、边界条件施加、源项处理等核心物理约束且单次正演耗时可控千级网格点可在毫秒级完成。注意此处“mt”并非指“大地电磁法”magnetotellurics而是项目命名中的标识符与电法正演本身无物理关联——避免与MT数据处理混淆。2.2 拉普拉斯方程的离散化从连续场到稀疏矩阵直流电法在准静态假设下满足泊松方程∇·(σ∇φ) -Iδ(r−rₛ)其中 σ 为电导率S/mφ 为电位VI 为供电电流Aδ 为狄拉克函数。对二维水平层状介质忽略地形在均匀网格步长 dx, dy上采用中心差分节点 (i,j) 处的离散形式为σᵢⱼ(φᵢ₊₁,ⱼ − 2φᵢ,ⱼ φᵢ₋₁,ⱼ)/dx² σᵢⱼ(φᵢ,ⱼ₊₁ − 2φᵢ,ⱼ φᵢ,ⱼ₋₁)/dy² -Iδᵢ,ⱼ整理后得线性系统Aφ b其中 A 为 N×N 稀疏矩阵N 为网格点总数b 为源项向量仅供电电极位置非零。关键在于电导率 σ 随位置变化时相邻节点间需采用调和平均处理界面电导如 σᵢ₊₀.₅,ⱼ 2σᵢ,ⱼσᵢ₊₁,ⱼ/(σᵢ,ⱼσᵢ₊₁,ⱼ)否则电流法向连续性被破坏。2.3 Python 实现从网格定义到电位求解的完整流程import numpy as np from scipy.sparse import diags, csr_matrix from scipy.sparse.linalg import spsolve def build_fdm_matrix(nx, ny, dx, dy, sigma): 构建二维有限差分稀疏矩阵 A sigma: (ny, nx) 形状的电导率数组 返回: csr_matrix 格式的 A 和源项向量 b N nx * ny # 初始化对角线存储 main_diag np.zeros(N) off_diag_x np.zeros(N-1) # x方向相邻 off_diag_y np.zeros(N-nx) # y方向相邻 for i in range(ny): for j in range(nx): idx i * nx j # 中心点系数 coeff 0.0 # x方向左邻j0 if j 0: sigma_harm 2 * sigma[i,j] * sigma[i,j-1] / (sigma[i,j] sigma[i,j-1]) coeff sigma_harm / dx**2 off_diag_x[idx-1] -sigma_harm / dx**2 # x方向右邻jnx-1 if j nx-1: sigma_harm 2 * sigma[i,j] * sigma[i,j1] / (sigma[i,j] sigma[i,j1]) coeff sigma_harm / dx**2 off_diag_x[idx] -sigma_harm / dx**2 # y方向下邻i0 if i 0: sigma_harm 2 * sigma[i,j] * sigma[i-1,j] / (sigma[i,j] sigma[i-1,j]) coeff sigma_harm / dy**2 off_diag_y[idx-nx] -sigma_harm / dy**2 # y方向上邻iny-1 if i ny-1: sigma_harm 2 * sigma[i,j] * sigma[i1,j] / (sigma[i,j] sigma[i1,j]) coeff sigma_harm / dy**2 off_diag_y[idx] -sigma_harm / dy**2 main_diag[idx] coeff # 组装稀疏矩阵 diagonals [main_diag, off_diag_x, off_diag_x, off_diag_y, off_diag_y] offsets [0, -1, 1, -nx, nx] A diags(diagonals, offsets, shape(N,N), formatcsr) # 构造源项 b供电电极在中心回路电极在无穷远Dirichlet边界 b np.zeros(N) src_i, src_j ny//2, nx//2 # 中心供电 b[src_i * nx src_j] 1.0 # 单位电流源 return A, b # 参数设置 nx, ny 100, 50 dx dy 1.0 sigma np.ones((ny, nx)) * 0.01 # 均匀半空间电导率 0.01 S/m A, b build_fdm_matrix(nx, ny, dx, dy, sigma) phi spsolve(A, b).reshape((ny, nx)) # 输出电位分布用于后续视电阻率计算 print(f电位矩阵形状: {phi.shape}, 最小值: {phi.min():.6f}, 最大值: {phi.max():.6f})提示此代码未显式施加边界条件实际应用中需将边界行/列设为零电位Dirichlet或添加镜像电荷Neumann否则矩阵奇异。spsolve要求 A 必须可逆因此必须处理边界。该代码输出phi是归一化电位场单位为 V/A。其物理意义是当供电电流为 1A 时各网格点的电位响应。后续计算视电阻率时需从中提取测量电极对M,N处的电位差 Δφ并代入公式 ρₐ k·Δφ/I其中 k 为装置系数如温纳装置 k2πaa 为电极间距。3. 从电位场到视电阻率装置响应计算与常见陷阱排查3.1 四极装置响应的三种计算路径及适用场景正演输出的是电位场 φ(x,y)但野外实测的是视电阻率 ρₐ二者通过装置几何关系转换。常见路径有方法计算方式优点缺点适用阶段解析插值法在 φ 矩阵中定位 M/N 电极坐标双线性插值得 φₘ, φₙ计算 Δφ实现简单精度依赖网格密度电极位置非网格点时插值误差大快速验证、粗网格测试有限差分卷积法将电极对视为局部电流源在 φ 上做加权平均如 M 点取 3×3 区域均值抑制网格阶梯效应需经验确定窗口大小中等精度要求子网格积分法对每个电极覆盖区域如半径 r 的圆在 φ 上做面积分再除以面积物理意义最严谨抗网格敏感计算开销增加 3~5 倍发表级结果、高精度对比对于mt_电法勘探编程_电法勘探正演_这类编程实践推荐从解析插值法起步待流程跑通后再升级。3.2 温纳装置正演响应的完整计算示例def calculate_wenner_apparent_resistivity(phi, dx, a, src_pos(0,0)): 计算温纳装置AMNB电极距 a的视电阻率 phi: (ny, nx) 电位矩阵 a: 电极间距米 src_pos: 供电电极 A 坐标全局坐标系原点 返回: ρ_a 数组长度 测点数 ny, nx phi.shape # 温纳装置电极坐标A-M-N-B 共线间距均为 a # 设 A 在 (0,0)则 M 在 (a,0), N 在 (2a,0), B 在 (3a,0) # 注意phi 网格原点在左下角y 向上为正需坐标映射 x_coords np.array([0, a, 2*a, 3*a]) y_coord 0.0 # 地表y0 对应 phi 最上一行索引 ny-1 # 插值获取各电极电位 phi_vals [] for x in x_coords: # 将物理坐标 x 映射到网格索引 jx 方向 j_float x / dx j_low, j_high int(np.floor(j_float)), int(np.ceil(j_float)) if j_low 0: j_low, j_high 0, 0 if j_high nx: j_high nx-1 # 双线性插值简化为 x 方向线性插值y 固定为地表行 row_idx ny - 1 # 地表对应最后一行 if j_low j_high: phi_elec phi[row_idx, j_low] else: weight j_float - j_low phi_elec phi[row_idx, j_low] * (1-weight) phi[row_idx, j_high] * weight phi_vals.append(phi_elec) phi_A, phi_M, phi_N, phi_B phi_vals delta_phi phi_M - phi_N # M-N 电位差 k_factor 2 * np.pi * a # 温纳装置系数 rho_a k_factor * delta_phi # 单位Ω·m return rho_a # 示例调用 a_spacing 5.0 # 电极距 5 米 rho_obs calculate_wenner_apparent_resistivity(phi, dx, a_spacing) print(f温纳装置首点视电阻率: {rho_obs:.4f} Ω·m)参数说明dx网格步长米必须与建模时一致a_spacing电极间距直接影响装置系数k_factorsrc_pos供电电极全局坐标用于确定 M/N 位置此处简化为原点phi由spsolve得到的电位矩阵已归一化I1A。注意若地下存在高阻薄层电位场在层界面处梯度剧变此时双线性插值会低估 Δφ导致 ρₐ 偏低。此时应启用子网格积分法或加密网格dx 减小至 a/5 以下。3.3 三类高频报错及其物理根源报错现象控制台输出物理原因排查指令LinAlgError: Matrix is singularspsolve报错边界未处理拉普拉斯算子无零空间约束检查A.diagonal().min()是否接近 0强制设置边界行/列为 0 后重建 Arho_a 为负值print(rho_obs)输出负数电位插值顺序错误如 M/N 坐标颠倒或供电方向与装置定义相反打印phi_M, phi_N值确认phi_M phi_N电流从 A→BM 更近 Aρₐ 曲线在浅部剧烈震荡绘图显示高频毛刺网格过粗dx a/3无法分辨电极尺度下的电位变化执行np.gradient(phi, axis1)查看 x 方向电位梯度若出现 1e4 的突变值需减小 dx这些错误均非代码 bug而是物理建模失配的信号——正演的本质是让数学模型忠实地成为地质假设的“数字孪生”。4. 多层介质正演加速利用对称性压缩矩阵规模与内存占用4.1 层状介质的天然对称性如何削减 75% 计算量当模型为水平层状σ 仅随深度 z 变化时电位场 φ 在 x 方向呈偶对称关于供电轴y 方向呈奇对称若考虑地形则失效但地表平坦时成立。这意味着只需计算半个剖面x≥0另一半通过镜像获得。更进一步若采用柱坐标系r,z则问题退化为二维轴对称未知数数量从 Nₓ×N_z 降至 N_r×N_z典型减少 3~10 倍。mt_电法勘探编程_电法勘探正演_中若涉及多层模型如 5 层必须启用此优化否则 200×100 网格将生成 20000×20000 矩阵内存超限。4.2 柱坐标有限差分矩阵的构建要点在 (r,z) 坐标下拉普拉斯方程变为∇·(σ∇φ) (1/r)∂/∂r(rσ∂φ/∂r) ∂/∂z(σ∂φ/∂z) -Iδ(r)δ(z-zₛ)离散时需注意r0 处的奇异性第一行r0的差分格式需特殊处理避免 1/r 除零网格非均匀r 方向常采用对数间隔rᵢ r₀·qⁱ以提高浅部分辨率系数含 r矩阵非对称A[i,j] ≠ A[j,i]但仍是稀疏的。def build_cylindrical_matrix(nr, nz, r_edges, z_edges, sigma_rz): 构建柱坐标系下稀疏矩阵nr: r向节点数, nz: z向节点数 r_edges: (nr1,) 数组r方向单元边界 sigma_rz: (nz, nr) 电导率对应单元中心 N nr * nz data, rows, cols [], [], [] for i in range(nz): # z索引 for j in range(nr): # r索引j0为r0轴 idx i * nr j r_center 0.5 * (r_edges[j] r_edges[j1]) dr r_edges[j1] - r_edges[j] dz z_edges[i1] - z_edges[i] # r0 处j0使用半圆对称仅保留 ∂²φ/∂z² 项 if j 0: # 中心点-σ∂²φ/∂z² ≈ -σ*(φ[i1,0]-2φ[i,0]φ[i-1,0])/dz² # 系数2σ/dz²对角-σ/dz²上下 if i 0: rows.append(idx); cols.append(idx-1); data.append(-sigma_rz[i,j]/dz**2) rows.append(idx); cols.append(idx); data.append(2*sigma_rz[i,j]/dz**2) if i nz-1: rows.append(idx); cols.append(idx1); data.append(-sigma_rz[i,j]/dz**2) else: # 一般点含 r 项 r_left r_edges[j] r_right r_edges[j1] sigma_r_avg 2*sigma_rz[i,j]*sigma_rz[i,j-1]/(sigma_rz[i,j]sigma_rz[i,j-1]) if j0 else 0 # r方向项系数(1/r)(∂/∂r)(rσ∂φ/∂r) 离散为... # 此处省略详细推导核心是引入 r_center 加权 # 实际代码需按标准文献如Zhdanov, 2002实现 pass return csr_matrix((data, (rows, cols)), shape(N,N))提示柱坐标实现比直角坐标复杂 3 倍但对层状模型是刚需。若仅需快速验证可先用直角坐标对称边界条件设置 x0 区域 φφ(|x|,z)近似牺牲 10% 精度换取 50% 开发时间。4.3 内存与速度的量化平衡网格密度选择指南场景推荐 dx/dz网格点数典型单次正演耗时i7-11800H适用性教学演示均匀半空间2.0 m50×25 125010 ms✅ 快速验证流程三层模型精度验证0.5 m200×100 20000~120 ms✅ 论文图表生成含薄层厚度2m0.1 m1000×500 5000003 s需迭代法⚠️ 改用 AMG 预条件子或 GPU 加速关键原则dx 应小于最薄目标层厚度的 1/3且小于最小电极距 a 的 1/5。例如探测 1m 厚煤层a2m则 dx ≤ 0.4m若用温纳装置 a10m则 dx ≤ 2m —— 这就是为何野外大间距测量无法识别薄层正演同样受此物理限制。5. 验证正演正确性的三个硬性指标从理论解到野外数据拟合5.1 均匀半空间解析解的自动比对脚本任何电法正演代码上线前必须通过均匀半空间σ₀的解析解验证。直流点源在半空间中电位为 φ(r,z) I/(2πσ₀√(r²z²))。编写自动比对脚本计算数值解与解析解的相对误差def validate_uniform_halfspace(): # 设置均匀模型 sigma_uniform 0.01 nx, ny 200, 100 dx dy 0.5 sigma np.full((ny, nx), sigma_uniform) # 数值求解 A, b build_fdm_matrix(nx, ny, dx, dy, sigma) phi_num spsolve(A, b).reshape((ny, nx)) # 解析解供电在 (0,0)地表 z0 r_grid, z_grid np.meshgrid( np.arange(nx)*dx - (nx//2)*dx, # x 从 -50 到 49.5 np.arange(ny)*dy # z 从 0 到 49.5 ) phi_ana 1.0 / (2 * np.pi * sigma_uniform * np.sqrt(r_grid**2 z_grid**2 1e-12)) # 取地表z0和浅部z2m两层比对 err_surf np.abs(phi_num[-1,:] - phi_ana[-1,:]) / (phi_ana[-1,:]1e-12) err_shallow np.abs(phi_num[4,:] - phi_ana[4,:]) / (phi_ana[4,:]1e-12) # z2m print(f地表最大相对误差: {err_surf.max():.2%}) print(f2m深最大相对误差: {err_shallow.max():.2%}) assert err_surf.max() 0.05, 地表误差超 5%模型有误 assert err_shallow.max() 0.10, 浅部误差超 10%需检查插值或网格 validate_uniform_halfspace()若地表误差 5%说明边界条件或源项处理错误若浅部误差 10%大概率是网格过粗或插值方法不当。这是mt_电法勘探编程_电法勘探正演_项目交付前不可绕过的“出厂测试”。5.2 三层模型响应的特征指纹识别正演正确性不仅看数值精度更要看是否复现经典地球物理特征。三层模型ρ₁/ρ₂/ρ₃的视电阻率曲线有明确指纹当 ρ₂ ρ₁ 且 ρ₂ ρ₃低阻层ρₐ 曲线呈“W”形极小值对应层厚当 ρ₂ ρ₁ 且 ρ₂ ρ₃高阻层ρₐ 曲线呈“驼峰”极大值位置与层顶深相关若正演曲线无此特征即使误差 1%也说明电导率赋值或层界面位置错误。实践中用matplotlib绘制不同层厚下的 ρₐ 曲线族观察形态演化比单纯看误差更有诊断价值。5.3 与野外实测数据的最小二乘拟合残差分析最终验证是拟合真实数据。取一段温纳装置测量剖面100 个点用正演引擎生成不同层参数h₁,h₂,ρ₁,ρ₂,ρ₃的响应计算 RMS 残差# 假设 obs_rho 为实测视电阻率数组100 点 def objective(params): h1, h2, rho1, rho2, rho3 params # 构建三层 sigma 矩阵 sigma_model build_layered_sigma(ny, nx, h1, h2, rho1, rho2, rho3) A, b build_fdm_matrix(nx, ny, dx, dy, sigma_model) phi spsolve(A, b).reshape((ny, nx)) rho_syn calculate_wenner_apparent_resistivity(phi, dx, a_spacing) return np.sqrt(np.mean((rho_syn[:100] - obs_rho)**2)) # RMS # 使用 scipy.optimize.minimize 调参 result minimize(objective, x0[5,10,100,10,100], methodL-BFGS-B) print(f最优层厚: h1{result.x[0]:.1f}m, h2{result.x[1]:.1f}m)若最优残差 RMS 5% 且拟合曲线形态匹配尤其拐点位置则证明正演引擎可投入实际解释。此时mt_电法勘探编程_电法勘探正演_已完成从代码到生产力的闭环——它不再是一个练习而是你手里的地质透视镜。本文还有配套的精品资源点击获取