地球物理学报  2018, Vol. 61 Issue (4): 1549-1562   PDF    
复杂地下异常体的可控源电磁法积分方程正演
汤井田1,2,3, 周峰1, 任政勇1,2,3 , 肖晓1,2,3, 邱乐稳1, 陈超健1, 陈煌1     
1. 中南大学地球科学与信息物理学院, 长沙 410083;
2. 中南大学有色金属成矿预测与地质环境监测教育部重点实验室, 长沙 410083;
3. 有色资源与地质灾害探查湖南省重点实验室, 长沙 410083
摘要:可控源电磁法具有分辨率高及抗干扰能力强等特点,是一种重要的地电磁勘探方法.目前,可控源电磁法的高精度正演计算一直是其核心研究问题之一.传统积分方程法一般采用近似积分公式、简单矩形网格和近似的奇异性体积分计算技术,制约了体积分方程法处理复杂地下异常体的能力,降低了计算精度.针对上述问题,本文基于完全积分公式、四面体非结构化网格和奇异体积分的精确解析解来高精度求解复杂可控源电磁模型的正演响应.首先,从电场积分公式出发,推导了可控源电磁问题满足的积分方程;其次,借助于非结构化四面体网格离散技术,实现了地下复杂异常体的有效模拟.最后,利用散度定理把强奇异值体积分转换为一系列弱奇异性的面积分公式,并通过推导获得了这些弱奇异性的面积分公式的解析解,从而最终实现三维可控源电磁问题的高精度积分求解.以块状低阻体地电模型为测试模型,采用本文提出的积分方程方法获得的数值解与其他公开数值算法解进行对比分析,其对比结果具有高度的吻合性,验证了算法的正确性;同时,设计了球状及复杂地电模型进行算法收敛性测试,进一步验证算法的正确性以及能够处理地下复杂模型的能力.
关键词: 可控源电磁法      积分方程      非结构化      正演      奇异值积分     
Three-dimensional forward modeling of the controlled-source electromagnetic problem based on the integral equation method with an unstructured grid
TANG JingTian1,2,3, ZHOU Feng1, REN ZhengYong1,2,3, XIAO Xiao1,2,3, QIU LeWen1, CHEN ChaoJian1, CHEN Huang1     
1. School of Geosciences and Info-Physics of Central South University, Changsha 410083, China;
2. Key Laboratory of Metallogenic Prediction of Nonferrous Metals and Geological Environment Monitoring of Ministry of Education, Central South University, Changsha 410083, China;
3. Key Laboratory of Non-Ferrous Resources and Geological Hazard Detection, Changsha 410083, China
Abstract: The controlled source electromagnetic method (CSEM) is characterized by high resolution and strong anti-interference ability, which is an important tool in geo-electromagnetic exploration. Inversion is a key step of data processing and interpretation in this method, while the forward modeling is the foundation of inversion. Therefore, searching for a high-accuracy forward algorithm is one of the core research questions to interpret CSEM data. The traditional volume integral equation formula is successfully applied to compute the electromagnetic response of the 3D CSEM as a semi-analytical solution. This method often adopts the approximate integral formula, a regular hexahedron grid and the approximate singular value integral processing technique, which restricts the ability of the volume integral equation method to deal with anomaly bodies with arbitrary complex geometry and reduces its calculation precision. To solve these problems, a new integral strategy is proposed to accurately calculate the 3D controlled-source electromagnetic forward response based on the complete integration formula using a tetrahedral unstructured grid and singularity-free analytical solution for the singular volume integral. Firstly, the integral equation of the CSEM problem is deduced from the formula of the electric field integral. Then, the underground complex abnormal body is discretized by the latest unstructured discrete technique based on a tetrahedral grid. Using the divergence theorem, we transform the strong singular value volume integral into a series of weak singularity integral formulas. And we obtain the analytic solutions of these weak singularity integral formulas by vector-scalar identity, and finally the new singular integral techniques are successfully applied to compute the electromagnetic response of 3D CSEM with high precision. At last, for a conductive block buried in a less conductive half-space with a 100 m grounded wire, the total and secondary electric fields calculated by our algorithm are compared with those calculated by the integral equation method based on the secondary electric field, the finite element method based on magnetic vector potential and the DC resistivity forward modeling code (DCIP3D), respectively. The results show that four numerical solutions coincide well each other, and the algorithm suggested by this work is correct. Meanwhile, we conduct tests of this method on a sphere model and a complex geoelectric model, demonstrating that it is effective and capable of dealing with complicated subsurface anomaly bodies.
Key words: CSEM    IE    Unstructured    Forward modeling    Singular value integral    
0 引言

可控源电磁法(CSEM)是一种重要的地电磁勘探方法,通过引入人工场源从而大大改善了观测信号质量,使得其具有分辨率高及抗干扰能力强等特点.目前广泛用于石油、天然气、金属矿、地热、水文、环境等勘探领域(Kalscheuer et al., 2012; Li et al., 2017).

