1 引言

近十年来，海洋可控源电磁法(MCSEM)已经成为实现寻找油气田、研究海洋深部构造等勘探目的一种行之有效的勘探方法(Constable 2007，2010；Weitemeyer 2010).由于数据解释的需要，海洋可控源三维电磁数值模拟技术一直是国内外学者研究的热点.一般而言，对于海洋可控源电磁法正演研究有两种方式.一种是直接求解电磁场各分量，另外一种是通过求解电磁场所满足的位场方程，将矢量位、标量位求解出来后，再通过对位求导，便可得到相应的电磁场各分量.在直接求解电磁场方面，Li和Key(2007)开展了海洋可控源电磁法二维自适应非结构化有限单元数值模拟，取得了较好的计算结果；Key和Ovall(2011)实现了海洋电磁法2.5维自适应有限元并行计算研究.Schwarzbach等(2011)采用高阶矢量有限元法实现了电导率各向同性介质的MCSEM三维数值模拟.Cai(2014)应用六面体网格单元进行了MCSEM三维矢量有限元正演研究.陈桂波(2009)应用压缩映射的积分方程法，实现了各向异性介质的MCSEM响应研究.殷长春等(2014)进行了三维任意各向异性介质海洋可控源电磁有限差分数值模拟.杨波等(2014)采用规则网格有限体积法开展了考虑海底地形的海洋可控源电磁响应特征研究.杨军等(2015)实现了海洋可控源电磁三维非结构化矢量有限元数值模拟；在通过计算矢量位、标量位的方式来间接求解电磁场研究方面，汤井田等(2007)详细地给出了Coulomb规范下磁矢量势-电标量势的hp型自适应有限元分析的基本理论，为提高电磁场数值模拟的精度提供了一种新的途径. Pyzyrev等(2013)以海水全空间为背景层二次场算法实现了各向异性介质海洋可控源电磁三维正演，但其仅考虑三主轴电导率各向异性介质.Ansari(2014)实现了电导率各向同性介质的可控源电磁法三维矢量有限元数值模拟.Jahandari和Farquharson(2015)基于非结构化网格开展了各向同性介质的可控源电磁法三维有限体积正演研究.蔡红柱(2015)在Puzyrev等(2013)的基础上，同样仅考虑电导率三主轴各向异性的条件下，以水平层状介质为背景层实现了电导率各向异性的海洋可控源三维有限元数值模拟.

由地质运动造成的地层的隆起、翻转等引起的电导率各向异性往往并不局限于三主轴，甚至可能更为复杂，并且在实际探测中，复杂的海底地形与电导率各向异性的影响对实际数据的解释可能会产生一定的偏差.因此，研究复杂地形及电导率任意各向异性对于海洋可控源电磁法探测的影响是十分有必要的.相比于规则网格的积分方程法和有限差分法，非结构化网格的有限单元法更有利于精确模拟复杂地质模型.有限单元法又分为矢量有限元法和节点有限元，相比于节点有限元，矢量有限元形成的有限元方程过程相对复杂，且其形成的系数矩阵的条件数往往比较大，采用迭代算法求解收敛比较慢(Jin, 2002).由于对Maxwell方程进行直接求解，会出现解不稳定的问题，并且一次场/二次场分离算法与直接计算总场算法相比具有更高的计算精度和更好的稳定性，同时，若将复杂的海底地形视为异常体，那么一次场/二次场分离算法便同样适用于复杂的海底地形条件下的电磁场数值模拟.因此，本文采用一次场/二次场分离算法对海洋可控源电磁响应进行模拟.使用Galerkin加权余量法推导了Coulomb规范下的二次场磁矢量位、标量位在电导率任意各向异性条件下所满足节点有限元方程.同时，采用Krylov子空间迭代法中的IDR(n)迭代算法(Sonneveld et al., 2008)结合不完全LU分解预处理方法求解有限元方程组.并设计两个典型的地电模型对本文算法的正确性和有效性进行了验证，并将IDR(s)迭代算法与Krylov子空间迭代法中的QMR、BICGSTAB算法效率进行对比.

2 正演理论

地球物理勘探中，往往是采用低频电磁信号，忽略位移电流, 采用正弦谐变时间因子e^-iwt，Maxwell′s方程可表示如下形式:

(1)

(2)

其中，w为角频率，μ₀为真空中的磁导率.J_s为场源电流分布. 为任意各向异性电导率张量：

(3)

为了计算电导率张量，通常事先设定一个参考电导率σ_c(其中三个对角线的元素表示三个主轴的电导率)：

(4)

正如图 1所示，接着通过三次欧拉旋转，便可得到任意各向异性的电导率张量：

(5)

其中旋转矩阵为

引入基于Coulomb规范下的磁矢量位A、标量位Ψ来表示电场、磁场，如下所示：

(6)

(7)

将方程(6)、(7) 代入方程(1)、(2) 得：

(8)

(9)

在二次场算法中，总磁矢量位和标量为可分解成二次场与背景场之和：

