基于振幅谱互相关的薄层厚度计算

时间：2024-05-22

刘道理，刘军，杨登锋，魏旭旺，刘琼

(中海石油(中国)有限公司深圳分公司研究院，深圳 518000)

0 引言

随着油气勘探的不断深入，岩性圈闭逐渐成为勘探开发的重要领域之一。相对于常规构造圈闭，岩性圈闭的隐蔽性强、储层薄。薄层砂体厚度的有效识别和砂体尖灭的准确刻画是岩性圈闭勘探的重点和难点[1-3]。

为了能够准确可靠地识别薄层，很多专家在该领域做了大量的研究,Widess[4]首次利用楔状体模型研究薄层反射与地层厚度的关系并给出了薄层的定义，将厚度小于入射子波八分之一主波长的地层定义为薄层;Neidell等[5]定义四分之一主波长对应的地层厚度为调谐厚度,在实际应用中受子波和噪声影响，业界将厚度小于调谐厚度的地层定义为薄层;Kallweit[6-9]进一步研究了地层厚度与振幅之间的关系，得出了当地层厚度等于四分之一主波长时，反射波振幅达到最大，随着地层厚度变薄，反射振幅也相应减小的结论，从而为利用地震波振幅信息识别薄层奠定了理论基础；凌云等[10]利用振幅的调谐作用探测地层厚度小于1/4波长地质目标；李伟等[11]针对薄层砂体与厚层砂体不同的地震振幅响应特征，采用分频段地震资料优选地震属性，对秦皇岛32-6油田砂体进行了预测。另外，不少的研究同样发现地震波的频谱特征与薄层厚度密切相关，即当地层厚度小于四分之一波长时，地层变薄，其能量变小，频率有所提高[12-18]。根据这一特性，可以利用峰值频率来计算薄层的厚度。Liu等[19]利用线性拟合方法近似预测薄层厚度与峰值频率的定量关系；Sun等[20]在泰勒展开的基础上推导出利用峰值频率求解薄层厚度的解析式；Partyka等[21]提出谱分解技术，利用局部地震频谱信息反演得到薄层反射系数，该项技术的缺点是在有限地震频带内难以准确拾取频陷周期；Portniaguine等[22]提出了谱反演的方法，可直接反演出反射系数序列，在薄层预测中成功应用；Puryear等[23]对传统谱反演方法进一步改进，将反射系数分解成奇偶分量，利用偶分量在厚度趋于零时的有效干涉提高了薄层反演的稳定性和分辨率；王贻朋等[24]基于匹配追踪算法保证了高精度时频分辨率的同时，建立地震信号偶分量的峰值频率-薄层厚度关系曲线模板，实现谱分解薄层预测；张军华等[2]将压缩感知技术引入谱反演中反演反射系数，结果证实分辨率明显提升。基于稀疏脉冲约束的常规反演技术在薄层预测中有着广泛应用，特别是利用自然伽马等其他测井曲线对声波测井曲线进行曲线重构的拟声波反演在薄层识别上效果比较理想；在井网密集区，地质统计学反演可用于薄层预测[2,25-26]。在无井或少井区这些常规的地震反演技术对薄砂体的刻画分辨率有限。此外，基于多属性模式识别技术同样广泛应用于薄砂体检测中，模式识别方法有聚类分析、人工神经网络、支持向量机、主成分分析等[27-29]。其基本原理是优选提取多种敏感地震属性，根据样本学习，在井约束的条件下可以实现薄层厚度预测，合理可靠的地震属性优选和足够多的学习样本是该类方法应用的前提。

Yang Senlin等[30]提出了利用振幅谱互相关预测地震品质因子，基于该思路，笔者提出了一种利用振幅谱互相关函数计算薄层厚度的新方法。新方法首先构建目的层的振幅谱与子波振幅谱乘以反射系数振幅谱之间的相关系数函数。然后，利用高斯迭代方法使得相关系数函数达到最大值，此时对应的厚度就是薄层的厚度。通过模型测试和实际资料的应用证实新方法是可行、有效的。相较于峰值频率方法，新方法能够更精确地计算薄层厚度。

1 基本理论

当地震波在地层中传播时，地层反射系数可以看成一个滤波器，地震子波为该滤波器的输入。对于单层厚度很薄的泥包砂储层，地震子波在其中传播后产生的地震响应振幅谱A(f)可以表示为式(1)[12,31]。

(1)

其中：f为频率；W(f)为子波振幅谱；r1和r2分别是薄层顶底反射系数；Δτ0为子波在该薄层中传播的双程旅行时。方便起见，定义薄层反射响应滤波器如式(2)所示。

(2)

其中：Δτ为子波在任意薄层中传播过程中所需双程旅行时。将薄层响应滤波器与震源子波相乘，得到式(3)。

(3)

