CN102323618A - 基于分数阶傅里叶变换的相干噪声抑制方法 - Google Patents

Abstract

提供一种基于分数阶傅里叶变换的相干噪声抑制方法，所述方法包括：在监测的地震数据中估计相干噪声的频带和视速度范围，并提取所述频带内的地震数据；在视速度范围内，针对提取的所述频带内的地震数据来确定与视速度对应的视速度方向；基于所述视速度和对应的视速度方向来确定视倾角；根据视倾角对提取的所述频带内的地震数据执行分数阶傅里叶变换；对分数阶傅里叶变换后的地震数据执行相干预测，以得到预测的相干噪声；以及从监测的地震数据中消除预测的相干噪声。

Description

基于分数阶傅里叶变换的相干噪声抑制方法

技术领域

本发明属于石油地震勘探数据处理领域，涉及用于消除噪声干扰的技术，更具体地说，涉及一种用于利用分数阶傅里叶变换来抑制相干噪声的方法。

背景技术

利用地下介质弹性和密度的差异，通过观测和分析大地对人工激发地震波的响应，推断地下岩层的性质和形态的地球物理勘探方法称为地震勘探。地震勘探是钻探前勘测石油与天然气资源的重要手段，在煤田和工程地质勘查、区域地质研究和地壳研究等方面，也得到广泛应用。在地震勘探中，最为基础和重要的工作是如何采集到野外地震数据并对其进行有效的处理和解释。为了能够获得有效的地震数据，需要在野外设置大量的监测仪器，其中，每一个监测仪器采集到的数据称作一道地震数据。

在采集到的地震资料中，不可避免地存在着一定数量的噪声干扰。对于山地地震资料而言，其噪声来源十分复杂，噪声对有效波有较强的干扰，尤其是地震资料中普遍存在着相当严重的线性干扰，如浅层强能量的地表多次折射波和面波，中层、深层的强声波干扰，以及频带很宽的地表直达波等。一般情况下，这些干扰破坏了有效反射信号，严重时覆盖了整个地震记录，完全淹没了有效信号，大大降低了地震资料的信噪比。对这些干扰，常规的处理方法是采用f-k滤波法、τ-p变换法、中值滤波法、速度陷波法等。这些方法虽然对相干干扰的消除有一定的作用，但它们本身都存在一定的局限性。例如，f-k滤波技术要求均匀的空间采样，而陆地数据尤其3D数据常常难以满足这一苛刻条件；经典的τ-p变换也存在同样的问题；速度陷波滤波器则要求满足噪声为线状噪声的这一条件。这些方法都对不同道间的地震数据的振幅变化相当敏感，即使细微的变化也会使滤波的效果变得极差。所以如果这些相干噪声不能得到很好的抑制，则最终会影响地震记录的叠加成像质量及地震剖面的横向分辨率。

发明内容

本发明的目的在于提供一种能够基于分数阶傅里叶变换来抑制地震资料中的相干噪声的方法，从而提高在较复杂的地质环境中获得的地震数据(尤其是低信噪比的地震数据)的处理质量。

根据本发明的一方面，提供一种用于在监测的地震数据中抑制相干噪声的方法，所述方法包括：在监测的地震数据中估计相干噪声的频带和视速度范围，并提取所述频带内的地震数据；在视速度范围内，针对提取的所述频带内的地震数据来确定与视速度对应的视速度方向；基于所述视速度和对应的视速度方向来确定视倾角；根据视倾角对提取的所述频带内的地震数据执行分数阶傅里叶变换；对分数阶傅里叶变换后的地震数据执行相干预测，以得到预测的相干噪声；以及从监测的地震数据中消除预测的相干噪声。

估计相干噪声的频带和视速度范围的步骤可包括：人工估计监测的地震数据中相干噪声的频带，并估计出相干噪声的视速度范围。

提取所述频带内的地震数据的步骤可包括：在估计出相干噪声的频带之后，构建相应的带通滤波器，应用所述带通滤波器对监测的地震数据进行滤波，以得到滤波后的所述频带内的地震数据。

确定视速度方向的步骤可包括：在估计出的视速度范围内，选取一个或多个视速度来对提取的所述频带内的地震数据依次进行扫描，并由此确定与所述一个或多个视速度对应的视速度方向。

