基于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.

1. 引言

目前，直流电阻率法三维有限元数值模拟已经取得了显著的进展。阮百尧等[1] [2]给出了电导率分块均匀及电导率分块连续变化的三维地电断面电阻率测深的有限元数值模拟方法。在处理边值问题方面，黄俊革[3]简化异常电位算法的边界条件以节约计算时间，汤井田[4]提出有限元–无限元耦合数值模拟方法以解决传统有限元截断边界所引起的问题。在优化网格剖分方面，任政勇等[5]实现了对复杂模型的完全非结构化四面体全自动剖分，给出了三维直流电阻率模拟中四面体网格的质量评价标准和最优指标，杨志龙[6]基于梯度恢复的局部网格自适应加密技术实现了网格合理剖分，显著节约计算内存并减少正演时间，赵宁等[7]应用h型自适应加密算法结合高阶形函数提高有限元解的精度并节约计算资源。在提高求解效率方面，刘斌等[8]采用预条件共轭梯度法求解大型稀疏线性方程组以提高正演速度，张钱江等[9]提出了一种近似边界条件方法将与场源位置相关的边界系数矩阵从刚度矩阵中分离出来以实现线性方程组的快速回代求解，欧东新等[10]利用改进的对称超松弛预条件共轭梯度法进行线性方程组求解，朱姣[11]利用直接求解器MUMPS对大型稀疏矩阵进行快速分解及线性方程组求解。

Python作为一种强大的科学计算平台，其具备丰富的科学计算和数据处理库，如NumPy，SciPy，Matplotlib，Pandas和PyPardiso等，能极大简化数据操作、数值计算与可视化工作。NumPy提供了高效的数组操作，SciPy则扩展了这些能力至科学计算领域，而PyPardiso专门优化稀疏矩阵的求解过程，显著提高处理大型稀疏线性代数方程组的计算效率。Matplotlib则是一个在数据分析和科学计算中功能强大的可视化工具。值得一提的是，在地球物理学领域，Python平台同样拥有众多出色的数据处理库。如用于地球物理正反演的开源库SimPEG [12]，其支持电阻率法、激发极化法和电磁法等多种物探方法；专注于地震学数据处理和分析的Pyrocko与ObsPy [13]；以及支持地球物理建模与反演的pyGIMLi [14]等。

基于Python丰富成熟的科学生态系统，本文将尝试利用Python进行直流电阻率法三维有限元数值模拟。为保证数值模拟的精度与效率，将测试对比SciPy与PyPardiso求解器的性能及准确度，同时分析在两种求解器中采取CSC与CSR两种不同的方式存储大型稀疏系数矩阵对其性能与精准度的影响，以挑选出更适合求解地球物理有限元正演中大型稀疏线性系统的科学计算工具。将水平层状介质模型的滤波算法计算结果与本文算法数值模拟结果对比，平均相对误差1.54%，再对单一低阻异常体与组合异常体模型进行数值模拟，验证了本文方法的有效性。

2. 边值问题

三维地电断面双电源总电位的边值问题可归纳为[2]：

$\nabla \cdot \left(\sigma \nabla U\right)=-2I\delta \left(A\right)+2I\delta \left(B\right)$ , $\in \Omega$ (1)

$\frac{\partial U}{\partial n}=0$ ，$\in {\Gamma }_s$ (2)

$\frac{\partial U}{\partial n}+\frac{1}{r_B-r_A}\left[\frac{r_B}{r_A}\cos \left(r_A,n\right)-\frac{r_A}{r_B}\cos \left(r_B,n\right)\right]U=0$ , $\in {\Gamma }_{\infty }$ (3)

式中，Ω是地面边界${\Gamma }_s$ 及无穷边界${\Gamma }_{\infty }$ 所围的区域，U是总场电位，$\sigma$ 是介质的电导率，$\delta \left(A\right)$ ，$\delta \left(B\right)$ 分别表示以电源点A，B为中心的$\delta$ 函数。$r_A$ ，$r_B$ 分别是测点至电源点A，B的距离，$\cos \left(r_A,n\right)$ ，$\cos \left(r_B,n\right)$ 分别是矢径$r_A$ ，$r_B$ 与外法向$n$ 之间夹角的余弦。

