麻昌英, 柳建新, 郭荣文, 等. 2018. 耦合有限单元法扩边的直流电阻率勘探无单元Galerkin法正演. 地球物理学报, 61(6): 2578-2588, doi: 10.6038/cjg2018L0492 . Ma C Y, Liu J X, Guo R W, et al. 2018. Element-free Galerkin forward modeling of DC resistivity using a coupled finite element method with extended the boundaries. Chinese J. Geophys . (in Chinese), 61(6): 2578-2588, doi: 10.6038/cjg2018L0492 . 2018-03-12 收修定稿 第一作者简介: 麻昌英, 男, 1988年生, 博士研究生, 主要从事电磁法正反演和应用研究.E-mail: [email protected] 摘要 ：本文分析了目前直流电阻率正演模拟中的无单元Galerkin法（EFGM）和有限单元法（FEM）的优缺点，针对采用第一类边界条件需要足够大的计算域时EFGM计算成本高的问题，在计算域外围区域采用FEM扩边，提出了直流电阻率的无单元Galerkin-有限单元耦合法（EFG-FE）.采用具有Kronecker delta函数性质的径向基点插值法（RPIM）构造EFGM形函数，在外围区域将EFGM与FEM直接耦合，无需其他处理手段，消除了传统EFGM与FEM耦合中存在的界面耦合困难.EFG-FE将模型计算域分割为EFGM区域和FEM区域，模型核心区域采用EFGM计算，发挥EFGM灵活性、适应性强和高精度的优点，使得模型建立简单方便，对任意复杂地电模型适应性强，同时获得高精度模拟结果.在模型计算域外围采用快速扩展的FEM单元网格进行剖分，利用其数值稳定性和高效性，使用少量FEM节点和单元网格将计算域大范围扩大满足第一类边界条件，同时不大幅增加计算成本，进而提高计算效率.最后，通过不同正演方法的模型算例的模拟结果对比，验证了本文提出的EFG-FE有效可行，其模拟结果具有很高的模拟精度，且相比于采用第三类边界条件的EFGM提高了计算效率，具有更好的模拟性能.

Key words : Element-free Galerkin-finite element coupling method RPIM shape function Element-free Galerkin method Finite element method Direct current resistivity

