CN105093294A - 基于可变模态分解的地震波衰减梯度估计方法 - Google Patents

Abstract

本发明公开了一种基于可变模态分解的地震波衰减梯度估计方法，该方法包括：利用测井、录井和合成地震记录标定目标层；对地震数据逐道进行可变模态分解，结合最小二乘法提取模态函数的衰减梯度，并对所述模态函数的衰减梯度进行加权求和获得各地震道的衰减梯度值；利用地震衰减梯度特征体确定目标层的岩性及烃类性质。本发明能够更有效地适用于低信噪比及频谱波动较大的地震信号的衰减梯度估计，因此，能够提高吸收衰减梯度指示油气的有效性，准确用于识别岩性、烃类性质等。

Description

基于可变模态分解的地震波衰减梯度估计方法

技术领域

本发明涉及石油地球物理勘探技术领域，尤其涉及一种基于可变模态分解的地震波衰减梯度估计方法。

背景技术

作为介质内在的属性，地震衰减是地震波在地下介质传播过程中发生的能量损失。除了在岩石中地震波传播期间的波前发散，岩石的非完全弹性特性也可能导致地震能量的衰减。由于地震波的弹性能量被不可逆地转换成热能然后消散，地震波的总能量从而发生衰减。影响地震波衰减的因素很多，但在一个相对稳定的地层构造，垂直和水平方向岩性变化不大的情况下，地震波的衰减主要是由流体特性导致的。在含油气层的内部及其下部，地震波的能量会发生明显的高频衰减。地层吸收的不一致性表现在地震信号的频谱上，体现出了低频强能量，高频弱能量的特征。利用地震信号高频衰减异常从地震反射数据进行地质解释及油气指示的衰减梯度估计技术，根据地震波能量及对应的频率拟合出能量与频率的衰减梯度，求出振幅衰减梯度因子，从而进行吸收衰减分析。常规的能量吸收分析方法采用基于短时傅里叶变换、小波变换等的时频分析方法及两点斜率或线性拟合的方法。由于短时傅里叶变换等常规时频分析方法受测不准原理的限制，时频分辨率较低，且由于各自方法的一些缺点导致方法的适用性受到了很大的限制，如，短时傅里叶变换很难保证在时、频域都有足够高的分辨率；小波分析的有效性取决于小波函数的选择；S变换的小波基函数是固定的，不能满足实际数据处理的需求。此外，利用两点斜率或线性拟合的方法求取能量与频率的衰减梯度，仅适用于高信噪比的地震信号或者具有光滑频谱的地震波，对于低信噪比的地震信号及频谱波动较大的地震波效果变差，尤其是当高频段噪声能量较强时，吸收衰减梯度不能完全指示油气。

发明内容

为解决上述技术问题，本发明的目的是提供一种基于可变模态分解的地震波衰减梯度估计方法，通过方法能够更有效地适用于低信噪比及频谱波动较大的地震信号的衰减梯度估计，因此，能够提高吸收衰减梯度指示油气的有效性，准确用于识别岩性、烃类性质等。

本发明的目的通过以下的技术方案来实现：

基于可变模态分解的地震波衰减梯度估计方法，包括：

A利用测井、录井和合成地震记录标定目标层；

B对地震数据逐道进行可变模态分解，结合最小二乘法提取模态函数的衰减梯度，并对所述模态函数的衰减梯度进行加权求和获得各地震道的衰减梯度值；

C利用地震衰减梯度特征体确定目标层的岩性及烃类性质。

与现有技术相比，本发明的一个或多个实施例可以具有如下优点：

1、使用了较经验模态分解等方法更具抗噪声性的可变模态分解方法进行地震道的分解。各个模态函数是具有不同频带宽度的窄带信号，其频谱相对原始地震信号更为平滑，有利于衰减梯度的准确计算，能够反映特定频带处的衰减情况。

2、对地震道进行可变模态分解产生的不同模态函数分别结合希尔伯特变换通过可变长度频率域窗函数利用最小二乘法计算衰减梯度，再通过相关加权系数对各个模态计算的衰减梯度进行加权求和，加强了地震信号中的主要贡献分量，能够反映出地震信号的微弱变化，加强了油气地震响应特征，加强了储层与非储层的区别，有利于识别岩性及烃类性质。

3、地震衰减梯度特征体的计算可适用于二维或三维数据，计算方式灵活多样，可以根据实际需求计算时间切片，地层切片或沿层切片进行分析。

4、本发明较常规衰减梯度估计方法具有更高的分辨率，也适用于低信噪比或者频谱波动较大的地震信号。

附图说明

图1是基于可变模态分解的地震波衰减梯度估计方法流程图；

图2是本发明实施例提供的一条二维叠后偏移过井剖面(目标层段)图；

图3是过井道及其VMD分解结果图；

