基于形状约束的多道地震反演方法对空腔水库的表征 李凯、 黄旭瑞、 周丹、 陈晓春、 牛丽萍

《Eng》：Characterization of Cavity Reservoirs Based on a Shape-Constrained Multi-Trace Seismic Inversion Method Kai Li, Xuri Huang, Dan Zhou, Xiaochun Chen and Liping Niu

　　摘要：中国西北部的超深海洋碳酸盐岩洞穴储层具有强烈的异质性和复杂的几何形态。传统的地震反演方法会产生过于平滑的结果，模糊地质边界，并且存在严重的非唯一性问题，这使得准确识别低阻抗洞穴异常变得困难。为了解决这些问题，本研究开发了一种基于Mumford–Shah泛函的形状约束多迹地震反演方法。为了满足计算要求，采用了多迹反演框架，并引入了Ambrosio–Tortorelli近似方法将不光滑的Mumford–Shah泛函转换为可求解的平滑形式。所提出的方法实现了声学阻抗和地质界面的联合反演，隐式地鼓励出现分段常数区域，并保持了清晰的地质边界，有效缓解了反演的非唯一性问题。数值实验和现场应用验证了该方法即使在有限的先验信息下也能提供高横向分辨率和边界精度，并且能够可靠地描绘出与钻井和测井数据高度一致的离散“串珠状”洞穴几何形态，为复杂碳酸盐岩洞穴储层的精细表征提供了可靠的解决方案。

1. 引言
碳酸盐岩储层在中国西北部的塔里木盆地广泛分布，这些储层的埋藏深度超过6500米[1]。沿着走滑断层带，多阶段构造作用促进了碳酸盐岩层中洞穴、孔隙和裂缝型储层的发展[2,3]。特别是洞穴储层，给勘探带来了巨大的挑战。这类洞穴储层是实现高产油气井的关键钻探目标[4]。经过强烈的后期改造，碳酸盐岩洞穴储层具有强烈的异质性、多样的储层类型、离散的空间分布以及各种组合[5]。这些特征导致目标层的地震响应非常复杂，给储层解释带来了巨大困难。
利用一种称为“串珠响应（SBRs）”的特征地震反射模式，地震成像可以识别碳酸盐岩层中的喀斯特洞穴发育范围[6]。然而，由于地震子波的带宽有限以及反射事件之间的干扰，地震成像缺乏关于洞穴形态、边界和组合类型的详细信息。在勘探的早期阶段，广泛使用了地震属性技术[7]。这种方法通常结合应用多种地震属性，包括地震相干性、反跟踪属性、梯度结构张量（GST）、反射强度和均方根（RMS）振幅，并在断层控制型洞穴储层预测中取得了许多成功的现场应用[8,9,10,11]。属性方法需要更高质量的地震数据（信噪比和分辨率），主要用于定性描述断层控制型储层中裂缝区的范围。此外，基于振幅的属性无法解决地震响应中由子波引起的模糊效应，而储层边界的识别需要基于通过钻井验证的已知储层的统计属性值来设定阈值[12]。
基于地震数据的一维稀疏脉冲反演和先验约束已被广泛用于估计地下阻抗模型，并也应用于碳酸盐岩储层。然而，传统的地震反演方法主要依赖于相对均匀的低频模型来预测储层参数[13,14,15]。在实际勘探中，构建可靠的碳酸盐岩洞穴储层先验模型存在重大困难。碳酸盐岩洞穴储层的强烈异质性导致其地震响应中的横向能量变化显著。此外，在钻井作业过程中，当钻柱穿过目标层中的储层时，经常会发生严重的循环损失。这些井下异常不仅影响测井绳的操作稳定性，还导致用于储层评估的有效测井曲线数据严重不足[16,17]。这些因素导致构建的先验模型不准确。当应用传统反演方法时，如果先验约束不足，反演剖面中会出现大量孤立且横向突变的条带状异常，这是由于单迹地震反演方法的固有限制所致。相反，基于不可靠先验模型的强约束可能会掩盖低阻抗异常，从而降低油气勘探的效果[18]。后续的地质统计反演和相控反演也应用于碳酸盐岩洞穴储层。然而，这些方法仍然严重依赖于地震属性或地震相的质量来约束反演，导致反演结果的一致性较差[19,20]。此外，全波形反演（FWI）和人工智能的最新发展也推动了地震反演和裂缝洞穴表征的进步。然而，FWI在计算上非常耗时，特别是对于超深碳酸盐岩储层，因为高分辨率反演需要极大的计算资源[21]。尽管基于机器学习的反演方法在传统储层方面取得了令人鼓舞的改进，但由于缺乏足够的实际现场数据用于训练，这仍然是限制其应用于超深洞穴储层的关键因素[12,22,23]。因此，这些方法在实际超深碳酸盐岩洞穴数据集上的性能仍需进一步研究和开发。
为了解决单迹反演的局限性，最近的研究越来越多地关注多迹阻抗反演。与单迹地震反演不同，多迹同步反演的固有特性避免了反演结果中的虚假横向异常[24,25]。同时，多迹反演将2D或3D模型作为输入和输出，使我们能够结合空间约束或结构正则化来实现所需的反演性能[26]。近年来，多迹反演方法在解决非平稳地震效应和正则化挑战方面取得了显著改进。先进的多迹反演结合了结构约束，实现了衰减参数和反射率的同时估计，提高了时间分辨率和空间结构连续性[27]。还开发了自适应乘法正则化策略用于多迹反演，实现了正则化强度的自动调整，避免了低效的手动参数调整[28]。这些发展促进了在复杂地质环境中获得更稳定、高分辨率和结构一致性的阻抗结果。
从方法论本质的角度来看，地质统计反演侧重于再现统计空间变异性，相控反演强调在预定义的相边界下的地质约束一致性，而多迹反演旨在通过利用相邻迹之间的结构相似性来增强横向结构连续性。尽管这三类方法提高了反演结果的稳定性和地质合理性，但它们本质上追求连续或平滑变化的参数场，这与碳酸盐岩洞穴储层的离散、孤立和强烈异质性特征本质上是不符的。
对于洞穴型储层，我们的目标是突出与离散高能量“串珠响应”相对应的低阻抗异常，并恢复其准确的几何形态。这项任务在数学上类似于图像分割，后者旨在识别分段常数区域和不连续边界[29]。与传统反演范式相比，图像分割提供了一种自然且强大的替代方法，完美符合表征清晰、孤立洞穴储层的目标。
作为图像分割的经典变分模型，Mumford–Shah（MS）泛函可以同时实现区域内的平滑和区域间的边界保持，从而在地震反演和洞穴形状表征之间提供了严格的理论桥梁。图像分割是计算机视觉中的一个基本课题，通常使用模型参数的梯度信息来划定不同结构单元的边界。作为自由不连续问题的代表性变分模型，MS泛函通过捕捉内部结构和界面特征实现了图像分割、去噪和重建[30]。这种类比激励我们将MS泛函引入多迹地震反演。
在本文中，我们提出了一种基于MS泛函的形状约束多迹地震反演方法。形状约束反演方法的关键创新在于同时反演界面和参数值。基于MS泛函框架实现了隐式形状表示。在这个框架内，隐式形状表示仅在真实界面处形成清晰的过渡层，从而避免了传统反演方法为减少数据不匹配而牺牲结构清晰度的常见缺陷。通过隐式鼓励出现分段常数区域，该框架有效地减少了优化过程中的过度自由度。
通过施加形状约束，强烈非均匀储层的反演结果实现了更高的空间分辨率和储层预测精度。我们首先使用合成数据验证了所提出的方法。数值测试表明，当先验信息有限时，形状约束反演在结构重建方面具有显著优势，并且也能提供更高的数值精度。我们进一步将该方法应用于中国西北部碳酸盐岩洞穴储层上采集的实际地震数据。现场应用表明，所提出的方法比传统方法能够更好地表征洞穴型储层。
为了明确总结本工作的主要贡献，以下突出了三个关键创新：首先，提出了一种嵌入MS泛函的形状约束多迹反演策略，以实现地质界面和地下参数的联合反演。其次，引入了隐式形状表示，以应对碳酸盐岩洞穴储层的强烈非均匀性，确保了边界清晰度的同时降低了优化复杂性。第三，所提出的方法有效地避免了传统方法中存在的结构清晰度与数据不匹配之间的权衡，提高了高分辨率储层表征的实用性和可靠性。

