1 引言

海洋可控源电磁法(Marine controlled-source electromagnetic method, MCSEM)勘探是评价海底地层中是否含有油气的重要地球物理工具(Edwards, 2005)，它通过观测油气储层与低阻沉积围岩之间电阻率差异引起的电磁异常信号确定油气储层的分布范围，油气储层的电阻率通常是海底地层电阻率的数倍甚至几十倍，并借助电阻率与储层含油饱和度的密切关系直接用于探测地层的含油气性(Constable and Srnka, 2007; 何展翔等, 2009; Constable, 2010; Zhdanov, 2010).

油气勘探MCSEM方法通常采用船载可移动水平电偶极子源激发0.1Hz到5Hz的低频电磁信号，置于海底上方100 m范围内的接收器接收来自海底地层反射和折射的电磁场信号，数据处理解释方法通常是绘制发射频率对应的电磁场振幅和相位随收发距变化的MVO (Magnitude versus offset)曲线和PVO (Phase versus offset)曲线，其振幅和相位依赖于海底介质的电导率.目前，基于电导率各向同性MCSEM数值模拟方法研究成果较多(Mitsuhata, 2000; Li and Key, 2007; Key, 2009; 杨波等, 2012; Grayver et al., 2013; 杨军等, 2015)，然而，由于受应力、热对流、温差以及物质迁移等多种因素的影响，往往会引起地质体电导率各向异性，同时，随着海洋可控源电磁装备野外观测精度的提高，加上在实际数据处理解释中发现很多现象仅仅利用电导率各向同性理论难以作出合理的解释，电导率各向异性MCSEM正演研究逐渐引起了人们的重视.Everett和Constable (1999)实现了电导率单轴各向异性一维MCSEM电磁响应计算，得到了与各向同性介质解释方法不同的结果；Løseth和Ursin (2007)提出了电导率任意各向异性介质一维MCSEM电磁响应计算方法；Streich和Becken (2011)实现了不同发射源、电导率垂直各向异性(VTI)介质一维MCSEM正演计算方法；刘颖和李予国(2015)利用双重傅里叶变换和矢量势函数，推导出了层状电导率VTI介质中任意取向电偶源的电磁场计算公式；罗鸣和李予国(2015)在LØseth算法的基础上，实现了电导率任意各向异性一维层状介质海洋可控源电磁场计算方法.Kong等(2008)对基于二次场电导率各向异性二维海洋可控源电磁法有限元正演算法进行了研究；Wiik等(2013)实现了电导率VTI介质二维MCSEM积分方程数值模拟计算方法；Li等(2011, 2013)提出了基于二次场电导率各向异性介质MCSEM响应自适应有限元数值模拟方法.电导率各向异性介质三维海洋可控源电磁数值模拟方法也已取得重要进展(陈桂波等, 2009; Newman et al., 2010; Puzyrev et al., 2013; 殷长春等, 2014; 周建美等, 2014; 蔡红柱等, 2015).一般来说，基于积分方程的二维MCSEM数值模拟方法对复杂地质的模拟普遍存在一定难度；有限元法适用于模拟复杂电性结构，且计算精度较高，但基于二次场MCSEM有限元数值模拟方法需要额外计算一次场，计算量较大；而三维正演算法很难合理地平衡计算速度和精度，要达到较高的精度，往往要花费大量的计算时间；另外，这些一维、二维和三维MCSEM正演研究大多假定电导率具有垂直对称轴的横向同性介质(VTI介质)或电导率具有水平对称轴的横向同性介质(Azimuthal anisotropy介质)，这些介质都仅有一个旋转对称轴.

本文在前人研究成果的基础上，以电导率具有双轴对称性的正交各向介质为模型，研究基于总场电导率正交各向异性2.5维海洋可控源电磁有限元数值模拟方法.为了模拟起伏地形和地下复杂电性结构，采用任意四边形单元对研究区域剖分，在单元分析中对电磁场采用双二次插值，并采用伪δ函数等效源项提高数值模拟精度，通过与一维解析解、二维模型已有文献自适应非结构有限元计算结果对比，验证了方法的正确性，通过对二维地电模型的有限元计算，分析电导率各向异性对海洋可控源电磁响应的影响.

2 2.5维海洋可控源电磁问题的偏微分方程

