基于局部加密等级网格的2.5D直流电法有限元模拟

胡宏伶¹，肖晓²，潘克家^{3, 4}，汤井田²，谢维²

(1. 湖南师范大学 数学与计算机科学学院，湖南 长沙，410081;

2. 中南大学 地球科学与信息物理学院，有色金属成矿预测教育部重点实验室，湖南 长沙，410083;

3. 中南大学数学与统计学院, 湖南 长沙，410083;

4. 成都理工大学 油气藏地质及开发工程国家重点实验室, 四川 成都，610059)

摘要：合理截取半圆形计算区域，采取局部加密的l-等级网格，结合对称行索引存贮格式(CSR)及并行稀疏直接求解器PARDISO，提出一种高效、高精度的2.5D直流电阻率法有限元正演方案，并编制相应的Fortran程序，对具有解析解的3个典型地电模型进行计算与分析。研究结果表明：圆形截断边界不仅便于在径向方向上采取l-等级网格剖分，而且能大大简化有限元模拟中单元刚度矩阵的计算；结构化的等级网格避开了通常非结构化网格有限元计算时繁琐的网格剖分及总体刚度阵的集成过程，且能在不增大问题规模的前提下，显著提高2.5D直流电法正演源点附近的模拟精度；Intel MKL的PARDISO求解器能在普通PC机上5 s内求解电法正演有限元离散得到的100万阶稀疏线性方程组，可广泛用于各种地球物理正演问题。

关键词：直流电阻率；等级网格；有限元；局部加密；并行稀疏直接求解器；行索引存贮格式

中图分类号：P631             文献标志码：A         文章编号：1672-7207(2014)07-2259-09

Finite element modeling of 2.5D DC resistivity based on locally refined graded mesh

HU Hongling¹, XIAO Xiao², PAN Kejia^{3, 4}, TANG Jingtian², XIE Wei²

(1. College of Mathematics and Computer Science, Hunan Normal University, Changsha 410081, China;

2. Key Laboratory of Metallogenic Prediction of Nonferrous Metals, Ministry of Education,

School of Geosciences and Info-Physics, Central South University, Changsha 410083, China;

3. School of Mathematics and Statistics, Central South University, Changsha 410083, China;

4. State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation, Chengdu University of Technology,

Chengdu 610059, China)

Abstract: Truncating semi-circular computational domain reasonably, and combining the compressed sparse row (CSR) format and the parallel direct sparse solver (PARDISO), an efficient and high-precision finite element method for 2.5D DC resistivity modeling was proposed based on locally refined graded mesh, and the corresponding FORTRAN program was developed. The numerical results and analysis were given for three typical geoelectric models which have analytical solutions. The results show that the circular truncation boundary is not only easy to conduct a l-graded mesh partition in the radial direction, but also can significantly simplify the calculation of the element stiffness matrix in finite element simulations. The structured locally refined graded meshes can not only avoid the cumbersome processes of mesh generation and assembly of global stiffness matrix in finite element computations with usual unstructured meshes, but also significantly improve the accuracy near the source of 2.5D forward modeling of direct current method under the premise of not increasing the size of the problem. The Intel MKL PARDISO solver can solve, a sparse linear system of millions of unknowns arising from finite element discretization of DC modeling on an ordinary PC within 5 s, and can be widely used for a variety of the forward modeling problems in geophysics.

Key words: direct current resistivity; graded mesh; finite element; local refinement; PARDISO; compressed sparse row

