0 引言

砂砾岩体是一种广泛存在于中国东部的碎屑岩储集体，主要发育于陆相断陷湖盆地中的陡坡带^[1]。由于其侧向上紧邻深洼烃源岩，油气成藏条件优越，可发育岩性油气藏、构造―岩性油气藏和构造油气藏等^[2-3]，勘探潜力巨大^[4]，一直以来都是常规油气藏勘探的重要领域之一^[5-6]，同时也是非常规油气勘探的重要研究方向^[7]。但砂砾岩体主要由多期砂砾沉积物滑塌、搬运、快速堆积和叠置而成^[8]，其沉积厚度大，埋深跨度大^[9]，横向变化快，沉积类型多样，空间展布及叠置关系复杂，具有强非均质性特征。与此对应，砂砾岩体成像结果的分辨率和信噪比偏低，顶、底界面之间的地震反射能量较弱，对内部结构特征认识不清^[10]，所以地震资料解释和储层预测难度较大，严重制约了砂砾岩体油气藏的勘探与开发。

除上述问题以外，砂砾岩体的岩性致密，粒度较粗，整体速度普遍较高，且明显不同于围岩的速度，而在地震数据处理过程中容易忽视砂砾岩体速度的空间变化，尤其是时间域砂砾岩体成像结果普遍存在构造假象^[11]，油水关系不清，圈闭难以落实。前人分析了时间域假象的成像机理，并对比了时、深域的成像差异，认为与时间域成像相比，深度域成像能够更好地刻画砂砾岩体产状^[12]。总体而言，为提高砂砾岩体的速度模型精度，需要将处理技术与解释信息相结合，开展有针对性的处理工作^[13]。

位于济阳凹陷的BS工区构造特征复杂，其中的陡坡带发育砂砾岩体，岩体与围岩速度差异较大，常规网格层析反演的结果很难反应出这种速度差异。本文在常规网格层析反演的基础之上，引入解释所得的砂砾岩体信息，通过施加约束条件控制岩体速度和空间分布，提高反演速度模型的精度，进而提高该工区砂砾岩体的成像质量。

1 成像域网格层析反演

如何精确计算地震波速度，一直是地震资料处理和解释中的难题。当前生产应用的主流深度域建模方法是成像域网格层析反演方法^[14-16]。该方法主要由构造倾角扫描、剩余曲率分析、反向射线追踪、构建方程组并求解更新量等步骤构成。对应的目标函数为

$ O=\frac{1}{2}\left({‖{\Delta }\boldsymbol{T}-\boldsymbol{G}{\Delta }\boldsymbol{m}‖}_{2}^{2}+{\varepsilon }_{\mathrm{s}}{‖\boldsymbol{L}{\Delta }\boldsymbol{m}‖}_{2}^{2}\right) $

(1)

式中：ΔT表示由成像道集拾取的剩余深度差转化得到的走时残差；G为由射线路径构建的敏感度矩阵；Δm表示待求解的慢度更新量；ε_s为平滑项权重系数；L表示由拉普拉斯算子构成的平滑矩阵，实际中可结合构造倾角使用^[17]。式(1)对应的反演方程组为

$ \left[\begin{array}{c}\boldsymbol{G}\\ \sqrt{{\varepsilon }_{\mathrm{s}}}\boldsymbol{L}\end{array}\right]{\Delta }\boldsymbol{m}=\left[\begin{array}{c}{\Delta }\boldsymbol{T}\\ 0\end{array}\right] $

(2)

该方程组可应用最小平方正交分解(LSQR)算法求解^[18]。配合叠前深度偏移使用，通过迭代更新可拉平成像道集，聚焦反射波能量，提高成像质量。

2 砂砾岩体特征约束层析反演

式(1)对应方法所得的速度模型较为平滑，分辨率低。引入倾角约束后，分辨率有所提高，结果与构造趋势更一致，但仍无法适用于岩体或地层速度突变的情况。对岩体速度建模的实现方法以融合建模^[19-20]和层位约束层析反演^[21-24]为主。前者直接向模型中引入岩体空间分布和速度信息，但难以拉平成像道集；后者可更新速度，但只引入岩体空间分布信息，没有对速度进行限制。

借鉴层位约束的思想，本文进一步在速度变化范围和空间分布区域两方面添加约束条件，以控制反演过程，获得符合砂砾岩体特征的速度模型。

在速度变化范围方面，本文在层析反演过程中应用不等式约束，将砂砾岩体的成像速度限制在测井速度的预设区间之内，以提高速度层析建模的精度和稳定性。罚函数法^[25]是目前实施不等式约束的常用方法之一，通过向目标函数中引入罚函数项，可将原有问题转化为无约束问题。本文采用内罚函数法，通过向目标函数中引入对数型障碍函数，实现不等式约束层析。在实施反演之前，需要通过融合建模方法将岩体速度赋值为测井速度，以保证速度在预设区间内。

