一种地震记录迭代型组合域去噪方法及系统与流程

本发明属于地震探测数据处理技术领域，尤其涉及一种地震记录迭代型组合域去噪方法及系统。

背景技术：

目前，业内常用的现有技术是这样的：

地震勘探一直是石油与天然气资源探测、矿产资源勘查、煤田和工程地质勘查、区域地质及地壳研究的重要物探方法，如反射波、折射波及瑞雷面波法等在探测不同深度的各类地质目标体方面有着不可取代的作用。然而，随着勘探目标要求越来越精细，如何在强能量的各种复杂干扰噪声湮没环境中提取微弱的目标信号是一项艰巨的任务。

当前，研究与开发的用于地震记录去噪方法技术颇多，理论与应用研究亦比较广泛。如常规的f-k(频率-波数)域滤波方法(wardenetal.，2012)、压制随机噪声的f-x(频率-空间)域预测反褶积方法(张军华等，2006)、高精度randon变换及其变形变换(线性的、多项式的)去噪方法(牛滨华等，2001；黄新武等，2003；wardenetal.，2012)、用于随机噪声聚类分析去噪的k-l(karhunen-lòeve)变换(ribaricandfratric，2005)、用于线性噪声压制的τ-p变换(王维红等，2006)等等。但是上述传统地震记录去噪技术多基于简单的最优化或正交的变换算法，其形成的是固定滤波窗口，具有特定噪声去噪的局限性。同时，固定式阀值滤波不适应于微弱有效波的保真，亦容易去噪过度，导致地震信号失真(shanetal.，2009)。因此，有必要对传统的地震记录去噪处理解释技术做进一步的改进创新。

“保真去噪”是地震记录数据处理的目标，即最大限度的保留目标信号的强度、形态及特征，同时尽可能去除干扰解释的噪声信号。傅里叶分析及窗口傅里叶变换提供了信号的整体频谱，并克服了时频局部化方面的缺陷，但仍使用非自适应时频窗口。wavelet(小波)理论的提出一定程度上弥补了上述缺点，即可根据信号的特点提供自适应的时频窗口，并提供多频带、多尺度的信号分析能力(saraccoetal.，2007；maandplonka，2010；maurietal.，2012)。对于地震记录去噪而言，每一个地震道可视为一个具有伸缩和平移特性的一维wavelet基函数，通过在特定频段的稀疏表示以达到保真去噪的目的。然而，一维wavelet变换具有明显的方向选择性，其没有充分利用二维地震剖面特性的优势，故去噪效果局限于单个数据道(roueffetal.，2006；shahetal.，2009)。二维wavelet变换的出现一定程度上改善了wavelet变换去噪的效果，其本质上可表述为一维离散wavelet变换的矢量积分，因而二维wavelet变换具有水平、对角及垂向上的去噪能力(xueetal.，2015；vuongetal.，2015)。当地震记录数据中含有多方向的线状特征信息，wavelet变换不能提供有效的表征和检测，无法对原始信号进行最优化非线性重建(ventosaetal.，2012；chopraandmarfurt，2015)。

wavelet域高阶相关统计叠加技术的提出和发展，为消除地震记录中的噪声并提高地震剖面的信噪比提供了新的技术手段(张三宗和徐义贤，2006；lietal.，2013)。该算法认为有用信号主要分布在wavelet域小尺度范围，而噪声信号主要由大尺度范围控制，通过不同尺度分离后叠加；同时采用高阶相关系数作为叠加权系数，达到突显有用信号而压制噪声信号的目的。相关统计叠加技术在地震记录去噪的应用，从二阶发展到高阶理论，获得广泛的应用并取得较好的去噪效果。但是，地震记录中有用信号的高阶相关系数的获得取决于原始记录数据具有良好的同相轴一致性，对于缺失或者受破坏的实测地震剖面，有效地震信号的高阶系数不明显，导致统计叠加后的地震剖面形成空白区域，影响后续准确的地质解释。因而，对于复杂地震剖面的“保真去噪”，wavelet域高阶相关统计叠加技术可有效地用于去噪，而无法高效地保真(张恒磊等，2010)。

