CN110865410A - 一种基于nar-tfpf压制地震勘探随机噪声的方法 - Google Patents

Abstract

本发明属于地震勘探技术领域，具体涉及一种基于NAR‑TFPF压制地震勘探随机噪声的方法。本发明基于NAR‑TFPF压制地震勘探随机噪声的方法包括以下步骤：对地震含噪信号进行分段；根据地震资料中有效信号和随机噪声非线性特性有很大的差异，通过非线性自回归模型(NAR)区分信号段和噪声段；采用时频峰值滤波(TFPF)进行去噪处理，噪声段选取大窗长进行噪声压制，信号段选取小窗长进行保幅，本发明可以在压制强随机噪声的同时保留有效信号，提高地震资料的信噪比，能更准确地用于地质解释。

Description

一种基于NAR-TFPF压制地震勘探随机噪声的方法

技术领域

本发明属于地震勘探技术领域，具体涉及一种地震勘探随机噪声压制的方法，尤其涉及一种基于NAR-TFPF压制地震勘探随机噪声的方法。

背景技术

地震勘探随机噪声始终伴随着有效信号出现在地震资料中，能量较大的随机干扰会直接影响动、静校正分析以及最终的成像效果，给后续的地质解释带来不利影响。为此，专家学者们设计、改进了很多用于处理随机噪声的算法。其设计思路大致上可以分为信号增强和噪声压制两类，目的是尽可能地提高地震资料的信噪比，例如f-x反褶积，多项式拟合等属于信号增强算法，中值约束下的矢量分解算法属于噪声压制方法。这些算法已经被广泛使用，且其应用效果相对能够满足生产要求。与此同时，这些算法在应用过程中受到一些假设条件的限制，直接影响了滤波效果。譬如这些算法在处理随机噪声时，以有效信号空间相关而随机噪声空间不相关为依据，且随机噪声是时空域二维的随机过程，只研究空间方向上的随机噪声是远远不够的。

Boashash和Mesbah于2004年提出了一年种新的信号增强算法——时频峰值滤波。基于时频分析理论，先将含噪信号调制解析信号，求出解析信号的维格纳时频分布(WVD)，通过估计WVD的最大值获取瞬时频率，达到信号增强的目的。TFPF最先应用于新生儿的脑电分析中，2005年被首次应用到地震勘探随机噪声的压制中，根据TFPF的滤波理论，它能够从强高斯白噪声中恢复线性信号，但会对非线性信号造成失真，而地震信号尤其是频率较高的信号都具有非线性特性。为了避免这一问题，用伪维格纳分布(PWVD)来替代WVD，以此来确保窗内瞬时频率的线性条件。改进后的TFPF对强随机噪声背景下地震信号的恢复取得了有效的结果，但窗长的选取对于TFPF压制随机噪声和保留有效信号波形至关重要。为了获得更好的滤波效果，在过去的几年里，林红波等人利用窗长，采样频率和信号主频之间的关系，提出了变窗长TFPF。该方法先将含噪信号分成任意长度的时间序列，在不同时间序列选择不同的窗长，以此来解决由于TFPF窗长不合适导致的高频分量失真。刘彦萍等人采用经验模式分解(EMD)来改进传统的TFPF。利用EMD的分解特性将含噪信号分解为从高频到低频的几种模态，并为不同的频率模态选择不同的TFPF窗长，通过识别信号分量来实现降噪以及保幅。对含噪信号直接进行时间段分段处理，无法准确区分噪声段和信号段，且受时间序列长度影响，不能采取最合适的窗长；对含噪信号进行EMD分解，容易出现模态混叠的问题。

近年来，非线性是随机信号分析中的重要标准之一，并被广泛应用于临床诊断，股票预测，机械故障诊断等领域。地震信号与随机噪声的非线性程度不同，通过非线性度区分地震数据中的信号部分和噪声部分，并在噪声部分选择长窗长，信号部分选择短窗长，同时增强信噪比和信号保幅度。当一道含噪信号使用固定窗长进行滤波时，选择长窗长能够较好地抑制噪声，但同时有效信号幅度也被衰减；选择短窗长对有效信号幅值有较好的保幅度，但不能更有效地压制噪声，滤波后的地震数据信噪比较低。如何能在压制噪声的同时保持有效信号的幅度，是我们亟待解决的问题。

