地球物理学报  2020, Vol. 63 Issue (9): 3520-3533   PDF    
引用本文
张继锋, 刘寄仁, 冯兵, 等. 2020. 三维陆地可控源电磁法有限元模型降阶快速正演. 地球物理学报, 63(9): 3520-3533, doi: 10.6038/cjg2020N0452.  
Zhang J F, Liu J R, Feng B, et al. 2020. Fast forward modeling of the 3D land controlled-source electromagnetic method based on model reduction. Chinese J. Geophys. (in Chinese), 63(9): 3520-3533, doi: 10.6038/cjg2020N0452.  
三维陆地可控源电磁法有限元模型降阶快速正演
张继锋, 刘寄仁, 冯兵, 郑一安     
长安大学地质工程与测绘学院, 西安 710054
收稿日期 2019-12-21, 2020-03-29 收修定稿
基金项目: 国家重点研发计划项目（2017YFC0602200）子课题"航空地球物理综合处理解释方法研究及软件开发"（2017YFC0602202）资助
第一作者简介: 张继锋, 男, 1978年生, 副教授, 主要从事电磁法理论及正反演研究.E-mail:zjf0201@126.com
通讯作者: 冯兵, 男, 1962年生, 副教授, 主要从事电磁勘探理论及应用研究.E-mail:fb8005@126.com
摘要：三维陆地可控源电磁法有限元快速正演的主要瓶颈在于多频率大型稀疏方程组求解问题.本文引入一种基于模型降阶的Krylov子空间投影算法，推导了有限元刚度矩阵的模型降阶形式，构建了频率域传递函数；采用标准正交向量序列，构建一个远远小于有限元刚度矩阵维度的矩阵，该矩阵与频率无关，通过一次模型降阶即可实现多频点有限元方程快速求解.采用基于电场的变分方程，加入散度校正条件，以消除伪解；引入伪δ函数，消除了源点的奇异性，可适用于复杂背景模型三维有限元数值模拟，并为多源的求解奠定了基础；以层状介质模型解析解为标准，通过和基于Pardiso直接求解器的有限元算法（3DFEM）进行比较，模型降阶法计算时间小于前者的1/10，平均相对误差在1.72%，在满足精度要求下，实现了高效率三维有限元数值求解；分别设计了横向高低阻模型和纵向高低阻模型，分析了从近区到远区电场和卡尼亚视电阻率的变化规律，假极值的表现特征，阴影效应的影响等，从而也验证了该算法的正确性.最后，建立了一个地层陷落柱模型，通过模型降阶有限元正演模拟，发现视电阻率断面图在陷落柱上方出现"凹陷"，与模型设计吻合，表明该算法对复杂地层模拟具有同样的适用性.
关键词: 可控源电磁法      模型降阶      有限元      电磁测深     
Fast forward modeling of the 3D land controlled-source electromagnetic method based on model reduction
ZHANG JiFeng, LIU JiRen, FENG Bing, ZHENG YiAn     
School of Geology Engineering and Geomatics, Chang'an University, Xi'an 710054, China
Abstract: The main bottleneck of fast forward modeling of 3-D land controlled-source electromagnetics based on the finite element method lies in solving multi-frequency large-scale sparse equations. In this paper, a Krylov subspace projection algorithm based on model reduction is introduced, a model reduction form of the finite element stiffness matrix is derived, and a frequency domain transfer function is constructed. A standard orthogonal vector sequence is used to construct a matrix that is much smaller than the finite element stiffness matrix, which is unrelated to frequency. Through model reduction, it is possible to quickly solve these multiple frequency finite element equations. The variational equation based on the electric field is used, and divergence correction conditions are added to eliminate spurious solutions. The pseudo-delta function is introduced to eliminate the singularity of the source point, which can be applied to 3D finite element numerical simulation of complex models. Such an algorithm may be regarded as the basic principle for a solution to the problem of multiple sources. Using the analytical solution of the layered medium model as the standard, comparison shows that the averaged percentage computing error of the proposed method is 1.72%. Compared with the 3DFEM algorithm using the Pardiso direct solver, the CPU time consumption of the proposed method is reduced by 1/10 of the original one. Therefore, with a required accuracy, it is possible to achieve highly efficient 3-D finite element numerical solutions. In the numerical simulation, with the horizontal and vertical high-low resistance models, we analyze the changes of the electric field and the Cagniard apparent resistivity from the near to the far region, the performance characteristics of false extremes, and the effects of shadow effects, which verify the correctness of the proposed algorithm. Finally, we apply the proposed algorithm to simulation of a stratum collapse column model. The results show that the apparent resistivity cross section looks like a "sag" above the collapse column, which is consistent with the model design, and further indicates that the proposed algorithm can be applied to real complex geoelectric structure.
Keywords: Controlled-source electromagnetic method    Model reduction    Finite element method    Electromagnetic sounding    
0 引言

为了克服了大地电磁场源信号微弱难以观测的缺点, Strangway和Myron Goldstein在70年代提出可控源音频大地电磁法.由于场源的加入, 该方法抗干扰能力大大增强, 在矿产资源、天然气、地热、水文工程、环境保护、煤田采空区富水性调查等各个方面(汤井田和何继善, 2005)得到了广泛的应用, 成为一种重要的地球物理勘探技术(彭荣华等, 2018).陆地可控源电磁法常采用接地电偶极子作为发射源, 通过发射一系列不同频率电磁波, 在远区观测相互正交的电磁场分量, 沿用大地电磁测深中卡尼亚电阻率计算公式, 达到电磁测深目的.海洋可控源一般是通过拖曳电偶极源方式进行海上作业, 由于海水电导率很低, 一般发射频率范围在0.1~10 Hz(Constable, 2010), 频点相对也较少, 本文主要以陆地可控源宽频电磁方法作为主要研究内容.

