去除地震信号中的随机噪声的方法和装置与流程

本公开涉及地震信号处理领域，更具体地，涉及一种去除地震信号中的随机噪声的方法和一种去除地震信号中的随机噪声的装置。

背景技术：

目前有多种压制地震信号中的随机噪声的处理技术，大致可分为时空域和变换域两类方法。时空域随机噪音衰减方法主要包括叠加处理、褶积滤波、多项式拟合技术、中值滤波、奇异值分解、混沌振子检测技术等。变换域随机噪声衰减方法主要包括LK-变换、XF-域滤波、τ-p变换、时频峰值滤波、小波变换去噪技术、稳定反Q滤波谱矩阵技术、变分数维方法、模型约束方法、复数道分析、零相位谱增强技术、非物理可实现空间预测滤波等技术、小波变换等。

现有的这些方法都难以在保留边缘、纹理等高频信息的同时有效地去除其中的随机噪声。

技术实现要素：

本公开提出了一种方法，其可在较好地保留边缘、纹理等高频信息的同时有效地去除其中包含的随机噪声。本公开还提出了相应的装置。

根据本公开的一方面，提出了一种去除地震信号中的随机噪声的方法，包括：对地震信号进行多尺度Contourlet变换，得到在多尺度多方向的Contourlet系数，Contourlet系数cl,k(m)表示为：其中，l表示尺度编号，k表示方向编号，0≤k<2^l-1；m、n表示在尺度l的二维离散希尔伯特空间L²(Z²)的位置信息，Z²表示二维离散信号的整数网格，L²表示平方可积空间；bl(n)表示在尺度l上进行DFB分解前的高频子带信号；Sl,k表示在尺度l、方向k的采样系数，hl,k(Sl,km-n)表示在尺度l的DFB分解中在方向k的楔形滤波器系数，是二维离散希尔伯特空间L²(Z²)中的一组基；得到第一子带系数P1(cl,k,m)：

其中H表示预设的硬阈值；基于第一子带系数P1(cl,k,m)重构所述地震信号的高频子带信号。

根据本公开的另一方面，还提出了一种去除地震信号中的随机噪声的装置，包括：Contourlet变换单元，用于对地震信号进行多尺度Contourlet变换，得到在多尺度多方向的Contourlet系数，Contourlet系数cl,k(m)表示为：其中，l表示尺度编号，k表示方向编号，0≤k<2^l-1；m、n表示在尺度l的二维离散希尔伯特空间L²(Z²)的位置信息，Z²表示二维离散信号的整数网格，L²表示平方可积空间；bl(n)表示在尺度l上进行DFB分解前的高频子带信号；Sl,k表示在尺度l、方向k的采样系数，hl,k(Sl,km-n)表示在尺度l的DFB分解中在方向k的楔形滤波器系数，是二维离散希尔伯特空间L²(Z²)中的一组基；硬判断单元，用于得到第一子带系数P1(cl,k,m)：

其中H表示预设的硬阈值；信号重构单元，用于基于第一子带系数P1(cl,k,m)重构所述地震信号的高频子带信号。

本公开的各方面通过采用Contourlet变换分离出地震信号中的Contourlet系数(也称为高频子带系数)，然后利用地震信号中有效信号与随机噪声信号的Contourlet系数的差异，有效地去除高频信息中的随机噪声。本公开尤其适用于去除随机噪声。

附图说明

通过结合附图对本公开示例性实施方式进行更详细的描述，本公开的上述以及其它目的、特征和优势将变得更加明显，其中，在本公开示例性实施方式中，相同的参考标号通常代表相同部件。

图1示出了Contourlet变换的分解流程图。

图2示出了Contourlet合成(即Contourlet变换的逆过程)的流程图。

图3示出了根据本公开的实施例的去除地震信号中的随机噪声的方法的流程示意图。

图4示出了原始地震信号和在尺度4(也称为4级分解)分解得到的Contourlet系数的波形示意图。

图5(a)示出了一个被加入随机噪声的合成地震记录。

图5(b)示出了应用本公开去除图5(a)中的随机噪声后的示意图。

图5(c)示出了应用本公开所去除的随机噪声的示意图。

图6(a)示出了一个实际地震记录。

