CN114563824A - 一种二阶多重同步挤压多项式chirplet变换薄储层识别方法 - Google Patents

Abstract

本发明公开了一种二阶多重同步挤压多项式chirplet变换薄储层识别方法，包括以下步骤：S1、输入待分析的地震数据；S2、对多项式chirplet变换依次求取时间和频率方向的偏导，得到二阶瞬时频率估计算子；并将该估计算子进行多次迭代得到一种新的二阶多重瞬时频率估计函数；S3、利用二阶多重同步挤压多项式chirplet变换(MSSPCT)对地震数据进行时频变换，得到最终的时频结果；S4、对上述时频结果取模，得其高精度时频谱；S5、采用傅里叶变换获取过井道的地震道振幅谱，将振幅最大值设置为低频段，将振幅相对较低处设置为高频段；S6、从MSSPCT得到的时频谱中分别抽取高频段和低频段对应的频谱，通过两个共频率剖面中地震信号的衰减情况来表征该地震资料中的薄储层展布特征。

Description

一种二阶多重同步挤压多项式chirplet变换薄储层识别方法

技术领域

本发明属于地震信号处理与地震储层识别领域，具体涉及一种二 阶多重同步挤压多项式chirplet变换薄储层识别方法。

背景技术

时频分析是挖掘地震信号瞬时谱信息的重要方法，可用于储层的 精细识别。根据发展历程，时频分析方法可分为传统时频分析方法和 高精度时频分析方法两大类。传统时频分析方法包括短时傅里叶变换 (STFT)、小波变换(WT)、S变换(ST)及广义S变换(GST)等，已被广泛应用于储层预测、烃类检测等领域；相较于传统方法，高精 度时频分析方法显著提高了时频分辨率，主要包括经验模态分解 (EMD)、同步挤压方法(SST)、同步提取方法(SET)等，在对实 际地震资料进行分频处理时，对信号刻画的精度更高，可以更为有效 地预测和识别薄储层。

除了前面提到的时频分析方法以外，chirplet变换(CT)是另一种 参数化的时频分析方法。CT对于线性调频信号能够获得很好的分析 效果。然而，当待分析信号瞬时频率为非线性函数时，例如强调频信 号，CT就无法准确获取该类信号的瞬时频率估计。为弥补这一缺点， 多项式Chirplet变换(PCT)通过引入频率旋转因子和频率平移因子， 将CT中的线性核函数改变为多项式非线性核函数，就可以很好的分 析瞬时频率非线性信号，产生较好的时频表征结果。然而，PCT的时 频分辨率也受到测不准原理的影响，时频谱分辨率达不到最优的效果。

发明内容

本发明的目的在于针对现有技术中的上述不足，提供的一种二阶 多重同步挤压多项式chirplet变换薄储层识别方法。本发明能够显著 提高时频分辨率，对薄储层的识别提供强有力的帮助。

为了达到上述发明目的，本发明采用的技术方案为：

S1、输入待分析信号x(t)；

S2、在该信号模型下，对多项式chirplet变换先进行时间方向的偏 导，再对其沿频率方向求偏导，从而得到二阶瞬时频率的初步估计 算子；根据上述瞬时频率估计算子，并以此类推，得到多次迭代后 一种新的二阶多重瞬时频率估计函数；

S3、根据同步挤压原理，在时频域上构建一种以信号瞬时频率为 中心的同步挤压算子，用于挤压原时频谱瞬时频率曲线处的时频系 数，得到二阶多重同步挤压多项式chirplet变换值；

S4、对二阶多重同步挤压多项式chirplet变换的时频结果取模，得 到每一点处的时频能量，从而得其高精度时频谱；

S5、采用傅立叶变换获取过井道的地震道振幅谱，将振幅最大值 处设置为低频段，将振幅相对较小值处设置为高频段；

S6、从二阶多重同步挤压多项式chirplet变换得到的时频谱中分别 抽取低频段和高频段对应的频谱，得到两个共频率剖面；通过对比两 个共频率剖面中地震信号的衰减情况识别该地震资料中的薄储层展 布特征。

优选地，步骤S2中对多项式chirplet变换先进行时间方向求偏 导，得到:

其中，

表示求关于时间变量t的偏导数，

为在窗口 g(t)和信号s(u)下的PCT，t为时间变量，ω为频率变量，A(t)和

分别表示信号的瞬时振幅和瞬时相位函数，A⁽¹⁾(t)为一阶导数， A⁽²⁾(t)为二阶导数，

为在窗口tg(t)和信号s(u)下的PCT。

为了得到

对

沿频率方向求偏导：

其中，