薄层振幅谱和震源子波薄层响应振幅谱之间的互相关系数函数为式(4)。

(4)

为了求取相关系数C(Δτ)的最大值，需要求解如下方程：

(5)

图1 振幅谱互相关法计算薄层厚度流程图Fig.1 The flow chart of thin bed thickness calculation from spectral correlation

图2 合成地震记录及其振幅谱Fig.2 The synthetic data and its amplitude spectrum(a)合成地震记录；(b)其振幅谱

式(5)是一个非线性方程，其解可以通过牛顿迭代法求取，具体迭代公式为式(6)。

(6)

其中：n表示迭代步数。初始迭代变量Δτ1的选取原则是：双程旅行时大于“0”小于λ/4。迭代收敛条件为

|Δτn+1-Δτn|<ε

(7)

其中：ε为任意小的常数。

2 模型试算

图3 振幅谱互相关函数系数曲线Fig.3 The amplitude spectral correlation function coefficient curves

图4 楔状体模型Fig.4 The wedge model

为了检验振幅谱互相关法估计薄层厚度的可行性和有效性，首先选择三个不同厚度水平层状薄层单砂体模型进行试算，三套单砂体对应厚度分别为5 m、10 m、20 m。砂体对应速度为3 500 m/s，密度为2 500 kg/m3；上下泥岩速度为3 000 m/s，密度为2 250 kg/m3。选用主频40 Hz雷克子波，与模型褶积合成地震记录，如图2(a)所示；其中1、2、3道分别对应的砂体厚度为5 m、10 m和20 m。图2(b)为每一道合成记录所对应的振幅谱。基于该振幅谱且根据式(4)构建每一道所对应的互相关系数函数，其所对应函数曲线如图3所示；其中蓝、绿、红曲线分别对应5 m、10 m、20 m单砂体厚度。由图3可以看出，对于某一特定厚度的薄层，其互相关系数函数都有一个极大值，其所对应的双程旅行时Δτ就是为所要求取的薄层厚度所对应的双程旅行时。例如，图3中红色互相关函数曲线，在图3上双程旅行时约为11.4 ms处互相关函数出现极大值，当砂体速度为3 500 m/s，该极大值对应的地层厚度约为20 m。这里对比分析振幅谱互相关法和峰值频率法两种方法薄层估计的精度。利用振幅谱互相关法估计的三套薄层砂岩的厚度分别为4.99 m、9.99 m、20 m；对于峰值频率法估计的厚度分别为5.03 m、10.63 m、24.30 m。从模型试算结果可见，振幅谱互相关法估计薄层厚度是可行有效的。相较于峰值频率法，新方法估计薄层厚度更精确，这是由于新方法不需要泰勒展开近似。

图5 基于楔状体模型的合成地震记录Fig.5 The synthetic data based on the wedge model

图6 合成地震记录所对应的振幅谱Fig.6 The amplitude spectrum of the synthetic data

图7 利用振幅谱互相关法和峰值频率法估计的薄层厚度Fig.7 The estimated thin bed thickness using the spectral correlation method and the peak frequency method

为了进一步验证新方法的可行性，设计一个楔状体模型，如图4所示。同样选用主频40 Hz雷克子波，与该模型褶积合成地震记录,如图5所示。图5中每一地震道对应振幅谱如图6所示。基于振幅谱互相关法和峰值频率法计算薄层厚度如图7所示。从图7中可以发现，新方法计算的砂体厚度曲线与实际厚度曲线叠合的非常好；砂体厚度越大，峰值频率法计算的砂体厚度偏离实际值越大，计算的误差越大。相对于峰值频率法，新方法计算精度明显改善。该模型测试进一步证实新方法识别薄层是可行有效的。

图8 研究区已钻井岩性柱状图频率法估计的薄层厚度Fig.8 The lithology histograms of drilled wells in research area(a)H1井；(b) H2井；(c) H3井

图9 过研究区某一地震剖面频率法估计的薄层厚度Fig.9 A seismic section in research area

图10 S30层附近地震振幅谱频率法估计的薄层厚度Fig.10 The amplitude spectrum nearby S30

图11 提取的地震子波及其振幅谱Fig.11 The abstracted wavelet and its amplitude spectrum(a)提取的地震子波;(b)振幅谱

3 应用实例和效果

