CN107315193B - 一种利用地震瞬时质心频率进行含气性检测的方法 - Google Patents

Abstract

本发明涉及天然气地球物理勘探领域，是一种利用地震数据来提取地震波瞬时质心频率，进而利用地震瞬时质心频率可靠地指示天然气存在的方法，具体为：(1)综合利用地质、测井及合成地震记录资料准确标定目标层；(2)利用基于CEEMD的地震瞬时质心频率方法沿着目标层开时窗提取地震瞬时质心频率属性体；(3)利用地震瞬时质心频率属性体确定目标层的含气性；本发明通过对地震波瞬时质心频率图像的分析，可以有效地进行含气性检测。本发明对于含气性检测的结果具有较高的时空分辨率，同时避免了常规的基于短时傅里叶分析、小波变换的谱分解方法需要大量分频剖面进行分析然后利用最佳的分频剖面进行烃类解释的情况。

Description

一种利用地震瞬时质心频率进行含气性检测的方法

技术领域

本发明涉及天然气地球物理勘探领域，是一种利用地震数据来提取地震波瞬时质心频率，进而利用地震瞬时质心频率可靠地指示天然气的存在。

背景技术

在天然气勘探中，如何从地震信号的瞬时属性中提取更多的有用信息，尤其是利用频率异常信息，并结合地质、测井等资料，来寻找有意义的天然气储集带是石油物探研究人员一直以来的追求目标，同时也是难点问题。含流体的岩石地层会造成地震波在传播过程中发生能量损失，在含气层的内部及其下部，地震波的能量会发生明显的高频衰减。谱分解技术是目前利用地震信号高频衰减异常从地震反射数据进行地质解释及油气指示的一种常用的有效的含气性预测技术。地震数据体经过谱分解可以产生一系列分频剖面，每个分频剖面都是所有地震数据体在某个特定频率或频率段处瞬时振幅的一种反映。不同频率处的振幅剖面能够体现不同尺度处的地质体的不同的响应特征，而含气性信息可能在某些频率处的分频剖面中得到强化反映，更容易被识别。通常含气区域会体现出“低频强能量，高频弱能量”的衰减特征。但是谱分解技术需要利用一系列地震分频剖面进行分析，然后选取合理的分频剖面进行解释，工作量较大。

此外，谱分解技术的核心是利用时频分析方法如短时傅里叶变换、S变换、小波变换等。相比于传统的地震信号时频分析方法，适用于非线性非平稳信号处理的基于经验模态分解(Empirical mode decomposition,EMD)的时频分析方法具有一些独特的优势，如不受测不准原理的限制，时频分辨率更高等。目前基于EMD的时频分析技术已经被广泛应用于天然气检测中，并取得了很好的效果。EMD可以把地震信号从高频到低频分解为有限个单频或窄带本征模态函数(Intrinsic mode function,IMF)的和。Huang等人(1998)认为只有本征模态函数的瞬时频率才有物理意义。虽然基于EMD时频分析方法的含气性检测技术较基于传统时频分析方法的含气性检测技术具有更好的时空分辨率和精确度，但是由于EMD方法本身存在模态混叠效应等问题，会导致EMD分解出的IMF缺乏物理意义和地质意义，从而降低含气性检测的准确性。

完备聚合经验模态分解(Complete Ensemble Empirical Mode Decomposition，CEEMD)是Torres等人(2011)提出的一种克服EMD方法中模态混叠效应的改进方法，它可以把地震信号分解的更彻底，误差在数值上可以忽略，得到的IMF分量更具有物理意义和地质意义。

发明内容

本发明的目的在于解决上述现有技术中存在的难题，提供一种新的地层含气性信息的频率异常检测方法，采用的是基于CEEMD的地震瞬时质心频率技术。通过定义一种新的地震属性：地震瞬时质心频率，进行天然气储层含气性检测，提高现有技术的精确度。

为实现上述目的，本发明采用的技术方案是：一种利用地震瞬时质心频率进行含气性检测的方法，包括如下步骤：

A、综合利用地质、测井及合成地震记录资料准确标定目标层；

B、对各个地震道进行CEEMD分解产生不同的IMF分量，对每个IMF分量，结合希尔伯特变换计算瞬时振幅和瞬时频率，再利用地震瞬时质心频率计算公式计算各个IMF分量的瞬时质心频率，通过相关加权系数对计算得到的各个IMF分量瞬时质心频率进行加权求和，获得相关地震道的瞬时质心频率值，逐道计算目标层范围内的地震数据的瞬时质心频率，获得相关工区目标层范围内的地震瞬时质心频率特征体,具体为：