2. 理论与方法
在本章中，我们首先介绍了多迹地震反演的基本原理。在此基础上，我们将Mumford–Shah泛函作为正则化项，建立了形状约束多迹反演的目标函数。我们还介绍了相应的解决方法以及应用所提出方法的工作流程。

2.1. 多迹阻抗反演
单迹叠后地震迹遵循卷积模型[31]：
(1)
其中 是单迹地震记录的向量， 是小波矩阵， 是反射率序列的向量， 表示噪声信号。方程中的符号*表示卷积操作。我们将时间域中某一层的阻抗表示为 ，其相邻下一层的阻抗表示为 ，然后反射率 可以表示为[32]：
(2)
一般来说，相邻层之间的声学阻抗对比度较弱。为了便于方程的线性化，可以近似阻抗与反射系数之间的关系[33]：
(3)
因此，在无噪声条件下，方程(1)可以重写为矩阵形式：
(4)
其中 是地震振幅， 表示小波， 是小波的长度， 是前向建模算子， 是阻抗向量。方程(4)是一个单迹前向函数，我们可以通过将向量和 展开为二维矩阵来获得多迹前向函数：
(5)
其中 是迹数， 是数据的第m个地震信号迹， 是模型的第m个阻抗迹。因此， 是地震信号矩阵， 是阻抗模型矩阵，它们都包含 个迹。基于多迹前向建模算子，我们可以通过最小二乘法轻松推导出多迹反演的目标函数：
(6)
其中 是观测地震数据， 表示L2范数。梯度通过对方程求导得到：
(7)
在本研究中，采用多迹反演有两个主要原因。首先，同时进行多迹反演可以提高反演结果的横向一致性。其次，随后引入的形状约束反演方法无法在单迹框架中有效捕获模型的结构信息。因此，需要采用多迹反演策略。