表示先求关于时间变量t的一阶偏导数，再求关于频率 变量ω的二阶偏导数。

优选地，步骤S2中得到二阶瞬时频率估计算子

为：

n＝1时：

其中，

表示二阶瞬时频率的初步估计算子，

表示求关于频率变量ω的偏导数，

为一重二阶瞬时频率估计算 子，I{·}表示取复数的虚部。

优选地，步骤S2中得到n次迭代后一种新构建的二阶多重瞬时 频率估计函数

其中，

表示为n重的二阶瞬时频率估计函数；

特别地，当n＝1时，一重二阶瞬时频率估计值即为所述步骤S2 中得到的二阶瞬时频率的初步估计值。

优选地，步骤S3中得到的与二阶多重瞬时频率估计函数相对应 的重排方法，即二阶多重同步挤压多项式chirplet变换值 MSSPCT^[n](t,η)为：

其中，η为频率变量，

表示MSSPCT为了进一步 提高PCT的能量集中度所采用的瞬时频率估计函数，使其越来越接 近真实频率。

优选地，所述步骤S3中所提出的二阶多重同步挤压多项式 chirplet变换是无损逆变换，不存在信息丢失，可以重构信号x(t)，其 逆变换表达式为：

本发明提供的一种二阶多重同步挤压多项式chirplet变换薄储层 识别方法，具有以下有益效果：

本发明的二阶多重同步挤压多项式chirplet变换相比于多项式 chirplet变换，是在利用时间和频率的二阶偏导数对同步挤压多项式 chirplet变换瞬时频率进行修正的基础上，采用迭代重排的方式逐步 将模糊的时频能量聚焦到相应的时频脊线，使信号的真实瞬时频率得 以凸显，且更加精细。在多项式chirplet变换的基础上能够显著提高 时频能量聚焦性，提高薄储层的识别精度。

附图说明

图1为本发明流程图；

图2为合成信号；

图3为对合成信号进行短时傅里叶变换(STFT)得到的时频谱；

图4为对合成信号进行多项式chirplet变换(PCT)得到的时频 谱；

图5为对合成信号进行一重二阶同步挤压多项式chirplet变换 (MSSPCT)得到的时频谱；

图6为对合成信号进行五重同步挤压变换(MSST)得到的时频 谱；

图7为对合成信号进行五重二阶同步挤压多项式chirplet变换得 到的时频谱；

图8为四川某区域地震剖面；

图9为过井道的地震道振幅谱；

图10(a)和(b)分别为STFT对实际地震数据处理后生成的低频和 高频剖面图，(c)和(d)分别为MSST对实际地震数据处理后生成的低 频和高频剖面图，(e)和(f)分别为SPCT对实际地震数据处理后生成的 低频和高频剖面图，(g)和(h)分别为MSSPCT对实际地震数据处理后 生成的低频和高频剖面图；

图11为图10低频剖面白色方框的局部放大图。

具体实施方式

下面将结合附图对本发明作进一步说明。

实施例1：参见图2-图7，合成信号是检验时频分布的时频聚焦 性能好坏的公认模型，采用图2展示的合成信号来进行试验，包括以 下步骤：

(1)输入待分析的合成信号；

(2)在该信号模型下，对多项式chirplet变换先进行时间方向的 偏导，再对其沿频率方向求偏导，得其二阶瞬时频率的初步估计算子；

其中，对时间方向求偏导，得到:

其中，

表示求关于时间变量t的偏导数，

为在窗口 g(t)和信号s(u)下的PCT，t为时间变量，ω为频率变量，A(t)和

分别表示信号的瞬时振幅和瞬时相位函数，A⁽¹⁾(t)为一阶导数， A⁽²⁾(t)为二阶导数，

为在窗口tg(t)和信号s(u)下的PCT。

为了得到

对

沿频率方向求偏导：

其中，

表示先求关于时间变量t的一阶偏导数，再求关于频率 变量ω的二阶偏导数。

联合公式(7)和公式(8)得到二阶瞬时频率的初步估计算子

为：

n＝1时：

其中，

表示二阶瞬时频率的初步估计算子，

表示求关于频率变量ω的偏导数，

为一重二阶瞬时频率估计算 子，I{·}表示取复数的虚部。

(3)根据上述瞬时频率估计算子，并以此类推，得到n次迭代 后一种新构建的二阶多重瞬时频率估计函数：

其中，

表示为n重的二阶瞬时频率估计函数。

特别地，当n＝1时，一重二阶瞬时频率估计值即为所述步骤(2) 中得到的二阶瞬时频率的初步估计值。