自1971年Coggon^[1]将有限单元法(finite element method，FEM)应用到直流电阻率法数值模拟以来，FEM以其理论完备、边界处理能力强和通用性强等优点，在直流电法正演中得到了广泛应用^[2-18]。直流电法数值模拟问题具有如下特点：地下往往存在复杂的电性不均匀结构，电性差异可达几个数量级；研究区域具有无界性，需要进行截断边界处理；点电源处电位为无穷大，具有奇性。这些都将导致传统有限元法模拟精度及效率降低。为解决上述问题，提高有限元的计算效率及点源处的模拟精度，汤井田等^[19]提出一种3D直流电阻率有限元—无限元耦合数值模拟方法，有效地降低了截断边界的影响。为提高源点附近的精度，陈小斌等^[20]采用组合网格技巧，对源点附近进行局部加密，并利用有限元直接迭代算法求解线源频率域大地电磁正演问题。汤井田等^[21-24]采用非结构化网格的自适应有限元，通过加密奇点(点电源)处网格密度以减小有限元的计算误差，提高了有限元模拟的速度和精度。采用局部加密的非结构化网格，能提高奇点附近有限元的模拟精度，但网格剖分及总体刚度矩阵集成繁琐，编程较复杂，且网格剖分具有很大的不确定性。非结构化网格破坏了FEM 固有的许多超收敛结构，在提高局部精度的同时降低了FEM的整体收敛速度^[25]。另一种处理方式是采用较密的均匀剖分，并利用多网格法求解超大规模有限元方程。为解决通常几何多重网格法粗网格难以模拟复杂电性差异的问题，Moucha等^[26]提出一种电导率的6参数化表示法，并给出基于9点格式有限差分的多网格法，求解2D直流电法正演问题；鲁晶津等^[27]提出直流电阻率3D正演的代数多重网格方法。潘克家等^[28-29]引入外推瀑布式多网格法求解2.5D/3D直流电阻率法数值模拟问题，大大提高了直流电法正演的计算速度。Li等^[30]提出一种3D直流电法的多源系统有限元模拟方案。Wang等^[31]给出非结构化网格下各向异性3D直流电阻率法有限元模拟。在此，本文作者拟结合局部加密策略和结构化网格的优点，提出一种基于λ-等级网格的2.5D直流电法有限元模拟方案。

1  电阻率法2.5D有限元逼近

1.1  点源稳定电场边值问题

假设地下为2D构造，以电源点为坐标原点建立合适坐标系使z轴平行于构造走向，结合z方向的Fourier变换，点源波数域中电位函数在内满足如下2D变系数非齐次亥姆霍茨(Helmholtz)方程边值问题^[16]：

      

其中：为地表以下求解区域；为地表-空气界面；为截断边界；，为点源到边界上的距离；为电导率；I为电流；n为边界的外法向；为狄拉克函数；K₀和K₁分别为第二类零阶、一阶修正贝塞尔函数。

考虑到本文将采用局部加密的等级网格，截断边界取为如图1所示以原点为圆心、为半径的半圆，则电位u满足的椭圆边值问题可进一步简化为

       

图1  2D结构模型

Fig. 1  Model of 2D structure

其中：正常数。对微分方程边值问题，可导出其等价的变分问题：

       

1.2  电阻率法有限元逼近

1.2.1  分块等级网格剖分

直流电阻率法正演实际上为求解具有点源奇点的椭圆边值问题。陈传淼^[25]认为解决奇点问题最好的办法是局部加密网格法。

考虑半圆型求解区域 。首先，将分为k个子扇形，确保每个子扇形的圆心角。然后，在每个子扇形上，以r方向的剖分为基准，进行局部加密的三角形剖分。

设先将半径方向作均匀剖分 (k=0, 1, …, N)，步长。为了在O点附近局部加密，可取区间，并在其上取m个等级网格节点。

             

其中：等级网格参数＞1。显然，。而小区间的长度为

          

为了使整个网格逐步加密，而不是突然变小，可希望接近点的单元长度大致接近h₀，即，因此，要求。

重新排列剖分节点为 ，并对扇形进行局部加密的三角剖分。将半径为r_j的圆弧作j等分得到j+1个节点，并将它们连接为折线。然后，用直线依次连接相邻两圆弧上的节点，得到1个三角形网格。几个子区域上的网格就构成了所需的分块局部加密等级网格，见图2。为了简便，本文实际计算中直接取l=N，即在整个求解区域采取-等级网格。

1.2.2  单元分析

在单元e(逆时针编号依次为i, j, m)上进行三角形线性插值：

         

图2  l=2时的分块局部加密等级网格

Fig. 2  Piecewise locally refined graded mesh when l=2

其中：

    

；

；

             

首先将式中的积分分解成各单元e上的积分。若单元e内电导率为常数，利用式可得

      

其中：，为单元e的u值列向量；为3阶对称矩阵；

    

类似地，可得积分项

      

其中：

        

考虑式中的线积分，若单元e的边落在截断边界上，则

         

其中：

          

l为边的长度。

1.2.3  总刚合成

设区域剖分为n个节点，将矩阵K_1e，K_2e和K_3e相加可得3阶单元刚度矩阵K_e，并扩展成n阶方阵，然后叠加得

        

其中：；荷载向量

                

只有电源点对应位置非零(坐标原点为第1个网格节点)；U为未知的电位向量。