如图 1a所示，二维地电模型由空气、海水、异常体和大地组成，设y方向与走向平行，并且该方向上的电性参数(包括电导率、介电常数和磁导率)不变，而在x-z平面内这些电性参数可以变化.假定模型中的磁导率和介电常数分别是自由空间磁导率μ和介电常数ε，任意方向的电性源置于海水的任意位置，电磁场是以e^iωt变化的谐变场，则二维地电模型中任意点的电场强度E和磁场强度H满足频率域Maxwell方程式：

(1)

(2)

式中, ω是角频率；i是虚数单位；J _s是源电流密度；是具有双轴对称性正交各向异性介质的电导率，是电阻率(x、y、z)中表示为

(3)

电导率各向同性介质(σ_x=σ_y=σ_z)、电导率具有垂直对称轴的横向同性介质(又称垂直各向异性介质或VTI介质，σ_x=σ_y≠σ_z)和电导率具有水平对称轴的横向同性介质(又称方位各向异性介质，σ_x≠σ_y=σ_z)是电导率正交各向异性介质的特例.

将(2)式代入(1)式并写成场分量的形式，可以得到关于E、H及其导数的方程，有

(4)

(5)

(6)

(7)

(8)

(9)

图 1a电导率分布是二维的，而发射源、电磁场是三维的，构成二维模型三维源、三维电磁场即2.5维MCSEM电磁响应问题，通过将(4)式-(9)式中E和H的每一个分量及源项，沿构造走向y方向作一维Fourier变换，从空间域转换到波数域，就能将三维的微分方程(1)式和(2)式转换为二维的情况.一维Fourier变换的表达式为

(10)

由Fourier变换的导数特性以及(10)式，有

(11)

(12)

(13)

其中，代表波数域中的傅氏场变量 _x、_y、_z、 _x、_y和_z；k_y为走向方向上的波数；F代表空间域中的场变量E_x、E_y、E_z、H_x、H_y和H_z.

利用(11)式-(13)式，空间域三维问题(4)式-(9)式转换得到以下波数域(k_y域)的二维微分方程：

(14)

(15)

(16)

(17)

(18)

(19)

其中，_sx、_sy和_sz代表波数域(ω，k_y)中场源电流密度沿直角坐标系坐标轴分量.

将(16)式和(17)式联立，有

(20)

(21)

其中，C=iωμ，Q_x=(σ_x+iωε)，k_ex²=k_y²+CQ_x.

将(14)式和(19)式联立，有

(22)

(23)

其中，Q_z=(σ_z+iωε)，k_ez²=k_y²+CQ_z.

将(21)式和(23)式代入(18)式，有

(24)

其中，Q_y=(σ_y+iωε).

将(20)式和(22)式代入(15)式，有

(25)

式(24)、式(25)是电导率正交各向异性2.5维海洋可控源电磁问题的偏微分方程，它是由_y和_y控制的耦合方程组，求解这个耦合方程组就可以获得波数域(ω, k_y)中的_y和_y.

3 偏微分方程相应的有限元方程

根据伽里金(Galerkin)加权余量法中余量的概念(徐世浙，1994; Jin, 2002)，由(24)式可得到电场_y控制方程的余量R^e为

(26)

由Galerkin加权余量法原理，要迫使余量R^e在包含N_e个单元的整个区域里面的加权积分为零，选择形函数N^e作为积分的权函数，有

(27)

式中，Ω_e表示第e个单元的求解区域，N_e表示全区单元总数，N^e表示第e个单元的形函数.

将(26)式代入(27)式，并利用下面的二维Green公式对(27)式中的被积函数做降阶处理，

(28)

(29)

式中，n_x、n_z表示曲面的外法线与坐标轴之间夹角的三角函数.

假设外法线与坐标轴x之间的夹角为α，有n_x=cosα、n_z=sinα，假设磁场的切向分量为_l，有_l=-_zcosα+_xsinα，则(27)式最终化简为

(30)

同理，由(25)式，有

(31)

(30)式和(31)式右端项中的曲线积分，是关于各个子单元区域边界Ω_e的积分.由于_l和_l在单元边界上是连续的，因此，此曲线积分在单元公共边界上相互抵消，只有在整个求解区域的边界上有值.为了简化计算，采用狄利克莱(Dirichlet)边界条件使得_y和_y在边界上都等于零，消去(30)式、(31)式中的环路积分项，这样，就获得2.5维海洋可控源电磁问题的偏微分方程相应的有限元方程.