2.2. 形状约束反演
正则化对于稳定地震反演和施加地质合理的地下参数分布约束至关重要。Tikhonov L2正则化结合初始模型约束是最广泛使用的正则化策略之一[34]：
(8)
其中 是模型参数， 是先验模型， 是控制约束强度的权重因子。
从物理上讲，这个项迫使反演模型符合先验模型并抑制数值振荡。然而，它在整个模型域上强制全局平滑，这与超深碳酸盐岩洞穴储层的物理特征直接矛盾：离散的空间分布、洞穴边界处的强烈阻抗不连续性以及非层状沉积环境。此外，在超深裂缝碳酸盐岩中，井壁坍塌是常见的现象，导致井壁粗糙度显著增加或井径明显扩大。这些不利的钻井条件导致密度、压缩波地震和中子测井数据缺失，使得模型缺乏足够的先验信息[18]。
为了解决这种物理不匹配，我们引入了MS泛函——这是一种最初为图像分割开发的变分模型，它自然地保留了结构边界处的不连续性。从数学和地球物理学的角度来看，碳酸盐岩洞穴与围岩之间的清晰阻抗边界类似于图像处理中的界面。MS泛函将地下结构划分为由清晰界面分隔的独立地质体，在每个地质体内强制参数变化平滑，同时允许界面之间的强烈对比。这种设计与实际腔体储层的分布在物理上是一致的。MS泛函定义为[30]：（9）其中是待更新的数据，是观测数据，是目标模型的空间域，是界面集。MS泛函可以分解为三个贡献项。方程（9）中的第一项用于减少模拟数据与观测数据之间的差异。第二项是对模型空间中除界面集之外的参数变化的惩罚。第三项积分项测量弧长并惩罚界面的振荡。两个系数和作为权重因子，用于平衡最后两项的贡献。MS泛函将模型划分为不同的子区域，强制每个区域内的参数平滑，并规范界面的长度[30]。换句话说，当应用于地下参数反演时，MS泛函有助于将域分割成多个独立的区域及其相应的界面。关键的是，它允许不同区域之间的参数有显著的空间变化，从而防止与腔体储层相关的低阻抗异常被抑制。然而，由于界面集未知，最小化方程（9）是困难的。尽管Mumford–Shah泛函为形状约束反演和边界保持提供了一个理想的数学框架，但其原始形式涉及未知的界面集K且不光滑，这使得在基于梯度的多道反演框架内直接求解变得困难。对于多道地震反演，我们需要一个光滑、可微分且数值稳定的泛函，可以使用基于梯度的迭代方法高效优化。为了克服这一困难并在多道框架中实现实际应用，我们引入了Ambrosio–Tortorelli近似。这种近似将非光滑的自由不连续问题转化为光滑的椭圆变分问题，允许在多道环境中联合迭代更新阻抗和界面。没有这种近似，Mumford–Shah泛函无法在多道计算流程中稳定且高效地最小化。

Ambrosio和Tortorelli用界面的光滑估计替换了界面集，并基于方程（9）开发了这种近似，称为Ambrosio–Tortorelli近似[35]。根据Ambrosio–Tortorelli近似的理论框架，我们假设MS泛函中涉及的Hausdorff长度可以在欧几里得范数下定义，基于此我们推导出所提出的形状约束正则化项[35,36]：（10）当的值接近0时，它对应于地质界面的位置；当的值接近1时，它对应于界面之间的光滑区域。因此，它将非光滑泛函转化为适合基于梯度求解的光滑椭圆泛函。是一个小的正常数，用于控制界面的光滑度。显然，当对所有时，方程（10）中的正则化项简化为Tikhonov正则化。方程（10）关于和分别是凸的。因此，在每次迭代中更新时，可以通过微分明确获得的形式解：（11）此外，界面的输出可以通过参数、和来控制：控制地质体内的内部平滑强度；控制反演界面的尺度和复杂性；控制边界过渡的锐度。

图1a展示了一个简单模型。该模型由两个正方形子单元组成。一个正方形的内部参数值为2，另一个为1，背景的值为10。这两个正方形部分重叠，它们界面处的数值差异小于与背景参数的差异。为了将这个模型分割成子区域和界面，我们使用方程（10）中描述的形状约束正则化。图1b显示了在最优参数设置下获得的界面结果。随着的增加，界面集的长度变小，即模型内部的分割尺度变大。随着的增加，恢复的界面变得更加模糊。通过选择适当的参数，我们可以获得一种自然的多尺度策略，以提高反演稳定性。图1. 简单的数值模型和用不同权重参数获得的界面集：(a) 简单的数值模型；(b) 用小值和获得的界面集；(c) 用小值和大值获得的界面集；(d) 用大值获得的界面集。一般来说，如果我们希望获得更多的模型细节（小尺度信息），我们需要保持和的值较小；然而，应该注意的是，过小的也会影响分割的稳定性。

