0 引言

海洋可控源电磁法(Marine Controlled Source Electromagnetic, 简称MCSEM)在海洋领域广泛地应用于洋中脊、大陆边缘地质构造以及海底油气资源和天然气水合物储层探测(Constable, 2010; Key, 2012)研究.正演采用结构化网格很难模拟地形和复杂地质体.通过使用非结构化网格，可以很好地适应起伏地形和地下复杂电性结构(杨军等，2015).但是，在大规模电磁数值模拟中，如何通过增加适当的网格来适应地下电性结构是提高计算效率的关键.面向目标的自适应算法与人工对目标体，激发源和接收站位置进行局部加密的算法相比，前者使用后验误差判断网格是否需要被优化，可以自动地对网格进行剖分，对计算区域的划分更合理.近些年，自适应有限元方法针对以上问题给出了很好的解决方案.其中，二维电磁场数值模拟中基于非结构化网格的自适应有限元法已经取得了很好的应用效果(Li and Key, 2007; Key and Ovall, 2011a; Li and Dai, 2011).针对三维大地电磁正演问题，国内外学者引入后验误差估计方法，得到了精确的电磁模型(Ren et al., 2013; Grayver and Kolev, 2014; 殷长春等, 2017).针对三维海洋电磁问题，基于背景/异常场方法实现了可控源电磁场三维矢量有限元正演模拟(李勇等, 2017; 刘颖等，2017).基于直接法求解系数矩阵和数据空间分解的方法实现了三维矢量有限元正演模拟(Zhang et al., 2015).

在海洋电磁模拟中，电磁场数值模拟策略分为背景/异常场和总场算法.为克服场源的奇异性问题，背景/异常场算法需要利用一维解析公式计算背景场.总场算法需要计算源项(Haber, 2014)，并通过靠近源的网格细化方法减少源的奇异性影响(Plessix et al., 2007).一般来说，对于简单模型，背景/异常场算法的计算精度要高，对于复杂模型，由于很难确定背景模型，总场方法的适应性更强(Jahandari and Farquharson, 2015).

有限元的条件数主要依赖于电导率对比度(Schwarzbach et al., 2011; Da Silva et al., 2012; Um et al., 2013).其中，空气层电导率比海水电导率低得多，电导率对比度为10⁸甚至更大.因此，三维海洋可控源电磁问题中系数矩阵的条件数通常很大，数值求解此类问题主要是采用直接法或迭代法.如采用直接法，由于计算的需求与问题大小呈非线性增长(Operto et al, 2007; Pardo et al，2012)，直接法更适合应用于中小规模矩阵的MCSEM问题(Streich, 2009; da Silva et al., 2012; Puzyrev et al., 2013).对于较大规模问题，想要构造一个高效的预条件使得方程快速收敛是很困难的.目前，快速辅助空间预处理技术可以满足三维大规模电磁正演的求解需求(Hiptmair and Xu, 2007; Grayver and Bürg, 2014).

本文针对复杂模型的三维海洋可控源电磁问题，采用基于总场的方法使得适应性更强；采用非结构化的六面体网格和八叉树的加密方式可以更容易在下一阶段自适应有限元反演中实现正反演网格的分离技术；采用迭代法可以满足三维大规模电磁正演问题中的求解需求.通过对比三维程序和半解析代码(Key, 2009)在层状模型的可控源电磁响应，验证了数值算法的精度.对于一个1300多万自由度的三维电性结构模型，分析了总体用时、内存使用以及不同自适应加密次数下相对误差变化情况.其中，线性方程组求解采用辅助空间迭代法，内存消耗仅线性增加，验证了算法的可扩展性.MCSEM大规模有限元计算采用Serial Graph Partitioning and Fill\|reducing Matrix Ordering(METIS)进行均衡负载，展示了不同迭代次数下网格划分情况.

1 正演理论 1.1 可控源电磁的控制方程

三维各向同性电导率模型中，假设电磁场时谐因子为e^-iωt，MCSEM控制方程为

(1)

(2)

从法拉第定律(1)和安培定律(2)，得到电场的赫姆霍兹方程(Ward and Hohmann, 1988)：

(3)

E为电场，H为磁场，I表示电流幅度，x₀表示源位置，ds表示在x方向电流源长度，ω为角频率，σ为电导率，δ(x, y, z)为狄拉克函数.更加复杂的源可以通过电偶源的组合得到(Haber, 2014).为了保证解的唯一性，在表面上强加dirichlet边界条件：

(4)

其中是上指向外的单位法向量.

1.2 矢量有限元分析

从得到的方程(3)出发，为避免节点有限元带来的伪解现象(Jin, 2002; Monk, 2003)，我们使用矢量形函数，并将形函数定义在网格棱边上，e_i为形函数(如图 1).