4 等参有限元法

求解2.5维海洋可控源电磁问题有限元方程步骤如下.

4.1 网格剖分

首先对研究区域进行二维网格剖分.为了能更好地模拟起伏地形和地下复杂电性结构，采用任意四边形单元对研究区域进行剖分，如图 1b所示，研究区域分成两部分，一部分是目标区域，另一部分是网格外延区域(网格外延区域主要是为了模拟无穷边界)，目标区域和网格外延区域均采用不等距网格，目标区域为地质体、地形和发射源赋存区域，也是数据采集区域.

为了便于计算，对剖分网格中的任意四边形单元进行等参变换，将直角坐标系内形成的不规则四边形网格单元(如图 2a所示的子单元)转换成自然坐标系下边长为2的正方形网格单元(如图 2b所示的母单元)，并定义正方形单元的中点为ξ、η坐标轴的原点，直角坐标系下任意四边形单元的八个节点分别与自然坐标系下正方形单元的八个节点对应，这样，子单元中的面积元dxdz与母单元之中的面积元dξdη的关系为

(32)

式中，| |是雅克比变换行列式.

4.2 双二次插值

在任意单元e中，对_y和_y采用双二次插值，则在每个单元内的坐标x、坐标z、_y和_y满足如下的表达式：

(33)

式中，j=1, 2, …, 8；x_j、z_j是子单元八个节点的坐标值； _yj、_yj是子单元八个节点上的电磁场值；N_j^e(ξ, η)为母单元八个节点在自然坐标系ξ、η下的双二次插值形函数，且有

(34)

(35)

(36)

式中，

4.3 单元积分

将单元离散变量(33)式代入(30)式可得：

(30)式左边第一项单元积分：

(37)

(30)式左边第二项单元积分：

(38)

(30)式左边第三项单元积分：

(39)

(30)式左边第四项单元积分：

(40)

(30)式左边第五项单元积分：

(41)

4.4 源项加载

为了消除电偶极子源处的奇异性，传统MCSEM有限元计算是通过二次场方法的计算实现的，将总电磁场值分为一次场和二次场，一次场选用均匀介质或一维层状介质的解析解，通过有限元方法求解二次场，这种计算策略需要额外计算一次场，计算量较大，另外，地下介质的电性结构复杂难以确定一次场计算的电性参数也让这种方法的适用性变差.本文采用具有一定面积的伪δ(r)函数表达源电流分布，避免源处电磁场的奇异性，从而直接求解总电磁场(Mitsuhata, 2000；沈金松和孙文博，2008；范翠松，2013；张继锋等，2013).伪δ(r)函数的表达式为

(42)

(42)式中参数τ是控制源分布宽度和幅值的参数.实际上伪δ(r)函数起到了一个低通滤波的作用，τ取值越大，源分布就越宽，消除奇异能力越好，等效于偶极子的作用越差，而τ取值越小，源分布就越窄，消除奇异能力越差，等效于偶极子的作用越好.所以伪δ函数必须选取合适的τ值来等效偶极子作用，才能达到较好的效果.在文中作者采用源附近对称的16个单元等间距剖分，其τ值取网格间距.这样既能保证源项非零元素的个数，又能使二阶插值函数适应源的分布，从而提高计算精度.

电偶极子处于坐标系不同方向时源项的加入也有所差异.当电偶极子平行于x轴方向时，在单元e内的离散形式如下：

(43)

当电偶极子平行于y轴方向时，在单元e内的离散形式如下：

(44)

当电偶极子平行于z轴方向时，在单元e内的离散形式如下：

(45)

式中，x₀、z₀为电偶极子源的中心坐标，I是电流强度，dl是电偶极子强度.

则有，

(30)式右边第一项单元积分：

(46)

(30)式右边第二项单元积分：

(47)

(30)式右边第三项单元积分：

(48)

将(37)式-(41)式、(46)式-(48)式代入(30)式，针对任意剖分单元e，有

(49)

同理，任意剖分单元e内(31)式离散为

(50)