在这里，我们介绍了参数值选择的策略。第一个参数涉及，它本质上是对分割区域内部应用的L2范数约束项的权重系数。其选择可以完全借鉴Tikhonov正则化参数调整的经验。当地震振幅归一化到-1到1的范围时，最优的alpha通常落在到区间内。与L1型正则化相比，L2范数约束的权重参数不太敏感，消除了对极细参数调整的需求。对于Ambrosio–Tortorelli（AT）近似，通常需要大于模型网格间隔，以确保优化过程中的数值稳定性。然而，在波阻抗反演中，梯度差异通常是在不考虑物理距离维度的情况下计算的，因此可以将模型网格间隔视为1。因此，建议将设置在1到10的范围内。除非反演问题严重不适定，否则大于10的ε值会导致边界过于模糊，使得形状约束无效。对于参数的选择，其值通常与模型的大小和地震数据的大小无关。它控制施加在边界长度上的惩罚。参数的敏感性在0.0001到0.001的范围内最高。当超过0.01时，敏感性降低，进一步增加只会导致结构特征的过度平滑和模糊。最后，应该注意的是，在选择这些参数时，方程10中第一项和第二项的平均能量应尽可能保持平衡。第一项的能量过大会导致过度平滑的反演结果，而第二项的能量过大则会导致由噪声引起的过多虚假界面。

MS泛函的核心优势在于它允许在这些边界处存在不连续性，而传统的Tikhonov正则化则强制解决方案的全局平滑性。通过上述介绍，我们可以得出结论，将这种约束纳入反演中可以实现对地质边界（如腔体和断层）的多尺度表征。最后，通过结合多道地震反演和形状约束正则化，我们得到目标函数：（12）为了最小化方程（12），可以使用共轭梯度、牛顿法和最速下降法等优化方法进行优化，从而求解和。通过方程（11）在给定迭代中获得的界面构成了下一次迭代中更新阻抗模型的隐式形状表示。

这种策略的核心物理优势在于它将反演从基于像素的网格优化重新定义为基于地质体的反演。地质体的形状通过MS泛函解决的界面集在数学上表示。由于形状不是显式参数化的，而是由边界隐式定义的，所提出的框架可以定义为隐式形状约束反演。形状约束多道地震阻抗反演的工作流程可以概括如下（图2）：图2. 形状约束地震反演的工作流程。(1)地震小波的估计。在合成记录的情况下，我们可以简单地使用用于观测数据建模的小波，例如Ricker小波。然而，在实际情况下，我们通常在目标层使用统计小波或通过井到地震校准来估计小波。(2)初始模型构建。使用先验信息（如测井和地震解释）来构建一个平滑的低频模型。(3)界面集的初始化。通常建议将界面集全部初始化为0，以确保稳定的更新。然而，如果有可靠的先验模型可用，初始界面也可以根据地质理解来设置。(4)求解目标函数。首先，设置较大的值并迭代求解以获得大尺度反演结果，然后基于大尺度解，给出较小的值以获得精细解。应该注意的是，方程（12）是双凸的，因此需要交替迭代来分别更新模型和界面。

3. 合成数据示例
在本节中，我们将形状约束多道反演方法（SCI）应用于一个数值模型案例，以验证所提方法的有效性。数值模型（图3a）是从Marmousi模型的一小部分提取的，包含复杂的陡峭结构，这可以更好地展示形状约束的优势[37]。无噪声的合成地震数据是使用前一节中的方程（5）和25 Hz峰值频率的Ricker小波计算得出的。然后我们在无噪声数据中添加10%的高斯随机噪声以获得观测数据（图3b）。模型和合成数据都包含201道迹线，每道迹线有51个样本，采样率为2毫秒。我们选择了真实模型的第121道迹线作为伪测井数据，然后对其进行平滑并横向扩展以构建初始模型（图3c）。为了与SC-I进行比较，图3d显示了使用先验模型约束的传统多道反演（PCI）的结果。用于反演的小波与正向小波一致，并采用了共轭梯度方法，进行了10次迭代。可以观察到，传统PCI方法得到的结果表现出强烈的平滑效果，陡峭结构的分层细节明显模糊，这对结构分析和解释不利。图3. 数值模型及其传统的先验约束反演结果：(a) 数值阻抗模型；(b) 带有10%高斯随机噪声的合成记录；(c) 通过伪井数据横向外推获得的初始模型；以及(d) 通过传统先验约束反演获得的阻抗模型。第二种方法是本研究中提出的形状约束地震反演（SCI），通过图2所示的工作流程同时获得阻抗和界面。在SCI过程中，小波和迭代次数与PCI保持一致。具体来说，大尺度SCI迭代了五次，小尺度SCI又迭代了五次，总共十次迭代。

为了确保反演的稳定性，我们首先进行大尺度反演，并将形状约束正则化的权重因子设置为（由地质特征的特征尺度指导）。大尺度反演得到的阻抗和界面集结果如图4所示。图4. 大尺度SCI得到的阻抗和界面集结果：(a) 通过大尺度形状约束地震反演得到的阻抗模型；(b) 通过大尺度形状约束地震反演得到的界面集。在图4a中，具有大梯度的界面非常清晰。然而，由于我们对总界面长度的严格约束（由权重因子控制），具有小梯度的阻抗界面被平滑了。图4b展示了反演后的界面集，与阻抗模型吻合良好。在这个过程中，界面集不仅作为目标模型，还作为阻抗更新的约束。

