基于聚合型代数多重网格法的三维直流电法自适应有限元正演

摘要

在各向异性、起伏地形、真实地质模型电法模拟中，经自适应有限元离散后形成的大型稀疏线性系统存在内存消耗高、求解效率低等缺陷。为此，提出了聚合型代数多重网格（AGMG）法与自适应有限元法的联合算法，在提高正演精度的同时提升计算效率，实现复杂模型三维直流电法大规模正演模拟。对于三维直流电法满足的二阶椭圆边值问题，采用非结构化四面体网格的有限元法离散，并通过自适应策略进行局部加密，再利用AGMG法求解离散形成的大规模稀疏线性方程组。最后，通过复杂地电模型和实际地质模型验证了联合算法的有效性。在千万级自由度的求解中，联合算法比传统迭代法快了20多倍，比代数多重网格法快了近10倍，随着模型复杂度的提高，联合算法的效率优势更加明显。

关键词

代数多重网格法; 自适应有限元; 各向异性; 直流电阻率; 起伏地形

直流电法的数值模拟方法有很多^［1-5］，其中有限元法（FEM）理论基础完善，能够有效处理复杂地形，得到了广泛应用^［6-10］。自适应有限元法是提高正演精度的一种有效方法，在电位、电磁场变化较为剧烈的区域，如异常体边界、间断面、测点等位置，以迭代的方式实现网格自动调节，进而提高解的精度和分辨率。Zienkiewicz等^［11］提出了一种基于梯度恢复的后验误差估计算法，该算法简单、高效且不依赖特定问题。在此基础上，Ren等^［12］将该算法与四面体网格相结合并应用于直流电法正演中，极大地改善了直流电法正演的精度。随后，Ren等^［13］研究了面向目标的自适应有限元法，解决了带地形任意复杂直流电阻率各向异性高精度正演这一难题。殷长春等^［14］对任意各向异性复杂介质进行研究，解决了传统网格模拟非规则异常体精度不足的问题。尽管自适应有限元法提高了解的精度，但是在大规模问题中，直流电法边值问题仍存在计算量大、求解效率低等缺点。

多重网格（MG）作为一种快速求解器，在数学上已被证实是求解椭圆型偏微分方程的最有效方法之一，现已被广泛应用于地球物理领域^［15-19］。MG法可分为几何多重网格（GMG）法和代数多重网格（AMG）法。GMG法操作简单，效率也高于AMG法，但GMG法需要一系列嵌套的粗细网格，而且在各向异性及复杂模型问题中缺乏一定的通用性。为了解决粗网生成问题，基于MG法的无网格法也受到了广泛关注^［20-21］。此外，经典AMG法虽然不依赖求解问题的几何信息，但是存在以下缺陷：一方面，基于AMG法的直流电法大多采用结构化正交网格剖分，在非结构化网格中的应用很少，同时适用于各向异性、真实地形的求解器更是少之又少；另一方面，已有的正演方程组的计算规模普遍偏小。近年来，基于聚类思想的代数多重网格（AGMG）法^［22-23］受到了越来越多的关注。Chen等^［24］采用AGMG法快速求解由有限差分法离散三维直流电法边值问题形成的线性系统，极大地提高了求解效率。陈辉等^［25］利用AGMG法加速了基于交错有限差分格式的大地电磁正演。上述文献都是基于结构化正交六面体网格，没有充分发挥AGMG法不依赖几何网格的固有优势；对于结构化正交网格离散问题的求解，采用GMG法^［16-17］加速效果更好。

综上所述，针对复杂模型（含各向异性、起伏地表、复杂异常体等）正演效率低的问题，提出了AGMG法和自适应有限元法的联合算法。采用非结构四面体网格对计算区域进行剖分，利用基于非结构化网格的自适应有限元法对控制方程进行离散；结合边界条件生成大型稀疏线性方程组，引入AGMG法高效求解该线性系统。最后，通过理论模型与实际地质模型验证AGMG法的有效性、稳定性和收敛性等。

1 正演理论

1.1　椭圆方程边值问题及弱形式