图6(b)示出了应用本公开去除图6(a)中的随机噪声后的示意图。

图6(c)示出了应用本公开所去除的随机噪声的示意图。

具体实施方式

下面将参照附图更详细地描述本公开的优选实施方式。虽然附图中显示了本公开的优选实施方式，然而应该理解，可以以各种形式实现本公开而不应被这里阐述的实施方式所限制。相反，提供这些实施方式是为了使本公开更加透彻和完整，并且能够将本公开的范围完整地传达给本领域的技术人员。

此处，先对本公开的基本原理进行简单介绍。

Contourlet变换是一种用于图像处理领域的二维表示方法，具有多分辨率、局部定位、多方向性、近邻界采样和各向异性等性质，其基函数分布于多尺度、多方向上，少量系数即可有效地捕捉图像中的边缘轮廓，而边缘轮廓正是自然图像中的主要特征。Contourlet变换的基本思想是首先用一个类似小波的多尺度分解捕捉边缘奇异点，再根据方向信息将位置相近的奇异点汇集成轮廓段。

图1示出了Contourlet变换的分解流程图。在每一尺度，先通过拉普拉斯塔式滤波器(LP)将原始图像分解成低通逼近图像和差值图像。该低通逼近图像可被称为低频子带信号，该差值图像可被称为高频子带信号。该高频子带信号被输入二维方向滤波器(DFB)，分解得到在各个方向的一系列Contourlet系数，其可将奇异点连成线形结构，从而捕获图像中的轮廓。LP和DFB结合形成的双层滤波器组结构，被称为塔形方向滤波器组(PDFB)。

而在每一尺度分解得到的低频子带信号则可乘以一采样矩阵，然后可被送入下一级的PDFB。依此类推，可进行多尺度LP分解，在每一尺度，可由DFB将经LP分解得到的低频子带信号进一步分解为在多方向的一系列Contourlet系数。

图2示出了Contourlet合成(即Contourlet变换的逆过程)的流程图。

发明人经过深入研究发现，地震信号在经Contourlet变换后，其中表示边缘、纹理等信息的高频有效信号的能量积聚在有限的Contourlet系数区域内，而随机噪声的能量分布范围很广、幅值很小，特别是高斯白随机噪声，其变换后能量在所有Contourlet系数上均匀分布，其幅值接近于零。因此，可考虑基于Contourlet系数的绝对值的大小，来提取有效信号。

实施例1

图3示出了根据本公开的实施例的去除地震信号中的随机噪声的方法的流程示意图，该方法包括：

S1，对地震信号进行多尺度Contourlet变换，得到在多尺度多方向的Contourlet系数，Contourlet系数cl,k(m)表示为：

其中，l表示尺度编号，k表示方向编号，0≤k<2^l-1；m、n表示在尺度l的二维离散希尔伯特空间L²(Z²)的位置信息，Z²表示二维离散信号的整数网格，L²表示平方可积空间；bl(n)表示在尺度l上进行DFB分解前的高频子带信号；Sl,k表示在尺度l、方向k的采样系数，hl,k(Sl,km-n)表示在尺度l的DFB分解中在方向k的楔形滤波器系数，是二维离散希尔伯特空间L²(Z²)中的一组基。

S2，得到第一子带系数P1(cl,k,m)：

其中H表示预设的硬阈值。

在本公开的一些实施例中，所有尺度的硬阈值都相同。本领域技术人员可根据经验来设置该硬阈值。

S3，基于第一子带系数P1(cl,k,m)重构上述地震信号的高频子带信号。

例如，重构的高频子带信号可表示为：

其中是的对偶基。

通过上述方法，可在重构地震信号中的边缘、纹理等信息的同时去除随机噪声(特别是高斯白随机噪声)。

在本公开所公开的上述方法中，硬阈值是否恰当会影响去噪效果。若硬阈值过大，可能会损伤有效信号；若硬阈值过小，则可能会产生随机噪声残留。