在微分方程的两边分别乘以权函数，本文取权函数与矢量形函数一致并积分，得到控制方程的弱解形式为

(5)

通过分部积分，可以获得微分方程的泛函形式：

(6)

采用第一格林函数理论，i为单元数，方程(6)变化为以下形式(Jin, 2014; Ansari and Farquharson, 2014):

(7)

1.3 方程求解

计算区域剖分为有限多个六面体单元网格，在求解区域内每个小单元上对上述泛函形式积分，并施加边界条件，可以获得一个大型复系数稀疏线性方程组：

(8)

其中K为复系数矩阵，u为待求解，s为右端项, 任意一点的场值E可表示为

(9)

线性方程组的求解采用迭代算法，迭代算法的效率主要取决于所使用的预条件.对于由Maxwell方程组离散得到的复数方程组，很难构造一个高效的预条件使迭代算法可以快速收敛.我们将式(8)分解为两个实数(u_r, u_i)形式进行求解.

(10)

其中K_r为矩阵K的实部，K_i为矩阵K的虚部.从数值计算角度，为了使系数矩阵具有对称性，我们将方程(10)的第二行乘上-1，方程(10)进一步可表示为

(11)

为简便起见，我们将式(11)中的系数矩阵记为A.如果直接采取Krylov子空间方法求解方程(11)，收敛速度会变得很慢，为此，采用预处理技术来克服这个问题，将线性方程组Ax= b转换为另一具有更优良谱性质的方程组，然后采用FGMRES迭代求解.其中预条件P为

(12)

B = K_r+ K_i是对称正定且可逆的.条件数(P^-1A)具有上界并且与矩阵的规模、频率和相对电导率无关(Chen et al., 2010, 2007).

实对称线性方程组的求解采用FGMRES迭代算法(Saad, 2003)，在使用P作为预条件的情况下，需要计算P^-1与任意向量的乘积.因为P具有对角特性，P^-1与任意向量的乘积可转化为求解两次如下的线性方程组:

(13)

其中，B的求解使用AMS预条件的共轭梯度法计算(Kolev and Vassilevski, 2009; Baker et al., 2012).

1.4 面向目标误差估计

在获得微分方程的解之后，我们需要对其进行后验误差估计来指导网格加密.全局自适应算法除了加密对观测点模拟误差较大的单元外，对其他单元也进行了加密，从而浪费了大量的计算资源.面向目标自适应方法可以根据解的性态来进行网格点的自动分布，克服一致加密引起计算机物理内存和CPU计算时间的剧烈增加.本文自适应方法中误差指示子η_K(E)的物理意义在Grayver & Bürg(2014)得到了充分讨论：

(14)

在单元K上的残差估计值可表示为

(15)

本文利用单元i界面处磁场的切向连续性和电流密度的法向连续性条件估计后验误差：

(16)

对于需要加密的网格，将其分为8个，使用基于八叉树的数据结构来存储，如图 2.在海洋可控源电磁正演模拟中，观测点往往仅包含于少量单元内，这些观测点处的数值模拟精度是人们关心的焦点.在测点处假设源的存在，从而获得一个对偶问题的解E^D，误差估计子可进一步表示为

(17)

在本文中，每次加密所有误差单元前10%的网格，直到达到迭代次数或自由度数后停止.

(18)

2 数值模拟算例

为了检验程序的正确性，我们设计了层状模型算例，一维层状模型是有解析解的，是检验数值解正确性的最佳模型.自适应加密策略的鲁棒性和可扩展性的验证，采用三维地电结构模型进行模拟.我们将程序部署在工作站上，系统配Intel Xeon E5\|2680V4 2.4 GHz十四核处理器，128GB内存.在以下算例中，FGMRES迭代截止条件是相对残差下降到10^-8.

2.1 正确性验证

通过比较一维和三维程序对于层状模型(如图 3)的响应特征以及相对误差随自适应迭代次数的变化规律，从而验证本文算法的精度.层状模型包含空气层，1 km厚的海水层，0.1 km厚的储层，储层从3 km开始，电导率为0.1 S·m^-1，空气电导率为10^-6 S·m^-1.间距0.5 km，接收站从0 km到15 km处，我们采用Kerry key公开的一维正演程序与面向目标自适应有限元程序进行模拟并对比响应曲线.图 4给出了电磁场六个分量第11次自适应加密的数值解与解析解的结果.可以看出，本文模拟结果与一维解析解模拟结果吻合较好.图 5展示了两种算法的相对误差曲线，从图 5给出的相对误差曲线可以发现，面向目标自适应算法在自适应加密初期相对误差达到95%，随着自适应次数的增加，电磁场水平方向相对误差可下降到1%以下，而垂直分量相位误差略大.证明了面向目标自适应算法具有较高的计算精度.发射源放置在距离海底0.1 km处，发射频率为0.25 Hz，具有45°旋转角，从而可以对比电磁场六个分量数据.