三维正演主要以数值计算方法为主, 目前应用于频率域电磁法的三维正演方法主要是有限体积法(张烨等, 2012; 周建美等, 2018; 彭荣华等, 2016, 2018), 有限差分法(Haber et al., 2000; Newman and Alumbaugh, 2002; Hou et al., 2006; Streich, 2009), 积分方程法(Zhdnov et al., 2006; 任政勇等, 2017; 汤井田等, 2018)和有限元法(张继锋等, 2009; 蔡红柱等, 2015; Ren et al., 2013, 2014; 李建慧, 2016; Um et al., 2017; 刘颖等, 2017; 王亚璐等, 2017).以上几种方法各有优缺点, 具体请参考Avdeev (2005)和Everett (2011), 其中有限元方法因其网格剖分灵活, 可模拟任意复杂形状地电体, 近年来得到了广泛的研究和应用.

由于源的加载, 可控源大地电磁法为有限元数值模拟带来了一定困难.目前主要有两种方法克服源的奇异性问题, 第一种方法是基于二次场法(张林成等, 2017; 王亚璐等, 2017; 彭荣华等, 2018; Cai et al., 2014), 该方法在计算背景场时耗费大量的时间和计算机内存, 且很难确定复杂模型的背景场.另一种方法是采用总场法(邱长凯等, 2018; Cai et al., 2017), 可直接把电偶极源用伪δ函数近似(Mitsuhata, 2000; 张继锋, 2009; 2013), 在场源附近可适当加密网格以减少场源引起的奇异性.该方法不需要计算背景场, 可适用于任意复杂地电结构模型和任意场源形式.

由于电磁场的水平分量连续, 法向分量不连续, 导致直接采用节点有限元会出现伪解.一般有三种处理方法.一是基于矢量位和标量位的间接计算方法(Haber et al., 2000; Badea et al., 2001; Mitsuhata and Uchida, 2004), 由于矢量位和标量位在求解区域连续, 符合节点有限元的前提条件, 但该方法通过位量数值求导计算电场和磁场, 导致误差增大; 另一种方法是基于棱边的矢量有限元(Nédélec, 1986; Hu et al., 2015), 该方法由于直接将自由度赋予单元的棱边上, 避开了电场法向不连续性问题, 且电磁场散度为零自动满足, 近年来受到了极大关注.第三种方法是采用传统节点有限元方法加散度校正(金建铭, 1998), 强制电磁场的散度为零, 该方法沿用传统的节点有限元, 不用修改形函数, 计算效率高.

频域可控源电磁有限元正演最后都会形成一个大型稀疏的复系数方程组, 如何快速求解该方程组成为提高计算速度的关键.对于多个发射源问题, 会形成多个右端项, 常采用直接方法, 主要有基于LU矩阵分解法和Open MP的Pardiso求解器(Schenk and Gartner, 2004), 基于多波前法和Mumps求解器(Da Silva et al., 2012), 只需要一次因式分解就可以获得多个场源的解.另外一种方法就是Krylov子空间迭代法, 该方法和直接法相比, 占用内存小, 计算速度快, 成为非常受欢迎的一种求解大型稀疏线性方程组方法(Börner et al., 2015).此外, 有限元-无限元相结合(张林成等, 2017)、多重网格技术(Mulder, 2008)和区域分解法(Zyserman, 2000; Ren et al., 2014)也成为减少网格和计算量, 实现快速正演的一种方法.

对于陆地多频点可控源电磁法三维正演, 采用基于模型降阶的Krylov子空间投影方法更具有吸引力, 它通过构建传递函数(蒋耀林, 2010), 一次求解即可得到一定范围内所有频率的数值解, 同时保证较高精度的计算结果.模型降阶方法最早应用于地球物理瞬变电磁场的计算(Börner et al., 2008), 通过选择合适的参考频率, 用Krylov子空间投影作为模型降阶技术(Antoulas, 2005), 产生一个传递函数的理想近似, 构建一个Krylov子空间正交基, 计算完成所有频率域解后, 用汉克尔变换得到瞬变电磁场(Newman et al., 1986), 该方法进一步扩展到直接在时间域计算瞬变场(Börner et al., 2015).周建美等(2018)采用基于Yee氏交错网格的拟态有限体积算法, 把m维单重复极点的有理函数Krylov子空间模型降阶方法应用于方程组求解, 大大提高了效率.

本文基于Krylov子空间技术, 采用单极点模型降阶算法, 把大型有限元刚度矩阵降阶为小维度矩阵, 采用直接求逆的方法求解小规模方程.模型降阶前后的矩阵与发射频率无关, 频率不再是影响计算速度的主要因素.通过一维解析解验证了程序的正确性, 给出了与3DFEM有限元计算精度和效率的比较结果.最后, 通过多个三维地电模型的正演计算, 说明程序的通用性和适用性.

1 有限元分析 1.1 控制方程

在均匀导电介质中, 当传导电流占支配地位时, 即σ≫ωε时, 称为准静态极限, 通常可以忽略位移电流∂D/∂t的作用.由于在通常情况下, 大地电阻率不大于1000 Ωm, 电磁波频率不高于100 kHz时, 可视为准静态极限(汤井田和何继善, 2005).在CSAMT中, 假设谐变因子为e－iωt, 那么在准静态极限下的电场E和磁场H所满足的频率域方程为(Cai et al., 2017)