令式变分为0，得到线性方程组

KU=P                       

由式，及知总体刚度矩阵K为高度稀疏的对称正定矩阵，可利用PARDISO等直接稀疏求解器快速求解得到波数域电位。得到波数域中对应一系列不同波数k_i时的电位后，利用Fourier逆变换即可得到感兴趣的过点电源剖面(z=0)上的总电位：

          

其中：k_i和g_i分别为离散波数及相应的权系数。采用文献[32]中由差分进化算法^[33]计算得到离散波数，见表1。利用单源电场的电位，通过组合不同的采集装置，即可得到不同装置对应的视电阻率曲线，以用于电法勘探解释。

与通常采取的矩形计算区域相比，本文选取半圆型截断边界，简化了原问题中截断边界上的混合边界条件，使得贝塞尔函数不必代入式的边界积分进行计算，而只需计算和各1次，大大简化了刚度矩阵的计算过程。另一方面，半圆形求解区域非常适合如图2所示的结构化分块等级网格剖分，且能自动对奇点(点电源)局部加密，在保证精度的同时提高计算速度。若对矩形区域采取通常的结构化局部加密策略，即x方向和y方向分别采取由密到疏的剖分(如对数均匀或-等级网格)，则会出现如图3所示非常扁长的矩形单元，网格质量下降，严重影响有限元模拟精度^[25]。

图3  矩形域上局部加密网格

Fig. 3  Locally refined mesh on rectangular domain

1.3  压缩存储格式及总刚直接集成

1.3.1 行索引存贮格式(CSR)

此种存储格式按照行的顺序，依次将对称稀疏矩阵的上三角非零元素存放在一段连续存储空间上。假设A为n阶对称方阵，m表示其上三角部分非零元的个数，数组下标从1开始，则这种存储结构由以下3个数组构成：

双精度数组a，即按行存储A的上三角部分的非零元素，长度为m；

整型数组J_a，即依次存储A中元素的列号，长度为m；

整型数组I_a，即存储A每行第1个非零元在数组a中的位置，长度为n+1。其中，I_a(n+1)=m+1，即第i行非零元个数为I_a(i+1)-I_a(i)。

图4所示为1个对称CSR格式存储的示例。此种存储格式可用于Intel MKL的PARDISO直接稀疏求解器，并且存储效率非常高。对图3所示的矩形网格采用矩形双线性元，有限元离散形成的总体刚度阵每行最多5个非零元。而对图2所示的等级网格，采用三角形线性元，总体刚度阵每行最多4个非零元。

1.3.2  CSR格式总刚直接集成

下面基于CSR存储格式，集成全局刚度矩阵，即将单元e中第i号和第j号节点对应的“单胞”累加到总刚的位置。

表1  离散波数及相应的权系数

Table 1  Discrete wavenumbers and corresponding weights

图4  CSR格式存储示例

Fig. 4  Example of CSR storage format

分析刚度矩阵的组装过程得知：只有2个节点相互关联(出现在同一单元中)，才会对总体刚度矩阵有贡献。在总刚K的第i行，非零元的列标正好为与节点i关联的节点编号。为此，为每个节点i设计1个名为“关联节点”的数据结构，用于存储非零元的列标。

整数N_i，为与节点i关联的节点个数；

整型数组C_i，为按升序依次存储关联节点编号。

考虑到总体刚度矩阵的对称性，实际上关联节点中不必存储小于i的节点编号。图5所示为1个基于2D三角形剖分的关联节点示例。用l _ij表示“单胞” 被累加到a中的位置，d_ij表示j在数组C_i中的序号，则可按如下规则直接集成CSR格式的总体刚度阵：

           

其中：

          

图5  关联节点示意图

Fig. 5  Schematic diagram of associated node

2  数值计算结果与分析

本文程序在Windows 7(Intel Core^TM i5-760，2.8 GHz，4.0 GB内存)下运行，软件平台为Intel Parallel Studio XE 2011，包含Intel Visual Fortran Compiler XE 12.0以及Intel Math Kernel Library 10.3。

首先考虑扇形区域上如下拉普拉斯方程边值问题：

             