发明内容

针对现有技术中存在的无法同时有效压制噪声和保持信号幅度的问题，本发明提供了一种基于NAR-TFPF压制地震勘探随机噪声的方法。

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

一种基于NAR-TFPF压制地震勘探随机噪声的方法，包括以下步骤：

步骤1，对地震数据进行PCNN分段处理；

步骤2，对每段时间序列进行NAR建模，计算NAR的最佳阶数估计其非线性度；

步骤3，为每段时间序列选择合适的窗长进行TFPF滤波，完成处理。

进一步地，所述步骤1对地震数据进行PCNN分段处理具体操作是：选择脉冲耦合神经网络(PCNN)对地震数据进行分段，将时间序列任意分割成几个片段。TFPF滤波窗长的选择是影响滤波效果的一个关键参数。使用同一窗长处理信号部分和噪声部分是不合理的，因此，我们可以将含噪地震信号分成几个时间段。由于常用的时间序列分割方法，如定长分割，快速分割，模糊分割等方法可能会将一个周期或一个完整的序列小波分成两部分，其结果受到序列的长度和振幅的影响，而PCNN用于时间序列分割时可以避免这些问题。如图5所示，图5(a)为只有一个子波的含噪信号，图5(b)为定长分割结果，将一个完整子波分成两部分，图5(c)为PCNN分割结果，避免了图5(b)中的问题。

再进一步地，所述步骤2中的进行NAR建模，计算NAR的最佳阶数估计其非线性度具体操作是：

y_n＝G[y_n-1,y_n-2,…,y_n-k] 公式(1)

其中y是系统响应，k是记忆因子，n是离散时间，G是非线性方程；

当G是混沌映射时，系统输出是混沌序列；当时间序列{y_n}，n＝1,2，…N是已知的，则它可以由方程G建模，得多项式G，多项式G包括所有用来获得通用非线性方程的递增阶数；

y_n是历史状态{y_n-1+y_n-2+…+y_n-k}0～d次幂的总和，其表达式为

其中

表示预测的系统响应，k是记忆因子，d是阶数，TL是项数，a_m，m＝0,1,2…k,k+1,k+2…TL-1是通过任意曲线拟合方法获得的系数；

公式2也可以为

确定记忆因子k和阶数d的值，根据公式(2)建立模型后按升幂排列所有项；截取建立模型的前r项得到一个多项式，系数a_m通过最小二乘法拟合求出，然后计算预测的时间序列

通过公式(2)得到预测值

预测值

和公式(1)中所示原始值{y_n}，n＝k+1，k+2…N.之间的均方误差为：

其中ε表示预测误差，

是均值并且满足

计算预测误差ε，改变r的取值并重复公式(2)～公式(3)的计算过程，绘制ε(r)-r曲线；

当预测误差ε最小时d的取值为最佳阶数，可记为d_optM，最佳阶数越大，表示该时间序列的非线性程度越高；

所述NAR模型用于对分段后的时间序列进行建模。图1为一道高斯白噪声，对其进行NAR建模得到的ε(r)-r曲线如图2所示，可以看出当项数r＝28时预测误差最小，此时d＝3，即图1所示的噪声序列最佳阶数d_optM＝3。图3所示是含有一个主频为30Hz地震子波的时间序列，对其进行NAR建模，ε(r)-r曲线如图4所示。可以看出，当r＝6时预测误差最小，此时d＝2，即图3所示的子波序列最佳阶数d_optM＝2。噪声序列的非线性程度高于有效信号，因此可以通过对其进行NAR建模来区分一道含噪信号中的信号部分和噪声部分。

更进一步地，所述步骤3中为每段时间序列选择合适的窗长进行TFPF滤波，完成处理的具体操作是：通常情况下，含噪信号表达式为：

其中s(t)为含噪信号，x_g(t)为有效信号，n_g(t)为加性噪声，x_kg(t)为有效信号的分量；

将含噪信号s(t)看作调频信号的瞬时频率，通过对s(t)调频得到解析信号，如下式所示：

其中z_s(t)为解析信号，exp是以自然常数e为底的指数函数，μ为尺度参数，j为纯虚数，π为圆周率，λ为积分变量,s(λ)为积分变量λ的含噪信号；