发明人进一步研究发现，地震信号经Contourlet变换后，在各个方向中，有效信号的Contourlet系数仍保持很强的相关性，而随机噪声的Contourlet系数则不具备该相关性，例如可参见图4。图4示出了原始地震信号和在尺度4(也称为4级分解)分解得到的Contourlet系数的波形示意图。其中左侧的一幅大图示出了原始的地震信号，右侧的16幅小图示出了经4级Contourlet变换后分解得到的全部16个方向的Contourlet系数。从图4中明显可以看出，在各个方向其有效信号都得以保留。因此，在根据本公开的一些实施例中，可在重构高频子带信号前基于上述特性对第一子带系数P1(cl,k,m)进行进一步处理，以降低对硬阈值的要求，这可获得更好的去噪效果。下文公开了本公开的一些实施例中对第一子带系数P1(cl,k,m)进行进一步处理的示意性方法。

S301，在每一尺度的每个方向，可以以互相关值为判断标准对各个点对应的第一子带系数P1(cl,k,m)进行倾角扫描，以及可以基于每个点在各个倾角方向上的互相关值判断该点是否是强有效信号。

根据Contourlet变换的理论可知，每一尺度可以包括多个方向，每一方向可以包括多个道，每一道可以包括多个点，每个点可以对应于一个Contourlet系数cl,k(m)。倾角扫描是地震信号处理领域常见的一种数据处理方法，可采用本领域技术人员已知的任意手段来执行该倾角扫描。而针对每个点，可以得到其在各个倾角方向上的互相关值中的最大值和第二大值，如果该最大值与该第二大值间的差大于互相关阈值，则可判断该点是强有效信号。该互相关阈值可以是根据经验设置的。

S302，在每一尺度的每个方向，可以基于在该尺度该方向上的强有效信号对应的第一子带系数P1(cl,k,m)确定在该尺度该方向的软阈值。

发明人研究发现，在每一尺度的每个方向上，可以认为任意强有效信号的能量都大于任意随机噪声的能量，因此，可以考虑基于在一定尺度的一定方向上的各个强有效信号对应的各个第一子带系数P1(cl,k,m)的绝对值中的最小值确定在该尺度该方向的软阈值。发明人经过深入研究进一步发现，随机噪声的幅值至少比强有效信号的幅值小一个数量级，因此，优选地，可将一定尺度的一定方向上的软阈值设置为该尺度该方向上的各个强有效信号对应的各个第一子带系数P1(cl,k,m)的绝对值中的最小值的10％，既可有效地滤除随机噪声，又可充分保留所需的轮廓、纹理等信息。

S303，可以得到第二子带系数P2(cl,k,m)：

其中Fl,k表示在尺度l、方向k的软阈值。

S304，可以基于第二子带系数P2(cl,k,m)重构所述地震信号的高频子带信号。

例如，重构的高频子带信号可表示为：

其中是的对偶基。

应用上述方法，可降低对硬阈值的设置要求。例如，实施过程中，可以将硬阈值设置地稍微小一点，随后可通过自适应软阈值进一步去除在硬阈值处理过程中没有去除干净的随机噪声。同时，先通过硬阈值处理过程，既可直接去除部分干扰，也可减少后续软阈值处理过程中的计算量。这些实施例很好地结合了硬阈值判断和自适应软阈值判断的优点，并补偿了其不足，有利于获得更好的去噪效果。

实施例2

本公开还提供了一种去除地震信号中的随机噪声的装置，包括Contourlet变换单元、硬判断单元和信号重构单元。

Contourlet变换单元，用于对地震信号进行多尺度Contourlet变换，得到在多尺度多方向的Contourlet系数，Contourlet系数cl,k(m)表示为：

其中，l表示尺度编号，k表示方向编号，0≤k<2^l-1；m、n表示在尺度l的二维离散希尔伯特空间L²(Z²)的位置信息，Z²表示二维离散信号的整数网格，L²表示平方可积空间；bl(n)表示在尺度l上进行DFB分解前的高频子带信号；Sl,k表示在尺度l、方向k的采样系数，hl,k(Sl,km-n)表示在尺度l的DFB分解中在方向k的楔形滤波器系数，是二维离散希尔伯特空间L²(Z²)中的一组基。

硬判断单元，用于得到第一子带系数P1(cl,k,m)：

其中表示预设的硬阈值。

信号重构单元，用于基于第一子带系数P1(cl,k,m)重构所述地震信号的高频子带信号。