(10)

(11)

(12)

将方程(10)—(12) 代入方程(8)、(9) 可得关于电磁场二次位表达式：

(13)

(14)

通过观察方程(13)、(14) 可知，方程右端项是关于一次位，而不含电流项.因此可避免源点的奇异性的影响.而一次位可通过层状介质或者均匀半空间为背景层求得.考虑如下关系：

(15)

方程(13)、(14) 可改写成如下形式：

(16)

(17)

这样就避免求解一次场标量位的梯度，提高一次场的精度.

3 有限单元分析

将方程(16)、(17) 按照x, y, z轴依次展开，可得如下方程：

(18)

(19)

(20)

(21)

利用伽辽金方法对方程(18)—(21) 进行加权，结合矢量恒等式以及散度定理，最后可得关于二次位的体积分方程组：

(22)

(23)

(24)

(25)

本文研究中，将研究区域剖分成有限多个四面体单元e，在各个单元中，电导率是固定不变的，单元中的二次矢量位和标量位呈线性变化：

(26)

线性插值函数定义如下所示(Jin, 2002)：

(27)

其中，V^e为单元体积.

将(26)、(27) 方程代入(22)—(25) 式最后可得离散方程组：

(28)

其中，

为了求解方程(29)，还需要加入相应的边界条件.为了简便起见，本文采用狄利克雷边界条件：

(29)

4 有限元方程组的求解

对于所形成的有限元线性方程组，一般有两种求解方法.一种是直接解法，另一种是迭代求解方法.与迭代法相比，直接法求解大型方程组具有精度高、稳定性好的优点，但其内存需求大、求解效率低，在一般的微机上实现较为困难.因此，对于大型、稀疏的线性方程组的求解，大多是采用Krylov子空间迭代法.本文引入Krylov子空间迭代法中的IDR(s)算法系统求解方程组，IDR(S)是近几年来才发展起来的比较有效的一种短递归求解线性方程组算法，该方法可以使得迭代求解残差限制在不断减小的子空间中，相对于常用的QMR、BI-CGSTAB迭代算法在每次迭代过程中需要计算2N次矩阵向量相乘，IDR(s)只需要计算N+N/S次矩阵向量相乘(其中N为矩阵的大小，s为固定子空间的维数，可取1~8)，因此求解效率可大幅度提升.特别是当系数矩阵阶数特别大时，IDR(s)优势更为明显.针对系数矩阵的稀疏性特点，本文将不完全LU分解预处理算法与IDR(s)相结合来加快求解有限元线性方程组的收敛速度.为了进一步说明IDR(s)算法，这里给出IDR(s)法具体实现过程：

给定系数矩阵A ∈ С^N×N，x₀，b ∈ С^N×1，P ∈ С^N×s，TOL ∈ (0, 1)，MAXIT>0，

计算r₀=b-Ax₀，

  for n=0 : s-1 do

   v =Ar_n，ω=v^Hr_n/v^Hv，Δx_n=ωr_n，

   Δr_n=-ωv，

   r_n+1=r_n+Δr_n，x_n+1=x_n+Δx_n，

  end for

ΔR_n+1=(Δr_n…r₀)，ΔX_n+1=(Δx_n…x₀)，

n=s，

While ‖ r_n ‖>TOL and n<MAXIT

  for k=0 : s

   Solve c from P^HΔR_nc=P^Hr_n，

   v =r_n-ΔR_nc，

   if k=0 then

    t=Av，ω=(t^Hv)/(t^Ht)，

    Δr_n=-ΔR_nc- ω t，

    Δx_n=-ΔX_nc+ ω v，

   else

    Δx_n=-ΔX_nc+ ω v，

    Δr_n=-AΔx_n，

   end if

    r_n+1=r_n+Δr_n，x_n+1=x_n+Δx_n，

    n=n+1，

    ΔR_n=(Δr_n-1…r_n-s)，

    ΔX_n+1=(Δx_n-1…x_n-s).

  end for

end while.

5 加权移动最小二乘法求导

当求解方程(28) 结束后得到二次场矢量位A_s、标量位Ψ_s.可通过以下方程求解相应的电磁场各分量：

(30)

(31)

通过观察方程(30)、(31) 可以发现，要求得电磁场各分量，必须要对二次磁矢量位A_s、标量位Ψ_s进行求导.然而，对于完全非结构化网格节点上的二次位，若采用常规的求导方法对其进行求导数，往往会引起较大的误差.因此为了提高计算精度，文中采用指数函数加权的移动最小二乘法对二次位进行求导(Tabbara et al., 1994).具体过程如下：

首先设单元中的二次位的线性描述为

(32)

显而易见，系数T₁、T₂、T₃即为待求解的二次位的导数.因篇幅有限，在此直接给出指数函数加权的移动平均最小二乘法所形成的方程：

(33)

这里：

权函数W形式如下：

这里β为权函数系数，一般取1.F为待求二次位导数点附近N个节点的二次位值.

