地球物理学报  2015, Vol. 58 Issue (3): 1059-1071   PDF    
复杂场源形态的海洋可控源电磁三维正演
韩波1,2, 胡祥云1, Adam SCHULTZ2, 左博新1, 李建慧1, 蔡建超1    
1. 中国地质大学(武汉)地球物理与空间信息学院, 武汉 430074;
2. 俄勒冈州立大学地球与海洋与大气科学学院, 科瓦利斯, 美国
摘要:在使用电偶极发射源的可控源电磁法(CSEM)勘探中,发射源的方位、长度、形状等对观测数据有重要的影响,然而现有的大部分三维数值模拟方法没有全面地将这些因素考虑进来,很多都只能应对非常简单的场源形态,例如单一方位的点电偶极子,这有可能显著降低模拟结果的准确性.本文实现了基于交错网格有限体积(FV)离散的海洋CSEM三维正演算法,能够模拟形态相对复杂的场源,包括任意方位的有限长直导线和弯曲导线发射源.该算法使用一次场/二次场方法,只需对二次场使用FV法求解,避免了场源的奇异性问题;一次场的计算为一维正演问题,使用准解析法求解,并且只要在计算一次场时考虑复杂的场源形态便可以实现同样场源的三维正演.通过与一维理论模型的解析解对比验证了三维程序的准确性,并针对三维理论模型进行了一系列正演测试,初步考察了场源形态对三维正演结果的影响.
关键词海洋可控源电磁法     三维正演     有限长线源     场源形态    
Three-dimensional forward modeling of the marine controlled-source electromagnetic field with complex source geometries
HAN Bo1,2, HU Xiang-Yun1, Adam SCHULTZ2, ZUO Bo-Xin1, LI Jian-Hui1, CAI Jian-Chao1    
1. Institute of Geophysics and Geomatics, China University of Geosciences, Wuhan 430074, China;
2. College of Earth, Ocean, and Atmospheric Sciences, Oregon State University, Corvallis, OR, 97331-5503, USA
Abstract: In controlled-source electromagnetic (CSEM) surveys that utilize electric dipole sources, observations can be greatly influenced by the orientation, length and the shape of the source. However, most of the current 3D CSEM modeling algorithms do not take all these factors into account, and many are only able to deal with very simple source geometries such as a point dipole with a fixed orientation, which may cause significant inaccuracies in the modeling result. A 3D forward modeling code for marine controlled-source electromagnetic (CSEM) surveys has been developed, which permits simulations of complex source geometries.
The solution is based on a staggered-grid finite volume (FV) discretization of the second-order Maxwell's equations. The primary/secondary field approach, in which only the secondary fields need to be solved using the FV method, is employed to avoid the source singularity problem. To accurately compute the primary fields over an 1D layered background model, the open source software Dipole1D is integrated into the 3D code. Due to Dipole1D's versatility and the fact that the primary field calculation is the only part of the 3D modeling that directly handles the source geometry, the 3D code is capable of simulating arbitrarily orientated finite-length wire sources. The resulting linear system of FV equations is solved by the quasi-minimum residual (QMR) method.
To validate the 3D code, we first benchmarked FV solutions against quasi-analytic solutions for a canonical 1D reservoir model. The numerical solutions show excellent agreement with analytic solutions for various transmitter orientations for both a point dipole source and a 300 m wire source. Then we calculated responses of a 3D reservoir model similar to a well-studied one, varying horizontal sizes of the resistor. Reasonable results for the 3D model, which are consistent with published work, can be observed.
Further, a suite of numerical experiments on the 3D model were conducted to investigate the influences of source geometries, and the observations include: (1) The full lateral extent of the 3D resistor can be well estimated from the responses normalized by responses of the background model without the resistor, by moving and rotating the transmitter. (2) Electromagnetic fields excited by wire sources with wire lengths ranging from 100 m to 500 m are compared with those by a point dipole, showing that the longer the source, the greater the differences, and the EM field amplitudes of 500 m wire source can be up to more than 4.5 times larger than that of a point dipole. Prominent differences occur not only in the region immediately around the source, but also in four narrow regions at oblique angles to the source. (3) While the 300 m straight wire source is slightly bended in four different ways, either in the horizontal plane or in the vertical plane, one or several components of the EM field change significantly.
From results of all the numerical experiments, we can conclude that the presented 3D FV CSEM modeling code is accurate and versatile, capable of modeling relatively complicated source geometries, such as finite-length straight/bent wire sources with arbitrary orientations. In CSEM surveys, source geometries, i.e. the orientation, the length and the shape of the source can have remarkable influence on the observed data, thus accounting for source geometry effects is necessary.
Key words: Marine CSEM     3D modeling     Finite-length source     Source geometry    
1 引言