通过估计解析信号z_s(t)的WVD峰值可以将有效信号从背景噪声中恢复出来，

其中

为滤波结果，W_Zs(t,f)为解析信号z_s(t)的WVD时频分布，arg max为求W_Zs(t,f)最大值的函数；

根据TFPF滤波理论，当有效信号随时间线性变化，且背景噪声是高斯白噪声时，通过公式(4)～公式(6)得到滤波后的信号是含噪信号s(t)的无偏估计，因此基于WVD的TFPF滤波方法关键是保证有效信号的线性特性；而地震信号尤其是主频较高的地震信号一般是非线性的，因而采用加窗WVD分布即PWVD代替WVD，从而保证窗内信号的近似线性特性；

伪维格纳分布PWVD定义式如下：

其中h(ττ)为窗函数，Z_s为公式(5)中的解析信号，

为Z_s的共轭，j为纯虚数，f为频率变量；

用伪维格纳分布PWVD的最大值代替WVD的最大值；

适合地震子波的TFPF窗函数取值为：

其中f_s是信号采样频率，f_d是地震子波的主频，WL是窗长且只能取奇数；窗长只能选奇数且窗长太长会造成信号畸变。TFPF滤波窗长的选择是影响滤波效果的一个关键参数。使用同一窗长处理信号部分和噪声部分是不合理的。

非线性程度高为噪声段，选择较长的TFPF窗长压制噪声；

非线性程度低为信号段，选择较短的TFPF窗长对信号保幅。

一、噪声的非线性程度比信号高，因此根据非线性程度来判断噪声段和信号段；二、窗长越大，时间序列幅度衰减的越厉害，噪声段就选大窗长，噪声衰减的严重，信号段就选小窗长，信号幅度衰减小，从而提高信噪比。

与现有技术相比本发明具有以下优点：

针对地震资料信噪比极低，TFPF窗长不可变时不能在有效压制噪声的同时保持信号幅度，本发明提出通过NAR建模寻找TFPF最佳滤波窗长，提高地震资料的信噪比和分辨率。整个数据处理过程中，根据NAR对每段数据进行建模，从而识别和检测有效信号，不需要对算法有条件假设；每段采用不同的窗长进行TFPF滤波，实现TFPF窗长自适应滤波，不会在压制噪声的同时消减有效信号，能够从强背景噪声中有效地恢复出有效信号，增大了地震数据的使用性，为后续地质解释提供更多有利的信息。

附图说明

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

图1是地震随机噪声时间序列；

图2是地震随机噪声序列的lnε-r曲线；

图3地震有效信号雷克子波；

图4是地震有效信号雷克子波的lnε-r曲线；

图5是单道地震含噪记录及其分段结果示意图，其中，(a)表示含噪时间序列，(b)表示固定时间分段结果，(c)表示PCNN分段结果；

图6是本发明实施案例1，将其应用到一道地震含噪信号处理中，其中(a)为PCNN分段后的含噪序列，(b)表示TFPF去噪结果，(c)表示NAR-TFPF去噪结果；

图7是本发明实施案例2，地震合成记录处理图，(a)纯净记录，(b)含噪记录，(c)TFPF噪声压制后的结果，(d)NAR-TFPF去噪结果；

图8是从图7所示的各记录的频域-波数谱，(a)为纯净记录(b)为含噪记录，(c)为TFPF噪声压制后的结果，(d)为NAR-TFPF去噪结果；

图9是部分野外采集数据处理图，(a)为含噪数据，(b)TFPF噪声压制后的结果，(c)为NAR-TFPF去噪结果；

图10是图9中方框中所示部分的放大图，(a)为含噪数据，(b)TFPF噪声压制后的结果，(c)为NAR-TFPF去噪结果。

具体实施方式

为了使本技术领域的人员更好地理解本发明的技术方案，下面将结合本实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述。

实施例1

将本发明应用于一道含有一个地震子波的含噪信号，采样频率为1000Hz，子波频率为20Hz，噪声为高斯白噪声，信噪比为-7dB，如图6(a)所示。

将该含噪信号经传统TFPF滤波处理，由公式(8)，滤波窗长WL≤19.2，窗长只能选奇数且窗长太长会造成信号畸变，这里选最佳滤波窗长为15，滤波结果如图6(b)所示。