当对提取的所述频带内的地震数据依次进行扫描时，沿着各个方向对所述频带内的地震数据进行采样，从而在每个方向上得到一组采样数据，然后对每个方向上的一组采样数据相加求和，将绝对值最大的求和值所对应的方向确定为相应视速度的视速度方向。

将视速度的方向确定为视倾角的方向。

确定视倾角的步骤可包括：对于二维叠前地震数据，基于视速度和监测的地震数据中的地震道间距来确定视倾角的大小。

视倾角的大小＝(1000×地震道间距)/视速度，其中，地震道间距的单位为米，视速度的单位为米/秒。

确定视倾角的步骤包括：对于三维叠前地震数据，基于视速度和监测的地震数据中的偏移距来确定视倾角的大小。

根据视倾角对提取的所述频带内的地震数据执行分数阶傅里叶变换的步骤可包括：将提取的所述频带内的地震数据与视倾角下的分数傅里叶变换核矩阵相乘，以得到分数阶傅里叶变换后的地震数据。

对分数阶傅里叶变换后的地震数据执行相干预测的步骤可包括：采用F-X域滤波算法对分数阶傅里叶变换后的地震数据执行相干预测。

对分数阶傅里叶变换后的地震数据执行相干预测的步骤可包括：利用在所述频带外的频段求取的预测因子对分数阶傅里叶变换后的地震数据执行相干预测。

所述频带外的频段是指相干噪声较少的频段。

所述相干噪声较少的频段为中频段。

利用在所述频带外的频段求取的预测因子对分数阶傅里叶变换后的地震数据执行相干预测的步骤可包括：用求取的预测因子与所述频带内的地震数据进行褶积，从而得到相应的相干预测滤波结果。

用求取的预测因子与所述频带内的地震数据进行褶积的步骤可包括：根据中频段与相干噪声的所述频带之间的频率关系，用求取的预测因子与所述频带内的地震数据进行褶积。

从监测的地震数据中消除预测的相干噪声的步骤可包括：将预测的相干噪声进行分数阶傅里叶逆变换，以得到时空域中的相干噪声，然后从监测到的地震数据中减去时空域中的相干噪声，从而得到抑制相干噪声后的有效地震数据。

附图说明

通过下面结合附图进行的对实施例的描述，本发明的上述和/或其它目的和优点将会变得更加清楚，其中：

图1是示出根据本发明示例性实施例的相干噪声抑制方法的流程图；

图2A示出根据本发明示例性实施例监测到的地震数据的示意图；

图2B示出根据本发明示例性实施例过滤出的地震数据的示意图；

图2C示出根据本发明示例性实施例对过滤出的地震数据执行分数阶傅里叶变换的示意图；

图2D示出根据本发明示例性实施例对过滤出的地震数据执行相干预测的示意图；

图2E示出根据本发明示例性实施例抑制相干噪声后的地震数据的示图；

图3A到图3C示出根据本发明示例性实施例针对叠前二维数据来抑制相干噪声的示图；以及

图4A到图4C示出根据本发明示例性实施例针对叠前三维数据来抑制相干噪声的示图。

具体实施方式

现将详细描述本发明的示例性实施例，所述实施例的示例在附图中示出，其中，相同的标号指示相同的部分。以下将通过参照附图来说明所述实施例，以便解释本发明。

图1是示出根据本发明示例性实施例的相干噪声抑制方法的流程图。如图1所示，根据本发明示例性实施例的相干噪声抑制方法通过基于相干噪声的视速度及其方向来确定其视倾角，并按照视倾角将估计的相干噪声数据转换到分数域中，从而在分数域执行相干噪声预测。上述方法能够在较低信噪比的地震资料中有效地预测出分数域的相干噪声，在将分数域的相干噪声变换到时空域之后，从原始的地震资料中减去变换到时空域的预测相干噪声，以起到抑制相干噪声的作用。

具体说来，在步骤S100，在监测的地震数据中估计相干噪声的频带和视速度范围，并提取所述频带内的地震数据。

