基于二阶挤压小波变换的地震资料时频分析方法与流程

【技术领域】

本发明属于地震勘探技术领域，涉及一种基于二阶挤压小波变换的地震资料时频分析方法。

背景技术：

由于地震波传播过程中受到地下介质的衰减和滤波作用，地震信号是典型的非平稳信号。传统的fourier变换不具备对信号频率的局部刻画能力，而时频分析可以用来描述地震信号随时间的变化规律，刻画地震信号的局部特征，从而反映这些特征所对应的地质结构和储层信息，因此时频变换方法是分析非平稳信号的理想工具。构造或者选择合适的时频分析工具是地震信号处理和属性提取的关键问题之一。随着近代信号处理理论的发展，涌现出大量的时频分析工具和时频分析方法，其中很多也被应用到地震信号分析中。

1946年，gabor提出了短时傅里叶变换(short-timefouriertransform,stft)，短时傅里叶变换的优点是算法实现容易，物理意义明确，反映整体时频趋势。加窗傅里叶变换中窗函数的选择是时间分辨率和频率分辨率的折中，一旦窗函数选定，时频图的时间分辨率和频率分辨率也随之确定。作为典型的非平稳信号，地震信号的频率随着时间不断变化，需要进行多尺度分析。为了解决加窗傅里叶变换存在的不足，小波变换应运而生。20世纪80年代后期，morlet等完善了小波变换理论。连续小波变换继承了短时傅里叶变换的优势，实现了变窗处理的要求，具有良好的时频局域化特性。stockwell等于1996年提出了s变换，s变换是一种线性、多分辨率、无损可逆的时频分析方法，综合了短时傅里叶变换和小波变换的优点且避免了他们的不足。但是，s变换的窗函数是以固定的趋势随频率变化的，不能根据应用需求而改变窗函数。通过对标准s变换中的窗函数按照一定的规则进行替换，就可以得到各种广义形式的s变换表达式。具有不同窗函数的广义s变换对不同信号的刻画能力不同，在解决包括地震资料处理在内的各种实际问题时得到了很好的应用。

但是，以上的线性时频分析方法由于受不确定性原理的制约，同时分离在时域上距离太小的高频分量及在频域上相距太近的低频信号仍很困难。为了克服不确定性原理的限制，许多学者将已有的时频分析方法进行发展和推广，提出了很多新的时频分析方法。其中应用比较广泛的为重排算法，该方法通过将每个时频系数重排至其对应的能量重心上，从而获得具有更高分辨率的时频结果，但此类算法的计算量比较大，而且重排以后的时频结果不能对原信号进行重构。daubechies等人在2011年提出了同步挤压变换(synchrosqueezingtransform),改方法在小波变换的基础上，引入了基于瞬时频率的重排算法，将扩散的小波系数在频率方向挤压到其对应的真实频率上，不仅大大提高了信号的时频分辨率，而且能对信号进行精确重构。但是该算法也存在一定的局限性，仅仅在分析由若干窄带信号叠加的多分量信号时可以取得比较好的结果。

技术实现要素：

针对现有技术的不足，本发明目的在于提供一种高精度的基于二阶挤压小波变换的地震资料时频分析方法，分析结果受噪声的影响较小，且易于实现，计算效率高。

为达到上述目的，本发明是采取如下技术方案予以实现的：

基于二阶挤压小波变换的地震资料时频分析方法，二阶同步挤压小波变换包括以下步骤：

1)对信号进行连续小波变换，得到小波变换谱w(a,b)，w(a,b)的时域和频域表达式分别为

其中，ψ(t)为母小波函数，a和b分别表示尺度和平移因子，x(ω)和ψ(ω)分别为x(t)和ψ(t)的傅里叶变换，*表示取共轭；基本小波ψ(t)平方可积，无直流分量，并且满足以下容许性条件：

2)计算每一个小波系数对应的参考频率

其中

3)将时间-尺度域的小波系数重排至时间-频率域，二阶同步挤压小波变换的形式为

