基于分区并行的2.5维海洋可控源电磁法反演研究

时间：2024-05-22

张迎利，李赓

(河南理工大学物理与电子信息学院，焦作 454000)

0 引言

海洋可控源电磁法(Marine Controlled Source Electromagnetic method, MCSEM)广泛地应用于海底油气资源和天然气水合物储层探测研究[1-2]。随着能源需求的不断增长和陆上油气资源储备量的日益减少，海洋油气资源探测逐渐成为能源争夺的焦点。近些年来，MCSEM被广泛应用于海底油气资源探测，并已经取得了明显的效果[3-4]。海洋电磁正演策略分为二次场和总场算法[5-6]，为克服场源奇异性问题，二次场算法需要利用一维解析公式计算背景场，而总场算法需要计算源项[7]。一般来说，对于简单模型而言，二次场算法的计算精度高，但对于复杂模型，由于二次场算法很难确定背景模型，因此总场算法的适应性更强[8-9]。

正演通过使用非结构化网格，可以很好地适应起伏地形和地下复杂电性结构[10]。但是在大规模电磁数值模拟中，通过增加适当的网格来适应地下电性结构是提高计算效率的关键。近些年，自适应有限元方法针对以上问题给出了很好的解决方案。其中二维电磁场数值模拟中基于非结构化网格的自适应有限元法已经取得了很好的应用效果[11-13]。

海洋可控源电磁反演研究中，Abubakar等[14]基于高斯-牛顿方法提出了一种快速、有效的2.5维反演算法;Chen等[15]采用非线性共轭梯度算法实现了MCSEM 2.5维反演;李刚等[16]采用有限差分法实现了2.5维频率域MCSEM正演，并通过直接法求解方程组;赵东东等[17]基于有限元方法对2.5维MCSEM进行了正演模拟，引入网格加密收缩算法提升了数值的计算效率;陈光源等[18]利用非线性共轭梯度法实现了2.5维海洋可控源电磁反演,并探讨了不同反演参数对计算速度的影响;Key等[12]实现了频率、发射源并行的2.5维自适应有限元OCCAM反演算法，提高了反演计算速度，但仍需要计算、存储整个反演域上的灵敏度矩阵。基于光滑模型约束的OCCAM反演方法，反演过程稳定，但结果不能清晰刻画地下电性结构的边界[19]。为了解决光滑模型约束的问题，Portniaguine等[20]提出了聚焦反演成像方法，其方法的优点在于能够得到较为清晰地地电结构图像，但反演过程不稳定，依赖初始电性模型的选择。在MCSEM实际资料解释过程中，反演计算量大、内存需求量高[21]。

笔者运用二维海洋可控源电磁资料处理，正演采用基于总场的自适应有限元方法，其中网格剖分采用非结构化三角单元。为提高计算速度，降低内存使用量，我们将反演区域进行分解计算。反演过程中利用地震资料和光滑模型进行约束反演，得到背景电阻率，然后将其作为初始模型进行模型空间约束反演，从而获得工区内高阻区域分布。

1 正演模拟

1.1 控制方程

在二维海洋地电模型中，设走向方向为x、y方向轴水平向右，z方向轴垂直向下，假设时谐因子为e-iωt，则电场和磁场满足的控制方程为[11]

(1)

式中:E为电场强度；ω为角频率；μ0为真空中的磁导率；H为磁场强度；σ为介质电导率；J为外加源电流密度。

将式(1)中电磁场分量沿着构造走向(x方向)做傅里叶变换[22]

(2)

(3)

(4)

(5)

(6)

(7)

(8)

1.2 自适应有限元

对于简单模型，依据经验对网格进行剖分，可以得到精度较高的数值解；而对于复杂模型来说，需要把有限元网格剖分得足够细，但网格过细会影响到计算效率。在自适应有限单元法中，通过后验误差估计方法获得每个网格单元上的误差分布，选择性的加密误差较大的部分单元，可以快速提高解的精度，同时避免全局网格加密所带来的计算时间和内存使用量的急剧增加。

(9)

式(9)的离散形式表示为

B(un,v)=F(v)

(10)

其中，un为有限元数值解，u为精确解，令全局误估计算子εn≈u-un。由于精确解u未知，因此εn不能直接求解，但是可以通过近似误差函数计算。

B(εn,v)=B(u-un)=F(v)-B(un,v)

(11)

定义误差函数[12]

(12)