3. Galerkin法建立控制方程

式(1)的余项可写为

$R=\nabla \cdot \left(\sigma \nabla U\right)+2I\delta \left(A\right)-2I\delta \left(B\right)$ ， (4)

这里，R是标量余项。由此，Galerkin的加权积分可表示成

${\displaystyle {\iiint }_{\Omega }{W}_{i}R\text{d}V}=0$ ，

将式(4)代入上式，同时借助Green第一定理

${\displaystyle {\iiint }_v\left(\varphi {\nabla }^2\psi +\nabla \varphi \cdot \nabla \psi \right)\text{d}V}={\displaystyle {∯}_S\varphi \frac{\partial \psi }{\partial n}\text{d}S}$ ，

可以得到

$\begin{array}{l}{\displaystyle {\iiint }_{\Omega }{W}_{i}R\text{d}V}\\ ={\displaystyle {\iiint }_{\Omega }{W}_i\left[\nabla \cdot \left(\sigma \nabla U\right)+2I\delta \left(A\right)-2I\delta \left(B\right)\right]\text{d}V}\\ =-{\displaystyle {\iiint }_{\Omega }\sigma \nabla U\cdot \nabla W_i\text{d}V}+{\displaystyle {∯}_{\Gamma }{W}_i}\sigma \frac{\partial U}{\partial n}\text{d}S+2{\displaystyle {\iiint }_{\text{Ω}}{W}_{i}I\delta \left(A\right)\text{d}V}-2{\displaystyle {\iiint }_{\text{Ω}}{W}_{i}I\delta \left(B\right)\text{d}V}\\ =0,\end{array}$

再将式(2)和(3)代入上式，得到控制方程

$\begin{array}{l}{\displaystyle {\iiint }_{\Omega }\sigma \nabla U\cdot \nabla W_i\text{d}V}+{\displaystyle {∯}_{{\Gamma }_{\infty }}{W}_i\sigma \frac{1}{r_B-r_A}\left[\frac{r_B}{r_A}\cos \left(r_A,n\right)-\frac{r_A}{r_B}\cos \left(r_B,n\right)\right]U\text{d}S}\\ =2{\displaystyle {\iiint }_{\Omega }{W}_{i}I\delta \left(A\right)\text{d}V}-2{\displaystyle {\iiint }_{\Omega }{W}_{i}I\delta \left(B\right)\text{d}V}.\end{array}$ (5)

4. 有限单元法

4.1. 求解区域离散化

将求解区剖分为一系列互不重叠的六面体单元(图1)，即可将式(5)中对区域和无穷边界积分分解为对各小单元的体积分与面积分，每个单元的顶点称为节点。

Figure 1. Schematic diagram of hexahedral mesh generation (a) and hexahedral element (b)

图1. 六面体网格剖分(a)与单元(b)示意图

在完成网格剖分的同时，于求解区形成若干节点。待定函数(如函数U)及其导数在这些节点上的数值，称为节点函数值。这样，便用求解区内有限个离散的待定节点函数值来近似表征待定函数在连续空间中的分布。这些待定的节点函数值$U_1$ 、$U_2$ 、……、$U_n$ ，组成一组独立变量。所谓求待定函数的数值解，就是确定这些节点函数值[15]。

4.2. 线性插值

通常地，利用各节点函数值在各单元内作线性内插来逐个单元计算参数U。采用六面体单元和三线性插值，图2(a)是母单元(边长为2的正方体)，图2(b)是子单元，两者的坐标关系为

$x=x_0+\frac{a}{2}\xi$ ，$y=y_0+\frac{b}{2}\eta$ ，$z=z_0+\frac{c}{2}\zeta$ ，

其中，$x_0$ 、$y_0$ 、$z_0$ 是子单元的中点，a、b、c是子单元的三个边长，两单元间的微分关系为