海洋可控源电磁法(controlled-source electromagnetic methods,CSEM)由于对高电阻率的油气储集层有很强的识别能力,在最近十多年引起了勘探地球物理学家们极大的关注,其仪器、数据采集技术、解释 技术等均获得了迅猛发展(Edwards,2005Constable,2010Key,2012a),并已经成为海洋油气勘探中一种不可或缺的地球物理技术.如今,海洋CSEM的反演解释正在迈入三维化阶段.正演是反演的引擎,反演解释的准确性与高效性很大程度上依赖于正演.CSEM中最常用的三维正演方法包括积分方程法(IE)(Avdeev and Knizhnik, 2009Zaslavsky et al., 2011),有限差分法(FD)(Newman and Alumbaugh, 1995Weiss and Newman, 2002沈金松,2003Streich,2009Mittet,2010)和有限单元法(FE)(Badea et al., 2001Mukherjee and Everett, 2011Schwarzbach et al., 2011da Silva et al., 2012Puzyrev et al., 2013).已有不少作者对这几种方法进行过总结(Avdeev,2005汤井田等,2007Börner,2010),一般来说,从IE到FD再到FE,所能模拟的地质构造的复杂性逐渐增强,但同时计算量也逐渐增加.若只进行极少次数的正演计算,具有灵活的离散网格的FE法或许是最佳选择;然而一次反演中往往包含成千上万次正演,FE法庞大的线性方程可能会导致无法接受的计算成本.此时可考虑使用计算量较小且精度也足够的FD法.在交错网格FD法基础上发展起来的有限体积法(finite-volume method,FV),由于对模型参数和电磁场在控制体积内进行积分而具有比FD法更高的精度(Haber and Ascher, 2001Weiss and Constable, 2006杨波等,2012李建慧等,2013);对FD/FV法进行一定的改进可模拟简单的起伏地形(杨波等,2012).

尽管现存的CSEM三维正演算法多种多样且对复杂介质的模拟能力越来越强,但其中全面考虑到场源形态复杂性(指场源的方位、长度、形状等)的还不多,多数将场源简化为导线长度无穷小的点电偶极子,有的甚至只允许场源固定在某个位置(如空气-大地界面)、为单一方向.陆地CSEM勘探中发射源接地导线长度通常在1 km以上,甚至可能达到4 km;海洋油气CSEM勘探中悬浮发射导线的典型长度为100~300 m(Constable and Srnka, 2007).Ward和Hohmann(1988)指出,只有当观测点到发射源中心的距离不小于发射导线长度的5~10倍时,发射源才能近似为点电偶极子.依照这个标准,很多实际的陆地CSEM勘探和海洋CSEM勘探的发射源都不应该被视为点电偶极子,而应该作为有限长度线源来对待.此外,由于地形起伏或海水的流动,发射导线可能产生方位上的偏离,自身也可能弯曲,几何形态变得复杂.场源的形态无疑会对CSEM的观测数据造成影响.Streich和Michael(2011)对点源与线源的电磁场作过较为全面的比较,认为发射导线的长度、方位等可能对勘探靶体的正演响应产生巨大的影响.因此,CSEM三维正演中考虑场源的复杂形态是很有必要的.

本文基于CSEM模拟中常用的二次场公式,使用有限体积法来求解二次场控制方程;使用一维准解析法计算一次场时,允许电偶极子为任意方位、位于模型中的任意位置,并对点源进行数值积分来获得线源的电磁场,对多个场源同时激发的电磁场进行线性叠加,最终实现了任意方位的有限长直导线和弯曲导线发射源激发下的海洋CSEM三维正演.通过理论模型的正演计算初步考察了场源形态对三维正演结果的影响.

2 方法理论 2.1 二次场方程

频率域CSEM属于典型的低频电磁法,在绝大部分已知的传播介质的电阻率范围内,可以忽略位移电流.则频率域Maxwell方程的微分形式为(取时谐因子为eiωt):

其中,H为磁场强度(A/m),E为电场强度(V/m),σ为电导率(S/m),ω为角频率,μ0为真空中的磁导率(H/m),J为外加场源的电流密度(A/m2).若直接对以上方程进行数值求解,由于J的存在,场源点的奇异性(J为无穷大)将给三维数值模拟带来很大的困难.为此,引入一个参考模型σp,在场源J激发下的电磁场同样满足上述形式方程:

将式(1)、(2)分别与式(3)、(4)相减,可得:

整理式(5)右端,重新写为

若将所计算模型的场(总场)看作参考模型所产生的一次场(primary field)与“剩余模型”所产生的二次场之和,即令二次场Hs=HHpE s=EEp;将式(7)右端第二项看作激发二次场的电流密度,即令Js=(σ-σp)E p,则可将式(6)和式(7)写成二次场方程的形式:

