CN106597539B - 针对黄土塬地区的曲波域Radon变换噪声压制方法 - Google Patents

Abstract

本发明公开了针对黄土塬地区的曲波域Radon变换噪声压制方法，包括：对原始地震数据进行曲波正变换，得到不同方向、不同尺度的曲波系数域的矩阵；对各矩阵逐个进行Radon正变换将曲波系数域的数据变换到Radon域，对各矩阵逐个进行阈值滤波，由于干扰波在Radon域表现为能量较强的点或能量团，随机噪音表现为不聚焦能量，而有效波虽然不聚焦，但仍然表现出较随机噪声较大的聚焦能量，通过阈值滤波，干扰波能得到有效的剔除；将数据通过Radon反变换从Radon域变换回曲波系数域；对曲波域的数据进行曲波反变换，得到最终数据。通过本发明所述方法能够很好的压制原始资料中强能量相干噪声，大大提高地震资料的信噪比，为后续地震解释人员提供有效的、更多的、高品质的地震信息。

Description

针对黄土塬地区的曲波域Radon变换噪声压制方法

技术领域

本发明属于油气勘探地震资料处理技术领域，特别涉及针对黄土塬地区的曲波域Radon变换噪声压制方法。

背景技术

在地震资料处理过程中噪声压制方法，主要可以分为三大类：1.基于时间—空间域噪声能量、波形特征、速度以及道间关系等采取的类似于外科手术切除方法，通过先识别噪声干扰，然后从原始地震记录中切除以达到压制的目的；2.基于变换域的方法，常见的有傅立叶变换、FK变换、Radon变换、小波变换等等，将原始信号从时间空间域变换到某个域以便于更好的分离出干扰噪声。3.基于波场变换的方法，常见的有SRME(自由表面多次波衰减)、波场延拓去噪等，通过已知速度模型，将波长延拓到高速层顶部并将时间为零的波场(表征噪声的波场)置零，最后延拓回基准面，从而消除噪声。第一类噪声压制方法去噪效果较好，具有高度保幅等优点，但一般需要较为复杂的交互操作，依赖于处理人员处理经验。第二类方法，一般来说要求处理参数较少且不需要太多的交互操作，也能达到较好的处理效果，因此在实际处理中运用广泛、方法种类颇多。但是有效信号和噪音在很多情况下通过这些变换也无法区分，造成去噪效果不理想。第三类噪声压制方法能够较为彻底的分离原始记录不同特征的噪声，但是它依赖于速度模型的精度，而建立精确的速度模型的在地震资料处理难以达到。

山前黄土塬地区地表条件复杂，山前黄土覆盖、地形起伏高差达、表层岩性变化剧烈、工区土质结构疏松，地震波路径复杂、低降速带厚度、速度变化大等因素，对地震资料运动学、动力学特征等带来了一系列影响。采用常规的上述三种地震资料处理过程中的噪声压制方法效果都不好。因此，如何消除山前黄土塬影响造成的干扰波问题，从而加深对山前黄土塬地区地下地质特征的认识，对山前黄土塬复杂构造区勘探具有重要意义。

发明内容

本发明的目的在于，针对现有技术的不足，提供一种针对黄土塬地区的曲波域Radon变换噪声压制方法。

本发明在用的技术方案如下。

针对黄土塬地区的曲波域Radon变换噪声压制方法，包括如下步骤：

步骤1：对原始地震数据进行曲波正变换，得到不同方向、不同尺度的曲波系数域的矩阵，此时，在曲波系数域中线性特征较好的干扰波在矩阵中仍然保持较好的线性性，而有效波则在矩阵中被分散开。

步骤2：对各尺度角度矩阵逐个进行Radon正变换将曲波系数域的数据变换到Radon域，对各尺度角度矩阵逐个进行阈值滤波；由于干扰波在Radon域表现为能量较强的点或能量团，随机噪音表现为不聚焦能量，而有效波虽然不聚焦，但仍然表现出较随机噪声较大的聚焦能量，通过阈值滤波，干扰波能得到有效的剔除。

步骤3：将步骤2中的数据通过Radon反变换从Radon域变换回曲波系数域；对曲波域的数据进行曲波反变换，得到最终数据。