无单元Galerkin法(Element-Free Galerkin Method, EFGM)是一种基于离散节点的数值模拟方法, 其利用局部支持域内的节点空间信息构造形函数, 通过Galerkin全局弱式将计算域上的弱式积分方程转化为离散系统方程, 具有模拟精度高, 预处理简单, 节点可根据需要任意布置, 适应性强等特点.起初, 随着滑动最小二乘(Moving Least Squares, MLS)被用于构造形函数, 使得无网格法(Meshless method)形函数的连续性问题而得以解决( Lancaster and Salkauskas, 1981 ), 在此基础上 Belytschko等(1994a) 改进了强加本质边界条件及数值积分处理等问题并提出EFGM, 随后EFGM在各个研究领域被广泛关注( Peng et al., 2011 ; Hajiazizi and Bastan, 2014 ; Hadinia and Jafari, 2015 ).EFGM已被证明是一种非常优秀的数值模拟方法, 在模拟精度、灵活性、便利性等方面均优于有限单元法(Finite Element Method, FEM), 不仅能应用到FEM的应用领域, 且可解决FEM由于网格的存在而难以处理的问题, 如场剧烈变化问题, 大变形问题, 复杂边界和形体问题等( Belytschko et al., 1994b ; Li et al., 2012 ).在地球物理勘探领域方面, EFGM已经得到了一些研究和应用, 贾晓峰等(2006) 对地震波场应用EFGM进行了正演模拟研究; 王月英等(2007) 探讨了EFGM应用于地震波正演模拟中的各种影响因素; 冯德山等(2013) 采用MLS构造形函数, 实现了探地雷达二维EFGM正演模拟; 严家斌和李俊杰(2014) 也采用MLS构造形函数, 将EFGM应用于大地电磁法二维模拟, 计算了水平地形下大地电磁二维模型响应; 嵇艳鞠等(2016) 利用EFGM计算了各向异性的2D大地电磁响应; 麻昌英等(2017) 实现了直流电阻率EFGM正演模拟, 并与FEM进行了比较, 展现了EFGM的高精度和高灵活性特点. 然而, EFGM采用高斯积分计算积分, 由于EFGM的近似函数一般为有理函数而不是多项式, 为了保证模拟精度须使用较多的积分点和较高阶的高斯积分点( Dolbow and Belytschko, 1999 )进行计算, 而EFGM计算量与积分点数成正比, 使得EFGM计算量大, 计算效率低.另外, 在直流电电阻率正演模拟中, 当采用第一类边界条件时, 计算域需要足够大, 此时在远离场源和异常体的计算域外围区域场值变化平缓, EFGM虽然对任意分布的节点具有很强的适应性, 但仍要求节点合理分布, 节点的不合理分布可能造成EFGM形函数构造失败, 因此在外围区域需要较多的节点覆盖以保证节点的合理分布, 这导致EFGM计算成本明显增加.虽然FEM由于单元网格的不可或缺而受到许多约束, 但FEM具良好的数值稳定性, 单元的形状及其分布变化情况不会造成模拟计算失败, 在某些特殊区域可充分利用该特点发挥FEM的性能, 例如在远离场源的场值变化平缓的区域, 可快速延伸单元扩大计算域以满足第一类边界条件, 快速延伸的单元虽然形状细长, 但FEM仍可顺利进行计算, 同时对模拟精度影响很小.针对上述直流电阻率EFGM正演中采用第一类边界条件时存在的问题, 本文尝试在外围区域使用FEM进行计算, 并尝试将EFGM与FEM直接耦合, 耦合法可将不同的数值模拟方法相结合, 发挥单一方法的优点的同时又能避免各自的缺点.EFGM与FEM的耦合已经得到了相关的研究和应用( Belytschko et al., 1995 ; Liu et al., 2006 ; Hadinia et al., 2016 ), 但在传统EFGM与FEM的耦合中, 由于传统EFGM基于MLS插值构造的形函数不具有Kronecker delta函数性质, 不能满足边界条件, 所以若将各自得到的代数方程直接耦合, 则会造成两个子域在耦合界面上场值不连续, 使得EFGM与FEM的耦合存在困难, 因此如何处理EFGM区域和FEM区域的耦合界面关系就成为EFGM与FEM耦合需要解决的首要问题, 在这方面已经展开了相关研究并取得了一定的成果( Kumar et al., 2014 ; Pathak et al., 2015 ; Samira et al., 2017 ), 但这些研究成果均为使用特殊处理手段实现EFGM与FEM耦合, 增加了算法复杂度. 本文从EFGM形函数方面进行研究, 尝试采用不同的形函数构造方法来解决EFGM与FEM耦合困难问题.地球物理研究领域中, 麻昌英等(2017) 利用径向基点插值法(Radial Point Interpolation Method, RPIM)构造EFGM形函数, 获得了具有Kronecker delta函数性质的RPIM形函数; 李俊杰和严家斌(2015) 初步展开了大地电磁法的EFGM与FEM耦合法的研究, 提出了有限元-径向基点插值法, 获得了良好的效果.RPIM以径向基函数(Radial Base Function, RBF)为基础, RBF具有优良数学特性使得RPIM形函数能满足EFGM形函数需要满足一些基本要求, 如适应性、稳定性、一致性、紧支性、相容性、Kronecker delta函数性质等( Liu and Gu, 2001 ).本文采用RPIM插值构造EFGM形函数解决EFGM与FEM耦合困难问题, 在计算域外围区域将EFGM与FEM直接耦合.同时, 针对采用第一类边界条件EFGM存在的问题, 基于矩形单元剖分的FEM, 采用快速扩展的单元网格覆盖远离场源的计算域外围区域, 以扩大计算域使得第一类边界条件得到满足.在此基础之上, 提出直流电阻率的无单元Galerkin-有限单元耦合法(EFG-FE).通过不同模拟方法对层状模型测深和均匀半空间点电源电位进行模拟对比, 验证了EFG-FE算法的有效性, 表明了EFG-FE可方便地通过外围FEM区域扩大计算域来满足第一类边界条件, 获得高精度的模拟结果.最后, 基于相同的EFGM计算域模型, 将EFG-FE与采用第三类边界条件的EFGM进行对比研究, 结果表明本文提出的EFG-FE相比于采用第三类边界条件的EFGM提高了计算效率, 具有更好的模拟性能. 1 直流电阻率边值问题 采用第一类边界条件, 2.5维稳定电流场波数域总电位函数 U ( x , λ )满足的边值问题及其对应的变分问题分别表示为( 阮百尧, 2001 )：