$\text{d}x=\frac{a}{2}\text{d}\xi$ ，$\text{d}y=\frac{b}{2}\text{d}\eta$ ，$\text{d}z=\frac{c}{2}\text{d}\zeta$ ，

$\text{d}x\text{d}y\text{d}z=\frac{abc}{8}\text{d}\xi \text{d}\eta \text{d}\zeta$ 。

三线性插值的形函数为

$N_i=\frac{1}{8}\left(1+{\xi }_i\xi \right)\left(1+{\eta }_i\eta \right)\left(1+{\zeta }_i\zeta \right)$ ， (6)

这里${\xi }_i$ 、${\eta }_i$ 、${\zeta }_i$ 是点$i\left(i=1,2,\cdots ,8\right)$ 的坐标。在单元中，U的插值函数是

$U={\displaystyle \sum _{i=1}^8{N}_i{U}_i}$ ，

式中$U_i$ 是六面体单元八个顶点处的待定U值。

Figure 2. Parent element (a) and subelement (b)

图2. 母单元(a)与子单元(b)

4.3. 单元分析

用形函数$N_i$ 作权函数$W_i$ ，将式(5)中的积分分解为各单元e上的积分。

4.3.1. 式(5)左端的体积分

${\displaystyle {\iiint }_{{}_e}\sigma \nabla U\cdot \nabla N_i\text{d}V}=\left(k_{ij}\right)U_e$ ， (7)

其中

$U_e={\left(U_1,U_2,\cdots ,U_8\right)}^{\text{T}}$ ，

$K_{1e}=\left(k_{ij}\right)$ ，

$k_{ij}={\displaystyle {\iiint }_{{}_e}\sigma \left(\frac{\partial N_i}{\partial \xi }\frac{\partial \xi }{\partial x}\cdot \frac{\partial N_j}{\partial \xi }\frac{\partial \xi }{\partial x}+\frac{\partial N_i}{\partial \eta }\frac{\partial \eta }{\partial y}\cdot \frac{\partial N_j}{\partial \eta }\frac{\partial \eta }{\partial y}+\frac{\partial N_i}{\partial \zeta }\frac{\partial \zeta }{\partial z}\cdot \frac{\partial N_j}{\partial \zeta }\frac{\partial \zeta }{\partial z}\right)\frac{abc}{8}\text{d}\xi \text{d}\eta \text{d}\zeta }$ ，

将式(6)对$\xi$ 、$\eta$ 、$\zeta$ 微商后，代入上式积分，可得$k_{ij}$ ，且$k_{ij}=k_{ji}$ 。$k_{ij}$ 的具体计算公式如下

$\left(\begin{array}{c}{k}_{11}\\ k_{21}\\ k_{31}\\ k_{41}\\ k_{51}\\ k_{61}\\ k_{71}\\ k_{81}\\ k_{22}\end{array}\right)=\left(\begin{array}{ccc}4& 4& 4\\ 2& -4& 2\\ -2& -2& 1\\ -4& 2& 2\\ 2& 2& -4\\ 1& -2& -2\\ -1& -1& -1\\ -2& 1& -2\\ 4& 4& 4\end{array}\right)\left(\begin{array}{c}\alpha \\ \beta \\ \gamma \end{array}\right)$ ，$\left(\begin{array}{c}{k}_{32}\\ k_{42}\\ k_{52}\\ k_{62}\\ k_{72}\\ k_{82}\\ k_{33}\\ k_{43}\\ k_{53}\end{array}\right)=\left(\begin{array}{ccc}-4& 2& 2\\ -2& -2& 1\\ 1& -2& -2\\ 2& 2& -4\\ -2& 1& -2\\ -1& -1& -1\\ 4& 4& 4\\ 2& -4& 2\\ -1& -1& -1\end{array}\right)\left(\begin{array}{c}\alpha \\ \beta \\ \gamma \end{array}\right)$ ，