进一步，在步骤1中，所述对原始地震数据进行曲波正变换包括：对工区具有代表性噪声的单炮进行分析及确定噪音主要的系数矩阵位置。单炮进行分析及确定噪音主要的系数矩阵位置的目的是便于后期进行处理时可以有选择性地进行操作。将工区不同位置单炮记录进行曲波正变换，并选取整条测线不同位置的具有代表性的单炮分析其噪声和有效信号在曲波域不同方向、不同尺度的特征。

进一步，步骤1中，所述曲波正变换时，首先，将曲波变换的系数定义为Curvelet系数C(i,j,k)；

c(i,j,k)＝∫f(w)U_j,l(w)dw，

其中：窗口函数定义为

尺度窗口定义为W_j(ω)＝W(2^-jω)，

方向窗口定义为V_j(ω)＝V(2^[j/2]ω₂/ω₁)；

对于一个二维数据来说，曲波正变换历以下五个步骤：

1)对数据进行二维傅立叶变换，将时间-空间域地震数据变换到频率-波数域；

2)对频率、波数域数据重采样或内插值；

3)将频率-波数域数据与窗函数U_j,l相乘；

4)绕原点Wrapping局部化；

5)进行二维傅立叶逆变换；得到Curvelet系数C(i,j,k)。

进一步，对于一个地震道集d(t,x)，所述Radon正变换、Radon反变换定义如下：

Radon正变换到τ-p域：

Radon反变换到t-x域：

把连续线性双曲线拉东变换积分变为有限区域内的相加，从而可以得到离散形式的线性拉东变换记作：

Radon正变换到τ-p域：

Radon反变换到t-x域：

其中：N_x----地震记录道数；

d_x---------------道间距；

dp-------拉东域采样间隔为；

Np---------拉东域记录道数；

对应的也可以写成矩阵形式：

Radon正变换到τ-p域：u＝Ld

Radon反变换到t-x域：d＝L^Hu

L和L^H为拉东变换算子，L为大小为(Np×Nt)×(Nx×Nt)的正变换算子，它将时域的抛物线变换到Radon域收敛为一个点；L^H为L的共轭转置算子，即反变换算子，将拉东域的一个点反变换到时间空间域；

通常为了提高拉东域的分辨率，在求解时，通常先定义拉东反变换，即使目标函数最小，对u求导得到最小平方解：u＝(LL^H+λ²I)^-1L^Hd；

其中λ是使求解稳定的范围在0.1-1的阻尼因子，I为单位矩阵。

进一步，在步骤2中，所述Radon正变换为对曲波域不同角度、不同尺度的的矩阵数据进行拉东变换，并分析噪声和有效信号在拉东域的对应关系及其特征。分析噪声和有效信号在拉东域的对应关系及其特征的目的是为后面设置滤波器提供参考。

进一步，在步骤2中，所述阈值滤波是指设置选取合适的阈值，将大于该阈值的数据置零而小于阈值的数据保持不变。

进一步，所述阈值滤波的具体步骤是：

初步选定一个阈值，将各矩阵中大于该阈值的数据置零而小于阈值的数据保持不变；对滤波结果进行检验；

当阈值滤波结果可以满足要求时，阈值滤波结束；

当阈值滤波结果不能满足要求时，改变选定的阈值的大小并按照重新设定的阈值完成对各矩阵新一轮的迭代计算，直到阈值滤波结果可以满足噪声压制要求时为止。

进一步，阈值滤波结果可以满足要求是指：计算相邻两次计算数据样点振幅差的平方和，如果当前计算平方和与上一次计算平方和差异小于1％,则退出迭代。

本发明的有益效果是：为专门针对于黄土塬低信噪比地区提出的噪声压制方法，同时也适用于地表条件复杂地区噪声的压制，通过该方法能够很好的压制原始资料中强能量相干噪声，大大提高地震资料的信噪比，为后续地震解释人员提供有效的、更多的、高品质的地震信息。

附图说明

图1为曲波变换示意图。

图2为Wrapping示意图。

图3为曲波域Radon变换去噪流程图。

图4为本文提出的曲波域Radon变换单炮去噪应用效果图。

图5为本文提出方法实际案例应用得到的叠加剖面图。

具体实施方式

下面，结合附图和实施例，对本发明做进一步说明。