比较二次场方程(8)、(9)与总场方程(1)、(2),由于二次场电流密度Js不存在奇异点,因此解二次场方程要比解总场方程容易得多.对式(9)取旋度并将式(8)代入其中,得到关于电场的二阶偏微分方程:

对方程(10)进行数值法求解得到二次场;对于一次场,一般将参考模型选为所需计算模型的一维“背景模型”,可为全空间或层状模型,因此可以通过解析法求得;一次场与二次场相加便得到总场.

2.2 有限差分法与有限体积法

求解类似于方程(10)的Helmholtz方程,基于Yee网格(Yee,1966)的交错网格有限差分法(staggered finite-difference,SFD)可能是被使用得最多的方法(Newman and Alumbaugh, 1995Weiss and Newman, 2002Streich,2009).将研究区域离散化为一系列规则长方体单元,然后在长方体单元上交错地对电场和磁场采样.这里选择如图 1所示的采样方式,将电场定义在单元棱边的中点,将磁场定义在单元面的中心点.

图 1 电磁场的交错采样Fig. 1 Sampling positions of EM fields on a staggered grid

对于有限差分离散方程的推导,很多文献中都 有详细的过程(例如Siripunvaraporn,1999),为了节约篇幅本文不再列出.使用SFD法对式(10)进行离散化直接得到的线性方程系数矩阵是不对称的,为了使其对称以便于方程的求解,每个采样点的线性方程两端需要同时乘以一个对应于该采样点的“尺度因子”V,V具有某种体积的意义,即SFD线性方程实际上对应的微分方程为:

有限体积法的电磁场采样方式与SFD完全一样,不同之处在于FV在每个采样点处定义一个包围该采样点的控制体积,如图 2所示.需要注意的是,控制体积与图 1中的网格单元是不同的,某采样点处的控制体积由共享该采样点所在棱边的相邻4个网格单元体积的1/4组成.对式(10)在某个体积单元Ω内求积分:

事实上,SFD对应的方程(11)中的尺度因子刚好是FV法中的控制体积值,因此对比式(11)和式(12),SFD法与FV法的区别在于各种变量是直接 取采样点处的值,还是取在一定范围内积分后的值.式(12)左端第一项可以完全展开(Weiss and Constable, 2006),并且与式(11)左端第一项具有完全相同的离散形式.对于式(12)的左端第二项,通常认为二次电场在控制体积范围内的变化不大,因此Es可提取到积分号以外;采样点处的电导率可与SFD中一样,使用采样点所在棱边的相邻4个单元电导率的体积加权平均值,也不必求积分.经过这样的简化,FV的左端项就与SFD完全相同,因此二者的线性方程系数矩阵相同;唯一不同的是右端项:FV法对一次场在控制体积内求积分,SFD法则简单地以采样点 处的场值来代表整个控制体积内的场值.显然,FV在电磁场变化剧烈的地方如场源附近,比SFD法更准确.

由于CSEM的场源是局部的,因此可让剖分域的边界离场源足够远,从而运用齐次的Dirichlet边界条件,即令边界上的切向电场为0.

图 2 三个电场分量采样点处的控制体积Fig. 2 Control volumes for sampling sites of three E components

FV法右端项的一次场体积积分使用高斯-勒让德积分公式计算(Davis and Rabinowitz, 2007);离散大型线性方程系数矩阵为对称且高度稀疏的复矩阵(仅对角元为复数),对其预优处理后使用QMR(拟最小残差法,Quasi-Minimal Residual)迭代解法求解(Saad,2003).

2.3 一次场计算——一维正演问题

根据前文,一次电场与“剩余模型”电导率的乘积具有二次电场的场源电流密度的物理意义,出现在FD/FV离散线性方程的右端项中.由于一次场为人工场源在一维背景模型中激发的场,所以一次场的计算实际上是一维正演问题,并且只要在一次场计算中考虑到人工场源的方位、长度等,就可以实现同样的场源激发下的三维正演.对于全空间模型这种最简单的情况,电偶极子源激发的电磁场可以用纯解析公式求得(Ward and Hohmann, 1988);对于水平层状模型(包括半空间),需要使用Hankel变换来求贝塞尔函数的积分,属于“准解析”方法.在本文的三维正演代码的实现中,植入了开源的Dipole1D程序(Key,2009)作为一维正演的核心.Dipole1D程序基于电磁场的势的公式,允许发射电偶极子位于层状介质中的任意位置,并可为任意方位,同时可对点电偶极子沿发射导线方向进行线积分,从而获得有限长线源的电磁场.这里简要介绍一下任意方位的电偶极子和有限长线源这两个特征的实现.给定直角参考坐标系,假设偶极子朝y轴正向偏离x轴正向的角度为φ,朝z轴正向偏离水平面的角度为θ,偶极矩为m,如图 3所示.

图 3 电偶极子的方位示意图Fig. 3 The orientation of an electric dipole