图1 自适应有限元正演算法流程图Fig.1 Adaptive finite element forward algorithm chart

图2 1D模型示意图Fig.2 1D model

1.3 正确性验证

为验证正演算法正确性，设计如图2所示海洋模型，使用本文正演算法计算其正演响应，将计算结果与解析解[24]进行对比。其中，海水层厚度为1 km，电阻率为0.3 Ω·m，在海底下方1 km处有一厚度为0.1 km的高阻储层，电阻率为100 Ω·m；发射源置于海底上方0.05 km处，坐标为(0 km,0 km,0.95 km)，发射电流为1 A，发射频率为0.3 Hz；30个海底接收点等间距分布于y=0.5 km～15 km的范围内。将本文正演计算的电磁响应(第六次自适应加密)与解析解进行对比，结果如图3(a)所示，可以看出数值解与解析解吻合较好。图3(b)为相对误差曲线，可以看到，在自适应加密初期相对误差很大，随着自适应次数增加，相对误差逐渐下降到1%以下，表明本文的正演方法具有较高的计算精度。

2 基于背景模型约束的OCCAM反演

反演过程中，如果知道某区域地球物理场的背景信息，那么可以将其作为约束增加到反演中[25]，但是背景信息往往不易直接获得。为了得到背景信息并且使其值更加逼近真实，利用已知的地震层位信息对反演区域进行剖分，将每个层位的电阻率作为一个参数参与反演，得到工区的电性参数背景(对应公式(13)，记为第一步)

图3 数值解与解析解对比Fig.3 Comparison of numerica and analytical solutions(a)电场Ey振幅响应曲线;(b)不同迭代次数下相对误差曲线

(13)

然后选择它与模型m之差的二范数定义为模型空间的目标函数，结合数据空间的目标函数可形成新的目标函数，对新的目标函数继续实施反演直至收敛(公式(14)，记为第二步)。

U=‖P(m-m*)‖2+

(14)

对于给定的初始模型mk，为了使总目标函数最小，需要求取目标函数的梯度，令mU=0，可得到新的迭代模型mk+1[18]：

对式(13)有

(15)

对式(14)有

mk+1=[λPTP+(WJk)TWJk]-1×

(16)

(17)

反演迭代过程中，通过计算模型的修改量以确定新的模型mk+1，从而计算、判断模型的正演响应与观测数据之间的拟合差是否达到设置要求，实现目标函数最小化求解。过程如下：首先利用层位信息建立反演初始模型，移除层位中模型粗糙度，并设置反演电阻率的左右边界范围为10-1Ω·m～103Ω·m，对其整体实施反演，直到RMS收敛为止。其次，用第一步反演结果作为初始模型，设置对应的P值，在目标区域范围内采用矩形剖分，扩展区三角剖分，最大纵横比设置为10。实施第二步反演，直到RMS不再下降为止。

3 分区并行方式

OCCAM反演方法需要计算和存储整个反演域上的灵敏度矩阵，从而耗费大量计算时间和内存。根据航空footprint概念，某测点响应大小与其下方一定范围内地电模型有关[26]。因此，为了提高反演计算效率，我们尝试将footprint引入到海洋可控源电磁中。如图4所示，计算了频率为0.3 Hz的10个接收机为一组和发射源(O)分别位于-20 km、-5 km、10 km处的归一化灵敏度示意图。通过对MCSEM影响范围进行分析可知，归一化灵敏度值主要集中测点下方所在的有效区域范围内，虚线范围外归一化灵敏度值相对很小。由于接收站附近“敏感区域”远小于整个测区，因此在海洋电磁数据反演时，可先基于“敏感区域”大小对测区进行单元划分，从而减小反演模型规模，在完成各测区单元的电磁数据反演后，组合得到反演区域内总的电性分布。

图4 归一化灵敏度示意图Fig.4 Normalized sensitivity diagram(a)-20 km;(b)-5 km;(c)10 km

图5 模型一Fig.5 Model I(a)模型示意图;(b)传统方法反演结果;(c)分区并行方式反演结果;(d)传统方法、分区并行方式反演RMS对比

4 反演算例

4.1 模型一

