0 引言

近十年来，海洋可控源电磁法因探测精度高、效率高等优点而被广泛用于海洋资源、深部构造勘探等方面，并取得了较好的成果(Constable and Srnka, 2007; Zhdanov, 2010).由于海洋环境特殊的沉积特性及地质运动等综合因素的影响，使得海底地下介质会呈现电阻率不均匀分布和复杂的各向异性特征，地电模型往往比较复杂.在这种情况下，如何实现高效、高精度的三维正演对于实现正确解释海洋可控源电磁数据有着重要的理论意义和实际价值(Newman et al., 2010; Ramananjaona et al., 2011).一直以来，海洋可控源电磁三维正演问题是国内外研究的热点问题.一般而言, 目前常用的数值模拟算法主要积分方程法(Wannamaker et al., 1984; Zhdanov et al., 2006)、有限差分法(Newman and Alumbaugh, 1995; Streich, 2009; Mittet, 2010; 殷长春, 2014)、有限单元法(汤井田, 2007; Da Silva et al., 2012; Ansari and Farquharson, 2014; 蔡红柱等, 2015; 陈汉波等, 2017)、有限体积法(Jahandari and Farquharson, 2014, 2015；杨波, 2014).有限差分法比较简单、直接, 但是采用的是规则网格对研究区域进行剖分，因此无法模拟比较复杂的地质模型(Mukherjee and Everett, 2010).积分方程法只需要对研究区域进行剖分，对于简单的模型，计算量比较小，精度较高.但是与有限差分法一样，通常采用的是规则网格剖分，不适合构建复杂的模型.并且其计算量会随着模型复杂度的增加或者规模的增大而显著增大(Cai et al., 2014).有限单元法分为矢量有限元法和节点有限元法，因为可以采用规则网格和非结构网格对研究区域进行剖分，灵活性较强，具有较强的模拟复杂的地形和地质体的能力(Jin, 2002).因此，有限单元法在地球物理数值模拟领域的应用越来越广泛(Chung et al., 2014; Cai et al., 2015).

与以上提到的数值模拟方法相比，谱元法(Spectral Element Method，SEM)是近些年来从其他研究领域中借鉴而来的将有限元法和伪谱法相结合的方法.该方法兼具有限元法处理边界和结构的灵活性以及伪谱法的精确性和快速收敛特性, 具有节省计算时间、存储空间等优点(Zhuang and Chen, 2017).该算法是Patera首先提出来用于求解Navier-Stokes方程(Patera, 1984).又经杜克大学的柳清伙教授推广用于电磁场数值模拟中去(Liu, 2002).经过近十年来的发展，谱元法已经广泛用于各种工程领域的求解中(Lee et al., 2006; Lee and Liu, 2007; Chen and Liu, 2009; Luo and Liu, 2009; Luo et al., 2010).谱元法的基本原理首先将整个正演区域进行离散成一系列小区域，然后在每个小区域内采用含有未知系数的已知的插值函数来近似表示子区域内的场的变化规律，该方法可视作为一种广义的有限元法，其与常规的有限元法既有相同点，又有不同点.相同的是，这两种方法均可采用结构网格和非结构网格对研究区域进行离散，可以模拟复杂的地形和地电结构，具有相同的离散过程.然而，谱元法在基函数的构建、节点的选择、空间映射和积分处理等方面却与有限元法有着显著的不同(Lee and Liu, 2005; Luo et al., 2009).例如：谱元法采用了能够改善电磁场插值精度的高效的插值函数(Liu et al., 2015).而有限元中通常采用的是二阶或者一阶的插值函数，因此，相比于传统的有限元法，谱元法具有较高的计算精度.并且，值得注意的是，谱元法可以采用稀疏的网格单元对研究区域进行网格剖分，即可通过较少的自由度获得与有限元法相同的精度，并可通过简单地提高插值函数阶数改善计算精度(Chung and Yu, 2015).

本文从麦克斯韦方程组方程出发，详细推导了海洋可控源电磁二次电场所满足的亥姆赫兹方程.再利用基于混合阶矢量基函数的谱元法在六面体网格中进行了求解，实现电导率任意各向异性介质中的二次电场的求解.设计典型的地电模型进行正演计算，并与一维解析解及有限单元法的数值解进行对比，对本文算法的有效性和可行性进行了验证.采用六面体网格单元进行网格剖分离散，有利于模拟复杂的海底地形和地电结构，便于研究海底地形对实际海洋可控源电磁勘探的影响.