实施例1。针对黄土塬地区的曲波域Radon变换噪声压制方法，包括如下步骤：

步骤1：对原始地震数据进行曲波正变换，得到不同方向、不同尺度的曲波系数域的矩阵，此时，在曲波系数域中线性特征较好的干扰波在矩阵中仍然保持较好的线性性，而有效波则在矩阵中被分散开；

步骤2：对各尺度角度矩阵逐个进行Radon正变换将曲波系数域的数据变换到Radon域，对各尺度角度矩阵逐个进行阈值滤波；由于干扰波在Radon域表现为能量较强的点或能量团，随机噪音表现为不聚焦能量，而有效波虽然不聚焦，但仍然表现出较随机噪声较大的聚焦能量，通过阈值滤波，干扰波能得到有效的剔除；

步骤3：将步骤2中的数据通过Radon反变换从Radon域变换回曲波系数域；对曲波域的数据进行曲波反变换，得到最终数据。

在步骤1中，所述对原始地震数据进行曲波正变换包括：对工区具有代表性噪声的单炮进行分析及确定噪音主要的系数矩阵位置，便于后期进行处理时可以有选择性地进行操作；将工区不同位置单炮记录进行曲波正变换，并选取整条测线不同位置的具有代表性的单炮分析其噪声和有效信号在曲波域不同方向、不同尺度的特征。

步骤1中，所述曲波正变换时，首先，将曲波变换的系数定义为Curvelet系数C(i,j,k)；

c(i,j,k)＝∫f(w)U_j,l(w)dw，

其中：窗口函数定义为

尺度窗口定义为W_j(ω)＝W(2^-jω)，

方向窗口定义为V_j(ω)＝V(2^[j/2]ω₂/ω₁)；

对于一个二维数据来说，曲波正变换历以下五个步骤：

1)对数据进行二维傅立叶变换，将时间-空间域地震数据变换到频率-波数域；

2)对频率、波数域数据重采样或内插值；

3)将频率-波数域数据与窗函数U_j,l相乘；

4)绕原点Wrapping局部化；

5)进行二维傅立叶逆变换；得到Curvelet系数C(i,j,k)。

对于一个地震道集d(t,x)，所述Radon正变换、Radon反变换定义如下：

Radon正变换到τ-p域：

Radon反变换到t-x域：

把连续线性双曲线拉东变换积分变为有限区域内的相加，从而可以得到离散形式的线性拉东变换记作：

Radon正变换到τ-p域：

Radon反变换到t-x域：

其中：N_x----地震记录道数；

d_x---------------道间距；

dp-------拉东域采样间隔为；

Np---------拉东域记录道数；

对应的也可以写成矩阵形式：

Radon正变换到τ-p域：u＝Ld

Radon反变换到t-x域：d＝L^Hu

L和L^H为拉东变换算子，L为大小为(Np×Nt)×(Nx×Nt)的正变换算子，它将时域的抛物线变换到Radon域收敛为一个点；L^H为L的共轭转置算子，即反变换算子，将拉东域的一个点反变换到时间空间域；

通常为了提高拉东域的分辨率，在求解时，通常先定义拉东反变换，即使目标函数最小，对u求导得到最小平方解：u＝(LL^H+λ²I)^-1L^Hd；

其中λ是使求解稳定的范围在0.1-1的阻尼因子，I为单位矩阵。

在步骤2中，所述Radon正变换为对曲波域不同角度、不同尺度的的矩阵数据进行拉东变换，并分析噪声和有效信号在拉东域的对应关系及其特征，为后面设置滤波器提供参考。

在步骤2中，所述阈值滤波是指设置选取合适的阈值，将大于该阈值的数据置零而小于阈值的数据保持不变。

所述阈值滤波的具体步骤是：初步选定一个阈值，将各矩阵中大于该阈值的数据置零而小于阈值的数据保持不变；对滤波结果进行检验；

当阈值滤波结果可以满足要求时，阈值滤波结束；

当阈值滤波结果不能满足要求时，改变选定的阈值的大小并按照重新设定的阈值完成对各矩阵新一轮的迭代计算，直到阈值滤波结果可以满足噪声压制要求时为止。

阈值滤波结果可以满足要求是指：计算相邻两次计算数据样点振幅差的平方和，如果当前计算平方和与上一次计算平方和差异小于1％,则退出迭代。