为了测试分区并行方式的有效性，设计了如图5(a)所示的二维地电模型。参数设置如下：空气电阻率为109Ω·m；海水层厚度为1 km，电阻率0.3 Ω·m；海底岩层电阻率1 Ω·m，在海底下方0.6 km处有一个8 km×0.2 km的高阻薄层，电阻率大小为50 Ω·m；电偶源布置在海底上方0.05 km处，沿测线方向等间距分布100个(-25 km ≤y≤25 km)，发射频率为0.3 Hz，发射电流1A；100个海底接收点等间距分布于y=-24 km～25 km的范围内。正演计算得到电场Ey分量的幅值、相位数据，加入4%高斯噪声，作为反演的拟合数据。反演初始模型由空气-海水-岩层组成，其中空气、海水层不参与反演，岩层初始电阻率为1 Ω·m。图5(b)、5(c)展示了单节点28个处理器反演迭代20次结果，由图可以发现二者反演结果很相似，异常体的位置与实际情况基本吻合，并且电阻率值和真实的电阻率值很接近。图5(d)展示了单节点28个处理器反演迭代20次RMS(数据拟合误差)曲线，由图可以得出，二者初始RMS均为17，经过20次迭代最终RMS均收敛到1.0左右。为了体现分区并行方式的优越性，计算了单节点28个处理器反演迭代20次的耗时，传统方法需要约29小时，分区并行方式反演需13.5小时左右，前者比后者的反演耗时增加了约两倍。此外，分区并行比传统方法反演的内存减少1/3左右。

图6 反演并行结果Fig.6 Inversion of parallel results(a)不同处理器核数反演耗时;(b)分区并行方式加速比曲线

加速比是衡量并行处理性能的重要指标之一，因此，首先定义加速比为

SP=T1/TP

(18)

其中:SP为加速比;p为并行计算所用处理器的核数；T1为串行计算所花时间；Tp为p个处理器并行计算所花时间。

为了表示效率的提升情况，在图5(a)所示模型基础上，计算了不同处理器核数(1、2、6、12、18、24、28)反演迭代20次消耗CPU时间以及相应的加速比。并行计算环境为一个拥有28个Intel(R) Xeon (R) CPU E5-2620 v4 @ 2.10GHz的处理器构成的集群上，操作系统为SUSE Linux Enterprise Server 11 (x86_64)，内存128GB。图6(a)为传统方法、分区并行方式反演迭代20次消耗时间对比。由图6可以看出，处理器核数相同时，传统方法比分区并行方式反演耗时增加了两倍左右。图6b展示了加速比与并行计算所用处理器核数之间的关系，随着计算规模的逐渐增加，计算效率有很大提升，同时取得了很好的加速比。

4.2 模型二

为了测试背景模型约束的OCCAM反演方法的有效性，设计了如图7(a)所示的模型。其中，海底地形和地层信息来自实际数据资料，每层的电阻率均以各向同性的方式填充。空气的电阻率为1×109Ω·m；海水电阻率为0.3 Ω·m；海底以下有一高阻薄层，宽度约6 km，深度为3.2 km～3.4 km，电阻率大小为60 Ω·m；电偶源布置在海底上方0.05 km处，发射频率为0.1 Hz，发射电流为1A；15个接收点布设在海底，范围为-13 km～6 km，间隔为1.35 km；正演计算电场Ey分量的幅值、相位数据，加入4%高斯噪声，作为反演的拟合数据；反演初始模型由空气-海水-岩层组成，空气、海水层不参与反演，岩层初始电阻率为1 Ω·m半空间。图7(b)展示了光滑反演(式(13))结果，可以看出，光滑反演结果中异常体边界不够清晰，并且异常体的纵向范围有所拉伸。图7(c)、图7(d)为背景模型约束的OCCAM反演结果，由图7(c)可以得出，整体电性结构与真实电性结构很接近，但靠近海底的两层信息没有明显反映出来，是因为它们的(海底～2.5 km、2.5 km～2.9 km)真实电阻率值很接近。由图7d可以看出，实施第二步反演后(式(14))，靠近海底的两层信息略微呈现了出来，将其与光滑反演结果(图7(b))对比可以发现，整体电性结构以及电阻率变化趋势与真实情况基本一致，异常体电性边界结构不仅更加清晰，而且位置更加准确、形态更加饱满。图7(e)、图7(f)为RMS随迭代次数的变化曲线，由图7(e)可得，初始RMS约为13，经过14次迭代RMS降到1.9左右；由图7(f)可以看出，经过第二步反演后RMS最终收敛到1.2左右。

图7 模型二Fig.7 Model II(a)模型示意图;(b)光滑反演结果;(c)第一步反演结果;(d)第二步反演结果;(e)、(f)RMS曲线