1 正演理论 1.1 控制方程

在地球物理电磁勘探领域，通常采用低频的电磁信号，忽略位移电流，采用正弦谐变时间因子e^iωt，Maxwell方程可写为(Zhdanov, 2009):

(1)

(2)

其中ω为角频率，μ₀为自由空间的磁导率，J_s为源中电流的分布，σ为任意各向异性的电导率张量，表达式为

(3)

为了得到任意各向异性电导率张量σ，我们首先定义一个参考电导率张量σ_e, 该参考电导率张量的主对角线上的元素对应各向异性介质的三个主轴的电导率值，公式为

(4)

为了获得任意方向的电导率张量，我们需要对参考电导率张量进行三重欧拉旋转，如图 1所示(沿着x、y、z三轴旋转一定的角度).那么任意方向的各向异性介质的电导率张量可通过方程(5)进行计算，公式为

(5)

其中：

ϕ、Ψ、θ分别为沿着x、y、z轴进行欧拉旋转的角度大小.

为了避免源的奇异性的影响，我们采用二次场算法开展可控源电磁三维正演，假定总的电导率张量为背景电导率张量和异常电导率张量之和，公式为

(6)

电磁场可类似的写成一次电场和异常电场之和(Zhdanov, 2009)，公式为

(7)

将方程(7)代入方程(1)、(2)可得异常电场矢量所满足的赫姆霍兹方程为

(8)

其中，E_p为背景介质的电场，可通过计算全空间的、层状介质为背景时的电场.当计算出异常电场E_s后，可通过方程(9)求解出异常磁场，公式为

(9)

1.2 网格离散

谱元法与有限元法类似，都需要对研究区域进行网格剖分，通常可以采用矩形块、六面体、四面体单元对研究区域进行剖分，为了简化计算和实现对复杂的地质模型进行模拟，本文采用任意六面体对研究区域进行剖分.然而，为了计算六面体单元的刚度矩阵，我们需要将网格中的任意六面体单元进行等参变换，即将直角坐标系统下的不规则的网格单元(图 2a)转换成自然坐标系统下的长方体网格(图 2b).

1.3 参考单元的混合阶矢量基函数

谱元法是以精度高而著称，其计算误差会随着基函数的阶数增加而呈指数减少.在这篇文章中，我们引入Gauss-Lobatto-Legendre(GLL)混合矢量基函数.对于一维等参单元而言，其N阶GLL基函数为(Lee et al., 2006)：

(10)

其中L_N(ξ_j)为勒让德多项式，L_N'(ξ_j)为L_N(ξ_j)的导数, ξ_j则为GLL节点.

对于三维等参单元而言，其N阶混合阶矢量GLL基函数为

(11)

(12)

(13)

其中Φ_rst^ξ(ξ)、Φ_rst^η(η)、Φ_rst^η(ζ)分别是等参单元的ξ、η、ζ方向的矢量基函数.N_ξ、N_η、N_ζ分别是等参单元的ξ、η、ζ方向的矢量基函数的阶数.为了简洁了解GLL混合阶矢量基函数，在此给出三维空间的GLL混合阶的矢量基函数空间分布，如图 3所示.

那么，在等参单元中的空间二次电场E_s(ξ, η, ζ)可以采用拉格朗日-勒让德多项式进行表示，公式为(Lee and Liu, 2007)：

(14)

其中N=N_ξ(N_η+1)(N_ζ+1)+(N_ξ+1)(N_η)(N_ζ+1)+(N_ξ+1)(N_η+1)(N_ζ)是参考单元的棱边数.E_s^ξ(ξ_r, η_s, ζ_t)、E_s^η(ξ_r, η_s, ζ_t)和E_s^ζ(ξ_r, η_s, ζ_t)是单元上节点(ξ_rη_s, ζ_t)上的二次电场分量.

1.4 谱元法分析

对方程(8)，应用迦里金法后可得方程(8)的弱解形式为

(15)

其中Φ_i是第i条边的矢量基函数.

将方程(14)代入方程(15)后，并加入狄里克雷边界条件后，可得离散系统方程为

(16)

其中i=1, 2, …, N, K是单元总数.

方程(16)可以进一步化简为矩阵方程：