在水平电偶极子和垂直电偶极子的计算已经实现的前提下,可将任意方位的电偶极子在水平面内和垂直面内投影,分解为一个水平电偶极子和一个垂直电偶极子,偶极矩分别为mcosθ和msinθ.分别计算水平电偶极子和垂直电偶极子的电磁场,并投影到参考坐标轴方向后进行线性叠加,便可得到任意旋转的电偶极子的电磁场.

对于有限长度线源,其激发的电磁场等于导线上点源的场沿导线方向的积分.依然使用高斯-勒让德求积公式来数值计算线积分,即:

其中F(L)表示长度为L的线源的场;f(x)表示导线上相对坐标为x的点源的场;n为高斯积分阶数;xk和Wk分别为高斯-勒让德求积公式的节点和系数(需要根据积分区间改变尺度).可以看出,这种方法的本质是对导线上有限个离散点处的点源激发的场进行加权求和来获得线源的场.离散点个数越多,积分阶数越高,线积分的结果就越准确,但同时计算量也越大. 3 算法正确性检验——一维层状模型的计算

通过存在准解析解的一维层状模型来验证FV三维正演程序的准确性.选用一个海洋CSEM中常 见的模拟高阻油气薄层的模型(Constable and Weiss, 2006Key,2009),如图 4所示,在背景电阻率值为1 Ωm的海底地层中包含一层100 Ωm的油气储层,位于海底以下1000~1100 m深度处;海水厚度1000 m,电阻率0.3 Ωm;空气的电阻率取1011 Ωm.观测系统设置为:发射电偶极位于海底上方50 m,取偶极中心在x-y平面上的投影为坐标原点,发射 频率为1 Hz;接收器位于海底表面,沿x轴方向布设.

图 4 海洋油气储层层状模型以及一个典型的海洋CSEM施工示意图Fig. 4 The 1D layered reservoir model and a typical marine CSEM survey geometry
3.1 点电偶极子

首先让点电偶极处于x轴正方向(φ=0°,θ=0°),然后在x-y平面内向y轴正方向旋转两次,每次旋转的角度为30°,即θ=0°不变,φ分别为30°和60°;紧接着让电偶极在垂直面内向z轴正向旋转两次,每次旋转的角度为30°,即φ=60°不变,θ分别为30°和60°.分别使用FV三维正演程序与Dipole1D程序对这5种不同方位的点电偶极的电磁场计算.三维正演时剖分域包含5层空气层(厚度为9 km);参考模型为空气-海水-海底地层3层模型(即去掉了原模型中的高阻油气层);离散方程右端项的一次场的高斯积分阶数取1;三维正演程序中的主场计算以及Dipole1D程序中,快速Hankel变换均使用Key的201点滤波器(Key,2012b).图 5a—5d给出了单位偶极矩激发的水平电场的振幅-收发距(MVO)曲线与相位-收发距(PVO)曲线.可以看出,对于所有5个方位的电偶极子,Ex分量的FV三维数值解与解析解高度吻合,而Ey的相位在收发距 较大(>3.5 km)时数值解与解析解出现了较明显的偏差,可能跟从网格节点到测点的场值插值方法有关.

图 5 不同方位的点电偶极子在图 4所示的层状模型上激发的水平电场的FV三维数值解(符号)与Dipole1D解析解(线条)对比Fig. 5 Horizontal electric fields for the layered model shown in Fig. 4,FV solution(symbols) and the Dipole1D analytic solution(lines)for a point dipole source with different orientations

由于接收器沿x方向布设且与发射偶极的中 心共线,根据电偶极子场的特征(Constable and Weiss, 2006Figure 2),x方向的偶极矩在测线上(此时为inline观测)激发的一次电场主要为x分量和z分量,y分量极微弱;y方向的偶极矩在测线上(此时为broadside观测)激发的一次电场主要为y分量,x分量和z分量极微弱;z方向的偶极矩在测线上激发的一次电场主要为x分量和z分量,y分量极微弱.因此,单独看图 5a—5d中偶极子在水平面内旋转时(仅φ变化)3个不同方位的电磁场,当偶极子从x轴逐渐向y轴旋转时,电偶极矩在x方向的投影越来越小,使得主要由x方向偶极矩产生的Ex分量振幅逐渐减小,直到偶极子完全处于y方向时(φ=90°)将会消失;相反,电偶极矩在y方向的投影越来越大,使得主要由y方向偶极矩产生的Ey分量振幅逐渐增强,直到偶极子完全处于y方向时将达到最大.但这3个方位的电场相位是重合的,这是因为此时Ex和Ey都只主要由某个单一方向的偶极矩产生,φ的改变只是改变了各自方向的偶极矩的幅度,导致电场的实部和虚部等比例地改变,并不影响复数电场的相位.若单独看偶极子在垂直面内旋转时(仅θ变化)3个不同方位的电磁场,可以发现当偶极子从水平方向逐渐向垂直方向旋转时,电偶极矩的水平分量越来越小,y分量越来越小,导致Ey振幅逐渐减小,直到偶极子完全垂直时(θ=90°)将会消失,但相位不会改变;而Ex的振幅不会随θ的变化而单向变化,且θ变化时相位也会改变,这是因为当θ>0°时偶极矩产生了z分量,Ex是偶极矩的x分量和z分量共同作用的结果,θ的增大会导致前者减小而后者增大,Ex的实部和虚部不会简单地单向且等比例变化.