其中：边界函数g由精确解确定，精确解在原点处不可导，具有奇性。采用如前所述的局部加密的等级网格，利用有限元方法求解椭圆边值问题，计算结果如表2所示(其中：N为径向剖分数；误差表示数值解与精确解之间的最大绝对误差)。从表2可看出，均匀网格下对此类奇异解问题，有限元只有半阶收敛；网格参数增加1，收敛阶增加半阶；固定N，随着增加，精度显著提高；取4时已达到线性有限元最高阶即2阶收敛，故没必要进一步增大。另外，此例说明对此类奇异解问题，均匀剖分收敛非常缓慢，而采用局部加密的等级网格可显著提高有限元解的精度，即使=2时128等级剖分的计算精度(4.69×10^-6)已与均匀1 024剖分时的精度(4.20× 10^-6)相当，而=4时128等级剖分的误差仅为2.45× 10^-7，也远比均匀1 024剖分时的4.20×10^-6小；同时也验证了本文方法的可行性和有限元程序的正确性。

表2  等级网格有限元解收敛性分析

Table 2  Convergence analysis of finite element solutions for graded meshes

需要指出的是：以上模型问题为调和方程边值问题，没有考虑物性的不均匀性以及附加源的影响；同时，精确解的奇性也比较弱，只在原点导数不存在，而函数仍是连续的。然而，在直流电法数值模拟中，地下存在复杂的电性不均匀结构，电性差异甚至可达几个数量级；同时电法正演问题为无界问题，需要进行截断边界处理，且在点电源处电位函数本身就不连续，为无穷大。这些都将导致有限元模拟精度及效率下降，因此，不能期望其达到如表2所示模型问题的计算精度。

下面对3个具有解析解的水平层状及垂直电性界面模型，均取截断半径 m，即计算区域为。网格采用局部加密的-等级网格，参数分别取为1，2和3，取1时的等级网格即退化为r方向均匀剖分。径向r剖分数N取为800，网格单元数为3×800²=1 920 000，网格节点数为962 001。

基于行索引稀疏存储模式，利用Intel MKL的PARDISO直接稀疏求解器求解有限元方程。因为基于Cholesky分解的直接求解器计算时间主要取决于问题的规模及系数矩阵的稀疏结构，故每次求解962 001阶有限元方程的时间差不多，约为4.5 s；而进行1次完整的正演计算(需要解5次方程组)，总耗时约为24.5 s，由此可知正演程序其它部分(包括5次单元刚度阵的计算及总体刚度矩阵的集成) 仅耗时2 s。这是因为局部加密的-等级网格为结构化网格，生成的刚度矩阵稀疏结构简单，且每行最多7个非零元，实际上考虑到对称性仅需存储4个元素，不仅减小计算机内存，而且可大大简化总体刚度矩阵的集成。而采取半圆形求解区域，简化了原边值问题截断边界上的混合边界条件，使得贝塞尔函数不必参与单刚的计算，提高了总刚的计算速度。

2.1  模型一^[28]

2层水平地层：第1层电阻率为10Ω·m，厚度为10 m；第2层电阻率为100Ω·m。采用单极供电，供电电极位于坐标原点。图6所示为两极装置得到的视电阻率曲线，其中精确解为线性滤波法计算结果^[34]。图7所示为主剖面上的相对误差分布。从图7可看出：相对误差关于y轴对称，且在y=10左右有明显突变，反映出介质水平分界面的位置。

图6  两层地电模型计算结果

Fig. 6  Numerical results for two-layer geoelectric model

图7  主剖面上的相对误差

Fig. 7  Relative errors on main cross-section

从图6可以看出：3种不同的网格剖分，在远离电源点时(r＞3)都能够得到非常精确的结果，相对误差控制在0.5%以内；径向采取均匀剖分时(=1)，数值解在点源附近的误差较大，超过4%；局部加密等级网格在不增加计算工作量的前提下，显著提高了奇点附近的精度，且2种不同的等级网格(=2, 3)得到几乎完全相同的结果。经计算，基于等级网格有限元解最大相对误差为0.65%，而平均相对误差仅为0.30%。从计算时间和模拟精度2个方面考虑，此例计算结果都与文献[28]中采用1 600×1 600(2 563 201个节点，2 560 000个单元)均匀矩形剖分时，外推瀑布式多网格法的计算结果相当。文献[28]中计算时间为28 s，最大相对误差0.87%，平均相对误差0.22%。

2.2  模型二^[5]

