CN110941015A

CN110941015A - 基于频率切片时频峰值滤波压制地震勘探随机噪声的方法

Info

Publication number: CN110941015A
Application number: CN201910367701.1A
Authority: CN
Inventors: 李光辉; 冯志强; 弓子卉; 梁美彦
Original assignee: Shanxi University
Current assignee: Shanxi University
Priority date: 2019-05-05
Filing date: 2019-05-05
Publication date: 2020-03-31
Anticipated expiration: 2039-05-05
Also published as: CN110941015B

Abstract

本发明涉及一种基于频率切片时频峰值滤波压制地震勘探随机噪声的方法，包括以下步骤：将时间域的含噪信号变换为频域空间的含噪信号，并对频域空间的含噪信号进行切片处理，将其分成若干段、任意长度的频率切片，然后确定每一段频率切片的主频和滤波窗长，之后分别对每一段频率切片进行时频峰值滤波，最后将滤波后的时间域信号分量进行叠加即得到滤波后信号。本发明将频率切片和时频峰值滤波TFPF结合起来，目的主要是通过切片处理能够有效帮助TFPF针对不同主频的信号选择最佳滤波窗长，解决了TFPF不能使用自适应窗长滤波的问题。本发明更能有效压制强随机噪声，有较高的信噪比和信号保幅效果，为进一步做出地质解释提供可靠依据。

Description

基于频率切片时频峰值滤波压制地震勘探随机噪声的方法

技术领域

本发明属于地震勘探技术领域，具体涉及一种压制地震勘探随机噪声的方法，尤其涉及一种基于频率切片时频峰值滤波压制地震勘探随机噪声的方法。

背景技术

地震勘探是油气勘探的主要手段，其主要过程包括三个环节：数据采集、数据处理和地质解释。地下构造和岩性有关的信息都包含在野外采集的地震数据中，有效信息与背景噪声叠加，被外界干扰因素所扭曲，被截断甚至直接被湮没，给地质解释带来了极大的干扰。根据噪声在地震剖面中出现的特征将其分为规则噪声和不规则噪声，规则噪声在时间上的出现有规律性，其波形具有一定的频率和视速度，相对比较容易压制；不规则噪声即随机噪声没有固定的频率，也没有特定的传播方向，在地震资料中形成杂乱无章的背景，噪声压制难度较大。

目前已有的滤波方法例如小波变换(Gao et al.2006)、中值滤波(Liu etal.2011)、Curvelet(Lari et al.2014)、多项式拟合(Yu et al.1989；Liu et al.2011)、奇异值分解(Kreimer et al.2012)、f-x反褶积(Aali et al.2005)等技术在压制随机噪声时常常不能达到令人满意的效果。时频峰值滤波(TFPF)是2004年由Boashash和Mesbah提出的，随后被应用到地震勘探随机噪声的压制当中。TFPF能够在极低信噪比下恢复出有效信号并具有良好的保幅性，取决于其具有能够对非平稳信号进行估计而不需要任何先验条件下这一优势。根据TFPF的滤波原理，窗长是其压制随机噪声恢复有效信号的一个关键因素，长窗长虽能有效平滑噪声，但是有效信号的幅度也被严重削弱，且窗长越长，高频信号的畸变也越严重；短窗长虽能提高有效信号的保幅度，但不能有效削减噪声，提高信噪比。因此，合理地选取滤波窗长是TFPF的一个关键问题。TFPF窗长在地震信号处理中的选取受信号主频和采样频率的影响，地震记录是不同主频信号在不同位置的叠加。目前，采用单一窗长进行滤波，信噪比和保幅度都不能得到有效提高。

此外，针对TFPF窗长的选取问题，林红波等人利用窗长、采样频率和信号主频之间的关系提出了最适合地震子波的TFPF窗长计算经验公式(Lin et al.2007；Li etal.2009)；刘彦萍等人(2014)提出经验模态分解时频峰值滤波方法，先利用经验模态(EMD)分解对含噪信号进行分解，对每一个分量选择合适的窗长对其进行滤波，然后滤波后的分量进行叠加得到去噪后的信号，这一方法取得了良好的滤波效果，但是存在EMD分解的缺陷——模态混叠的问题。

