地球物理学报  2015, Vol. 58 Issue (3): 1072-1087   PDF    
基于二次场电导率分块连续变化的三维可控源电磁有限元数值模拟
李勇1,2, 吴小平1,3, 林品荣2    
1. 中国科学技术大学地球和空间科学学院 地震与地球内部物理实验室, 合肥 230026;
2. 中国地质科学院地球物理地球化学勘查研究所, 廊坊 065000;
3. 深部岩土力学与地下工程国家重点实验室, 中国徐州 221008
摘要:从电偶源三维地电断面可控源电磁法的二次电场边值问题及其变分问题出发,采用任意六面体单元对研究区域进行剖分,并且在单元分析中同时对电导率及二次电场进行三线性插值,实现电导率分块连续变化情况下,基于二次场的可控源电磁三维有限元数值模拟.这个新的可控源电磁三维正演方法可以模拟实际勘探中地下任意形状及电性参数连续变化的复杂模型.理论模型的计算结果表明,均匀大地计算的视电阻率误差和相位误差分别为0.002%和0.0005°.分层连续变化模型的有限元计算结果表明,其与对应的分层均匀模型解析结果有明显差异.三维异常体组合模型以及倾斜异常体等复杂模型的有限元计算结果也有效地反映了异常形态.
关键词可控源电磁法     二次场     电导率分块连续变化     有限元法     三维    
Three-dimensional controlled source electromagnetic finite element simulation using the secondary field for continuous variation of electrical conductivity within each block
LI Yong1,2, WU Xiao-Ping1,3, LIN Pin-Rong2    
1. Laboratory of Seismology and Physics of Earth's Interior, School of Earth and Space Sciences, University of Science and Technology of China, Hefei 230026, China;
2. Institute of Geophysical and Geochemical Exploration, Chinese Academy of Geological Science, Langfang 065000, China;
3. State Key Laboratory for GeoMechanics and Deep Underground Engineering, Xuzhou 221008, China
Abstract: We focus on the forward problem of controlled source electromagnetic (CSEM) and present a new three-dimensional (3D) finite-element modeling (FEM) algorithm using the secondary field for continuous variation of electrical conductivity within each block. Usually delta sources are used in the governing partial differential equations and the electrical conductivity of each block is assumed to be a constant in traditional FEM for CSEM modeling, which produces near-field error due to the source singularity and boundary reflections among blocks. On the other hand, the traditional uniform hexahedral grids significantly limit their capacity to handle practical underground structures with complex geometry and continuous variation of electrical conductivity.
On the basis of Maxwell's equations, the electromagnetic field is separated into background field with analytical solution of a layered media and the secondary field caused by anomalous bodies, which is numerically computed by the finite element method, to avoid the source singularity. Firstly, the corresponding variation question of the three-dimension boundary value problem for CSEM is given according to the generalized variational principle. Secondly, the FEM is implemented to solve the variational equation. The computational area is divided into many arbitrary hexahedral element meshes, in which the tri-linear interpolation is performed on the electromagnetic field and conductivity parameter simultaneously to simulate the model with electrical conductivity of continuous variation. Thus the variational equation is converted into a linear equation system, and it is solved to obtain the secondary electric field at each node. In the process, the divergence condition is added to eliminate the spurious solution. Using the relation between the electric fields and the magnetic fields, the secondary magnetic field value of each node is also obtained. In combination with the background field, the total electromagnetic fields of each node can be calculated, resulting in apparent resistivity and impedance phase values on the ground surface. Finally, numerical modeling results for a series of models are carried out to validate the effectiveness of our 3D CSEM finite element algorithm.
The result of our method shows 0.002 percent error in apparent resistivity and 0.0005°error in impedance phase for a homogeneous half-space model. The FEM modeling for a continuous multilayer model shows obviously differences in comparison to analytic solutions for its corresponding layered model, which are up to 44.43% for apparent resistivity and 13.27° for impedance phase. For more complicated models, such as the combined model with multiple 3D anomaly bodies or sloping interface anomaly bodies, our FEM method works well and illustrates subsurface structures effectively.
A CSEM modeling method is developed using FEM to compute the electromagnetic responses for three-dimensional models. The main feature of our method is that the secondary electrical field is solved by FEM with arbitrary hexahedral element meshing, in which the tri-linear interpolation is performed on the electromagnetic field and conductivity parameter simultaneously. In addition, the divergence condition is used to overcome the spurious solution. Finally the total field is obtained by adding the secondary electrical field to the analytical background field. Numerical experiments show our new three-dimensional CSEM forward modeling algorithm ensures the accuracy of the solution and is able to simulate more complicated model with arbitrarily shaped structures and continuous variation of electrical conductivity.
Key words: Controlled source electromagnetic method     Secondary field     Continuous variation of conductivity within each block     Finite element method     Three-dimensional    
1 引言