直流电法勘探是基于地下岩石和矿石的电阻率差异进行的。直流电在地下传输时，会遇到不同电阻率地层，从而产生电位差。通过测量电极间电位差，可获得与地质体有关的直流电场分布，进而推断地下电性结构和岩层分布等信息。取$xoy$平面为水平地表，$z$轴垂直向下，三维直流电阻率正演问题可由以下二阶椭圆方程混合边值问题描述^［26］：

$$\left\{\begin{array}{l}\nabla \left(\mathit{\sigma }\nabla u\right)=-I\delta (\mathit{r}-{\mathit{r}}_{\mathrm{0}}),\mathrm{属于}\Omega \\ \mathit{\sigma }\frac{\partial u}{\partial \mathit{n}}=\mathrm{0},\mathrm{在}{\Gamma }_{\mathrm{1}}\mathrm{上}\\ \mathit{\sigma }\frac{\partial u}{\partial \mathit{n}}+\frac{\mathrm{c}\mathrm{o}\mathrm{s}(\mathit{r}-{\mathit{r}}_{\mathrm{0}},\mathit{n})}{\mathit{B}}u=\mathrm{0},\mathrm{在}{\Gamma }_{\mathrm{2}}\mathrm{上}\end{array}\right.$$

（1）

式中：$u$为待求的未知电位；${\mathit{r}}_{\mathrm{0}}=(x_{\mathrm{0}},y_{\mathrm{0}},z_{\mathrm{0}})$为电流源的位置，其电流大小为$I$；$\Omega$为有界求解区域且$\partial \Omega ={\Gamma }_{\mathrm{1}}\bigcup {\Gamma }_{\mathrm{2}},{\Gamma }_{\mathrm{1}}\bigcap {\Gamma }_{\mathrm{2}}=\mathrm{\varnothing }$，${\Gamma }_{\mathrm{1}}$、${\Gamma }_{\mathrm{2}}$上任一位置记为$\mathit{r}=(x,y,z)$；$\delta$为狄拉克函数；$\mathit{n}$为边界的法向量；$\mathit{\sigma }$为地下介质的电导率。若地下介质为均匀各向同性时，则$\mathit{B}={\sigma }^{-\mathrm{1}}\left|\mathit{r}-{\mathit{r}}_{\mathrm{0}}\right|$且$\sigma$为常数；若地下介质呈现各向异性时，$\mathit{\sigma }$为对称正定的电导率张量，且电导率张量与电阻率张量满足$\mathit{\sigma }={\mathit{\rho }}^{-\mathrm{1}}$，则$\mathit{B}=(\mathit{r}-{\mathit{r}}_{\mathrm{0}}{)}^{\mathrm{T}}\mathit{\rho }(\mathit{r}-{\mathit{r}}_{\mathrm{0}})$。对任一电导率张量，可由主轴电导率经过3次欧拉旋转计算得到，即：

$$\mathit{\sigma }=\left[\begin{array}{ccc}{\sigma }_{xx}& {\sigma }_{xy}& {\sigma }_{xz}\\ {\sigma }_{yx}& {\sigma }_{yy}& {\sigma }_{yz}\\ {\sigma }_{zx}& {\sigma }_{zy}& {\sigma }_{zz}\end{array}\right]=\mathit{R}\left[\begin{array}{ccc}{\tilde{\sigma }}_x& & \\ & {\tilde{\sigma }}_y& \\ & & {\tilde{\sigma }}_z\end{array}\right]{\mathit{R}}^{\mathrm{T}}$$

（2）

式中：${\tilde{\sigma }}_x\mathrm{、}{\tilde{\sigma }}_y\mathrm{、}{\tilde{\sigma }}_z$为主轴电导率；$\mathit{R}={\mathit{R}}_{\mathrm{S}}{\mathit{R}}_{\mathrm{D}}{\mathit{R}}_{\mathrm{L}}$，其中${\mathit{R}}_{\mathrm{S}}$、${\mathit{R}}_{\mathrm{D}}$、${\mathit{R}}_{\mathrm{L}}$分别为走向角${\alpha }_{\mathrm{S}}$、倾角${\alpha }_{\mathrm{D}}$、倾斜角${\alpha }_{\mathrm{L}}$的旋转矩阵，具体表达式可参考文献［26］。利用Galerkin有限元法对方程组（1）进行离散，得到大型稀疏线性系统。通过计算得到测量电极处的电位，由此可得不同排列装置的标量视电阻率或张量视电阻率^［27］。

1.2　自适应网格加密

参照文献［12］进行自适应加密，以某个内部单元为中心，将其相连的单元组成一个单元片，利用单元片计算得到数值解的恢复梯度$\mathit{G}\left(\nabla u\right)$。记$‖\mathit{e}‖$为所有子单元误差之和，定义全局误差$\beta =‖\mathit{e}\text{'}‖/‖\mathit{G}(\nabla u)‖$。设定自适应网格误差限为$\widehat{\beta }$。由恢复梯度与数值梯度之差等于梯度误差^［11］可知：${‖\mathit{G}\left(\nabla u\right)‖}^{\mathrm{2}}={‖\nabla u‖}^{\mathrm{2}}+{‖\mathit{e}\text{'}‖}^{\mathrm{2}}$。每个单元的允许误差不超过${‖\mathit{e}\text{'}‖}_m^{\mathrm{*}}=\widehat{\beta }{\left({‖\mathit{G}\left(\nabla u\right)‖}^{\mathrm{2}}/m\right)}^{\mathrm{1}/\mathrm{2}}$。若计算得到的单元误差比${‖\mathit{e}\text{'}‖}_m^{\mathrm{*}}$大，则进行局部加密。