式(49)和(50)为电导率正交各向异性2.5维海洋可控源电磁问题等参有限元法任意剖分单元e离散方程，方程中形如类型的二重积分，其被积函数f (ξ, η)是与四边形形态有关的复杂函数，不存在解析解，采用高斯数值积分计算(徐世浙, 1994).

4.5 方程组求解

按照剖分节点的总体序号，将各单元离散后的方程组(49)式和(50)式单元系数矩阵中的各元素进行扩展总体合成后得到总体系数矩阵，联合右端源项，得到大型复系数线性方程组：

(51)

其中，A是大小为2N的待求列向量(N为所有剖分单元节点总数)，用于存放各节点处的_y和_y；B为源项；K是大小为2N×2N的复系数稀疏矩阵，采用非零元素压缩存储.

代入Dirichlet边界条件，采用不完全Cholesky分解-双复共轭梯度法求解(张永杰和孙秦，2007)线性方程组(47)式，可以获得各个节点上的_y和_y；然后根据辅助场的表达式(20)-(23)获得波数域中的_x、_z、_x和_z；利用反傅里叶变换可以得到空间域电磁场值(底青云等，2004；张继锋等，2013).

在求解过程中，波数k_y选取是影响计算精度和计算效率的重要因素，由于各种地电断面的场函数不同，满足一种地电断面的波数组不一定适用于另一种地电断面；另外，增加一个波数求解线性方程组的次数增加，必然使得计算时间增加，倘若减少波数个数，虽然计算速度较快，但难以保证精度，许多学者均针对这些问题做了大量研究(Mitsuhata, 2000, 沈金松和孙文博，2008；张继锋等，2013).根据研究区域大小和网格剖分大小，确定最优的k_y值范围，同时，考虑到波数域电磁场的变换特征，在k_y低值端要适当加密取值，高值端可以抽稀，适当选取k_y的个数，以减少计算量并保证获得高精度计算结果的原则，本文采用15个最优化波数k_y(×10^-5)：1.78、3.16、5.63、10、17.8、31.6、56.3、80、100、178、316、563、1000、2000、5000(底青云等，2004).

上述理论是建立在电偶极源与坐标轴平行的基础上，然而，实际海洋电磁观测中由于海底洋流作用及拖船航线的影响，发射偶极源可能产生摆动、倾斜和旋转而为任意取向的电偶极源，可以根据坐标旋转公式将任意取向电偶极源转换到观测坐标系中，分解为三个相互正交等效电偶极源_sx、_sy和_sz，分别计算三个电偶极源的电磁响应，然后矢量相加即可得到任意取向电偶源的电磁场.另外，当发射源为有限长的线源激发时，可将线源激发下的电磁场值由多个电偶极源下的场值线性叠加.

5 模型算例

为了检验算法的正确性，我们设计了电导率各向异性一维、二维模型；然后对二维地电模型不同各向异性系数的情况进行正演计算，分析电阻率各向异性系数对海洋可控源电磁响应的影响；最后，建立海底起伏地形的油藏模型，考察电导率各向异性对地形引起的异常和油藏引起的异常的影响.在数值模拟中，考虑电阻率各向同性空气层的影响，设定电阻率为10¹²Ωm；x、z方向采用非均匀网格剖分，网格大小130×50.

5.1 电阻率各向异性一维模型算例

层状模型选用中间层为电阻率各向异性且相对高阻的油气薄层模型，如图 3所示，海水层厚度为950 m，电阻率ρ_x1=ρ_y1=ρ_z1=0.25 Ωm；覆盖层厚度700 m，电阻率ρ_x2=ρ_y2=ρ_z2=1 Ωm；高阻油气薄层厚度600 m，电阻率ρ_x3=ρ_y3=50 Ωm、ρ_z3=200 Ωm；基底为沉积层，其电阻率ρ_x4=ρ_y4=ρ_z4=1 Ωm；水平电偶极子Tx位于海底上方50 m，中心坐标为(0 m, 0 m, 900 m)，发射频率均为0.25 Hz，采用两种装置：(1) Inline装置，即单位电偶极子沿着x轴方向，且测线方向与x方向一致；(2) Broadside装置，即单位电偶极子沿着y轴方向，且测线方向与x方向一致.位于海底上方25 m的电磁采集站Re平行x轴在-8000m < x < 8000 m范围内按不同间距布设(Re的z坐标均为925 m).