发明内容

鉴于此，本发明的目的在于，提供一种基于频率切片时频峰值滤波压制地震勘探随机噪声的方法，旨在克服现有技术中存在的采用单一窗长滤波不能有效提高信噪比和保幅度的缺陷以及解决滤波过程中存在频率混叠的问题。

为了解决上述技术问题，本发明采用的方案如下：

一种基于频率切片时频峰值滤波压制地震勘探随机噪声的方法，包括以下步骤：

基于频率切片时频峰值滤波压制地震勘探随机噪声的方法，其特征在于，包括以下步骤：

S1：将时间域的含噪信号变换为频域空间的含噪信号

对时间域的含噪信号进行傅里叶变换，得到频域空间的含噪信号，变换公式如下：

其中，F(ω)为频域空间的含噪信号，s(t)为时间域的含噪信号，ω是角频率，j是虚数单位；

S2：对频域空间的含噪信号进行切片处理

将步骤S1中的频域空间的含噪信号沿着频率方向切成若干段、任意长度的频率切片，切片采用的公式如下：

其中，n为切片的数目，F_i(ω)为每一段频率切片；

S3：确定每一段频率切片的主频

将步骤S2中的每一段频率切片的峰值频率近似看作是该段频率切片的主频；

S4：确定每一段频率切片的滤波窗长

每一段频率切片的滤波窗长记作WL，窗长WL的确定函数如下：

其中，f_h为每一段频率切片的主频，f_s是采样频率，σ²为随机噪声方差，C是特征参数，并且窗长WL只能取奇数；

S5：分别对步骤S2中的每一段频率切片进行傅里叶逆变换后进行时频峰值滤波

S5-1：分别对步骤S2中的每一段频率切片进行傅里叶逆变换，得到相应含噪信号的时间域分量，公式如下：

其中，s_i(t)为每一段频率切片的傅里叶逆变换，也是含噪信号s(t)的时间域分量；

S5-2：将每一个含噪信号的时间域分量s_i(t)表示成有效信号与随机噪声的和，采用公式如下：

其中，x_g(t)是有效信号，x_kg(t)是有效信分量，N_s是分量个数，n_g(t)是加性随机噪声；

S5-3：将步骤S5-2中的时间域分量s_i(t)进行频率编码，变换为具有瞬时频率的解析信号z_s(t)，采用公式如下：

其中，μ为尺度参数，λ为第一积分变量，j为虚数单位；

S5-4：对步骤S5-3中的解析信号运行伪维格纳分布运算，伪维格纳分布

表示为：

其中，f表示频率，*表示共轭操作，ττ表示第二积分变量，h(ττ)为窗函数且窗函数h(ττ)记作窗长WL；

S5-5：对步骤S5-4中已运行伪维格纳分布运算的解析信号计算该解析信号的时频分布峰值频率作为有效信号的无偏估计，即获得滤波后的时间域信号分量，表达式如下：

其中，

即为滤波后的时间域信号分量，μ为尺度参数；

S6：将滤波后的时间域信号分量进行叠加即得到滤波后信号，公式如下：

式中，x(t)为滤波后的信号。

本发明的有益效果：针对使用单一窗长的时频峰值滤波无法同时有效削减噪声和保持信号幅度的问题，本申请提出将频率切片和TFPF结合起来来压制地震勘探中的随机噪声，本申请对每一道单道的地震记录进行频率切片，将其分成一系列不同频率范围的数据片段，根据每一频率片段的主频选择合适的滤波窗长进行TFPF滤波，将所有滤波后的时域数据片段叠加后得到滤波后的信号，从而实现了不同频率片段使用不同窗长滤波，达到了对地震记录进行自适应变窗长的TFPF滤波的目的，通过此种滤波方法能够有效压制随机噪声的基础上保留有效信号，使反射同相轴清晰连续，同时提高了地震数据的信噪比和保幅性。