δ(g)表示单位冲击函数；通过公式(4),在时间-尺度域所有和瞬时频率对应的小波系数沿着尺度方向进行纵向求和，然后将求和的结果放置在时间-频率域w所对应的位置；

基于以上二阶挤压小波变换的地震资料时频分析方法具体包括以下步骤：

步骤1：在三维数据体中的某一地震剖面选取典型地震道，并对该地震数据进行二阶挤压小波变换，并找出异常区域所对应的频率成分；

步骤2：对整个地震剖面逐道进行二阶同步挤压小波变换，提取出于异常区域相对应的频率切片；

步骤3:对整个三维数据进行二阶挤压小波变换，得到对应的频率数据体，根据解释人员提供的层位信息，制作沿层切片，供地质人员解释。

进一步，所述步骤3)通过步骤2)得到每个小波系数对应的瞬时频率后，通过能量重排，将时间-尺度域的小波系数挤压到时间-频率域中；

二阶同步挤压小波变换有以下几种形式：

分布形式

近似形式

在公式(5)中，h(g)表示无限光滑并带有局部支撑的函数，且ε表示为了防止出现数值计算不稳定所设置的阈值参数，通过公式(4)或(5),在时间-尺度域所有和瞬时频率对应的小波系数沿着尺度方向进行纵向求和，然后将求和的结果放置在时间-频率域w所对应的位置。

进一步，所述步骤3)二阶同步挤压小波变换的离散形式为

本发明具有以下有益效果：

本发明基于二阶同步挤压小波变换的地震信号时频分析方法，首次在连续小波变换的基础上，引入二阶同步挤压算子，该变换通过对小波变换的系数在时间-尺度域进行重排，从而压制了能量扩散，该方法可以得到一种具有更高时频分辨率的分析结果。相比于之前的同步挤压变换，改方法更加适合对地震信号这种带有较强频率调制的信号进行分析，进一步将其应用于实际地震数据的时频分析中，可以更加准确的界定地下碳氢储层的边界，更精细的刻画河道与断层等地质信息，从而有利于地震资料的进一步解释和井位确定。

本发明基于二阶同步挤压小波变换的地震信号时频分析方法，合成信号表明，相比于传统的时频分析方法，如小波变换、同步挤压小波变换，二阶同步挤压小波变换通过对小波系数的重排，获得了一种具有更高时频分辨率的时频分析结果，并且能够准确刻画信号的频率随时间的变化关系。通过对二维地震剖面提取单频切片，二阶同步挤压小波变换的结果可以更加清晰的刻画出储层的位置和边界。在三维地震数据体的沿层切片上，二阶同步挤压小波变换的结果可以更加清晰的刻画地下断层和河道的特征。

【附图说明】

图1为测试信号的几种时频分析结果；

(a)测试信号；(b)连续小波变换变换；(c)同步挤压小波变换；(d)二阶同步挤压小波变换

图2为某油气田二维剖面；

图3为图2二维剖面的35hz频率切片；

(a)小波变换结果；(b)同步挤压变换结果；(c)二阶同步挤压小波变换结果；

图4为某油田三维数据体35hz沿层切片示意图；

(a)同步挤压小波变换35hz沿层切片；(b)二阶同步挤压小波变换35hz沿层切片。

【具体实施方式】

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

本发明所提出的基于二阶同步挤压小波变换的地震资料时频分析方法具体包括以下步骤：

二阶同步挤压小波变换包括以下步骤：

步骤1：对信号进行连续小波变换，得到小波变换谱w(a,b)，w(a,b)的时域和频域表达式分别为

其中，ψ(t)为母小波函数，a和b分别表示尺度和平移因子。x(ω)和ψ(ω)分别为x(t)和ψ(t)的傅里叶变换，*表示取共轭。

步骤2：计算每一个小波系数对应的参考频率

其中

我们以线性chirp信号h(t)＝acos(2πf0t²)为例，可以得到当|w(a,b)|≠0时，因此，当分析信号为线性chirp信号是，经过公式(3)的操作，可以从小波变换的系数中将信号的瞬时频率解调出来。