式中, Ω 为计算域, Γ _S 为地表边界, Γ _∞ 为截断边界, σ =( x , z )为介质电导率, U ( λ , x , z )为波数域电位, λ 为波数, I ₀ 为电流, δ 为Kronecker delta函数, x =[ x , z ] ^T 为 Ω 上的任意一点, A 为场源点位置, 为二维哈密顿算子, n 为外法线单位向量.直流电阻率二维模型如 图 1 所示. 更高阶导数同理可得.由(13)—(14)式可知, 为了获得RPIM形函数及其导数均需要独立求解方程, 这是导致EFGM计算效率较低的重要因素. 3 直流电阻率的无单元-有限单元耦合法 如 图 2 所示, 假设计算域 Ω 被分割为EFGM区域 Ω ₁ 和FEM区域 Ω ₂ , Ω = Ω ₁ ∪ Ω ₂ , Ω ₁ 中包含 N ₁ 个EFGM节点, Ω ₂ 中包含 N ₂ 个FEM节点.其中EFGM区域 Ω ₁ 边界上的 N ₃ 个节点为EFGM和FEM共用节点, 则计算域 Ω 中节点总数为 N = N ₁ + N ₂ － N ₃ . Ω ₁ 为模型的核心区域, 采用EFGM进行正演计算. Ω ₂ 为外围扩边区域, 采用FEM进行正演计算, 采用快速延伸网格扩大计算域, 以满足第一类边界条件.特别指出的是, 如 图 2 所示EFGM区域被FEM区域包围, 与截断边界 Γ _∞ 无交集, 因此EFGM区域内无需做边界计算. 3.1 无单元Galerkin法 EFGM的求解过程如 图 3 所示.EFGM中首先使用一组总数为 N ₁ 的节点离散EFGM区域 Ω ₁ , 然后使用一组背景单元将 Ω ₁ 覆盖, 接着在每一个背景单元中布置高斯积分点计算积分, 将所有背景单元的积分求和得到 Ω ₁ 上的积分.在EFGM计算过程中, 对每一个高斯积分点需要独立构造局部支持域 Ω _q , 并利用支持域内所包含的节点信息构造无单元法形函数.特别指出的是, 在EFGM中节点与背景单元是相互独立的, 节点用于模拟地电模型为EFGM的关键模型信息, 背景单元用于计算积分, 与FEM中的单元概念是完全不同的.在保证数值积分精度的情况下, 如果可以直接在 Ω ₁ 上布置高斯积分点, 则可使得EFGM更加简便, 但通常这样很难控制积分计算精度, 因此引入背景单元便于布置高斯积分点进行积分计算. 其中 , a 和 b 为矩形单元的长和宽.矩形单元四个节点坐标分别表示为 x _l =( x _l , z _l ), l =1, 2, 3, 4, 下标1、2、3和4为矩形单元顶点对应的节点局部编号; U _e =( U ₁ , U ₂ , U ₃ , U ₄ ) ^T 为矩形单元节点上的波数域电位列向量. 将FEM区域 Ω ₂ 所有矩形形单元的子方程组(19)式组装得到直流电阻率FEM方程组为 式中 K ^FE 为由全部矩形单元子域 Ω _e 的 K _e 相加组成的 N ₂ × N ₂ 维的FEM刚度系数矩阵, U ^FE 为FEM区域内所有矩形单元节点上的波数域电位子向量 U _e 组合成的1× N ₂ 维列向量, S 为1× N ₂ 维的FEM方程右端项列向量. 3.3 无单元Galerkin法与有限单元法耦合 在无单元Galerkin-有限单元耦合法中, 耦合界面 Γ _int 上两种方法需满足电位相容性条件( Gu and Liu, 2005 ; Liu et al., 2016 )，公式为 传统的无单元法Galerkin法基于MLS形函数, 由于MLS形函数不具有Kronecker delta函数性质, 因此不具有相容性, 若将各自得到代数方程直接耦合, 则会造成两个子域在耦合界面上场值不连续, 产生较大误差.通常, 在传统EFGM与FEM耦合中在两种方法的子域之间设置连接区域，并引入斜坡函数等技术以克服相容性问题( Liu et al., 2006 ; Kumar et al., 2014 ; Pathak et al., 2015 ; Hosseini et al., 2017 ), 这些处理手段增加了耦合难度和算法的复杂性.本文采用具有Kronecker delta函数性质的RPIM形函数, 使得(26)式自然得到满足, 因此可以克服传统MLS形函数不具有Kronecker delta函数性质而产生的耦合困难.在耦合方法中, 相容性条件(26)式是相容性条件和自然边界条件两项要求中最重要的偶合条件, 其必须得到满足, 即使自然边界条件(27)式在一些耦合方法中不能完全满足, 也可得到良好的效果, 如果两种方法的形函数在耦合界面上满足Kronecker delta函数性质, 则这两种数值方法在耦合界面上可直接进行耦合( Gu and Liu, 2005 ).目前, Liu等(2016) 将基于RPIM形函数的局部Petrov-Galerkin无单元-有限单元耦合法应用于电磁计算领域中, 展示了两种方法子域耦合界面上RPIM形函数与FEM形函数的相容性, 并采用两种方法直接耦合的方式简化了耦合过程, 并获得了良好的效果. 事实上, EFGM和FEM主要不同之处在于采用不同的插值方法, 但采用单元插值的FEM和采用局部支持域插值的EFGM均为分片插值, 它们的形函数均满足紧支性, 即FEM形函数在单元外为零, EFGM形函数在局部支持域外为零.因此, 当EFGM形函数满足Kronecker delta函数性质, 且局部支持域限定在EFGM区域内时, 则与FEM单元连接的EFGM局部支持域可视为FEM的一个连续插值单元(广义单元). 贾亮等(2007) 等基于广义单元的概念, 将FEM单元和EFGM局部支持域视为广义单元, 利用在每一个单元内假设的形函数来分片地表示求解域上待求的未知场函数, 将FEM和EFGM直接耦合, 获得了良好的结果. 张琰等(2008) 利用RPIM方法近似函数满足插值条件, 将无单元法与有限单元法直接耦合, 应用于应力分析获得了良好的效果. 基于上述分析结合直流电阻率地质模型, 本文在模型核心区域采用EFGM计算, 发挥EFGM的灵活性、适应性强和高效性等特点.选择耦合面两侧介质电阻率相同的区域分割FEM区域和EFGM区域, 使得计算域内部介质界面上满足自然边界条件(27)式, 同时将EFGM局部支持域限定在EFGM区域, 将两种方法子域直接耦合.因此, 将(15)和(21)式根据节点整体编号组装在一起得到直流电阻率EFG-FE的离散方程为 式中 K = K ^EFG + K ^FE , U = U ^EFG + U ^FE , P = F + S , K 为 N × N 维的EFG-FE刚度系数矩阵.为方便两种方法的计算, 首先在两种方法子域中对节点分别采用局部编号.在FEM区域 Ω ₂ , 节点局部编号按常用的编号方法进行编号, 如 图 4a 所示, N _z 表示垂向方向FEM节点数.在EFGM区域往往采用不规则分布的场节点, 节点局部编号方法如 图 4b 所示, 先按 X 方向从左到右排列, 再按Z方向由上到下排列, 事实上除了节点分布不规则外, 这样的编号方式与FEM区域采用的节点编号方式是一样的.最后, 在全局域上全部节点(包含FEM节点和EFGM节点)也按先 X 方向从左到右排列, 再按 Z 方向由上到下排列进行整体编号, 这有助于将刚度系数矩阵非零元素带宽减小, K 与有限单元法形成的刚度系数矩阵相似, 它是稀疏对称正定矩阵, 使用定带宽的方式储存稀疏刚度矩阵( 徐世哲, 1994 ), 减少零元素的储存.采用MKL库中开源求解器Pardiso求解方程组(28)式, 得到计算域 Ω 上 N 个节点的波数域电位, 通过反傅立叶变换可得到节点的空间域电位( 潘克家等, 2012 ). 4.1 层状模型 分别应用EFG-FE和基于矩形和线性插值的FEM对如 图 5 所示的三层地电模型进行对称四极测深模拟, 验证EFG-FE程序的正确性.如 图 5 所示, 第一层电阻率 ρ ₁ =100 Ωm, 厚度 h ₁ =5 m, 第二层电阻率 ρ ₂ =50 Ωm, 厚度 h ₂ =5 m, 第三层电阻率 ρ ₃ =200 Ωm.建立均匀半空间模型水平方向宽度200 m, 垂直方向高100 m的EFGM区域.在该区域内, 地表附近和地层界面附近采用稠密的网格分布, 离场源较远的区域使用稀疏的网格分布, 节点总数15725.当使用FEM计算时, FEM单元数15456, 由于使用矩形单元, 因此采用的节点分布没有采用无单元法计算时分布灵活.当采用无单元法计算时, 为了避免节点不合理分布导致EFGM形函数构造失败, 在EFM计算时的节点分布基础上适当改变了一些节点的分布.在外围区域使用快速延伸的FEM网格扩大计算域, 该区域使用1062个FEM单元, 1190个节点, 使得计算域水平宽度超过4 km, 垂直高度超过2 km.以地面 X =0 m为测点, 模拟不同供电极距对称四极测深, 图 6 给出了两种方法模拟的视电阻率曲线及其相对误差.两种方法的模拟结果与解析解均吻合得很好, 相对误差均在1%以内, EFG-FE结果平均绝对相对误差为0.20%, FEM结果平均绝对相对误差为0.37%, 两种方法均具有很好的模拟精度, 表明了EFG-FE算法的有效性. 4.2 点电源均匀半空间电位模拟 如 图 5a 所示, 建立均匀半空间模型水平方向宽度120 m( X :-60~60 m), 垂直方向高60 m( Z :0~-60 m)的EFGM区域 Ω ₁ , 电阻率为100 Ωm, 点电源设置在地表 X =0 m处, 电流强度1 A.模型节点分布如 图 7a 所示, 为了提高模拟精度同时兼顾计算成本, 利用EFGM的灵活性设置节点以点电源位置为中心向外逐渐由稠密逐渐稀疏分布, 节点主要集中分布在电场强烈变化的场源附近, 共3057个节点.如 图 7b 所示, 在 Ω ₁ 外围建立FEM区域 Ω ₂ , FEM单元设置为向外成双倍扩展延伸8层, 水平范围 X :-2014~2104 m, 垂直范围 Z :0~-2104 m, 共790个节点, 702个单元. 首先, 对 图 7 所示的计算域 Ω ( Ω = Ω ₁ + Ω ₂ )采用第一类边界条件的EFGM(记为F-EFGM)进行点电源电位模拟, 此时在 Ω ₂ 区域采用向外逐渐稀疏的节点离散, 以保证节点合理分布, 使得EFGM形函数构造顺利, 计算域总节点数为7712.然后, 对 图 7a 所示的节点模型和EFGM计算域 Ω ₁ , 采用第三类边界条件的EFGM(记为T-EFGM)( 麻昌英等, 2017 )进行点电源电位模拟.最后, 将 图 7a 和 b 所示的EFGM区域 Ω ₁ 和FEM区域 Ω ₂ 相结合的, 采用第一类边界条件的EFG-FE耦合法进行点电源电位模拟.在1~20 m范围内选取离点电源位置不同距离的35个地面节点的模拟电位进行对比, F-EFGM、T-EFGM和EFG-FE的模拟结果如 图 8 所示.如 图 8 所示, 三种方法的模拟结果与解析解均吻合得很好, 具有很高的模拟精度, 进一步表明了本文的EFG-FE算法的有效性. 表 1 列出了在同一台计算机上三种正演模拟方法分别进行10次模拟统计的平均计算时间和35个采样点的模拟电位平均绝对相对误差, 三种方法均得到了良好的模拟结果, F-EFGM由于采用了更多的节点而模拟精度更好.由 表 1 分析, F-EFGM中在 Ω ₂ 中采用EFGM计算, 大大增加了EFGM计算区域, 使得EFGM计算时间大幅增加, 且为了保证EFGM计算稳定, 大量增加了节点, 增加了方程求解时间. 表 1 表明EFG-FE的计算时间相比于F-EFGM明显降低了, 表明本文提出的EFG-FE可有效解决采用第一类边界条件时由于扩大计算域导致EFGM计算量大量增加的问题.对于单个电极模拟, T-EFGM和EFG-FE计算时间很接近, 下文将展开对多电极供电观测的地电模型进行模拟, 进一步研究EFG-FE的模拟性能. 4.3 复杂地形起伏模型 本小节应用上节所用的T-EFGM和EFG-FE对复杂地形起伏模型进行视电阻率观测模拟.如 图 9 所示, 建立水平方向宽度80 m( X :-70~70 m)、垂直方向宽度40 m( Z :0~-40 m)的复杂地形起伏模型的EFGM区域 Ω ₁ , 在 X =-14 m处为山脊最高处(最高为7.6 m), X =10处为山谷最低处(最高为3.4 m).在直流电阻率EFGM正演模拟中, 可在局部域任意加密节点提高模拟精度, 如在场源附近加密节点可显著提高近场源模拟精度( 麻昌英等, 2017 ).因此, 如 图 9 所示, 针对该复杂地形起伏模型, 在电场场值变化剧烈的场源区域(即浅地表区域)和地形变化区域采用密集的节点分布以提高模拟精度, 同时利用EFGM的强适应性, 采用任意的节点分布以适应任意起伏的地形, 若有进一步加密节点需求, 可直接加密节点而无需其他处理, 相比于FEM中加入新的节点需要重新剖分单元网格, 这体现了EFGM的便利性和灵活性.随着与场源及地形距离的增加, 电场场值变化趋于平缓, 节点分布随之设置为逐渐稀疏, 以降低计算成本.EFGM区域 Ω ₁ 共使用7749个节点.外围采用与上节相同形式的FEM单元网格扩展方式建立EFM区域 Ω ₂ (如 图 7b 所示), 水平范围超过4 km宽, 垂直范围超过2 km高, 共740个节点, 657个单元.复杂地形起伏模型介质电阻率为100 Ωm.在地表 X :-58~58 m范围内等间隔2 m布置59根供电和观测电极, 对模型进行高密度电法温纳装置视电阻率观测正演模拟, 两种方法正演模拟结果及它们所花费的计算时间分别如 图 10 和 表 2 所示. Fig. 10 Pseudo-cross sections of apparent resistivity of a complex undulating terrain simulated model (a) T- EFGM; (b) EFG-FE. 将 图 10a 和 b 对比分析可知, 两种方法的模拟结果几乎一致, 进一步表明当采用相同的EFGM模型时, 本文的EFG-FE可获得与T-EFGM一致的模拟精度. 表 2 列出了两种方法对复杂地形起伏模型的正演计算时间, T-EFGM施加第三类界条件而增加了边界计算时间, 由于边界计算量与电极数、波数个数、边界段数、高斯点数等成正比, 因此在多电极供电观测条件下往往耗费较多计算时间.EFG-FE中由于FEM单元数很少, 且FEM计算效率高, 使得FEM计算时间很少, 但EFG-FE由于增加了FEM区域的节点需要花费更多方程求解时间, 然而EFG-FE中仅使用少量的FEM节点, 因此不会明显增加方程求解时间.根据 表 2 , 虽然EFG-FE需要花费更多的方程求解时间, 但T-EFGM需要较多的边界计算时间, 使得总体上EFG-FE比EFGM花费计算时间少.表明本文的EFG-FE相比于T-EFGM, 通过使用快速延伸的FEM网格扩展计算域, 仅少量增加节点的情况下使得第一类边界条件得到满足, 节省了边界计算时间, 从而提高计算效率.考虑到两者精度基本一致, 表明本文的EFG-FE具有更好的模拟性能.需要指出的是, 外围FEM单元剖分需适当选择, 若选择过于密集的FEM单元剖分, 可能造成FEM节点数量过多而导致EFG-FE方程求解时间明显增加, 反而使得EFG-FE比T-EFGM计算效率低.因此, EFG-FE中应使用合适的的网格剖分方式尽可能发挥FEM的高效性, 从而提高EFG-FE的模拟性能. 本文利用RPIM插值构造无单元法形函数, 在外围区域将EFGM与FEM直接耦合, 提出直流电阻率的无单元Galerkin-有限单元耦合法(EFG-FE).首先, EFG-FE实现了仅需节点离散的EFGM正演模拟, 使得任意复杂电模型可方便地使用任意分布的节点离散覆盖, 而无需网格剖分, 可根据模型和模拟需要任意加密或稀释节点, 利于使节点合理分布, 从而有助于提高模拟精度和控制计算成本.其次, EFG-FE利用FEM的数值稳定性和高效性, 基于矩形单元, 采用快速扩展的单元网格覆盖远离场源的计算域外围区域, 使用少量的FEM节点和单元网格的情况下实现计算域大范围扩大, 使得满足第一类边界条件, 且不造成EFG-FE计算成本大幅增加.最后, 通过多电极供电和观测的复杂模型视电阻率观测正演模拟对比, 结果表明选择适当的外围有限单元法区域及单元剖分, EFG-FE不仅可获得高精度的模拟结果, 且相比于采用第三类边界条件的EFGM计算效率高, 具有更好的模拟性能.虽然EFG-FE相比于传统的数值模拟方法(如FEM)计算效率仍然较低, 但本文的研究可为后续研究计算成本更低的高精度无单元模拟方法提供研究思路, 以及为直流电阻率正演方法的选择提供选项.

特别感谢4位审稿专家中肯的修改意见，使得本文内容更加完善.

Wang Y Y. 2007. Study of element-free galerkin method in the seismic forward modeling. Progress in Geophysics , 22(5): 1539-1544. Xu S Z. 1994. Finite Element Method in Geophysics . Beijing: Beijing Science Press. Zhang Y, Wang J G, Zhang B Y. 2008. Coupled FEM and meshless radial point interpolation method. Journal of Tsinghua University (Science and Technology) , 48(6): 951-954. Zhdanov M S, Cox L H, Endo M. 2015. Large-scale 3D inversion of airborne electromagnetic data based on the hybrid IE-FE method and the moving sensitivity domain approach. //International Workshop and Gravity, Electrical & Magnetic Methods and their Applications. Chenghu, China, 2015: 303-306.