基于第一子带系数P1(cl,k,m)重构所述地震信号的高频子带信号还可以包括：

每一尺度的每一方向可以包括多个道，每一道可以包括多个点，每个点可以对应于一个Contourlet系数cl,k(m)，相应地，每个点可以对应于一个第一子带系数P1(cl,k,m)，在每一尺度的每个方向，可以以互相关值为判断标准对各个点对应的第一子带系数P1(cl,k,m)进行倾角扫描，以及可以基于每个点在各个倾角方向上的互相关值判断该点是否是强有效信号；

在每一尺度的每个方向，可以基于在该尺度该方向上的强有效信号对应的第一子带系数P1(cl,k,m)确定在该尺度该方向的软阈值；

可以得到第二子带系数P2(cl,k,m)：

其中Fl,k可以表示在尺度l、方向k的软阈值；

可以基于第二子带系数P2(cl,k,m)重构所述地震信号的高频子带信号。

例如，针对每个点，可以得到其在各个倾角方向上的互相关值中的最大值和第二大值，如果该最大值与该第二大值间的差大于互相关阈值，则可以判断该点是强有效信号。本领域技术人员可根据经验来设置硬阈值。

例如，可以基于在一定尺度的一定方向上的各个强有效信号对应的各个第一子带系数P1(cl,k,m)的绝对值中的最小值确定在该尺度该方向的软阈值，例如，软阈值可以是该尺度该方向上的各个强有效信号对应的各个第一子带系数P1(cl,k,m)的绝对值中的最小值的10％。

应用示例

为便于理解本公开实施例的方案及其效果，以下给出两个具体应用示例。本领域技术人员应理解，该示例仅为了便于理解本公开，其任何具体细节并非意在以任何方式限制本公开。图5和图6中各幅图的横坐标均表示道号，纵坐标表示时间。

图5(a)示出了一个被加入随机噪声的合成地震记录。图5(b)示出了应用本公开去除图5(a)中的随机噪声后的示意图。图5(c)示出了应用本公开所去除的随机噪声的示意图。从图5(a)中可以看出，该合成地震记录中存在很强的随机噪声，资料的信噪比很低。如图5(b)所示，应用本公开去除随机噪声后，资料的信噪比得到了明显提高。图5(c)中未见有效信号。可见，本公开能有效抑制地震信号中的随机噪声，且保幅保真效果好。

图6(a)示出了一个实际地震记录。图6(b)示出了应用本公开去除图6(a)中的随机噪声后的示意图。图6(c)示出了应用本公开所去除的随机噪声的示意图。从图6(a)中可以看出，该地震资料受吸收衰减的影响，远偏移距资料能量较弱，压制干扰波能力较差，随机噪声严重，信噪比很低。应用本公开去除随机噪声后，如图6(b)所示，单炮记录的整体面貌有了较大改善，信噪比大幅提高，远偏移距低信噪比区域的有效反射信号得到增强(如图中箭头指向部位所示)。图6(c)中未见明显有效信号。

本公开可以是系统、方法和/或计算机程序产品。计算机程序产品可以包括计算机可读存储介质，其上载有用于使处理器实现本公开的各个方面的计算机可读程序指令。

计算机可读存储介质可以是可以保持和存储由指令执行设备使用的指令的有形设备。计算机可读存储介质例如可以是――但不限于――电存储设备、磁存储设备、光存储设备、电磁存储设备、半导体存储设备或者上述的任意合适的组合。计算机可读存储介质的更具体的例子(非穷举的列表)包括：便携式计算机盘、硬盘、随机存取存储器(RAM)、只读存储器(ROM)、可擦式可编程只读存储器(EPROM或闪存)、静态随机存取存储器(SRAM)、便携式压缩盘只读存储器(CD-ROM)、数字多功能盘(DVD)、记忆棒、软盘、机械编码设备、例如其上存储有指令的打孔卡或凹槽内凸起结构、以及上述的任意合适的组合。这里所使用的计算机可读存储介质不被解释为瞬时信号本身，诸如无线电波或者其他自由传播的电磁波、通过波导或其他传输媒介传播的电磁波(例如，通过光纤电缆的光脉冲)、或者通过电线传输的电信号。

