CN104730580B - 地震资料异常振幅压制方法 - Google Patents

Abstract

公开一种地震资料异常振幅压制方法，包括：提取地震道信号进行时频双域分解；分别在每个对应的频率下，根据有效波振幅分布，找到有效波振幅区；在有效波振幅区，计算有效振幅的能量指数作为平均振幅能量指数；计算采样点的能量指数，将其与平均振幅能量指数进行比较，若大于，则将采样点振幅划分为疑似异常振幅；对于疑似异常振幅，将其能量指数与平均振幅能量指数做差再与平均振幅能量指数相比，若该比值大于给定阀值则确认为该疑似异常振幅为异常振幅值；对异常振幅值进行衰减；以及将时频谱信号反变换到时间域，即完成异常振幅压制。上述方法克服常规方法异常振幅压制不彻底的问题保护有效信号，提高地震资料信噪比。

Description

地震资料异常振幅压制方法

技术领域

本发明涉及地震资料的数字信号处理方法，具体地，涉及一种地震资料异常振幅压制方法。

背景技术

在地震资料处理领域中，通常需要将异常振幅与有效波振幅进行分离。异常振幅通常具有随机，振幅强等特点，而且异常振幅频带通常与有效波频带有很大部分重叠，这使得压制异常振幅的同时保护有效信号的问题变得非常棘手，目前仍没有解决常规方法异常振幅压制不彻底的问题。

发明内容

本发明的目的克服常规方法异常振幅压制不彻底的问题，提供一种地震资料异常振幅压制方法。

为了实现上述目的，本发明提供一种地震资料异常振幅的压制方法，该方法包括：S101，提取地震资料的时域信号进行时频双域分解，以及分析分解出的时频谱以得出异常振幅占优频率；S102，对于任意一个所述异常振幅占优频率下的信号，进行有效波振幅能量指数分析寻找有效波振幅区，并且计算有效波振幅的能量指数作为平均能量指数；S103，在所述异常振幅占优频率下提取采样点，并计算所述采样点振幅的能量指数，当所述采样点振幅的能量指数大于所述平均能量指数时确定所述采样点具有疑似异常振幅；S104，将所述疑似异常振幅的能量指数数量级与所述平均能量指数数量级做差再与所述平均能量指数数量级相比以得到比值，当所述比值大于给定阈值时确定所述疑似异常振幅为异常振幅，并对所述异常振幅进行衰减；S105，对所有所述异常振幅占优频率下的时频谱信号进行S102、S103、以及S104步骤处理以得到压制异常振幅后的时频谱信号；以及S106，对所有所述压制异常振幅后的时频谱信号进行反变换以得到压制异常振幅后的地震资料。

优选地，所述时域信号进行时频双域分解的方法包括小波变换或S变换。

优选地，利用下述公式进行时域信号的时频双域分解：

S_k(t)为时域信号，W_gS_k(t,a)为时频谱，g(t)为复数母小波函数，(·)^*表示复共轭，a为频率尺度因子，其中，采用的复数Morlet母小波为

优选地，计算所述有效波振幅的能量指数的方法包括：对所述振幅能量指数取平均值或对所述振幅能量指数取均方根。

进一步优选地，利用下述公式计算得出所述有效波振幅的平均能量指数：

其中，THD为平均能量指数，W_gS_k(t,a)为时频谱，i为采样点，i＝1,2,3,…,N。

优选地，计算所述采样点振幅的能量指数的方法包括对振幅能量取自然对数或对振幅能量取常用对数。

进一步优选地，利用下述公式计算所述采样点振幅的能量指数：

EA_ik＝Ceil(lg(fabs(W_gS_k(i,a))))，

其中，EA_ik为所述采样点振幅的能量指数，W_gS_k(t,a)为时频谱，i为采样点，i＝1,2,3,…,N。

优选地，所述能量指数与有效波振幅的能量指数进行比较的方法包括：

EA_ik＞THD，则任意一个时频谱所在采样点振幅值为疑似异常振幅值；

EA_ik＜THD，则任意一个时频谱所在采样点振幅值为有效波振幅;

其中，EA_ik为所述采样点振幅的能量指数，THD为平均能量指数。