图4是利用本发明计算的时频谱和短时傅里叶变换、小波变换计算的时频谱的对比结果图；

图5是利用本发明计算的过井道的地震衰减梯度归一化曲线图；

图6是利用小波变换结合两点斜率的常规衰减梯度估计方法计算的过井道的地震衰减梯度归一化曲线图；

图7是利用小波变换结合最小二乘拟合法计算的过井道的地震衰减梯度归一化曲线图；

图8是利用本发明计算的归一化的地震衰减梯度特征剖面图；

图9是利用小波变换结合两点斜率的常规衰减梯度估计方法计算的归一化的地震衰减梯度特征剖面图；

图10是利用小波变换结合两点斜率的常规衰减梯度估计方法计算的归一化的地震衰减梯度特征剖面图；

图11是某一盆地某气田的一条二维叠后偏移过含气井剖面图；

图12是利用本发明技术计算的归一化的地震衰减梯度特征剖面图；

图13是利用小波变换结合最小二乘拟合法计算的归一化的地震衰减梯度特征剖面图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面将结合实施例及附图对本发明作进一步详细的描述。

如图1所示，为基于可变模态分解的地震波衰减梯度估计方法流程，所述方法包括如下步骤：

步骤10利用测井、录井和合成地震记录准确标定目标层；

步骤20对地震数据逐道进行可变模态分解，结合最小二乘法提取模态函数的衰减梯度，并对所述模态函数的衰减梯度进行相关加权求和获得各地震道的衰减梯度值；

利用基于可变模态分解的地震波衰减梯度估计方法沿着目标层开时窗提取地震衰减梯度。地震衰减梯度计算的核心问题是对地震记录逐道进行可变模态分解，对每一道地震记录生成的各个模态函数，结合希尔伯特变换生成相应的时频图，计算每个时间采样点处对应的对数频谱图，结合最小二乘拟合法计算各个模态函数在每一个时间采样点上的衰减梯度，利用相关加权法将各道模拟函数的衰减梯度相加，生成该地震记录的衰减梯度。

步骤30利用地震衰减梯度特征体确定目标层的岩性及烃类性质。

结合测井和地质等资料，利用井旁地震衰减梯度特征剖面确定不同岩性、流体等引起的地震衰减梯度特征体区别，再将其外推到无井区域，从而进行岩性、烃类检测。

本实施例的具体实现原理如下：

1、首先在利用测井、地质和合成地震记录等资料准确标定目标层的基础上，对工区内目标层的地震道逐道进行可变模态分解，确定最佳模态个数，计算生成的各个模态函数的时频图。可变模态分解是由DragomiretskiyandZosso于2014年提出的由下式表示的约束可变问题

$\underset{\left\{u_k\right\},\left\{{\mathrm{\ω}}_k\right\}}{\min }\left\{\underset{k}{\mathrm{\Σ}}{\left|\right|{\∂}_t[(\mathrm{\δ}\left(t\right)+\frac{j}{\mathrm{\πt}})*u_k\left(t\right)]e^{-j{\mathrm{\ω}}_{k}t}\left|\right|}_2^2\right\}---\left(1\right)$

其中，u_k是第k个模态函数，ω_k是第k个模态函数的中心频率。每个模态函数可表示为

${\hat{u}}_k\left(\mathrm{\ω}\right)=\frac{\hat{f}\left(\mathrm{\ω}\right)-{\mathrm{\Σ}}_{i\≠k}{\hat{u}}_i\left(\mathrm{\ω}\right)+(\widehat{\mathrm{\λ}}\left(\mathrm{\ω}\right)/2)}{1+2\mathrm{\α}{(\mathrm{\ω}-{\mathrm{\ω}}_k)}^2}---\left(2\right)$

对每一条地震道提取出来的各个模态分量u(t)，分别利用希尔伯特变换用下式提取瞬时幅度A(t)和瞬时频率ω(t)：

$\left\{\begin{array}{c}A\left(t\right)=\sqrt{u^2\left(t\right)+y^2\left(t\right)}\\ \mathrm{\φ}\left(t\right)=\arctan (y\left(t\right)/u\left(t\right))\\ \mathrm{\ω}\left(t\right)=\frac{1}{2\mathrm{\π}}\frac{\mathrm{d\φ}\left(t\right)}{\mathrm{dt}}\end{array}\right.---\left(3\right)$

其中，H[·]表示Hilbert变换，P为柯西主值。

为了避免(3)式中计算瞬时频率时相位解卷绕导致的模糊度，采用下式计算瞬时频率ω(t)：

$\mathrm{\ω}\left(t\right)=\frac{1}{2\mathrm{\π}}\frac{u\left(t\right)\frac{\mathrm{dy}\left(t\right)}{\mathrm{dt}}-\frac{\mathrm{du}\left(t\right)}{\mathrm{dt}}y\left(t\right)}{u^2\left(t\right)+y^2\left(t\right)}---\left(4\right)$

利用时间和瞬时频率和瞬时幅度定义一个三维空间[t,ω(t),A(t)]。令

H(ω,t)＝Re{A(t)e^i∫ω(t)dt}(5)

其中，Re表示取结果的实部。

于是，三维空间通过将函数H(ω,t)转变成三个变量的函数[t,ω,H(ω,t)]来实现，其中，A(t)＝H[ω(t),t)]。从而，可以获得某条地震道的各个模态函数的联合时频分布。