随后，ridgelet(脊波)变换被提出并结合radon变换和wavelet变换的优势，具有很强的方向选择性和识别能力，更为有效地改善了wavelet域去噪后的模糊现象。但对于地震记录曲线边缘的特征描述，其逼近程度低及计算效率低(zhangetal.，2009)。curvelet(曲波)变换可视为描述一系列ridgelet变换组合特征，除了具备wavelet域伸缩及平移特性，其通过增加一个方向参量而具有更好的方向识别能力。curvelet变换属于多尺度分析，可用于高维信号的稀疏表示，因而该方法自提出便迅速广泛应用于地球物理领域，特别是针对地震记录数据噪声压制、多次波分离等领域，研究成果不断涌现(candèsanddonoho，2004；candèsetal.，2006；张恒磊等，2010)。然而，curvelet变换去噪的关键点是获取较为准确的阀值参数，针对不同的噪声类型，所需阀值亦有所不同。尽管当前发展了诸多计算阀值的方法，如分段线性滤波方法、l2标准差法、curvelet正变换法、对角实数及复数阀值法等，但不同阀值法具有实用局限性，在实测地震记录中并不具普遍性。因而，curvelet变换总的来说达到了“保真去噪”的保真效果，但需要估计准确的阀值参数，无法实现普遍的去噪目的(仝中飞等，2008；shanetal.，2009)。

为此，人们自然而然提出了wavelet和curvelet多域组合滤波算法，即利用不同域的去噪算法，顺序迭代对地震记录进行处理，或迭代融合不同域去噪算法，期望达到去除噪声的同时，有效保持目标地震记录信息，实现地震记录的保真处理。多域组合滤波算法的发展较为成熟，研究成果颇多，如starck等提出用于图像去噪的多域组合滤波法、saevarsson等为sar数据设计的组合域去噪技术、马建伟等进一步将starck等提出的方法引入到地震记录去噪；韩立国等提出curvelet-中值滤波、curvelet-radon变换等组合域滤波算法并成功应用于金属矿勘探地震记录去噪；徐颖等将多域组合去噪技术用于油气勘探地震记录去噪，张恒磊等提出基于高阶相关的curvelet域和空间域的倾角扫描噪声压制方法。总体上，现有的多域组合去噪算法多为curvelet域与传统去噪算法、wavelet变换在处理顺序上的结合，从而导致整体的去噪效果接近于单个curvelet变换的去噪效果，无法满足包含复杂噪声的地震记录去噪需求(shanetal.，2009)。

综上所述，随着地震勘探目标越发复杂多变，地震记录含躁性越来越复杂，传统地震记录去噪方法技术已不能满足精细去噪要求，发展高效的“保真去噪”算法迫在眉睫。wavelet和curvelet组合域去噪技术是上述发展的方向之一，近年来研究成果颇丰，但多波域算法有机融合去噪技术的研究程度不高，仍需要加强此方向的理论与应用研究。

综上所述，现有技术存在的问题是：

(1)传统地震记录去噪方法技术已不能满足精细去噪要求；而多波域算法有机融合去噪技术的研究程度不高。

(2)地震记录中有用信号的高阶相关系数的获得取决于原始记录数据具有良好的同相轴一致性，对于缺失或者受破坏的实测地震剖面，有效地震信号的高阶系数不明显，导致统计叠加后的地震剖面形成空白区域，影响后续准确的地质解释。因而，对于复杂地震剖面的“保真去噪”，wavelet域高阶相关统计叠加技术可有效地用于去噪，而无法高效地保真。

(3)curvelet变换总的来说达到了“保真去噪”的保真效果，但需要估计准确的阀值参数，无法实现普遍的去噪目的。

(4)现有的多域组合去噪算法多为curvelet域与传统去噪算法、wavelet变换在处理顺序上的结合，从而导致整体的去噪效果接近于单个curvelet变换的去噪效果，无法满足包含复杂噪声的地震记录去噪需求。

解决上述技术问题的难度和意义：

针对地震记录wavelet高阶相关统计叠加去噪容易导致信号失真、curvelet变换阀值去噪局限于特定阀值函数的获取，而现有的多域组合去噪技术无法有效滤除复杂地震记录噪声，本研究拟将wavelet高阶相关统计叠加算法和curvelet变换阀值去噪技术有机结合，提出迭代型和交叉型地震记录组合域去噪技术，实现高效的地震记录“保真去噪”。