K型3层层状地层：第1层电阻率为50Ω·m，厚度为5 m；第2层电阻率为100Ω·m，厚度为10 m；第3层电阻率为20Ω·m。数值模拟结果见图8。从图8可看出：采用等级网格的有限元法能在200 m范围内得到高精度的计算结果，平均相对误差为0.29%；而200 m以外相对误差显著增加，这是因为Fourier离散波数是基于均匀地层精确解利用200 m内的极距序列最优化选取的，适用范围为1~200 m^[29]。除此之外，可得到与模型一完全类似的结论：(1) 局部加密等级网格可显著提高奇点附近的模拟精度(x=1附近相对误差从4.60%显著减为0.10%；(2) 网格参数=2和=3对应的视电阻率误差曲线几乎重合，即对此类电法正演问题网格参数取为2即可。

2.3  模型三^[21]

2D垂直典型界面，左边介质电阻率为1Ω·m，右边介质电阻率为100Ω·m；单位点电源位于左边介质上，距界面5 m。基于局部加密等级网格的有限元法得到的视电阻率及相对误差如图9所示，其中图9(b)中虚线表示分界面的位置。从图9可看出：采用均匀网格剖分，奇点及分界面附近的相对误差比较大；而等级网格相对误差在1~200 m范围内都比较小，控制在0.70%以内，与文献[21]中基于非结构化网格的自适应有限元计算结果精度相当，文献[21]中给出的误差限为0.69%。在介质分界面处(x=5 m)，这2种网格剖分下视电阻率的相对误差均有体现，具有明显的突变。值得注意的是：利用等级网格，x=1处的相对误差从3.94%骤降到0.12%，这进一步说明了等级网格的确能显著提高奇点附近的精度。

图8  3层地电模型计算结果

Fig. 8  Numerical results for three-layer geoelectric model

图9  垂直电性界面模型计算结果

Fig. 9  Numerical results of vertical electrical interface

3  结论

(1) 在无穷远处截取圆形边界，大大简化了2.5D直流电法模拟电位满足的椭圆边值问题及单元刚度矩阵的计算。

(2) 局部加密等级网格能在不改变问题规模的前提下，显著提高2.5D直流电法正演的有限元模拟精度。

(3) 对简单的模型问题，随着等级网格参数l的增大(小于4)，奇异解的计算精度显著增大，但对直流电法正演影响不大，取为2即可。

(4) 局部加密等级网格为结构化网格，避开了通常非结构化网格计算时繁琐的网格剖分及总刚集成过程，有限元程序实现简便。

(5) Intel MKL的PARDISO为一种高效、并行直接稀疏求解器，能在普通PC机上5 s内求解直流电阻率法正演有限元离散得到的100万阶稀疏线性方程组，可用于各种地球物理正演问题的计算。

参考文献：

[1] Coggon J H. Electromagnetic and electrical modeling by the finite element method[J]. Geophysics, 1971, 36(1): 132-145.

[2] Dey A, Morrison H F. Resistivity modeling for arbitrarily shaped two dimensional structures[J]. Geophysical Prospecting, 1979, 27(1): 106-136.

[3] 底青云, 王妙月. 稳定电流场有限元法模拟研究[J]. 地球物理学报, 1998, 41(2): 252-260.

DI Qingyun, WANG Miaoyue. The real-like 2D FEM modeling research on the field characteristics of direct electric current field[J]. Chinese Journal of Geophysics, 1998, 41(2): 252-260.

[4] LI Changwei, XIONG Bin, QIANG Jianke, et al. Multiple linear system techniques for 3D finite element method modeling of direct current resistivity[J]. Journal of Central South University, 2012, 19(2): 424-432.

[5] 李长伟, 阮百尧, 吕玉增, 等. 点源场井-地电位测量三维有限元模拟[J]. 地球物理学报, 2010, 53(3): 729-736.

LI Changwei, RUAN Baiyao, L Yuzeng, et al. Three- dimensional FEM modeling of point source borehole-ground electrical potential measurements[J]. Chinese Journal of Geophysics, 2010, 53(3): 729-736.

[6] Li Y G, Spitzer K. Three-dimensional DC resistivity forward modelling using finite elements in comparison with finite-difference solutions[J]. Geophysical Journal International, 2002, 151(31): 924-934.

[7] 罗延钟, 孟永良. 关于用有限单元法对二维构造作电阻率法模拟的几个问题[J]. 地球物理学报, 1986, 29(6): 613-621.

LUO Yanzhong, MENG Yongliang. Some problems on resistivity modeling for two-dimensional structures by the finite element method[J]. Chinese Journal of Geophysics, 1986, 29(6): 613-621.

[8] Mundry E. Geoelectrical model calculations for two-dimensional resistivity distributions[J]. Geophysical Prospecting, 1984, 32(1): 124-131.

[9] 强建科, 罗延钟. 三维地形直流电阻率有限元法模拟[J]. 地球物理学报, 2007, 50(5): 1606-1613.

QIANG Jianke, LUO Yanzhong. The resistivity FEM numerical modeling on 32D undulating topography[J]. Chinese Journal of Geophysics, 2007, 50(5): 1606-1613.

[10] 阮百尧, 徐世浙. 电导率分块线性变化二维地电断面电阻率测深有限元数值模拟[J]. 地球科学: 中国地质大学学报, 1998, 23(3): 303-307.

RUAN Baiyao, XU Shizhe. FEM for modeling resistivity sounding on 2D geoelectric model with line variation of conductivity within each block[J]. Earth Science: Journal of China University of Geoscience, 1998, 23(3): 303-307.

[11] 阮百尧, 熊彬. 电导率连续变化的三维护电阻率测深有限元模拟[J]. 地球物理学报, 2002, 45(1): 131-138.

RUAN Baiyao, XIONG Bin. A finite element modeling of three-dimensional resistivity sounding with continuous conductivity[J]. Chinese Journal of Geophysics, 2002, 45(1): 131-138.

[12] TANG Jingtian, WANG Feiyan, XIAO Xiao, et al. 2.5D DC resistivity modeling considering flexibility and accuracy[J]. Journal of Earth Science, 2011, 22(1): 124-130.

[13] WU Xiaoping. A 3D finite-element algorithm for DC resistivity modeling using shifted incomplete cholesky conjugate gradient method[J]. Geophysical Journal International, 2003, 154(3): 947-956.

[14] 吴小平, 汪彤彤. 利用共轭梯度算法的电阻率三维有限元正演[J]. 地球物理学报, 2003, 46(3): 428-432.

WU Xiaoping, WANG Tongtong. A 3D finite-element resistivity forward modeling using conjugate gradient algorithm[J]. Chinese Journal of Geophysics, 2003, 46(3): 428-432.

[15] 熊彬, 阮百尧. 电位双二次变化二维地电断面电阻率测深有限元数值模拟[J]. 地球物理学报, 2002, 45(2): 285-295.

XIONG Bin, RUAN Baiyao. A numerical simulation of 2D geoelectric section with biquadratic change of potential for resistivity sounding by the finite element method[J]. Chinese Journal of Geophysics, 2002, 45(2): 285-295.

[16] 徐世浙. 地球物理中的有限单元法[M]. 北京: 科学出版社, 1994: 131-205.

XU Shizhe. The finite element method in geophysics[M]. Beijing: Science Press, 1994: 131-205.

[17] 徐世浙, 刘斌, 阮百尧. 电阻率法中求解异常电位的有限单元法[J]. 地球物理学报, 1994, 37(S2): 511-515.

XU Shizhe, LIU Bin, RUAN Baiyao. The finite element method for solving anomalous potential for resistivity surveys[J]. Chinese Journal of Geophysics, 1994, 37(S2): 511-515.

[18] Zhou B, Greenhalgh S A. Finite element three-dimensional direct current resistivity modeling: Accuracy and efficiency considerations[J]. Geophysical Journal International, 2001, 145: 679-688.

[19] 汤井田, 公劲喆. 三维直流电阻率有限元-无限元耦合数值模拟[J]. 地球物理学报, 2010, 53(3): 717-728.

TANG Jingtian, GONG Jinzhe. 3D DC resistivity forward modeling by finite-infinite element coupling method[J]. Chinese Journal of Geophysics, 2010, 53(3): 717-728.

[20] 陈小斌, 胡文宝. 有限元直接迭代算法及其在线源频率域电磁响应计算中的应用[J]. 地球物理学报, 2002, 45(1): 119-130.

CHEN Xiaobin, HU Wenbao. Direct iterative finite element (DIFE) algorithm and its application to electromagnetic response modeling of line current source[J]. Chinese Journal of Geophysics, 2002, 45(1): 119-130.

[21] 汤井田, 王飞燕, 任政勇. 基于非结构化网格的2.5D直流电阻率自适应有限元数值模拟[J]. 地球物理学报, 2010, 53(3): 708-716.

TANG Jingtian, WANG Feiyan, REN Zhengyong. 2.5D DC resistivity modeling by adaptive finite-element method with unstructured triangulation[J]. Chinese Journal of Geophysics, 2010, 53(3): 708-716.

[22] 任政勇, 汤井田. 基于局部加密非结构化网格的三维电阻率法有限元数值模拟[J]. 地球物理学报, 2009, 52(10): 2627-2634.

REN Zhengyong, TANG Jingtian. Finite element modeling of 3D DC resistivity using locally refined unstructured meshes[J]. Chinese Journal of Geophysics, 2009, 52(10): 2627-2634.

[23] REN Zhengyong, TANG Jingtian. 3D direct current resistivity modeling with unstructured mesh by adaptive finite-element method[J]. Geophysics, 2010, 75(1): 7-17.

[24] TANG Jingtian, REN Zhengyong, WANG Feiyan, et al. 3D direct current resistivity forward modeling by adaptive multigrid finite element method[J]. Journal of Central South University of Technology, 2010, 17(3): 587-592.

[25] 陈传淼. 科学计算概论[M]. 北京: 科学出版社, 2007: 120-262.

CHEN Chuanmiao. Introduction to scientific computing[M]. Beijing: Science Press, 2007: 120–262.

[26] Moucha R, Bailey R C. An accurate and robust multigrid algorithm for 2D forward resistivity modelling[J]. Geophysical Prospecting, 2004, 52(3): 197-212.

[27] 鲁晶津, 吴小平, Spitzer K. 直流电阻率三维正演的代数多重网格方法[J]. 地球物理学报, 2010, 53(3): 700-707.

LU Jingjin, WU Xiaoping, Spitzer K. Algebraic multigrid method for 3D DC resistivity modeling[J]. Chinese Journal of Geophysics, 2010, 53(3): 700-707.

[28] 潘克家, 汤井田, 胡宏伶, 等. 直流电阻率法2.5D正演的外推瀑布式多重网格法[J]. 地球物理学报, 2012, 55(8): 2769-2778.

PAN Kejia, TANG Jingtian, HU Hongling, et al. Extrapolation cascadic multigrid method for 2.5D direct current resistivity modeling[J]. Chinese Journal of Geophysics, 2012, 55(8): 2769-2778.

[29] PAN Kejia, TANG Jingtian. 2.5D and 3D DC resistivity modelling using an extrapolation cascadic multigrid method[J]. Geophysical Journal International, 2014, 197(3): 1459-1470.

[30] LI Changwei, XIONG Bin, QIANG Jianke, et al. Multiple linear system techniques for 3D finite element method modeling of direct current resistivity[J]. Journal of Central South University, 2012, 19(2): 424-432.

[31] Wang W, Wu X, Spitzer K. Three-dimensional DC anisotropic resistivity modelling using finite elements on unstructured grids[J]. Geophysical Journal International, 2013, 193(2): 734-746.

[32] 潘克家, 汤井田. 2.5D直流电法正演中Fourier逆变换离散波数的最优化选取[J]. 中南大学学报(自然科学版), 2013, 44(7): 2819-2826.

PAN Kejia, TANG Jingtian. Optimized selection of discrete wavenumbers for inverse Fourier transform in 2.5D DC resistivity modeling[J]. Journal of Central South University (Science and Technology), 2013, 44(7): 2819-2826.

[33] 潘克家, 陈华, 谭永基. 基于差分进化算法的核磁共振T₂谱多指数反演[J]. 物理学报, 2008, 57(9): 5956-5961.

PAN Kejia, CHEN Hua, TAN Yongji. Multi-exponential inversion of T₂ spectrum in NMR based on differential evolution algorithm[J]. Acta Physica Sinica, 2008, 57(9): 5956-5961.

[34] Guptasarma D, Singh B. New digital linear filters for Hankel J₀ and J₁ transforms[J]. Geophysical Prospecting, 1997, 45(5): 745-762.

(编辑  陈灿华)

收稿日期：2013-07-30；修回日期：2013-09-18

基金项目：国家自然科学基金资助项目(41204082, 11301176, 41104071)；中国博士后科学基金特别资助项目(2013T60781)；高等学校博士学科点专项科研基金资助项目(20120162120036)；成都理工大学油气藏地质及开发工程国家重点实验室基金资助项目(PLC201305)

通信作者：肖晓(1981-)，男，湖北丹江口人，博士，从事电磁法正反演及电磁信号处理研究；电话：13786149397；E-mail: csuxiaox@csu.edu.cn