尽管大尺度反演模型不是最终输出，但这个大尺度反演步骤通常是必要的。直接进行小尺度反演可能导致反演结果不稳定和模型过度锐化，从而在最终输出中出现伪影。我们将大尺度反演结果作为初始模型，并进行小尺度SCI（Shape-Constrained Inversion），权重因子设置为特定值。经过五次迭代后，小尺度SCI的结果如图5所示。图5展示了通过小尺度SCI获得的阻抗和界面集结果：(a) 通过小尺度形状约束地震反演得到的阻抗模型；(b) 通过小尺度形状约束地震反演得到的界面集。图5a中显示的反演阻抗在结构上与真实模型非常接近，模型的小梯度细节也被很好地刻画出来。这是因为我们放宽了对界面长度欧几里得度量的惩罚，从而在小尺度反演中恢复了更多的界面细节，如图5b所示。与通过PIC（Particle Inversion）方法得到的阻抗模型（图3d）相比，SCI的反演结果几乎恢复了每个陡倾结构的阻抗界面，例如最左侧60毫秒处的薄层。同样，小尺度界面集也与阻抗模型一致，这也有助于解释模型结构。此外，尽管PCI方法包含了Tikhonov正则化项，但其反演结果在模型上方的高阻抗区域仍然产生了更多的人为伪影。这是因为在这个模型测试中提供的先验信息不可靠，且Tikhonov正则化项的权重设置得不够大。增加权重将进一步过度平滑模型，使其更接近初始模型中给出的错误结构信息。需要注意的是，PCI方法在此实验中采用了多道反演策略，如果使用更常见的单道反演，结果会更差。相比之下，模型约束正则化在保留反演界面集的结构信息的同时，确保了模型的内部平滑性。

为了进一步定量评估不同反演策略对声学阻抗恢复精度的影响，我们从模型中选取了60毫秒处的水平切片进行详细数值比较。图6中，红色曲线代表本文提出的形状约束反演结果，蓝色曲线对应于传统的先验约束反演结果，黑色曲线表示真实的声学阻抗作为参考基准。从图6可以看出，SCI方法不仅在结构信息上具有优势，而且在反演阻抗方面也提供了更高的数值精度。我们分别计算了大尺度SCI、小尺度SCI和PCI反演结果的相对误差，相应的结果见表1。图6展示了使用不同反演策略在60毫秒处得到的时间切片与真实阻抗的比较。基于此，我们还在一个简单的腔体水库数值模型上进一步测试了所提出的方法。在这个测试中，我们将SCI（形状约束反演）方法的腔体反演性能与总变分（TV）正则化反演进行了比较，后者在保持边缘的地震反演中广泛应用。如图7所示，各面板分别展示了真实阻抗模型、使用25赫兹峰值频率Ricker小波进行正向建模生成的合成地震数据、反演的初始模型，以及基于TV正则化的波阻抗反演得到的阻抗模型。从TV正则化反演结果可以看出，水库的内部细节被过度平滑，只有大尺度的地层界面被保留下来，并且引入了明显的带状伪影。TV正则化本质上倾向于分段常数解，产生与水平或垂直方向对齐的阶梯状边界，这使得它无法准确捕捉腔体水库特有的局部不规则阻抗对比。仔细观察放大的局部细节，可以发现腔体的内部边缘明显“模糊”，未能恢复清晰明确的边界特征。此外，由于该模型测试中提供的先验信息不可靠，且Tikhonov正则化项的权重设置得不够大，因此其反演结果在模型上方的高阻抗区域仍产生了更多的人为伪影。值得注意的是，PCI方法在此实验中采用了多道反演策略，如果使用更常见的单道反演，结果会更差。相比之下，模型约束正则化在保留反演界面集的结构信息的同时，确保了模型的内部平滑性。

为了进一步定量评估不同反演策略对声学阻抗恢复精度的影响，我们从模型中选取了60毫秒处的水平切片进行详细数值比较。图6中，红色曲线代表本文提出的形状约束反演结果，蓝色曲线对应于传统的先验约束反演结果，黑色曲线表示真实的声学阻抗作为参考基准。表1展示了使用不同反演方法得到的阻抗结果的相对误差比较。在此基础上，我们在一个简单的腔体水库数值模型上进一步测试了所提出的方法。在这个测试中，我们将SCI（形状约束反演）方法的腔体反演性能与总变分（TV）正则化反演进行了比较，后者在保持边缘的地震反演中广泛应用。如图7所示，各面板分别展示了真实阻抗模型、使用25赫兹峰值频率Ricker小波进行正向建模生成的合成地震数据、反演的初始模型，以及基于TV正则化的波阻抗反演得到的阻抗模型。从TV正则化反演结果可以看出，水库的内部细节被过度平滑，只有大尺度的地层界面被保留下来，并且引入了明显的带状伪影。TV正则化本质上倾向于分段常数解，产生与水平或垂直方向对齐的阶梯状边界，这使得它无法准确捕捉腔体水库特有的局部不规则阻抗对比。通过检查放大的局部细节，可以发现腔体的内部边缘明显“模糊”，未能恢复清晰明确的边界特征。此外，由于与裂缝-腔体体相关的地震道强梯度，高阻抗的围岩也被TV正则化强制分割成人工段。虽然TV正则化在处理大尺度、规则的地层界面方面表现出色（例如水平分层介质），但在应用于孤立的、局部的腔体水库时存在固有局限性。