图2A示出根据本发明示例性实施例监测到的地震数据的示意图。如图2A所示，在由多个监测仪器监测到的具有多个地震道的地震数据中，必然存在着一定数量的相干干扰。相干干扰如地表折射多次波、面波和声波及直达波等，由于传播机理各不相同，使得它们相互之间以及与有效反射波之间存在诸多方面的差异，其显著特征之一就是频谱的衰减规律有较大不同，也就是说，不同的干扰波，其优势频带各不相同。因此，针对不同的线性干扰，在各自的优势频带范围内进行识别，不仅可以更准确地识别出线性干扰波，而且可以保证其它频段内的有效信号不受影响，使处理结果具有更高的保真度，而相干噪声一般集中于地震数据所有频率中的低频部分。因此，作为示例，可人工识别地震数据中相干噪声的频带，并同时估计出相干噪声的视速度范围。本领域技术人员应理解：本发明并不受限于任何具体的用于估计相干噪声的频带和视速度范围的方式，任何能够有助于从监测到的地震数据中分离出相干噪声并估计出相干噪声的大致视速度范围的手段均可等同地应用于本发明。在估计出相干噪声的频带之后，构建相应的带通滤波器，应用所述带通滤波器对地震数据进行滤波，以得到滤波后的地震数据。图2B示出根据本发明示例性实施例过滤出的地震数据的示意图。如图2B所示，滤波后的地震数据中主要包括相干噪声，但同时也含有一定数量的有效波。为了进一步得到其中的有效波，有必要对滤波后的地震数据进行进一步的处理。

然后，在步骤S200，在视速度范围内，针对提取的所述频带内的地震数据来确定与视速度对应的视速度方向。

作为示例，在估计出的视速度范围内，可根据实际需要选取一个或多个视速度(例如，按照预定的速度间隔来选取一个或多个典型的视速度)来对过滤出的地震数据依次进行扫描，并由此确定与所述一个或多个视速度对应的视速度方向。在利用每个视速度进行扫描时，可沿着各个方向对所述过滤出的地震数据进行采样，这样，每个方向上可得到一组采样数据，然后对每个方向上的一组采样数据相加求和，根据相干噪声的特性可知：绝对值最大的求和值所对应的方向则为相应视速度的视速度方向，也是相应视速度下的相干噪声的视速度方向。

在得到了相干噪声的视速度和视速度方向之后，可在步骤S300基于所述视速度和对应的视速度方向来确定视倾角，其中，视倾角的方向即为视速度的方向。

关于如何确定视倾角的大小，作为示例，对于二维叠前数据而言，可认为其中的相干噪声满足以下两个条件：(1)干扰波沿同相轴呈现线性关系；(2)干扰噪声的视速度不同于有效波的视速度。因此，在这种情况下，可基于道间距和视速度来确定视倾角的大小。作为优选而非限定性的计算方式，可根据下面的等式来确定视倾角的大小a：

a＝(1000×地震道间距)/视速度。

其中，地震道间距的单位为m，视速度的单位为m/s。

作为另一示例，对于三维叠前数据而言，由于在进行三维数据的野外采集时，炮点不在检波线上，因此，地震道间距是变化的，以致在较远的偏移距上，面波和直达波不呈线性变化。在这种情况下，可基于偏移距和视速度来确定视倾角的大小。

在步骤S300确定视倾角之后，本发明的示例性实施例可在步骤S400根据视倾角对提取的所述频带内的地震数据执行分数阶傅里叶变换，从而得到分数域的相应数据。图2C示出根据本发明示例性实施例对滤波后的地震数据执行分数阶傅里叶变换的示意图。如图2C所示，当根据在步骤S300确定的视倾角对在步骤S100滤波后的地震数据执行分数阶傅里叶变换(将所述地震数据与视倾角下的分数傅里叶变换核矩阵相乘以得到分数域下的对应数据)时，可看出所述变换处理相当于将地震数据在时空域中的坐标轴沿视倾角的方向旋转相应的角度a。通过分数阶傅里叶变换，可以看出分数域中的相干噪声在视倾角的方向上呈现良好的相关性，而与有效波的相关性则比较差。

