CN102608658A - 强能量干扰抑制方法 - Google Patents

Abstract

提供了一种强能量干扰抑制方法，所述方法包括：将单炮地震数据沿地震道与时间的方向划分为多个窗口；分别对每个窗口的地震数据进行分数阶傅里叶变换，在分数频域中获得与每个窗口的地震数据对应的频域地震数据；基于获得的各个窗口的频域地震数据，为每个窗口建立强能量衰减函数；使用为每个窗口建立的强能量衰减函数，分别对每个窗口的频域地震数据进行衰减；分别对衰减后的每个窗口的频域地震数据进行分数阶傅里叶反变换，以获得去除了强能量干扰的地震数据。

Description

强能量干扰抑制方法

技术领域

本发明涉及一种地震数据处理方法，更具体地讲，涉及一种在处理地震数据时使用的强能量干扰抑制方法。

背景技术

近年来，随着油气勘探形势的发展，地震勘探的范围也越来越广泛，越来越深入，地震勘探工作也逐渐从东部转向西部，由平原转向沙漠、丘陵、山地等复杂地表区，在这些复杂的地貌地区采集的地震资料中，普遍存在着相当严重的强能量干扰，例如：声波干扰、交流电干扰、脉冲噪声及在记录深部出现的高频干扰，一般情况下这类噪声严重干扰了有效反射信号，极大地降低了地震资料的信噪比。

对于这些强能量干扰，目前常用的处理方法是滤波或道剔除，但无论是交互剔除废道还是人工编辑，都需要花费大量的时间，而且淹没在能量干扰中的有效信号也会损失，因此目前的干扰处理方法已不能满足日渐兴起的高分辨率勘探的精度要求，要解决这个问题，只有寻求自动识别和剔除这些异常能量的新方法。因此需要一种能够有效的去处采集的地震数据中包含的强能量干扰，同时提高大量地震数据处理效率的地震数据处理方法。

发明内容

根据本发明的一方面，提供了一种强能量干扰抑制方法，所述方法包括：(a)将单炮地震数据沿地震道与时间的方向划分为多个窗口；(b)分别对每个窗口的地震数据进行分数阶傅里叶变换，在分数频域中获得与每个窗口的地震数据对应的频域地震数据；(c)基于获得的各个窗口的频域地震数据，为每个窗口建立强能量衰减函数；(d)使用为每个窗口建立的强能量衰减函数，分别对每个窗口的频域地震数据进行衰减；(e)分别对衰减后的每个窗口的频域地震数据进行分数阶傅里叶反变换，以获得去除了强能量干扰的地震数据。

此外，可通过以下等式对每个窗口的地震数据进行分数阶傅里叶变换：

$f_k(x,u)=\underset{-\∞}{\overset{\∞}{\∫}}{K}_a(t,u)f_k(x,t)\mathrm{dt},$

其中，f_k(x，t)表示划分的多个窗口中的第k个窗口的地震数据，f_k(x，u)是频域中与所述第k个窗口的地震数据f_k(x，t)对应的频域地震数据，K_a(t，u)为分数阶傅里叶变换的核函数并且表示为：