快速及高精度的CSEM正演算法是其反演的前提条件之一.目前,可控源电磁法中只有简单地电模型存在解析表达式,对于复杂的地电模型则需要涉及数值计算方法(Avdeev, 2005).目前常用的数值解法包含有限元、有限差分(有限体积)、表面积分法和体积分方程法.有限单元法(Jin, 2002, Yoshimura and Oshiman, 2002, Ren and Tang, 2014)能够模拟任意复杂起伏地形及地下异常体,但是其需要全区域剖分,未知数巨大,如Li和Key(2007)张继峰等(2009)Ren等(2013, 2014)、蔡红柱等(2015)杨军等(2015)李建慧等(2016)Yin等(2016)实现了MT和CSEM正演模拟.有限差分和有限体积法(谭捍东等, 2003; 孙怀凤等, 2013; 张烨等, 2012; 杨波等,2012; Jahandari and Farquharson, 2015; 陈辉等, 2016; 彭荣华等, 2016; Du et al., 2016)同样需要全区域剖分,未知数较多,并且大多数采用规则的六面体单元无法模拟复杂地下异常体(除Farquharson采用的非结构化Voronoi网格的有限体积方法外).边界单元法(Ren et al., 2012)仅仅适合于非常简单的地电模型.综上所述,有限元及有限体积等数值方法虽然已逐渐成为主流数值方法,但是,积分方程法与其他方法相比却有其自身不可比拟的优点:1)只对异常体进行剖分,未知数相对较少;2)积分方程法的数值结果具有半解析解的精度,常用于测试新开发的算法的正确性及求解精度;3)基于积分方程的反演算法只需要研究异常体区域,具有非常高的反演效率和精度(Kruglyakov et al., 2016).因此,开发高精度的CSEM积分方程正演算法具有一定的研究价值以及意义.

积分方程法20世纪六七十年代被应用到电磁法数值模拟中,如Hohmann等(1975)给出了基于六面体的奇异核积分技术.Weidelt(1975)给出了层状介质张量格林函数的积分表达式,从而奠定了积分方程法的基础,之后积分方程法获得了大量的研究结果,主要集中在地下规则异常体和近似方法的求解.例如Wannamaker等(1984)应用积分方程法求解了地下规则的异常体产生的大地电磁响应.张辉等(2005, 2006)采用规则的六面体网格实现了电偶源三维电磁积分方程正演.陈桂波等(2009a, 2009b)采用规则的六面体网格实现了三维电磁积分方程各向异性的研究.另外,Habashy等(1993)Torres-Verdín和Habashy(1994)Zhdanov等(2000, 2002)、王若等(2009)采用近似求解方法(Born近似方法和拟线性近似等)提高了积分方程法计算效率,但降低了求解精度.Born近似方法适用电导率差异较小及频率低情况,难适用于电导率高对比情况.上述研究结果大多数都是基于结构化的六面体网格来模拟地下复杂异常体,会导致由于网格离散来模拟异常体不准确产生的系统误差.更为重要的是,上述大部分研究工作常采用挖点法计算奇异值体积分,从而进一步降低了计算精度.

因此,本文提出了一套新的积分方程求解策略,实现了CSEM问题的高精度积分方程法求解.首先,直接从体积分方程出发,避免了近似求解方法可能带来的精度降低问题;其次,采用Tetgen(Si, 2007)开发的非结构化的四面体离散技术对异常体单元进行离散,降低了由网格单元拟合目标体不准确带来的系统误差.同时,系统地分析了奇异值积分问题,推导出一种消除积分奇异性的解析表达式,解决了并矢格林函数强奇异性问题,大大提高了求解精度.最后,以块状低阻异常体为测试模型,验证了本文算法的正确性;以球状和组合等复杂地电模型为例,测试本文算法的收敛性和处理地下复杂模型的能力.

1 三维CSEM积分方程 1.1 基本原理

假设时间因子为e-iωt,在准静态条件下,麦克斯韦方程组如下(Nabighian, 1988):

(1)

(2)

(3)

(4)

式中,E为电场,H为磁场,Js为激发源,σ为电导率,μ为磁导率,ε为介电常数,ω=2πf为角频率,f为频率.

利用公式(1)和(2)可知:

(5)

利用电张量格林函数以及公式(5),可以得出传统的电磁积分方程公式:

(6)

其中,Ω为异常体区域,Ωs为场源区域的积分.(6)式既是观测点数据求取的方程亦是进行求解异常电场的方程.本文采用的并矢Green函数为均匀半空间下的频率域张量格林函数,其具体形式由Hohmann(1975)给出(见附录).从公式(6)可以看出,任意点的总场是由源Js激发的背景场和异常电导率(σ-σ0)产生的散射场组成,图 1a展示了3D CSEM工作原理.

图 1 模型示意图 (a) CSEM勘探原理;(b)四面体单元. Fig. 1 Sketch of the model (a) Principle of controlled source electromagnetic method (CSEM); (b) Tetrahedral element.

将地下异常体离散成M个四面体子单元(如图 1b所示),令每个单元内部电场为常值(即零次基函数),其值取为单元中心ri(i=1, 2, …, M)的值,并对公式(6)进行离散化,可得:

(7)

将观测点放在异常体区域,可得到3M个线性方程,公式(7)写成矩阵形式为:

(8)