接下来，我们将形状约束反演应用于中国西北部一个盆地的真实地震数据，以表征超深奥陶纪碳酸盐腔体水库。在这种情况下，目标层是中奥陶纪海洋碳酸盐地层，具有典型的喀斯特特征；目标层的平均深度高达6000米。受热液溶解和多相构造运动的影响，在目标层段形成了空间离散的腔体型水库[38,39]。如图9所示的目标区域局部地震剖面图，地震事件轴是不连续且无序的，这使得构建有效的先验约束变得困难。此外，该地区的井在钻探到腔体顶部时会出现钻井中断或循环丢失的情况，导致只有目标层段的上部有声波测井数据，这对构建初始模型构成了挑战。图9显示了研究区域目标层段的时间域地震剖面图。紫色曲线代表目标层段的顶部和底部。在此应用中，我们简单地垂直扩展低频测井阻抗曲线，并沿目标地层的顶部和底部进行插值，以构建反演的初始模型。

为了进行比较，我们首先使用工业应用中常见的稀疏尖峰反演方法计算了相对阻抗。其反演结果如图10所示。稀疏尖峰反演未能识别腔体，并显示出强烈的非唯一性。图中标记为“AI”的曲线代表该位置A15井的声学阻抗曲线；这口井并未包含在初始模型中。可以看出，稀疏尖峰反演的结果与测井曲线中低阻抗异常的位置存在差异。此外，由于横向平滑和初始模型约束的影响，整体阻抗变化较小，使用阈值方法难以区分腔体与其他地质特征。图10显示了使用单道稀疏尖峰反演方法得到的目标层段的时间域阻抗剖面图。标记为“AI”的曲线是A15井的声学阻抗测井曲线，其值使用与阻抗剖面相同的颜色图例填充。随后，我们使用形状约束策略对这个地震数据集进行反演。反演中使用的初始模型和小波与稀疏尖峰反演的一致。图11展示了SCI得到的阻抗反演结果，图12显示了界面集的反演结果。与稀疏尖峰反演相比，SCI结果有效地突出了对应于腔体水库的低阻抗异常，并抑制了传统反演中的额外分层伪影。SCI结果能够清晰地识别腔体几何形状，这与强幅度SBRs的空间离散分布一致。同样，反演得到的界面集也能更直观地划分出水库边界。此外，使用A15井的实测声学阻抗曲线进行验证表明，SCI的准确性高于传统反演方法。SCI反演识别出的低阻抗异常位置与测井数据的一致性更高。图11显示了使用形状约束反演方法得到的目标层段的时间域阻抗剖面图。标记为“AI”的曲线是A15井的声学阻抗测井曲线，其值使用与阻抗剖面相同的颜色图例填充。图12显示了使用形状约束反演方法得到的目标层段的时间域阻抗剖面图。图13展示了深度域中SCI阻抗体积的inline1225剖面图，其中蓝线表示A15井井路径在剖面上的投影。由于在井底发生循环丢失事件，A15井的钻探被暂停。在腔体水库中，钻探过程中的循环丢失表明井筒已经遇到了裂缝-腔体系统；因此，这类事件是识别水库的关键指标[18]。如图11所示，沿井轨迹的钻探暂停点（循环丢失事件的位置）与SCI方法表征的腔体顶部完全吻合。图13显示了通过时深转换转换为深度域的inline1225 SCI阻抗剖面图。蓝线表示A15井井路径在剖面上的投影，由于在轨迹末端发生循环丢失事件，钻探被终止。沿井轨迹的曲线是测量的孔隙度测井曲线。此外，我们还应用所提出的SCI方法对腔体水库模型的合成数据进行了反演，相应的反演结果如图8所示。与TV正则化反演结果相比，SCI方法得到的腔体水库阻抗模型有了显著改进。SCI结果中的水库边界清晰明确，与腔体的真实轮廓高度吻合，同时完全保留了腔体内的低阻抗特征，使得能够清晰识别水库的内部结构。这表明SCI方法有效地平衡了保留大尺度地层界面和局部异常的特征，实现了边缘保留与细节恢复之间的良好权衡。

此外，虽然TV正则化在噪声抑制过程中不可避免且不可逆地减弱了与水库边界和内部结构相关的高频信息，但SCI方法利用形状约束将反演的自由度集中在局部水库特征上。这使得该方法能够有效地恢复低阻抗腔体内的高频阻抗变化，从而保留了水库边界的清晰度和内部阻抗分布的真实性。这些结果全面展示了所提出的SCI方法在应用于离散分布的腔体水库时的优越性能。更重要的是，该方法有效突破了传统保持边缘反演技术在表征小尺度异常体方面的固有局限性。