$K_a(t,u)=\left\{\begin{array}{cc}\sqrt{1-i\cot (\mathrm{a\π}/2)e}*e^{\mathrm{i\π}[\cot (\mathrm{a\π}/2)-u^2-2\csc (\mathrm{a\π}/2)+\cot (\mathrm{a\π}/2)t^2]}& a\≠2n\\ \mathrm{\δ}(t-u)& a=4n\\ \mathrm{\δ}(t+u)& a=4n+2\end{array}\right.$

其中，n是整数。

此外，对于所述多个窗口中的第k个窗口的频域地震数据，建立强能量衰减函数的步骤可包括：将第k个窗口的频域地震数据f_k(x，u)划分为多个频段的频域地震数据f_k(x₁，u)、...、f_k(x_i，u)、...、f_k(x_m，u)，其中，m是划分的频段的数量，f_k(x_i，u)表示第k个窗口的第i个频段的频域地震数据，并且1≤i≤m；基于以下等式计算第k个窗口的每个频段的频域地震数据的振幅值，

A_k(x_i，u)＝|(f_k(x_i，u))|，

其中，A_k(x_i，u)表示第k个窗口的第i个频段的频域地震数据f_k(_xi，u)的振幅值；基于计算的第k个窗口的多个频段的频域地震数据的振幅值，计算第k个窗口的多个频段的频域地震数据的振幅值的加权中值M_k(u)；根据计算的第k个窗口的频域地震数据的振幅值的加权中值M_k(u)，计算第k个窗口的加权因子S_k(u)＝Thr×M_k(u)，其中，Thr是预设的阈值；基于以下等式为第k个窗口的每个频段的频域地震数据建立强能量衰减函数，

${\mathrm{Att}}_k\left(u_i\right)=\left\{\begin{array}{cc}p*(A_k(x_i,u)/M_k\left(u\right))& A_k(x_i,u)>S_k\left(u\right)\\ 1& A_k(x_i,u)\≤S_k\left(u\right)\end{array}\right.,$

其中，Att_k(u_i)表示针对第k个窗口的第i个频段的频域地震数据建立的强能量衰减函数，p是预设的衰减系数。

此外，对于所述多个窗口中的第k个窗口的频域地震数据，进行衰减的步骤可包括：根据以下等式对第k个窗口的各频段的频域地震数据进行衰减，

$f_k^{\′}(x_i,u)=f_k(x_i,u)/{\mathrm{Att}}_k\left(u_i\right)$

其中，

表示衰减后的第k个窗口的第i个频段的频域地震数据。

此外，对衰减后的每个窗口的频域地震数据进行分数阶傅里叶反变换的步骤可包括：对各窗口的衰减后的每个频段的频域地震数据进行组合，以获得衰减后的每个窗口的频域地震数据；分别对衰减后的各窗口的频域地震数据进行分数阶傅里叶反变换，以获得去除了强能量干扰的地震数据。

此外，计算第k个窗口的多个频段的频域地震数据的振幅值的加权中值的步骤可包括：为每个频段的频域地震数据的振幅值分别设置加权系数w₁、w₂、...、w_i、...、w_m，其中，w_i是为第i个频段的频域地震数据的振幅值A_k(x_i，u)设置的加权系数，并且加权系数w₁、w₂、...、w_i、...、w_m均为正整数；根据设置的加权系数，对振幅值序列中的多个振幅值按振幅值的大小进行升序/降序排序，其中，振幅值序列中的加权系数w_i表示存在w_i个A_k(x_i，u)；当

为奇数时，确定第k个窗口的振幅值的加权中值M_k(u)等于排序后的振幅值序列中的第

个振幅值，当

为偶数时，确定第k个窗口的振幅值的加权中值M_k(u)等于排序后的振幅值序列中的第

个振幅值与第

个振幅值的均值。

此外，可在可编程图形处理单元(GPU)平台下通过多线程执行本发明的强能量干扰抑制方法中的一个或多个步骤。

有益效果

本发明采用分数阶傅立叶变换来去除地震数据中的强能量干扰，能更好的分离出强能量干扰与有效信号，具有自适应性强和良好的振幅保真性的优点。另一方面，在本发明中还采用了GPU并行技术来对大量的地震数据进行并行处理，能够充分节约计算机资源，并且与现有技术相比，处理效率更高，运行速度更快，更能够满足实际生产的需求。因此，采用本发明能够在石油地震资料处理过程中发挥更一步的积极作用。

附图说明

通过下面结合附图进行的详细描述，本发明的上述和其它目的、特点和优点将会变得更加清楚，其中：

图1是示出根据本发明的实施例的应用于地震数据的强能量干扰抑制方法的流程图。

具体实施方式

以下，将参照附图更充分地描述本发明的示例性实施例，其示例性实施例在附图中示出。然而，可以以许多不同的形式实施示例性实施例，并且不应被解释为局限于在此阐述的示例性实施例。相反，提供这些实施例从而本公开将会彻底和完整，并将完全地将示例性实施例的范围传达给本领域的技术人员。

图1是示出根据本发明的实施例的应用于地震数据的强能量干扰抑制方法的流程图。

如图1中所示，在步骤101，可将单炮地震数据沿地震道与时间的方向划分为多个窗口。所述单炮地震数据是指包括多道地震数据的一个炮集的地震数据，其中，一道地震数据表示一个道上在不同时间点记录的值。可将一个炮集(例如，具有300个道且在每个道上记录2500个时间点的值的炮集)视为二维平面记录，其中，横向方向是道，纵向方向是时间点。例如，如果划分了一个32×64的窗口，则表示横向方向取了32道，而纵向方向取了64个时间点，从而随后可将该炮集按照给定横向大小(即，32道)与给定纵向大小(即，64个时间点)划分成多个窗口。在本发明中，可沿地震道与时间的方向对单炮地震数据进行任意划分，也可根据数据处理软件和硬件性能要求来对所述地震数据进行划分。

在步骤103，分别对每个窗口的地震数据进行分数阶傅里叶变换，从而获得在频域中与每个窗口的地震数据对应的频域地震数据。具体地讲，可通过下面的等式(1)对每个窗口的地震数据进行分数阶傅里叶变换：

$f_k(x,u)=\underset{-\∞}{\overset{\∞}{\∫}}{K}_a(t,u)f_k(x,t)\mathrm{dt}---\left(1\right)$

其中，f_k(x，t)表示划分的多个窗口中的第k个窗口的地震数据，f_k(x，u)是频域中与所述第k个窗口的地震数据f_k(x，t)对应的频域地震数据，K_a(t，u)为分数阶傅里叶变换的核函数并且表示为：

$K_a(t,u)=\left\{\begin{array}{cc}\sqrt{1-i\cot (\mathrm{a\π}/2)e}\×e^{\mathrm{i\π}[\cot (\mathrm{a\π}/2)-u^2-2\csc (\mathrm{a\π}/2)+\cot (\mathrm{a\π}/2)t^2]}& a\≠2n\\ \mathrm{\δ}(t-u)& a=4n\\ \mathrm{\δ}(t+u)& a=4n+2\end{array}\right.---\left(2\right)$

其中，n是整数，a是由用户设置的分数傅里叶变换的阶数并且在划分的每个窗口的分数阶傅里叶变换中具有相同的值。此外，当a取-a时，即得到分数阶傅氏逆变换，当a＝1时，为普通傅里叶正变换。当a＝-1时，为普通傅里叶反变换。

然后在步骤105，可基于获得的各个窗口的频域地震数据，为每个窗口建立强能量衰减函数。以下将以第k个窗口的频域地震数据为例详细说明为每个窗口建立强能量衰减函数的具体方法。

对于第k个窗口的频域地震数据，可首先将所述第k个窗口的频域地震数据f_k(x，u)划分为多个频段(例如，1～30Hz、30～60Hz、60～90Hz......)的频域地震数据f_k(x₁，u)、...、f_k(x_i，u)、...、f_k(x_m，u)，其中，m是划分的频段的数量，f_k(x_i，u)表示第k个窗口的第i个频段的频域地震数据，1≤i≤m。在对每个窗口的频域地震数据的频段划分中采用相同的频段划分方式。然后基于下面的等式(3)计算划分的每个频段的频域地震数据的振幅值：

A_k(x_i，u)＝|(f_k(x_i，u))|                 (3)，

其中，A_k(x_i，u)表示第k个窗口的第i个频段的频域地震数据f_k(x_i，u)的振幅值。

在获得第k个窗口的每个频段的频域地震数据的振幅值之后，可基于获得的振幅值来计算第k个窗口的多个频段的频域地震数据的振幅值的加权中值M_k(u)。在本发明中，可通过排序的方法来获得上述振幅值的加权中值M_k(u)。例如，对于第k个窗口的多个频段的频域地震数据，可首先为每个频段的振幅值分别设置加权系数w₁、w₂、...、w_i、...、w_m，其中，w_i是为第i个频段的频域地震数据的振幅值A_k(x_i，u)设置的加权系数，并且加权系数w₁、w₂、...、w_i、...、w_m均为正整数。然后，根据设置的加权系数，对振幅值序列

中的多个振幅值按振幅值的大小进行升序/降序排序，其中，振幅值序列中的加权系数w_i表示存在w_i个A_k(x_i，u)。在对上述振幅值序列中的多个振幅值完成排序之后，可从排序后的振幅值序列中获得加权中值M_k(u)。具体地讲，当

为奇数时，可确定第k个窗口的振幅值的加权中值M_k(u)等于排序后的振幅值序列中的第

个振幅值；当

为偶数时，可确定第k个窗口的振幅值的加权中值M_k(u)等于排序后的振幅值序列中的第

个振幅值与第个振幅值的均值。

应该理解，上述加权系数w₁、w₂、...、w_i、...、w_m的设置可根据强能量抑制的需求而设置，或者由用户根据经验而设置。此外，也可将加权系数全部设置为1，从而简化加权中值的计算。

在获得第k个窗口的多个频段的频域地震数据的振幅值的加权中值M_k(u)之后，可根据计算的加权中值M_k(u)来计算第k个窗口的加权因子S_k(u)＝Thr×M_k(u)，其中，Thr是预设的阈值，并可根据需要而由用户任意确定。然后，可基于以下等式(4)为第k个窗口的每个频段的频域地震数据建立强能量衰减函数：

${\mathrm{Att}}_k\left(u_i\right)=\left\{\begin{array}{cc}p*(A_k(x_i,u)/M_k\left(u\right))& A_k(x_i,u)>S_k\left(u\right)\\ 1& A_k(x_i,u)\≤S_k\left(u\right)\end{array}\right.---\left(4\right),$

其中，Att_k(u_i)表示针对第k个窗口的第i个频段的频域地震数据建立的强能量衰减函数，p是预设的衰减系数，并可根据用户的需求而设置。

在步骤107，使用为每个窗口建立的强能量衰减函数，分别对每个窗口的频域地震数据进行衰减。具体地讲，对于第k个窗口的频域地震数据，可根据以下等式(5)对每个频段的频域地震数据进行衰减：

$f_k^{\′}(x_i,u)=f_k(x_i,u)/{\mathrm{Att}}_k\left(u_i\right)---\left(5\right)$

其中，

表示衰减后的第k个窗口的第i个频段的频域地震数据。

然后，在步骤109，分别对衰减后的每个窗口的频域地震数据进行分数阶傅里叶反变换，以获得去除了强能量干扰的地震数据。具体地讲，可在频域中对各窗口的衰减后的每个频段的频域地震数据进行组合，以获得衰减后的每个窗口的频域地震数据，然后分别对衰减后的各窗口的频域地震数据进行分数阶傅里叶反变换，从而获得去除了强能量干扰的地震数据。继续以第k个窗口的频域地震数据为例说明，在完成步骤107之后，可通过对衰减后的每个频段的频域地震数据

进行组合来获得衰减后的第k个窗口的全频段的频域地震数据

然后，对获得的第k个窗口的频域地震数据

进行分数阶傅里叶反变换，以在时域中获得衰减后的第k个窗口的地震数据F_k(x，t)。在本发明中，可通过下面的等式(6)来对每个窗口的频域地震数据进行分数阶傅里叶反变换：

$F_k(x,t)=\underset{-\∞}{\overset{\∞}{\∫}}{K}_{-a}(t,u)f_k^{\′}(x,u)\mathrm{du}---\left(6\right)$

其中，K_-a(t，u)为分数阶傅里叶反变换的核函数且可表示为：

$K_{-a}(t,u)=\left\{\begin{array}{cc}\sqrt{1-i\cot (-\mathrm{a\π}/2)e}\×e^{\mathrm{i\π}[\cot (-\mathrm{a\π}/2)-u^2-2\csc (-\mathrm{a\π}/2)+\cot (-\mathrm{a\π}/2)t^2]}& a\≠2n\\ \mathrm{\δ}(t-u)& a=4n\\ \mathrm{\δ}(t+u)& a=4n+2\end{array}\right.---\left(2\right)$

其中，n是整数。

通过以上方法，由于通过加权因子以及衰减系数对各窗口的地震数据的每个频段的振幅值进行了不同程度的衰减，因此可对地震数据中的能量干扰部分进行衰减，从而有效地抑制地震数据中的强能量干扰。

在本发明中，考虑需要对大量的地震数据同时进行例如分数阶傅里叶变换/反变换等数据处理，因此，可在可编程图形处理单元(GPU，GraphicProcessing Unit)平台下使用多线程的方法来进行处理，这是因为GPU是一种高度并行化、多线程、多核的处理器，具有杰出的计算功率和极高的存储器带宽，能够适应大规模数据运算的需要。因此在大量地震数据的情况下，使用GPU能够有效地提高数据处理的效率。

在本发明中，可通过GPU平台来执行用于实现根据本发明的强能量抑制方法的程序。首先，可在GPU平台下对需要去除强能量干扰的地震数据进行初始化(即，在GPU上申请显存的临时空间，将数据从内存复制到显存)，减少临时空间的使用(例如，将临时空间进行合并或者重复使用同一块临时空间)，最大化并行度(例如，可根据本发明的强能量抑制方法中的窗口的划分来选取并行度，诸如可由每个线程读取一个窗口的数据，或者每个线程读取一个窗口数据)以获取最高的性能提升。然后，可在GPU平台下通过多线程的方式来完成图1中示出的方法的一个或多个步骤：例如，可在GPU平台下通过多线程同时进行多个窗口的地震数据的分数阶傅里叶变换，以获得所述多个窗口的频域地震数据；可在GPU平台下通过多线程同时进行多个窗口中的多个频段的频域地震数据的振幅值的计算；可在GPU平台下使用多线程通过例如双调排序算法对振幅值序列进行并行排序，从而求出多个频段的频域地震数据的振幅值的加权中值；可在GPU平台下通过多线程对多个窗口的多个频段的频域地震数据同时进行衰减；可在GPU平台下通过多线程来对衰减后的多个窗口的频域地震数据同时进行分数阶傅立叶反变换等。

应该理解，在图1中示出的方法中的一个或多个步骤均可在GPU平台下通过多线程的方式来实现，或者可由用户根据软硬件性能、处理程序等因素来确定哪些步骤需要通过GPU的多线程来完成，从而更快更好的对大量的地震数据进行有效处理。

以上方法可更好的分离出强能量与有效信号，并有自适应性强和良好的振幅保真性的优点。另一方面，在本发明中还采用了GPU并行技术来对大量的地震数据进行并行处理，能够充分节约计算机资源，并且与现有技术相比，处理效率更高、运行速度更快，更能够满足实际生产的需求。因此，采用本发明能够在石油地震资料处理过程中发挥更一步的积极作用

尽管已经参照其示例性实施例具体显示和描述了本发明，但是本领域的技术人员应该理解，在不脱离权利要求所限定的本发明的精神和范围的情况下，可以对其进行形式和细节上的各种改变。

Claims (7)

1.一种强能量干扰抑制方法，包括：

(a)将单炮地震数据沿地震道与时间的方向划分为多个窗口；

(b)分别对每个窗口的地震数据进行分数阶傅里叶变换，在分数频域中获得与每个窗口的地震数据对应的频域地震数据；

(c)基于获得的各个窗口的频域地震数据，为每个窗口建立强能量衰减函数；

(d)使用为每个窗口建立的强能量衰减函数，分别对每个窗口的频域地震数据进行衰减；

(e)分别对衰减后的每个窗口的频域地震数据进行分数阶傅里叶反变换，以获得去除了强能量干扰的地震数据。

2.如权利要求1所述的方法，其中，通过以下等式对每个窗口的地震数据进行分数阶傅里叶变换：

$f_k(x,u)=\underset{-\∞}{\overset{\∞}{\∫}}{K}_a(t,u)f_k(x,t)\mathrm{dt}$

其中，f_k(x，t)表示划分的多个窗口中的第k个窗口的地震数据，f_k(x，u)是频域中与所述第k个窗口的地震数据f_k(x，t)对应的频域地震数据，K_a(t，u)为分数阶傅里叶变换的核函数并且表示为：

$K_a(t,u)=\left\{\begin{array}{cc}\sqrt{1-i\cot (\mathrm{a\π}/2)e}\×e^{\mathrm{i\π}[\cot (\mathrm{a\π}/2)-u^2-2\csc (\mathrm{a\π}/2)+\cot (\mathrm{a\π}/2)t^2]}& a\≠2n\\ \mathrm{\δ}(t-u)& a=4n\\ \mathrm{\δ}(t+u)& a=4n+2\end{array}\right.$

其中，n是整数。

3.如权利要求2所述的方法，其中，对于所述多个窗口中的第k个窗口的频域地震数据，建立强能量衰减函数的步骤包括：

将第k个窗口的频域地震数据f_k(x，u)划分为多个频段的频域地震数据f_k(x₁，u)、...、f_k(x_i，u)、...、f_k(x_m，u)，其中，m是划分的频段的数量，f_k(x_i，u)表示第k个窗口的第i个频段的频域地震数据，并且1≤i≤m；

基于以下等式计算第k个窗口的每个频段的频域地震数据的振幅值，

A_k(x_i，u)＝|(f_k(x_i，u))|，

其中，A_k(x_i，u)表示第k个窗口的第i个频段的频域地震数据f_k(x_i，u)的振幅值；

基于计算的第k个窗口的多个频段的频域地震数据的振幅值，计算第k个窗口的多个频段的频域地震数据的振幅值的加权中值M_k(u)；

根据计算的第k个窗口的频域地震数据的振幅值的加权中值M_k(u)，计算第k个窗口的加权因子S_k(u)＝Thr×M_k(u)，其中，Thr是预设的阈值；

基于以下等式为第k个窗口的每个频段的频域地震数据建立强能量衰减函数，

${\mathrm{Att}}_k\left(u_i\right)=\left\{\begin{array}{cc}p*(A_k(x_i,u)/M_k\left(u\right))& A_k(x_i,u)>S_k\left(u\right)\\ 1& A_k(x_i,u)\≤S_k\left(u\right)\end{array}\right.,$

其中，Att_k(u_i)表示针对第k个窗口的第i个频段的频域地震数据建立的强能量衰减函数，p是预设的衰减系数。

4.如权利要求3所述的方法，其中，对于所述多个窗口中的第k个窗口的频域地震数据，进行衰减的步骤包括：

根据以下等式对第k个窗口的各频段的频域地震数据进行衰减，

$f_k^{\′}(x_i,u)=f_k(x_i,u)/{\mathrm{Att}}_k\left(u_i\right)$

其中，

表示衰减后的第k个窗口的第i个频段的频域地震数据。

5.如权利要求4所述的方法，其中，步骤(e)包括：

对各窗口的衰减后的每个频段的频域地震数据进行组合，以获得衰减后的每个窗口的频域地震数据；

分别对衰减后的各窗口的频域地震数据进行分数阶傅里叶反变换，以获得去除了强能量干扰的地震数据。

6.如权利要求3所述的方法，其中，计算第k个窗口的多个频段的频域地震数据的振幅值的加权中值的步骤包括：

为每个频段的频域地震数据的振幅值分别设置加权系数w₁、w₂、...、w_i、...、w_m，其中，w_i是为第i个频段的频域地震数据的振幅值A_k(x_i，u)设置的加权系数，并且加权系数w₁、w₂、...、w_i、...、w_m均为正整数；

根据设置的加权系数，对振幅值序列

中的多个振幅值按振幅值的大小进行升序/降序排序，其中，振幅值序列中的加权系数w_i表示存在w_i个A_k(x_i，u)；

当为奇数时，确定第k个窗口的振幅值的加权中值M_k(u)等于排序后的振幅值序列中的第

个振幅值，

当

为偶数时，确定第k个窗口的振幅值的加权中值M_k(u)等于排序后的振幅值序列中的第

个振幅值与第

个振幅值的均值。

7.如权利要求1所述的方法，其中，在可编程图形处理单元平台下通过多线程执行步骤(a)到(e)中的一个或多个步骤。