其中,E=[E1, E2, …, E3M]TEb=[E1b, E2b, …, E3Mb]T分别为由离散单元中心未知场以及背景场组成的3M列向量,T表示矩阵转置,矩阵A为离散单元的系数,表达式为:

(9)

其中,I表示单位张量.求解方程组(8)可以得到地下异常体各离散单元的场值,然后通过公式(7)可以获取空间内任意点场值.

1.2 并矢Green函数奇异核积分的新解析表达式

计算公式(8)中的矩阵A的系数时,需要涉及并矢格林函数积分计算.然而,当r=r′时,并矢Green函数具有很强的奇异性,导致系数矩阵A主对角元素不准确,从而影响计算结果的精度.首先,我们将矩阵A的系数分解成电流以及电荷两部分(Harrington,1968),矩阵A的系数的奇异积分项仅仅存在于电流及电荷的一次场项(见附录公式(A4)和(A5)),其他积分项不存在奇异性.传统的奇异值积分处理技术常常采用挖点法计算获得的奇异值体积分,从而降低了计算精度.对于四面体单元而言,该方法会造成严重的数值误差,因此基于作者在重力方法的最新研究成果(Ren et al., 2017a, 2017b, 2017c; Jiang et al., 2017),开发一种格林函数的奇异性去除办法,实现了任意多面体下格林函数积分解析计算,大大提高了奇异值积分的计算精度及实用性.

含奇异性及强奇异性积分项的电流及电荷的积分公式为:

(10)

(11)

其中,R=|r-r′|,Ω为积分区域.

r=r′时,公式以及存在强奇异性,为了能够处理此类强奇异体积分,对公式(10)及(11)进行如下变换,可得:

(12)

(13)

其中,, , , .

由于是连续光滑函数,因此,, ψ1是连续可导函数,其可以采用具有代数精度的高斯数值积分进行计算.因此,奇异积分项转移到式ψ2中,为了能够对其进行奇异性去除,即对式(12)、(13)分别进行矢量恒等式及梯度定理处理,得:

(14)

(15)

(16)

其中,Ti为四面体第i个面,为四面体第i个面的单位外法向矢量.令di=[(r-r′)·]表示场点到四面体单元(源点单元)各面中心有向距离.从包含奇异性积分公式(14)及(16)可知,其主要涉及到这两类奇异面积分计算.