步骤3：将时间-尺度域的小波系数重排至时间-频率域；

通过步骤2得到每个小波系数对应的瞬时频率后，通过能量重排，可以将时间-尺度域的小波系数挤压到时间-频率域中。

二阶同步挤压小波变换有以下几种形式：

分布形式

近似形式

在公式(4)中，δ(g)表示单位冲击函数；在公式(5)中，h(g)表示无限光滑并带有局部支撑的函数，且ε表示为了防止出现数值计算不稳定所设置的阈值参数。通过公式(4)或(5),在时间-尺度域所有和瞬时频率对应的小波系数沿着尺度方向进行纵向求和，然后将求和的结果放置在时间-频率域w所对应的位置。

在进行数值计算时，需要对尺度a进行离散，离散方式为aj＝a02^j/nv，a0＝1/t，t为信号的时间长度,nv一般取32或64.尺度间隔aj-aj-1＝(δa)j；然后对频率进行划分wj＝δw×j,其中δw＝1/δt/n，δt,n分别为信号的采样间隔和采样点数。二阶同步挤压小波变换的离散形式为

本发明的物质基础为三维地震数据体，采用都是逐道处理方法。

本发明基于二阶同步挤压小波变换的地震资料时频分析方法具体步骤为：

步骤1:在三维地震数据体中选取一些典型地震道,首先进行频谱分析，找到该地震数据体对应的主频。

步骤2:选取三维数据体的某一地震剖面，逐道进行二阶同步挤压小波变换，然后制作主频所对应的频率切片。

步骤3：对整个三维数据体逐道进行二阶同步挤压小波变换，提取主频附近的频率切片，得到频率数据体；

步骤4：根据地质人员解释的层位信息，制作沿层切片。

效果分析

一、数值实验

首先，对比二阶同步挤压小波变换与常用的小波变换和同步挤压小波变换对合成信号时频分析的结果。合成信号如图1(a),此信号由以下三个分量组成：

s1(t)＝sin(2π(330t+16cos(3πt)))

s2(t)＝asin(160πt)

图1(b)-(d)分别对应小波变换、同步挤压小波变换和二阶同步挤压小波变换的结果。尽管这三种变换的结果都可以把这两个分量，但是由于受到不确定性原理的限制，小波变换结果出现的明显的能量扩散。对比图1(c)和(d)可以看出,通过能量重排，同步挤压小波变换和二阶同步挤压小波变换的结果大大提高了时频分辨率。值得注意的是，在图1(c)和(d)中，在信号的瞬时频率变化较快的地方(黑色箭头)，二阶同步挤压小波变换的时频聚集性比同步挤压小波变换的结果好，从而能更加准确的刻画信号的时变特性。

二、实际地震资料

下面，我们将二阶同步挤压小波变换应用实际地震资料时频分析，并对比了小波变换及二阶同步挤压小波变换的结果。图(2)为某海上油气田的二维剖面，该剖面含有1000道，每道的采样点为751，采样间隔为2ms。对每道地震数据进行时频变换，从中提取出35hz的频率切片，如图(3)所示。由于二阶同步挤压小波变换有着更高的时频分辨率，可以更加准确的定位储层的位置。图(3)中用黑色方框指示了储层的位置。最后，对某工区的三维数据体进行时频分析。分别采用同步挤压小波变换和二阶同步挤压小波变换对该数据体进行分析，首先得到一个35hz的频率数据体，根据解释人员提供的层位数据，制作沿层切片。如图(4)所示，可以看出地下古河道的特征在同步挤压小波变换的沿层切片上变得比较模糊，由于具有更高的时频分辨率，在二阶同步挤压小波变换的沿层切片上，河道(黑色箭头)的特征得到十分清晰的刻画。

以上实施例仅用于说明本发明的技术方案而非对其限制，尽管参照上述实施例对本发明进行了详细说明，所属领域的普通技术人员依然可以对本发明的具体实施方案进行修改或者等同替换，而这些并未脱离本发明精神和范围的任何修改或者等同替换，其均在本发明的权利要求保护范围之内。