优选地，利用下述公式计算所述采样点的均值指数数量级振幅值：

W_gS_k(i,a)_av＝W_gS_k(i,a)·10^(THD-EA_ik)，

其中，THD为平均能量指数，EA_ik为采样点振幅的能量指数，W_gS_k(t,a)为时频谱，W_gS_k(i,a)_av为所述采样点的均指指数数量级振幅值。

优选地，所述确定异常振幅的步骤包括：当所述比值小于所述给定阈值，则确定该采样点振幅值为有效波振幅值；以及当所述比值大于所述给定阈值，则确定该采样点振幅值为异常振幅值。

通过上述技术方案，在优势频带，根据异常振幅强能量的特性进行异常振幅识别，联合异常振幅的能量指数特征，对异常振幅进行压制的办法，很好地解决了常规方法异常振幅压制不彻底的问题，能极大压制低信噪比资料中的异常振幅，同时保护与其频带重叠的有效波信号不被压制，提高地震资料信噪比，满足工业化生产的需求。

本发明的其他特征和优点将在随后的具体实施方式部分予以详细说明。

附图说明

附图是用来提供对本发明的进一步理解，并且构成说明书的一部分，与下面的具体实施方式一起用于解释本发明，但并不构成对本发明的限制。在附图中：

图1是本发明提供的一种地震资料异常振幅压制方法流程图；

图2（a）-图2（c）是利用本发明提供的方法对某地震资料异常振幅压制效果图；

图3（a）-图3（c）是利用现有技术中的方法对某地震资料异常振幅压制效果图；

图4（a）-图4（c）是利用本发明提供地震资料异常振幅压制方法对与图3中相同的地震资料异常振幅压制效果图；以及

图5是本发明提供的优选的实施方式的地震资料异常振幅压制方法流程图。

具体实施方式

以下结合附图对本发明的具体实施方式进行详细说明。应当理解的是，此处所描述的具体实施方式仅用于说明和解释本发明，并不用于限制本发明。

图1是本发明提供的一种地震资料异常振幅压制方法流程图。如图1所示，本发明提供一种地震资料异常振幅的压制方法，该方法包括：S101，提取地震资料的时域信号进行时频双域分解，以及分析分解出的时频谱以得出异常振幅占优频率；S102，对于任意一个所述异常振幅占优频率下的信号，进行有效波振幅能量指数分析寻找有效波振幅区，并且计算有效波振幅的能量指数作为平均能量指数；S103，在所述异常振幅占优频率下提取采样点，并计算所述采样点振幅的能量指数，当所述采样点振幅的能量指数大于所述平均能量指数时确定所述采样点具有疑似异常振幅；S104，将所述疑似异常振幅的能量指数数量级与所述平均能量指数数量级做差再与所述平均能量指数数量级相比以得到比值，当所述比值大于给定阈值时确定所述疑似异常振幅为异常振幅，并对所述异常振幅进行衰减；S105，对所有所述异常振幅占优频率下的时频谱信号进行S102、S103、以及S104步骤处理以得到压制异常振幅后的时频谱信号；以及S106，对所有所述压制异常振幅后的时频谱信号进行反变换以得到压制异常振幅后的地震资料。

其中，通过S101，可以实现对信号的时频分析，从而得出占优频率段，为下一步处理做准备；通过S102，得出占优频率段的有效波振幅能量指数作为平均能量指数，为下一步做准备；通过S103，分离出具有疑似异常振幅的采样点；通过S104，对于具有疑似振幅值的采样点，将其能量降为平均振幅能量级别，得到该采样点的平均能量指数振幅值，再通过一系列的比较变换可以确定其是否为异常振幅，并对异常振幅进行衰减以压制（例如，使其乘以衰减指数进行衰减）；通过S105，可以得到压制异常振幅后的时频谱信号；通过S106对衰减后的数据进行反变换，可以得到滤波结果。

优选地，所述时域信号进行时频双域分解的方法包括小波变换或S变换，两种方法都可以实现时频双域分解。

进一步优选地，以小波变换步骤为例，利用下述公式进行时域信号的时频双域分解：