(1)

(2)

式(1)中μ为磁导率, 本研究不考虑磁导率的变化, 故取真空磁导率μ0; ω为圆频率, f为频率, 式(2)中Js为传导电流密度, σ为电导率.

通过式(1)和式(2)可导出电场的双旋度方程为

(3)

其中, k为准静态极限下的波数, k2=iωμσ.

对式(3)采用广义变分原理, 可得到其变分方程为(金建铭, 1998)：

(4)

利用以上方程求解节点有限元矢量场时, 有时会得到伪解.这种解不满足散度条件, 主要由于在有限元数值模拟中, 离散计算区域只要求插值函数连续, 而未对其导数加以限制, 因此这种条件下解得的电场有时是错误的.为了减弱这种现象, 在原变分方程中加入一罚项强加散度条件(汤文武等, 2018)：

(5)

s为罚因子, 衡量罚项所占的比重, 通常取s=1.

1.2 有限单元分析

本文采用节点有限元方法, 六面体单元剖分整个区域, 将六面体单元经过等参变换为母单元, 母单元是一个边长为2的正方体, 取8个角点为节点, 母单元和子单元的编号及坐标见图 1, 两者之间的坐标关系为(徐世浙, 1994)

图 1 六面体单元 (a)母单元; (b)子单元. Fig. 1 Hexahedral element (a) Parent element; (b) Subelement.

(6)

其中x0, y0, z0为子单元的中点, a、b、c为子单元的三个方向的边长, 两个单元的微分关系为：

(7)

构造如下形函数：

(8)

ξi, ηi, ζi为母单元中节点, i(i=1, …, 8)

(9)

通过如上构造的形函数, 任意单元e中各个节点上的电场分量大小可表示为

(10)

若将区域剖分为n个单元, 离散化后的变分方程为

(11)

求该变分方程的极值问题, 相当于对多元函数F(Ex, Ey, Ez)取极值, 多元函数取极值应当满足以下条件(徐世浙, 1994)：

(12)

其中, i=1, 2, …, Nd, Nd为节点总数, n为单元总数.

经过单元分析, 离散方程组中与某一单元相关的内容可以表示为

(13)

其中; be为右端源项, 采用伪δ函数以消除源的奇异性.

按照单元节点的总体序号, 将单元系数矩阵K1e和K2e的各个元素分别放在总体刚度矩阵C和M的相应行与列的交叉位置上, 扩展成3Nd×3Nd的大型稀疏矩阵(徐世浙, 1994).同理将Ee和be扩展成成3Nd×1的矢量矩阵E和b.

(14)

最终离散化变分方程求极值问题转变为求解方程组式(14), 该方程组为大型稀疏复系数方程组.

电偶极子源在源点处表现出奇异性, 一种消除奇异性的较好的办法是采用伪δ函数, 将源点分布到一个小的区域, 消除源点奇异性(Mitsuhata, 2000).当场源为沿x方向的电偶源时, Js为

(15)

其中(x0, y0, z0)为场源坐标, I为电流, ds为电偶极子长度.

2 模型降阶算法 2.1 传递函数的构建

由最终的离散方程组式(14), 令A=M-1C, X=M-1b, 则式(14)变为

(16)

其中I为单位阵, 则E可看成ω的函数：

(17)

构建频率域的传递函数gω(A)(蒋耀林, 2010)：

(18)

其中, gω(A)=-iω(A+iωI)-1.

求解E(ω)的关键就是求解其传递函数gω(A).

在所有的迭代方法中, Krylov子空间迭代方法被认为是求解超大型线性方程组最有效的方法(张继锋, 2008).由于在CSAMT三维正演中, 右端源项是不随频率改变的, X=M-1b, 因此我们采用Krylov子空间投影算法, 在满足精度要求的前提下, 可大大提高运算效率.

2.2 模型降阶过程

传递函数gω(A)是一个大型稀疏的矩阵, Krylov子空间方法可用于构造此类矩阵的有效逼近rm(A), 使得

(19)

其中, rm∈qm－1pm, qm(z)=(z－ξ1)(z－ξ2)…(z－ξm), ξi(i=1, …, m)为m个极点, pm为m阶非零多项式.

因此求解rm(A)等价于求解A在m维Krylov子空间Qm+1=qm(A)－1κm+1(A, X)上的投影(周建美等, 2018).由于本文所得到的系数矩阵M并非单位阵, 若直接对式(18)的形式采用传统的Arnoldi算法来求解正交化矩阵会使计算变得复杂, 因此可以对Krylov子空间采用一种M正交Arnoldi算法, 即M正交运算, 进而得到标准列正交向量Vm+1(蒋耀林, 2010; Börner et al., 2015).最终的表达形式为

(20)

其中, =.

Am+1是A对于Vm+1的正交投影, 在Krylov子空间Qm+1上产生的m+1维的方阵.

从式(20)可看出, 模型降阶算法所构建的矩阵Am+1的维度远远小于原来的维度, 因此可以直接求逆, 由于模型降阶前后的矩阵均与源的频率无关, 故只要一次模型降阶便可实现多个频点的运算, 在Börner原有的研究成果可看出, 模型降阶的效率取决于多个不同极点所产生方程的求解.图 2给出了模型降阶的主要步骤, 其中, Lp为第p次迭代的Cholesky分解下三角阵, vj+1为正交基向量, w是中间向量.

图 2 有限元模型降阶算法流程图 Fig. 2 Flow chart of model reduction for finite element method