技术实现要素：

针对现有技术存在的问题，本发明提供了一种地震记录迭代型组合域去噪方法及系统。

本发明是这样实现的，一种地震记录迭代型组合域去噪方法，所述地震记录迭代型组合域去噪方法将wavelet域高阶相关统计叠加理论与curvelet域变换有机结合；利用基于wavelet和curvelet小波基函数的curvelet域阀值去噪算法，结合wavelet域高阶相关统计叠加技术，实现迭代型组合域去噪。

进一步，所述地震记录迭代型组合域去噪方法进一步包括：

利用curvelet变换算法，引入wavelet域高阶相关统计叠加技术进行curvelet不同尺度、不同角度的变换系数相关性扫描，实现迭代型组合域去噪。

进一步，wavelet高阶相关统计叠加技术地震记录去噪过程包括：

1)对地震记录进行一维或二维wavelet正变换，获取wavelet域参数；

2)利用高阶相关统计理论获取wavelet系数的相关参数；

3)根据已获取的相关系数在不同尺度叠加，增强目标wavelet系数的相关性；

4)采用一维或二维wavelet反变换重建去噪后的地震记录。

进一步，curvelet变换去噪算法为：对地震记录进行二维快速傅里叶正变换获取系数分布，形成尺度和角度窗口，根据给定阀值函数截断去噪，利用二维快速傅里叶反变换获取去噪结果；则尺度窗口w和角度窗口v表述为：

式中，j为尺度，m为旋转方向，

r和t分别为空间及时间域参数。角度窗口v为环形域，并满足|x|∈[2^j，2^j+1]及极坐标定义θj，m＝2πm·2^-[j/2]；而频率域的窗口统一定义为：

uj(r，θ)＝2^-3/4w(2^-jr)v(2[j/2]·θ/2π)；

对尺度j、旋转角度θj，m及空间位置则curvelet变换有：

相应的旋转角度矩阵为：

对f(x)∈l²(r²)，curvelet系数为：

对于合成地震记录，根据加入的噪声类型及强度给出准确的阀值函数；对于实测地震记录数据，采用l2标准差法估计阀值函数。

进一步，迭代型组合域去噪方法具体包括：

1)在小尺度顺序迭代使用curvelet基函数、wavelet基函数对地震记录进行curvelet变换阀值去噪，拟达到拾取地震记录复杂曲线、线状特征，并通过估计的阀值函数滤除噪声；

2)对curvelet变换阀值去噪后的地震记录，再顺序迭代使用wavelet高阶相关统计叠加算法去噪；迭代型组合域去噪为通过在小尺度中使用不同的基函数，curvelet变换阀值去噪在curvelet基函数的情况残留点状噪声成份，wavelet基函数的情况残留线状噪声成份；再引入高阶相关统计理论，拾取目标地震信号，去除噪声成份。

进一步，步骤2具体包括：

步骤一：设置curvele系数所有尺度j和选择角度θ的硬阀值函数值δj，θ；获取大尺度阀值ej，θ＝3δj，θ和小尺度阀值ej，θ＝4δj，θ；

步骤二：采用curvelet及wavelet基函数，对原始地震记录f做curvelet变换获取系数分布和使用阀值函数和采用curvelet反变换获取去噪的地震记录和

步骤三：fori＝1，…，mdo；

对和做wavelet变换，获取系数分布；

提取近似系数和详细系数；

对wavelet域中所有的尺度：

forjc＝1，…，jcdo；

计算相关系数并归一化

进行地震记录叠加

end；

采用wavelet反变换获取去噪结果；

i＝i+1；

end。

本发明的另一目的在于提供一种实现所述地震记录迭代型组合域去噪方法的计算机程序。

本发明的另一目的在于提供一种实现所述地震记录迭代型组合域去噪方法的信息数据处理终端。

本发明的另一目的在于提供一种计算机可读存储介质，包括指令，当其在计算机上运行时，使得计算机执行所述的地震记录迭代型组合域去噪方法。

本发明的另一目的在于提供一种实现所述地震记录迭代型组合域去噪方法的地震记录迭代型组合域去噪控制系统。