将含噪信号经NAR-TFPF处理，首先，将含噪序列经PCNN分段，如图6(a)所示。其次，对每一段时间序列进行NAR建模，由公式(2)-(3)，求得图6(a)中各段的最佳阶数分别为d_optM1＝3,d_optM2＝3,d_optM3＝2,d_optM4＝3,d_optM5＝3，第三段最佳阶数最小，即非线性程度最低，因此第三段为信号段。

第三，分别对每一段时间序列进行TFPF滤波，第1、2、4、5段为噪声段，选大一些的窗长来压制噪声，这里WL＝19；第3段为信号段，选短一些的窗长保持信号幅度，这里选WL＝11，滤波结果如图6(c)所示。可以看出，NAR-TFPF能够有效地压制背景噪声并较好地保持信号幅度。

实施例2

本发明利用一组40道并含有三条相交同相轴的人工合成的地震记录进行试验，地震子波主频分别为25Hz，40Hz和50Hz，采样频率为1000Hz，背景噪声为高斯白噪声，信噪比为-7dB，如图7所示。图7(a)为纯净记录，图7(b)为含噪记录。分别将图7(b)中的每一道含噪记录进行TFPF滤波，子波主频最高频率为50Hz，由公式(8)，最佳滤波窗长为7，整个记录经TFPF处理后如图7(c)所示。

对图7(b)中的每一道记录分别进行NAR-TFPF处理。首先对每一道进行PCNN分段处理。其次，分别对每一个时间序列片段进行NAR建模，由公式(2)-(3)，求得每一段时间序列的最佳阶数，判断其非线性特性，由此找出信号段和噪声段。最后，分别对每段采取合适的滤波窗长，噪声段采用大窗长WL＝17，信号段采用小窗长WL＝5，处理后的记录如图7(d)所示。可见以传统TFPF滤波处理，即使选择了最佳窗长，压制了大部分的随机噪声，同时有效同相轴也被部分衰减，导致这些同相轴的连续性不够好。NAR-TFPF能够抑制绝大部分的随机噪声，所有的反射同相轴比图5(c)中的更加清晰。

对图7(a)～(d)中所示的记录进行二维傅里叶变换，得到其频域-波数(f-k)谱，如图8所示。图8(a)为纯净记录的f-k谱，图8(b)为含噪信号的f-k谱，图8(c)是TFPF滤波后的f-k谱，图8(d)是NAR-TFPF处理后的f-k谱。可以看出，NAR-TFPF的处理结果更接近纯净记录。另外，与TFPF相比，NAR-TFPF能够抑制低频(0-10Hz)随机噪声，如图8中的黑圈内所示。

实施例3

将本发明应用于实际野外地区的地震勘探背景噪声的压制处理，截取某地区野外采集数据如图9(a)所示，可以看出共炮点记录中存在大量的噪声，信噪比较低。

使用TFPF对原始记录进行滤波，实际地震记录不能确定有效信号的主频，根据文献记载，一般地震勘探子波的主频范围大约为60Hz以内，由公式(8)，滤波窗长最小为7，在最大程度地压制噪声并保证有效信号不失真的情况下，不断调整滤波窗长，确定最佳滤波窗长为11，滤波结果如图9(b)所示。在图9(b)中，可以看到，TFPF可以抑制大量随机噪声，但不能很好地恢复同相轴。

对原始信号进行NAR-TFPF进行处理，首先对图9(a)中每一道进行PCNN分段处理。其次，分别对每一个时间序列片段进行NAR建模，由公式(,2)-(3)，求得每一段时间序列的最佳阶数，判断其非线性特性，由此找出信号段和噪声段。最后，分别对每段采取合适的滤波窗长，噪声段采用大窗长WL＝15，信号段采用小窗长WL＝7，处理后的记录如图9(c)所示。可以看出，NAR-TFPF除了压制大量的随机噪声，同时更多的反射同相轴被恢复出来。为了更清楚有效地对比滤波想过，将图9中框内的部分放大如图10所示。可以看出，图10(c)中的同相轴比图10(b)中的更清晰连续，说明通过NAR建模改善TFPF滤波窗长的选择，可以提高地震数据的信噪比和分辨率。