从正交基构建的过程可看出, 运算的速率取决于极点的个数, 而周建美等(2018)的研究表明, 当τ=iω为正虚数时, 存在一个最优化实数极点ξ0=, 使得在足够多的迭代次数下, 求解误差能快速收敛并趋于定值.对方程(14)进行单极点模型降阶, 只要使用稀疏矩阵求解器pardiso进行一次矩阵分解, m次矩阵回代即可完成正交基的构建, 因此通过这种重复极点的选取方法, 能够在满足精度条件下快速求解多个频点的正演问题.

在选定单个重复极点后, 可利用该极点进行C－ξpM的矩阵构建, 将单元矩阵组装为整体刚度矩阵, 由于矩阵是一个稀疏对称矩阵, 在整体刚度矩阵集成之前, 非零元素的个数和位置是未知的, 若直接存储将会占用很大的静态空间, 造成存储空间的浪费.采用CSR压缩存储非零元素可大大节省内存, 提高运算效率, 本次模拟采用动态的双向链表存储上三角部分, 利用MKL库函数Pardiso进行高精度求解.

通过模型降阶方法进行快速CSAMT正演, 可得到各个网格节点的电场三分量, 便可通过求导计算得到磁场Hy和卡尼亚视电阻率:

(21)

(22)

3 数值算例 3.1 程序可靠性检验

为检验程序的可靠性和适用性, 用含有解析解的模型进行误差分析, 本文采用三维有限元程序3DFEM和模型降阶算法分别计算了二层地电模型的水平电场和视电阻率, 方程求解器均采用IntelⓇ Math Kernel Library (IntelⓇ MKL)中的大型稀疏方程组求解器Pardiso, 比较了精度和运算速率.采用陈辉等(2016)的三维模拟结果对比分析, 验证了模型降阶程序的正确性.本次数值模拟使用普通的笔记本电脑, 处理器为Inter(R) Core(TM) i7-9750H CPU @ 2.60GHz.

3.1.1 二层模型

构建二层模型, 将坐标原点置于电性源中心, XYZ坐标轴按右手螺旋放置, Z轴向下为正, 电性源沿X轴放置, 发射电流1 A, 频点个数为12个, 最低频率2 Hz, 最高频率4096 Hz, 测点坐标为(-1500, -2400, 0).第一层地层电阻率500 Ωm, 厚度200 m, 第二层电阻率为50 Ωm, 选取单个重复极点; 将研究区域剖分为51×47×57的网格, 网格剖分边界为±20×104 m, 以保证边界场值为零, 在测点、场源处加密网格, 以减小场源的奇异性, 降低因等距插值而带入的误差.

图 3a为电场水平分量的对比图, 可看出拟合程度较好, 图 3b为相对误差图, 3DFEM与理论值相比, 高频和低频误差较大, 但大部分点相对误差小于5%, 模型降阶解同解析解间的平均相对误差为1.72%.图 4为卡尼亚视电阻率和相对误差对比图, 3DFEM同理论值之间的误差为2.4%, 模型降阶解同解析解之间的相对误差为2.1%, 误差有所增大, 主要是因为求解磁场需要进行数值求导, 会引进数值计算误差, 导致卡尼亚电阻率相对误差变大.值得一提的是模型降阶方法在高频段计算比较稳定, 误差要远小于3DFEM.

图 3 二层模型降阶解、3DFEM与解析解的电场水平分量对比 (a)电场水平分量; (b)相对误差. Fig. 3 Comparison of horizontal electric field of model reduction, 3DFEM and analytical solutions for two layers medium (a) Horizontal electric field; (b) Relative error.
图 4 二层模型降阶解、3DFEM与解析解的卡尼亚视电阻率对比 (a)卡尼亚视电阻率; (b)相对误差. Fig. 4 Comparison of Cargniard apparent resistivity of model order reduction, 3DFEM and analytical solution for two layers (a) Cargniard apparent resistivity; (b) Relative error.

该模型的运算效率对比如表 1所示, 在同一计算条件下, 3DFEM的运算时间大于模型降阶解的10倍以上.

表 1 运算速率对比表 Table 1 Comparison of calculation time

本文采用二层模型验证模型降阶程序的正确性, 可以看出, 模型降阶解的结果同3DFEM解的结果具有一致性, 甚至计算精度比3DFEM求解更高, 这主要是由于模型降阶解所带入的数值计算误差比直接求解大型稀疏矩阵的3DFEM所带入的误差角度看, 在误差允许范围内模型降阶求解大型稀疏矩阵是可行的.从运算速率角度考虑, 模型降阶只需要求解一次包含大型稀疏矩阵的方程, 多次矩阵回代即可构建所需的正交阵, 因此多个频点的求解变得快速, 而使用3DFEM直接求解时, 每个频点都需要求解一次包含大型稀疏矩阵的方程, 因此在运算速率方面, 远逊色于模型降阶算法.而基于频点并行的3DFEM求解, 只是将每个频点分配给不同的线程处理, 在CSAMT陆地三维正演求解中, 需要计算较多频点, 3DFEM并行也仅仅只能在增加线程的数量上来提高运算速率, 而在模型降阶中频点的个数对运算时间的影响并不是主要的, 因此将模型降阶运用到CSAMT陆地三维正演中是可行的.

3.1.2 三维模型

如图 5, 构建一个三维地质模型, 低阻异常体模型大小为800 m×800 m×500 m, 顶部埋深为500 m, 异常体电阻率为10 Ωm, 围岩电阻率为100 Ωm.计算频率为1~1000 Hz, 接地导线长为2000 m, 采用赤道模式观测, 收发距为6000 m.