可控源电磁法(Controlled source electromagnetic method,CSEM)是金属矿、石油天然气、地下水、地热等资源勘查的重要方法(何继善,1990石昆法,1999汤井田,2005底青云等,2008).由于数据解释的需要,二、三维CSEM数值模拟一直是国内外学者研究的热点.Coggon(1971)首次提出2.5维电磁波场问题的有限元理论;Unsworth等(1993)提出了基于二次场CSEM二维有限元数值模拟系统,并重点讨论有限长度电偶源激发场的特征;Mitsuhata(2000)借用地震数值模拟中震源的加载方法,使用伪delta函数的方法模拟了2.5维可控源电磁响应;Li和Key(2007)提出基于后验误差估计的自适应有限元算法并进行二维海洋可控源电磁法模拟;孟永良和罗延钟(1996)实现了基于总场的二维CSEM正演;陈小斌和胡文宝(2002)采用有限元直接迭代算法实现了线源频率域响应的二维正演计算;王若等(2006)张继锋等(2009a)实现了频率域线源大地电磁法有限元正演模拟;底青云(2004a2004b)建立了2.5维有限元CSAMT数值模拟方法,较准确地获得复杂介质结构的波场特征;雷达(2010)针对起伏地形下CSAMT二维地电模型,采用加权余弦数值积分法,进行了波数域电磁场二维有限单元法正演.三维CSEM正演模拟多采用有限差分法(沈金松,2003Weiss and Constable, 2006付长民等,2009邓居智等,2011)和积分方程法(Zhdanov, et al.,2006),三维CSEM有限元数值模拟的工作目前还很少.Pridmore等(1981)采用有限单元法求解三维地电结构的电磁响应;闫述和陈明生(2000)王若等(2007)发表了有关三维频率域电磁测深地电模型有限元正演;汤井田和何继善(2007)系统建立了Coulomb规范下磁矢量势-电标量势的hp型自适应有限元分析的基本理论,但没有给出计算结果;张继锋(2009b)薛云锋和张继锋(2011)在Mitsuhata加载源方法的基础上,实现了基于总场的三维可控源电磁法有限元数值模拟;徐志锋和吴小平(2010)采 用电磁场的磁矢量位和电标量势,实现了可控源电磁三维频率域有限元模拟方法,但没有给出磁场的计算结果,因而不能计算常用的卡尼亚视电阻率.

上述二、三维CSEM数值模拟方法均未考虑地下介质电性连续变化的情况,而实际地下岩矿石的组分及其所在环境的温度、湿度可能是渐变的,使得电导率参数往往连续变化.针对实际勘探中地下介质的这种连续变化,阮百尧和徐世浙(1998)阮百尧和熊彬(2000)实现了电阻率分块连续变化的二、三维直流电阻率法有限元数值模拟;李勇等(2011)实现了基于异常复电位复电导率分块连续变化的2.5维复电阻率有限元数值模拟;徐世浙等(1995)刘云和王绪本(2012)实现了电性参数分块连续变化的二维大地电磁(MT)有限元数值模拟.可控源电磁法方面,仅Li和Key(2007)在二维海洋可控源电磁法有限元模拟中考虑了电性参数分块连续变化,还未见三维情况下的相关工作成果发表.
另外,CSEM采用人工可控信号源,基于总场CSEM有限元数值模拟由于源的奇异性会带来较大的计算误差,因为其变分问题包含源项(张继锋,2009b),总场在源处是奇异的,线形插值函数不能精确描述源附近总场的快速衰减,从而带来较大的计算误差,通过对源附近的网格加密也只能有限地提高计算精度,但计算量将大量增加.基于二次场CSEM有限元数值模拟方法,源产生的一次场单独处理,其有解析解,因此在变分问题中不包含源项,计算精度会提高,且不增加有限元方法的计算工作量.本文研究新的基于二次场CSEM有限元数值模拟方法,采用任意六面体单元对研究区域进行剖分,并且在单元分析中同时对电导率及二次电场进行三线性插值,以模拟实际勘探中地下任意形状及电性参数连续变化的复杂模型.

2 二次电场变分问题

考虑如图 1a所示的三维地电模型,全空间由空气、异常体和大地组成,电导率σ是空间任意点的函数,任意方向的水平电偶极子源置于地表的任意位置. 假定大地中的磁导率和介电常数与自由空间磁导率μ和空气介电常数ε相同,取时谐因子e-iωt.根据电磁场理论,三维地电模型中任意点的电场强度E和磁场强度H满足如下的麦克斯韦方程式:

式中,ω是角频率,i是虚数单位,J是激励源,即水平电偶极子源Tx电流密度.

对(1)式两边取旋度,并将(2)式代入可导出电场强度E所满足的微分方程为

式中,k为波数,

根据叠加原理,总电场E分解为一次电场E0和二次电场Es,有

仿照(3)式推导,E0和Es分别满足下述方程:

式中,k202με+iωμσ0,k2s=iωμσs,k22με+iωμσ,σ0为均匀介质(背景)电导率,σs是异常电导率,σs=σ-σ0.