在空间分布区域方面，层析反演的重点控制在砂砾岩体所在区域之内。对于岩体以外的区域，应减小反演对围岩速度的改变。为实现上述目的，需要对当前的成像剖面进行解释，提取砂砾岩体空间分布区域并引入到反演之中，对区域内施加不等式约束，对区域外施加正则化约束，使外部的更新量趋于0。

基于以上分析，将式(1)的目标函数拓展为

$ \begin{align}O=&\frac{1}{2}\left[{‖{\Delta }\boldsymbol{T}-\boldsymbol{G}{\Delta }\boldsymbol{m}‖}_{2}^{2}+{\varepsilon }_{\mathrm{s}}{‖\boldsymbol{L}{\Delta }\boldsymbol{m}‖}_{2}^{2}+\right.\\ &\left.{\varepsilon }_{\mathrm{c}}{‖\left(\boldsymbol{I}-\boldsymbol{C}\right){\Delta }\boldsymbol{m}‖}_{2}^{2}\right]-{\varepsilon }_{\mathrm{p}}\sum _{i\in \boldsymbol{\varOmega }}\left[\mathrm{l}\mathrm{n}\left({m}_{i}+{\Delta }{m}_{i}-{m}_{i, \mathrm{m}\mathrm{i}\mathrm{n}}\right)+\right.\\ &\left.\mathrm{l}\mathrm{n}\left({m}_{i, \mathrm{m}\mathrm{a}\mathrm{x}}-{m}_{i}-{\Delta }{m}_{i}\right)\right] \end{align} $

(3)

式中：I为单位矩阵；ε_c为正则化系数；Ω表示砂砾岩体的空间分布区域；C为与Ω区域相对应的对角矩阵，由0或1构成，在区域之内的元素为1，之外的元素为0，该项约束岩体外的区域；对数函数项为施加的罚函数，ε_p为其权重系数；m_i为模型当前值，m_{i, min}、m_{i, max}分别为设置的下限和上限。对式(3)求导，并对结果中的罚函数部分作一阶近似，可得相应的最小二乘解

$ \begin{align}{\Delta }\boldsymbol{m}=&\left[{\boldsymbol{G}}^{\mathrm{T}}\boldsymbol{G}+{\varepsilon }_{\mathrm{s}}{\boldsymbol{L}}^{\mathrm{T}}\boldsymbol{L}+{\varepsilon }_{\mathrm{c}}{\left(\boldsymbol{I}-\boldsymbol{C}\right)}^{\mathrm{T}}\left(\boldsymbol{I}-\boldsymbol{C}\right)\right.+\\& {\begin{array}{cc}& \left.{\varepsilon }_{\mathrm{p}}\boldsymbol{C}\left({\boldsymbol{R}}_{\mathrm{m}\mathrm{i}\mathrm{n}}^{\mathrm{T}}{\boldsymbol{R}}_{\mathrm{m}\mathrm{i}\mathrm{n}}+{\boldsymbol{R}}_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{T}}{\boldsymbol{R}}_{\mathrm{m}\mathrm{a}\mathrm{x}}\right)\right]\end{array}}^{-1}\times \\& \begin{array}{cc}& \left[{\boldsymbol{G}}^{\mathrm{T}}{\Delta }\boldsymbol{T}+{\varepsilon }_{\mathrm{p}}\boldsymbol{C}\left({\boldsymbol{R}}_{\mathrm{m}\mathrm{i}\mathrm{n}}+{\boldsymbol{R}}_{\mathrm{m}\mathrm{a}\mathrm{x}}\right)\right]\end{array}\end{align} $

(4)

式中R_min、R_max为与所设的下限、上限有关的对角矩阵，即

