一种基于曲波变换-联合双边滤波的地震数据去噪方法

1.本发明涉及信号处理领域，特别是涉及一种传统曲波变换算法与联合双边滤波算法相结合的地震数据去噪方法。


背景技术：

2.随着我国地震勘探进程的不断发展，地震勘探项目涉及到的地区往往是复杂地区。在地震数据采集中，由于复杂地区的地表层激发条件不太理想，以及野外环境噪声的干扰，采集到的地震数据资料的分辨率和信噪比比较低，地震资料解释人员容易作出错误的判断，从而影响油层位置的勘探精度，增加石油开发成本，造成巨大经济损失。因此，地震资料去噪仍是地震资料处理分析的难点和重点之一。
3.现有的地震去噪的方法比较多，主要充分利用时频变换域内稀疏特点，对变换域中的系数进行阈值处理，只保留幅值较大的稀疏系数，幅值较小的稀疏系数将被置0，从而实现去噪的效果。但这种处理方法处理后，稀疏系数在阈值处不连续，导致去噪后在地震资料成像中出现部分弱同相轴模糊化，连续性降低，影响地震资料的保真度。


技术实现要素：

4.为了解决现有地震去噪方法的问题，本发明提供了一种基于曲波变换-联合双边滤波的地震数据去噪方法。通过改进曲波变换阈值去噪方法，实现在去噪的同时，保护弱同相轴的效果，在曲波域中，引入联合双边滤波算法，确保稀疏系数在阈值处连续，提高弱同相轴的连续性。
5.本发明所采用的具体技术方案是：首先，对含噪声的地震数据进行曲波变换，得到原始曲波系数。其次，对原始曲波系数进行尺度分析，将第一尺度曲波系数作为粗尺度信号，粗尺度信号主要包含主曲波系数，表征地震数据构造的特征；将最后两个尺度曲波系数作为细尺度信号，细尺度信号通常包含噪声信号；剩下的尺度曲波系数作为中间尺度信号，中尺度信号通常包含大量的弱同相轴信号。然后，由于粗尺度信号包含低频信息，而且同相轴基本连续。因此，对粗尺度信号使用双边滤波处理，得到新粗尺度曲波系数；中间尺度信号包含大量的弱同相轴信息，对不同尺度的中间尺度信号曲波系数设置不同的阈值，通过每个尺度的不同阈值筛选出最多25％较大幅值曲波系数并保留，得到阈值系数并作为引导曲波系数；将中间尺度信号曲波系数作为被处理的曲波系数，利用联合双边滤波算法；联合双边滤波处理后，将阈值系数和联合双边滤波后的曲波系数进行重组，得到新中间尺度曲波系数；细尺度信号主要包含噪声信号，对细尺度信号进行阈值处理，得到新细尺度曲波系数。最后，将新粗尺度曲波系数、新中间尺度曲波系数和新细尺度曲波系数进行重组，并用曲波逆变换处理得到去噪后的地震数据。
6.与现有方法相比，本发明中在曲波变换域中引入联合双边滤波算法，联合双边滤波具有引导图，结合中间尺度曲波系数及其相似性，对缺失的同相轴信息进行补全，能够提高弱同相轴的连续性。曲波变换阈值法能够较好的压制噪声。本发明在有效去除噪声的同
时，保护弱同相轴信号。
附图说明
7.图1为本发明的方法流程图；
8.图2为含噪地震数据的各个尺度分析的频谱图；
9.图3为基于曲波变换-联合双边滤波地震数据去噪方法的具体实施步骤流程图；
10.图4为实际地震数据噪声压制结果对比全局剖面图及其局部放大图；
11.图5为不同方法对实际地震数据噪声的压制结果对比剖面图及其局部放大图；
12.图6为不同方法对实际地震数据噪声的压制结果对比单道图；
13.图7为不同方法对实际地震数据噪声的压制结果对比频谱图。
具体实施方式
14.本发明公开了一种基于曲波变换-联合双边滤波的地震数据去噪方法，属于地震数据处理领域。本发明通过改进曲波变换阈值方法，针对曲波变换硬阈值法在阈值处不连续的问题，在曲波域中引入联合双边滤波，实现了地震数据去噪的同时，提高弱同相轴的连续性。下面通过介绍本发明具体实施步骤，结合实例和附图做更进一步的具体说明，根据以下说明，本发明的内容、特征、优点和效果等将会变得更加清楚。
15.图1为本发明的方法流程图，具体实施步骤如下：
16.步骤101：对含噪声的地震数据进行曲波正变换处理。实际采集到的地震数据包含有效信号与噪声，可用以下表达式描述地震数据的构成：
17.f(t)＝o(t) n(t)
18.其中，f(t)代表含噪声地震数据，o(t)代表无噪声的原始地震数据，n(t)代表噪声。定义母曲波其中，x为空间域因子，j为尺度因子，j＝1,2,
…
,r。总尺度数r由原始地震数据f(t)的地震道数m和采样点数n来确定，
[0019][0020]
其中，为向上取整运算。的傅里叶变换与傅里叶窗函数uj(ω)满足ω为频率域因子，所有曲波函数均可通过平移和旋转得到。引入位移因子序列k＝(k1,k2)∈z2，旋转角度和位置其中，0≤θ
l
＜2π，l为角度因子，且l＝0,1,
…
，为向下取整运算，则曲波函数可表示为：
[0021][0022]
式中，旋转矩阵
[0023][0024]
对于原始地震信号f(t)，在实数域中通过曲波正变换得到原始曲波系数c
j,l,k
，
[0025]
[0026]
式中，为f(t)的傅里叶变换。
[0027]
步骤102：对原始曲波系数进行尺度分析。对地震数据的每个尺度曲波系数分别进行曲波逆变换，得到各个尺度信息的频谱图如图2所示，可以看出，第一尺度只包含低频信息，表征地震数据构造的特征，归为粗尺度信号；第二尺度到倒数第三个尺度包含多个角度，这些角度包含了大量有效信号，归为中间尺度信号；最后两个尺度主要是由高频噪声信号组成，故归为细尺度信号。
[0028]
基于此分析，将原始曲波系数划分为粗尺度信号曲波系数q
j,l,k
，中间尺度信号成分的曲波系数s
j,l,k
和细尺度信号成分的曲波系数t
j,l,k
。
[0029][0030]
步骤103：对粗尺度信号进行双边滤波处理。粗尺度信号包含低频信息，而且同相轴基本连续，无需对缺失的同相轴进行补全。本文应用双边滤波处理粗尺度信号的曲波系数，弱同相轴曲波系数之间的相似度较大，噪声信号曲波系数之间的相似度较小。双边滤波根据曲波系数之间的相似度信息进行处理，保留弱同相轴曲波系数，缩小噪声信号曲波系数。对粗尺度信号的曲波系数q
j,l,k
进行双边滤波，得到新的曲波系数
[0031][0032]
其中，为双边滤波公式的归一化系数，
[0033][0034]
式中，k
q1
，k
q2
为粗尺度信号曲波系数q
j,l,k
在频率域空间的位移因子，kq＝(k
q3
,k
q4
)，且k
q3
，k
q4
为范围内粗尺度信号的曲波系数在频率域空间的位移因子；表示以(k
q1
,k
q2
)为中心，n1×
n1大小的双边滤波范围，n1为双边滤波边长且n1＝2a1 1，a1＝0,1,2
…
；
[0035]wd1
(k
q3
,k
q4
)为双边滤波的空间域权重，
[0036][0037]wr1
(k
q3
,k
q4
)为双边滤波的值域权重，
[0038][0039]
式中，为双边滤波算法几何距离的方差，为双边滤波灰度相似度的方差。
[0040]
步骤104：对中间尺度信号进行联合双边滤波处理。中间尺度信号曲波系数s
j,l,k
包含大量同相轴的信息，为提高曲波系数在阈值处的连续性和抑制噪声系数，对不同尺度的中间尺度信号曲波系数设置不同的阈值λj'，通过每个尺度的不同阈值筛选出最多25％较大幅值曲波系数并保留，得到阈值系数并作为引导曲波系数，
[0041][0042]
将中间尺度信号曲波系数s
j,l,k
作为被处理的曲波系数，利用联合双边滤波算法，结合曲波系数的相似性，对缺失的同相轴信息进行补全。得到联合双边滤波处理后的曲波系数
[0043][0044]
其中，为联合双边滤波公式的归一化系数，
[0045][0046]
式中，k
s1
，k
s2
为中间尺度信号曲波系数s
j,l,k
在频率域空间的位移因子，ks＝(k
s3
,k
s4
)，k
s3
，k
s4
为范围内的中间尺度信号曲波系数在频率域空间的位移因子；表示以(k
s1
,k
s2
)为中心，n2×
n2大小的联合双边滤波范围，n2为联合双边滤波的边长且n2＝2a2 1，a2＝0,1,2
…
；为联合双边滤波的空间域权重，
[0047][0048]
为联合双边滤波的值域权重，
[0049][0050]
式中，为联合双边滤波算法几何距离的方差，为联合双边滤波灰度相似度的方差。和分别为(k
s1
,k
s2
)点和(k
s3
,k
s4
)点对应的曲波系数。联合双边滤波处理后，将大曲波系数(阈值系数)和联合双边滤波后的曲波系数进行重组，得到去噪后的中间尺度信号曲波系数
[0051][0052]
步骤105：对细尺度信号进行硬阈值处理。细尺度中，绝大部分信号为噪声信号，因此对细尺度信号只用硬阈值处理即可，硬阈值公式如下所示，
[0053][0054]
步骤106：对新粗尺度信号，新中间尺度信号和新细尺度信号进行重组，得到新曲波系数。利用双边滤波、联合双边滤波和硬阈值处理对不同尺度的曲波系数进行去噪处理之后，构造去噪后的曲波系数
[0055][0056]
步骤107：对新曲波系数进行曲波逆变换，得到去噪后的曲波系数，曲波逆变换公式如下所示，
[0057][0058]
最后，本发明通过频谱图对去噪后的结果进行评价。
[0059]
下面结合具体实施例对本发明的实施过程进行演示，图3为本发明提供的具体实施步骤。
[0060]
步骤1：本示例采用实际地震数据cmp_inl2000_stk_org_gain.segy。图4为实际地震数据和本发明去噪后的地震数据剖面图。其中，(1-1)和(1-2)分别为全部的实际地震数据和红色方框中的实际地震数据。
[0061]
步骤2：将数据带入曲波变换-联合双边滤波算法编程。
[0062]
步骤3：将地震数据分为3个尺度部分。
[0063]
步骤4：对粗尺度信号用双边滤波处理。
[0064]
步骤5：先设置阈值系数并作为引导图，再结合中间尺度信号，使用联合双边滤波进行处理。
[0065]
步骤6：对细尺度信号用硬阈值去噪。
[0066]
步骤7：对分别用双边滤波、联合双边滤波和硬阈值处理后的结果进行重组，再用曲波逆变换得到去噪后的地震数据。
[0067]
步骤8：对比本发明和小波变换、曲波变换的去噪结果。图5为各个方法去噪后的地震数据剖面图。其中，(1-1)和(1-2)分别为小波变换去噪后的全部地震数据及其红色方框显示的局部放大图；(2-1)和(2-2)分别为曲波变换去噪后的全部地震数据及其红色方框显示的局部放大图；(3-1)和(3-2)分别为本发明去噪后的全部地震数据及其红色方框显示的局部放大图。从图5可以看出，采用小波变换，曲波变换后，弱信号被剔除，分辨率降低。本发明在去除噪声的同时，弱信号被保留和增强，尤其是深层弱信号信息丰富，为深层资料后续的建模和解释提供了高质量的地震资料。局部放大图和单道对比都表明了本文方法的有效性。其中，单道对比图如图6所示，(1-1)、(1-2)、(1-3)和(1-4)分别为实际含噪地震数据单道图、小波变换去噪后的地震数据单道图、曲波变换去噪后的地震数据单道图和本发明去噪后的地震数据单道图。
[0068]
步骤9：求整个含噪实际地震数据和去噪后的地震数据的频谱图。图7为各个方法去噪后的频谱图。从图7可以看出，对于高频部分，三种方法均能较好地压制噪声。对于低频部分，小波变换处理的结果中包含的能量信息较少，显然，大部分同相轴信号被压制；曲波变换处理结果优于小波变换，但与实际资料相比，仍有部分弱同相轴信号被压制；本文所提方法处理的结果明显优于小波和曲波的处理结果。
[0069]
通过上述过程可知，本发明提供了一种基于曲波变换-联合双边滤波地震去噪方法，解决了硬阈值法在阈值处不连续的问题，在去噪的同时，提高了弱同相轴的连续性，显著提高信噪比，其效果优于传统方法。