图 5 低阻模型 Fig. 5 Model of conductive body

图 6a和图 6b分别表示模型降阶结果同陈辉等文章模拟结果的对比, 图 6a表示卡尼亚视电阻率的对比, 图 6b表示相位的对比, 从中可以看出, 二者的变化趋势以及数值大小都吻合得较好, 因此可以进一步证明模型降阶方法在CSAMT陆地三维正演中的可行性.

图 6 低阻体正演结果对比 (a)卡尼亚视电阻率; (b)相位. Fig. 6 Comparison of the results for conductive body (a) Cagniard apparent resistivity; (b) Phase.
3.2 水平方向高阻-低阻模型

建立如图 7所示的水平高阻-低阻三维模型, 图 7a为平面, 图 7b为剖面图.电性源沿y轴放置, 发射电流1A, 源的中心在坐标轴原点.

图 7 水平方向高阻-低阻模型 (a)平面图; (b)剖面图. Fig. 7 Model of resistive and conductive bodies in horizontal direction (a) Plane view; (b) Cross-section.

我们分别选取高频4096 Hz和低频256 Hz, 绘制其平面图, 如图 8.在高频4096 Hz的平面图中, 异常电阻率与模型电阻率刚好相反, 这个就是电磁法勘探中常见的“假极值”问题; 随着频率的降低, 探测深度增加, 因而能看到正确反映目标体电阻率大小和目标体范围的视电阻率异常, 低阻体异常范围较高阻体大, 这与场源的“阴影效应”使低阻体产生更大的”阴影”吻合.在较低频率, 由于高阻体一侧距离场源更近, 因而高阻体频率先进入过渡带(如图 10中32 Hz、16 Hz切片), 随着频率的降低, 更远的地层也相继进入过渡带(如图 10中8 Hz、4 Hz切片).图 9是两个频率所对应的视电阻率和电场分量, “假极值”现象也可看到, 而且低阻体异常幅值远大于高阻体异常幅值, 可控源电磁法对低阻体的灵敏度要远大于高阻体.

图 8 (a) 4096 Hz切片图; (b) 256 Hz切片图 Fig. 8 (a) 4096 Hz slice; (b) 256 Hz slice
图 9 2700号测线在4096、256 Hz的卡尼亚视电阻率和电场水平分量对比 (a)卡尼亚视电阻率; (b)电场水平分量. Fig. 9 Comparison of Cargniard apparent resistivity and horizontal electric field for 4096 Hz, 256 Hz at the survey line 2700 (a) Cargniard apparent resistivity; (b) Horizontal electric field.
图 10 多频点卡尼亚视电阻率三维切片图 Fig. 10 Three-dimensional slice chart for multi-frequency Cagniard apparent resistivity
3.3 垂直方向高阻-低阻模型

设计如图 11所示的垂向高阻低阻模型, 图 11a为平面图, 图 11b为剖面图, 电性源沿y轴放置, 发射电流1A, 源的中心在坐标轴原点.

图 11 垂直方向高阻-低阻模型 (a)平面图; (b)剖面图. Fig. 11 Model of vertical resistive and conductive model (a) Plane view; (b) Cross-section.

图 12为卡尼亚电阻率纵向切片图, 展示了不同测线剖面的视电阻率变化, 横坐标测线上对应的测点的x坐标, 纵坐标为视深度, 在其中的测线2700、3000、3300、3600均能明显看到异常体的存在, 且低阻体产生明显的“阴影效应”, 异常范围远大于异常体本身.在中间几条主测线上, 可以看到上下高阻体和低阻体响应, 但距离异常体越远, 异常越弱, 特别是高阻体, 在旁测线几乎看不到异常响应.

图 12 卡尼亚视电阻率纵向切片图 Fig. 12 Vertical slice chart of Cargniard apparent resistivity
3.4 低阻陷落柱

为了验证模型降阶程序对复杂模型的适用性, 本文构建了一个低阻陷落柱模型.模型图如图 13所示, 该异常体中心坐标为(7000 m, 6000 m, 215 m).电性源沿y轴放置, 发射电流1A, 源的中心在坐标轴原点.

图 13 陷落柱模型示图 (a)平面图; (b)剖面图. Fig. 13 Conductive collapse column model (a) Plane view; (b) Cross-section.

图 14为中心测线剖面图, 由于异常体规模较小, 上覆地层厚, 在剖面图上并未出现圈闭, 但是可以明显看到异常体中心位置所对应的测点7000 m附近视电阻率出现了明显的下凹, 与模型位置吻合, 证明了低阻陷落柱的存在.

图 14 中心测线视电阻率断面图 Fig. 14 Apparent resistivity section for the center survey line

图 15分别展示了1024 Hz, 256 Hz, 128 Hz的横向切片图, 从中可以看出, 1024 Hz低阻体位置处呈现高阻, 而随着频率的降低, 勘探深度增加, 低阻异常体中心呈现与之对应的低阻异常, 这证明在高频段出现了“假极值”现象, 但其值相对于异常幅值比较小.

图 15 (a) 1024 Hz切片图; (b) 256 Hz切片图; 128 Hz切片图 Fig. 15 (a) 1024 Hz slice; (b) 256 Hz slice; (c) 128 Hz slice
4 结论

本文采用基于Krylov子空间投影算法的模型降阶方法, 采用单重复极点迭代, 实现了多频发射的陆地三维可控源电磁法的快速正演.推导了三维可控源电磁法有限元变分方程, 施加了散度条件, 消除了伪解; 采用基于总场的有限元方法, 通过伪δ函数, 减小了源的奇异性.