图 4给出了利用本文2.5维电导率正交各向异性海洋可控源电磁有限元计算的Inline装置E_x响应和H_y响应结果(2.5D, anisotropic)与解析解结果(1D, anisotropic)对比曲线，图 5给出了利用本文2.5维电导率正交各向异性海洋可控源电磁有限元计算的Broadside装置E_y响应和H_x响应结果(2.5D, anisotropic)与解析解结果(1D, anisotropic)对比曲线.由图 4、图 5可见，2.5维海洋可控源电磁MVO曲线和PVO曲线变化稳定，与解析解曲线吻合非常好，说明我们的算法具有较高的计算精度.为进一步分析高阻油气层电导率各向异性对MCSEM电磁响应的影响，图 4、图 5给出了高阻油气层电阻率为各向同性(ρ_x3=ρ_y3=ρ_z3=50 Ωm)时的电磁响应(1D, isotropic)，由图可见，随着收发距的增大，电阻率各向异性的MVO曲线、PVO曲线与电阻率各向同性的MVO曲线、PVO曲线分离，并且同一个收发距，电磁场振幅随着ρ_z3的增大而增大，而电磁场相位随着ρ_z3的增大而减小.将油气层电阻率各向异性与油气层电阻率各向同性时的电磁场振幅归一化(normalized responses)和相位归一化(phase difference)，曲线形态更加直观地显示了高阻油气层电导率各向异性对MCSEM电磁响应的影响.

5.2 电导率各向异性二维模型算例

为了进一步验证2.5D电导率正交各向异性MCSEM有限元数值模拟计算程序的正确性，设计了如图 6所示的二维地电模型，海水的电阻率是各向同性的，电阻率ρ_x1=ρ_y1=ρ_z1=0.3 Ωm，深度为1 km；顶界面距海底1000 m处有一个6000 m×100 m的高阻油气薄层，中心坐标为(0 m, 0 m, 2050m)，其电阻率是各向同性的，电阻率ρ_x2=ρ_y2=ρ_z2=50 Ωm；高阻油气薄层周围沉积物的电阻率是各向异性的，电阻率ρ_x3=ρ_y3=1 Ωm、ρ_z3=4 Ωm；沿x方向的水平单位电偶极子源发射频率为0.25 Hz，中心坐标为(0 m, 0 m, 950 m)；沿x方向布置的接收站置于海底.图 7是位于海底本文2.5维有限元数值模拟结果(Numerical solution from 2D IFEM)与海洋电磁二维自适应非结构有限元数值模拟结果(Numerical solution from 2D AFEM)的对比，海洋电磁二维自适应非结构有限元数值模拟结果采用加州大学圣地亚哥分校Scripps海洋研究所海洋电磁实验室的计算结果，由中国海洋大学罗鸣博士提供.从图 7对比曲线可以看出，本文有限元数值模拟程序计算结果和二维自适应非结构有限元计算结果十分吻合，这也进一步验证了基于总场电导率正交各向异性2.5维海洋可控源电磁等参有限元数值模拟计算程序的正确性.

海洋可控源电磁勘探常常采用多发射源的形式，接收点数量相对较少且放置在目标区域上方，对于图 6模型，令发射源沿x方向在-2000~2000 m之间按400 m间距移动，其余参数不变.此时为2.5维电导率正交各向异性多发射源海洋可控源电磁正演问题，取计算区域对所有发射源均满足截断边界条件，则多个发射源对应的控制方程的系数矩阵保持不变，只需要改变控制方程的右端项，其计算精度与单个发射源的计算精度相同，相应的数值模拟结果如图 8所示.

5.3 各向异性系数变化对海洋可控源电磁响应的影响

为了分析电导率正交各向异性对海洋可控源电磁响应的影响，用和分别表示垂直各向异性系数和水平各向异性系数，在中间层为电阻率为各向异性且相对高阻的油气薄层模型(如图 3所示)基础上，我们设计了如图 9所示的二维地电模型.深度为950 m的海水层电阻率是各向同性的，电阻率ρ_x1=ρ_y1=ρ_z1=0.3 Ωm.覆盖层下方的基岩电导率是各向同性的，其电阻率ρ_x2=ρ_y2=ρ_z2=1 Ωm.覆盖层下方的右侧基岩中存在电阻率各向同性且相对高阻的形态不规则油气层(resistor)，其电阻率ρ_x3=ρ_y3=ρ_z3=100 Ωm.假设测线方向的水平电偶极子源Tx位于海底正上方50 m处，中心坐标为(0 m, 0 m, 900 m)，发射频率为0.25 Hz.位于海底上方25 m的电磁采集站Re平行x轴在-8000 m < x < 8000 m范围内按不同间距布设(Re的z坐标均为925 m).覆盖层(cover rock)的厚度为700 m，其电阻率正交各向异性，为.