图 1 三维模型及离散化示意图
(a)三维模型;(b)三维模型离散化.
Fig. 1 Schematic diagram for a three-dimensional model and the model meshing
(a)Three-dimensional model;(b)The 3-D model meshing.

方程(5)是均匀半空间中水平电偶极子激发产生的一次电场,可用解析法求解.方程(6)是二次电场Es所满足的矢量波动方程,取外边界Γ足够大,使异常体产生的二次电场Es衰减到外边界处趋于0,可获得足够精确的第一类边界条件,构成三维地电模型CSEM二次电场的边值问题.

根据广义变分原理(徐世浙,1994),上述二次电场边值问题对应的变分问题为

式中,V为三维研究区域,F(Es)是Es的泛函,δF(Es)为F(Es)的变分.

3 有限单元法

求解CSEM三维地电模型二次场变分问题的有限元方法步骤如下:

3.1 网格剖分

首先对研究区域进行三维网格剖分.为了模拟不规则形状的复杂地电结构,在x、y方向,采用任意形状的六面体单元将区域划分成若干个曲面层,每个曲面层包括目标区域和网格外延区域(网格外 延区域主要是为了模拟无穷边界),目标区域x、y方向采用等间隔网格,网格外延区域采用不等距网格;z方向(包括空气层)采用不等距网格.目标区域为地质体赋存区域,也是数据采集区域,地表的源可以在目标区也可以在网格外延区域,如图 1b所示.为了便于计算,对剖分网格中的任意六面体单元进行等参变换,将直角坐标系内不规则的网格单元(如图 2a所示的子单元)转换成自然坐标系下的正方形网格(如图 2b所示的母单元),这样,式(7)在整个研 究区域V的积分可分解为各六面体等参单元的积分:

而子单元中的体积元dxdydz与母单元中的体积元dζdηdγ的关系为

图 2 任意形状的六面体单元示意图
(a)子单元;(b)母单元.
Fig. 2 Schematic diagram for an arbitrary hexahedral element
(a)Sub-element;(b)Parent element.
3.2 线性插值

将任意六面体的8个角点取为节点,其节点的编号及坐标如图 2所示.在六面体单元内一次电场Ee0、二次电场Ees、电导率σ和异常电导率σs均呈三线性变化,则在每个单元内Ee0,Ees,σ和σs满足如下的表达式:

式中,i=1,2,…,8;Ee0i为单元节点i上的一次电场值;Eesi为单元节点i上的待定二次电场值,σi和σsi为单元节点i上的电导率和异常电导率,Ni为形函数,且有

其中,ζi,ηi,γi是自然坐标系ζ,η,γ下节点i的坐标.

3.3 单元分析

式(8)的第1项单元积分:

其中,为单元内x、y、z方向的二次电场.
将式(12)写成矩阵形式,有

其中

N=[N1 N2 N3 N4 N5 N6 N7 N8]T
式(8)的第2项单元积分

将式(14)写成矩阵形式,有

其中,

式(8)的第3项单元积分

将式(16)写成矩阵形式,有

其中,

3.4 总体合成及求变分

在单元e内,将(13)、(15)、(17)式的积分结果相加,再扩展成全部单元的相加,有

以上步骤将连续函数Es的泛函F(Es)离散成各节点Es值的多元函数,对多元函数求极值即可获得泛函的极值,有

将(18)式代入(19)式,得线性方程组


代入边界条件,解线性方程组(20)式,得各节点上的二次电场Es,再加上一次电场E0之后获得各节点上的总电场E.

求解过程中,常规立方体剖分时单元积分得到矩阵中的各个元素为常数(徐世浙,1994).而采用任 意形状的六面体剖分,其单元积分需要用等参变 换来处理,把式(9)代入式(20)将形成形如类型的积分,被积函数f(ζ,η,γ)是与六面体曲面形态有关的复杂函数,需采用高斯数值积分计算(徐世浙,1994).

4 散度条件的施加

计算结果表明,采用上述节点有限元离散方程(6)式来计算二次电场,计算结果存在伪解,有资料表明,这种解不满足散度条件(金建铭,1998底青云和王若,2008).为了消除伪解,在(7)式的变分问题中增加一个罚项来强加散度条件,有

则有

其中,

则强加散度条件后,线性方程组式(20)式转换为

5 磁场、视电阻率和阻抗相位的计算

将式(1)展开,二次磁场Hsy的表达式为

其中,

式中,Esx1为地面上的节点x方向二次电场值;Esx2,Esx3,Esx4分别为地面下第一、二、三个节点x方向二次电场值,z为z方向的网格间距;Esz1为地面上的所求偏导数对应的节点的垂直分量的二次电场值;Esz2,Esz3,Esz4分别为在x方向与之相邻的第一、二、三个节点的垂向二次电场值,x为x方向的网格间距.

求得 Esx z和Esz x 后代入式(25)便可计算出地面上节点的二次磁场值Hsy,加上由解析解计算的一次磁场值H0y获得各节点上的总磁场Hy,便可用下式分别计算视电阻率ρxy和相位φxy

6 模型计算与结果分析