3.2 有限长线源

下面在正演模拟中考虑发射导线的长度.图 6显示了对应于3.1节中的每种情况下的300 m长线源的正演响应,对场源的线积分使用10个高斯点.可以看出,FV三维数值解与解析解仍然吻合得非常好,说明本文的三维正演程序对任意方位的有限长线源的计算是准确的;线源旋转时电磁场的变化规律与点源类似.

图 6 不同方位的线源(300 m)在图 4所示的层状模型上激发的水平电场的FV三维数值解(符号)与Dipole1D解析解(线条)对比Fig. 6 Horizontal electric fields for the layered model shown in Fig. 4,FV solution(symbols) and the Dipole1D analytic solution(lines)for a 300 m long wire source with different orientations

Ex为例来简单对比一下该层状模型中点源和线源激发的电磁场.图 7为使用点源对长度在50~500 m之间的线源产生的Ex分量进行归一化(振幅相比,相位相减)的结果,所有场源均为x方向.可以看出,导线长度越大,其正演响应与点电偶源的差别越大,尤其是在场源附近,100 m以上长度的导线与点源的差别非常明显,500 m长导线的水平电场振幅可达点电偶源的4.5倍以上;随着收发距的增大,二者差别逐渐减小,收发距达到2 km时500 m长的导线与点源已经几乎看不出差别.

图 7 使用点源对不同长度的线源在图 4所示模型上激发的Ex归一化后的结果
(a)为振幅比,(b)为相位差.场源均为x方向.
Fig. 7 Ex for different wire lengths normalized by point dipole for the model displayed in Fig. 4
(a)Amplitude ratio and (b)phase difference. All sources are x-directed.
4 三维模型算例

在第3节层状模型的基础上,将高阻储层在水平方向的边界截断,则得到更接近真实情况的三维油气储层模型,如图 8所示. 类似的模型曾被多次研究(Constable and Weiss, 2006Weiss and Constable, 2006Streich,2009),这些作者均使高阻体在x-y平面内的截面为圆形.为了避免模拟弧形边界导致的大量网格,这里让高阻体在x-y平面内的截面为正方形.

图 8 三维油气储层模型Fig. 8 The 3D reservoir model

分别对高阻体截面边长(设为a)为1 km,2 km和5 km的三种情况进行了正演.参数设置为:场源 位于高阻体x负向边界中点的正上方,海底上方50 m,场源坐标为(0,0,-50 m);发射频率为1 Hz;接收器分布于海底y=0的测线上(inline观测),x坐标范围为-1000~8000 m,测站间距50 m,则总共有181个测站;发射偶极长度100 m,使用4点高斯积分;构建FV体积积分右端项使用1个高斯积分点.所得到的水平电场Ex分量如图 9所示.图 9中还给出了不含高阻体的半空间(a=0)和高阻体在水平方向无限延伸(a=∞)这两种极限情况的响应(使用解析法获得).

图 9 三维储层模型的储层截面边长变化时的Ex振幅(a)及相位(b)曲线Fig. 9(a)Ex amplitude and (b)phase for the 3D reservoir model shown in Fig. 8, with various horizontal edge of the reservoir

可以看出,在小收发距(r<1 km)内,经海水传 播的直达波能量占主导地位,地层几乎毫无影响,电磁场振幅大约按1/r3衰减,类似于陆地CSEM中的“近场效应”;随着收发距的增大,电磁波在海水中大大衰减,地层波开始起主导作用,含高阻层(a>0)与不含高阻层(a=0)的电磁场逐渐区分开来,前者幅值更强,这是由于电磁波在地层中发生了长距离折射,折射波在高阻介质中衰减较慢而在低阻介质中衰减较快;电磁场曲线在对数坐标系中的斜率趋于常数,说明电磁场趋于按指数衰减,衰减因子取决于折射层电阻率.a=1 km与半空间模型的电磁场曲线几乎重合,说明该尺寸的异常体难以分辨;a= 2 km与半空间的电磁场虽然在大收发距时的衰减 速度保持一致,但二者的绝对差值较为明显,说明该尺寸的异常体容易分辨;a=5 km与半空间的电磁场差别非常明显,且在收发距大于高阻体边长时与高阻体无限延伸(a=∞)的电磁场区分开来,衰减速度逐渐趋于半空间的衰减速度.