此外，将含噪信号在频域范围内进行切片处理，不受含噪信号的时长、频率、统计特性等影响，可以任意选择特定的频率范围，通过计算信号主频选择合适的窗长进行有针对性地滤波，而数据的其余部分不受影响，因此，在滤波过程中不存在频率混叠问题。

附图说明

构成本申请的一部分的附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：

图1是本发明基于频率切片时频峰值滤波压制地震勘探随机噪声的方法中所涉及的频率切片流程图；

图2是本发明实施例1中所涉及的一道含有三个地震子波的含噪信号与有效信号的示意图，其中，(a)表示纯净信号，(b)表示含噪信号；

图3是图2所示的含噪信号图2(b)经傅里叶变换后得到含噪信号的频谱图；

图4是图3所示的频域空间的含噪信号进行切片处理后的频谱切片图；

图5是图2所示的含噪信号进行滤波后的地震信号图，其中，(a)表示噪声段的放大部分，(b)表示信号峰值的放大部分；

图6是地震合成记录处理图，(a)纯净记录，(b)含噪记录，(c)TFPF处理结果，(d)FS-TFPF处理结果；

图7是地震记录频率-波数谱图，(a)纯净记录，(b)含噪记录，(c)TFPF处理结果，(d)FS-TFPF处理结果；

图8是部分野外采集数据处理图，(a)含噪数据，(b)TFPF处理结果，(c)FS-TFPF处理结果。

具体实施方式

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

根据TFPF的滤波理论，其滤波窗长是影响滤波效果的直接因素。长窗长能够有效压制噪声同时也削弱信号幅值，甚至会引起信号波形畸变，尤其是高频信号；短窗长信号保幅度好但不能有效去除噪声。由于地震记录是各种主频信号和噪声的综合记录，因此在滤波时采用单一窗长并不合适。

基于此，本发明在处理地震勘探随机噪声时，针对TFPF的最佳滤波窗长取决于信号的主频率这一关键参数对滤波效果的影响，考虑将含噪信号在频域范围内进行切片处理，可以任意选择特定的频率范围，通过计算信号主频选择合适的窗长进行有针对性地滤波，而数据的其余部分不受影响。

本发明通过频率切片能够有效帮助TFPF针对不同主频的信号选择最佳滤波窗长，解决了TFPF不能使用自适应窗长滤波的问题。

下面将结合附图说明来具体解释本发明的技术方案。

本发明技术方案的具体实施方式为：一种基于频率切片时频峰值滤波压制地震勘探随机噪声的方法，包括以下步骤：

S1：将时间域的含噪信号变换为频域空间的含噪信号；

S2：对频域空间的含噪信号进行切片处理；

S3：确定每一段频率切片的主频；

S4：确定每一段频率切片的滤波窗长；

S5：分别对步骤S2中的每一段频率切片进行傅里叶逆变换后进行时频峰值滤波和

S6：将滤波后的时间域信号分量进行叠加即得到滤波后信号。

接下来，具体阐述压制地震勘探随机噪声的方法中的每个步骤。

在步骤S1中，将时间域的含噪信号变换为频域空间的含噪信号，具体的变换操作为：

其中，F(ω)为频域空间的含噪信号，s(t)为时间域的含噪信号，ω是角频率，j是虚数单位。

同时，为了更好的理解上述变换，可结合图1中A-B的流程示意。在图1的A中，通过横测线、纵测线和时间建立起来的立体坐标轴形象的表示了时间域的输入数据，将A中的时间域的输入数据通过傅里叶变换变换成频率空间的数据，即B中由横测线、纵测线和频率建立起来的立体坐标轴。在该B中立体坐标轴体现的频率空间的数据分为实部与虚部，正好与上述的公式(1)对应。