B1、将各个地震道进行CEEMD分解，得到IMF分量

1)确定CEEMD分解过程中添加的高斯白噪声、集合成员的个数和最佳最大筛选迭代次数，添加的高斯白噪声的幅值a和集合成员的个数N应满足：

其中，ε_n是最终的误差标准偏差，高斯白噪声的幅值a取值范围为0.1～0.4，集合成员的个数N＞100；

2)CEEMD分解：在CEEMD分解过程中，使用EMD经过I次分解信号x(t)+ε₀w_i(t)(i＝1,2,Λ,I)以获得第一个模态，并计算：

其中，w_i(t)(i＝1,2,Λ,I)是不同的高斯白噪声序列，x(t)是原始地震信号；

在第一个阶段(k＝1)中，计算第一个余量r₁(t)：

分解r₁(t)+ε₁E₁(w_i(t))(i＝1,2,Λ,I),直到产生它们的第一个IMF分量，然后定义第二个IMF分量为：

E_j(·)表示产生第j个模态，

对于k＝2,Λ,K,计算第k次的余量：

分解r_k(t)+ε_kE_k(w_i(t))(i＝1,2,Λ,I)直到获得它们的第一个EMD模态，然后定义第(k+1)个模态为：

重复前述步骤计算第k个模态；

B2、计算每个地震道分解产生的IMF分量的瞬时属性

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

其中，H[·]表示Hilbert变换，P为柯西主值；为了避免(7)式中计算瞬时频率时相位解卷绕导致的模糊度，采用下式计算瞬时频率ω(t)：

B3、计算各条地震道的瞬时质心频率

瞬时质心频率f定义为：

其中，ω_k为每个IMF分量在每个时间采样点处的瞬时频率；A_k为每个IMF分量在每个时间采样点处的瞬时振幅；

其中，首先计算每条地震道各个IMF分量的瞬时质心频率，然后，采用相关加权系数对计算得到的地震道中各个IMF分量的瞬时质心频率进行加权求和，将加权求和结果作为该地震道的瞬时质心频率；其中，相关加权方案如下：

利用各条地震道分解后生成的IMF分量与原始地震道的相关系数R的大小对各个IMF分量获得的瞬时质心频率进行加权，相关加权系数Rc定义为：

式(10)对具有强相关的IMF分量计算得到的瞬时质心频率保持不变，对具有较强相关的IMF分量计算得到的瞬时质心频率进行10^-1的衰减，对具有弱相关的IMF分量计算得到的瞬时质心频率进行10^-2的衰减；

对计算出的最终该条地震道的瞬时质心频率进行归一化，采用如下归一化计算公式：

x2＝(x-min(x))/(max(x)-min(x))……(11)

其中，x为该条地震道的瞬时质心频率，x2为该条地震道的归一化瞬时质心频率，将结果归一化到区间[0,1]范围内，min(·)表示取数据的最小值，max(·)表示取数据的最大值；

C、结合地质、测井和试油气信息相关资料，利用井旁地震瞬时质心频率特征剖面确定不同岩性、流体引起的地震瞬时质心频率属性体区别，再将其特征外推到无井区域，从而进行岩性、含气性检测。

本发明的有益技术效果是：本发明的基于CEEMD的地震瞬时质心频率提取方法，具有如下特点，主要表现为：

(1)使用了克服模态混叠效应的CEEMD方法进行地震道的分解。通过CEEMD方法获得的各个IMF分量是具有不同频带宽度的窄带信号，且较常规EMD方法更具物理意义和地质意义。准确的IMF分量保证了通过它计算获得的瞬时频率具有物理意义和更为准确、明确的地质意义。

(2)对地震道进行CEEMD分解产生的各个不同IMF分量分别结合希尔伯特变换提取瞬时振幅和瞬时频率，再计算瞬时质心频率，通过相关加权系数对各个IMF分量计算的瞬时质心频率进行加权求和，弱化了地震信号中的次要贡献成分，突出了主要贡献成分，能够更有效地反映地震信号的微弱变化，加强了天然气地震响应特征，有利于含气性检测。

(3)地震瞬时质心频率属性体的计算可适用于二维或三维数据的剖面分析、沿层切片或者时间切片等分析，计算方式灵活多样。

(4)基于CEEMD的地震瞬时质心频率属性体估计方法较常规谱分解技术和衰减分析技术更为简单方便，同时具有较高的分辨率。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为含气地质衰减模型；