(17)

其中：

(18)

(19)

(20)

值得注意的是，方程(17)—(20)是在物理空间中的表达式，然后，我们需要计算参考单元下的相应表达式，在此，我们引入协变映射表达式为(Peterson et al., 1997)：

(21)

(22)

其中Φ和分别是物理空间和参考单元的矢量阶基函数，J表示雅克比矩阵，表达式为(Jin，2014)：

(23)

因此，体积微元dv在ξηζ坐标下可表示为

(24)

方程(16)中的参考单元的形函数和权函数的旋度可通过公式(25)—(27)进行计算，即：

(25)

(26)

(27)

那么方程(18)—(19)则可转化为

(28)

(29)

(30)

其中：

(31)

这里的V_i是第i个矢量基函数的旋度分量.

为了求解方程(17)，我们需要计算方程(28)—(31)，在此我们采用数值计算方法进行计算，具体方法可以参照文献(Cohen, 2013).由于我们所选的SEM单元基函数是正交的，所以最后生成的总体系数矩阵是对角矩阵，这是SEM重要的特性之一.

1.5 边界条件

为了使得方程(17)的解是唯一的，我们需要加入合适的边界条件.对于一次场/二次场分离算法而言，通常设定区域边界远离发射源，以致边界上的二次电磁场可忽略不计.因此，我们可在边界处加入狄利克雷边界条件：

(32)

其中∂Ω表示边界区域.

1.6 线性方程组的求解

对于方程组(17)的求解，常用有两种方法：一种就是采用直接求解算法，该方法比较简单直接, 精度高，但是内存需求高；另一种就是迭代求解算法，相比于直接求解算法，迭代求解算法内存消耗较少，易于实现.因此，本文采用不完全LU分解的IDR(s)迭代算法对方程组进行求解(Sonneveld and Van Gijzen, 2008).

2 数值模拟算例

为了检验本文算法的正确性和有效性，设计了三个典型地电模型开展正演计算.具体计算平台参数为：2个16核CPU(Intel(R) Xeon(R)E5-2630), 主频2.40 GHz, 内存64 G, Win7 64位操作系统.

2.1 算法验证

为了验证本文算法的正确性，设置一个各向异性的水平层状介质，如图 4所示.第一层为海水层, 厚度是300 m，电导率为3.2 S·m^-1.第二层为沉积层, 电导率为1 S·m^-1.有一各向异性介质薄层赋存于沉积层中，厚度为100 m, 埋深为1300 m.参考电导率为σ_e=diag(0.01, 0.01, 0.025)，并令各向异性电导率张量沿着z轴旋转30°，即φ=0，Ψ=0，θ=30°.通过给定的参考电导率张量σ_e和三次欧拉旋转角度φ=0、Ψ=0、θ=30°，利用方程(5)即可计算沿着z轴旋转30°时的电导率张量.采用水平电偶极子源进行激发，源点处于海底水平面以上30 m, 具体坐标是(0, 0, 270).频率为0.25 Hz.该模型是参考Loseth和Ursin(2007)中的模型，利用本文提出的算法计算的结果将与Loseth和Ursin(2007)的结果进行对比.电场幅值对比结果如图 5所示, 从图 5可知, SEM数值解与解析解比较吻合, 图 6为数值解与解析解的相对误差, 从图 6可知, 相对误差均低于1.2%.数值解与解析解对比结果表明本文算法是正确的.

2.2 算例2