以层状地层模型解析解为标准, 比较了3DFEM和模型降阶的有限元算法, 结果表明：在不降低精度条件下, 模型降阶方法求解时间不到3DFEM的1/10, 极大地提高了计算效率, 为进一步精确反演提供了条件.

分别对水平方向高阻-低阻模型, 垂直方向高阻-低阻模型以及复杂地层低阻陷落柱模型进行了正演计算, 异常体的空间位置与卡尼亚视电阻率断面图结果吻合, 进一步证明该算法的正确性和可实用性; 分析了由源所产生的阴影效应, “假极值”现象和过渡带、近区响应的特征, 在实际勘探中对这些现象和特征要仔细分析辨识, 才能更好地进行处理解释.

基于模型降阶的快速正演能够大幅度提高计算效率, 频点个数不再是影响计算时间的主要因素, 在本质上改变了求解方式, 因而该算法可运用于三维陆地多频点CSAMT正演模拟, 为进一步实现三维快速反演奠定基础.

致谢  感谢MKL高性能计算库提供的Pardiso求解器, 为本文的完成提供了方便.
References
Antoulas A C. 2005. Approximation of Large-Scale Dynamical Systems. Philadelphia, PA: SIAM Publications.
Avdeev D B. 2015. Three-dimensional electromagnetic modelling and inversion from theory to application. Surveys in Geophysics, 26(6): 767-799.
Badea E A, Everett M E, Newman G A, et al. 2001. Finite-element analysis of controlled-source electromagnetic induction using Coulomb-gauged potentials. Geophysics, 66(3): 786-799. DOI:10.1190/1.1444968
Börner R U, Ernst O G, Spitzer K. 2008. Fast 3-D simulation of transient electromagnetic fields by model reduction in the frequency domain using Krylov subspace projection. Geophysical Journal International, 173(3): 766-780. DOI:10.1111/j.1365-246X.2008.03750.x
Börner R U. 2010. Numerical modeling in geo-electromagnetics:Advances and challenges. Surveys in Geophysics, 31(2): 225-245. DOI:10.1007/s10712-009-9087-x
Börner R U, Ernst O G, Güttel S. 2015. Three-dimensional transient electromagnetic modelling using Rational Krylov methods. Geophysical Journal International, 202(3): 2025-2043. DOI:10.1093/gji/ggv224
Cai H Z, Xiong B, Han M R, et al. 2014. 3D controlled-source electromagnetic modeling in anisotropic medium using edge-based finite element method. Computers & Geosciences, 73: 164-176.
Cai H Z, Xiong B, Zhdnov M. 2015. Three-dimensional marine controlled-source electromagnetic modelling in anisotropic medium using finite element method. Chinese Journal of Geophysics (in Chinese), 58(8): 2839-2850. DOI:10.6038/cjg20150818
Cai H Z, Hu X Y, Li J H, et al. 2017. Parallelized 3D CSEM modeling using edge-based finite element with total field formulation and unstructured mesh. Computers & Geosciences, 99: 125-134.
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 (in Chinese), 59(8): 3087-3097. DOI:10.6038/cjg20160831
Constable S. 2010. Ten years of marine CSEM for hydrocarbon exploration. Geophysics, 75(5): 75A67-75A81. DOI:10.1190/1.3483451
Da Silva N V, Morgan J V, MacGregor L, et al. 2012. A finite element multifrontal method for 3D CSEM modeling in the frequency domain. Geophysics, 77(2): E101-E115.
Everett M E. 2011. Theoretical developments in electromagnetic induction geophysics with selected applications in the near surface. Surveys in Geophysics, 33(1): 29-63.
Haber E, Ascher U M, Aruliah D A, et al. 2000. Fast simulation of 3D electromagnetic problems using potentials. Journal of Computational Physics, 163(1): 150-171.
Hou J S, Mallan R K, Torres-Verdín C. 2006. Finite-difference simulation of borehole EM measurements in 3D anisotropic media using coupled scalar-vector potentials. Geophysics, 71(5): G225-G233. DOI:10.1190/1.2245467
Hu Y C, Li T L, Fan C S, et al. 2015. Three-dimensional tensor controlled-source electromagnetic modeling based on the vector finite-element method. Applied Geophysics, 12(1): 35-46.
Jiang Y L. 2010. Model Reduction Method (in Chinese). Beijing: Science Press.
Jin J M. 1998. Fnite Element Method of Electromagnetic Field (in Chinese). Wang J G Trans. Xi'an: Xi'an Electronic University Press.
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 (in Chinese), 59(4): 1521-1534. DOI:10.6038/cjg20160432
Liu Y, Li Y G, Han B. 2017. Adaptive edge finite element modeling of the 3D CSEM field on unstructured grids. Chinese Journal of Geophysics (in Chinese), 60(12): 4874-4886. DOI:10.6038/cjg20171227
Mitsuhata Y. 2000. 2-D electromagnetic modeling by finite-element method with a dipole source and topography. Geophysics, 65(2): 465-475. DOI:10.1190/1.1444740
Mitsuhata Y, Uchida T. 2004. 3D magnetotelluric modeling using the T-Ω finite-method. Geophysics, 69(1): 108-119.
Mulder W A. 2008. Geophysical modelling of 3D electromagnetic diffusion with multigrid. Computing and Visualization in Science, 11(3): 129-138. DOI:10.1007/s00791-007-0064-y
Nédélec J. 1986. A new family of mixed finite elements in R3. Numerische Mathematik, 50(1): 57-81. DOI:10.1007/BF01389668
Newman G A, Alumbaugh D. 2002. Three-dimensional induction logging problems, Part 2:A finite-difference solution. Geophysics, 67(2): 484-491. DOI:10.1190/1.1468608
Pan K J, He D D, Hu H L, et al. 2017. A new extrapolation cascadic multigrid method for three dimensional elliptic boundary value problems. Journal of Computational Physics, 344: 499-515. DOI:10.1016/j.jcp.2017.04.069
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 (in Chinese), 59(10): 3927-3939. DOI:10.6038/cjg20161036
Peng R H, Hu X Y, Li J H, et al. 2018. 3-D finite-volume forward modeling of wide-field EM using scattered potentials. Chinese Journal of Geophysics (in Chinese), 61(10): 4160-4170. DOI:10.6038/cjg2018L0363
Puzyrev V, Koldan J, De la Puente J, et al. 2013. A parallel finite-element method for three-dimensional controlled-source electromagnetic forward modelling. Geophysical Journal International, 193(2): 678-693. DOI:10.1093/gji/ggt027
Qin C. 2015. Study of three dimensional staggered-grid finite difference modeling and source effect of CSAMT[Master's thesis] (in Chinese). Chengdu: Chengdu University of Technology.
Qiu C K, Yin C C, Liu Y H, et al. 2018. 3D forward modeling of controlled-source audio-frequency magnetotellurics in arbitrarily anisotropic media. Chinese Journal of Geophysics (in Chinese), 61(8): 3488-3498. DOI:10.6038/cjg2018L0326
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, Chen C J, Tang J T, et al. 2017. A new integral equation approach for 3D magnetotelluric modeling. Chinese Journal of Geophysics (in Chinese), 60(11): 4506-4515. DOI:10.6038/cjg20171134
Ruhe A. 1994. Rational krylov algorithms for nonsymmetric eigenvalue problems.//Golub G, Luskin M, Greenbaum A eds. Recent Advances in Iterative Methods. New York, NY: Springer, 149-164.
Schenk O, Görtner K. 2004. Solving unsymmetric sparse systems of linear equations with PARDISO. Future Generation Computer Systems, 20(3): 475-487.
Streich R. 2009. 3D finite-difference frequency-domain modeling of controlled-source electromagnetic data:direct solution and optimization for high accuracy. Geophysics, 74(5): F95-F105. DOI:10.1190/1.3196241
Tang J T, He J S. 2005. Controlled-source Audio Magnetotelluric Method and Its Application (in Chinese). Changsha: Central South University Press.
Tang J T, Zhou F, Ren Z Y, et al. 2018. Three-dimensional forward modeling of the controlled-source electromagnetic problem based on the integral equation method with an unstructured grid. Chinese Journal of Geophysics (in Chinese), 61(4): 1549-1562. DOI:10.6038/cjg2018L0303
Tang W W, Liu J X, Ye Y X, et al. 2018. Comparison of 3D controlled-source electromagnetic forward modeling based on the nodal finite element and the edge-based finite element. Oil Geophysical Prospecting (in Chinese), 53(3): 617-624.
Um E S, Harris J M, Alumbaugh D L. 2010. 3D time-domain simulation of electromagnetic diffusion phenomena:a finite-element electric-field approach. Geophysics, 75(4): F115-F126. DOI:10.1190/1.3473694
Um E S, Kim S S, Fu H H. 2017. A tetrahedral mesh generation approach for 3D marine controlled-source electromagnetic modeling. Computers & Geosciences, 100: 1-9.
Wang Y L, Di Q Y, Wang R. 2017. Three-dimensional modeling of controlled-source audio-frequency magnetotellurics using the finite element method on an unstructured grid. Chinese Journal of Geophysics (in Chinese), 60(3): 1158-1167. DOI:10.6038/cjg20170326
Weiss C J. 2013. Project APhiD:A Lorenz-gauged A-Φ decomposition for parallelized computation of ultra-broadband electromagnetic induction in a fully heterogeneous Earth. Computers & Geosciences, 58: 40-52.
Xu S Z. 1994. The Finite Element Method in Geophysics (in Chinese). Beijing: Science Press.
Yan S, Xue G Q, Qiu W Z, et al. 2017. Interpretation of CSAMT single-component data. Chinese Journal of Geophysics (in Chinese), 60(1): 349-359. DOI:10.6038/cjg20170129
Yin C C, Zhang B, Liu Y H, et al. 2017. A goal-oriented adaptive algorithm for 3D magnetotelluric forward modeling. Chinese Journal of Geophysics (in Chinese), 60(1): 327-336. DOI:10.6038/cjg20170127
Zhang J F. 2008. Three dimensional controlled source electromagnetic numerical simulation based on electric field double curl equation using finite element method[Ph. D thesis] (in Chinese). Changsha: Central South University.
Zhang J F, Tang J T, Yu Y, et al. 2009. Three dimensional controlled source electromagnetic numerical simulation based on electric field vector wave equation using finite element method. Chinese Journal of Geophysics (in Chinese), 52(12): 3132-3141. DOI:10.3969/j.issn.0001-5733.2009.12.023
Zhang J F, Zhi Q Q, Li X, et al. 2013. 2.5D finite element numerical simulation for electric dipole source on ridge terrain. The Chinese Journal of Nonferrous Metals (in Chinese), 23(9): 2540-2550.
Zhang L C, Tang J T, Ren Z Y, et al. 2017. Forward modeling of 3D CSEM with the coupled finite-infinite element method based on the second field. Chinese Journal of Geophysics (in Chinese), 60(9): 3655-3666. DOI:10.6038/cjg20170929
Zhdnov 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. DOI:10.1190/1.2358403
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 J. Geophys. (in Chinese), 55(6): 2141-2152. DOI:10.6038/j.issn.0001-5733.2012.06.036
Zhou J M, Liu W T, Liu H, et al. 2018. Research on rational Krylov subspace model order reduction algorithm for three-dimensional multi-frequency CSEM modelling. Chinese Journal of Geophysics (in Chinese), 61(6): 2525-2536. DOI:10.6038/cjg2018L0311
Zyserman F I, Santos J E. 2000. Parallel finite element algorithm with domain decomposition for three-dimensional magnetotelluric modelling. Journal of Applied Geophysics, 44(4): 337-351. DOI:10.1016/S0926-9851(00)00012-4
蔡红柱, 熊彬, Zhdnov M. 2015. 电导率各向异性的海洋电磁三维有限单元法正演. 地球物理学报, 58(8): 2839-2850. DOI:10.6038/cjg20150818
陈辉, 殷长春, 邓居智. 2016. 基于Lorenz规范条件下磁矢势和标势耦合方程的频率域电磁法三维正演. 地球物理学报, 59(8): 3087-3097. DOI:10.6038/cjg20160831
蒋耀林. 2010. 模型降阶方法. 北京: 科学出版社.
金建铭. 1998.电磁场有限元方法.王建国译.西安: 西安电子科技大学出版社.
李建慧, Farquharson C G, 胡祥云, 等. 2016. 基于电场总场矢量有限元法的接地长导线源三维正演. 地球物理学报, 59(4): 1521-1534. DOI:10.6038/cjg20160432
刘颖, 李予国, 韩波. 2017. 可控源电磁场三维自适应矢量有限元正演模拟. 地球物理学报, 60(12): 4874-4886. DOI:10.6038/cjg20171227
彭荣华, 胡祥云, 韩波, 等. 2016. 基于拟态有限体积法的频率域可控源三维正演计算. 地球物理学报, 59(10): 3927-3939. DOI:10.6038/cjg20161036
彭荣华, 胡祥云, 李建慧, 等. 2018. 基于二次耦合势的广域电磁法有限体积三维正演计算. 地球物理学报, 61(10): 4160-4170. DOI:10.6038/cjg2018L0363
秦策. 2015.可控源音频大地电磁法三维有限差分正演模拟及场源效应研究[硕士论文].成都: 成都理工大学.
邱长凯, 殷长春, 刘云鹤, 等. 2018. 任意各向异性介质中三维可控源音频大地电磁正演模拟. 地球物理学报, 61(8): 3488-3498. DOI:10.6038/cjg2018L0326
任政勇, 陈超健, 汤井田, 等. 2017. 一种新的三维大地电磁积分方程正演方法. 地球物理学报, 60(11): 4506-4515. DOI:10.6038/cjg20171134
汤井田, 何继善. 2005. 可控源音频大地电磁法及其应用. 长沙: 中南大学出版社.
汤井田, 周峰, 任政勇, 等. 2018. 复杂地下异常体的可控源电磁法积分方程正演. 地球物理学报, 61(4): 1549-1562. DOI:10.6038/cjg2018L0303
汤文武, 柳建新, 叶益信, 等. 2018. 基于节点有限元与矢量有限元的可控源电磁三维正演对比. 石油地球物理勘探, 53(3): 192-199.
王亚璐, 底青云, 王若. 2017. 三维CSAMT法非结构化网格有限元数值模拟. 地球物理学报, 60(3): 1158-1167. DOI:10.6038/cjg20170326
徐世浙. 1994. 地球物理中的有限单元法. 北京: 科学出版社.
闫述, 薛国强, 邱卫忠, 等. 2017. CSAMT单分量数据解释方法. 地球物理学报, 60(1): 349-359. DOI:10.6038/cjg20170129
殷长春, 张博, 刘云鹤, 等. 2017. 面向目标自适应三维大地电磁正演模拟. 地球物理学报, 60(1): 327-336. DOI:10.6038/cjg20170127
张继锋. 2008.基于电场双旋度方程的三维可控源电磁法有限单元法数值模拟[博士论文].长沙: 中南大学.
张继锋, 汤井田, 喻言, 等. 2009. 基于电场矢量波动方程的三维可控源电磁法有限单元法数值模拟. 地球物理学报, 52(12): 3132-3141. DOI:10.3969/j.issn.0001-5733.2009.12.023
张继锋, 智庆全, 李貅, 等. 2013. 起伏地形电偶源2.5维有限元数值模拟. 中国有色金属学报, 23(9): 2540-2550.
张林成, 汤井田, 任政勇, 等. 2017. 基于二次场的可控源电磁法三维有限元-无限元数值模拟. 地球物理学报, 60(9): 3655-3666. DOI:10.6038/cjg20170929
张烨, 汪宏年, 陶宏根, 等. 2012. 基于耦合标势与矢势的有限体积法模拟非均匀各向异性地层中多分量感应测井三维响应. 地球物理学报, 55(6): 2141-2152. DOI:10.6038/j.issn.0001-5733.2012.06.036
周建美, 刘文韬, 刘航, 等. 2018. 多频可控源电磁法三维有理函数Krylov子空间模型降阶正演算法研究. 地球物理学报, 61(6): 2525-2536. DOI:10.6038/cjg2018L0311