1 引言

直流电阻率法早已广泛应用于铁路、公路、水利、城建、地质和矿产勘察等生产实践当中，而现代电子技术的发展使得高效率的三维数据采集成为可能.海量三维数据的反演依赖于强大的正演计算，因此，如何在所需的计算精度范围内加速正演计算成为了一个极有价值的研究课题.对于目前广泛使用的有限单元法，传统的混合边界条件可在合理的计算范围内得到较高的正演精度，但是由于刚度矩阵与电源位置有关，每移动一次电源位置，就需重新生成刚度矩阵和解方程组.对于需要多次移动电源的探测方法，例如偶极-偶极方法，一条断面的测量就涉及数十次的电源移动，使三维地电断面的正演模拟十分耗时，以此为基础的反演更是没有实用价值.目前广泛采取的解决方案是用齐次边界条件（Neumann边界条件）^[1]或Dirichlet边界条件来取代混合边界条件^{[2, 3]}，即在半无限边界上强制∂u/∂n=0或u=0，其中u为电位，n为边界外法相单位矢量；或者认为边界积分函数中cos（r，n）/r为常数，提至积分号以外而不考虑其变化^{[4, 5]}.采用上述这些方案必须将计算范围取得更加大，才能减少边界对计算结果的影响，这必定增加网格节点数，从而增加计算量.

为解决上述问题，本课题组在传统有限元数值模拟的基础上引入了无限单元法(Infinite Element Method)，使得刚度矩阵与电源位置不再相关，并且保持计算区域相对较小.目前，无限元大体可分为Bettess映射无限元^{[6, 7]}，Astley波包映射无限元^{[8, 9]}(Mapped wave envelope elements)和Burnett多极展开(Multipole expansion)无限元^{[10, 11]}.其中Astley波包映射无限元应用最为广泛.其基本思想是坐标变换和位场采用不同的插值函数，通过坐标映射实现无限域的积分，并在映射后的局部坐标里建立位场插值形函数，使位场在有限元网格和无限元网格结合处保持一致性，而在无限远处衰减为零^[12].无限元目前已较多的应用于声学、电磁学、岩土力学等领域^[13~17].在地球物理领域，Fu和wu^[18]在弹性波方面应用了无限元来处理吸收边界条件；Blome等^[19]应用Astley无限元处理直流电法边界问题，缩小了计算范围，减少了节点数，并且避免了采用Dirichlet或混合边界条件，但是其采用的无限元形函数较为复杂，且由于引入了Astley无限元中特有的额外权重因子，使整体系数矩阵失去了对称性，本文的实验还表明该种无限元会在某些情况下导致线性方程组不收敛.国内，张玉娥等^[20]用于研究地铁区间隧道地震反应分析；姜忻良^[21]等用于研究交叉隧道地震反映；朱军等在三维电测井中成功应用无限元来改进计算速度.在直流电法领域，孙跃^[22]提出了用有限元、无限元相结合来计算三维稳定电流场的方法，效果优于Dirichlet边界条件，但仅给出了基于六面体单元的基本公式和简单结果，对无限元映射函数和形函数均未作详细说明.

本文由三维点源电场边值问题出发，利用Galerkin方法简单推导了有限元基本公式.随后推导了三维情况下无限元坐标映射函数.之后，提出了一种全新的无限元形函数并应用于本文提出的有限元-无限元耦合法中.它避免了传统Astley无限元法中引入额外权重因子导致系数矩阵失去对称性的缺点，并且与多种常见无限元形函数的对比显示，基于此种无限元形函数的算法在精度和时间消耗上均占优势.最后，对称四极测深和偶极-偶极装置模型验证了本算法的有效性.一系列的数值模拟实验表明，本文提出的有限元-无限元耦合法具备以下优点; （1）可以在近似测区大小的计算范围内得到与混合边界条件相当的计算精度，并优于同等计算范围下齐次边界条件的解，有利于减少节点数，从而减小计算量；（2）电位在整个半空间区域内分布，无需再考虑截断边界条件，刚度矩阵与电源位置无关，移动电源后不需改变刚度矩阵，有利于加速反演计算；（3）由于采用非结构化四面体有限元剖分，在边界外耦合的无限单元数量不多，并且不引入额外的权重因子，系数矩阵可保持稀疏对称性，有利于用迭代法快速求解.