2 AGMG求解器

AGMG法^［22-23］是一种特殊的AMG法，利用矩阵图论的思路，从局部搜索出发，将连接关系强的节点尽量聚集在一起，形成只与该节点相关的集合。AGMG法分为初始阶段和求解阶段2个过程。初始阶段包含代数磨光、粗细节点划分和插值算子构建三部分，求解阶段主要利用MG法的V-cycle形式嵌套迭代进行求解。

2.1　初始阶段

2.1.1　插值算子构建

考虑$n\times n$的大型稀疏线性方程组

$$\mathit{A}\mathit{u}=\mathit{b}$$

（3）

式中：$\mathit{A}$、$\mathit{b}$分别为联合算法离散后的系数矩阵和右端项；$\mathit{u}$为待求未知数，粗节点集合为$\mathit{u}$的子集。假设细节点的个数为$n$（即$\mathit{u}$的个数），粗节点的个数为$n_{\mathrm{c}}$，矩阵$\mathit{P}$为由细节点至粗节点的限制算子（即一个$n\times n_{\mathrm{c}}$的矩阵）。通过Galerkin形式得到粗节点下的矩阵

$${\mathit{A}}_{\mathrm{H}}={\mathit{P}}^{\mathrm{T}}\mathit{A}\mathit{P}$$

（4）

${\mathit{P}}^{\mathrm{T}}$为一个插值算子，限制算子的转置。

由于AGMG法粗化中形成的聚类$C_i$是互不相交的子集，因此矩阵$\mathit{P}$可以被定义为布尔矩阵