这里所描述的计算机可读程序指令可以从计算机可读存储介质下载到各个计算/处理设备，或者通过网络、例如因特网、局域网、广域网和/或无线网下载到外部计算机或外部存储设备。网络可以包括铜传输电缆、光纤传输、无线传输、路由器、防火墙、交换机、网关计算机和/或边缘服务器。每个计算/处理设备中的网络适配卡或者网络接口从网络接收计算机可读程序指令，并转发该计算机可读程序指令，以供存储在各个计算/处理设备中的计算机可读存储介质中。

用于执行本公开操作的计算机程序指令可以是汇编指令、指令集架构(ISA)指令、机器指令、机器相关指令、微代码、固件指令、状态设置数据、或者以一种或多种编程语言的任意组合编写的源代码或目标代码，所述编程语言包括面向对象的编程语言—诸如Smalltalk、C++等，以及常规的过程式编程语言—诸如“C”语言或类似的编程语言。计算机可读程序指令可以完全地在用户计算机上执行、部分地在用户计算机上执行、作为一个独立的软件包执行、部分在用户计算机上部分在远程计算机上执行、或者完全在远程计算机或服务器上执行。在涉及远程计算机的情形中，远程计算机可以通过任意种类的网络—包括局域网(LAN)或广域网(WAN)—连接到用户计算机，或者，可以连接到外部计算机(例如利用因特网服务提供商来通过因特网连接)。在一些实施例中，通过利用计算机可读程序指令的状态信息来个性化定制电子电路，例如可编程逻辑电路、现场可编程门阵列(FPGA)或可编程逻辑阵列(PLA)，该电子电路可以执行计算机可读程序指令，从而实现本公开的各个方面。

这里参照根据本公开实施例的方法、装置(系统)和计算机程序产品的流程图和/或框图描述了本公开的各个方面。应当理解，流程图和/或框图的每个方框以及流程图和/或框图中各方框的组合，都可以由计算机可读程序指令实现。

这些计算机可读程序指令可以提供给通用计算机、专用计算机或其它可编程数据处理装置的处理器，从而生产出一种机器，使得这些指令在通过计算机或其它可编程数据处理装置的处理器执行时，产生了实现流程图和/或框图中的一个或多个方框中规定的功能/动作的装置。也可以把这些计算机可读程序指令存储在计算机可读存储介质中，这些指令使得计算机、可编程数据处理装置和/或其他设备以特定方式工作，从而，存储有指令的计算机可读介质则包括一个制造品，其包括实现流程图和/或框图中的一个或多个方框中规定的功能/动作的各个方面的指令。

也可以把计算机可读程序指令加载到计算机、其它可编程数据处理装置、或其它设备上，使得在计算机、其它可编程数据处理装置或其它设备上执行一系列操作步骤，以产生计算机实现的过程，从而使得在计算机、其它可编程数据处理装置、或其它设备上执行的指令实现流程图和/或框图中的一个或多个方框中规定的功能/动作。

附图中的流程图和框图显示了根据本公开的多个实施例的系统、方法和计算机程序产品的可能实现的体系架构、功能和操作。在这点上，流程图或框图中的每个方框可以代表一个模块、程序段或指令的一部分，所述模块、程序段或指令的一部分包含一个或多个用于实现规定的逻辑功能的可执行指令。在有些作为替换的实现中，方框中所标注的功能也可以以不同于附图中所标注的顺序发生。例如，两个连续的方框实际上可以基本并行地执行，它们有时也可以按相反的顺序执行，这依所涉及的功能而定。也要注意的是，框图和/或流程图中的每个方框、以及框图和/或流程图中的方框的组合，可以用执行规定的功能或动作的专用的基于硬件的系统来实现，或者可以用专用硬件与计算机指令的组合来实现。

以上已经描述了本公开的各实施例，上述说明是示例性的，并非穷尽性的，并且也不限于所披露的各实施例。在不偏离所说明的各实施例的范围和精神的情况下，对于本技术领域的普通技术人员来说许多修改和变更都是显而易见的。本文中所用术语的选择，旨在最好地解释各实施例的原理、实际应用或对市场中的技术改进，或者使本技术领域的其它普通技术人员能理解本文披露的各实施例。