1. 引言
海洋地震是海洋油气勘探的传统方法。然而,海洋地震探测海底油气存在固有的缺陷——它在探测出可能的储油构造后无法区分其含油/含水特性,并且勘探成本高。海洋可控源电磁法(Marine Controlled Source Electromagnetic Method, MCSEM)则可以很好地弥补它的不足。MCSEM是一种主动源电磁勘探技术,该方法通过在海洋环境中激发低频率的电磁场,利用接收器采集经过地下介质传播后的响应信号,进而通过数据处理与反演有效获取海底地下介质的电性分布信息[1]-[4]。MCSEM依托电阻率能良好反映油气、水合物以及矿产等地质体与其围岩的差异,可以有效地提高海洋油气勘探的成功率[5],减小勘探成本,因而在深海油气资源、天然气水合物、海底矿产以及构造演化研究中得到了广泛应用[6]。
当前,MCSEM最常规的观测方式是“拖曳发射源 + 海底布设接收器”,即通过调查船拖曳电磁发射源(多为水平电偶极子源)在目标海域航行,同时将大量接收器预先布设在海底固定点位,形成接收阵列[7]。作业时,拖曳在海水中的发射源向海底发射特定频率的电磁信号,信号穿透海水与地下介质后,由海底接收器记录电磁场响应,待完成一定区域的测量后,再通过调查船回收接收器并提取数据。由于接收器是直接置于海底,远离海面波浪以及船体噪声,信噪比相对较高,但局限性也较为明显,一方面,接收器的海底布设与回收需耗费大量时间,尤其是在深海区域,单次布放可能需要数小时甚至数天,且受海况、海底地形影响极大,不仅拖慢了作业进度,还可能因设备丢失、损坏导致数据缺失;另一方面,在接收器投放时,由于洋流等因素导致接收器未能投放到预定位置;此外这种观测模式下的勘探效率较低,难以实现大范围、高密度的连续测量,这种时间跨度的作业模式进一步推高了成本,限制了海洋电磁法在大规模资源勘探中的应用普及。正是这些固有局限,催生了对更高效、集成化作业模式的探索,为深拖式发射接收装置的研发与应用奠定了现实基础。如Scripps研究所开发的Deep-Towed EM系统、挪威EMGS的SeaBed Logging系统、WesternGeco的Stingray系统[8]-[10]等都是该技术的探索研究。
目前常用的电磁三维正演方法常用有限单元法(FE)、有限差分法(FD)、积分方程法(IE)以及有限体积法(FV)。FE可以用灵活的单元拟合复杂几何异常体和地形,精度可控但计算成本高[11];FD数学思想简单,易于编程实现、效率高,但依赖于结构化网格,难以处理复杂几何形状和边界[12] [13];IE的核心是借格林函数建立目标异常区电流与电磁场的积分关系、仅离散异常区域而非全空间以简化求解,虽具计算效率较高、天然满足无穷远边界条件无需额外处理吸收边界的优势,但存在对复杂地形(如剧烈起伏海底)和强电性差异结构模拟精度有限、积分核计算易产生大规模稠密矩阵导致内存占用高而受限大规模三维问题应用的缺点[14];FV是将研究区域划分为一系列互不重叠的控制体积,并在每个控制体积上对控制方程的积分形式进行近似求解。通过这种方式,可以在保证局部和整体守恒性的同时,将连续的偏微分方程转化为离散代数方程。与有限元或有限差分方法相比,有限体积法天然满足物理量的局部守恒性,因此在电磁场数值模拟中具有显著优势[15] [16]。由于MCSEM的勘探规模,既要兼顾方法适用性,又要考虑在有限的计算成本下兼顾计算效率,本文采用拟态有限体积法(Mimetic Finite Volume, MFV)的离散策略[17] [18]。
MCSEM一、二维正反演理论已经很成熟了,并且在实际生产中得到了广泛的应用。一、二维方法虽能简化计算、快速获取地层垂向或剖面电性特征,但难以准确刻画实际地质体(如倾斜油气藏、不规则矿体、非均质地层界面)的三维空间形态与电性参数分布,而这正是实际生产所追求的。MCSEM的三维正演理论研究已经走向成熟。当然,MCSEM的三维反演研究也涌现了很多,如Newman等基于非线性共轭梯度实现了横向各向异性的三维MCSEM反演[19];A.V. Grayver等使用了高斯牛顿法实现了三维MCSEM反演算法[20];Zach等实现了有限内存的BFGS (L-BFGS)的三维MCSEM反演算法[21]。然而,目前大多数的三维MCSEM反演研究都是基于“拖曳发射源 + 海底布设接收器”这种观测方式。因高效率、低成本而研发的深拖式海洋电磁装备在三维MCSEM反演研究还是比较少。Michael S. Zhdanov等采用重加权正则化共轭梯度(Reweighted Regularized Conjugate Gradient, RRCG)算法,对北海北部Troll油田的拖曳电磁数据开展了三维反演研究[22] [23];Jenny-Ann Malmberg等利用高斯牛顿法实现了拖曳电磁数据的三维反演[24];Joel Skogman等同样基于高斯牛顿法,对Barents海的拖曳电磁数据进行了三维反演分析[25];此外,Johan Mattsson等也开展了基于高斯牛顿算法的拖曳电磁数据三维反演工作[26]。
综上所述,现有的MCSEM三维反演研究在拖曳方式、算法优化以及观测系统等方面已取得了重要进展[27],但多数工作仍以海底接收器为主要观测模式,对深拖式MCSEM反演研究相对较少。深拖式方案在观测效率、作业灵活性以及连续测线覆盖方面具有显著优势,能够有效提升数据采集与横向成像效率。因此,开展针对深拖式频率域MCSEM的三维反演研究,不仅有助于完善该观测系统的理论基础和方法体系,也为深水油气资源勘探提供了新的技术支撑与实践参考。为此,本文开展了基于simpeg开源库[28] [29]的三维频率域MCSEM近似高斯牛顿法反演研究。首先介绍了基于树网格剖分的三维MCSEM拟态有限体积正演理论;随后提出了反演中采用的目标函数构建方法、正则化参数的自适应计算策略、步长因子的搜索方案;最后通过理论算例,验证了所构建反演程序的正确性与可行性。
2. 正演方法
2.1. 电磁场的控制方程
在频率域形式下,采用复指数时间因子
,以便将时间微分转化为频率算子
,针对非均质各向异性介质中时谐电磁场的三维可控源电磁正演问题,此时麦克斯韦方程组可表示为:
(1)
(2)
其中,
和
分别为电场强度与磁场强度,
与
分别为磁感应强度与电流密度;
表示磁源项,对应源磁通密度;
表示电源项,对应源电流密度。
为建立电磁场与介质参数之间的关系,引入本构方程:
(3)
(4)
其中,
为磁导率,
为电导率,
为外加体电流源。在MCSEM应用中,由于海洋电磁法的典型频率和介质条件下,
,这时位移电流
相对于传导电流可忽略,从而电磁场主要表现为扩散传播。因此常采用准静态近似,仅保留
。
2.2. 有限体积法
有限体积法,通常采用矩形剖分(图1),在矩形单元内,梯度(grad)、旋度(curl)、散度(div)等算子满足相应的积分恒等式和守恒律。同时,将电场通常定义在网格的边(edge)上,磁感应强度定义在网格的面(face)上,电流密度与磁场的散度则与体积(cell center)相关[30] [31]。这种基于边–面–体的物理量配置方式与麦克斯韦方程的微分算子结构相匹配,使得离散后的系统矩阵不仅稀疏,而且能够较好地保持稳定性与数值精度。
Figure 1. Location of variables within a cell (Lindsey J. Heagy et al., 2017)
图1. 单元内的变量位置(Lindsey J. Heagy等,2017)
为了对方程进行离散化,引入矢量测试函数
,并对方程(1)~(4)在计算域Ω上做内积。
(5)
(6)
(7)
(8)
其中(6)式中的曲面积分项在自然边界条件下为零,即
根据原则:离散电场
定义在网格边上,离散磁通密度
定义在网格面上。因此
必须定义在边上,
必须定义在面上。
是与电场
同离散位置(网格边的中点,“边变量”);
是与磁感应强度
同离散位置(网格面的中点,“面变量”),可得以下离散内积关系:
(9)
(10)
(11)
(12)
其中:
为离散旋度算子;
与
分别为积分形式的磁源与电源项;
为边上的内积矩阵;
为面上的内积矩阵;
为投影到边上的电导率内积矩阵;
为投影到面上的磁导率倒数的内积矩阵。
利用(9)~(12)式,整理得到求解电场的标准形式
(13)
其中
同理,可以得到磁场的标准形式
(14)
其中
式中,
为投影到面上的电阻率的内积矩阵;
为投影到边上的磁导率内积矩阵。
在求解大型稀疏矩阵系统(13)、(14)式时,为了在保证精度的同时提升求解效率,本文采用了基于直接法的高性能稀疏矩阵求解器PARDISO [32]。
3. 反演理论
3.1. 目标函数
为了获得海底地质结构的电性分布特征,反演是数据处理解释中关键环节。MCSEM的三维反演问题,本质上求解式(15)的目标函数。
(15)
式中:右边第一项为数据拟合项,用于衡量观测数据与模型正演响应之间的拟合程度。
是数据协方差矩阵,是一个对角矩阵,其元素等于
,
表示对第i个观测数据的标准差的估计值,它通常用来刻画第i个观测值的不确定性或误差大小。
是由反演模型通过正演得到的数据;
是真实模型的观测数据。第二项为模型约束项,可引入先验信息,对模型施加合理约束以提高反演稳定性。
是一个“正则化矩阵”,它不是单一的矩阵,而是有几个部分堆叠组合而成,这些部分分别代表
,
,
,
是模型在x,y,z方向上的一阶导数,
是单位矩阵,各部分前面的
,
,
,
是权重系数,用来平衡这些正则化项的重要性;
为参考模型。
为正则化参数。
依据实测数据采用高斯牛顿法求解式(15),使得目标函数最小,则可以反演得到模型的电阻率结构。
3.2. 正则化参数β的计算策略
正则化系数
的合理选取对于平衡数据拟合精度与模型稳定性至关重要,当
过大时,反演结果会过度依赖先验模型,而观测数据的影响被削弱;当
过小时,则可能引入不稳定或错误的解。本文通过Tikhonov正则化求取
。
具体方法为,先利用式(17)估计一个初始的正则化权重
,随后依据式(18)动态调整正则化权重。
(17)
其中:
是可以自定义的目标比率。
(18)
其中:
为冷却因子,它决定了每次冷却时
的缩减幅度;r为冷却频率,它控制每隔多少次迭代执行一次冷却操作。
4. 数值算例
为评估本文所提出三维MCSEM反演方法的可靠性与适用性,本节设计并开展数值算例研究。算例模型依照由简至繁的思路构建,检验算法的成像能力与稳定性。首先选取几何结构相对简单的单异常体模型,对反演流程的精度、参数设置及收敛特性进行基础验证;在此基础上,再引入多异常体模型,以进一步考察方法在相对复杂的条件下的适应性和稳健性。对于算例模型,利用电场分量Ex的实部和虚部进行反演,并在正演得到的数据中加入5%的高斯噪声。所有计算均使用处理器2.60 GHz、运行内存32 GB的计算机。
为在保证正演与反演计算精度的同时有效控制计算规模,本文采用三维八叉树(TreeMesh)网格剖分策略对研究区域进行离散化。首先构建基础网格,单元初始尺寸设定为x,y,z = 100,100,100 m,沿x、y、z三个方向分别划分为[64, 64, 64]个基元单元,共57,268个六面体单元,并以模型中心为原点(x0 = “𝐶𝐶𝐶”),从而获得覆盖整个研究域的均匀基础网格(如图2)。利用refine_tree_xyz()算法对全域实施一次全局细化,以形成具备多层次分辨率的初始网格体系,为后续局部加密提供基础。在基础网格的框架上,为兼顾近源异常体与远场关键层位的分辨需求,进一步采用区域约束的局部加密策略。其中,首先在[−1800, 1800] × [−1800, 1800] ×[−1000, −200] m的空间范围内执行一次局部细化,以提升模型中心及潜在异常体邻近区域的网格分辨率,从而增强电磁响应的刻画能力;然后,在[−3200, 3200] × [−3200, 3200] × [0, 50] m的范围内,对位于海底及其邻近区域的网格实施进一步加密。同时,x方向往外扩充2400 m (两端各1200 m),y方向往外扩充2400 m (两端各1200 m),z方向往外扩充4200 m (上端2200 m,下端2000 m)。
Figure 2. Model grid
图2. 模型网格
本次调查共布设11条测线,测线沿y方向分布在−1500至1500 m范围内,线距为300 m。每条测线上布设16个水平偶极源(源长为200 m,频率为1.0 Hz,z坐标固定为50 m),测量分两次进行:一次沿x正方向,源中心位置从0 m至2000 m,间隔200 m;另一次沿x反方向,源中心位置从0 m至−2000 m,间隔同为200 m。源牵引9个接收器,第一个接收器距源中心400 m,其余接收器间间隔200 m,接收器z坐标固定为30 m,如图3。整个观测中,共获得观测数据2178个。
Figure 3. Towed survey track
图3. 拖曳轨迹
4.1. 单异常体模型
设置空气层厚度为2200 m,电导率为1.0 × 10−12 S∙m−1;海水层厚度为1000 m,电导率为3.33 S∙m−1;异常体顶部距离海底400 m,电导率为0.25 S∙m−1,大小为1600 m × 1000 m × 200 m;海底背景电导率为1.0 S∙m−1 (图4)。
Figure 4. Horizontal seafloor model
图4. 水平海底模型
将海底的背景模型作为初始模型,采用本文方法反演的结果如图5所示。
图5清晰反映了正演模型与添加5%高斯噪声后数据反演结果的对比。在水平方向(x, y)及垂直方向(z)上,反演结果较好地恢复了异常体的位置及其主要电性分布,表明该方法具备较好的三维空间分辨能力。从反演收敛过程来看(图6),均方根误差(RMS)在迭代初期快速下降,随后趋于稳定,最终收敛至约1.4,接近理想拟合水平1.0。该过程表明反演算法具有良好的收敛性。进一步对异常体区域进行定量评价,反演电阻率值的平均相对误差为39.8%,表明在海水高导电性及有限观测几何条件约束下,反演结果能够有效刻画异常体的空间位置,但在电阻率幅值恢复方面仍存在一定偏差。
Figure 5. Single anomaly model (left: forward model, right: inverted model; a: slice at y = 0, b: slice at x = 0, c: slice at z = −500 m)
图5. 单异常体模型(左侧为正演模型,右侧为反演模型,其中a为y = 0的切片,b为x = 0的切片,c为z = −500的切片)
Figure 6. Data misfit
图6. 数据拟合差
4.2. 多异常体模型
在实际海洋地质环境中,地下介质通常具有复杂的非均匀性,常表现为多个空间位置的地质体共存。为更真实地模拟此类地质条件,并检验反演方法在复杂电性结构下的分辨能力与稳定性,本文设计了包含两个异常体的理论模型。
Figure 7. Multiple-anomaly model (left: forward model, right: inverted model; a: slice at y = 0, b: slice at z = −300 m, c: slice at z = −700 m)
图7. 多异常体模型(左侧为正演模型,右侧为反演模型,其中a为y = 0的切片,b为z = −300的切片,c为z = −700的切片)
模型设置一个异常体顶部距离海底200 m,另一个异常体顶部距离海底600 m,大小都为900 m × 1000 m × 200 m,电导率都为0.25 S∙m−1,其他设置与上一节模型均相同(图7左)。
反演结果如图7所示。
图中对比了真实模型与添加5%高斯噪声后的反演结果在Y = 0 m的垂直剖面以及Z = −300 m和Z = −700 m水平切片上的电性分布特征。整体而言,反演结果能够较好地恢复低电导率异常体的空间位置,表明反演方法在异常体空间位置重构方面具有一定有效性。在异常体电导率和边界刻画方面,反演结果仍与真实模型存在差异。在垂直剖面(Y = 0 m)上,反演能够识别异常体的埋深及横向分布,但较深的异常体电导率恢复明显减弱,边界呈现平滑扩散特征,并在异常体下方局部区域产生高电导率伪异常,反映出垂向分辨率随深度增加而降低。异常体区域电阻率值的平均相对误差为56.2% (较浅异常体区域电阻率值的平均相对误差为50.0%,较深异常体区域电阻率值的平均相对误差为62.4%),定量评价进一步说明反演在深部异常体幅值恢复上存在更大偏差。综合分析认为,该反演结果在异常体定位方面具有一定可靠性,但对深部异常体电性幅值及精细结构的分辨能力仍然不足,这体现了三维海洋电磁反演在有限观测条件下的固有分辨率限制。
Figure 8. Data misfit
图8. 数据拟合差
从图8可知,RMS在迭代后期逐渐至1.4,趋近于理论理想值1.0,说明了反演具有良好的收敛性。整个收敛过程平稳可控,进一步验证了该反演算法在多异常体复杂模型中依然具有一定的可靠性与实用性。
4.3. 结论
本文通过三维深拖式海洋电磁正演与反演模拟,对单异常体与双异常体模型条件下的反演成像能力进行了对比分析。结果表明,该反演方法能够有效恢复异常体的三维空间位置与整体形态,在由单异常体向双异常体变化的情况下仍保持较好的稳定性与成像能力。
对于单异常体模型,反演结果在异常体定位、尺度及电性分布方面与真实模型高度一致,显示出较高的反演精度;而在双异常体模型中,反演能够成功区分多个异常体的空间展布,但异常体电导率值有所削弱,边界趋于平滑,异常体下方区域局部出现伪影,反映出复杂电性结构条件下电性结构刻画能力的降低。
反演过程中均方根误差(RMS)在迭代初期快速下降并最终稳定收敛,表明算法具有良好的收敛性。总体而言,该方法在异常体几何成像方面表现可靠,但对深部及复杂结构的精细电性刻画仍有改进空间,后续可通过多频联合反演与观测系统优化进一步提升成像质量。
本文在数值反演中假设介质均为各向同性,忽略了实际中可能存在的各向异性。这一假设可能导致反演异常体电阻率值和形态与实际存在一定偏差,尤其在层状或裂隙导电方向明显的区域。未来研究可考虑引入各向异性模型以更准确地恢复地下介质结构。
NOTES
*通讯作者。