在山前黄土源复杂地表条件下的地震勘探中，首先由于这些地区表层被黄土覆盖，黄土结构疏松，厚度在纵向和横向上变化剧烈，地震波的吸收和衰减较为严重。其次加上其复杂的地形条件，使得得到的地震资料混有大量的噪音，干扰类型多，能量强，线性规律差，这些因素使得从野外采集的地震资料信噪比较低并且强的噪音干扰会使地震资料的处理难以进行。除常见的面波、多次波、声波、线性干扰、相干干扰及背景随机噪声干扰外，强能量、分布范围广的线性特征的散射干扰显得尤为突出，往往使得有效信号淹没其中。因此采用合适的噪声压制方法来提高该地区地震资料信噪比显得尤为重要。由于黄土塬地区噪声干扰的复杂性，难以从单一的频率、速度、时间、空间等方面去识别去除噪声，因此目前商业软件中常规的去噪方法难以有效的压制噪声。为了适应地震勘探的需求，必须开展针对黄土塬地区新的去噪方法研究。曲波变换理论提出早期主要应用于图像处理中，由于它具有多分辨率、多尺度、多方向等诸多有点，在地震资料处理中得到了快速的发展，目前曲波变换运用到噪声压制中，主要为压制随机噪声，也有压制面波的相关文献，但是整体来说方法不成熟。实践表明，Radon变换在压制线性干扰、相干干扰方面具有较好的效果。而对于针对黄土塬、山前低信噪比地区曲波域Radon换噪声压制方法目前国内外还没有相关方面文献资料对其进行介绍、阐述。发明人经过研究发现，曲波变换是一种充分考虑了地震记录波场特征，能够从频率、视速度、空间、时间各个方向进行信号分离的一种变换手段。经过曲波变换后，获得的曲波系数能够在多个维度上进行信号和噪音的分离。拉东变换可以将符合一定传播规律的信号聚焦，在曲波系数域使用拉东变换可以将信号和噪音进一步分离，从而获得高保真的信噪分离方法。在本实施例，针对黄土源低信噪比地区噪声传播特征提出适应于该地区的噪声压制方法，通过线性Radon变换，在去除地震面波、直达波等线性干扰中达到了较好的效果，其基本原理是将数据中的线性信息聚焦后进行噪音剔除。由于曲波变换具有较好的方向、尺度和位置分辨率，可以将各种噪音在曲波系数域里得到较好的区分。本实施例结合曲波变换与Radon变换的这种特性，将二者结合起来进行去噪，去噪效果得到大幅度提升，为进一步的处理解释工作提供更可靠的信息。

实施例2。针对黄土塬地区的曲波域Radon变换噪声压制方法，包括如下步骤：

①首先对原始地震数据进行曲波正变换，得到不同方向、不同尺度的矩阵。在曲波域中线性特征较好的干扰波在矩阵中仍然保持较好的线性性，而有效波则被分散开；

②分别对不同的小矩阵作Radon变换，将曲波域的数据变换到Radon域，并设置阈值，将大于阈值的数据置零，小于阈值的数据保持不变。由于干扰波在Radon域表现为能量较强的点或能量团，随机噪音表现为不聚焦能量，而有效波虽然不聚焦，但仍然表现出较随机噪声较大的聚焦能量，通过处理，干扰波能得到有效的剔除；

③将Radon域阈值处理后的数据变换回曲波系数域，最后进行曲波反变换，得到最终数据。具体处理流程见图2(曲波域Radon变换去噪流程示意图)。

线性Radon变换，在去除地震面波、直达波等线性干扰中具有较好的效果，其基本原理是将数据中的线性信息聚焦后进行噪音剔除。而曲波变换具有较好的方向、尺度和位置分辨率，可以将各种噪音在曲波系数域里得到较好的区分。本实施例结合曲波变换与Radon变换的这种特性，将二者结合起来进行去噪。

曲波变换、拉东变换的基本原理如下。

一、曲波变换

自从上世纪90年代Curvelet变换理论提出以来，Curvelet变换经历第一代Curvelet变换、第二代Curvelet发展。Curvelet变换具有多尺度、多分辨率、多方向等特点并能够很好的稀疏表达地震信号,在地震处理中压制相干干扰、随机噪声、提高分辨率、地震数据插值等方面都得到了广泛应用。