模拟实验和实际数据处理结果均可看出，本发明提出的NAR-TFPF压制地震勘探随机噪声的方法，能够实现TFPF的可变窗长，寻找TFPF最佳窗长进行滤波。除了可以抑制绝大部分的随机噪声，NAR-TFPF在分辨反射同相轴，保持信号幅度和同相轴的的连续性方面都有更好的处理结果，使地震数据更可靠，方便作进一步的地质解释。

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

Claims (4)

1.一种基于NAR-TFPF压制地震勘探随机噪声的方法，其特征在于：包括以下步骤：

步骤1，对地震数据进行PCNN分段处理；

步骤2，进行NAR建模，计算NAR的最佳阶数估计其非线性度；

步骤3，为每段时间序列选择合适的窗长进行TFPF滤波，完成处理。

2.根据权利要求1所述一种基于NAR-TFPF压制地震勘探随机噪声的方法，其特征在于：所述步骤1对地震数据进行PCNN分段处理具体操作是：选择脉冲耦合神经网络PCNN对地震数据进行分段，将时间序列任意分割成几个片段。

3.根据权利要求1所述一种基于NAR-TFPF压制地震勘探随机噪声的方法，其特征在于：所述步骤2中的进行NAR建模，计算NAR的最佳阶数估计其非线性度具体操作是：

y_n＝G[y_n-1,y_n-2,…,y_n-k] 公式(1)

其中y是系统响应，k是记忆因子，n是离散时间，G是非线性方程；

当G是混沌映射时，系统输出是混沌序列；当时间序列{y_n}，n＝1,2，…N是已知的，则它可以由方程G建模，得多项式G，多项式G包括所有用来获得通用非线性方程的递增阶数；

y_n是历史状态{y_n-1+y_n-2+…+y_n-k}0～d次幂的总和，其表达式为

其中

表示预测的系统响应，k是记忆因子，d是阶数，TL是项数，a_m，m＝0,1,2…k,k+1,k+2…TL-1是通过任意曲线拟合方法获得的系数；

确定记忆因子k和阶数d的值，根据公式(2)建立模型后按升幂排列所有项；截取建立模型的前r项得到一个多项式，系数a_m通过最小二乘法拟合求出，然后计算预测的时间序列

通过公式(2)得到预测值

预测值

和原始值{y_n}，n＝k+1，k+2…N.之间的均方误差为：

其中ε表示预测误差，

是均值并且满足

计算预测误差ε，改变r的取值并重复公式(2)～公式(3)的计算过程，绘制ε(r)-r曲线；

当预测误差ε最小时d的取值为最佳阶数，可记为d_optM，最佳阶数越大，表示该时间序列的非线性程度越高，为噪声段；最佳阶数越小，表示该时间序列的非线性程度越低，为信号段；

所述NAR模型用于对时间序列进行建模。

4.根据权利要求1所述一种基于NAR-TFPF压制地震勘探随机噪声的方法，其特征在于：所述步骤3中为每段时间序列选择合适的窗长进行TFPF滤波，完成处理的具体操作是：通常情况下，含噪信号表达式为：

其中s(t)为含噪信号，x_g(t)为有效信号，n_g(t)为加性噪声，x_kg(t)为有效信号的分量；

将含噪信号s(t)看作调频信号的瞬时频率，通过对s(t)调频得到解析信号，如下式所示：

其中z_s(t)为解析信号，exp是以自然常数e为底的指数函数，μ为尺度参数，j为纯虚数，π为圆周率，λ为积分变量,s(λ)为积分变量λ的含噪信号；

通过估计解析信号z_s(t)的WVD峰值可以将有效信号从背景噪声中恢复出来，

其中

为滤波结果，W_Zs(t,f)为解析信号z_s(t)的WVD时频分布，arg max为求W_Zs(t,f)最大值的函数；

伪维格纳分布PWVD定义式如下：

其中h(ττ)为窗函数，z_s为公式(5)中的解析信号，

为z_s的共轭，j为纯虚数，f为频率变量；

用伪维格纳分布PWVD的最大值代替WVD的最大值；

适合地震子波的TFPF窗函数取值为：

其中f_s是信号采样频率，f_d是地震子波的主频，WL是窗长且只能取奇数；

非线性程度高为噪声段，选择较长的TFPF窗长压制噪声；

非线性程度低为信号段，选择较短的TFPF窗长对信号保幅。

Images