为了验证算法的正确性,以及对电导率连续变化的模型进行有限元计算,本文设计了五个理论地电模型.其所有模型的三维有限元网格剖分均如图 1b所示,x、y方向目标区(-3000 m≤x≤3000 m、-3000 m≤y≤3000 m)采用均匀网格剖分,网格间距100 m,网格节点数均为61个,网格外延区域采用稀疏网格,其左、右、前、后网格节点数均为13个,网格间距逐渐向边界扩展,分别为200、400、400 m、 800、800、1600、1600、3200、3200、6400、12800、25600 m; z方向采用非均匀网格剖分,网格节点由空气层的网格节点和大地的网格节点组成,空气层的网格节点 数是18个,其网格节点坐标从上到下分别为-50 km、 -30、-10、-8、-5、-2、-1.5 km和-1 km; -500、-300、-200、-175、-150、-125、-100 m、-75、-50 m和-25 m,大地的网格节点数61个,网格节点坐标从0~-50 km.值得注意的是,为模拟电导率分块连续变化,本文与传统分块均匀模型将各单元的电导率赋予一个常数的剖分不同,电导率值是赋予在剖分单元的节点上,通过三线性插值的方法实现电导率分块连续变化.计算频率为f=2nHz,n取值从-6到12,共19个频点,如表 1中所示.

表 1 均匀半空间模型三维有限元数值模拟解与解析解对比 Table 1 Comparison between analytical and 3-D FE solutions for homogeneous half-space model
6.1 均匀半空间模型三维有限元数值模拟精度检验

均匀半空间模型如图 3所示,电阻率为100 Ωm,发射源长度和电流强度分别为1 m和1A的电偶极子Tx位于目标区域,Tx中点与直角坐标系坐标原点0重合,其方向沿着x方向.该均匀半空间模型的视电阻率和相位存在解析解.点A(0,1500,0)m处有限元三维数值模拟计算的不同频率视电阻率和相 位与解析结果的比较见表 1,相应的视电阻率和相 位曲线如图 4所示,解析结果采用中国地质科学院地球物理地球化学勘查研究所研发的的电法勘探工作站WEM2.5计算.用公式

表示有限元法计算的视电阻率ρxy与解析解比较的相对误差,其中ρ′xy和ρxy分别是解析法和有限元法计算的视电阻率,n是频点数.用公式

表示相位的计算误差,其中φ′ xy和φxy分别是解析法和有限单元法计算的相位.

图 3 均匀半空间模型示意图 Fig. 3 chematic diagram of homogeneous half-space model

图 4 均匀半空间模型三维有限单元计算结果与解析解结果对比
(a)视电阻率;(b)相位.
Fig. 4 Comparison of the analytical and 3-D FEM solutions for homogeneous half-space model
(a)Apparent resistivity;(b)Phase.

由公式(28)、(29)计算得到视电阻率相对误差Xρxy和相位误差Xφxy分别为0.002%和0.0005°,三维有限元计算结果与解析解在很宽的频段上均吻合非常好.图 5还给出频率2 Hz的电场Ex和磁场Hy振幅全空间分布,可以直观地看到,地表的 Ex和Hy 呈4个花瓣状的图案,4个花瓣被4个低值带所分开,为可能的三维全空间测量提供参考.

图 5 模型一电磁场空间分布示意图
(a)电场振幅;(b)磁场振幅.
Fig. 5 Space distribution of the electrical and magnetic fields for homogeneous half-space model
(a)The amplitude of the electrical field;(b)The amplitude of the magnetic field.
6.2 水平层状连续介质模型

水平层状一维连续介质模型如图 6所示.发射源长度和电流强度分别为1000 m和30 A的电偶 极子Tx位于网格外延区域,Tx中点坐标为(0、-8000、 0)m,其方向平行于x轴.第一层电阻率为10 Ωm,厚度为500 m,第三层电阻率为500 Ωm,第二层为电阻率连续变化介质,厚度为500 m,电阻率随着深度 线性增加,其电阻率变化如图 7a所示,该模型目前 未见到解析解,可由本文的有限元方法计算.图 7b是对应的三层电阻率均匀模型,有解析解,其第一层和第三层电阻率与图 7a相同,中间层厚度仍为500 m,但电阻率为第一和第三层电阻率的对数平均值,以近似图 7a的线性变化.分层连续变化模型在原点的视电阻率和相位与三层均匀模型的解析结果比较见表 2,相应的视电阻率和相位曲线如图 8所示.可以看出,分层连续变化介质模型与均匀模型的视电阻率曲线和相位曲线形态基本一致,但在某些频点上差异较大,视电阻率差异最大可达44.43%,相位差可达13.27°,而且在低频的视电阻曲线上产生了 系统偏差,表明电阻率连续变化相对电阻率均匀层 状模型有明显影响.如果实际地下岩石电阻率为连续变化,而依然采用传统的分层均匀模型进行解释,显然会造成较大的误差,值得引起重视.

图 6 水平层状一维连续介质模型示意图
Fig. 6 Schematic diagram of a stratified model with continuous media