综上所述，本发明的优点及积极效果为：

本发明将wavelet域高阶相关统计叠加理论与curvelet域变换有机结合，提出一种迭代型组合域地震记录去噪技术，即基于wavelet和curvelet小波基函数的curvelet域阀值去噪算法，拾取地震记录复杂曲线、线状特征，并通过阀值函数滤除噪声；结合wavelet域高阶相关统计叠加技术，根据点状特征(椒盐式)噪声并不具备高阶相关性质，而有效地震记录具备较好的高阶相关性质，从而利用高阶相关值的大小拾取目标地震波有效信号，而剔除点状特征噪声信号，提出并实现迭代型组合域去噪技术；利用合成含躁地震记录及实测地震记录去噪实验验证本发明提出的迭代型组合域去噪算法的有效性及可行性。项目研究成果为复杂含躁地震记录的去噪处理提供了新的技术手段，为组合域去噪算法的实用化奠定了理论基础。

附图说明

图1是本发明实施例提供的地震记录迭代型组合域去噪方法流程图。

图2是本发明实施例提供的含随机噪声地震记录组合域算法与单独域算法去噪结果图。

图中：(a)含躁数据、(b)curvelet小波基去噪、(c)curvelet曲波基去噪、(d)组合域去噪、(e)小波变换去噪。

图3是发明实施例提供的含相关噪声地震记录组合域算法与单独域算法去噪结果图。

图中：(a)含躁数据、(b)curvelet小波基去噪、(c)curvelet曲波基去噪、(d)组合域去噪、(e)小波变换去噪。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

如图1，本发明实施例提供的地震记录迭代型组合域去噪方法包括：

将wavelet域高阶相关统计叠加理论与curvelet域变换有机结合；利用基于wavelet和curvelet小波基函数的curvelet域阀值去噪算法，结合wavelet域高阶相关统计叠加技术，实现迭代型组合域去噪。

所述地震记录迭代型组合域去噪方法进一步包括：

利用curvelet变换算法，引入wavelet域高阶相关统计叠加技术进行curvelet不同尺度、不同角度的变换系数相关性扫描，实现迭代型组合域去噪。

下面结合具体实施例对本发明的应用作进一步描述。

1.wavelet高阶相关统计叠加技术、curvelet变换地震记录去噪实用性评价：

本发明拟结合wavelet高阶相关统计叠加技术与curvelet变换算法开展地震记录组合域去噪技术研究与开发，因而首先需明确用于组合的单独算法地震记录去噪的实用性及适用范围，以取长补短，研发更有效的地震记录组合域去噪算法。

wavelet高阶相关统计叠加技术地震记录去噪过程主要为：1)对地震记录进行一维或二维wavelet正变换，获取wavelet域参数；2)利用高阶相关统计理论获取wavelet系数的相关参数；3)根据已获取的相关系数在不同尺度叠加，增强目标wavelet系数的相关性；4)采用一维或二维wavelet反变换重建去噪后的地震记录。其中，未偏移的三阶相关函数表示为(张三宗和徐义贤，2006；lietal.，2013)：

归一化为：

则加权的wavelet系数为：

上式中，i为地震剖面道数，t为时间采样点，p为平移因子，一般取值为1，fi^k为第i道地震记录、第k个wavelet域尺度正变换系数。

本发明拟采用第二代截断离散curvelet变换去噪算法(candèsanddonoho，2004；candèsetal.，2006)开展组合域去噪技术研究。其主要思路为：对地震记录进行二维快速傅里叶正变换获取系数分布，形成尺度和角度窗口，根据给定阀值函数截断去噪，利用二维快速傅里叶反变换获取去噪结果。则尺度窗口w和角度窗口v可表述为：

式中，j为尺度，m为旋转方向，r和t分别为空间及时间域参数。角度窗口v为环形域，并满足|x|∈[2^j，2^j+1]及极坐标定义θj，m＝2πm·2^-[j/2]。而频率域的窗口统一定义为：

uj(r，θ)＝2^-3/4w(2^-jr)v(2[j/2]·θ/2π)(6)

对尺度j、旋转角度θj，m及空间位置则curvelet变换可有：

相应的旋转角度矩阵为：

对f(x)∈l²(r²)，curvelet系数定义为：