图2为模型的地震响应；

图3为模型的基于CEEMD的地震瞬时质心频率剖面；

图4为四川盆地某地的须家河组一条二维叠后偏移过含气井剖面；

图5为利用本技术计算的该地震剖面基于CEEMD的地震瞬时质心频率属性剖面；

图6为利用基于短时傅里叶变换的谱分解方法计算的该地震剖面的分频剖面的对比结果；

图7为利用基于小波变换的谱分解方法计算的该地震剖面的分频剖面的对比结果；

图8为四川盆地某地的须家河组另一条二维叠后偏移过气水井剖面；

图9为利用本技术计算的该叠后偏移过气水井剖面的地震瞬时质心频率属性剖面。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

实施例1

一种利用地震瞬时质心频率进行含气性检测的方法，包括如下步骤：

A、综合利用地质、测井及合成地震记录资料准确标定目标层；

B、对各个地震道进行CEEMD分解产生不同的IMF分量，对每个IMF分量，结合希尔伯特变换计算瞬时振幅和瞬时频率，再利用地震瞬时质心频率计算公式计算各个IMF分量的瞬时质心频率，通过相关加权系数对各个IMF分量计算的瞬时质心频率进行加权求和，获得相关地震道的瞬时质心频率值，逐道计算目标层范围内的地震数据的瞬时质心频率，获得相关工区目标层范围内的地震瞬时质心频率特征体；具体为：

B1、将各个地震道进行CEEMD分解，得到IMF分量

1)确定CEEMD分解过程中添加的高斯白噪声、集合成员的个数和最佳最大筛选迭代次数，添加的高斯白噪声的幅值a和集合成员的个数N应满足：

其中，ε_n是最终的误差标准偏差，白噪声的幅值a取值范围为0.1～0.4，集合成员的个数N＞100；

CEEMD方法使用特定的噪声添加在每次分解的过程中，计算唯一的剩余量以获得每个IMF分量，最终分解获得的IMF分量是完全的，具有数值可忽略的误差；通常，添加的白噪声幅度不宜过大，否则会在分解过程中引入虚假的IMF分量；添加的白噪声分量的幅度也不宜太小，否则它会导致模态混叠现象不易消除。

2)CEEMD分解：在CEEMD分解过程中，使用EMD经过I次分解信号x(t)+ε₀w_i(t)(i＝1,2,Λ,I)以获得第一个模态，并计算：

其中，w_i(t)(i＝1,2,Λ,I)是不同的高斯白噪声序列，x(t)是原始地震信号；

在第一个阶段(k＝1)中，计算第一个余量r₁(t)：

分解r₁(t)+ε₁E₁(w_i(t))(i＝1,2,Λ,I),直到产生它们的第一个IMF分量，然后定义第二个IMF分量为：

E_j(·)表示产生第j个模态，

对于k＝2,Λ,K,计算第k次的余量：

分解r_k(t)+ε_kE_k(w_i(t))(i＝1,2,Λ,I)直到获得它们的第一个EMD模态，然后定义第(k+1)个模态为：

重复前述步骤计算第k个模态；

B2、计算每个地震道分解产生的IMF分量的瞬时属性

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

其中，H[·]表示Hilbert变换，P为柯西主值；为了避免(7)式中计算瞬时频率时相位解卷绕导致的模糊度，采用下式计算瞬时频率ω(t)：

B3、计算各条地震道的瞬时质心频率

瞬时质心频率f定义为：

其中，ω_k为每个IMF分量在每个时间采样点处的瞬时频率；A_k为每个IMF分量在每个时间采样点处的瞬时振幅；

本发明中，首先计算每条地震道各个IMF分量的瞬时质心频率，然后，采用相关加权系数对计算得到的地震道中各个IMF分量的瞬时质心频率进行加权求和，将加权求和结果作为该地震道的瞬时质心频率；其中，相关加权方案如下：

利用各条地震道分解后生成的IMF分量与原始地震道的相关系数R的大小对各个IMF分量获得的瞬时质心频率进行加权，相关加权系数Rc定义为：

式(10)对具有强相关的IMF分量计算得到的瞬时质心频率保持不变，对具有较强相关的IMF分量计算得到的瞬时质心频率进行10^-1的衰减，对具有弱相关的IMF分量计算得到的瞬时质心频率进行10^-2的衰减；该操作可以加强该地震道的主要贡献成分，同时减弱次要贡献分量的信息成分。

对计算出的最终该条地震道的瞬时质心频率进行归一化，采用如下归一化计算公式：