为了去除单元积分的奇异性,构建一个局部坐标系,如图 2所示.点o为点r(观测点)在三角单元Tj上投影点,β(o)=∑β(o)j为投影点o与每一边Cj两端点的夹角之和,分别为边Cj的单位法向和切向. hi为点r到多边形Tj的距离,即hi=(r-r′)· 为多边形Tj的面法向,mj是投影点o到边Cj的垂向距离,即mj=(o-r′)·ρ表示点o到边Cj上任一点距离,即ρ=|o-r′|,v0v1分别表示边Cj两个端点,假设中间变量s=(r′-os0s1分别是边Cj顶点v0v1中间变量,R0R1分别表示观测点到边Cj两个端点的距离,R表示观测点与点o到边Cj的投影点r之间的距离,即R表示观测点到原点的距离,即.根据上述公式导出了关于的积分恒等式为:

(17)

图 2r和多边形Tj的边Cj的关系 Fig. 2 llustration of the geometrical relationship between the observation site r and an Cj of a polygon Tj. The point o is the projection point of r onto the polygon Tj

其中,∇′s为2D局部坐标下面散度算子.假设(u′, v′)为r′的局部坐标系,o为2D局部坐标中心点(0, 0),.其相关证明如下:

(18)

为了处理R→0弱奇异积分,将任意多边形Tj的积分分解为两个部分,一部分为奇异部分,即以o点为中心的无限小区域Oεε→0;另一部分为除奇异部分的其他区域Ωi-Oε,此时的积分公式可写为如下:

(19)

式中,β(o)为投影点o在多边形内Tj的固体角,固体角采用下面公式进行计算:

(20)

式中为点o到点r的单位矢量.从公式(20)可知,当点o位于多边形内时β(o)=2π;若在多边形边上,β(o)=π;若在多边形拐角点上,β(o)=Θ;若在多边形外部,β(o)=0.

q=±1时,公式(19)中的Ajq+2的解析表达式的具体形式如下表示.当q=-1时,

(21)

q=1时,

(22)

同理,对于任意多边形,可推导出的积分计算,其表达式为:

(23)

其中,.

公式(23)中同样涉及到q=±1时的积分计算,当q=-1时,

(24)

q=1时,

(25)

上述的公式推导得出了并矢格林函数所包含的强奇异性积分公式解析表达式,该表达式能够降低由强奇异积分带来的数值积分误差等问题.

2 算例分析

本文采用的测试平台内存为7.7 G,intelR CoreTM i5-4590 CPU@ 3.30 GHz.

2.1 奇异值积分正确性验证

为了验证本文推导解析表达式的正确性,设计如图 3所示的四面体单元,针对矢量χ2(k1=i),观测点为四面体单元的重心,对比分析本文推导出的解析表达式与高阶高斯数值积分(N为阶数)计算结果(标量解).从表 1中可知,随着阶数N的增加,χ2的实部值在小数点后7位变化相对较小,而虚部达到了小数点后12位,说明高斯数值积分是逐渐收敛的.与解析表达式结果对比,实部值的相对误差为1.860/0000,虚部值的精度更高.结果表明本文推导出的解析表达式计算结果正确且精度高,并为本文后续工作奠定扎实的基础.

图 3 四面体单元 Fig. 3 Tetrahedral element
表 1 χ2高斯数值积分与解析表达式对比 Table 1 Comparison between Gauss′s numerical integration and analytic expression of χ2
2.2 算法验证以及收敛性测试 2.2.1 与其他积分方程和有限元结果对比

首先,设计块状低阻异常体进行程序正确性测试.图 4所示的块状低阻异常体模型,电导率为0.2 S·m-1的异常体被置于均匀半空间模型中,背景电导率为0.02 S·m-1,异常体尺寸为120 m×200 m× 400 m,中心点坐标为(1000 m, 0 m, 300 m).沿着x方向布设有限长导线源(其可以看成若干个偶极子的叠加),源的长度为100 m,源的中心坐标为(50 m, 0 m, 0 m),发射电流1 A,发射频率为3 Hz.沿着x方向布设一条测线,测线起始位置(400 m, 0 m, 0 m),终点位置为(1400 m, 0 m, 0 m).采用开源程序Tetgen(Si, 2007)非结构四面体网格离散技术进行单元剖分,剖分单元数为1155,未知数为3465,内存消耗为189 M,计算时间2143 s(直接求解).

图 4 块状低阻异常体模型 Fig. 4 Massive model of low-resistivity anomalies. A 100 m grounded wire source is on the air-ground interface. The size of the conductive prism is 120 m×200 m×400 m in the x-, y-, and z-directions

将本文的电场Ex分量的总场及二次场计算结果与采用六面体网格(512单元)的积分方程结果(Farquharson and Oldenburg, 2002)、非结构化有限元(713542单元)的结果(Ansari and Farquharson, 2014)和程序DCIP3D(Li and Oldenburg, 1999)的计算结果进行对比.图 56分别表示Ex分量的总场及二次场对比曲线.从图 5中可以看出,本文的结果与前人的结果高度吻合, 验证了本文算法的正确性.同样地,图 6二次电场结果可知,本文的结果与基于磁矢势的Helmholtz方程的有限元的结果具有高度吻合性,与Farquharson的积分方程结果存在一定误差,但误差相对较小,进一步验证本文算法的正确性.图 5的结果还可以了解到在异常体上方电场的实部与虚部存在下凹的现象,这说明异常体存在使得电流发生明显的变化,使得电流趋于低阻体方向流动,导致流经地表电流密度降低.

图 5 块状低阻体置于均匀半空间下的地电模型的Ex总场分量 (a)实部; (b)虚部. Fig. 5 Components of the total Ex field from a geoelectric model with a massive low-resistivity body in a uniform half-space (a) Real part; (b) Imaginary part.
图 6 块状低阻体置于均匀半空间下的Esx二次场分量 (a)实部; (b)虚部. Fig. 6 Components of the secondary field Esx from a geoelectric model with a massive low-resistivity body in a uniform half-space (a) Real part; (b) Imaginary part.
2.2.2 算法收敛性测试

设置如图 7a所示的球状异常体模型分析本文算法收敛性情况.球体的半径为225 m,球体的中心位置为(0 m, 0 m, 450 m).沿x方向放置长接地导线,导线长度为1 km,发射电流为1 A,源的中心坐标为(0 m, 6000 m, 0 m).以坐标原点为中心,沿着x轴布设19个测点,测点x方向的坐标范围为(-1000~1000 m),yz等于0 m.发射频率为16 Hz.球状异常体被置于均匀半空间中,电导率为1 S·m-1,背景电导率为0.01 S·m-1. 图 7b展示了不同异常体单元剖分数视电阻率曲线变化图.从图 7b结果中可知,视电阻率曲线能够较明显体现低阻异常响应,同时随着网格数的增加,视电阻率值未发现明显的变化,说明了本文编制的3D CSEM程序收敛且正确.

图 7 (a) 球状异常体模型以及(b)不同网格剖分收敛性分析 Fig. 7 (a) Sphere model of an anomaly body and (b) convergence analysis for different grid subdivision
2.3 倾斜板状的影响

设计如图 8所示的倾斜板状异常体模型,异常体被置于电导率为0.01 S·m-1的均匀半空间中,异常体的电导率分别为0.1 S·m-1和0.333 S·m-1,异常区尺度为400 m×400 m×800 m,x方向的坐标范围为:顶部[-200 m, 200 m], 底部[-400 m, 0 m];y方向的坐标范围为:[-200 m, 200 m];z方向的坐标范围为:[100 m, 900 m].沿着x方向布设电偶源,偶极子源长度为10 m,源的中心坐标为(0 m, 1000 m, 0 m),发射电流100 A,激发频率分别为1 Hz和32 Hz, 沿着y方向布设五条测线(y=100 m, 0 m, -100 m, -200 m, -300 m),每条测线的长度为2000 m,每条测线均29个测点,具体模型如图 8所示.从图 9所示的结果中可以看出,电场分量Ex随着收发距增大而逐渐减小,这一现象也满足电磁场能量随距离增大而逐渐衰减的规律.地下不存在板状异常体(均匀半空间)时,电场Ex振幅关于y轴对称;当地下存在倾斜板状体时,电场Ex振幅不关于y轴对称,其明显受到倾斜高导板状异常体的影响.从不同高导异常体电场Ex振幅等值线图 9b9c9e9f可知,当电导率为0.333 S·m-1电场Ex振幅等值线凹陷程度明显大于电导率为0.1 S·m-1.

图 8 倾斜板状低阻异常体 (a) x-z平面断面; (b) x-y平面断面. Fig. 8 Model of an inclined plate-like low-resistivity body (a) x-z cross section; (b) x-y cross section.
图 9 倾斜板状异常体模型的电场总场Ex分量振幅 1)激发频率为1 Hz: (a)均匀半空间; (b) 0.1 S·m-1的倾斜板状异常体; (c) 0.333 S·m-1的倾斜板状异常体. 2)激发频率为32 Hz: (d)均匀半空间; (e) 0.1 S·m-1的倾斜板状异常体; (f) 0.333 S·m-1的倾斜板状异常体. Fig. 9 X-component amplitude of the total electric field Ex for a model of an inclined plate-like anomaly body 1) Excitation frequency 1 Hz: (a) Homogeneous half-space; (b) Inclined plate-like anomaly body of 0.1 S·m-1; (c) Inclined plate-like anomaly body of 0.333 S·m-1. 2) Excitation frequency 32 Hz: (d) Homogeneous half-space; (e) Inclined plate-like anomaly body of 0.1 S·m-1; (f) Inclined plate-like anomaly body of 0.333 S·m-1.
2.4 异常体组合模型