2 有限单元法基本方程 2.1 三维点源场边值问题^[23]

点源置于三维无限半空间Ω的表面Γ_s时，电位满足的微分方程

(1)

其中σ为介质电导率，u为电位，I为点源供电电流，δ（A）为关于场源所在位置A点的δ函数，其满足

(2)

在地面Γ_s上，电位满足的边界条件为

(3)

在半无限边界上Γ_∞上，电位满足混合边界条件

(4)

其中n为边界外法向单位矢量，r为电源点距边界的距离.区域Ω内部的电阻率边界为自然边界条件，可以不予考虑.

2.2 有限单元法

采用四面体单元和线性插值，得到单元中的插值函数;

(5)

其中N_i为插值形函数，与四面体自然坐标L_i相同，u_i为节点上的电位值.

基于变分原理，可得到与上节中边值问题等价的变分问题^[23];

(6)

基于Galerkin加权余量法，也可得到相同的方程形式.

令

(7)

其中ω_i为权函数，在Galerkin法中一般取权函数为形函数N_i（i=1，…，4）.根据Green定理，有

(8)

其中为边界Γ的外法向单位矢量.将（3），（4）和（5）式代入（8）式，得

(9)

令R=0，则

(10)

对四面体单元的四个节点进行组合，可得到单元积分的矩阵型式.其中

(11)

这里

由全局坐标和局部坐标的转换关系;

(12)

(13)

其中[J]为雅克比矩阵

(14)

这样，全局坐标下的单元积分就可转化到局部坐标下的积分.将全部单元相加，即可得到整体矩阵K₁.（10）式中的其他积分均可按上述思路转换，混合边界条件产生的积分;

(15)

点源积分

(16)

将上述积分相加，得线性代数方程组

(17)

解方程组，即可得到各节点上的点位u.

3 无限元与有限元结合 3.1 无限元几何映射

本文主要讨论Astley无限单元^{[8, 9, 13, 14]}.在有限元网格的边界Γ_∞外增加无限单元，通过坐标映射将单元在某一方向扩展到无限远处，这样就可实现无限远处的积分，也就无需再考虑边界条件.与有限元法不同，其坐标变换和位场插值采用不同的插值函数，分别称为映射函数（Mappingfunction）和形函数（Shapefunction）.一维情况下的映射原理如图 1所示.