首先，Curvelet变换定义为c(i,j,k)＝∫f(w)U_j,l(w)dw，其中窗口函数定义为其中尺度窗口定义为W_j(ω)＝W(2^-jω)，方向窗口定义为V_j(ω)＝V(2^[j/2]ω₂/ω₁)。

Wrapping快速算法对于一个二维数据来说，第二代Curvelet变换一般需要经历以下五个步骤：

(1)对数据进行二维傅立叶变换(2DFFT)，变换到频率、波数域。

(2)对频率、波数域数据重采样或内插。

(3)与窗函数U_j,l相乘，如图所示。

(4)绕原点Wrapping局部化，如图2所示。

(5)二维傅立叶逆变换(2DIFFT)得到Curvelet系数C(i,j,k)。

二、拉东变换

对于一个地震道集d(t,x)，对于线性拉东正反变换定义如下：

正变换到τ-p域：

反变换到t-x域：

把连续线性双曲线拉东变换积分变为有限区域内的相加，从而可以得到离散形式的线性拉东变换记作：

正变换到τ-p域：

反变换到t-x域：

其中：N_x----地震记录道数；

d_x---------------道间距；

dp-------拉东域采样间隔为；

Np---------拉东域记录道数。

对应的也可以写成矩阵形式：

正变换到τ-p域：u＝Ld

反变换到t-x域：d＝L^Hu

L和L^H为拉东变换算子，L为大小为(Np×Nt)×(Nx×Nt)的正变换算子，它将时域的抛物线变换到Radon域收敛为一个点。L^H为L的共轭转置算子，即反变换算子，将拉东域的一个点反变换到时间空间域。

通常为了提高拉东域的分辨率，在求解时，通常先定义拉东反变换，即使目标函数最小，对u求导得到最小平方解：u＝(LL^H+λ²I)^-1L^Hd

其中λ是使求解稳定的范围在0.1-1的阻尼因子，I为单位矩阵。

实施例3。本实施例以一个二维地震资料为目标靶区，应用本方法对该资料进行处理，以验证本方法的效果。该实际资料采用二维单线采集，观测系统为6090-110-20-110-6090，地震资料时间长度6000ms，时间采样间隔为2ms，采样点数3000，每炮道数为408道。采用上述方法对该资料进行处理。

(1)首先步骤1，对工区典型单炮分析及选取合适的处理参数。

将工区不同位置单炮记录进行曲波正变换(图1为曲波变换示意图，其中左图为曲波变换过程中将地震信号在多个维度上进行细致的划分，其中四边代表尺度、角度代表维度方向示意图)，并选取整条测线不同位置的具有代表性噪音的记录分析其噪声和有效信号在曲波域不同方向、不同尺度的特征。接下来对曲波域不同方向、不同尺度再进行拉东变换并分析噪声和有效信号在拉东域的对应关系及其特征，最后选取合适的拉东域阈值，迭代计算，从而选取合适的阈值以及迭代处理次数作为下一步批量处理的输入参数。图4为本文提出的曲波域Radon变换单炮去噪应用效果，其中左图为去噪前单炮记录，右图为去噪后单炮记录。

(2)步骤2，将步骤1优选的曲波域拉东变换处理参数应用于整个工区。如图5所示，其中左图为原始叠加剖面(去噪前叠加剖面)，右图为利用本文提出的适应于噪声干扰复杂地区的曲波域Radon变换去噪方法处理后得到的叠加剖面(去噪后叠加剖面)。可以看出处理后叠加剖面效果得到了有效的改善。实际处理效果证明，曲波域Radon变换噪声压制方法在黄土塬低信噪比地区能够取得较好的效果。

最后说明的是，以上优选实施例仅用以说明本发明的技术方案而非限制，尽管通过上述优选实施例已经对本发明进行了详细的描述，但本领域技术人员应当理解，可以在形式上和细节上对其作出各种各样的改变，而不偏离本发明权利要求书所限定的范围。

Claims (3)

1.针对黄土塬地区的曲波域Radon变换噪声压制方法，其特征在于，包括如下步骤：