图 7 水平层状一维介质模型电阻率分布
(a)介质分层连续变化模型;(b)三层均匀介质模型.
Fig. 7 Distribution of resistivity for 1D stratified models
(a)Continuous variation of resistivity;(b)Step variation of resistivity.

图 8 水平层状一维介质模型三维有限单元计算结果与解析解结果对比
(a)视电阻率;(b)相位.
Fig. 8 Comparison between the analytical and 3-D FEM solutions for 1D stratified models
(a)Apparent resistivity;(b)Phase.

表 2 水平层状一维介质模型三维有限元数值模拟解与解析解对比 Table 2 Comparison between the analytical and 3-D FEM solutions for 1D stratified models
6.2 两个三维异常体模型

为了模拟电导率连续变化三维地电模型地表二次电磁场的分布,设计了两个异常体模型,如图 9所示,电阻率为1000 Ωm的均匀半空间中存在两个电 阻率为100 Ωm的低阻异常体(异常体大小500 m×500 m×1000 m,埋深500 m).用x-y平面图和x-z断面图表示的网格离散化后空间电阻率的分布示意图如图 9(a,b).注意到,异常体与围岩过渡区域的电导率是连续变化的,如图 9(a,b)的色标渐变所示.发射源Tx(长度和电流强度分别为1500 m和 20 A)位于网格外延区域,其方向与x轴平行,Tx中点坐标为(0、-8000、0)m.该模型不同频率地表二次电磁场、视电阻率和阻抗相位数值模拟计算结果如图 10所示,从图中可以看出,地表二次电磁场等值线图较好地显示了两个异常体的位置,磁场的异常分辨率大于电场的异常分辨率,另外,从视电阻率和相位等值图可以看出,频率为128 Hz时异常反应比较明显,这主要是与异常体的埋藏深度有密切的关系,随着频率的降低,视电阻率和相位异常分辨率随之降低.

图 9 与围岩过渡区域电导率连续变化的两个三维异常体模型示意图
(a)平面图;(b)断面图.
Fig. 9 Schematic diagram for the model with two 3D anomaly bodies around which shows continuous variation of electrical conductivity

图 10 两个三维异常体模型不同频率地表二次电磁场、视电阻率和相位等值线图
(a1)、(a2)、(a3)和(a4)分别为1024Hz地表二次电场振幅、二次磁场振幅、视电阻率和相位;
(b1)、(b2)、(b3)和(b4)分别为128Hz地表二次电场振幅、二次磁场振幅、视电阻率和相位;
(c1)、(c2)、(c3)和(c4)分别为16Hz地表二次电场振幅、二次磁场振幅、视电阻率和相位.
Fig. 10 Contour plot of the secondary fields,apparent resistivities and phases with different frequencies on the ground for the model with two 3D anomaly bodies
Fig.(a1),(a2),(a3) and (a4)show the amplitudes of the secondary electric and magnetic fields,the apparent resistivity and phase at 1000Hz on the surface,respectively; Fig.(b1),(b2),(b3) and (b4)show the amplitudes of the secondary electric and magnetic fields,the apparent resistivity and phase at 128Hz on the surface,respectively; Fig.(c1),(c2),(c3) and (c4)show the amplitudes of the secondary electric and magnetic fields,the apparent resistivity and phase at 16Hz on the surface,respectively.

6.3 “V”字型异常体三维模型

为了模拟电导率连续变化倾斜界面的异常特征,设计“V”字型异常体,如图 11所示,电阻率为1000 Ωm的均匀半空间中含有电阻率为100 Ωm、厚度为100 m的“V”字型异常体,“V”字型异常体在地表投影的中心剖面为y=0 m,其电阻率分布断面图如图 11a所示,“V”字型异常体沿着中心剖面向两侧分别沿伸400 m,网格离散后电阻率的空间分布如图 11b所示,与图 9类似异常体周边的电阻率也是渐变的,发射源Tx(长度和电流强度分别为1000 m和30 A)位于网格外延区域,其方向与x轴平行,Tx中点的坐标为(0,-15000,0)m.

图 11图 12分别是地表不同剖面y=-600 m、 -300、0、300 m和600 m的视电阻率等值线图和相位等值线图,纵坐标为频率,横坐标x为测点坐标,测点间距100 m,x=-3000~3000 m.由图 12图 13可以看出,视电阻率和相位等值线很好地反应了“V”字型异常形态.

图 11 与周边电性渐变的V字型异常体模型示意图
(a)断面图;(b)电阻率分布图.
Fig. 11 Schematic diagram for anomaly body shaped as V and continuous variation of electrical conductivity is shown around it.
(a)Section view;(b)The spatial distribution of resistivity.

图 12 不同剖面视电阻率等值线图 Fig. 12 Contour plot of apparent resistivity in x planes

图 13 不同剖面相位等值线图 Fig. 13 Contour plot of phase in x planes
6.4 电导率连续变化的复杂模型