根据以上分析并且与公开发表的类似模型的结果对比,不仅能证明本文的三维正演算法对三维模型的计算是准确的,而且可以再次确认Constable和Weiss(2006)的结论:对于此类水平薄板状高阻体,只有当其水平尺寸至少为埋深的2倍时,才可以利用CSEM将其分辨出来.

4.1 多方位观测

上述分析表明inline观测方式对inline方向的电性边界有灵敏的反应.为了将三维储层在水平方向的延伸范围更完整地反映出来,仅仅做一次观测是不够的,通常需要移动并旋转场源进行多个方位的观测.对以上储层边长为2 km和5 km的两个模型,让发射源从储层边缘正上方(图 10a,10e)移动到储层中心正上方(图 10b,10f),然后朝y轴正向分别旋转45°(图 10c,10g)和90°(图 10d,10h).计算了-2 km≤x≤8 km,-4 km≤y≤4 km的测网上的电磁场,并将水平分量投影到沿场源方向.为了突出储层的响应,使用无储层的半空间模型进行归一化,得到的归一化的水平电场振幅如图 10所示.可以看到,以broadside测线(穿过偶极中心且垂直于偶极方向)为轴的一定旋转角度范围内是一个“盲区”,电性边界几乎没有任何反映,而沿发射偶极方向的电性边界都被准确地反映出来.通过这种多方位的观测,最终准确地圈定了储层在水平方向的延伸范围,说明进行不同方位的场源的模拟是非常重要的.

图 10 储层边长为2 km(a—d)和5 km(e—h)的模型在场源移动并旋转时的归一化水平电场分量Fig. 10 Normalized horizontal electric filed amplitude with the source moving and rotating for the model depicted in Fig. 8, with the reservoir edge length of 2 km(a—d) and 5 km(e—h)
4.2 场源长度的影响

以上均为100 m长线源的计算结果,图 11图 12则分别显示了点源和多个长度的线源的正演结果.这里使用储层边长为5 km的模型,场源位于储层中心正上方(图 10f).图 11a—11d分别为点源和100 m、300 m、500 m长的线源激发的Ex振幅在二维测网上的分布,从中可以粗略看出发射源长度的影响:点源的场的极大值区域聚焦于一个点,线源的场的极大值区域则聚焦于一条线,且场源长度越长,这条线的长度越长.图 12a图 12b分别显示了一个相对较小的测网上300 m线源被点源归一化的Ex振幅和相位,可以看出,线源与点源的场差异最大的地方除了场源附近的一小片区域外,还包括一个延伸得较远的“X”型区域,在“X”的4个端点附近场的差异甚至与场源处相当.图 12c图 12d分别为inline和broadside测线上的线源与点源的Ex振幅比,可以看出,在inline测线上场源中心点附近点源的场幅值强于线源的场,而当远离场源中心点一定距离时,则是线源的场较强,这是由于点源的能量都集中于一点,而线源的能量则均匀分布在导线上;对于broadside测线,在收发距不太大时点源的场幅值始终强于线源的场,且发射导线越长场越弱,这是因为点源始终处在导线上离broadside测点最近的地方.

图 11 三维储层模型在点源与不同长度的线源激发下的Ex振幅(以10为底的对数)Fig. 11 Ex amplitude(on a log10 scale)for wire sources with various lengths for the 3D reservoir model

图 12 三维储层模型线源被点源归一化的正演响应
(a)(b)分别为归一化的300 m线源的Ex振幅和相位在平面区域的分布,(c)(d)分别为归一化的不同长度线源的Ex振幅在inline和broadside测线上的分布.
Fig. 12 Forward responses of wire sources normalized by point dipole source for the 3D reservoir model Normalized Ex
(a)amplitude and (b)phase for a 300 m-long wire, and normalized Ex amplitude for wire sources with various lengths at(c)inline sites and (d)broadside sites.
4.3 场源发生弯曲的影响

海洋CSEM勘探经常将发射导线悬浮在海水中,很容易受海水流动的影响,使导线不仅产生方位上的偏离,还会弯曲.对曲线的模拟无非是使用多个分段直线来近似.由于我们使用了一次场/二次场方法,因此只需在一次场计算中(一维正演)对多个同时发射的不同方位的直线源的场叠加便可以构建出曲线源的方程右端项,实现曲线源的三维正演模拟.