$${\mathit{P}}_{ij}=\left\{\begin{array}{l}\mathrm{1},i\in C_j\\ \mathrm{0},\mathrm{其他}\end{array}\right.$$

（5）

每行只有一个非零元素$(\mathrm{1}\le i\le n,\mathrm{1}\le j\le n_{\mathrm{c}})$。

2.1.2　粗细节点划分

定义1  　记集合$U$为全体节点集，集合$F$为细节点集，集合$C$为粗节点集。

定义2  　对某一特定节点$i$，若对任意给定的节点$j$满足

$$K_i=\left\{j\in U:j\ne i,a_{ij}\ne \mathrm{0}\right\}$$

（6）

则称$K_i$为节点$i$在集合$U$上的邻接集，其中$a_{ij}$为矩阵$\mathit{A}$中的元素。

定义3  　对某一特定节点$i$，若对任意给定的节点$j$满足

$$S_i=\left\{j\in K_i:a_{ij}<-\varsigma \underset{a_{ik}<\mathrm{0}}{\mathrm{m}\mathrm{a}\mathrm{x}}\left|a_{ik}\right|\right\}$$

（7）

则$S_i$为节点$i$的强耦合集，其中$\varsigma$为阈值，在［0，1］之间（通常$\varsigma$取0.25最优）。记${\lambda }_i^{}$为节点$i$的测度（即集合$S_i$中元素的个数）。

以图1为例，初始节点共有12个，由定义3可得每个节点的测度，图1a显示每个节点的测度大小。选择测度最小的节点$i$作为初始粗节点，在所有节点中选择$a_{ij}$最小的节点$j$作为待定点，若与节点$i$有强连接关系，则化为聚类（如图1b，细节点$j$可能不止一个），否则单独作为粗节点。修正与节点$j$强耦合的其他节点的测度，再选择测度最小的节点作为下一个粗节点。重复此过程，直至所有节点全部标记，图1f验证了聚类之间是互不相交的。图1中，黑色节点为粗节点，灰色节点为细节点。

图1  AGMG粗化

Fig.1  AGMG coarsening

2.2　AGMG法与经典AMG法的区别

AGMG法与经典AMG法在粗化方式和插值算子的构建上存在以下不同：

（1）在粗化方式上，AGMG粗化以聚类方式划分粗细节点，聚类$C_i$之间互不相交，如图1所示。经典AMG粗化较为复杂，每个细节点会与多个粗节点及细节点建立联系。

（2）在插值算子上，AGMG法的插值算子为布尔矩阵，即每行只有一个非零元素（见式（7））；经典AMG法的插值算子构造非常复杂，以加权Jacobi迭代得到的粗网格上的残差方程$\mathit{A}\mathit{e}\approx \mathrm{0}$为例（$\mathit{e}$为真解与近似解的差），对每个细节点$i$，其插值算子由三部分构成，分别是与节点$i$相关的粗节点集$C_i$、与节点$i$强连接的细节点集$F_{i,\mathrm{S}}$、与节点$i$弱连接的细节点集$F_{i,\mathrm{W}}$，即：

$$a_{ii}{e}_i\approx -\sum _{j\in C_i}{a}_{ij}{e}_j-\sum _{j\in F_{i,\mathrm{S}}}{a}_{ij}{e}_j-\sum _{j\in F_{i,\mathrm{W}}}{a}_{ij}{e}_j$$

（8）

经过一系列转化可得到$e_i={\omega }_{ij}{e}_j$，其中${\omega }_{ij}$为相应权重。以自由度为8 024 268的方程组为例，分别采用经典AMG法和AGMG法进行粗化。由图2可知：AGMG法的粗化次数比AMG法少，同等粗化次数下，AGMG法的粗化率低于AMG法，即AGMG法的自由度更少。经计算，AGMG法的初始阶段时间为26.3 s，总计算时间为57.14 s，而经典AMG法的初始阶段时间和总计算时间分别为221.81、571.59 s。由此可见，在大规模数值计算中，AGMG法的计算效果要优于经典AMG法。

图2  AGMG法与经典AMG法的粗化对比

Fig.2  Comparison of coarsening between AGMG and classical AMG methods

3 数值实验

在数值实验中，将自适应有限元法和AGMG法联合，完成三维直流电法正演模拟，并验证联合算法在复杂地形和各向异性模型下的有效性。对求解区域分别采用TetGen和ParaView软件进行四面体网格生成与可视化，利用自适应有限元离散生成大 型线性方程组，并采用AGMG、Hypre‒AMG^［28］、Hypre‒AMG‒CG^［28］和SSOR‒CG等方法进行求解。程序运行环境为带有mpich编译器的64位Linux系统，硬件为Intel Xeon Gold 6248R CPU（48核），内存为192 GB。

3.1　精度验证与效率分析

3.1.1　球体模型

模型一为球体模型^［12］，如图3所示。模型区域为（2 000，2 000，1 200） m，取向下方向为$z$轴正方向，均匀半空间的电阻率${\rho }_{\mathrm{0}}$为1 $\mathrm{\Omega }\cdot \mathrm{m}$；异常体为球形区域，半径$R=$2.25 m，距地面深度$H=\mathrm{2}R$，电阻率${\rho }_{\mathrm{r}}$为10 $\mathrm{\Omega }\cdot \mathrm{m}$。采用pole-pole测量装置，电源中心为$A$（-10，0，0） m，从电源中心通入大小为1 A的电流，以坐标原点为中心沿$x$轴长度为10 m的位置处每隔0.25 m布设测点进行测量。

图3  球体异常体模型

Fig.3  Spherical anomaly model

利用自适应加密技术对计算区域进行剖分，在球体及源电极附近网格越来越密，剖分后的网格图由ParaView软件绘制。图4a所示为由3 441个节点和18 155个单元组成的初始网格；图4b为第1次加密后的截面，共有39 248个节点和232 712个单元；图4c、d为第2次、第3次加密的示意图，分别有 65 765个节点、396 036个单元以及206 069个节点、1 265 607个单元。

图4  四面体网格自适应加密示意图

Fig.4  Schematic diagram of adaptive tetrahedral mesh refinement

对不同的自适应网格误差，分别采用4种不同的求解器计算最密网格下的大规模线性系统，相同停机标准（相对残差小于10^-8）下计算结果如表1所示。当自由度为13 568 104时，AGMG法求解时间约为SSOR‒CG法的1/23，约为Hypre‒AMG法的1/10，证实了AGMG法在求解此类问题时的高效性；另外，AGMG法的迭代次数远远小于SSOR‒CG法，且迭代次数几乎不随自由度的增加而增加。

表1  不同求解器的数值结果（球体模型）

Tab.1  Numerical results of different solvers（sphere model）

自由度	SSOR‒CG法			Hypre‒AMG‒CG法			Hypre‒AMG法			AGMG法
	迭代次数	时间/s	内存/GB	迭代次数	时间/s	内存/GB	迭代次数	时间/s	内存/GB	迭代次数	时间/s	内存/GB
206 069	205	6.55	0.5	7	7.42	0.5	19	7.13	0.5	12	0.98	0.1
1 945 747	441	154.30	5.2	8	113.97	4.3	24	166.64	4.9	12	11.67	1.3
8 024 268	719	1 152.23	14.8	8	554.69	17.4	25	571.59	17.6	13	57.14	5.4
11 686 326	774	1 951.66	31.2	8	851.59	24.7	23	849.75	24.4	13	88.08	7.2
13 568 104	808	2 396.10	36.8	8	991.79	27.9	25	1 018.51	28.5	13	102.16	9.2

图5a为求解时间随自由度变化的曲线。可知，AGMG法拟合线的斜率最小，说明AGMG法具有更好的可扩展性。图5b为AGMG法收敛曲线，残差随着迭代次数单调递减，且在不同规模下有近乎相同的下降趋势。

图5  不同自由度下的数值对比

Fig.5  Numerical comparisons between different degrees of freedom

图6为模型经过3次加密达到自适应误差限（此处设定为3%）时的视电阻率变化曲线。在用AGMG

图6  不同算法视电阻率计算结果

Fig.6  Apparent resistivity calculation results of different algorithms

法求解过程中，将每一层网格的视电阻率数值解与精确解进行对比，随着网格不断加密，数值解越来越逼近精确解。图6b中，4种求解器在最密层的视电阻率数值解与精确解有相同趋势，在异常体区域与精确解存在一定偏差。

3.1.2　各向异性层状模型

各向异性层状模型（2 000 m×2 000 m×2 000 m）如图7所示，第1层厚度为20 m，采用pole-pole测量装置，沿$x$轴正向布设20个间距为2 m的测点，电极位于坐标原点。

图7  各向异性层状模型

Fig.7  Anisotropic layered model

与球体模型类似，4种求解器计算结果如表2所示。可见，在8 193 852自由度下，Hypre‒AMG‒CG法的求解时间约为SSOR‒CG法的1/6，而AGMG法相比Hypre‒AMG法计算速度提高了近9倍，仅仅是SSOR‒CG法的1/54。随着自由度的增加，AGMG法的迭代次数保持不变。

表2  不同求解器的数值结果（各向异性层状模型）

Tab.2  Numerical results of different solvers (anisotropic layered model)

自由度	SSOR‒CG法		Hypre‒AMG‒CG法		Hypre‒AMG法		AGMG法
	迭代次数	时间/s	迭代次数	时间/s	迭代次数	时间/s	迭代次数	时间/s
217 400	248	16.54	7	8.34	20	8.28	15	1.34
897 180	389	136.36	8	55.87	22	54.98	15	7.12
1 693 175	480	655.72	8	125.07	23	126.10	15	14.40
8 193 852	814	6 197.55	9	1 001.81	26	1 007.83	15	113.79

在自适应加密过程中，将初始网格、第2次加密和第4次加密的视电阻率数值解与精确解进行对比。图8横坐标代表$x$轴方向电极测量的位置，纵坐标分别代表视电阻率和相对误差。随着网格不断加密，数值解误差越来越小，在第4次加密网格中最大相对误差仅为0.284%。

图8  不同网格视电阻率计算结果

Fig.8  Apparent resistivity calculation results of different meshes

3.2　起伏地形

3.2.1　带地形的复杂模型

考虑如图9所示带地形的复杂模型^［29］，取地表水平方向为$xoy$方向，垂直向下为$z$轴正向。ABCD是棱台的一个切面，AB水平长度为45 m，A点到地面的垂直高度为45 m。模型内部有一个正六棱柱异常体，该异常体埋于AB下方1 m的位置，且正六棱柱的棱边长$a$为20 m，高$h$为40 m，电阻率${\rho }_{\mathrm{1}}$为100 $\mathrm{\Omega }\cdot \mathrm{m}$；在正六棱柱右侧相隔20 m处，有一个正六面体异常体，边长$b$为40 m，电阻率${\rho }_{\mathrm{2}}$为1 $\mathrm{\Omega }\cdot \mathrm{m}$，正六面体的上表面与正六棱柱的下表面平行，且两异常体的上下表面均与水平方向平行；均匀介质的电阻率${\rho }_{\mathrm{0}}$为10 $\mathrm{\Omega }\cdot \mathrm{m}$。$O$点坐标为（0，0，-45） m，从$O$点注入1 A的电流，在$y$轴方向-150 m至150 m的位置上，每隔20 m放置一个测量电极。

图9  复杂地形结构示意图

Fig.9  Illustration of complex terrain structure

如表3所示，在8 864 454自由度下，AGMG法的求解时间是SSOR‒CG法的1/16，Hypre‒AMG‒CG法的求解时间约为SSOR‒CG法的1/2，证实了在起伏地表复杂模型中AGMG法也可以完成大规模线性方程组的快速求解。

表3  不同求解器的数值结果（起伏地形）

Tab.3  Numerical results of different solvers (rolling terrain)

自由度	SSOR‒CG法		Hypre‒AMG‒CG法		Hypre‒AMG法		AGMG法
	迭代次数	时间/s	迭代次数	时间/s	迭代次数	时间/s	迭代次数	时间/s
105 025	177	2.32	7	2.92	18	2.87	12	0.45
156 545	203	4.56	7	5.01	20	5.04	12	0.75
894 432	374	58.66	7	42.46	22	44.16	17	5.84
5 060 274	648	608.57	8	324.44	23	325.80	19	42.52
8 864 454	791	1 445.82	8	705.06	24	751.81	19	88.49

在自适应加密过程中，绘制沿$yoz$平面的地表视电阻率切面图（见图10）。整个过程共加密3次，由上至下网格不断加密，在50 m测点位置处有清晰的异常体等值线分布；在源电极附近存在高阻异常体，该区域的视电阻率值偏高。初始网格的误差为28.32%，最密网格的误差为4.97%，验证了AGMG法在复杂模型中保持高效的同时，结合自适应有限元法保持了高精度。

图10  起伏地形下视电阻率剖面

Fig.10  Apparent resistivity profile of rolling terrain

3.2.2　各向异性起伏地形模型

图11为一个4 400 m×2 400 m×2 200 m的山脉峡谷模型^［30］，山脉顶部与峡谷底部均距地面10 m，且地形区域为220 m×120 m×80 m。采用pole-pole测量装置，于A点处注入电流，沿$x$轴正向从M点到N点依次布设19个间距为10 m的测点。山脉峡谷及半空间视电阻率均为倾斜各向异性，视电阻率具体参数如图11所示。

图11  各向异性起伏地形模型

Fig.11  Anisotropy model within rolling terrain

不同自适应加密环境下，最密网格中的大规模线性方程组在不同迭代方法上的计算结果如表4所示。在1 229 565自由度下，AGMG法的计算效率比SSOR‒CG法提高了约37倍，迭代次数仅有15次；与Hypre‒AMG法和Hypre‒AMG‒CG法相比，求解速度也快了近10倍。此外，与球体模型类似，在不同自由度下AGMG法每次迭代的相对残差迅速减小。

表4  不同求解器的数值结果（各向异性起伏地形）

Tab.4  Numerical results of different solvers (anisotropic rolling terrain)

自由度	SSOR‒CG法		Hypre‒AMG‒CG法		Hypre‒AMG法		AGMG法
	迭代次数	时间/s	迭代次数	时间/s	迭代次数	时间/s	迭代次数	时间/s
206 802	252	14.13	7	6.41	27	7.92	15	1.12
284 386	285	29.13	7	10.00	27	11.67	15	1.67
670 087	376	116.38	7	31.08	28	36.68	15	4.65
1 229 565	465	350.15	8	70.64	29	85.29	15	9.22

3.3　真实地质模型

为了进一步验证联合算法处理复杂问题的能力，以安徽省泥河铁矿真实地质模型^［31］为基础进行数值模拟。真实地质模型的经度和纬度分别为117.250 0°、31.083 3°，经已有的地质、钻孔和岩石物性测量得到的4个地层和2个矿体的视电阻率值如表5所示。在2个矿体正上方的地形布置21个观测剖面，每个观测剖面上有21个接收器，观测剖面面积为3 000 m×3 000 m，线间距为150 m。3个电源分别距观测剖面中心5 000、7 000、9 000 m，源的长度为500 m，电流大小为50 A。

表5  地质单元的视电阻率

Tab.5  Apparent resistivity of geologic unit

岩性单元	视电阻率/（$\mathrm{\Omega }\cdot \mathrm{m}$）
上双庙组（USF）	552.0
下双庙组（LSF）	495.0
砖桥组（BBF）	624.0
闪长玢岩（DP）	800.0
磁铁矿（MG）	12.5
黄铁矿（PR）	8.0

真实地质模型计算结果如表6所示。可以看出，AGMG法的迭代次数并没有随着自由度的增加而增加；在40 896 890自由度下，AGMG法求解时间约为SSOR‒CG法的1/23，与Hypre‒AMG法相比，计算效率同样提高了近10倍。

表6  不同求解器的数值结果（真实地质模型）

Tab.6  Numerical results of different solvers (realistic terrain model)

自由度	SSOR‒CG法		Hypre‒AMG‒CG法		Hypre‒AMG法		AGMG法
	迭代次数	时间/s	迭代次数	时间/s	迭代次数	时间/s	迭代次数	时间/s
146 618	148	3.62	7	5.27	17	4.93	12	0.73
1 208 988	301	68.85	8	69.40	18	60.77	12	6.82
9 153 622	557	1 057.09	8	657.38	20	593.40	12	67.49
40 896 890	828	8 905.53	8	3 456.72	21	3 257.59	13	372.51

图12为控制单元误差2.5%情况下自适应加密后不同网格的视电阻率剖面。

图12  真实地形下视电阻率剖面

Fig.12  Apparent resistivity profile of real topography

可以看出，随着网格的不断加密，视电阻率等值线所围成的区域越来越清晰。同时，随着剖面深度增加，低阻部分区域逐渐清晰起来，在z=-450 m处可清晰地看到低阻的等值线，与真实地质模型的2个矿床相吻合。

4 结论

（1）当正演得到的线性方程组未知数达千万量级时，AGMG法比SSOR‒CG法快20多倍，比

Hypre‒AMG‒CG法快近10倍。

（2）不论是各向异性模型，还是起伏地形或真实地质模型，随着AGMG法迭代的进行，相对残差快速下降，展现出很强的鲁棒性和可扩展性。

（3）AGMG法继承了GMG法的优点，不依赖实际几何网格，有望在其他地球物理大规模正反演问题中得到广泛应用。

作者贡献声明

潘克家：研究选题，研究思路提供，论文撰写指导。

王鹏德：数值计算，论文整体撰写。

胡双贵：论文整体结构安排及修改。

王晋轩：算法对比实验指导。

邱乐稳：提供层状各向异性模型及起伏各向异性模型等实验数据。

汤井田：提供真实地质模型等实验数据。