为了模拟电导率连续变化的复杂地质情况异常特征,图 14为电阻率100 Ωm均匀半空间中有电阻率分别为20、500、20、500、20 Ωm 5个相间的异常体,在x、y、z方向的长度均分别为400、800、400 m,埋深均为200 m,其周边介质的电阻率同样显示为 渐变.图 14a表示地下异常体在地表的投影,图 14b表 示 y=0断面(x=-3000~3000 m,z=-2000~0 m)的电阻率分布,发射源Tx(长度和电流强度分别为1000 m和30 A)位于网格外延区域,其方向与x轴平行,中心点的坐标为(0、-12000、0)m.图 15是电偶极子发射频率为8 Hz时地表y=0剖面(x=-3000~3000 m,z=0)的视电阻率曲线和相位曲线,图 16是电偶极子发射频率为8 Hz时地表的视电阻率等值线和相位等值线,由图可以看出,数值模拟结果反映了异常体高低阻相间的特性,与地下异常体的分布相吻合.

图 14 电导率连续变化的复杂模型示意图
(a)平面图;(b)断面图.
Fig. 14 Schematic diagram for combined complex model with multiple 3D anomaly bodies and continuous variation of electrical conductivity around them
(a)Plan view;(b)Section view.

图 15 8 Hz时正演模拟结果
(a)视电阻率曲线图;(b)相位曲线图.
Fig. 15 The apparent resistivity and phase curve for combined model with multiple 3D anomaly bodies
(a)Apparent resistivity;(b)Phase.

图 16 率为8 Hz的地表视电阻率和相位等值线图
(a)视电阻率(单位:Ωm);(b)相位(单位:(°)).
Fig. 16 Contour plots of the apparent resistivity and phase on the ground at 8 Hz
(a)Apparent resistivity;(b)Phase.
7 结论

本文采用任意六面体单元剖分,并且在单元分析中同时对电导率及二次电场进行三线性插值,实现了基于二次场电导率分块连续变化的三维可控源电磁有限元数值模拟新方法.同时,在变分问题中增加一个罚项来强加散度条件,克服了有限单元求解二次电场伪解的出现,均匀大地模型理论计算结果表明,本文的三维可控源电磁有限元数值模拟取得了很高的计算精度.进一步在模拟实际勘探中地下任意形状及多异常体组合、电性参数连续变化等复杂模型的计算中,也均取得可靠的结果.

分层连续变化模型的有限元计算结果表明,其与对应的层状模型解析结果有明显差异.实际数据解释中,应避免用简单的分层模型代替可能的地下介质连续变化模型,否则会造成较大的偏差.

致谢 感谢两位匿名评审专家提出了非常有价值的修改意见,对论文总体质量的提高帮助很大.中国地质科学院地球物理地球化学勘查研究所王书民博士审阅了修改稿,在此表示致谢!