$\left(\begin{array}{c}{k}_{63}\\ k_{73}\\ k_{83}\\ k_{44}\\ k_{54}\\ k_{64}\\ k_{74}\\ k_{84}\\ k_{55}\end{array}\right)=\left(\begin{array}{ccc}-2& 1& -2\\ 2& 2& -4\\ 1& -2& -2\\ 4& 4& 4\\ -2& 1& -2\\ -1& -1& -1\\ 1& -2& -2\\ 2& 2& -4\\ 4& 4& 4\end{array}\right)\left(\begin{array}{c}\alpha \\ \beta \\ \gamma \end{array}\right)$ ，$\left(\begin{array}{c}{k}_{65}\\ k_{75}\\ k_{85}\\ k_{66}\\ k_{76}\\ k_{86}\\ k_{77}\\ k_{87}\\ k_{88}\end{array}\right)=\left(\begin{array}{ccc}2& -4& 2\\ -2& -2& 1\\ -4& 2& 2\\ 4& 4& 4\\ -4& 2& 2\\ -2& -2& 1\\ 4& 4& 4\\ 2& -4& 2\\ 4& 4& 4\end{array}\right)\left(\begin{array}{c}\alpha \\ \beta \\ \gamma \end{array}\right)$ ，

这里

$\alpha =\frac{\sigma }{36}\frac{bc}{a}$ ，$\beta =\frac{\sigma }{36}\frac{ca}{b}$ ，$\gamma =\frac{\sigma }{36}\frac{ab}{c}$ 。

4.3.2. 式(5)左端的面积分

若单元e的一个面$S_{3478}$ 落在无穷远边界上，则边界积分

${\displaystyle {∯}_e{W}_i}{N}_i\sigma \frac{1}{r_B-r_A}\left[\frac{r_B}{r_A}\cos \left(r_A,n\right)-\frac{r_A}{r_B}\cos \left(r_B,n\right)\right]U\text{d}S=K_{2e}{U}_e$ ， (8)

其中

$K_{2e}=\frac{bc}{36}\cdot \frac{\sigma }{r_B-r_A}\left[\frac{r_B}{r_A}\cos \left(r_A,n\right)-\frac{r_A}{r_B}\cos \left(r_B,n\right)\right]\left(\begin{array}{cccccccc}0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 4& 2& 0& 0& 2& 1\\ 0& 0& 2& 4& 0& 0& 1& 2\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 2& 1& 0& 0& 4& 2\\ 0& 0& 1& 2& 0& 0& 2& 4\end{array}\right)$ 。

同理，可得出其他落在无穷边界上的面积分。

4.3.3. 式(5)右端的积分

根据$\delta$ 函数的积分性质，即