x2＝(x-min(x))/(max(x)-min(x))……(11)

其中，x为该条地震道的瞬时质心频率，x2为该条地震道的归一化瞬时质心频率，将结果归一化到区间[0,1]范围内，min(·)表示取数据的最小值，max(·)表示取数据的最大值。

C、结合地质、测井和试油气信息相关资料，利用井旁地震瞬时质心频率特征剖面确定不同岩性、流体引起的地震瞬时质心频率属性体区别，再将其特征外推到无井区域，从而进行岩性、含气性检测。

应用实施例1

(1)利用地质、测井、速度等资料结合弥散黏滞方程一个建立地质衰减模型。表1是依据鄂尔多斯盆地某气田的资料建立的含气衰减地质模型的参数。其中，层④为含气层，层③为干层(地质模型如图1所示)。

表1模型的地质参数

其中，V_p是纵波速度，ρ是密度，ζ是弥散系数，η是黏滞系数。Q是衰减系数。

模型含气层厚度为75m，子波频率为40Hz，采样频率为512Hz。

(2)图2为模型的地质响应。

(3)图3为模型对应的基于CEEMD的瞬时质心频率剖面。如图所示，本技术很好的检测到了含气区域。

(4)图4为四川盆地某地的须家河组一条二维叠后偏移过井剖面。该区域为致密砂岩储层。图中上部的椭圆所示区域为含气区域，下部的椭圆所示区域为含水区域。含气井命名为well。

(5)利用本技术计算的基于CEEMD的地震瞬时质心频率剖面如图5所示。从图中可以看到，在上部分的椭圆所示的含气区域中存在强烈的瞬时质心频率异常值。在下部分的椭圆所示的含水区域有较强的瞬时质心频率异常值。排除岩性等其他地层影响因素，所提方法给出了一个含烃类的统计性解释结果。同时，所提方法对于气、水的衰减程度在检测上有所差异，含气区域的异常值更大，水层异常值较小。气层、水层区分较大。图像时间分辨率和空间分辨率较强。

(6)图6为地震剖面的利用基于短时傅里叶变换的谱分解方法计算的分频剖面的对比结果。图6a和6b所示为利用短时傅里叶变换分别提取的低频分频剖面和高频分频剖面。该数据主频为50Hz。这里，低频分频剖面为45Hz，高频分频剖面为65Hz。从图中可以看出，基于短时傅里叶变换的低频分频剖面在含气区域和含水区域都检测到了较强的强振幅异常，在高频分频剖面中在含气区域和含水区域强能量减弱，都表现出了“低频强能量、高频弱能量”的特征，排除岩性等其他地层影响因素，都给出了一个含烃类的统计解释结果。但是分频剖面时空分辨率较低，对气、水区域差异不大。

(7)图7为利用基于小波变换的谱分解方法计算的地震剖面的分频剖面的对比结果。图7a和7b所示为利用小波变换分别提取的低频分频剖面和高频分频剖面。这里，低频分频剖面为45Hz，高频分频剖面为65Hz。从图中可以看出，基于小波变换的低频分频剖面和高频分频剖面在含气区域和含水区域都呈现出了“低频强能量、高频弱能量”的特征，都给出了一个含烃类的统计解释结果。与图6基于短时傅里叶变换提取的分频剖面相比，基于小波变换的分频剖面的时空分辨率较高，但是从图6a可以看出，气、水区域差异不大。对比图6、图7利用常规谱分解技术与图5所提技术得到的结果，本技术给出的结果时空分辨率较好，对储层信息体现较好，对气、水的区分能力更好。同时，本技术不需要类似常规谱分解技术需要分析大量的分频剖面。

(8)图8为四川盆地某地的须家河组另一条二维叠后偏移过气水井剖面。图中的椭圆所示区域为含气水区域。

(9)图9为利用本技术计算的该叠后偏移过气水井剖面的地震瞬时质心频率属性剖面。从图中可以看到，在含气水区域存在较强的瞬时质心频率异常值。排除岩性等其他地层影响因素，所提方法给出了一个含烃类的统计性解释结果。

本技术领域技术人员可以理解，除非另外定义，这里使用的所有术语(包括技术术语和科学术语)具有与本发明所属领域中的普通技术人员的一般理解相同的意义。还应该理解的是，诸如通用字典中定义的那些术语应该被理解为具有与现有技术的上下文中的意义一致的意义，并且除非像这里一样定义，不会用理想化或过于正式的含义来解释。

