海洋可控源电磁法(Marine Controlled-Source Electromagnetic,简称MCSEM)是近年来快速发展起来的一种应用于海洋油气勘探的地球物理方法.自2000年挪威国家石油公司Statoil利用频率域海洋可控源电磁在Angola近海实现SBL(Sea Bed Logging)海底勘探以来,MCSEM作为一种重要的钻前储层评价手段可减小勘探风险,正逐渐广泛应用于海洋油气勘探领域中(Ellingsrud等,2002).MCSEM的发射机采用近海底拖曳的水平电偶极子源,向海水中发送0.1~10 Hz的大功率电磁脉冲,接收机沿测线阵列布设于海底,通过综合分析不同位置观测的电场和磁场信号,获得海底三维电性分布特征.随着工业勘探不断广泛的开展,众多学者对MCSEM正反演理论展开了研究(Constable,2010;Key,2012),精确有效的三维MCSEM正反演成为当今电磁勘探领域的研究热点.在过去的二十几年里,MCSEM数值模拟取得了巨大的进步.Zhdanov(2006)利用积分方程法对非均匀背景介质中含三维异常体模型进行正演模拟,取得了非常好的效果;陈桂波等(2009)进一步用积分方程法模拟了各向异性介质中三维海洋可控源电磁响应,分析了各向异性对MCSEM响应的影响特征.Li(2007)利用基于非结构网格的自适应有限元算法实现了复杂地形情况下二维海洋可控源正演模拟;杨军等(2015)利用矢量有限元法实现了非结构单元剖分海洋可控源电磁三维数值模拟.殷长春等(2014)利用有限差分法实现三维任意各向异性介质海洋可控源电磁正演模拟.韩波等(2015)利用交错网格有限体积法实现了基于并行化直接解法的频率域海洋可控源电磁三维正演.现有的三维MCSEM正演方法主要基于积分方程法、有限元法、有限差分法和有限体积法,这些数值方法各有优劣,并在应用中得到不断改进和完善.Avdeev(2005)和Brner(2010)对这些方法的特点和现状进行了较为全面的评述.然而,如何获得高精度和快速三维正演模拟结果目前仍然是海洋可控源电磁的一个技术难点,在一定程度上制约了MCSEM技术的应用和发展.因此,本文尝试将一种新型的数值计算方法—谱元法引入到海洋可控源电磁三维正演模拟中.
谱元法是一种求解偏微分方程的高精度数值方法.它最早由Patera提出并应用于流体力学(Patera,1984).由于谱元法具有很好的精度和收敛特性(Maday等,1987),国内外许多学者对其应用展开了系列研究.Seriani等(1991)与Komatitsch和Tromp(1999)分别将基于Chebyshev和Legendre基函数的谱元法引入到地震波场传播模拟中;王秀明等(2007)引入基于预先条件下共轭梯度法的元到元算法,提高了谱元法的计算效率和精度;Komatitsch等人开发了谱元软件SPECFEM3D_Cartesian,在地震波场模拟中取得了很好的应用效果(Tape等,2010);刘有山(2014)等利用Koornwinder-Dubiner正交多项式和Fekete积分节点,实现了三角网格剖分下的地震波场谱元模拟.李琳等(2014)采用Cohen节点及对应的Lagrange形函数实现三角谱元法在地震正演模拟中的应用.谱元法也是计算电磁中应用较多的技术,Lee等(2006)将基于Gauss-Lobatto-Legendre(GLL)多项式的混合阶谱元法用于求解三维矢量电磁波动方程和三维瞬变电磁问题中,相对传统高阶有限元显著减少了CPU耗时;Ren等(2015)发展了时域间断伽辽金谱元法求解一阶麦克斯韦方程组,使得未知数个数和计算量大大减少.然而,在地球物理电磁正演模拟中,除了Yin and Huang(2017)利用基于GLL基函数的谱元法实现了频率域航空正演模拟外,谱元法在其它领域尚未得到足够的重视和广泛的应用.
与有限元类似,谱元法将模型区域剖分为互不重叠的基本单元,在每个单元内对未知量(如电磁场值)在配置点处用类似谱元法的N阶Chebyshev或Legendre正交多项式进行插值逼近,再借助Galerkin加权残差法形成单元矩阵,进而合成总体大型线性方程组,最后求解得到全局的近似解.因此,谱元法结合了有限元方法良好的模型适应性和谱方法高精度的优点,可通过增加单元数(h型)或提高插值函数阶数(p型)两种方式提高收敛性.单元剖分使得谱元法相对于谱方法可以更好的模拟介质模型;在每个单元内采用正交完备的Chebyshev或Legendre多项式作为基函数,使得随阶数增长正演结果呈指数收敛,因而谱元法较之于传统p型有限元更具优势.在基函数选取上,可以采用Gauss-Lobatto-Legendre(GLL)或Gauss-Lobatto-Chebyshev(GLC)多项式.通过降阶的GLL积分,基于GLL多项式的谱元法能够得到对角的质量矩阵.然而,Chebyshev级数在逼近函数时是一致收敛的,理论上相对于其他多项式(如Legendre正交多项式)应当得到更精确的解,并且使用Chebyshev多项式在计算单元矩阵时可以直接求得解析表达式,无需使用数值积分,因而在方程组求解之前保持很高的精度.此外,Chebyshev节点可以保证在边界处配置更密的节点间隔,有效地避免了传统高阶插值可能产生的Runge现象(Runge,1901).鉴于上述基于GLC基函数谱元法的独特优势,我们将其应用于本文的三维海洋可控源电磁正演模拟中.下面的讨论中,我们由海洋可控源电磁三维频率域控制方程出发,重点介绍谱元法的基本原理和实施步骤,利用一维模型数值结果验证本文方法的准确性,并通过与有限差分方法对比检验方法的有效性.
1 理论从微分形式的麦克斯韦方程组出发,采用正时谐因子ejωt,可得频率域麦克斯韦方程组为
(1) |
(2) |
其中E和H分别为电场和磁场,μ为磁导率,ε为介电常数,σ为电导率,ω为角频率,
为了避免源的奇异性,将电磁场分解成背景场和二次场,则由(1)和(2)式得到频率域下二次电场Es满足的矢量亥姆赫兹方程
(3) |
其中
(4) |
式中,σp和Ep分别为背景电导率和背景电场.我们利用开源程序Dipole1D(Key,2009)计算背景场.
定义函数空间H(curl, Ω)={V∈L2(Ω), ∇×V∈L2(Ω)},其中L2中函数平方可积,内积定义为
(5) |
其中V和W是属于H(curl, Ω)空间的矢量函数.
为了应用谱元法,我们用六面体单元离散整个模型区域Ω.在单元内将电场关于高阶矢量插值基函数
(6) |
其中Nx、Ny、Nz分别为插值基函数沿三个坐标方向的阶数,高阶插值相当于对单元物理棱边进行了再次剖分.图 1给出Nx=3,Ny=3,Nz=3的单元剖分. Nb=Nx(Ny+1)(Nz+1)+(Nx+1)Ny(Nz+1)+(Nx+1)(Ny+1)Nz为单元棱边总数,
对(3)式的亥姆赫兹方程应用伽辽金加权余量法,并运用第一矢量格林定理和Dirichlet边界条件,可得
(7) |
其中,K为剖分的单元总数.将(7)式改写成矩阵形式
(8) |
其中A为刚度矩阵,B为质量矩阵,b为右端源向量,k02=jωμσ-ω2με.
谱元法中的基函数GLL或GLC正交多项式是定义在区间[-1, 1]上的,首先我们需在单元上通过坐标映射将各物理单元(x, y, z)∈[xe, xe+1]×[ye, ye+1]×[ze, ze+1]转换为标准参考单元(ξ, η, ζ)∈[-1, 1]×[-1, 1]×[-1, 1],进而在参考单元内构造插值基函数Φ=[Φξ, Φη, Φζ].物理单元的插值基函数
(9) |
(10) |
其中J为雅可比矩阵,其表达式为
(11) |
由此,参考坐标系下的单元刚度矩阵、单元质量矩阵和单元源向量可表示为
(12) |
(13) |
(14) |
其中 |J|为雅可比行列式,Ejp为第j条单元棱边对应的背景场,Φj为参考单元内第j个插值基函数,Vj为参考单元内第j个基函数的旋度.
(15) |
(16) |
至此,我们通过坐标映射将物理单元下的关于(x, y, z)坐标的离散方程(7)全部转换为参考单元下关于(ξ, η, ζ)坐标的形式.在参考单元内,本文通过三个坐标轴方向GLC多项式的混合张量积构造插值基函数 Φ =[Φξ, Φη, Φζ],保证了切向电场的连续性,其形式如下
(17) |
其中 eξ,eη和eζ表示沿ξ, η, ζ方向的单位向量,hk(N)(ξ)为N阶GLC插值多项式
(18) |
其中k, t=0, 1, 2, …, N,Tt(ξ)=cos(t·arccos(ξ))为t阶Chebyshev多项式,ξk=-cos(kπ/N)为N阶Chebyshev多项式的第k个极值点.
由于Chebyshev多项式有确定的极值点和简单的表达式,当选用GLC多项式作为插值基函数时,(12)—(14)式中的单元积分可以得到解析表达式,在方程组求解之前极大地减小了数值误差(王秀明等,2007).例如,单元质量矩阵包含积分项:
(19) |
其中
(20) |
类似地,我们可以解析计算单元刚度矩阵中的积分项.获得各单元矩阵元素后,通过单元连接关系合成整体系数矩阵和右端源向量.针对发射源随船移动的海洋可控源电磁问题,对于某一频率正演模拟中,仅源项变化而系数矩阵不变,使用直接求解器求解时仅需进行一次矩阵分解,通过不断回带源项实现多源问题快速求解.因此,本文使用并行多波前稀疏矩阵求解器MUMPS(Amestoy等,2001)求解方程(8)获得各棱边电场值,进而插值得到接收点处的二次电场,接收点处的总场值由该二次场加上背景场获得.谱元法的具体计算流程如图 2所示.
为了检验本文谱元算法的准确性,我们建立如图 3所示的海底三层介质模型.海水深度1000 m,电阻率0.3 Ωm,海底围岩电阻率1.0 Ωm,高阻层位于海底1000 m处,电阻率为80 Ωm,厚度为100 m.我们将接收机布设在海底,水平电偶极子发射源随船沿x轴方向拖曳,电偶极子长100 m,距海底高50 m,发射频率和发射电流分别为1 Hz和1 A.为满足Dirichlet边界条件,整个模型区域为[-100 km, -100 km]×[-100 km, -100 km]×[-100 km, 50 km],我们设计了粗、细两套网格,粗网格将模型区域剖分成30×20×17个六面体单元,最小网格尺寸为500 m×500 m×100 m;细网格将区域剖分成111×52×39个六面体单元,最小网格尺寸为100 m× 100 m×20 m.我们在粗网格上应用二阶谱元法而在细网格上应用一阶谱元法.图 4为两种网格下所得同线Ex分量振幅和相位随偏移距变化曲线,以及相对于Dipole1D结果(Key, 2009)的相对误差.从图中可以看出,两种网格下谱元法计算的响应曲线与Dipole1D结果都能很好地吻合.误差曲线显示高阶谱元法应用于粗网格与低级谱元法应用于细网格可以获得相同的计算精度.为了进一步检验谱元阶数对计算精度的影响,我们在图 5中给出了粗网格下不同阶数(N=2, 3, 4)时Ex振幅和相位相对于Dipole1D的相对误差.从图中可以看出,随着多项式阶数的增长,正演精度显著提高,每提高一阶误差减小约一个数量级.N=4时谱元法获得的振幅和相位响应最大相对误差分别为0.026%和0.01%.通过对图 5分析可以得出如下结论:谱元法可以在稀疏网格下达到很高的精度,这是因为谱元法通过高阶多项式对各剖分单元内场的变化特征进行有效地刻画.相比之下,传统有限元方法由于采用线性基函数,为了有效刻画单元内场的特征,需要将网格细化到一定程度使得单元内的场近似为线性变化,因此模拟算法的有效性对网格的依赖性较大.
我们进一步设计了如图 6所示的三维地电模型,并与殷长春等(2014)给出的有限差分(finite-difference method,简称FDM)结果进行对比.模型海水深度1000 m,电阻率0.3 Ωm,高阻异常体电阻率100 Ωm,大小为6000 m×6000 m×100 m,顶部埋深为1000 m,中心在海底投影点坐标为(5000 m, 0 m, 1000 m).发射机为沿x方向的水平电偶源,发射频率0.25 Hz,发射电流1 A,距海底高30 m.我们设计了三套网格:网格1将模型区域划分为46×20× 19个网格,除扩边网格外,区域网格尺寸为500 m×500 m×200 m,异常体所在位置z方向网格长度细分为100 m;网格2将模型区域划分为106×32×33个网格,除扩边网格外,区域网格尺寸为200 m×200 m×100 m,异常体所在位置z方向网格长度细分为25 m;网格3将模型区域划分为206×52×53个网格,除扩边网格外,区域网格尺寸为100 m×100 m×50 m,异常体所在位置z方向网格长度细分为25 m.我们首先在粗网格1上应用本文二阶谱元法,并在不断细化的三种网格上应用殷长春等(2014)给出的有限差分算法.图 7展示了计算的模型MCSEM响应.从图中可以看出,随着网格的细化,有限差分响应结果逐渐向二阶谱元响应结果收敛,当采用网格3时,有限差分结果与谱元结果吻合较好,说明二阶谱元法在粗网格1上获得了有限差分在细网格3上同样的计算精度.表 1给出粗网格上谱元法和三种网格下有限差分法关于自由度和计算时间的对比结果.从表中可以看出:相对于细网格下有限差分法,在获得相同计算精度的前提下,粗网格下谱元法产生更少的未知数,计算时间也大大减少.
本文最后设计了如图 8所示的三维海洋地电模型,以分析异常体埋深和海水深度对MCSEM响应特征的影响.假设模型原点位于海面正中心,海水深度h0,高阻储油层电阻率为80 Ωm,嵌入于1 Ωm的海底围岩中,储油层大小为3000 m×3000 m×100 m,顶面中心坐标为(2500 m, 0, h0+h1).接受点位于(0,0,h0)处,水平电偶极子发射源沿x轴方向拖曳,系统参数及网格剖分情况与图 3中粗网格相同.
假设海水深度h0=1000 m,我们分别用3阶和4阶谱元法计算了储油层不同埋深情况下(h1=100 m/500 m/1000 m) Ex分量振幅和相位曲线(图 9)及归一化振幅和相位差曲线(图 10).从图中可以看出:3阶和4阶谱元法计算的MCSEM响应完全相同,归一化振幅和相位差曲线结果重合;由电流密度的连续性可知,在高阻储油层位置电场增大,Ex振幅和相位在油层上方都明显增大,且异常体埋深越浅,响应增幅越大.图 11和图 12展示了油层埋深为h1=500 m时不同海水深度(h0=100 m/500 m/1000 m)条件下Ex响应及归一化振幅和相位差曲线.从图中可以看出:不同海水深度的响应尾支趋于平行,表明空气波占主导地位.随海水深度减小,空气波出现越早,影响范围越大.当海水深度为100 m时,在2000 m的收发距下就观测到空气波,导致高阻层中的导波响应很小,归一化幅值接近于1.随着海水深度的增加,空气波在海水中垂向传播时消耗大部分能量,使得导波在接收信号中占的比例增大.当海水深度为1000 m时,归一化幅值约为17,相位差约为220°,明显反映高阻异常体的存在.
(1) 本文将谱元法成功应用于海洋可控源电磁三维正演模拟.与传统的有限元方法不同,利用高阶GLC正交多项式构造插值基函数,能够有效地刻画单元内场的变化特征,使得在稀疏网格下也能获得高精度的正演结果.数值实验结果表明谱元阶数提高对改善计算精度效果明显.和有限差分法的对比进一步表明谱元法对网格的依赖程度较之于传统方法较弱.因此,谱元法是三维海洋电磁正演模拟的有效方法.海底高阻储油层可以通过观测场与背景场的差异进行有效识别.海水越深或异常体埋藏越浅,响应差异越大,高阻储油层越容易被识别出来.
(2) 然而,本文仅初步尝试将基于GLC多项式的谱元法应用到电磁勘探领域,数值模拟是在比较简单的结构化网格下进行的.为了模拟复杂地质情况,可以考虑使用任意六面体单元剖分实现谱元算法,这将是我们未来的研究方向.
致谢感谢编辑部和审稿专家对本文提出的宝贵的意见和建议,同时感谢吉林大学“千人计划”电磁研究团队成员的积极支持和帮助.
Amestoy P R, Duff I S, L'Excellent J Y, et al.
2001. A fully asynchronous multifrontal solver using distributed dynamic scheduling. SIAM Journal on Matrix Analysis and Applications, 23(1): 15-41.
DOI:10.1137/S0895479899358194 |
|
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 |
|
Börner R U.
2010. Numerical modelling in geo-electromagnetics:advances and challenges. Surveys in Geophysics, 31(2): 225-245.
DOI:10.1007/s10712-009-9087-x |
|
Chen G B, Wang H N, Yao J J, et al.
2009. Modeling of electromagnetic responses of a 3-D electrical anomalous body in a layered anisotropic earth using integral equations. Chinese J. Geophys, 52(8): 2174-2181.
DOI:10.3969/j.issn.0001-5733.2009.08.028 |
|
Constable S.
2010. Ten years of marine CSEM for hydrocarbon exploration. Geophysics, 75(5): 75A67-75A81.
DOI:10.1190/1.3483451 |
|
Ellingsrud S, Eidesmo T, Johansen S, et al.
2002. Remote sensing of hydrocarbon layers by seabed logging (SBL):Results from a cruise offshore Angola. The Leading Edge, 21(10): 972-982.
DOI:10.1190/1.1518433 |
|
Han B, Hu X Y, Huang Y F, et al.
2015. 3-D frequency-domain CSEM modeling using a parallel direct solver. Chinese J. Geophys, 58(8): 2812-2826.
DOI:10.6038/cjg20150816 |
|
Hesthaven J S, Warburton T.
2002. Nodal high-order methods on unstructured grids:Ⅰ.Time-domain solution of Maxwell's equations. Journal of Computational Physics, 181(1): 186-221.
DOI:10.1006/jcph.2002.7118 |
|
Key K.
2009. 1D inversion of multicomponent, multifrequency marine CSEM data:Methodology and synthetic studies for resolving thin resistive layers. Geophysics, 74(2): F9-F20.
DOI:10.1190/1.3058434 |
|
Key K.
2012. Marine electromagnetic studies of seafloor resources and tectonics. Surveys in Geophysics, 33(1): 135-167.
DOI:10.1007/s10712-011-9139-x |
|
Komatitsch D, Tromp J.
1999. Introduction to the spectral-element method for 3-D seismic wave propagation. Geophysical journal international, 139: 806-822.
DOI:10.1046/j.1365-246x.1999.00967.x |
|
Lee J H, Xiao T, Liu Q H.
2006. A 3-D spectral-element method using mixed-order curl conforming vector basis functions for electromagnetic fields. IEEE Transactions on Microwave Theory and Techniques, 54(1): 437-444.
DOI:10.1109/TMTT.2005.860502 |
|
Li L, Liu T, Hu T Y.
2014. Spectral element method with triangular mesh and its application in seismic modeling. Chinese J. Geophys, 57(4): 1224-1234.
DOI:10.6038/cjg20140419 |
|
Li Y, Key K.
2007. 2D marine controlled-source electromagnetic modeling:Part 1-An adaptive finite-element algorithm. Geophysics, 72(2): WA51-WA62.
DOI:10.1190/1.2432262 |
|
Liu Y S, Teng J W, Xu T, et al.
2014. Numerical modeling of seismic wavefield with the SEM based on Triangles. Progress in Geophysics, 29(4): 1715-1726.
DOI:10.6038/pg20140430 |
|
Maday Y, Patera A T, Ronquist E M. 1987. A well-posed optimal spectral element approximation for the Stokes problem.
|
|
Patera A T.
1984. A spectral element method for fluid dynamics:laminar flow in a channel expansion. Journal of Computational Physics, 54(3): 468-488.
DOI:10.1016/0021-9991(84)90128-1 |
|
Ren Q, Tobón L E, Sun Q, et al.
2015. A new 3-D nonspurious discontinuous Galerkin spectral element time-domain (DG-SETD) method for Maxwell's equations. IEEE Transactions on Antennas and Propagation, 63(6): 2585-2594.
DOI:10.1109/TAP.2015.2417891 |
|
Runge C.
1901. Vber empirische Funktionen und die Interpolation zwischen äquidistanten Ordinaten. Zeitschrift für Mathematik und Physik, 46: 224-243.
|
|
Seriani G, Priolo E.
1991. High-order spectral element method for acoustic wave modeling. SEG Technical Program Expanded Abstracts 1991. Society of Exploration Geophysicists: 1561-1564.
DOI:10.1190/1.1888989 |
|
Sun X Y, Nie Z P, Zhao Y W, et al.
2008. The electromagnetic modeling of logging-while-drilling tool in tilted anisotropic formations using vector finite element method. Chinese J. Geophys, 51(5): 1600-1607.
|
|
Tape C, Liu Q, Maggi A, et al.
2010. Seismic tomography of the southern California crust based on spectral-element and adjoint methods. Geophysical Journal International, 180(1): 433-462.
DOI:10.1111/j.1365-246X.2009.04429.x |
|
Wang X M, Seriani G, Lin W J.
2007. Some theoretical aspects of elastic wave modeling with a recently developed spectral element method. Science in China Series G:Physics Mechanics and Astronomy, 50(2): 185-207.
DOI:10.1007/s11433-007-0022-1 |
|
Yang J, Liu Y, Wu X P.
2015. 3D simulation of marine CSEM using vector finite element method on unstructured grid. Chinese J. Geophys, 58(8): 2827-2838.
DOI:10.6038/cjg20150817 |
|
Yin C C, Ben F, Liu Y H, et al.
2014. MCSEM 3D modeling for arbitrarily anisotropic media. Chinese J. Geophys, 57(12): 4110-4122.
DOI:10.6038/cjg20141222 |
|
Yin C, Huang X, Liu Y, et al.
2017. 3-D Modeling for Airborne EM using the Spectral-element Method. Journal of Environmental and Engineering Geophysics, 22(1): 13-23.
DOI:10.2113/JEEG22.1.13 |
|
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.
DOI:10.1190/1.2358403 |
|
陈桂波, 汪宏年, 姚敬金, 等.
2009. 用积分方程法模拟各向异性地层中三维电性异常体的电磁响应. 地球物理学报, 52(8): 2174–2181.
DOI:10.3969/j.issn.0001-5733.2009.08.028 |
|
韩波, 胡祥云, 黄一凡, 等.
2015. 基于并行化直接解法的频率域可控源电磁三维正演. 地球物理学报, 58(8): 2812–2826.
|
|
李琳, 刘韬, 胡天跃.
2014. 三角谱元法及其在地震正演模拟中的应用. 地球物理学报, 57(4): 1224–1234.
|
|
刘有山, 滕吉文, 徐涛, 等.
2014. 三角网格谱元法地震波场数值模拟. 地球物理学进展, 29(4): 1715–1726.
|
|
王秀明, 林伟军.
2007. 利用谱元法计算弹性波场的若干理论问题. 中国科学:G辑, 37(1): 41–59.
|
|
杨军, 刘颖, 吴小平.
2015. 海洋可控源电磁三维非结构矢量有限元数值模拟. 地球物理学报, 58(8): 2827–2838.
|
|
殷长春, 贲放, 刘云鹤, 等.
2014. 三维任意各向异性介质中海洋可控源电磁法正演研究. 地球物理学报, 57(12): 4110–4122.
|
|