设计如图 10所示的组合异常体模型,异常体被置于电导率为0.01 S-1·m的均匀半空间中,图 10a中从左到右的异常体的电导率分别为0.001 S-1·m, 0.1 S-1·m, 0.002 S-1·m, 异常区尺度为600 m×200 m×400 m,x方向的坐标范围为:[-300 m, 300 m];y方向的坐标范围为:[-100 m, 100 m];z方向的坐标范围为:[200 m, 600 m].在x轴线以及y轴线上分别布设有限长导线源,源长为1000 m,源的坐标分别为(5000 m, 0 m, 0 m),(0 m,5000 m, 0 m),激发频率分别为1 Hz,2 Hz和4 Hz,激发电流为1 A,测量了主剖面y=0 m,共计算了53测点.图 11, 12分别表示源点坐标为(5000 m, 0 m, 0 m)激发时,电场Ex总场与二次场的实部与虚部曲线.由图 11可以看出,在异常体正上方曲线发生明显弯曲的变化,并在低阻异常体正上方曲线向下弯曲,体现了低阻异常体吸引电流的现象,而在高阻异常体的正上方曲线有明显向上凸的趋势,说明了高阻异常体具有排斥电流的现象.同时随着测点距离源越远,电场幅值逐渐衰减.而从纯二次电场Esx曲线(图 12)中可知,在异常区上方出现明显的三个极值点,极值点的位置区域能够很好地对应各异常体区的正上方.从图 12中可知,受场源的影响,在异常体电导率为0.002 S-1·m的正上方电场峰值明显要大于其他两个区域,这一现象与源点位于(0 m,5000 m, 0 m)激发产生的二次电场Esx的曲线存在明显的不同.图 13中可知,异常体区域同样出现明显的三个异常极值点,且可与异常体相对应,与图 12不同之处为异常体的电导率在0.001 S-1 · m得到的幅值最大,然后逐渐依次减小.

图 10 组合模型示意图 (a) x-z平面断面; (b) x-y平面断面. Fig. 10 Sketch of combined anomaly model (a) x-z cross section; (b) x-y cross section.
图 11Js的坐标为(5000 m, 0 m, 0 m)激发时,组合异常体Ex总场分量 (a)实部; (b)虚部. Fig. 11 Components of the total electric field Ex of a combined anomaly body with source Js at (5000 m, 0 m, 0 m) (a) Real part; (b) Imaginary part.
图 12Js的坐标(5000 m, 0 m, 0 m)为激发时,组合异常体二次场Esx分量 (a)实部; (b)虚部. Fig. 12 Components of the secondary electric field Esx of a combined anomaly body with source Js at (5000 m, 0 m, 0 m) (a) Real part; (b) Imaginary part.
图 13Js的坐标(0 m, 5000 m, 0 m)为激发时,组合异常体二次场Esx分量 (a)实部; (b)虚部. Fig. 13 Components of the secondary electric field Esx of a combined anomaly body with source Js at (0 m, 5000 m, 0 m) (a) Real part; (b) Imaginary part.
3 结论

本文采用非结构化四面体网格对地下异常体进行离散,推导出了全新的任意非规则四面体(适合任意多面体)单元的张量格林函数的积分解析解表达式,实现了3D CSEM高精度的正演模拟.通过设计一系列地电模型进行测试,获得不同地电模型下三维可控电磁法正演响应,同时与公开的程序结果进行对比,得出以下几点认识:

(1) 本文给出任意多面体的张量格林函数奇异性去除的解析表达式与高阶(400阶)高斯积分结果具有高度吻合,误差精度达到万分之一数量级,高斯数值积分随着积分阶数的增加求解效率将大大降低,而本文的解析解表达式不存在此类问题.同时,相比于规则六面体网格奇异性处理技术(扣点法)传统积分方程法,本文方法更具有一般性且求解精度较高.

(2) 文中采用了非结构化四面体的网格剖分技术离散复杂异常体,弥补了传统规则网格不能准确模拟复杂异常体的不足,提高了计算精度.

(3) 本文的研究成果为加速积分方程法的求解速度、精确计算地下复杂模型的CSEM响应等后续研究提供了基础;同时,文中的奇异值积分技术还可以扩展到大地电磁、直流电阻率及地震波的积分求解中;另外,本文所开发的长导线CSEM程序可为同行测试新开发的算法提供有力的参考.

致谢

感谢纽芬兰纪念大学地球科学系Colin.G Farquharson教授课题组的Hormoz Jahandari与SeyedMasoud Ansari两位博士提供算法验证数据,同时特别感谢两位审稿专家中肯的修改意见,使得本文内容更加完善.

附录

频率域半空间张量Green函数的表达式可以写成一次场GP和二次场GS和的形式(Hohmann, 1975; 鲍光淑等, 1999):

(A1)

(A2)

(A3)

其中,

(A4)

(A5)

(A6)

(A7)

这里, , , , , ,

, 其中,r(x, y, z)为场点坐标,r′(x′, y′, z′)为源点坐标.