${\displaystyle {\iiint }_{\text{Ω}}\delta \left(A\right)\text{d}V}=\{\begin{array}{l}0,A\notin \Omega \\ \frac{1}{2},A\in \Gamma \end{array}$ ，

容易写出式(5)右端的体积分

$2{\displaystyle {\iiint }_{\Omega }{N}_{i}I\delta \left(A\right)\text{d}V}-2{\displaystyle {\iiint }_{\Omega }{N}_{i}I\delta \left(B\right)\text{d}V}=I_A-I_B$ ， (9)

仅与两电源的$I_A$ ，$I_B$ 有关。

4.4. 总体合成

在单元内，将式(7)~(9)相加，然后扩展成全体节点组成的矩阵或列阵，再把全部单元相加，得到

$\begin{array}{c}F\left(U\right)=\sum F_e\left(U\right)=\sum \left(K_{1e}+K_{2e}\right)U_e-I_A+I_B\\ =\sum K_e{U}_e-I_A+I_B=\sum {\overline{K}}_{e}U-P\\ =KU-P,\end{array}$ (10)

其中$K_e=K_{1e}+K_{2e}$ ，${\overline{K}}_e$ 是$K_e$ 的扩展阵，$K=\sum {\overline{K}}_e$ ，$P={\left(0,\cdots ,I_A,\cdots ,I_B,\cdots ,0\right)}^{\text{T}}$ 。

令式(10)为零，得到线性代数方程组

$KU=P$ ， (11)

解方程组，得到各节点的电位函数值U。

5. 求解器与存储方式

基于结构化网格剖分的有限单元法中，由于单元的连通性是固定且有规律的，每个节点只与邻近的节点有联系，即只有邻近的节点会出现在相同的单元内，因此，在系数矩阵中对应的非零元素大多集中在矩阵的对角线上或对角线附近，形成了条带状的非零元素区域。

Figure 3. Distribution map of non-zero elements in a structured grid’s coefficient matrix

图3. 结构化网格系数矩阵中非零元素分布图

以六面体网格剖分为例，如图1(a)，网格剖分的节点数量为11 × 11 × 6 = 726，其全体节点组成的系数矩阵(726 × 726)非零元素分布如图3，其非零元素占比2.42%。随着网格节点数量的增加，系数矩阵的稀疏程度将越来越大，可通过定带宽[16]或者变带宽[17]-[20]的存储方式以减少存储量。而定带宽或变带宽的存储方式依然存储了很多零元素，采用稀疏存储方式可进一步压缩存储量[21]。

基于Python平台，分别采用压缩稀疏列矩阵(CSC)与压缩稀疏行矩阵(CSR)存储总体系数矩阵，对比分析SciPy及PyPardiso两种科学计算库的求解器速度。其中，SciPy采用专门处理稀疏线性系统的scipy.sparse.linalg.spsolve [22]方法进行求解；PyPardiso则采用pypardiso.spsolve [23]方法。本文的有限元数值模拟程序均在同一计算机上运行及测试，操作系统Windows10 64位，CPU型号Intel (R) Core (TM) i5-8500，内存48 GB (DDR4 2666 MHz)。如表1，在求解地球物理正演中的大型稀疏线性方程组时，以CSC或CSR形式存储的大型稀疏线性系统，在同一求解器中两者求解速度相差微乎其微；随着网格单元节点总数的增加，SciPy与PyPardiso两种求解器用时均有所增长；当节点总数超过5000后，PyPardiso求解效率明显高于SciPy；当节点总数达37万时，PyPardiso求解时间约17 s，而SciPy耗时急剧增长，约1130 s，两者相差约66倍，这是由于SciPy求解线性系统的未知向量非稀疏时会提高计算成本[22]，而直流电阻率三维有限元正演中待求解的全体节点电位为稠密向量，与基于SuperLU分解默认单线程运行的SciPy而言，PyPardiso是专门针对稀疏矩阵优化的求解器，且支持多进程并行计算，具有高性能、稳健、节约内存等优点。

Table 1. Solving speed of SciPy and PyPardiso with different storage methods

表1. SciPy与PyPardiso中不同存储方式的求解速度

网格剖分节点总数	求解用时(s)
	SciPy		PyPardiso
	CSC	CSR	CSC	CSR
021 × 21 × 11 = 004851	0.14	0.14	0.14	0.14
073 × 23 × 17 = 028543	5.13	5.11	0.38	0.38
135 × 35 × 23 = 108675	61.55	61.96	1.83	1.87
141 × 45 × 28 = 177660	191.64	191.74	3.37	3.32
141 × 59 × 35 = 291165	816.95	812.41	11.74	11.82
181 × 59 × 35 = 373765	1129.30	1132.46	16.95	16.97

6. 数值算例

6.1. 模型1：水平层状介质模型

首先利用水平层状介质模型的对称四极测深数据检验本文算法的有效性。模型1为H型水平层状介质(图4)，第1层厚度为2 m，电阻率为100 Ω∙m，第2层厚度为2 m，电阻率为10 Ω∙m，第3层电阻率为200 Ω∙m。

Figure 4. Model 1: H-layered media

图4. 模型1：H型层状介质

采用140点一阶Hankel滤波系数的滤波算法[24]得到的正演结果，作为模型1的解析解，同时对比分析本文算法中SciPy与PyPardiso在CSC与CSR存储方式下的精度。模型1的视电阻率测深曲线如图5(a)，可见SciPy与PyPardiso求解器在不同存储方式下的正演结果与解析解拟合较好。两种求解器在不同存储方式下数值模拟结果的相对误差曲线如图5(b)，容易看出它们精度一致，由于总场算法电源点附近电位的奇异性，距场源最近的浅部测深点相对误差为3.26%，其他位置最高相对误差均低于3.00%，平均相对误差为1.54%，表明本文算法行之有效。

Figure 5. Model 1 sounding curve of apparent resistivity (a) and relative error curves for two solvers with different storage methods (b)

图5. 模型1视电阻率测深曲线(a)及两种求解器在不同存储方式下的相对误差曲线(b)

Figure 6. Schematic diagram of Model 2 (a) and contour section of apparent resistivity for dipole device (b)

图6. 模型2示意图(a)及偶极装置视电阻率等值线断面图(b)

6.2. 模型2：单一异常体模型

模型2为均匀大地中单一低阻块状体，如图6(a)，其沿x，y，z方向的规模为20 m × 30 m × 20 m，电阻率为10 Ω∙m，顶面埋深20 m，围岩电阻率为100 Ω∙m，以低阻体顶面中心在地面上的投影作为坐标原点。采用偶极装置，偶极长度AB = MN = 10 m，分别在y = 30 m，15 m，0 m，−15 m，−30 m测线上进行偶极测深的数值模拟结果如图6(b)。在主测线上(y = 0 m)，低阻体异常最明显，从低阻体埋深处开始，视电阻率呈现极小值，随着深度的增加，视电阻率极小值圈向两侧扩展形成“八”字型双谷。当测线远离主测线，视电阻率极小值逐渐趋于背景场，低阻体产生的异常随之减弱。各断面上的视电阻率等值线图能较好表示出低阻体产生的视电阻率异常。

Figure 7. Schematic diagram of Model 3 (a) and contour section of apparent resistivity for dipole device (b)

图7. 模型3示意图(a)及偶极装置视电阻率等值线断面图(b)

6.3. 模型3：组合异常体模型

模型3为均匀大地中平行排列的高阻与低阻组合异常体，如图7(a)，两者规模、走向及顶面埋深与模型2中单一低阻块状体一致，横向间距20 m。其中，低阻块状体电阻率为10 Ω∙m，高阻块状体电阻率为1000 Ω∙m，围岩电阻率为100 Ω∙m。采用与模型2一致的排列方式进行数值模拟，得到多条测线上视电阻率等值线断面图，见图7(b)。在主测线上(y = 0 m)，高阻体与低阻体产生的异常相互叠加，视电阻率极小值圈顶、极大值圈顶大致对应低阻体、高阻体埋深位置。当测线远离主测线，高、低阻体产生的异常随之减弱，视电阻率极小值逐渐增大趋于背景场，视电阻率极大值逐渐减小趋于背景场。各断面上的视电阻率等值线断面图能较好表示出组合异常体产生的视电阻率异常。

7. 结论

1) 本文基于Python平台，实现了电导率分块均匀的直流电阻率法三维有限元数值模拟，通过对比水平层状介质模型对称四极视电阻率测深曲线的解析解与有限元数值解，除近场源处相对误差达3.26%，其他位置相对误差均低于3.00%，平均相对误差1.54%，表明本文算法准确、有效。

2) 通过对比分析SciPy与PyPardiso求解器性能、精度的试验发现：两者计算精度相当，而PyPardiso效率更高，当方程组未知个数达37万时，PyPardiso比SciPy求解速度快约66倍；以CSC或CSR形式存储的大型稀疏线性系统，在同一求解器中计算效率相当。因此，基于Python平台求解直流电阻率法三维有限元数值模拟中的大型稀疏线性方程组时，优选支持多进程并行计算的PyPardiso高性能求解器可有效保证计算效率。

基金项目

广西自然科学基金(2024GXNSFBA010257)；广西人才基地专项基金(桂科AD23026231)；江西省地质局青年科学技术带头人培养计划项目(2024JXDZKJRC10)资助。