整体坐标下的x∈ [a，∞）被映射到了局部坐标下的v∈[-1,1]有限范围内，映射关系为

(18)

其中M₁（v）和M₂（v）为映射函数;

(19)

这种映射具备惟一性，并且在二维和三维情况下均适用.由（18）式和（19）式可以得到;

(20)

上式表明v呈1/x型衰减.对于电位形函数，设v=-1，0，1处的电位值分别为u₁，u₂，u₃，假使采用传统有限元二次插值形函数，即

(21)

将（20）式代入（21）式，并考虑到在无限元处u₃=0，得

(22)

上式表明u呈2型衰减，符合点电源在无限半空间中衰减的物理规律，因此将此种映射应用于点源直流电法问题具备合理性.

将上述映射推广到三维情况，采用非结构化四面体有限单元与六节点无限单元^{[12, 19]}相结合，其单元映射示意图如下;

映射原点P为所有无限单元的映射出发点，一般可取在上表面的中心点处.以P点为原点，经过有限单元在外边界上的三角形面的三个顶点，向外引出三条射线，形成无限单元的几何形式.点P到点4、5、6的距离分别为点P到点1、2、3的距离的两倍.在无限远处还有无限单元的三个节点，它们使得势场分布扩展到无限远处，但由于无限远处电位为0，故此三点不参与插值运算，因此也就没有必要考虑.在局部坐标中，无限单元被映射为一个有限大小的三棱柱，从而可以进行积分运算.其中ζ方向表示无限远方向，其映射方式与前文所述的一维无限单元映射方式相同，电位插值则采用特殊的形函数；在ξ-η平面内，无限单元与有限单元采用相同的映射形式和形函数，即三角形面单元映射和面积坐标线性插值^[23].

坐标映射为

(23)

其中L_i为三角形面积坐标.结合图 2，有

(24)

再结合ζ方向的坐标映射关系式（18），得到六节点无限单元的结点坐标映射函数;

(25)

其中

(26)

经过上述坐标映射之后，电位借助于无限单元延伸至无限远处，并且在无限远处衰减为零.此时便无需再考虑边界条件，即计算中无需考虑边界积分（15）式.无限单元的单元分析与1.2节中的有限元单元分析相类似，只是四面体单元变为了六节点无限元单元.

3.2 无限元形函数

在ξ-η平面内，电位形函数与坐标映射函数相同，均为三角形面积坐标，见式（24）.在图 2中的ζ方向，无限单元采用特殊的电位插值形函数.该方向上形函数的选择一直是研究和争论的热点，只有根据不同的场源性质选择合适的形函数才能达到较好的效果.李录贤等^{[13, 14]}对多种常见的无限元的形函数做了简介.国内岩土力学领域^{[12, 24, 25]}常通过2.1节所述的坐标映射后直接采用对应的有限元形函数，如式（21）所示；或者采用Lagrange插值函数与某一衰减函数的乘积构造“衰减型”无限元，但其坐标映射方式与前文所述不同.Dreyer和Estorff^[26]在三维开域声学（Exterioracoustics）有限元-无限元仿真中基于Helmholtz方程的求解比较了多种形函数形成的系数矩阵的条件数，这些形函数分别基于Lagrange多项式、Legendre多项式和Jacobi多项式，最终认定基于Jacobi多项式P_i^2，0的形函数产生的系数矩阵性质优良，适合迭代法求解.Blome等^[19]在三维直流电法数值模拟中采用了上述形函数，并取得了较好的结果.值得指出的是上述方法采用加权余量法推导无限元积分表达式，并且在引入权函数时加入了额外的权重因子，最终形成了系数矩阵为非对称矩阵，而且实验表明在某些模型中还会出现方程组不收敛的情况.

由文献[13]，声辐射问题的声压p由Helmholtz方程

(27)

控制，并在远场满足Sommerfeld辐射条件;

(28)

其中，k为波数，依据Astley-Leis映射无限元相关文献中的表述^[27]，无限远方向的无限元形函数应为如下形式

(29)

其中v=1，…，n，表示该方向上的n个插值点r₁，…r_n，r为整体坐标中计算点离无限元发散原点的距离.根据坐标映射关系式（18）

(30)

所以有

(31)

F_v（）为一个n-1阶的多项式.它通常可取以下几种形式;

(32)

(33)

其中L_v（a/r）为基于插值点r₁，…r_n的n-1阶Lagrange插值多项式

(34)

注意到整体坐标下最后一个插值点在无限远处，即r_n→ ∞，上式与局部坐标ζ表示下的Lagrange插值多项式

是完全等价的.

(35)

其中P_v-1^s（a/r）为v-1阶Legendre多项式或Jacobi多项式.在第一个点，即v=1，ζ=-1时，F₁=1；在其他插值点上则需要通过+1或-1来使得F_v（a/r）=0，即在局部坐标ζ=-1时其他插值点上插值权重均为0，从而使得无限单元与有限单元在结合点处保持一致性.

通常Astley型无限元的单元积分都是通过Petrov-Galerkin加权余量法推导的.为保证单元积分可积，其权函数取额外的权重因子（a/r）²乘以形函数的共轭;

(36)

将（29）式和（36）式转换到局部坐标下，有

(37)

(38)

对于直流电场问题，电位由（1）式控制，不含波数k，因此指数部分可以不予考虑.

本文基于半空间中的三维球体模型，比较了采用多种无限元形函数情况下的计算精度和计算时间，最终选取一种最优的无限元形函数，形成有限元-无限元耦合法.这些形函数分别为

（1）国内岩土力学领域常用的传统有限元二次Lagrange插值函数

(39)

C₃作用于无限单元的无限远节点处，此处电位已衰减为零，因此插值中可以不用考虑.结合ξ-η平面内形函数（24）式，可得无限单元电位插值形函数

(40)

采用此形函数直接应用于无限元单元的单元分析，不需考虑指数项和权重因子等，最终形成的系数矩阵为对称矩阵.

（2）改进型二次Lagrange插值函数

在无限元形函数的构建中，F_v（a/r）前的系数a/r对整个形函数起到修正的作用，保证r→ ∞时形函数趋于0，由式（31）已知在局部坐标中该系数为（1-ζ）/2.直流电场问题中没有波数k，不需考虑表示相位的指数部分.额外的权系数（1-ζ）²/2是采用加权余量法时引入的，在使用变分原理或权系数与形函数相同的经典Galerkin加权余量法推导时不存在这个系数.因此我们尝试在传统有限元二次Lagrange插值函数前乘以修正系数（1-ζ）/2，再组合单元矩阵，这样得到的系数矩阵依然为对称矩阵.此时无限单元电位插值形函数为

(41)

（3）改进型四次Lagrange插值函数

与改进型二次Lagrange插值函数构建方式类似，只是将二次插值换为四次插值，主要目的是比较插值次数的增加对精度有无改善.在局部坐标中插值点分别取，构造方式为

(42)

此时无限单元电位插值形函数为

(43)

（4）考虑额外加权因子的改进型二次Lagrange插值函数

按式（37）和式（38）分别构造插值形函数和权函数，考虑额外的权重因子（1-ζ）²/2但不考虑指数项.电位插值形函数与式（41）相同.而单元积分采用加权余量法的形式，此时系数矩阵元素应为

(44)

其中M_i为无限元形函数；D为权重因子（1-ζ）²/2.由于引入了额外权重因子D，形成的系数矩阵为非对称矩阵.

（5）考虑额外加权因子的改进型四次Lagrange插值函数

电位插值形函数与式（43）相同，而单元积分采用加权余量法的形式，权函数为额外的权重因子（1-ζ）²/2与形函数的乘积.此时单元积分的系数矩阵元素形式与式（44）相同.由于引入了额外权重因子D，形成的系数矩阵为非对称矩阵.

（6）考虑额外加权因子的五次Jacobi多项式

此种形函数为Blome所采用^[19]，其构造方式遵从式（35）.Jacobi多项式选用P_n^2，0，其具体表达式可由下式获得

(45)

此时无限单元电位插值形函数为

(46)

单元积分采用加权余量法推导的形式，考虑额外权重因子，系数矩阵元素表达式同式（44）.由于引入了额外权重因子D，形成的系数矩阵为非对称矩阵.

4 算例及分析 4.1 多种形函数计算效果比较

基于半空间中的三维球体模型，我们对上节中所述六种形函数所构建的有限元-无限元耦合方程分别进行了计算，并对比了结果精度.模型中球体半径为2.25m，中心埋深5m，置于上表面中心点正下方.球体电阻率为1Ωm，无限半空间电阻率为10Ωm.采用AM装置，点源A置于坐标（0，-5，0）处，40个测点沿y轴在球体正上方的地面上排列，点距0.25m，由（0，-4.75，0）至（0，5，0）.整体模型如图 3所示.计算程序基于本文自行编制的三维有限元软件^[28]，采用非结构化四面体单元进行网格剖分，再在外边界处耦合无限元单元.通过不断缩小无限半空间的范围，即有限元网格剖分范围，与解析解对比计算结果.

计算模型中无限半空间的范围（长×宽×高）分别为50m×50m×25m，20m×20m×15m，15m×15m×15m和12m×12m×10m.基于第（1）至第（3）种无限元形函数的计算方案将得到对称的全局矩阵，因此采用不完全LU预处理（iLU precondi-tioned）的CG（ConjugateGradient）迭代法求解；基于第（4）至第（6）种无限元形函数的计算方案将得到非对称的全局矩阵，因此采用不完全LU预处理的QMR（Quadratic Minimum Residual）迭代法求解.与解析解对比得到的相对误差分析如表 1所示.

对于范围大小不同的每个计算模型而言，虽然采用的无限元形函数不同，但模型剖分情况是相同的，计算条件也是相同的，只是对于对称矩阵和非对称矩阵分别采用了不同的迭代求解方法，因此平均相对误差限、最大相对误差限等参数具有相对可比性.由表 1中数据可以清楚地看出，基于第（2）种无限元形函数（即改进型二次Lagrange插值函数）的计算方案在各个模型中均表现良好，无论是计算精度还是计算时间，大都领先于其他形函数的方案；增加Lagrange插值阶数或考虑额外的权重因子D并没有明显提高计算精度，反而因为增加了自由度而增加了计算量和计算时间；基于Jacobi多项式的形函数计算精度并不高，且在15 m×15 m×15 m和12m×12m×10m两个模型中出现了不收敛的情况（迭代次数达到了设定的上限5000）.虽然Dreyer等的研究表明基于Jacobi多项式的形函数形成的系数矩阵具有更好的条件数，有利于快速收敛，但此研究是基于开域声学中的Helmholtz方程讨论的，本实验表明该结论并不适用于三维直流电法问题.整体而言，这些算例表明无限元方法适用于处理三维直流电法边界问题，并且可以在仅包含测区的计算范围内达到符合要求的计算精度，有利于减少节点数，加快计算速度；其中第（2）种无限元形函数，即修正项（1-ξ）/2与二阶Lagrange插值多项式相乘而得到的形函数，在计算精度和速度上都最具优势，因此在我们提出的有限元-无限元耦合法中均采用此种无限元形函数.

4.2 其他算例

为了验证本文提出的有限元-无限元耦合计算方法的可靠性和实用性，我们基于第（2）种无限元形函数对其他几种三维模型进行了计算.模型描述及结果如下.

4.2.1 三层水平大地模型

K型三层水平大地，第一层电阻率为50 Ωm，厚度为5m；第二层电阻率为100Ωm，厚度为5m；第三层电阻率为10Ωm.模型大小为160m×160m×50m，有限元网格采用非结构化四面体网格剖分.采用对称四极测深装置，测点M、N置于地表中心处，极距MN=1m，逐点移动A、B进行测深.计算结果与数值滤波解析解相对比，计算结果及相对误差分布如表 2所示.

由表 2可以看出，在电源点A、B距边界较远时，即AB/2较小时，传统混合边界条件有限元法、有限元-无限元耦合法和齐次边界条件（Neumann）有限元法得到的结果大致相同，精度也都较高.当电源点A、B距边界越来越近时，即AB/2逐渐增大时，混合边界条件和齐次边界条件解的误差均呈明显的上升趋势.其中，齐次边界条件下，即不考虑边界积分项时，当AB/2=33m时误差就超过了5%，随后急剧增大.但是有限元-无限元耦合法在距边界最近的点处（AB/2=65 m时）仍能得到较好的解，甚至优于混合边界条件.总体而言，本文提出的有限元-无限元耦合法适用于层状介质上的对称四极测深，最远测点距边界仅为15 m时仍能得到合理的解，且精度优于相同计算范围时齐次边界条件下的解.

4.2.2 垂直两地层含一六面体模型

本模型较为复杂，为垂直两层大地模型并在一侧内埋藏一立方体异常矿体.Hvozdara^[29]曾用边界积分方程法对其进行过求解.模型如图 4所示，左边地层电阻率为10Ωm，右边地层电阻率为100Ωm，电阻率为10 Ωm的立方体位于右边地层中.采用Cartesian坐标系，两地层的分界面与地表面的交线与x轴重合，原点位于交线中点.立方体边长为2m，中心点坐标位置为（0，3，3）.采用对称四极测深装置，M、N坐标分别为（0，-3.4，0）和（0，-2.6，0），点源A、B沿y轴分布.

边界积分方程法、有限元-无限元耦合法、混合边界条件有限元法与齐次（Neumann）边界条件有限元法的计算结果对比如图 5所示.

由图 5可以得到与3.2.1节相类似的结论.混合边界条件有限元法和有限元-无限元耦合法得到的结果十分相近，并且与边界积分方程法得到的结果拟合得很好；齐次边界条件有限元法在电源距边界较远时，即AB/2较小时与前述两种方法的结果相近，但是随着AB/2的增大，由于忽略边界积分而引起的误差逐渐加大.本例再次证明了有限元-无限元耦合法的实用性和可靠性，其精度与混合边界条件相当，优于相同计算范围下的齐次边界条件有限元法.

4.2.3 半空间下含一立方体

本模型来源于文献[30].电阻率1Ωm，边长为4m，顶部埋深3m的低阻立方体位于电阻率100Ωm的水平大地中.坐标原点位于立方体正上方的地表面上，采用偶极-偶极装置，点距1 m，仅布置一条测线，测线沿x轴方向分布，共21个极点.模型如图 6所示.

采用有限元-无限元耦合法计算，对模型范围为50m×50m×50 m和24 m×24 m×20 m的两个模型进行了计算，得到的视电阻率剖面图如图 7所示.

可以看出，计算范围不同的两个模型都反映出了“八字形”的低阻异常，并且与文献[30]中的计算结果吻合得很好，足以说明计算结果的正确性.图 7a模型在低阻立方体上方和下方视电阻率均趋于围岩电阻率100Ωm；对于图 7b这样小的模型，虽然在低阻立方体下方视电阻率与围岩实际电阻率相差较大，但也足以与低阻体引起的异常相区分.整体而言，有限元-无限元耦合法适用于偶极-偶极观测装置的正演模拟，并且可以在与探测区域相当的计算范围内得到合理的解.

5 结论

本文从三维直流电法边值问题出发，利用Galerkin加权余量法推导了基于总场的有限元方程，再引入Astley型映射无限元，将电位分布扩展到无限远处，从而在方程中不用再考虑与点源位置有关的混合边界条件.利用具有理论解的三维球体模型，对比了多种无限元形函数的计算效果，最终认定基于修正项（1-ξ）/2与二阶Lagrange插值多项式乘积的形函数在计算效率和精度上都占有优势.将其与有限元方法相耦合，形成了适用于三维直流电法正演模拟的有限元-无限元耦合法.若干计算实例表明，本文提出的有限元-无限元耦合法可以很好地解决三维直流电法正演模拟中半无限区域边界问题.无论对于单点AM装置，AMNB对称四极测深装置还是偶极-偶极装置，有限元-无限元耦合法都可以在近似测量区域的计算范围得到合理的解，且精度与混合边界条件解相当，优于相同计算范围下的齐次边界解，有利于在反演计算中减少计算量，加快计算速度.

在有限元-无限元耦合法中，iLUCG法被用于求解稀疏对称的线性方程组，一般经过几十次至几百次迭代即可得到收敛的解，单个方程组求解总时间一般不超过1s.对于具有多条剖面的偶极-偶极方案，如果采用直接解法求解方程，则在反演中可能更具速度优势.因为有限元-无限元耦合法已不需考虑边界条件，在电源位置改变时方程组的右端矩阵是不变的，在直接法中只需对全局矩阵进行一次分解，对于不同的右端项，通过简单的回代即可得到计算结果.但是直接法的主要缺点是矩阵分解耗时较多和内存占用巨大，应用于节点数较多的大规模正演问题还存在困难.

本文提出的有限元-无限元耦合法是基于总场的算法，相对于二次场算法而言，在电源点附近由于奇异性会产生较大的误差.但是总场算法更具普遍性，可以用于计算任意地形条件下的地电模型，电源点附近还可以利用非结构化网格加密的方法尽量减少误差，使得计算误差缩小到可接受的范围内.无限元的引入相当于将有限元计算范围扩展到了无限远，是对传统的混合边界条件的取代，因此有限元-无限元耦合法继承了有限元法的所有优点，并且解决了有限元法截断边界容易引起误差的缺点.

致谢

感谢国家自然科学基金对本研究课题提供资金资助.作者感谢任政勇博士（ETH Zurich，Switzerland）给予的卓有成效的讨论.