S_k(t)为时域信号，W_gS_k(t,a)为时频谱，g(t)为复数母小波函数，(·)^*表示复共轭，a为频率尺度因子，其中，采用的复数Morlet母小波为

优选地，计算所述有效波振幅的能量指数的方法包括：对所述振幅能量指数取平均值或对所述振幅能量指数取均方根。

进一步优选地，以平均值为例，利用下述公式计算得出所述有效波振幅的平均能量指数：

其中，THD为平均能量指数，W_gS_k(t,a)为时频谱，i为采样点，i＝1,2,3,…,N。

优选地，计算所述采样点振幅的能量指数的方法包括对振幅能量取自然对数或对振幅能量取常用对数。

进一步优选地，以利用常用对数为例，下述公式计算所述采样点振幅的能量指数：

EA_ik＝Ceil(lg(fabs(W_gS_k(i,a))))，

其中，EA_ik为所述采样点振幅的能量指数，W_gS_k(t,a)为时频谱，i为采样点，i＝1,2,3,…,N。

优选地，所述能量指数与有效波振幅的能量指数进行比较的方法包括：

EA_ik＞THD，则任意一个时频谱所在采样点振幅值为疑似异常振幅值；

EA_ik＜THD，则任意一个时频谱所在采样点振幅值为有效波振幅;

其中，EA_ik为所述采样点振幅的能量指数，THD为平均能量指数。

优选地，利用下述公式计算所述采样点的均值指数数量级振幅值：

W_gS_k(i,a)_av＝W_gS_k(i,a)·10^(THD-EA_ik)，

其中，THD为平均能量指数，EA_ik为采样点振幅的能量指数，W_gS_k(t,a)为时频谱，W_gS_k(i,a)_av为所述采样点的均指指数数量级振幅值。

优选地，所述确定异常振幅的步骤包括：当所述比值小于所述给定阈值，则确定该采样点振幅值为有效波振幅值；以及当所述比值大于所述给定阈值，则确定该采样点振幅值为异常振幅值。

在不同的实施方式中，根据本发明的方法可以包括如下步骤：A、提取地震道信号进行时频双域分解；B、分别在每个对应的频率下，根据有效波振幅分布，找到有效波振幅区；C、在有效波振幅区，计算有效振幅的能量指数作为平均振幅能量指数；D、计算采样点的能量指数，将其与平均振幅能量指数进行比较，若大于，则将采样点振幅划分为疑似异常振幅；E、对于疑似异常振幅，将其能量指数与平均振幅能量指数做差再与平均振幅能量指数相比，若该比值大于给定阀值则确认为该疑似异常振幅为异常振幅值；F、对异常振幅值进行衰减；G、将时频谱信号反变换到时间域，即完成异常振幅压制。

上述方法中，方法所针对的原始数据为：一个采样点数为N的地震道，记为S_k(t)，其中k＝1,2,3,…,K。

以下，结合本发明提供的方法对于对Sk(t)进行面波压制的步骤具体描述如下：

A.对时域信号S_k(t)进行时频双域分解，得时频谱W_gS_k(t,a)。以小波变换为例，采用复数母小波函数g(t)，取其连续小波变换，为：

其中，a为尺度因子（频率），(·)^*表示复共轭，采用的复数Morlet母小波为

B.对时频谱W_gS_k(t,a)进行有效波范围分析（例如，找出有效波振幅区），计算出有效波的平均振幅能量指数THD；

其中，i为时间采样点，i＝1,2,3,…,N

C.根据采样点振幅能量指数与平均振幅能量指数比较大小的结果，将时采样点分为有效振幅和疑似异常振幅。对于某个频率下的信号，统计出采样点的振幅能量指数EA_ik。其计算公式如下：

EA_ik＝Ceil(lg(fabs(W_gS_k(i,a))))

D.如若：EA_ik＞THD，则判断采样点振幅值为疑似异常振幅值，将其指数数量级降为时窗内均值指数数量级，易得到该采样点的均值指数数量级振幅值：

W_gS_k(i,a)_av＝W_gS_k(i,a)·10^(THD-EA_ik)