因此，本发明的示例性实施例在步骤S500对分数阶傅里叶变换后的地震数据执行相干预测，以得到预测的相干噪声。根据分数阶傅里叶变换后的地震数据的特性(即，相干噪声呈现良好的相关性，而又与有效波不相关，其中，随机噪声不可预测)，可通过执行相干预测来得到相干噪声。图2D示出根据本发明示例性实施例对滤波后的地震数据执行相干预测的示意图。

以下，首先结合最简单的情况来说明相干预测的基本概念：

假设对于一组信道时差为Δt₁的相邻地震道，其中，第一地震道上的地震波数据可表示为w₁(t)，则在第一地震道上经过时间t₁的反射波为w₁(t-t₁)，对其进行傅里叶变换，可得到对于某一个频率f，则相应地构成了一个复数序列：

$W_1\left(f\right)e^{-i2\mathrm{\π}{\mathrm{ft}}_1},$ $W_1\left(f\right)e^{-i2\mathrm{\π}{\mathrm{ft}}_1}\·e^{-i2\mathrm{\πf\Δ}{t}_1},$ $W_1\left(f\right)e^{-i2\mathrm{\π}{\mathrm{ft}}_1}\·e^{-i2\mathrm{\πf}2\mathrm{\Δ}{t}_1},\·\·\·$

其中：

是对反射波w₁(t-t₁)进行傅里叶变换的输出结果，f为频率，t₁为时间，对应的离散Z变换形式为：

$S_1\left(Z\right)=\underset{n=0}{\overset{\∞}{\mathrm{\Σ}}}{W}_1\left(f\right)e^{-i2\mathrm{\π}{\mathrm{ft}}_1}\·e^{-i2\mathrm{\πfn\Δ}{t}_1}\·Z^n=\frac{W_1\left(f\right)e^{-i2\mathrm{\π}{\mathrm{ft}}_1}}{1-e^{i2\mathrm{\πf\Δ}{t}_1}Z}---\left(1\right).$

对于假设的最简单情况，可明显看出预测因子就是也就是说，前一地震道上的地震波数据乘上

后，正好是下一地震道上的地震数据。由此，可将预测出的地震数据表示为：

$W_{m+1}\left(f\right)=W_m\left(f\right)e^{-i2\mathrm{\πf\Δ}{t}_1}---\left(2\right).$

其中，m表示傅里叶变换后的第m道地震数据。

基于上述原理，实际中，可采用常规F-X域滤波算法来执行相干预测。此外，由于常规F-X域滤波算法对实际信号与理想信号之间的差异比较敏感，因此，可能会在剖面上出现蚯蚓化或细微断层模糊不清的现象。为了改善这种不稳定的预测偏差，作为优选方式，可采用改进的F-X域滤波算法，即，考虑到预测因子中频率与信道时差的对偶关系，利用相干噪声所在的频带外的预测因子来执行频带内的相干噪声预测，从而避免了如果直接利用频带内的预测因子则会受到频率低且能量较强的相干噪声影响的情况。

在选择用于获取预测因子的具体频带时，作为优选方式，可选择中频带作为用于计算预测因子的频带，这是因为：高频带信噪比较低，能量弱、求取的预测因子受噪声的影响较大，从而使得滤波去噪处理后的高频带有效信号严重畸变，不利于进一步提高分辨率。而中频带的信号中无干扰波影响，且能量比较集中，信噪比较高，同相轴的连续性很强，这个频带范围内求取的预测因子很适合F-X预测滤波的条件。

以下将说明选择中频带来计算预测因子，并利用该预测因子在滤波后的低频地震数据中执行相干预测的处理。

首先，可利用复数维纳滤波来求解预测因子OP(f，l)，其中，f为频率，l为预测步长。设原始的地震记录为s(f，x)，其中，f为频率，x为空间位置，那么对于某一频率f₀，预测误差能量E(f₀)为：

$E\left(f_0\right)=\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{[\underset{l=1}{\overset{L}{\mathrm{\Σ}}}s(f_0,x_{n-l})\·\mathrm{OP}(f_0,l)-s(f_0,x_l)]}^2---\left(3\right).$

这里，根据误差能量最小原则，本发明示例性实施例可求出与最小误差能量相应的预测误差因子OP(f₀，l)，l＝0，1，2，…，L，作为用于在过滤出的低频地震数据中执行相干预测的预测因子，其中，n表示时窗内的第n地震道，N表示时窗内的地震道数，其中l为(空间)预测步长，L表示空间预测长度。