为了检验本文算法的效率，设计一个金字塔起伏地形条件下的主轴各向异性介质模型，即不考虑电导率张量发生旋转的情况.金字塔地形含异常体模型如图 7所示.金字塔的底面范围为Ω={-1000 m≤x, y≤1000 m}, z=1000 m, 顶面范围为Ω={-500 m≤x, y≤500 m}, z=800 m.空气和海水深度均为1000 m.电导率分别为10^-6 S·m^-1和3.3 S·m^-1.沉积岩为各向异性介质，电导率张量为σ=diag(1.0, 1.0, 0.8).基岩的电导率为0.1 S·m^-1.异常体为各向同性介质，电导率为0.05 S·m^-1, 其范围为Ω_a={-1000 m≤x, y≤1000 m, 2400≤z≤2600 m}.采用水平电偶极子进行激发，源点坐标为(-3000, 0, 950).激发频率为0.5 Hz.接收站均放置于海底平面以上50 m处.研究区域设计为Ω={-4 km≤x, y≤4 km, -1 km≤z≤3 km}.谱元法的矢量基函数阶数为N_ξ=N_η=N_ζ=4.将谱元法的数值计算结果与有限单元法数值解进行对比，检验本文算法的效率.图 8和图 9是谱元法和有限单元法计算的异常电场E_x和异常磁场H_y对比图.从这两个图上可以观察到，SEM和FEM计算结果基本一致.从图 8和图 9，我们还发现，海底地形的可对海洋可控源电磁响应产生比较明显的影响, 进一步验证本文算法的有效性.图 10为各向同性和各向异性沉积岩下的二次电场E_x分量幅值对比，从图 10可知，电导率各向异性同样对电场E_x分量产生显著的影响.表 1为这两种算法所消耗的计算时间和内存之比.从表 1可以看到，在相同自由度(Degree of freedom，DOF)的条件下，SEM的计算消耗比FEM要小.

2.3 算例3

为了进一步验证本文提出的算法，设计一个海洋沉积模型如图 11和图 12所示.空气层厚度为1000 m，电导率为10^-6 S·m^-1.海洋厚度为1000 m，电导率为3.3 S·m^-1.各向异性异常体大小规模为2000×2000×500 m，各向异性电导率为σ=diag(0.01, 0.01, 0.025)S·m^-1.沉积层为各向异性介质，参考电导率为σ_e=diag(2, 1, 1).模型区域大小定义为Ω={-5000 m≤x, y≤5000 m, -1000 m≤z≤4000 m}.采用x方向的水平电偶极子进行激发，频率为0.25 Hz，源点坐标为(-3000，0，950).接收站均处于海底平面(z=1000 m).我们分别考虑两种沉积层电导率旋转情况，即分别沿着y、z轴旋转30°、45°、60°(φ=0，Ψ=30°, 45°, 60°，θ=0和φ=0，Ψ=0，θ=30°, 45°, 60°)，研究海洋可控源电磁电场异常响应特征.谱元法的矢量基函数阶数为N_ξ=N_η=N_ζ=4.

图 13为沉积层电导率张量不发生旋转时的海地平面电场分布.图 14—15分别为沉积层电导率向y轴和z轴旋转30°、45°、60°时的海底电场分布.对比图 13和图 14—15.我们可以发现, 随着电导率张量发生旋转，电场三分量的分布均发生了明显的变化.从图 14可以看到，随着向y轴旋转的角度越大，电场E_x分量幅值变大，从中心位置向外扩展分布.而E_y分量的变化趋势与电场E_x分量相反，即幅值变小，向中心位置收缩.电场E_z分量变得明显地不对称，由一开始的左边的响应值大于右边，变化成右边的响应值大于左边.而且从E_z的分布中存在一个低幅值区域，可以看出与异常体在海底平面的投影位置一致.从图 15可以看出, 随着沿着z轴旋转的角度增大，电场E_x分量由近似对称分布变为不对称.E_y分量也亦如此.E_z分量，随着沿着z轴旋转的角度增大，变为沿着z轴旋转的分布特征.因此，可知电场的分布能够反映出海底各向异性介质的特征，有助于识别海底各向异性介质分布规律.

3 结论

本文采用谱元法结合混合阶的矢量基函数，实现了电导率呈任意各向异性介质的海洋可控源电磁三维数值模拟研究.采用一次场/二次场分离算法，相对于总场算法而言，可以有效回避引入场源所带来的奇异性.设计水平层状介质模型进行正演计算，计算结果与解析解进行对比, 验证了本文算法的正确性和有效性.通过设计典型的三维各向异性介质模型进行正演计算，并与有限单元结果进行对比, 对比结果展示了谱元法的优势.深入分析任意各向异性介质下的电场分布特征，有利于识别地下介质各向异性特征分布规律.本文算法具有良好的通用性，适用于井中电磁、海洋电磁、环境地球物理等非均匀介质中的电磁感应基础研究.笔者在后续研究中，将会以本文提出的正演方案为基础，进一步实现各向异性介质的海洋可控源电磁三维反演.

致谢  感谢石会彦、李少朋、王恒、李根在论文准备过程中给予的帮助，特别感谢两名匿名审稿人和编辑对本文提出宝贵的修改意见，才使得文章得以进一步完善.