E.对于疑似异常振幅，则进一步判断是否超过一定限度，以确认是否为异常振幅。判断标准可以为：该疑似异常振幅与均值振幅指数数量级的差值和均值振幅指数数量级之比是否超过给定的阀值STD，即达到足够大，才判定为异常振幅进行衰减。即：（为简洁，以下表述均假设振幅值为正值。）

(lg(W_gS_k(i,a))-lg(W_gS_k(i,a)_av))lg(W_gS_k(i,a)_av)＞STD

F.识别出异常振幅后，乘以给定的衰减指数AttenScale加以衰减。

W_gS_k(i,a)＝W_gS_k(i,a)_av·AttenScale

G.分别对每个频率下的时频谱信号进行B、C、D三步处理，得到压制异常振幅后的时频谱信号W′_gS_k(t,a)；

F.对压制异常振幅后的时频谱信号W′_gS_k(t,a)进行反变换，重构得到压制异常振幅后的地震道信号S'_k(t)。

通过本发明提供的方法，能够有效压制地震资料中的异常振幅，保护有效信号，提高地震资料的信噪比。

下面结合附图对使用根据本发明实施方式提供的地震资料异常振幅压制方法的效果进行说明。图2（a）-图2（c）是本发明提供的某块地震资料异常振幅压制效果图。如图2（a）-图2（c）所示，其中，图2（a）为原始地震剖面，图2（b）为滤波后地震剖面，图2（c）为滤波去除的异常振幅，从图中我们可以看出，通过本发明提供的方法确实能最大限度的压制异常振幅信号，提高了地震资料的信噪比。而且，从图2（b），图2（c）对比可以看出，在压制异常振幅的同时，有效信号被保护得很好，这说明本发明提供的方法确实能在不伤害有效信号的情况下对异常振幅进行很好的压制。

在另一个示例中，利用本发明提供的方法对地震资料异常振幅压制效果与使用现有技术中方法的压制效果对比图如图3和图4所示，其中，图3（a）-图3（c）是利用现有技术中的方法对某地震资料异常振幅压制效果图；图4（a）-图4（c）是利用本发明提供地震资料异常振幅压制方法对与图3中相同的地震资料异常振幅压制效果图，其中，图3（a）以及图4（a）为原始地震剖面；图3（b）以及图4（b）为滤波后结果；图3（c）以及图4（c）压制掉的异常振幅。上述资料均采用四川某区块的地震资料。从上述对比中可以看到，本发明能在保护有效信号的前提下最大限度去除异常振幅，效果优于现有技术中的常规方法。

图5是本发明提供的优选的实施方式的地震资料异常振幅压制方法流程图。在一个优选的实施方式中，本发明提供的地震资料异常振幅压制方法具体的实施步骤可以包括：

第一步：对地震信号进行时频分解，在异常振幅占优的频率段进行滤波处理。

第二步：在各占优频率中，选定有效波区域，统计有效波振幅能量指数的平均值。

第三步：逐个采样点计算能量指数，与有效波振幅能量指数进行比较，若小于不处理，若大于则判定为疑似振幅值。

第四步：对于疑似振幅值，将其能量降为平均振幅能量级别，得到该采样点的平均能量指数振幅值。

第五步：对于疑似振幅值，进一步判断是否为异常振幅。判断标准为：该疑似异常振幅与均值振幅指数数量级的差值和均值振幅指数数量级之比是否超过给定的阀值，即达到足够大，才判定为异常振幅进行衰减。

第六步：衰减后的数据进行反变换，得到滤波结果，从而得到地震道信号。

通过上述步骤可以得到振幅压制后的地震道。

以上结合附图详细描述了本发明的优选实施方式，但是，本发明并不限于上述实施方式中的具体细节，在本发明的技术构思范围内，可以对本发明的技术方案进行多种简单变型，这些简单变型均属于本发明的保护范围。

另外需要说明的是，在上述具体实施方式中所描述的各个具体技术特征，在不矛盾的情况下，可以通过任何合适的方式进行组合，为了避免不必要的重复，本发明对各种可能的组合方式不再另行说明。

此外，本发明的各种不同的实施方式之间也可以进行任意组合，只要其不违背本发明的思想，其同样应当视为本发明所公开的内容。