选取珠江口盆地惠州凹陷南部某区块进行应用分析。研究区内，在地震剖面上无法直接分辨，厚度小于四分之一波长的地层都称之为薄层。研究区珠江组S系列沉积体系以三角洲前缘近端或远端沉积为主，砂体沉积厚度薄，岩性圈闭发育。S系列沉积的薄层砂体厚度可靠估计对于该区域油气勘探和开发至关重要。例如S30层是该地区的主要目的层，已钻井结果表明目的层岩性为粉砂岩，其上下地层岩性都为泥岩；地层组合特征为厚层泥岩夹单层薄砂岩，属于典型的泥包砂类型，如图8所示。其中，H1井S30层砂体厚度为20 m，H3井S30层砂体厚度为7.2 m。图9位该地区原始地震资料某一任意地震剖面，其受噪声干扰弱、信噪比高，同相轴连续性好；在S30层附近振幅谱如图10所示，主频约为40 Hz，有效频带宽度为10 Hz ～75 Hz，通过已钻井分析目的层砂岩速度约为3 500 m/s，目的层砂体沉积厚度在四分之一波长以内，在地震剖面上难以直接分辨。S30层这套砂体沉积在地震剖面上表现为波谷反射，见图9中蓝色虚线；在剖面上很难直接识别其尖灭点及计算其沉积厚度。这里用振幅谱互相关法预测薄层厚度，新方法可靠识别薄层厚度的难点主要有：①砂体厚度越薄，预测难度越大，从图6中可以看出，地层厚度越薄，振幅谱的差异越小，增加了预测难度；②需构建可靠的地震子波；③需要估算出可靠的地层速度。速度精度会影响厚度估算的精度，在无井区需通过对地震道集进行速度分析得到目的层段速度；在已钻井区，通过声波测井资料获取目的层段速度。新方法具体步骤如下：

1)利用商业软件在目的层段提取统计子波，提取的子波波形和振幅谱特征，如图11所示。

2)假设薄层砂岩顶底反射强度相等，极性相反，得到泥包砂类型沉积的薄层反射响应滤波器。

3)将提取子波的振幅谱和薄层反射响应滤波器相乘，得到不同地层厚度下提取子波的响应振幅。

4)在三维地震数据体上对S30反射地震同相轴进行层位解释追踪，再沿着所解释的层位上下15 ms开时窗，对S30层地震反射数据做短时傅里叶变换，得到S30薄层反射振幅谱。

5)利用实际S30薄层反射振幅谱和提取子波响应振幅谱构建得到关于地层厚度Δτ的相关系数函数。

6)求取互相关函数最大值对应的双程旅行时。

7)由双程旅行时和S30薄层砂体速度，计算出薄层砂体厚度。

图12 利用振幅谱互相关法计算的S30砂体厚度与其等值线叠合图Fig.12 The overlayed plots for the estimated thickness of S30 using the spectral correlation method and the thickness contour

图13 S30砂体厚度等值线图Fig.13 The thickness contour of S30

图14 沿S30层提取的均方根振幅平面分布图Fig.14 The RMS amplitude calculated around S30

图15 沿S30层提取的峰值频率平面分布图Fig.15 The peak frequency calculated around S30

由图12、图13可以看出，结果显示沉积砂体厚度普遍小于25 m，图中暖色调区(偏红色)指示沉积砂体厚度较厚，冷色调区(偏蓝色)指示沉积砂体厚度较薄。在图12砂体厚度平面图的东边，可以看出发育一套北北东向条带状较厚砂体沉积。为了进一步验证所估计砂体厚度的可靠性，同样沿着S30层提取了均方根振幅和峰值频率属性，分别如图14和图15所示。对比图13、图14和图15可以看出，它们之间具有很好的对应性。砂体沉积厚的区域对应于强振幅、低峰值频率；砂体沉积薄的区域对应于弱振幅、高峰值频率。在图14和图15中同样存在一套北北东向条带状沉积特征。

结合振幅、频率属性图和厚度图综合分析S30砂体沉积微相，该套砂体地层主要有两大沉积体系，西侧靠近物源方向为三角洲前缘近端或远端河口坝-席状砂沉积体系；在三角洲前缘外的东南侧由一系列北东向分布的条带状陆架砂脊沉积体系。三角洲前缘沉积不断向海盆中心推进过程中，受波浪作用将三角洲前缘的粗粒沉积物搬运到远离前三角洲区沉积，在沉积物供给充足的条件下，波浪作用会继续搬运同期沉积的粗粒沉积物，从而在三角洲前缘外侧形成一系列平行于岸线的砂脊，另外受后期北东-南西向潮流的冲刷和改造，导致这些砂脊条带状分明[32]。从沉积微相上综合分析认为砂体沉积厚度分布较为合理可靠。此区域三口已钻井H1、H2和H3井，应用新方法计算得到的三口井处S30层砂体厚度分别为H1井20.3 m、H2井16.5 m、H3井7.6 m；实钻结果三口井处S30层砂体厚度分别为H1井20 m、H2井16 m、H3井7.2 m。对比计算砂体厚度和实钻砂体厚度，可以发现计算结果与实钻厚度非常接近，误差很小。通过实际地震资料应用，证实了振幅谱互相关法计算薄层砂体厚度是可行、有效的。