4. 实际数据应用
接下来，我们将形状约束反演应用于中国西北部一个盆地的真实地震数据，以表征超深奥陶纪碳酸盐腔体水库。在这种情况下，目标层是中奥陶纪海洋碳酸盐地层，具有典型的喀斯特特征；目标层的平均深度达到6000米。受热液溶解和多相构造运动的影响，在目标层段形成了空间离散的腔体型水库[38,39]。如图9所示的目标区域局部地震剖面图，地震事件轴是不连续且无序的，这使得构建有效的先验约束变得困难。此外，该地区的井在钻探到腔体顶部时会出现钻井中断或循环丢失的情况，导致只有目标层段的上部有声波测井数据，这对构建初始模型构成了挑战。图9显示了研究区域目标层段的时间域地震剖面图。紫色曲线代表目标层段的顶部和底部。在此应用中，我们简单地垂直扩展低频测井阻抗曲线，并沿目标地层的顶部和底部进行插值，以构建反演的初始模型。

为了进行比较，我们首先使用工业应用中常见的稀疏尖峰反演方法计算了相对阻抗。其反演结果如图10所示。稀疏尖峰反演未能识别腔体，并显示出强烈的非唯一性。图中标记为“AI”的曲线代表该位置A15井的声学阻抗曲线；这口井并未包含在初始模型中。可以看出，稀疏尖峰反演的结果与测井曲线中低阻抗异常的位置存在差异。此外，由于横向平滑和初始模型约束的影响，整体阻抗变化较小，使用阈值方法难以区分腔体与其他地质特征。图10显示了使用单道稀疏尖峰反演方法得到的目标层段的时间域阻抗剖面图。标记为“AI”的曲线是A15井的声学阻抗测井曲线，其值使用与阻抗剖面相同的颜色图例填充。随后，我们使用形状约束策略对这个地震数据集进行反演。反演中使用的初始模型和小波与稀疏尖峰反演的一致。图11展示了SCI得到的阻抗反演结果，图12显示了界面集的反演结果。与稀疏尖峰反演相比，SCI结果有效地突出了对应于腔体水库的低阻抗异常，并抑制了传统反演中的额外分层伪影。SCI结果能够清晰地识别腔体几何形状，这与强幅度SBRs的空间离散分布一致。同样，反演得到的界面集也能更直观地划分出水库边界。此外，使用A15井的实测声学阻抗曲线进行验证表明，SCI的准确性高于传统反演方法。SCI反演识别出的低阻抗异常位置与测井数据的一致性更高。图11显示了使用形状约束反演方法得到的目标层段的时间域阻抗剖面图。标记为“AI”的曲线是A15井的声学阻抗测井曲线，其值使用与阻抗剖面相同的颜色图例填充。图12显示了使用形状约束反演方法得到的目标层段的时间域阻抗剖面图。图13展示了深度域中SCI阻抗体积的inline1225剖面图，其中蓝线表示A15井井路径在剖面上的投影。由于在井底发生循环丢失事件，A15井的钻探被暂停。在腔体水库中，钻探过程中的循环丢失表明井筒已经遇到了裂缝-腔体系统；因此，这类事件是识别水库的关键指标[18]。如图11所示，沿井轨迹的钻探暂停点（循环丢失事件的位置）与SCI方法表征的腔体顶部完全吻合。图13显示了通过时深转换转换为深度域的inline1225 SCI阻抗剖面图。蓝线表示A15井井路径在剖面上的投影，由于在轨迹末端发生循环丢失事件，钻探被终止。沿井轨迹的曲线是测量的孔隙度测井曲线。此外，我们还提取了A15井轨迹上的SCI阻抗进行与测井曲线的比较（图14）。图14中，MD表示井的测量深度，LN_RT表示真实电阻率的自然对数，POR代表孔隙度。在超深碳酸盐腔体水库中，低声学阻抗表明裂缝和洞穴发育广泛，这通常对应于较高的孔隙度，因为地层内的空隙空间增加。同时，这样的腔体区域通常充满烃类或低导电性流体，导致高电阻率响应。因此，低阻抗、高孔隙度和高电阻率之间存在明显的对应关系，这可以作为识别有效腔体水库的可靠标准。图14显示了SCI阻抗曲线与A15井轨迹上测井曲线的比较。MD表示测量深度，LN_RT表示真实电阻率的自然对数，POR表示测量孔隙度，SCI阻抗表示沿井轨迹提取的反演结果。可以观察到，SCI阻抗曲线中的低阻抗区间（对应于图13中特征为SCI的腔体顶部）与测量孔隙度对数中的高孔隙度区间相吻合。同时，LN_RT对应于高电阻率区域，这些区域被解释为有效的储层。同样，我们还展示了inline 1714的深度域SCI阻抗剖面，如图15所示，其中蓝线代表Well A12的井轨迹。在Well A12的钻井过程中，在井轨迹的末端发生了循环损失和钻井中断，表明井底与腔体储层相连。井轨迹的末端正好位于剖面中低阻抗异常的顶部，这与我们的特征化结果一致。图16显示了沿Well A12轨迹提取的SCI阻抗与测井曲线之间的比较。反演结果中的低阻抗区间与测量孔隙度曲线的高值区间非常吻合。图15显示了通过时深转换得到的inline1714的SCI阻抗剖面。蓝线表示Well A12的井轨迹在剖面上的投影，其中由于循环损失事件，钻井在轨迹末端终止。沿井轨迹的曲线是测量孔隙度对数。图16显示了沿Well A12轨迹的SCI阻抗曲线与测井曲线的比较。MD表示测量深度，LN_RT表示真实电阻率的自然对数，POR表示测量孔隙度，SCI阻抗表示沿井轨迹提取的反演结果。