在上述说明中，涉及的频率空间的数据与频域空间的含噪信号为同一意思表示。

在步骤S2中，对频域空间的含噪信号进行切片处理。具体操作为：参见图1中的C图，将步骤S1中的频域空间的含噪信号沿着频率方向切成若干段、任意长度的频率切片，切片采用的公式如下：

其中，n为切片的数目，F_i(ω)为每一段频率切片。

在步骤S3中，确定每一段频率切片的主频。这里的每一段频率切片的峰值频率被近似看作是该段频率切片的主频。

这是由于，对于地震子波来说，不同的子波主频与峰值频率之间的关系不同，例如Ricker子波的峰值频率约为主频的0.77倍；正弦指数衰减子波的峰值频率约为主频的0.9倍；带通Zinc子波不存在峰值频率，其主频取决于高截频或上限频率，且为高截频的0.7倍。考虑到实际地震资料不能确定地震子波类型，因此，可以将地震记录的峰值频率近似看作是地震子波的主频，尽管不够准确和严格，但是作为保守估计基本上是合理的。

在步骤S4中，确定每一段频率切片的滤波窗长。这里每一段频率切片的滤波窗长记作WL，窗长WL的确定函数如下：

其中，f_h为每一段频率切片的主频，f_s是采样频率，σ²为随机噪声方差，C是特征参数，由信号和噪声的统计特性来描述，并且窗长WL只能取奇数。

根据TFPF的滤波理论显示窗长是TFPF的一个关键参数，直接影响了TFPF的滤波效果。若将含噪信号由相同的窗长进行平滑，窗长较小，能够得到较好的有效信号保幅效果，但同时噪声压制能力弱，信噪比不能得到有效提高；窗长较大，噪声压制能力强，同时有效信号的幅度衰减严重，甚至会造成波形畸变。因此，先将含噪信号进行频率切片，由步骤S4中，确定每一个片段的滤波窗长，噪声段采取大窗长压制，信号段采取小窗长保幅，极大地提高地震资料的信噪比和分辨率。

在步骤S5中，分别对步骤S2中的每一段频率切片进行傅里叶逆变换后进行时频峰值滤波，该步骤与图1中的D图相匹配。这一步骤具体可分为5个分步骤。即：

其中，μ为尺度参数，λ为第一积分变量，j为虚数单位；

S5-4：对步骤S5-3中的解析信号运行伪维格纳分布PWVD运算，伪维格纳分布

表示为：

其中，

即为滤波后的时间域信号分量，μ为尺度参数；

在本发明的步骤S5-4中，依据传统TFPF理论，在时频域内，若解析信号的瞬时频率为关于时间的线性函数，解析信号的维格纳分布WVD能量主要集中在瞬时频率，解析信号的时频峰值即其瞬时频率的无偏估计。实际地震反射信号是关于时间的非线性函数，通过WVD得到的无偏估计存在误差，尤其是地震子波的波峰波谷处其非线性特性更强。因此，加窗WVD即伪维格纳分布PWVD计算解析信号的时频分布

在WVD的基础上加一个窗函数进行计算，在窗长范围内，解析信号的瞬时频率可近似为线性的，计算解析信号PWVD的峰值，从而获得有效信号的无偏估计。

在步骤S6中，将滤波后的时间域信号分量进行叠加即得到滤波后信号，公式如下：

式中，x(t)为滤波后的信号，

为滤波后的时间域信号分量。

本发明的滤波方法针对使用单一窗长的时频峰值滤波无法有效削减噪声的缺陷，提出结合频率切片和TFPF的方法，该方法通过对每一道单道的地震记录分成一系列不同频率范围的数据片段，根据每一频率片段的主频选择合适的滤波窗长进行TFPF滤波，从而解决了TFPF不能使用自适应窗长滤波的问题，从而实现了反射同相轴的清晰连续，地震数据较高信噪比的目的，具有实用性和有效性，为进一步做出地质解释提供可靠依据。这里涉及的数据片段与频率切片为同一意思表示。