参考文献
[1] Chen X B, Hu W B. 2002. Direct iterative finite element (dife) algorithm and its application to electromagnetic response modeling of line current source. Chinese J. Geophys. (in Chinese), 45(1): 119-130.
[2] Coggon J H. 1971. Electromagnetic and electrical modeling by the finite element method. Geophysics, 36(1): 132-155.
[3] Deng J Z, Tan H D, Chen H, et al. 2011. CSAMT 3D modeling using staggered-grid finite difference method. Progress in Geophysics (in Chinese), 26(6): 2026-2032, doi: 10.3969/j.issn. 1004-2903.2011.06.017.
[4] Di Q Y, Unsworth M, Wang M Y. 2004a. 2. 5D CSAMT modeling with finite element method. Progress in Geophysics (in Chinese), 19(2): 317-324,10.3969/j.issn. 1004-2903.2004.02.018.
[5] Di Q Y, Unsworth M, Wang M Y. 2004b. 2. 5D CSAMT modeling with finite element over 2D complex earth media. Chinese J. Geophys. (in Chinese), 47(4): 723-730.
[6] Di Q Y, Wang R. 2008. CSAMT Forward Modeling and Inversion and Its Application (in Chinese). Beijing: Science Press.
[7] Fu C M, Di Q Y, Wang M Y. 2009. 3D numeric simulation of marine controlled source electromagnetics (MCSEM). Oil Geophysical Prospecting (in Chinese), 44(3): 358-363.
[8] He J S. 1990. Controlled Source Audio Frequency Magnetotellurics Method (in Chinese). Changsha: Central South University of Technology Press.
[9] Jin J M. 1998. Finite Element Method of Electromagnetic Field (in Chinese). Xi'an: Xi'an Electronic Science and Technology University Press.
[10] Lei D. 2010. Studies and applications of 2-D CSAMT modeling and inversion with a dipole source and topography. Chinese J. Geophys. (in Chinese), 53(4): 982-993, doi: 10.3969/j.issn.0001-5733.2010.04.023.
[11] Li Y, Lin P R, Li T L, et al. 2011. Finite element method for solving anomalous complex potential of 2. 5-D complex resistivity. Journal of Jilin University (Earth Science Edition) (in Chinese), 41(5): 1596-1604.
[12] Li Y G, Key K. 2007. 2D marine controlled-source electromagnetic modeling: Part 1—An adaptive finite-element algorithm. Geophysics, 72(2): WA51-WA62.
[13] Liu Y, Wang X B. 2012. The FEM for modeling 2-D MT with continuous variation of electric parameters within each block. Chinese J. Geophys. (in Chinese), 55(6): 2079-2086, doi: 10.6038/j.issn..0001-5733.2012.06.029.
[14] Meng Y L, Luo Y Z. 1996. Two-dimensional bipolar source CSAMT finite element forward (in Chinese). //Collected Works Geophysical Exploration Geochemistry. No. 20 Electrical Album. Beijing: Geological Publishing House.
[15] Mitsuhata Y J. 2000. 2-D electromagnetic modeling by finite-element method with a dipole source and topography. Geophysics, 65(2): 465-475.
[16] Pridmore D F, Hohmann G W, Ward S H, et al. 1981. An investigation of finite-element modeling for electrical and electromagnetic data in three dimensions. Geophysics, 46(7): 1009-1024.
[17] Ruan B Y, Xu S Z. 1998. FEM for modeling resistivity sounding on 2-D geoelectric model with line variation of conductivity within each block. Earth Science-Journal of China University of Geosciences (in Chinese), 23(3): 303-307.
[18] Ruan B Y, Xiong B. 2000. A finite element modeling of three-dimensional resistivity sounding with continuous conductivity. Chinese J. Geophys. (in Chinese), 45(1): 131-138.
[19] Shen J S. 2003. Modeling of 3-D electromagnetic responses in frequency domain by using straggered-grid finite difference method. Chinese J. Geophys. (in Chinese), 46(2): 281-288.
[20] Shi K F. 1999. The Theory and Application of Controlled-Source Magnetotellurics (in Chinese). Beijing: Science Press.
[21] Tang J T, He J S. 2005. Controlled Source Audio Frequency Magnetotellurics Method and Its Applications (in Chinese). Changsha: Central South University Press.
[22] Tang J T, Ren Z Y, Hua X R. 2007. Theoretical analysis of geoelectromagnetic modeling on Coulomb gauged potentials by adaptive finite element method. Chinese J. Geophys. (in Chinese), 50(5): 1584-1594.
[23] Unsworth M J, Travis B J, Chave A D. 1993. Electromagnetic induction by a finite electric dipole source over a 2-D earth. Geophysics, 58(2): 198-214.
[24] Wang R, Wang M Y, Di Q Y. 2006. Electromagnetic modeling due to line source in frequency domain using finite element method. Chinese J. Geophys. (in Chinese), 49(6): 1858-1866.
[25] Wang R, Wang M Y, Lu Y L. 2007. Preliminary study on 3D CSAMT method modeling using finite element method. Progress in Geophysics (in Chinese), 22(2): 579-585, doi: 10.3969/j.issn.1004-2903.2007.02.035.
[26] Weiss C J, Constable S. 2006. Mapping thin resistors and hydrocarbons with marine EM methods, Part II: Modeling and analysis in 3D. Geophysics, 71(6): G321-G332.
[27] Xu S Z. 1994. The Finte Element Method in Geophysics (in Chinese). Beijing: Science Press.
[28] Xu S Z, Yu T, Li Y G, et al. 1995. The finite element method for modeling 2-D MT field on a geoelectrical model with continuous variation of conductivity within each block(Ⅰ). Geological Journal of Universities (in Chinese), 1(2): 65-73.
[29] Xu Z F, Wu X P. 2010. Controlled source electromagnetic 3-D modeling in frequency domain by finite element method. Chinese J. Geophys. (in Chinese), 53(8): 1931-1939, doi: 10.3969/j.issn.0001-5733.2010.08.019.
[30] Xue Y F, Zhang J F. 2011. Three dimensional controlled source electromagnetic numerical simulation based on the rock properties of the west line of South-to-North Water Diversion Project using Finite Element Method. Chinese J. Geophys. (in Chinese), 54(8): 2160-2168, doi: 10.3969/j.issn.0001-5733.2011.08.024.
[31] Yan S, Chen M S. 2000. Finite element solution of three dimensional geoelectric models in frequency electromagnetic sounding excited by a horizontal electric dipole. Coal Geology & Exploration (in Chinese), 28(3): 50-56.
[32] Zhang J F, Tang J T, Yu Y, et al. 2009a. Finite element numerical simulation on line controlled source based on quadratic interpolation. Journal of Jilin University (Earth Science Edition) (in Chinese), 39(5): 929-935.
[33] Zhang J F, Tang J T, Yu Y, et al. 2009b. Three dimensional controlled source electromagnetic numerical simulation based on electric field vector wave equation using finite element method. Chinese J. Geophys. (in Chinese), 52(12): 3132-3141, doi: 10.3969/j.issn.0001-5733.2009.12.023.
[34] Zhdanov M S, Lee S K, Yoshioka K. 2006. Integral equation method for 3D modeling of electromagnetic fields in complex structures with inhomogeneous background conductivity. Geophysics, 71(6): G333-G345.
[1] 陈小斌, 胡文宝. 2002. 有限元直接迭代算法及其在线源频率域电磁响应计算中的应用. 地球物理学报, 45(1): 119-130.
[2] 邓居智, 谭捍东, 陈辉等. 2011. CSAMT三维交错采样有限差分数值模拟. 地球物理学进展, 26(6): 2026-2032, doi: 10.3969/j.issn.1004-2903.2011.06.017.
[3] 底青云, Unsworth M, 王妙月. 2004a. 有限元法2. 5维CSAMT数值模拟. 地球物理学进展, 19(2): 317-324, doi: 10.3969/j.issn.1004-2903.2004.02.018.
[4] 底青云, Unsworth M, 王妙月. 2004b. 复杂介质有限元法2. 5维可控源音频大地电磁法数值模拟. 地球物理学报, 47(4): 723-730.
[5] 底青云, 王若. 2008. 可控源音频大地电磁数据正反演及方法应用. 北京: 科学出版社.
[6] 付长民, 底青云, 王妙月. 2009. 海洋可控源电磁法三维数值模拟. 石油地球物理勘探, 44(3): 358-363.
[7] 何继善. 1990. 可控源音频大地电磁法. 长沙: 中南工业大学出版社.
[8] 金建铭. 1998. 电磁场有限单元法. 西安: 西安电子科技大学出版社.
[9] 雷达. 2010. 起伏地形下CSAMT二维正反演研究与应用. 地球物理学报, 53(4): 982-993, doi: 10.3969/j.issn.0001-5733.2010.04.023.
[10] 李勇, 林品荣, 李桐林等. 2011. 基于异常复电位2.5维CR有限元数值模拟. 吉林大学学报(地球科学版), 41(5): 1596-1604.
[11] 刘云, 王绪本. 2012. 电性参数分块连续变化二维MT有限元数值模拟. 地球物理学报, 55(6): 2079-2086, doi: 10.6038/j.issn.0001-5733.2012.06.029.
[12] 孟永良, 罗延钟. 1996. 双极源CSAMT法二维正演的有限单元算法. //勘查地球物理勘查地球化学文集. 第20集电法专辑. 北京: 地质出版社.
[13] 阮百尧, 徐世浙. 1998. 电导率分块线性变化二维地电断面电阻率测深有限元数值模拟. 地球科学—中国地质大学学报, 23(3): 303-307.
[14] 阮百尧, 熊彬. 2000. 电导率连续变化的三维电阻率测深有限元模[CM)][LL] 拟. 地球物理学报, 45(1): 131-138.
[15] 沈金松. 2003. 用交错网格有限差分法计算三维频率域电磁响应. 地球物理学报, 46(2): 281-288.
[16] 石昆法. 1999. 可控源音频大地电磁法理论与应用. 北京: 科学出版社.
[17] 汤井田, 何继善. 2005. 可控源音频大地电磁法及其应用. 长沙: 中南大学出版社.
[18] 汤井田, 任政勇, 化希瑞. 2007. Coulomb规范下地电磁场的自适应有限元模拟的理论分析. 地球物理学报, 50(5): 1584-1594.
[19] 王若, 王妙月, 底青云. 2006. 频率域线源大地电磁法有限元正演模拟. 地球物理学报, 49(6): 1858-1866.
[20] 王若, 王妙月, 卢元林. 2007. 三维三分量CSAMT法有限元正演模拟研究初探. 地球物理学进展, 22(2): 579-585, doi: 10.3969/j.issn.1004-2903.2007.02.035.
[21] 徐世浙. 1994. 地球物理中的有限单元法. 北京: 科学出版社.
[22] 徐世浙, 于涛, 李予国等. 1995. 电导率分块连续变化的二维MT有限元模拟(Ⅰ). 高校地质学报, 1(2): 65-73.
[23] 徐志锋, 吴小平. 2010. 可控源电磁三维频率域有限元模拟. 地球物理学报, 53(8): 1931-1939, doi: 10.3969/j.issn.0001-5733.2010.08.019.
[24] 薛云锋, 张继锋. 2011. 基于南水北调西线工程岩性特征的CSAMT法有限元三维数值模拟研究. 地球物理学报, 54(8): 2160-2168, doi: 10.3969/j.issn.0001-5733.2011.08.024.
[25] 闫述, 陈明生. 2000. 电偶源频率电磁测深三维地电模型有限元正演. 煤田地质与勘探, 28(3): 50-56.
[26] 张继锋, 汤井田, 喻言等. 2009a. 基于二次插值的线源可控源有限元数值模拟. 吉林大学学报(地球科学版), 39(5): 929-935.
[27] 张继锋, 汤井田, 喻言等. 2009b. 基于电场矢量波动方程的三维可控源电磁法有限单元法数值模拟. 地球物理学报, 52(12): 3132-3141, doi: 10.3969/j.issn.0001-5733.2009.12.023