用以上求取的预测因子OP(f₀，l)与每一地震道的地震数据进行褶积，就得到相应的相干预测滤波结果，即：

$S(f_i,x_n)=\underset{l=1}{\overset{L}{\mathrm{\Σ}}}\mathrm{OP}(f_0,l)S(f_i,x_{n-l})---\left(4\right).$

上式中，f_i表示执行相干预测的频带内的频率，l为预测步长，L为预测长度，n表示时窗内的第n地震道。

以下考虑到预测因子中频率与信道时差的对偶关系，利用相干噪声所在的频带外的预测因子来执行频带内的相干噪声预测。

首先，假设地震数据的频率分布为(f₁，f_N)，其中f₁为直流分量，f_N为Nyquist频率；在步骤S100可通过带通滤波来过滤出频率范围为(f_l1，f_l2)的低频相干噪声数据。此外，本发明示例性实施例还可过滤出信噪比较高的中频带数据，其频率范围为(f_d1，f_d2)，由于中频带(f_d1，f_d2)可较好地满足相干预测的假设条件且基本无噪声干扰，因此，在这个中频带范围内按照上述方式求取预测因子OP(f_i，x)。

然后考虑到中频带与相干噪声的低频带之间的关系(例如，假设预测因子的具体形式表示为则可将其看作频率为f、地震道间时差为Δt₁的预测因子，也可将其看作频率为f/2、道间时差为2Δt₁的预测因子)，将与最小误差能量相应的预测因子OP(f_i，x)按照频率关系应用于在步骤S100过滤出的频带(f_l1，f_l2)内的地震数据进行褶积滤波，即：

$S(f_i,x_n)=\underset{l=1}{\overset{L}{\mathrm{\Σ}}}S(f_i,x_{n-2l})\mathrm{OP}(2f_i,l),f_{l1}\≤f_i\≤f_{l2}---\left(5\right).$

由于将信噪比相对较高的中频带(f_d1，f_d2)(基本不含有相干噪声)中求取的预测因子直接应用到含有相干噪声主要频率成分的频带(f_l1，f_l2)上，从而避免了相干噪声对求取预测因子过程的不利影响，从而能有效地预测出地震数据上的相干干扰波，使得干扰波和有效信号完全分离开来。

本领域技术人员应理解：当对在步骤S100过滤出的地震数据执行相干噪声预测时，并不需要受限于上述给出的任何预测方式，即，任何现有的相干噪声预测方式及其改进方式均可等同地应用于本发明。虽然利用不同频带来计算预测因子的改进方式能够进一步消除噪声的影响，但是这种改进并不构成对本发明的限制。

然后，在步骤S600，从监测的地震数据中消除预测的相干噪声。具体说来，可将在步骤S500预测出的相干噪声进行分数阶傅里叶逆变换，以得到时空域中的相干噪声，然后从监测到的原始地震数据中减去这部分相干噪声，从而得到抑制相干噪声后的有效地震数据。图2E示出根据本发明示例性实施例抑制相干噪声后的地震数据的示图。

图3A到图3C示出根据本发明示例性实施例针对叠前二维数据来抑制相干噪声的示图，图4A到图4C示出根据本发明示例性实施例针对叠前三维数据来抑制相干噪声的示图。其中，图3A和图4A分别示出原始的实际地震资料，图3B和图4B示出根据本发明示例性实施例预测出的相干噪声，图3C和图4C示出根据本发明示例性实施例获取的有效地震信号。

由此可见，根据本发明，可有效地抑制相干噪声，并相应地得到更为精确的有效地震波信号，从而很好地反映了实际地震记录中的地层信息和地质特征。利用分数阶傅里叶变换来去除相干噪声具有较好的准确性和可靠性，从而为后续处理提供了信噪比和保真度均较高的地震数据。此外，本发明示例性实施例在分数域执行相干噪声预测时，采用了改进的滤波算法，使得求取预测因子的过程不必受噪声干扰，从而进一步提高了抑制相干噪声的可靠性。

本发明的以上各个实施例仅仅是示例性的，而本发明并不受限于此。本领域技术人员应该理解：在不脱离本发明的原理和精神的情况下，可对这些实施例进行改变，其中，本发明的范围在权利要求及其等同物中限定。

