基于Python的直流电阻率法三维有限元数值模拟

Three-Dimensional Finite Element Numerical Simulation of Direct Current Resistivity Method Based on Python

摘要: 利用Python平台强大的科学计算库，拟实现电导率分块均匀的直流电阻率法三维有限元数值模拟。基于Galerkin法给出稳定电流场三维双电源总场电位边值问题的控制方程，采用六面体单元剖分和三线性插值进行有限单元分析，将控制方程化为线性方程组，分别采用压缩稀疏列矩阵(CSC)与压缩稀疏行矩阵(CSR)存储总体系数矩阵，对比科学计算库SciPy及PyPardiso的求解速度与准确度，试算层状介质模型及单一、组合异常体模型，结果表明：H型层状介质的视电阻率测深解析解与有限元数值解最大相对误差低于3.00%，平均相对误差1.54%，本文算法有效；在Python平台求解地球物理正演中大型稀疏线性方程组时，PyPardiso与SciPy精度一致，当方程未知个数达37万时，PyPardiso比SciPy求解速度快约66倍；以CSC与CSR形式存储的大型稀疏线性系统，在同一求解器中求解效率相当。

Abstract: Utilizing the powerful scientific computing libraries on the Python platform, this study aims to implement a three-dimensional finite element numerical simulation for direct current (DC) resistivity with uniform conductivity blocks. By applying the Galerkin method, we derived control equations for the boundary value problem of a steady-state electric field with dual-source potential. Hexahedral meshing and trilinear interpolation were used for finite element analysis, converting the control equations into linear algebraic systems. The global coefficient matrix was stored in Compressed Sparse Column (CSC) and Compressed Sparse Row (CSR) formats to evaluate the speed and accuracy of solvers in SciPy and PyPardiso. Test cases included layered media models and single and combined anomaly models. Results show that, for H-type layered media, the maximum relative error between the apparent resistivity depth profiling analytical solution and the finite element numerical solution was below 3.00%, with an average relative error of 1.54%, validating the algorithm. When solving large-scale sparse linear systems in geophysical forward modeling on Python, PyPardiso and SciPy demonstrated consistent accuracy; for systems with up to 370,000 unknowns, PyPardiso achieved approximately 66 times faster performance than SciPy. For large sparse systems stored in either CSC or CSR formats, both achieved comparable efficiency with the same solver.