2、计算各条地震道各个模态函数的衰减梯度。由于各个模态函数是窄带信号且具有不同的频带宽度，其频谱相对原始地震信号更为平滑，有利于衰减梯度的准确拟合，也有利于反映特定频带处的衰减情况。对某条地震道的每一个模态函数，沿着每个时间采样点从时频图中提取相应的频率-幅度谱，取对数，在对应的频率-对数幅度谱中开一个可变长度频率域窗函数W(f)，可变长度的频率域窗函数W(f)定义为

$W\left(f\right)=,\left\{\begin{array}{c}40,|f_{\mathrm{zr}}-f_{\max }|>40\\ (f_{\mathrm{zr}}-f_{\max }),20<|f_{\mathrm{zr}}-f_{\max }|<40\\ 20,|f_{\mathrm{zr}}-f_{\max }|<20\end{array}\right.---\left(6\right)$

其中，f_max为某个时间采样点处频率-对数幅度谱中最大幅度处的频率值；f_zr为频率-对数幅度谱中的第一个幅度过零点处的频率值。即如果最大幅度处的频率值与第一个幅度过零点处的频率值的差小于20，窗函数长度取为20；如果最大幅度处的频率值与第一个幅度过零点处的频率值的差大于20同时小于40，那么窗函数长度为他们两者之间的距离，反之，窗函数长度取为40。这里，频率域窗函数长度20和40为一个经验值，根据实际情况可调节。

然后利用最小二乘法拟合对数能量和频率的斜率，即衰减梯度。

3、采用相关加权系数对某条地震道计算的各个模态函数对应的衰减梯度进行加权求和，结果作为该地震道的衰减梯度。首先计算各条地震道分解后生成的模态分量与原始地震道的相关系数R，如果相关系数大于a(a>0.1,可根据实际情况选择较强相关阈值)，则加权系数为1；相关系数大于0.1小于a，则加权系数为10^-1；相关系数小于0.1，则加权系数为10^-2。即相关加权系数Wc定义为

$\mathrm{Wc}=\left\{\begin{array}{c}1,\left|R\right|\&GreaterEqual;a\\ {10}^{-1},0.1\&le;\left|R\right|<a\\ {10}^{-2},\left|R\right|\&le;0.1\end{array}\right.---\left(7\right)$

该操作主要目的在于识别出该地震道信号的主要贡献分量并加强该地震道的主要成分，减弱次要贡献分量的信息成分。

4、逐道逐点计算出整个地震数据体的衰减梯度特征体，确定目标层的岩性及烃类性质。结合测井和地质等资料，利用井旁地震衰减梯度剖面特征确定不同岩性、流体等引起的地震衰减梯度特征体的区别，再将其外推到无井区域，对于三维数据可以结合沿层切片、等时切片等进行岩性、烃类检测。

本实施例利用基于可变模态分解和希尔伯特变换的高分辨率时频分析方法通过最小二乘拟合法逐道提取地震波高频段能量与其对应频率的拟线性关系来计算各个时间采样点处的地震波能量衰减梯度。含油气层在地震波衰减梯度图像上表现为强幅值异常。通过对地震波衰减梯度图像的分析，可以有效地检测地层中是否存在烃类。

本实施例以四川盆地某地的须家河组一条二维叠后偏移过井剖面处理为例：该区域为致密砂岩储层。

图2为四川盆地某地的须家河一条二维叠后偏移过井剖面(目标层段)。其中，AC为测井声波时差曲线。图中A为井位置。这里，仅显示了目标层段的地震剖面和对应的测井曲线。测井曲线解释该层段存在一段含气储层和两段含水储层，测井曲线解释储层位置如图中所示。

图3为过井道及其VMD分解结果。

图4为利用本技术计算的时频谱和短时傅里叶变换、小波变换计算的时频谱的对比结果。其中，本技术计算的时频谱是经过了5*15高斯平滑的结果，平滑使得本技术计算的时频谱分辨率下降，但是与常规方法更好对比。从中可以看到，本技术计算的时频谱在时间上和频率上的分辨率都较高，对储层信息体现较好。

图5为利用本技术计算的过井道的地震衰减梯度归一化曲线。从中可以看出，在研究层位之间，气层衰减梯度值最大，水层有相对较大的衰减梯度值。从中可以看到，时间分辨率较高。

图6为利用小波变换结合最小二乘拟合法计算的过井道的地震衰减梯度归一化曲线。从中可以看出，在研究层位之间，该方法给出的统计性解释结果与本技术方法类似，气层衰减梯度值最大，水层有相对较大的衰减梯度值，气层和水层差异相对本技术图5中给出的结果偏小。该方法的时间分辨率相比本技术得到的时间分辨率较差。

图7为利用小波变换和两点斜率的常规衰减梯度估计方法计算的过井道的地震衰减梯度归一化曲线。从中可以看出，在研究层位之间，气层和水层存在较大的衰减梯度值，气层和水层差异不明显。在储层外存在较强衰减梯度值，不利于储层解释。该算法相对本技术得到的结果图5和目前存在的较好方法图6，效果最差。

图8为利用本技术计算的地震衰减梯度特征剖面。从中可以看出，气层所在区域衰减梯度较强，水层不是很明显。气层、水层区分较大。图像时间分辨率和空间分辨率较强。

图9为利用小波变换结合最小二乘拟合法计算的归一化的地震衰减梯度特征剖面。从中可以看出，气层和水层所在区域衰减梯度较强。气层、水层区分不大。该方法相对利用小波变换和两点斜率的常规衰减梯度估计方法对于烃类的解释效果有所提高，但相对本技术提出的衰减梯度估计方法，图像时间分辨率和空间分辨率较差，且气水区分效果差。

图10为利用小波变换和两点斜率的常规衰减梯度估计方法计算的归一化的地震衰减梯度特征剖面。从中可以看出，气层和水层所在区域存在较强的衰减梯度，气层和水层区分不大。该方法由于利用两点斜率的计算方法估计衰减梯度，准确性和烃类检测效果相对最差。

图11为苏里格盆地某气田的一条二维叠后偏移过含气井剖面。图中，well为含气井所在位置。椭圆所示区域为含气层。

图12为利用本技术计算的归一化的地震衰减梯度特征剖面。从中可以看出，气层所在区域衰减梯度最强，该地震衰减特征剖面给出了含气性的有利解释结果。

图13为利用小波变换结合最小二乘拟合法计算的归一化的地震衰减梯度特征剖面。从中可以看出，气层所在区域衰减梯度较强。与本技术计算的地震衰减特征剖面图12相比，图13给出的含气性统计性解释结果较差，能量聚集性及时空分辨率比本技术所得结果都较差。

虽然本发明所揭露的实施方式如上，但所述的内容只是为了便于理解本发明而采用的实施方式，并非用以限定本发明。任何本发明所属技术领域内的技术人员，在不脱离本发明所揭露的精神和范围的前提下，可以在实施的形式上及细节上作任何的修改与变化，但本发明的专利保护范围，仍须以所附的权利要求书所界定的范围为准。

Claims (5)

1.基于可变模态分解的地震波衰减梯度估计方法，其特征在于，所述方法包括：

A利用测井、录井和合成地震记录标定目标层；

B对地震数据逐道进行可变模态分解，结合最小二乘法提取模态函数的衰减梯度，并对所述模态函数的衰减梯度进行加权求和获得各地震道的衰减梯度值；

C利用地震衰减梯度特征体确定目标层的岩性及烃类性质。

2.如权利要求1所述的基于可变模态分解的地震波衰减梯度估计方法，其特征在于，所述步骤B具体包括：

对每一道地震记录进行可变模态分解后生成的各个模态函数，利用希尔伯特变换计算瞬时振幅和瞬时频率，进一步生成各个模态函数的时频图；

计算各模态函数每个时间采样点处对应的对数频谱图，进一步通过最小二乘拟合法计算各个模态函数在每一个时间采样点上的衰减梯度；

采用加权系数对各个模态函数的衰减梯度进行相关加权求和获得各地震道的衰减梯度值。

3.如权利要求2所述的基于可变模态分解的地震波衰减梯度估计方法，其特征在于，对地震道进行可变模态分解产生的不同模态函数分别结合希尔伯特变换并通过可变长度频率域窗函数利用最小二乘法计算衰减梯度。

4.如权利要求1所述的基于可变模态分解的地震波衰减梯度估计方法，其特征在于，所述步骤C具体包括：

根据测井和地质资料信息，并利用井旁地震衰减梯度剖面特征，确定不同岩性、流体引起的地震衰减梯度特征体的区别，再将不同的岩性、流体引起的不同的地震衰减梯度特征外推到无井区域，从而进行岩性、烃类检测。

5.如权利要求4所述的基于可变模态分解的地震波衰减梯度估计方法，其特征在于，采用井标定，结合沿层切片、等时切片或地层切片进行分析，确定岩性及烃类性质。