Claims (17)

1.一种用于在监测的地震数据中抑制相干噪声的方法，所述方法包括：

在监测的地震数据中估计相干噪声的频带和视速度范围，并提取所述频带内的地震数据；

在视速度范围内，针对提取的所述频带内的地震数据来确定与视速度对应的视速度方向；

基于所述视速度和对应的视速度方向来确定视倾角；

根据视倾角对提取的所述频带内的地震数据执行分数阶傅里叶变换；

对分数阶傅里叶变换后的地震数据执行相干预测，以得到预测的相干噪声；以及

从监测的地震数据中消除预测的相干噪声。

2.如权利要求1所述的方法，其中，估计相干噪声的频带和视速度范围的步骤包括：人工估计监测的地震数据中相干噪声的频带，并估计出相干噪声的视速度范围。

3.如权利要求1所述的方法，其中，提取所述频带内的地震数据的步骤包括：在估计出相干噪声的频带之后，构建相应的带通滤波器，应用所述带通滤波器对监测的地震数据进行滤波，以得到滤波后的所述频带内的地震数据。

4.如权利要求1所述的方法，其中，确定视速度方向的步骤包括：

在估计出的视速度范围内，选取一个或多个视速度来对提取的所述频带内的地震数据依次进行扫描，并由此确定与所述一个或多个视速度对应的视速度方向。

5.如权利要求4所述的方法，其中，当对提取的所述频带内的地震数据依次进行扫描时，沿着各个方向对所述频带内的地震数据进行采样，从而在每个方向上得到一组采样数据，然后对每个方向上的一组采样数据相加求和，将绝对值最大的求和值所对应的方向确定为相应视速度的视速度方向。

6.如权利要求1所述的方法，其中，将视速度的方向确定为视倾角的方向。

7.如权利要求6所述的方法，其中，确定视倾角的步骤包括：对于二维叠前地震数据，基于视速度和监测的地震数据中的地震道间距来确定视倾角的大小。

8.如权利要求7所述的方法，其中，视倾角的大小＝(1000×地震道间距)/视速度，其中，地震道间距的单位为米，视速度的单位为米/秒。

9.如权利要求6所述的方法，其中，确定视倾角的步骤包括：对于三维叠前地震数据，基于视速度和监测的地震数据中的偏移距来确定视倾角的大小。

10.如权利要求1所述的方法，其中，根据视倾角对提取的所述频带内的地震数据执行分数阶傅里叶变换的步骤包括：将提取的所述频带内的地震数据与视倾角下的分数傅里叶变换核矩阵相乘，以得到分数阶傅里叶变换后的地震数据。

11.如权利要求1所述的方法，其中，对分数阶傅里叶变换后的地震数据执行相干预测的步骤包括：采用F-X域滤波算法对分数阶傅里叶变换后的地震数据执行相干预测。

12.如权利要求1所述的方法，其中，对分数阶傅里叶变换后的地震数据执行相干预测的步骤包括：利用在所述频带外的频段求取的预测因子对分数阶傅里叶变换后的地震数据执行相干预测。

13.如权利要求12所述的方法，其中，所述频带外的频段是指相干噪声较少的频段。

14.如权利要求13所述的方法，其中，所述相干噪声较少的频段为中频段。

15.如权利要求14所述的方法，其中，利用在所述频带外的频段求取的预测因子对分数阶傅里叶变换后的地震数据执行相干预测的步骤包括：用求取的预测因子与所述频带内的地震数据进行褶积，从而得到相应的相干预测滤波结果。

16.如权利要求15所述的方法，其中，用求取的预测因子与所述频带内的地震数据进行褶积的步骤包括：根据中频段与相干噪声的所述频带之间的频率关系，用求取的预测因子与所述频带内的地震数据进行褶积。

17.如权利要求1所述的方法，其中，从监测的地震数据中消除预测的相干噪声的步骤包括：将预测的相干噪声进行分数阶傅里叶逆变换，以得到时空域中的相干噪声，然后从监测到的地震数据中减去时空域中的相干噪声，从而得到抑制相干噪声后的有效地震数据。

Images