假定覆盖层水平各向异性系数λ²=0.5，垂直各向异性系数η²分别取值为1.0、2.5、5.0、7.5和10.0，相应的电阻率参数值ρ_x、ρ_y和ρ_z见表 1.图 10为数值模拟结果，同时，为了考察高阻油气层的影响，给出了η²=1.0且不存在油气层的数值模拟结果(η²=1.0, without resistor).从图 10a的MVO曲线和图 10b的PVO曲线看，电场振幅随着收发距的增大而衰减，大收发距时电场相位随着收发距的增大而减小，这主要是电磁能量向海水和海底地层扩散引起，而中间小收发距(发射源附近)时，电场相位随着收发距的增大而增大，随后随着收发距的增大而减小，主要是由于小收发距时，电磁场仅与海水的电导率有关，受海底地层介质的电导率影响甚小；电场振幅随着垂直各向异性系数η²的增大而增大，电场相位随着垂直各向异性系数η²的增大而减小，存在高阻含油气层的电场振幅高于不存在高阻含油气层的电场振幅，存在高阻含油气层的电场相位小于不存在高阻含油气层的电场相位，这主要是低阻层对电磁场能量的强衰减作用引起，垂直各向异性系数η²的减小和不存在异常体相当于平均电阻率降低.

为了更直观地显现电导率垂直各向异性系数和高阻油气层对MCSEM电磁响应的影响，图 10c和图 10d分别给出了电导率垂直各向异性系数变化的电场振幅归一化曲线和电场相位归一化曲线，其中电场振幅归一化响应用| E_i|/|E₀ |计算，电场相位归一响应化用φ_i-φ₀计算，|E_i|、φ_i分别代表覆盖层电导率垂直各向系数η²变化时的电场振幅和相位，| E₀|和φ₀分别代表覆盖层电导率垂直各向系数η²=1且不存在油气层resistor时的电场振幅和相位.从曲线可以看出，电导率垂直各向异性系数对海洋可控源电磁数据产生严重的影响，左侧虽然没有油气层存在，但由于电导率各向异性的影响出现明显的电场归一化异常和相位差异常.

同样，假定覆盖层垂直各向异性系数η²=1，水平各向异性系数λ²分别取值为1.0、2.0、5.0、和10.0，相应的电阻率参数值ρ_x、ρ_y和ρ_z见表 2.图 11所示为数值模拟结果，由图可以看出，电场振幅随着水平各向异性系数λ²的增大而增大，电场相位随着水平各向异性系数λ²的增大而减小，曲线右侧明显地显现了高阻油气层的影响，同时，由于电导率各向异性的影响，在图 11c和图 11d的左侧出现明显的电场归一化异常和相位差异常.与图 10相比较可见，电导率水平各向异性系数对MCSEM电磁响应的影响小于电导率垂直各向异性系数对MCSEM电磁响应的影响.

5.4 起伏地形模型

为了了解起伏地形条件下电导率正交各向异性对MCSEM勘探实测数据的影响，设计二维起伏海底地形含油气构造模型，如图 12a所示.深度为950 m的海水电阻率是各向同性的，电阻率为0.25 Ωm，水平电偶极子源Tx位于海底正上方50 m处，中心坐标为(0 m, 0 m, 900 m)，发射频率为0.25 Hz，位于海底上方25 m的电磁采集站Re平行x轴在-8000 m < x < 8000 m范围内按不同间距布设.发射电偶极子左侧和右侧分别有一个凸地形和凹地形，右侧距离海底700 m存在1000 m×600 m大小的高阻油气藏，高阻油气藏和海底的电导率呈正交各向异性，其海底的电导率为ρ_x2、ρ_y2和ρ_z2，而高阻油气藏的电导率分别为ρ_x3、ρ_y3和ρ_z3.该模型按MOD1、MOD2、MOD3和MOD4四种情况分别计算，如图 12b所示，MOD1为无高阻油气藏、海底电导率为各向同性(ρ_x2=ρ_y2=ρ_z2=1 Ωm)的水平层状模型，MOD2为存在电导率各向同性高阻油气藏(ρ_x3=ρ_y3=ρ_z3=100 Ωm)、海底电导率为各向同性和存在起伏地形模型，MOD3为存在电导率各向同性高阻油气藏(ρ_x3=ρ_y3=ρ_z3=100 Ωm)、海底电导率为各向异性(ρ_x2=0.5 Ωm、ρ_y2=1 Ωm、ρ_z2=1.2 Ωm)和存在起伏地形模型，MOD4为存在电导率各向异性高阻油气藏(ρ_x3=100 Ωm、ρ_y3=100 Ωm、ρ_z3=100 Ωm)、海底电导率为各向同性和存在起伏地形模型.