设计了几种简单的弯曲导线,如图 13所示(黑色实线),其中导线长度为300 m,平均分成3段发生偏离,灰色虚线表示未弯曲的导线;图 13a—13d为导线在水平面内弯曲,图 13e—13h为导线在垂直面内弯曲.仍然使用储层边长为5 km的模型,场源位于储层中心正上方(如图 10f).图 13显示的是弯 曲导线与直导线的电场振幅比.单独来看,图 13a,13b表明这种弯曲方式会使导线弯向的一侧场值变强,另一侧场值变弱,显然这是由于跟直线源相比一侧离场源的距离减小而另一侧离场源的距离变大了;Ey分量在inline测线两侧会发生从偏大到偏小或从偏小到偏大的突变,说明如果接收器没有严格地布设在inline测线上,则场源的这种弯曲会对观测结果造成很大的影响;正如所预期的,水平电场的差异是关于y轴对称的.在图 13c,13d的这种弯曲方式下,Ex和Ey分量的差异不再关于任何轴对称,而是近似关于场源中心点旋转对称;Ey分量的差异在inline测线和broadside测线两侧均会突变.当导线在垂直面内发生如图 13e,13f所示时的小角度弯曲时,电场差异非常小,这是因为导线的第一段和第三段在方位上的偏离相反,互相抵消了一部分影响.图 13g,13h的这种弯曲方式下Ez分量的差异会在broadside测线两侧发生突变.若将图 13e,13f与13g,13h对比,将图 13a,13b与图 13c,13d对比,可以发现当导线端点的连线发生方位上的偏离时,其影响明显比未偏离时大.

总的来说,场源的弯曲形式不同,影响各不相同,但基本上都会对某一种或几种电磁场分量造成非常显著的影响,因此在利用多分量的数据进行解释时,应该将场源弯曲的因素考虑进来.

图 13 三维储层模型弯曲线源与直线源的电场振幅比Fig. 13 The ratio of electric fields′ amplitude of bent wire sources to that of a straight wire source for the 3D reservoir model
5 结论与讨论

本文基于一次场/二次场方法,实现了考虑复杂场源形态的海洋CSEM三维有限差分/有限体积法正演.在一次场计算中,将水平电偶极子和垂直电偶极子矢量合成来模拟任意方位的电偶极子,对点源进行高斯-勒让德数值积分来获得线源的电磁场,对多个同时发射的直线场源的场进行线性叠加来获得弯曲场源的场.通过与一维理论模型的解析解对比证明了三维正演算法的准确性.通过三维理论模型初步考察了场源形态对正演响应的影响,结果表明:要想探测出完整的电性边界,必须使场源位于多个不同的方位进行观测;场源的长度对收发距不太大的区域的影响较大;场源发生弯曲时一般会对某一种或几种电磁场分量造成非常显著的影响,因此在利用多分量的数据进行解释时,应该考虑到场源的弯曲.总而言之,考虑复杂场源形态的三维正演对提高CSEM勘探的准确性是很重要的.

与海洋CSEM相比,在陆地CSEM勘探中考虑复杂的场源形态或许更具有实际意义.正如文中所提到的,典型的陆地频率域CSEM勘探如CSAMT法所使用的发射导线长度通常为1~4 km,而接收点受制于多种因素经常无法布置在离场源足够远的地方,不仅平面波假设的条件难以满足,点电偶极子近似也非常粗糙;此外,受地形条件的影响,发射导线弯曲并且偏离所设计的方向是很容易产生的状况.本文的三维正演算法框架完全适用于接地线源的陆地频率域CSEM方法,但由于少了能够迅速吸收电磁波能量的海水层的存在,在一些细节上还需改动,比如边界条件的选取,线性方程组的求解等等.这也是我们接下来将会考虑的问题.