本领域普通技术人员可以理解：实施例中的装置中的部件可以按照实施例描述分布于实施例的装置中，也可以进行相应变化位于不同于本实施例的一个或多个装置中。上述实施例的部件可以合并为一个部件，也可以进一步拆分成多个子部件。

最后所应说明的是：以上实施例仅用以说明而非限制本发明的技术方案，尽管参照上述实施例对本发明进行了详细说明，本领域的普通技术人员应该理解：依然可以对本发明进行修改或者等同替换，而不脱离本发明的精神和范围的任何修改或局部替换，其均应涵盖在本发明的权利要求范围当中。

Claims (1)

1.一种利用地震瞬时质心频率进行含气性检测的方法，其特征在于，包括如下步骤：

A、综合利用地质、测井及合成地震记录资料准确标定目标层；

B、对各个地震道进行CEEMD分解产生不同的IMF分量，对每个IMF分量，结合希尔伯特变换计算瞬时振幅和瞬时频率，再利用地震瞬时质心频率计算公式计算各个IMF分量的瞬时质心频率，通过相关加权系数对计算得到的各个IMF分量瞬时质心频率进行加权求和，获得相关地震道的瞬时质心频率值，逐道计算目标层范围内的地震数据的瞬时质心频率，获得相关工区目标层范围内的地震瞬时质心频率特征体,具体为：

B1、将各个地震道进行CEEMD分解，得到IMF分量

1)确定CEEMD分解过程中添加的高斯白噪声、集合成员的个数和最佳最大筛选迭代次数，添加的高斯白噪声的幅值a和集合成员的个数N应满足：

其中，ε_n是最终的误差标准偏差，高斯白噪声的幅值a取值范围为0.1～0.4，集合成员的个数N＞100；

2)CEEMD分解：在CEEMD分解过程中，使用EMD经过I次分解信号x(t)+ε₀w_i(t)(i＝1,2,Λ,I)以获得第一个模态，并计算：

其中，w_i(t)(i＝1,2,Λ,I)是不同的高斯白噪声序列，x(t)是原始地震信号；

在第一个阶段(k＝1)中，计算第一个余量r₁(t)：

分解r₁(t)+ε₁E₁(w_i(t))(i＝1,2,Λ,I),直到产生它们的第一个IMF分量，然后定义第二个IMF分量为：

E_j(·)表示产生第j个模态，

对于k＝2,Λ,K,计算第k次的余量：

分解r_k(t)+ε_kE_k(w_i(t))(i＝1,2,Λ,I)直到获得它们的第一个EMD模态，然后定义第(k+1)个模态为：

重复前述步骤计算第k个模态；

B2、计算每个地震道分解产生的IMF分量的瞬时属性

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

其中，H[·]表示Hilbert变换，P为柯西主值；为了避免(7)式中计算瞬时频率时相位解卷绕导致的模糊度，采用下式计算瞬时频率ω(t)：

B3、计算各条地震道的瞬时质心频率

瞬时质心频率f定义为：

其中，ω_k为每个IMF分量在每个时间采样点处的瞬时频率；A_k为每个IMF分量在每个时间采样点处的瞬时振幅；

其中，首先计算每条地震道各个IMF分量的瞬时质心频率，然后，采用相关加权系数对计算得到的地震道中各个IMF分量的瞬时质心频率进行加权求和，将加权求和结果作为该地震道的瞬时质心频率；其中，相关加权方案如下：

利用各条地震道分解后生成的IMF分量与原始地震道的相关系数R的大小对各个IMF分量获得的瞬时质心频率进行加权，相关加权系数Rc定义为：

式(10)对具有强相关的IMF分量计算得到的瞬时质心频率保持不变，对具有较强相关的IMF分量计算得到的瞬时质心频率进行10^-1的衰减，对具有弱相关的IMF分量计算得到的瞬时质心频率进行10^-2的衰减；

对计算出的最终该条地震道的瞬时质心频率进行归一化，采用如下归一化计算公式：

x2＝(x-min(x))/(max(x)-min(x))……(11)

其中，x为该条地震道的瞬时质心频率，x2为该条地震道的归一化瞬时质心频率，将结果归一化到区间[0,1]范围内，min(·)表示取数据的最小值，max(·)表示取数据的最大值；

C、结合地质、测井和试油气信息相关资料，利用井旁地震瞬时质心频率特征剖面确定不同岩性、流体引起的地震瞬时质心频率属性体区别，再将其特征外推到无井区域，从而进行岩性、含气性检测。