curvelet变换地震记录去噪技术的关键是获取较为准确的阀值函数。对于合成地震记录，可根据加入的噪声类型及强度给出准确的阀值函数；对于实测地震记录数据，本发明拟采用l2标准差法(candèsetal.，2006)估计阀值函数。

2构建地震记录迭代型组合域去噪技术：如图1迭代型组合域去噪的主要过程为两步：

s101：在小尺度顺序迭代使用curvelet基函数、wavelet基函数对地震记录进行curvelet变换阀值去噪，拟达到拾取地震记录复杂曲线、线状特征，并尽可能通过估计的阀值函数滤除噪声；

s102：对curvelet变换阀值去噪后的地震记录，再顺序迭代使用wavelet高阶相关统计叠加算法去噪。迭代型组合域去噪核心思想是通过在小尺度中使用不同的基函数，curvelet变换阀值去噪在curvelet基函数的情况残留点状噪声成份，而wavelet基函数的情况残留线状噪声成份；再引入高阶相关统计理论，拾取目标地震信号，去除噪声成份。

其基本步骤如下：

步骤一：设置curvele系数所有尺度j和选择角度θ的硬阀值函数值δj，θ；

获取大尺度阀值ej，θ＝3δj，θ和小尺度阀值ej，θ＝4δj，θ；

步骤二：采用curvelet及wavelet基函数，对原始地震记录f做curvelet变换

获取系数分布和使用阀值函数和采用curvelet反变换获取去噪的地震记录和

步骤三：fori＝1，…，mdo.

对和做wavelet变换，获取系数分布；

提取近似系数和详细系数；

对wavelet域中所有的尺度：

iforjc＝1，…，jcdo.

按公式(1)计算相关系数并按公式(2)归一化

ii按公式(3)进行地震记录叠加

end。

iii采用wavelet反变换获取去噪结果。

i＝i+1；

end。

下面结合实验结果对本发明的应用作进一步描述。

通过采用上述方法开展地震数据含随机噪声(图2)、含相关噪声(图3)所示。

图2为本发明实施例提供的含随机噪声地震记录组合域算法与单独域算法去噪结果。(a)含躁数据、(b)curvelet小波基去噪、(c)curvelet曲波基去噪、(d)组合域去噪、(e)小波变换去噪。

图3本发明实施例提供的含相关噪声地震记录组合域算法与单独域算法去噪结果。(a)含躁数据、(b)curvelet小波基去噪、(c)curvelet曲波基去噪、(d)组合域去噪、(e)小波变换去噪。

本发明完成基于wavelet域高阶相关统计叠加技术、xurvelet域阀值去噪算法，开展地震合成含躁记录、实测含躁记录数值去噪效果对比分析，明确组合波域去噪技术在地震记录处理中的有效性及适用范围。结合wavelet高阶相关统计叠加技术与curvelet变换算法开展地震记录组合域去噪技术分析与开发，明确用于组合的单独算法地震记录去噪的实用性及适用范围，以取长补短，研发更有效的地震记录组合域去噪算法。

在上述实施例中，可以全部或部分地通过软件、硬件、固件或者其任意组合来实现。当使用全部或部分地以计算机程序产品的形式实现，所述计算机程序产品包括一个或多个计算机指令。在计算机上加载或执行所述计算机程序指令时，全部或部分地产生按照本发明实施例所述的流程或功能。所述计算机可以是通用计算机、专用计算机、计算机网络、或者其他可编程装置。所述计算机指令可以存储在计算机可读存储介质中，或者从一个计算机可读存储介质向另一个计算机可读存储介质传输，例如，所述计算机指令可以从一个网站站点、计算机、服务器或数据中心通过有线(例如同轴电缆、光纤、数字用户线(dsl)或无线(例如红外、无线、微波等)方式向另一个网站站点、计算机、服务器或数据中心进行传输)。所述计算机可读取存储介质可以是计算机能够存取的任何可用介质或者是包含一个或多个可用介质集成的服务器、数据中心等数据存储设备。所述可用介质可以是磁性介质，(例如，软盘、硬盘、磁带)、光介质(例如，dvd)、或者半导体介质(例如固态硬盘solidstatedisk(ssd))等。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。