(4)根据同步挤压原理，在时频域上构建一种以信号瞬时频率 为中心的同步挤压算子，用于挤压原时频谱瞬时频率曲线处的时频系 数，得到二阶多重同步挤压多项式chirplet变换值：

其中，η为频率变量，

表示MSSPCT为了进一步 提高PCT的能量集中度所采用的瞬时频率估计函数，使其越来越接 近真实频率。

通过对图3、图4、图5、图6、图7进行对比分析发现，STFT 表征的时频谱能量聚焦不高；与STFT相比，PCT的能量聚焦性更好， 且两个分量信号之间的干涉效应对其影响更小。

其次，多重同步挤压变换作为高精度时频分析方法，它通过重新 分配比例系数或频率系数来规避常规线性变换中固有的不确定性关 系，从而提高时频分析的清晰度和集中度，但其仍受强调频信号限制， 在强时变位置存在一定程度的发散。

最后，由图可见，二阶多重同步挤压多项式chirplet变换相比前 三种方法具有更高的时频分辨率，时频聚焦性更好，可以准确地刻画 时变信号的频率分布和变化情况。此外纯谐波信号的时频聚焦性比强 调频信号的聚焦性更好更准确。

实施例2：参见图8-图11，为验证本发明的实用性，实施例2使 用不同时频分析方法对四川某区域地震资料进行分析。该气藏共包含 304条地震道，采样点数为1126个，时间采样间隔为2ms。具体实施 步骤如下：

(1)输入待分析的地震数据，得到如图8所示的原始地震剖面 图；

(2)在该信号模型下，对多项式chirplet变换先进行时间方向的 偏导，再对其沿频率方向求偏导，从而得到二阶瞬时频率的初步估计 算子；

其中，对时间方向求偏导，得到:

其中，

表示求关于时间变量t的偏导数，

为在窗口 g(t)和信号s(u)下的PCT，t为时间变量，ω为频率变量，A(t)和

分别表示信号的瞬时振幅和瞬时相位函数，A⁽¹⁾(t)为一阶导数， A⁽²⁾(t)为二阶导数，

为在窗口tg(t)和信号s(u)下的PCT。

为了得到

对

沿频率方向求偏导：

其中，

表示先求关于时间变量t的一阶偏导数，再求关于频率 变量ω的二阶偏导数。

联合公式(12)和公式(13)得到二阶瞬时频率的初步估计算子

为：

n＝1时：

其中，

表示二阶瞬时频率的初步估计算子，

表示求关于频率变量ω的偏导数，

为一重二阶瞬时频率估计算 子，I{·}表示取复数的虚部。

(3)根据上述瞬时频率估计算子，并以此类推，得到n次迭代 后一种新构建的二阶多重瞬时频率估计函数：

其中，

表示为n重的二阶瞬时频率估计函数。

特别地，当n＝1时，一重二阶瞬时频率估计值即为所述步骤(2) 中得到的二阶瞬时频率的初步估计值。

(4)根据同步挤压原理，在时频域上构建一种以信号瞬时频率 为中心的同步挤压算子，用于挤压原时频谱瞬时频率曲线处的时频系 数，得到二阶多重同步挤压多项式chirplet变换值：

其中，η为频率变量，

表示MSSPCT为了进一步 提高PCT的能量集中度所采用的瞬时频率估计函数，使其越来越接 近真实频率。

(5)对地震数据进行二阶多重同步挤压多项式chirplet变换的时 频结果取模，得到每一点时频能量，从而得其高精度时频谱；

(6)采用傅立叶变换获取过井道的地震振幅谱，将振幅最大值 处设置为低频段，将振幅相对较小值处设置为高频段，如图9所示；

(7)从二阶多重多项式chirplet变换得到的时频谱中分别抽取低 频段和高频段对应的频谱，得到两个共频率剖面；通过对比两个共频 率剖面中地震信号的衰减情况识别该地震资料中的薄储层展布特征；

通过对比不同时频分析方法的高低频剖面可得，低频剖面上能量 分布较强，而在高频剖面能量明显衰减。符合地震波在油气藏中传播 的特征。图10中可以看到在白色方框内，低频段区域能量更强大， 高频段区域有明显的能量衰减。这种现象是符合客观事实的。

综合以上分析表明，四种方法均可用于识别储层。但在刻画薄层 含气性方面存在一定差异。为了突出显示图10中第一列低频剖面白 色方框区域内的更多细节，我们将其放大如图11所示。图11中可以 看出，STFT可观察到异常衰减，但由于分辨率较低，对于薄储层只 能提供较为粗糙的时频信息。与STFT相比，作为高精度时频分析方 法的MSST、SPCT和MSSPCT在储层垂直圈定的刻画精度方面均得 到了提高，如图10(b)，(c)和(d)所示。与其他三种时频分析方法相比， MSSPCT可以更有效地定位薄层，对薄层的上下边界描述更加清晰，分 辨率精度明显提高，能量也更加聚焦。总之，该方法更有效的刻画了 薄储层。