参考文献
Ansari S, Farquharson C G. 2014. 3D finite-element forward modeling of electromagnetic data using vector and scalar potentials and unstructured grids. Geophysics, 79(4): E149-E165. DOI:10.1190/geo2013-0172.1
Avdeev D B. 2005. Three-dimensional electromagnetic modelling and inversion from theory to application. Surveys in Geophysics, 26(6): 767-799. DOI:10.1007/s10712-005-1836-x
Bao G S, Zhang B X, Jing R Z, et al. 1999. Numerical modeling for 3-D EM responses using integral equation solution. Journal of Central South University, 30(5): 472-474.
Cai H Z, Xiong B, Zhdanov M. 2015. Three-dimensional marine controlled-source electromagnetic modelling in anisotropic medium using finite element method. Chinese Journal of Geophysics, 58(8): 2839-2850. DOI:10.6038/cjg20150818
Chen G B, Wang H N, Yao J J, et al. 2009a. Three-dimensional numerical modeling of marine controlled-source electromagnetic responses in a layered anisotropic seabed using integral equation method. Acta Physica Sinica, 58(6): 3848-3857.
Chen G B, Wang H N, Yao J J, et al. 2009b. Modeling of electromagnetic responses of 3-D electrical anomalous body in a layered anisotropic earth using integral equations. Chinese Journal of Geophysics, 52(8): 2174-2181.
Chen H, Yin C C, Deng J Z. 2016. A finite-volume solution to 3D frequency-domain electromagnetic modelling using Lorenz-gauged magnetic vector and scalar potentials. Chinese Journal of Geophysics, 59(8): 3087-3097. DOI:10.6038/cjg20160831
Du H K, Ren Z Y, Tang J T. 2016. A finite-volume approach for 2D magnetotellurics modeling with arbitrary topographies. Studia Geophysica et Geodaetica, 60(2): 332-347. DOI:10.1007/s11200-014-1041-9
Farquharson C G, Oldenburg D W. 2002. An integral equation solution to the geophysical electromagnetic forward-modelling problem. Methods in Geochemistry and Geophysics, 35: 3-19. DOI:10.1016/S0076-6895(02)80083-X
Habashy T M, Groom R W, Spies B R. 1993. Beyond the Born and Rytov approximations:A nonlinear approach to electromagnetic scattering. Journal of Geophysical Research:Solid Earth, 98(B2): 1759-1775. DOI:10.1029/92JB02324
Harrington R. 1968. Field Computation by Moment Methods. Macmillan.
Hohmann G W. 1975. Three-dimensional induced polarization and electromagnetic modeling. Geophysics, 40(2): 309-324. DOI:10.1190/1.1440527
Jahandari H, Farquharson G C. 2015. Finite-volume modelling of geophysical electromagnetic data on unstructured grids using potentials. Geophysical Journal International, 202(3): 1859-1876. DOI:10.1093/gji/ggv257
Jiang L, Zhang J Z, Feng Z B. 2017. A versatile solution for the gravity anomaly of 3D bodies with density-contrast variation with depth. Geophysics, 82(4): G77-G86. DOI:10.1190/geo2016-0394.1
Jin J M. 2002. The Finite Element Method in Electromagnetics. (2nd ed). New York: Wiley.
Kalscheuer T, Hübert J, Kuvshinov A, et al. 2012. A hybrid regularization scheme for the inversion of magnetotelluric data from natural and controlled sources to layer and distortion parameters. Geophysics, 77(4): E301-E315. DOI:10.1190/geo2012-0018.1
Kruglyakov M, Geraskin A, Kuvshinov A. 2016. Novel accurate and scalable 3-D MT forward solver based on a contracting integral equation method. Computers & Geosciences, 96: 208-217.
Li G, Xiao X, Tang J T, et al. 2017. Near-source noise suppression of AMT by compressive sensing and mathematical morphology filtering. Applied Geophysics, 14(4): 581-590. DOI:10.1007/s11770-017-0645-6
Li J H, Farquharson C G, Hu X Y, et al. 2016. A vector finite element solver of three-dimensional modelling for a long grounded wire source based on total electric field. Chinese Journal of Geophysics, 59(4): 1521-1534. DOI:10.6038/cjg20160432
Li Y, Oldenburg D W, Shekhtman E. 1999. DCIP3D: A program library for forward modelling and inversion of DC resistivity and induced polarization data over 3D structures. Report of the INDI Consortium, University of British Columbia.
Li Y G, Key K. 2007. 2D marine controlled-source electromagnetic modeling:Part 1-An adaptive finite-element algorithm. Geophysics, 72(2): WA51. DOI:10.1190/1.2432262
Nabighian M N. 1988. Electromagnetic methods in applied geophysics: theory. Tulsa, USA: Society of Exploration Geophysicists. http://ci.nii.ac.jp/ncid/BA29889651
Peng R H, Hu X Y, Han B, et al. 2016. 3D frequency-domain CSEM forward modeling based on the mimetic finite-volume method. Chinese Journal of Geophysics, 59(10): 3927-3939. DOI:10.6038/cjg20161036
Ren Z Y, Kalscheuer T, Greenhalgh S, et al. 2012. Boundary element solutions for broad-band 3-D geo-electromagnetic problems accelerated by an adaptive multilevel fast multipole method. Geophysical Journal International, 192(2): 473-499.
Ren Z Y, Kalscheuer T, Greenhalgh S, et al. 2013. A goal-oriented adaptive finite-element approach for plane wave 3-D electromagnetic modelling. Geophysical Journal International, 194(2): 700-718. DOI:10.1093/gji/ggt154
Ren Z Y, Kalscheuer T, Greenhalgh S, et al. 2014. A finite-element-based domain-decomposition approach for plane wave 3D electromagnetic modeling. Geophysics, 79(6): E255-E268. DOI:10.1190/geo2013-0376.1
Ren Z Y, Tang J T. 2014. A goal-oriented adaptive finite-element approach for multi-electrode resistivity system. Geophysical Journal International, 199(1): 136-145. DOI:10.1093/gji/ggu245
Ren Z Y, Chen C J, Pan K J, et al. 2017a. Gravity anomalies of arbitrary 3D polyhedral bodies with horizontal and vertical mass contrasts. Surveys in Geophysics, 38(2): 479-502. DOI:10.1007/s10712-016-9395-x
Ren Z Y, Tang J T, Kalscheuer T, et al. 2017b. Fast 3-D large-scale gravity and magnetic modeling using unstructured grids and an adaptive multilevel fast multipole method. Journal of Geophysical Research:Solid Earth, 122(1): 79-109. DOI:10.1002/2016JB012987
Ren Z Y, Zhong Y Y, Chen C J, et al. 2017c. Gravity anomalies of arbitrary 3D polyhedral bodies with horizontal and vertical mass contrasts up to cubic order. Geophysics, 83(1): G1-G13.
Si H. 2007. TetGen: A quality tetrahedral mesh generator and three-dimensional delaunay triangulator, http://tetgen.berlios.de, accessed 1 January 2014.
Sun H F, Li X, Li S C, et al. 2013. Three-dimensional FDTD modeling of TEM excited by a loop source considering ramp time. Chinese Journal of Geophysics, 56(3): 1049-1064. DOI:10.6038/cjg20130333
Tan H D, Yu Q F, Booker J, et al. 2003. Magnetotelluric three-dimensional modeling using the staggered grid finite difference method. Chinese Journal of Geophysics, 46(5): 705-711.
Torres-Verdín C, Habashy T M. 1994. Rapid 2.5-dimensional forward modeling and inversion via a new nonlinear scattering approximation. Radio Science, 29(4): 1051-1079.
Wang R, Di Q Y, Wang M Y, et al. 2009. Research on the effect of 3D body between transmitter and receivers on CSAMT response using integral equation method. Chinese Journal of Geophysics, 52(6): 1573-1582.
Wannamaker P E, Hohmann G W, SanFilipo W A. 1984. Electromagnetic modeling of three-dimensional bodies in layered earths using integral equations. Geophysics, 49(1): 60-74. DOI:10.1190/1.1441562
Weidelt P. 1975. Electromagnetic induction in three-dimensional structures. J. Geophys., 41: 85-109.
Yang B, Xu Y X, He Z X, et al. 2012. 3D frequency-domain modeling of marine controlled source electromagnetic responses with topography using finite volume method. Chinese Journal of Geophysics, 55(4): 1390-1399. DOI:10.6038/j.issn.0001-5733.2012.04.035
Yang J, Liu Y, Wu X P. 2015. 3D simulation of marine CSEM using vector finite element method on unstructured grids. Chinese Journal of Geophysics, 58(8): 2827-2838. DOI:10.6038/cjg20150817
Yin C C, Zhang B, Liu Y H, et al. 2016. A goal-oriented adaptive finite-element method for 3D scattered airborne electromagnetic method modeling. Geophysics, 81(5): E337-E346. DOI:10.1190/geo2015-0580.1
Yoshimura R, Oshiman N. 2002. Edge-based finite element approach to the simulation of geoelectromagnetic induction in a 3-D sphere. Geophysical Research Letters, 29(3): 9-1.
Zhang H, Li T L, Dong R X, et al. 2005. Computation of Green's tensor integrals for electromagnetic problems using Gaussian quadrature and continued fraction. Progress in Geophysics, 20(3): 667-670.
Zhang H, Li T L, Dong R X. 2006. Modeling 3-D electromagnetic responses of the electric dipole using volume integral equation method. Progress in Geophysics, 21(2): 386-390.
Zhang J F, Tang J T, Yu Y, a l. 2009. Three dimensional controlled source electromagnetic numerical simulation based on electric field vector wave equation using finite element method. Chinese Journal of Geophysics, 52(12): 3132-3141.
Zhang Y, Wang H N, Tao H G, et al. 2012. Finite volume algorithm to simulate 3D responses of multi-component induction tools in inhomogeneous anisotropic formation based on coupled scalar-vector potentials. Chinese Journal of Geophysics, 55(6): 2141-2152. DOI:10.6038/j.issn.0001-5733.2012.06.036
Zhdanov M S, Dmitriev V I, Fang S, et al. 2000. Quasi-analytical approximations and series in electromagnetic modeling. Geophysics, 65(6): 1746-1757. DOI:10.1190/1.1444859
Zhdanov M S, Portniaguine O, Hursan G. 2002. Compression in 3-D integral equation modeling. Methods in Geochemistry and Geophysics, 35: 21-42. DOI:10.1016/S0076-6895(02)80084-1
鲍光淑, 张碧星, 敬荣中, 等. 1999. 三维电磁响应积分方程法数值模拟. 中南大学学报(自然科学版), 30(5): 472–474.
蔡红柱, 熊彬, ZhdanovM. 2015. 电导率各向异性的海洋电磁三维有限单元法正演. 地球物理学报, 58(8): 2839–2850. DOI:10.6038/cjg20150818
陈桂波, 汪宏年, 姚敬金, 等. 2009a. 各向异性海底地层海洋可控源电磁响应三维积分方程法数值模拟. 物理学报, 58(6): 3848–3857.
陈桂波, 汪宏年, 姚敬金, 等. 2009b. 用积分方程法模拟各向异性地层中三维电性异常体的电磁响应. 地球物理学报, 52(8): 2174–2181.
陈辉, 殷长春, 邓居智. 2016. 基于Lorenz规范条件下磁矢势和标势耦合方程的频率域电磁法三维正演. 地球物理学报, 59(8): 3087–3097. DOI:10.6038/cjg20160831
李建慧, FarquharsonC G, 胡祥云, 等. 2016. 基于电场总场矢量有限元法的接地长导线源三维正演. 地球物理学报, 59(4): 1521–1534. DOI:10.6038/cjg20160432
彭荣华, 胡祥云, 韩波, 等. 2016. 基于拟态有限体积法的频率域可控源三维正演计算. 地球物理学报, 59(10): 3927–3939. DOI:10.6038/cjg20161036
孙怀凤, 李貅, 李术才, 等. 2013. 考虑关断时间的回线源激发TEM三维时域有限差分正演. 地球物理学报, 56(3): 1049–1064. DOI:10.6038/cjg20130333
谭捍东, 余钦范, BookerJ, 等. 2003. 大地电磁法三维交错采样有限差分数值模拟. 地球物理学报, 46(5): 705–711.
王若, 底青云, 王妙月, 等. 2009. 用积分方程法研究源与勘探区之间的三维体对CSAMT观测曲线的影响. 地球物理学报, 52(6): 1573–1582.
杨波, 徐义贤, 何展翔, 等. 2012. 考虑海底地形的三维频率域可控源电磁响应有限体积法模拟. 地球物理学报, 55(4): 1390–1399. DOI:10.6038/j.issn.0001-5733.2012.04.035
杨军, 刘颖, 吴小平. 2015. 海洋可控源电磁三维非结构矢量有限元数值模拟. 地球物理学报, 58(8): 2827–2838. DOI:10.6038/cjg20150817
张辉, 李桐林, 董瑞霞, 等. 2005. 利用高斯求积和连分式展开计算电磁张量格林函数积分. 地球物理学进展, 20(3): 667–670.
张辉, 李桐林, 董瑞霞. 2006. 体积分方程法模拟电偶源三维电磁响应. 地球物理学进展, 21(2): 386–390.
张继锋, 汤井田, 喻言, 等. 2009. 基于电场矢量波动方程的三维可控源电磁法有限单元法数值模拟. 地球物理学报, 52(12): 3132–3141. DOI:10.3969/j.issn.0001-5733.2009.12.023
张烨, 汪宏年, 陶宏根, 等. 2012. 基于耦合标势与矢势的有限体积法模拟非均匀各向异性地层中多分量感应测井三维响应. 地球物理学报, 55(6): 2141–2152. DOI:10.6038/j.issn.0001-5733.2012.06.036