5. 讨论
传统的基于像素的反演方法存在结构不一致性和非唯一性问题。相比之下，形状约束的地震反演能够很好地保留地质界面处的参数突变，从而突出异常地质体。这一特性使得形状约束反演方法特别适合于空间隔离的腔体储层的特征化。在塔里木盆地的现场数据应用中，反演结果清晰地揭示了喀斯特储层的空间分布模式，这与区域地质情况和井数据高度吻合。与传统阻抗反演方法相比，这种方法有效降低了储层解释的不确定性。作为形状约束反演的输出产品之一，界面集本身可以在合理的权重参数设置下作为地震属性用于储层特征化，无需进行阈值处理来划分三维储层体。对于实际的地质建模，反演得到的三维界面集可以作为相控约束用于建模。此外，从地震反演中得到的分层结构界面可以重新用作精细的地层边界，以支持地层划分和后续的储层建模。我们完成了前述工作区域的完整三维反演，图17显示了通过形状约束反演得到的目标层界面集的平面视图。可以看出，喀斯特储层界面在工作区域的北部广泛分布，而在南部则主要沿走滑断层呈带状分布。工作区域的北部和南部分别属于埋藏丘陵喀斯特区和断层控制喀斯特区。这是因为在覆盖层的尖灭边界（绿色曲线）以北，目标层在埋藏丘陵喀斯特区形成了大规模的洞穴[40,41]。这与该地区的先验地质知识相符。图17显示了整个调查范围内目标层界面集的切片。其中，红线代表目标层的走滑断层系统，绿色曲线表示覆盖层的尖灭线。尽管所提出的方法已经展示了其优势和适用性，但仍有一些实际方面和方法论上的局限性需要进一步讨论。三维地震数据的多道反演以及界面集的迭代更新和存储引入了额外的计算和内存需求。在现场数据计算中，整个数据集包含950个测线，每个测线包含1180个道，垂直反演维度包括400个采样点。计算是在配备有64 GB内存和32个CPU核心的服务器上进行的，总计算时间为9小时。值得注意的是，如果计算仅限于单道基础，则无法应用形状约束反演框架。在我们的实际应用中，我们采用了分块处理策略来缓解这个问题，并将在未来的工作中继续研究提高计算效率的方法。其次，形状约束正则化项包含了三个权重参数，每个参数控制不同模型特征的惩罚权重。整体行为由它们的相对比率而不是绝对值决定。在未来的工作中，我们将分别研究大规模和小规模反演的自适应参数选择。第三，在现有研究的基础上，我们计划将多源地质先验信息（包括沉积相、断层系统和喀斯特古地貌）整合到反演框架中，以进一步提高反演结果的地质合理性。此外，我们将开发一个深度域反演框架，以减少与时深转换相关的误差。为了解决深度变量小波带来的挑战，我们计划引入一个深度域广义S变换或点扩散函数，以估计每个深度采样点的波数信息并构建一个空间变化的小波[42]。这将实现深度域反演，使反演结果能够更可靠地整合到三维地质建模工作流程中。

6. 结论
在本文中，我们提出了一种基于形状约束地震反演的腔体储层特征化方法，并进行了相应的数值实验。在此基础上，我们将所提出的方法应用于真实的地震场数据。本研究的核心创新在于将MS函数引入到多道堆叠地震阻抗反演中。基于MS函数框架实现了隐式的形状表示。它在真实界面处只形成狭窄的过渡层，避免了在传统反演中为了减少数据不匹配而牺牲异常体的结构清晰度。数值结果证实，与现有技术相比，形状约束反演能够更好地恢复地质结构。值得注意的是，即使在先验信息有限的情况下，它也能显著提高阻抗值的准确性和横向分辨率。这是因为反演框架隐式地促进了分段常数区域，大大减少了优化过程中的自由度。实际现场的案例研究表明，与传统方法相比，所提出的方法在表征腔体型储层方面有显著改进。通过将腔体隐式地识别为分段常数区域而不是孤立的网格单元，我们大大减少了非物理自由度，保留了真实腔体边界处的清晰阻抗对比，并避免了过度平滑低阻抗储层异常。这使得该方法特别适合于表征具有清晰边界的离散碳酸盐岩腔体储层。通过与钻井和测井数据的相关性分析，我们确认了井中解释的腔体型储层与我们的结果一致，从而验证了腔体特征化的可靠性和真实性。这些验证结果表明，所提出的方法能够有效地表征地下异常体的结构形态，为在高度异质性的储层（如碳酸盐岩）中的油气勘探提供了宝贵的支持。