6 算例 6.1 算法验证

为了验证文中算法的正确性.设计一个三层水平层状沉积模型，如图 2a所示.采用海水-沉积层两层介质为背景层.设上半空间中的空气电导率为10^-12S·m^-1.海水电导率为3.2 S·m^-1，深度为0.3km.距离海水底面深度1 km有一水平方向无限延伸的各向异性高阻体，其厚度为100 m，参考电导率为σ_c=diag(0.01, 0.01, 0.025) S·m^-1，向y轴旋转30°.沉积层的电导率为1 S·m^-1.采用沿着x方向的水平电偶极子为激发源，其坐标位于(0, 0, 970).发射频率为0.25 Hz.偶极矩为1.沿着z=970处进行观测.网格剖分如图 2b所示，剖分区域大小(-4 km≤x, y≤4 km, -1 km≤z≤5 km)，网格单元总数1467492.在AMD Athlon(tm)×463 quad-core processor 2.6 GHz, 内存12 G计算平台上执行.将计算结果与一维拟解析解(Løseth and Ursim, 2007)进行对比.从图 3可知，本文计算的数值解与准解析解基本一致，验证了本文算法的正确性.图 4、表 1为IDR迭代算法与多种迭代算法在不同预处理条件下的计算效率对比，可以发现，IDR迭代算法在迭代次数与其余两种算法相近的情况下，计算时间明显减少.

6.2 沉积层电导率沿着y轴/z轴旋转0°、30°、45°、60°

为了进一步验证文中算法的有效性，设计一个沉积层电导率任意各向异性海洋地质模型.如图 5所示.上半空间空气电导率设为10×^-12，海水深度为1 km，电导率为3.3 S·m^-1.沉积层中有一高阻体，顶面距离海水底面1 km，大小为2 km×2 km×0.1 km，中心坐标为(2000, 0, 2050)，电导率为0.01 S·m^-1.沉积层的参考电导率设为σ_c=diag(2, 1, 1).激发源为沿着x方向的水平电偶极子，发射频率为0.25 Hz，偶极距为100，源点坐标(0, 0, 950).接收机位于海水深度z=950处.剖分区域大小(-4 km≤x, y≤4 km, -1 km≤z≤4 km).网格单元总数为2930304.令沉积层参考电导率分别沿着y轴、z轴旋转0°、30°、45°、60°后研究电场三分量幅值的异常响应特征.

从图 6可知，参考电导率张量沿着y轴旋转的角度为30°、45°时，电场E_x分量的幅值在x、y方向都有所延伸，特别在y方向的延展更为明显；而电场E_y分量正好相反.相比于E_x、E_y分量，E_z分量幅值变化更明显，由一开始的对称到明显的不对称，且右侧的响应值比左侧的大.当参考电导率张量沿着y轴旋转60°时，电场E_x分量幅值分布有所收缩；电场E_y分量变化不大；而电场E_z幅值分布明显向外扩展，且左侧的幅值变大，但右侧幅值依然比左侧大.通过观察图 7，可知，参考电导率沿着z轴旋转30°、45°时，电场E_x分量分布发生倾斜向外扩张.E_y分量的幅值随着旋转的角度增大而有所增大，且在x、y方向分布也有所延伸；E_z分量幅值分布发生明显的倾斜变化.当参考电导率沿着z轴旋转60°时，E_x分量分布进一步发生扩展；E_y分量的幅值分布有所收缩；E_z分量幅值分布变成了斜对称.因此，可知不同方向的电导率各向异性所引起的电场异常响应与三主轴电导率各向异性的响应特征有明显的区别.

7 结论

本文提出了一种基于Coulomb规范下的磁矢量位、电标量位的电导率呈任意各向异性条件下的海洋可控源电磁正演算法，通过对典型的地电模型的试算和分析，得到了如下结论：

(1) 海洋层状沉积模型理论计算结果表明：本文提出的电导率呈任意各向异性条件下的海洋可控源电磁有限元数值模拟算法是正确的，并且具有较高的计算精度；该算法是基于二次场矢量位、标量位以及非结构化网格剖分技术，既可以有效避免源点奇异性的影响，又能满足节点有限元对于场的连续性要求，因此具有良好的通用性，可用于几何复杂、电导率任意各向异性介质的航空电磁、井中电磁等电磁数值模拟研究.

(2) 引入IDR(s)迭代求解算法，并将其与不完全LU分解预条件因子相结合来求解有限元方程，并与不同预条件因子的Krylov子空间迭代算法效率进行对比，结果表明ILU-IDR(s)迭代算法加快了求解方程的收敛速度，有效减少计算时间，为进一步实现电磁三维快速反演奠定基础.

致谢

我们向吴宇豪、王少博、薛园兵、唐钊、王恒、黄廷哲在论文准备过程中给予的帮助表示感谢.特别感谢两名匿名审稿人对本文提出宝贵的修改意见，才使得文章得以进一步完善.