$ \left\{\begin{array}{l}{\boldsymbol{R}}_{\mathrm{m}\mathrm{i}\mathrm{n}}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(\frac{1}{{m}_{i}-{m}_{i, \mathrm{m}\mathrm{i}\mathrm{n}}}\right)\\ {\boldsymbol{R}}_{\mathrm{m}\mathrm{a}\mathrm{x}}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(\frac{1}{{m}_{i}-{m}_{i, \mathrm{m}\mathrm{a}\mathrm{x}}}\right)\end{array}\right. $

(5)

在应用LSQR算法求解时，要输入的反演方程组为

$ \left(\begin{array}{c}\boldsymbol{G}\\ \sqrt{{\varepsilon }_{\mathrm{s}}}\boldsymbol{L}\\ \sqrt{{\varepsilon }_{\mathrm{c}}}\left(\boldsymbol{I}-\boldsymbol{C}\right)\\ \sqrt{{\varepsilon }_{\mathrm{p}}}\boldsymbol{C}{\boldsymbol{R}}_{\mathrm{m}\mathrm{i}\mathrm{n}}\\ \sqrt{{\varepsilon }_{\mathrm{p}}}\boldsymbol{C}{\boldsymbol{R}}_{\mathrm{m}\mathrm{a}\mathrm{x}}\end{array}\right){\Delta }\boldsymbol{m}=\left(\begin{array}{c}{\Delta }\boldsymbol{T}\\ 0\\ 0\\ \sqrt{{\varepsilon }_{\mathrm{p}}}\boldsymbol{C}\\ \sqrt{{\varepsilon }_{\mathrm{p}}}\boldsymbol{C}\end{array}\right) $

(6)

由式(6)可知，当模型当前值逐步趋于区间边界时，R矩阵中对应元素的值将急剧增大，从而减小模型更新的步长，并会产生反向改变的效果；当模型当前值远离区间边界时，R矩阵中对应元素的值将减小，结合较小的权重系数，则不会对问题的解产生显著影响。

3 实际应用

BS工区位于东营凹陷北部，含油层系多，油气藏类型丰富，有较大的勘探潜力。该地区构造特征较为复杂，浅部地层平缓且发育多条断层，深部为太古界潜山，地震速度高，内幕不清晰。潜山侧壁的陡坡带有砂砾岩体发育，其分布形态不规则，沉积较为杂乱。该工区地震数据经过常规处理，有声波测井、地震成果剖面、速度模型等资料，其中速度模型由常规网格层析反演获得，道集中岩体上覆地层的同相轴已基本拉平。通过分析测井数据得知砂砾岩体速度较为稳定，整体在4.6±0.4 km/s内变化，而其上覆地层的速度小于3.7 km/s，下伏潜山的速度超过5.5 km/s，均与砂砾岩体有较大速度差异。

为提取到较准确的岩体空间分布信息，本文基于原有成果资料开展测试应用工作，并采用原有成果资料作参照，以便对比分析。具体采用的技术流程如图 1所示。首先对地震资料进行偏移前的预处理工作，并基于原有成果剖面识别刻画砂砾岩体空间分布，通过填充、融合方式建立初始速度模型；然后进行叠前深度偏移，生成新的成像剖面和道集；之后基于成像道集拾取剩余曲率，基于成像剖面开展倾角分析，并识别、刻画出剖面中的砂砾岩体边界，再进行反向射线追踪，构建方程组；基于测井分析结果设置速度区间，基于识别、刻画结果设置约束区域，向原方程组中添加平滑约束、正则化约束、不等式约束等条件方程组，最后统一求解并更新速度模型，再进行下一次深度偏移，直到道集拉平、成像结果改善为止，输出最终结果。本文此次应用中，速度区间设为[4.2 km/s, 5 km/s]，整个过程进行1次反演计算。

图 2为解释人员从成像剖面中刻画的砂砾岩体三维空间分布，岩体深度范围为[1.2 km, 5.6 km]，x方向跨度超过20 km，y方向跨度接近10 km，最大厚度近2 km。图 3为原有速度模型、融合速度模型、反演速度模型的对比。由图可见，三个模型整体上没有明显差别，差别主要存在于砂砾岩体区域(图中箭头所指位置)。与原有速度模型(图 3a)相比，反演速度模型(图 3c)中的砂砾岩体速度与前者有明显不同，速度突变更为清晰；与融合速度模型(图 3b)相比，反演速度模型能够一定程度上反映出砂砾岩体内部的速度变化。通过对比分析可知，反演中引入的砂砾岩体信息起到了局部约束的效果。

图 4为不同速度模型和成像剖面的叠合显示。对比可知，原有速度模型(图 4a)与剖面中的砂砾岩体(蓝色虚线所示位置)分布并不吻合，而反演速度模型(图 4b)的吻合程度更高，且其中的岩体速度、上覆地层速度、潜山速度与测井分析结果基本一致。上述对比结果说明与常规网格层析的结果相比，本文结果更符合地质认识，证明了约束反演的有效性。

为进一步验证反演的效果，图 5对比了融合速度模型与反演速度模型的偏移结果。与融合速度模型的成像剖面(图 5a)相比，反演速度模型的成像剖面(图 5b)整体质量相当，但对比箭头所指之处，反演速度模型的成像剖面中的砂砾岩体底界面反射更连续。以上结果证明，与融合建模技术相比，本文结果在一定程度上提高了砂砾岩体的成像精度。

4 结束语

本文应用约束层析反演技术解决砂砾岩体速度建模问题，实际应用的结果表明：

与原有速度模型相比，结合砂砾岩体特征的约束层析反演结果与成像剖面的构造吻合程度更高，砂砾岩体速度与测井资料更符合，反演计算并未对其他区域的速度产生显著影响。

与融合速度模型相比，反演速度模型能在一定程度上反映出砂砾岩体内部的速度变化，二者的成像剖面质量整体相当，但后者的砂砾岩体底界成像连续性有所提高。

约束层析反演技术可以有效提高砂砾岩体速度模型的精度，本文处理流程有望推广应用到其他地质体的地震速度建模。