致谢     感谢QMR程序和Dipole1D程序的开发者们公开他们的源代码,没有他们无私的贡献,本文的研究将难以完成;感谢两位匿名审稿专家提出的宝贵意见使得本文的质量得到显著提升.
参考文献
[1] Avdeev D B. 2005. Three-dimensional electromagnetic modelling and inversion from theory to application. Surveys in Geophysics, 26(6): 767-799.
[2] Avdeev D, Knizhnik S. 2009. 3D integral equation modeling with a linear dependence on dimensions. Geophysics, 74(5): F89-F94.
[3] 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.
[4] Börner R U. 2010. Numerical modelling in Geo-Electromagnetics: advances and challenges. Surveys in Geophysics, 31(2): 225-245.
[5] Constable S, Srnka L J. 2007. An introduction to marine controlled-source electromagnetic methods for hydrocarbon exploration. Geophysics, 72(2): WA3-WA12.
[6] Constable S. 2010. Ten years of marine CSEM for hydrocarbon exploration. Geophysics, 75(5): A67-A81.
[7] Constable S C, Weiss C J. 2006. Mapping thin resistors and hydrocarbons with marine EM methods: Insights from 1D modeling. Geophysics, 71(2): G43-G51.
[8] 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.
[9] Davis P J, Rabinowitz P. 2007. Methods of Numerical Integration. 2nd ed. New York: Dover Publications.
[10] Edwards N. 2005. Marine controlled source electromagnetics: principles, methodologies, future commercial applications. Surveys in Geophysics, 26(6): 675-700.
[11] Haber E, Ascher U M. 2001. Fast finite volume simulation of 3D electromagnetic problems with highly discontinuous coefficients. SIAM Journal on Scientific Computing, 22(6): 1943-1961.
[12] 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.
[13] Key K. 2012a. Marine electromagnetic studies of seafloor resources and tectonics. Surveys in Geophysics, 33(1): 135-167.
[14] Key K. 2012b. Is the fast Hankel transform faster than quadrature? Geophysics, 77(3): F21-F30.
[15] Li J H, Hu X Y, Zeng S H, et al. 2013. Three-dimensional forward calculation for loop source transient electromagnetic method based on electric field Helmholtz equation. Chinese J. Geophys. (in Chinese), 56(12): 4256-4267, doi: 10.6038/cjg20131228.
[16] Mittet R. 2010. High-order finite-difference simulations of marine CSEM surveys using a correspondence principle for wave and diffusion fields. Geophysics, 75(1): F33-F50.
[17] Mukherjee S, Everett M E. 2011. 3D controlled-source electromagnetic edge-based finite element modeling of conductive and permeable heterogeneities. Geophysics, 76(4): F215-F226.
[18] Newman G A, Alumbaugh D L. 1995. Frequency-domain modelling of airborne electromagnetic responses using staggered finite differences. Geophysical Prospecting, 43(8): 1021-1042.
[19] Puzyrev V, Koldan J, Puente J, et al. 2013. A parallel finite-element method for three-dimensional controlled-source electromagnetic forward modelling. Geophys. J. Int., 193(2): 678-693.
[20] Saad Y. 2003. Iterative Methods for Sparse Linear Systems. 2nd ed. New York: Society for Industrial and Applied Mathematics.
[21] Schwarzbach C, Borner R U, Spitzer K. 2011. Three-dimensional adaptive higher order finite element simulation for geo-electromagnetics—a marine CSEM example. Geophys. J. Int., 187(1): 63-74.
[22] Shen J S. 2003. Modeling of 3-D electromagnetic responses in frequency domain by using staggered-grid finite difference method. Chinese J. Geophys. (in Chinese), 46(2): 281-289.
[23] Siripunvaraporn W. 1999. An efficient data-subspace two-dimensional magnetotelluric inversion and its application to high resolution profile across the San Andreas Faults at Parkfield, California. Corvallis: Oregon State University.
[24] 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.
[25] Streich R, Michael B. 2011. Electromagnetic fields generated by finite-length wire sources: comparison with point dipole solutions. Geophysical Prospecting, 59(2): 361-374.
[26] Tang J T, Ren Z Y, Hua X R. 2007. The forward modeling and inversion in geophysical electromagnetic field. Progress in Geophysics (in Chinese), 22(4): 1181-1194, doi: 10.3969/j.issn.1004-2903.2007.04.025.
[27] Ward S H, Hohmann G W. 1988. Electromagnetic theory for geophysical applications. // Nabighian M ed. Electromagnetic Methods in Applied Geophysics. SEG, 131-311.
[28] Weiss C J, Newman G A. 2002. Electromagnetic induction in a fully 3-D anisotropic earth. Geophysics, 67(4): 1104-1114.
[29] Weiss C J, Constable S C. 2006. Mapping thin resistors and hydrocarbons with marine EM methods, Part Ⅱ — Modeling and analysis in 3D. Geophysics, 71(6): G321-G332.
[30] 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 J. Geophys. (in Chinese), 55(4): 1390-1399, doi: 10.6038/j.issn.0001-5733.2012.04.035.
[31] Yee K S. 1966. Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media. IEEE Trans. Ant. Prop., 14(3): 302-307.
[32] Zaslavsky M, Druskin V, Davydycheva S, et al. 2011. Hybrid finite-difference integral equation solver for 3D frequency domain anisotropic electromagnetic problems. Geophysics, 76(2): F123-F137.
[33] 李建慧, 胡祥云, 曾思红等. 2013. 基于电场Helmholtz方程的回线源瞬变电磁法三维正演. 地球物理学报, 56(12): 4256-4267, doi: 10.6038/cjg20131228.
[34] 沈金松. 2003. 用交错网格有限差分法计算三维频率域电磁响应. 地球物理学报, 46(2): 281-289.
[35] 汤井田, 任政勇, 化希瑞. 2007. 地球物理学中的电磁场正演与反演. 地球物理学进展, 22(4): 1181-1194, doi: 10.3969/j.issn.1004-2903.2007.04.025.
[36] 杨波, 徐义贤, 何展翔等. 2012. 考虑海底地形的三维频率域可控源电磁响应有限体积法模拟. 地球物理学报, 55(4): 1390-1399, doi: 10.6038/j.issn.0001-5733.2012.04.035.