本发明提出的二阶多重同步挤压多项式chirplet变换(MSSPCT)， 在利用时间和频率的二阶偏导数对同步挤压多项式chirplet变换的瞬 时频率计算式进行修正的基础上，采用迭代重分配的方式逐步将模糊 能量集中到相应的时频脊线上。

本发明以多重同步挤压变换与同步挤压多项式chirplet变换为基 础，结合了多重同步挤压变换与同步挤压多项式chirplet变换的优点， 具有更高的时频分辨率。

虽然结合附图对发明的具体实施方式进行了详细地描述，但不应理解 为对本专利的保护范围的限定。在权利要求书所描述的范围内，本领 域技术人员不经创造性劳动即可做出的各种修改和变形仍属本专利 的保护范围。

Claims (6)

1.一种二阶多重同步挤压多项式chirplet变换薄储层识别方法，其特征在于，包括以下步骤：

S1、输入待分析信号x(t)；

S2、在该信号模型下，对多项式chirplet变换先进行时间方向的偏导，再对其沿频率方向求偏导，从而得到二阶瞬时频率的初步估计算子；根据上述瞬时频率估计算子，并以此类推，得到多次迭代后一种新的二阶多重瞬时频率估计函数；

S3、根据同步挤压原理，在时频域上构建一种以信号瞬时频率为中心的同步挤压算子，用于挤压原时频谱瞬时频率曲线处的时频系数，得到二阶多重同步挤压多项式chirplet变换值；

S4、对二阶多重同步挤压多项式chirplet变换的时频结果取模，得到每一点处的时频能量，从而得其高精度时频谱；

S5、采用傅立叶变换获取过井道的地震道振幅谱，将振幅最大值处设置为低频段，将振幅相对较小值处设置为高频段；

S6、从二阶多重同步挤压多项式chirplet变换得到的时频谱中分别抽取低频段和高频段对应的频谱，得到两个共频率剖面；通过对比两个共频率剖面中地震信号的衰减情况识别该地震资料中的薄储层展布特征。

2.根据权利要求1所述的一种二阶多重同步挤压多项式chirplet变换薄储层识别方法，其特征在于，所述步骤S2中对多项式chirplet变换先进行时间方向求偏导，得到:

其中，

表示求关于时间变量t的偏导数，

为在窗口g(t)和信号s(u)下的PCT，t为时间变量，ω为频率变量，A(t)和

分别表示信号的瞬时振幅和瞬时相位函数，A⁽¹⁾(t)为一阶导数，A⁽²⁾(t)为二阶导数，

为在窗口tg(t)和信号s(u)下的PCT；

为了得到

对

沿频率方向求偏导：

其中，

表示先求关于时间变量t的一阶偏导数，再求关于频率变量ω的二阶偏导数。

3.根据权利要求2所述的对多项式chirplet变换先进行时间方向求偏导，再对其沿频率方向求偏导，其特征在于，根据权利要求2所述,联合公式(1)和公式(2)得到二阶瞬时频率的初步估计算子

为：

n＝1时：

其中，

表示二阶瞬时频率的初步估计算子，

表示求关于频率变量ω的偏导数，

为一重二阶瞬时频率估计算子，I{·}表示取复数的虚部。

4.根据权利要求3所述的二阶瞬时频率的初步估计算子

其特征在于，根据权利要求3所述得到多次迭代后一种新的二阶多重瞬时频率估计函数

其中，

表示为n重的二阶瞬时频率估计函数；

特别地，当n＝1时，一重二阶瞬时频率估计值即为权力要求3中得到的二阶瞬时频率的初步估计值。

5.根据权利要求1所述的一种二阶多重同步挤压多项式chirplet变换薄储层识别方法，其特征在于，所述步骤S3中得到的与二阶多重瞬时频率估计函数相对应的重排方法，即二阶多重同步挤压多项式chirplet变换值MSSPCT^[n](t,η)为：

其中，η为频率变量，

表示MSSPCT为了进一步提高PCT的能量集中度所采用的瞬时频率估计函数，使其越来越接近真实频率。

6.根据权利要求1所述的一种二阶多重同步挤压多项式chirplet变换薄储层识别方法，其特征在于，所述步骤S3中所提出的二阶多重同步挤压多项式chirplet变换是无损逆变换，不存在信息丢失，可以重构信号x(t)，其逆变换表达式为：

Images