图 13a和图 13b分别为MOD1、MOD2、MOD3和MOD4四种不同情况电磁采集站电场的振幅和相位响应.图 13c和图 13d分别为MOD1、MOD2、MOD3和MOD4四种不同情况的电场振幅归一化| E′|/|E₁|响应曲线和相位归一化φ′-φ₁响应曲线，其中| E′|、φ′分别代表MOD1、MOD2、MOD3和MOD4的电场振幅和相位，| E₁|和φ′分别代表MOD1的电场振幅和相位.从图 13c和图 13d可以直观地看出，相对于MOD1，MOD2在凸地形上方出现高振幅、高相位异常，在高阻油气藏上方也出现高振幅、高相位异常，而在凹地形上方出现低振幅、低相位异常，这主要是因为凸地形相当于高阻体，凹地形属于低阻体，而地质体低电阻率会对电磁场能量产生强衰减作用；与MOD2的MVO曲线和PVO曲线类比，由于背景电阻率各向异性的影响，MOD3在凸地形上方出现的高振幅异常更强，在凹地形上方也出现了高振幅异常(将凹地形引起的实际低振幅异常掩盖)，因此在野外实际勘探中应注重凹地形下电导率各向异性的影响，另外，相对于MOD1，MOD3左侧不存在油气藏，但由于背景电导率各向异性的影响出现假异常；与MOD2类比，由于高阻油气藏电导率各向异性的影响，MOD4油气藏上方的异常被消弱，因此，油气藏电导率各向异性对MVO曲线和PVO曲线的影响要引起重视，否则会得到错误的地电模型.

6 结论

针对电导率具有双轴对称性的正交各向异性介质模型，本文提出了基于总场电导率正交各向异性2.5维海洋可控源电磁等参有限元数值模拟方法.从麦克斯韦方程组出发，推导出了电导率正交各向异性波数域电磁场耦合微分方程；利用伽里金加权余量法给出了相应的有限元方程，采用任意四边形单元对研究区域进行剖分，以模拟起伏地形和地下复杂电性结构；在单元分析中对电磁场采用双二次插值，并采用伪δ函数等效源项提高数值模拟精度；采用等参变换解决任意四边形单元积分计算问题.与一维模型解析解以及已有二维模型自适应非结构有限元计算结果对比验证了方法的正确性.进一步在模拟实际勘探中电导率正交各向异性地下任意形状的异常体以及地形模型的等参有限元计算中，也取得了可靠的结果.

水平海底二维地电模型的数值结果表明，电场振幅随着收发距的增大而衰减，电场振幅随着垂直各向异性系数和水平各向异性系数的增大而增大，电场相位随着垂直各向异性系数和水平各向异性系数的增大而减小，电导率水平各向异性系数对MCSEM电磁响应的影响小于电导率垂直各向异性系数对MCSEM电磁响应的影响.海底起伏地形二维地电模型的数值结果表明，电导率各向异性会引起假异常，同时会消弱甚至掩盖凹地形和高阻油气藏引起的MVO异常和PVO异常，在实际资料的处理解释中应引起重视.

致谢

中国海洋大学罗鸣博士提供了文中二维模型海洋可控源电磁自适应非结构有限元数值模拟结果，中国地质科学院地球物理地球化学勘查研究所方慧博士审阅了修改稿，两位审稿专家提出了宝贵的修改意见和建议，在此一并表示感谢.