下面将本发明涉及的基于频率切片时频峰值滤波压制地震勘探随机噪声的方法应用于一道含有三个地震子波的含噪信号和对一组40道并含有三条同相轴的人工合成的地震记录进行试验以及应用于实际野外地区的地震勘探随机噪声的压制处理。

实施例1

将本发明应用于一道含有三个地震子波的含噪信号，信噪比为-6dB，采样频率为1000Hz，子波频率分别为20Hz，40Hz和60Hz。此时的含噪信号如图2(b)所示，有效信号即纯净信号如图2(a)所示。可见在图2(b)中，有效信号几乎被噪声湮没。

将时间域含噪信号图2(b)经傅里叶变换后得到含噪信号的频谱图，如图3所示，可以看出信号能量集中在150Hz以内。然后针对图3所表示的频域空间的含噪信号进行切片处理，具体为在频率范围内沿着频率方向进行任意长度的切片，以形成若干个切片片段即频率切片，如图4所示。随后计算每一个切片片段的峰值频率以确定每一段频率切片的主频。

最后确定每一个切片片段的滤波窗长。在此基础上，分别对每一个切片片段进行傅里叶反变换后分别进行TFPF滤波，最后将所有滤波后的信号分量叠加。滤波后的信号如图5所示，在图5中还分别给出了噪声段和信号段的放大部分，具体可参见图5(a)和图5(b)。由此可以看出由频率切片-时频峰值滤波FS-TFPF处理后的信号噪声幅值小于时频峰值滤波TFPF滤波结果，而有效信号幅值大于TFPF结果。

实施例2

本发明利用一组40道并含有三条同相轴的人工合成的地震记录进行试验，信噪比为-6dB，地震子波主频分别为20Hz，40Hz和60Hz，如图6所示，图6(a)为纯净记录，图6(b)为含噪记录，图6(c)为TFPF滤波后的记录，图6(d)为FS-TFPF滤波后的记录。

可以看出几乎被噪声湮没的纯净记录经过TFPF处理后，除去了大部分噪声，但信号幅度尤其是高频信号受到了一定的衰减，出现了同相轴集中性和连续性相对较差的问题。图6(d)显示的经过FS-TFPF处理后的记录，可以看出不同的频率段采用不同的滤波窗长，在噪声得到有效压制的情况下，也能保留有效信号的幅度，同相轴清晰连续。

图7显示了图6中各个地震合成记录的频率-波数谱。图7(a)～(d)分别为为纯净记录、含噪记录、TFPF处理结果和FS-TFPF处理结果的频率-波数谱。可以看出，TFPF能够压制大部分高频噪声，但是对低频噪声的压制和有效信号的保幅都不是特别理想，而FS-TFPF不仅能压制绝大部分高频噪声，低频部分也能有效压制，如图7(a)和(d)中蓝色方框内所示，有效信号的能量也更加集中。

实施例3

将本发明应用于实际野外地区的地震勘探随机噪声的压制处理，截取某地区野外采集数据如图8(a)所示，可以看出共炮点记录中存在大量的随机噪声，信噪比较低。图(b)和图(c)分别为TFPF和FS-TFPF处理后的数据，通过比较可以看出，经TFPF处理后的数据，虽能去除大部分噪声，信噪比得到了较大的提高，但是基于单一窗长滤波后的记录中，有的反射同相轴并未能清晰连续地显现出来；由FS-TFPF滤波后的数据，可以基于不同频率段的数据选取适合的窗长，在有效压制噪声的同时，反射同相轴也清晰连续地显现出来。

模拟实验和实际数据处理结果均可看出，相较于TFPF，本发明提出的FS-TFPF更能在有效压制随机噪声的基础上保留有效信号，使反射同相轴清晰连续，同时提高地震数据的信噪比和保幅性。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。