步骤1：对原始地震数据进行曲波正变换，得到不同方向、不同尺度的曲波系数域的矩阵，此时，在曲波系数域中线性特征较好的干扰波在矩阵中仍然保持较好的线性性，而有效波则在矩阵中被分散开；

步骤2：对各尺度角度矩阵逐个进行Radon正变换将曲波系数域的数据变换到Radon域，对各尺度角度矩阵逐个进行阈值滤波；由于干扰波在Radon域表现为能量较强的点或能量团，随机噪音表现为不聚焦能量，而有效波虽然不聚焦，但仍然表现出较随机噪声较大的聚焦能量，通过阈值滤波，干扰波能得到有效的剔除；

步骤3：将步骤2中的数据通过Radon反变换从Radon域变换回曲波系数域；对曲波域的数据进行曲波反变换，得到最终数据；

在步骤1中，所述对原始地震数据进行曲波正变换包括：对工区具有代表性噪声的单炮进行分析及确定噪音主要的系数矩阵位置；将工区不同位置单炮记录进行曲波正变换，并选取整条测线不同位置的具有代表性的单炮分析其噪声和有效信号在曲波域不同方向、不同尺度的特征；

所述曲波正变换时，首先，将曲波变换的系数定义为Curvelet系数C(i，j，k)；

c(i,j,k)＝∫f(w)U_i,j(w)dw，

其中：窗口函数定义为U_i,j(w)＝W_i(w)V_j(w)，

尺度窗口定义为W_i(w)＝W(2^-iw),

方向窗口定义为V_j(ω)＝V(2^[j/2]ω₂/ω₁)；

对于一个二维数据来说，曲波正变换历以下五个步骤：

1)对数据进行二维傅立叶变换，将时间-空间域地震数据变换到频率-波数域；

2)对频率、波数域数据重采样或内插值；

3)将频率-波数域数据与窗口函数相乘；

4)绕原点Wrapping局部化；

5)进行二维傅立叶逆变换；得到Curvelet系数C(i,j,k)；

对于一个地震道集d(t,x)，所述Radon正变换、Radon反变换定义如下：

Radon正变换到τ-p域：

Radon反变换到t-x域：

把连续线性双曲线拉东变换积分变为有限区域内的相加，从而可以得到离散形式的线性拉东变换记作：

Radon正变换到τ-p域：

Radon反变换到t-x域：

其中：N_x----地震记录道数；

d_x---------------道间距；

dp-------拉东域采样间隔为；

Np---------拉东域记录道数；

对应的也可以写成矩阵形式：

Radon正变换到τ-p域：u＝Ld

Radon反变换到t-x域：d＝L^Hu

L和L^H为拉东变换算子，L为大小为(Np×Nt)×(Nx×Nt)的正变换算子，它将时域的抛物线变换到Radon域收敛为一个点；L^H为L的共轭转置算子，即反变换算子，将拉东域的一个点反变换到时间空间域；

为了提高拉东域的分辨率，在求解公式时，先定义拉东反变换，即使目标函数最小，对u求导得到最小平方解：u＝(LL^H+λ²I)^-1L^Hd；

其中λ是使求解稳定的范围在0.1-1的阻尼因子，I为单位矩阵；

在步骤2中，所述Radon正变换为对曲波域不同角度、不同尺度的的矩阵数据进行拉东变换，并分析噪声和有效信号在拉东域的对应关系及其特征；

所述阈值滤波是指设置选取合适的阈值，将大于该阈值的数据置零而小于阈值的数据保持不变。

2.如权利要求1所述的针对黄土塬地区的曲波域Radon变换噪声压制方法，其特征在于：

所述阈值滤波的具体步骤是：

初步选定一个阈值，将各矩阵中大于该阈值的数据置零而小于阈值的数据保持不变；对滤波结果进行检验；

当阈值滤波结果可以满足要求时，阈值滤波结束；

当阈值滤波结果不能满足要求时，改变选定的阈值的大小并按照重新设定的阈值完成对各矩阵新一轮的迭代计算，直到阈值滤波结果可以满足噪声压制要求时为止。

3.如权利要求2所述的针对黄土塬地区的曲波域Radon变换噪声压制方法，其特征在于：阈值滤波结果可以满足要求是指：计算相邻两次计算数据样点振幅差的平方和，如果当前计算平方和与上一次计算平方和差异小于1％,则退出迭代。