2.2 三维起伏地形模型及误差指示子分布

前面简单模型证明了程序的准确性和自适应矢量有限元的有效性.为了考察起伏地形对于海洋可控源电磁场的影响情况，在本节展示了一个更加复杂的模型，通过等参变换来模拟地形，给出模型如图 6所示，山脊顶部面积为1 km×1 km，底部矩形面积为3 km×3 km，山脊高度为0.3 km.发射源放置在距离海底0.1 km处，发射频率为0.25 Hz.采集站主剖面放置在X＝0测线处，沿-15 km到15 km处扩展，间距0.5 km.三维模型的空气，海水电导率和2.1节中模型一致.

本文算法对接收点附近的起伏地表进行了很好的网格加密，面向目标的自适应算法在观测点附近将网格划分得很密，而在远离观测点处，网格剖分比较稀疏.

表 1展示的是每次迭代自由度信息和求解方程耗费的时间，总体用时以及内存使用.其中，自由度数量为复数线性方程组的自由度数量.FGMRES迭代次数不依赖于自由度数量的增长，CG迭代次数随着问题规模的增加有轻微的增长，体现了算法的鲁棒性.

图 7展示的是平地形均匀半空间模型与带地形均匀半空间模型的对比结果，其中左边为振幅，右边为相位.从给出的模拟结果可以看出，本文算法的响应曲线光滑，场值在起伏地形区域被扭曲.

自适应策略中一开始进行网格加密的误差是最大的，进行误差估计以后，误差具有显著的减小.其中，暖色代表误差指示子大的区域.图 8a展示了初始网格的误差分布，在电导率对比度大的区域，误差指示值最大.网格细化之后，图 8b展示了十二次自适应加密后的误差指示，整个区域的误差显著下降.通过面向目标的自适应方法把更多的自由度放在那些解函数剧烈变化的地方, 使得较小的计算量获得较大的计算精度.

2.3 三维模型计算区域均衡负载

为了测试解决方案的可扩展性，我们构造了一个波浪地形模型，发射源放置在距离海底0.1 km处，发射频率为0.25 Hz，采集站放置在Y=X处，沿-15 km到15 km处扩展，间距0.5 km，如图 9a和9b.图 9c和图 9d分别展示的是第1次和第10次网格细化，在网格细化过程中，大多数细化集中于接收点、发射源和异常体的区域，避免了采用结构化网格易导致计算界面处出现扁平单元和狭长单元的现象.本文算法不仅在观测点附近的起伏地表处对网格进行了很好的加密，而且对具有一定埋深的异常体表面也进行了网格加密.此外，在异常体的顶部，网格加密程度要高于底部，这由于随着埋深的增加，介质对观测点影响的程度逐渐变小，因而减小了误差指示子.为提高三维大规模电磁并行数值求解的速度，关键因素之一是处理好进程间的数据分布，通过将数据分布存储在各进程中，以保证负载平衡.传统有限元程序在计算过程中，仅需在开始做一次静态划分，数据分布保持不变.然而对自适应有限元程序，每个进程的负载是在不断地变化的，需要动态地调整数据分布以保持负载平衡.本文采用开源的METIS(Karypis and Kumar 1998)区域分解函数库在计算过程中动态地对目标区域计算任务进行分配，使得三维复杂模型的每个子区域计算任务量大致相等.图 10展示了一个动态负载平衡进行并行计算的例子，初始计算区域通过METIS被分为14个子区域, 如图 10a.通过8次自适应迭代后，计算量较大的区域更加集中在发射源附近(如图 10b)，总的单元数为3974826，每一个区域单元数大体为总单元数的1/14, 即283916.因此，计算区域通过频繁的划分以满足动态负载平衡.

本文分别用1、2、4、8和14个CPU核心对该模型进行正演模拟，以比较本文使用并行策略的效率, 计算时间如图 11所示.可以看出，本文获得了不错的加速比.

3 结论

本文实现了一种三维复杂结构的面向目标自适应海洋可控源电磁有限元方法.通过自适应网格细化来提高有限元解决方案的准确性，对于大规模三维MCSEM数据，我们使用了METIS来进行网格计算任务量的划分，使并行计算达到均衡负载的效果.通过对比一维解析解与三维自适应MCSEM计算结果，验证了程序的正确性.三维复杂模型在不同迭代条件下的网格划分和误差估计子的计算，证明了程序的鲁棒性和可扩展性.

致谢  感谢河南理工大学高性能计算中心支持.特别感谢审稿专家和编辑的建议.