Claims (9)

1.一种地震资料异常振幅的压制方法，其特征在于，该方法包括：

S101，提取地震资料的时域信号进行时频双域分解，以及分析分解出的时频谱以得出异常振幅占优频率；

S102，对于任意一个所述异常振幅占优频率下的信号，进行有效波振幅能量指数分析寻找有效波振幅区，并且计算有效波振幅的能量指数作为平均能量指数；

S103，在所述异常振幅占优频率下提取采样点，并计算所述采样点振幅的能量指数，当所述采样点振幅的能量指数大于所述平均能量指数时确定所述采样点具有疑似异常振幅；

S104，将所述疑似异常振幅的能量指数数量级与所述平均能量指数数量级做差再与所述平均能量指数数量级相比以得到比值，当所述比值大于给定阈值时确定所述疑似异常振幅为异常振幅，并对所述异常振幅进行衰减；

S105，对所有所述异常振幅占优频率下的时频谱信号进行S102、S103、以及S104步骤处理以得到压制异常振幅后的时频谱信号；以及

S106，对所有所述压制异常振幅后的时频谱信号进行反变换以得到压制异常振幅后的地震资料。

2.根据权利要求1所述的方法，其特征在于，所述时域信号进行时频双域分解的方法包括小波变换或S变换。

3.根据权利要求2所述的方法，其特征在于，利用下述公式进行时域信号的时频双域分解：

$W_g{S}_k(t,a)=<g_{t,a},S_k>=\underset{-\mathrm{\&infin;}}{\overset{+\mathrm{\&infin;}}{\&Integral;}}\frac{1}{a}{g}^{*}\left(\frac{\mathrm{\&tau;}-t}{a}\right)S_k\left(\mathrm{\&tau;}\right)d\mathrm{\&tau;},a\&Element;R,t\&Element;R$

S_k(t)为时域信号，W_gS_k(t,a)为时频谱，g(t)为复数母小波函数，(·)^*表示复共轭，a为频率尺度因子，τ为信号中的变量，k为分量的数量，其中，采用的复数Morlet母小波为

4.根据权利要求1所述的方法，其特征在于，计算所述有效波振幅的能量指数的方法包括：对振幅能量指数取平均值或对所述振幅能量指数取均方根。

5.根据权利要求4所述的方法，其特征在于，利用下述公式计算得出所述有效波振幅的平均能量指数：

$THD=\frac{1}{NK}\underset{k=1}{\overset{K}{\&Sigma;}}\underset{i=1}{\overset{N}{\&Sigma;}}Ceil\&lsqb;\lg \left(fabs\right(W_g{S}_k\left(i,a\right)\left)\right)\&rsqb;,$

其中，THD为平均能量指数，i为采样点，i＝1,2,3,…,N，k＝1,2,3,…,K，W_gS_k(i,a)为采样点i的时频谱，k为分量的数量，。

6.根据权利要求5所述的方法，其特征在于，计算所述采样点振幅的能量指数的方法包括对振幅能量取自然对数或对振幅能量取常用对数。

7.根据权利要求6所述的方法，其特征在于，利用下述公式计算所述采样点振幅的能量指数：

EA_ik＝Ceil(lg(fabs(W_gS_k(i,a))))，

其中，EA_ik为所述采样点振幅的能量指数，i为采样点，i＝1,2,3,…,N，W_gS_k(i,a)为采样点i的时频谱，k为分量的数量，。

8.根据权利要求7所述的方法，其特征在于，所述能量指数与有效波振幅的能量指数进行比较的方法包括：

EA_ik＞THD，则任意一个时频谱所在采样点振幅值为疑似异常振幅值；

EA_ik＜THD，则任意一个时频谱所在采样点振幅值为有效波振幅；

其中，EA_ik为所述采样点振幅的能量指数，THD为平均能量指数。

9.根据权利要求1所述的方法，其特征在于，所述确定异常振幅的步骤包括：

当所述比值小于所述给定阈值，则确定该采样点振幅值为有效波振幅值；以及

当所述比值大于所述给定阈值，则确定该采样点振幅值为异常振幅值。