基于多级中值的双边滤波算法的地震数据去噪研究

2022-04-24魏亚杰朱煜嘉曹静杰

河北地质大学学报 2022年2期

魏亚杰,朱煜嘉,曹静杰

河北地质大学 a.京津冀城市群地下空间智能探测与装备重点实验室、b.河北省战略性关键矿产资源重点实验室,河北石家庄 050031

0 引言

在地震勘探中,地震数据的精确处理可以直接影响后续的反演和偏移成像精度。地震数据去噪是提高地震信号信噪比(SNR) 一个不可缺少的环节,在地震AVO 分析、地震属性分析和微地震监测等许多地震勘探技术中采用高信噪比地震数据是至关重要的。在野外地震数据采集过程中,通常采用扩宽信号频带宽度来获取更加丰富的地震反射波信息,这样,在获取有效波的同时,不可避免地会记录各种噪声干扰,从而降低地震数据的信噪比,因此压制随机噪声是提高地震数据信噪比的关键。

随机噪声压制方法有很多,传统的去噪方法主要有空间去噪算法、变换域去噪算法和综合去噪算法。空间去噪算法常用的有中值滤波 (Huo 等,2012;Chen 等,2015; Gan 等,2016a)[1-3]、f-k 滤波(Spitz等,2008; Kim 等,2009)[4,5]等。一系列稀疏变换算法也获得了很好的应用,包括Fourier 变换(Liu 和Sacchi,2004; Xu 等,2005)[6,7]、不同种类的Radon变换(Verschuur 等,2012)[8]、Curvelet 变换 (Cao等,2019; 2020)[9,10]、Shearlet 变换(Guo 和 Labate,2007; 2013)[11,12],此外还有seislet 变换(Gan 等,2016b)[13]、focal 变换 (Kontakis 等,2016)[14]等。Bonar 和Sacchi 提出了一种非局部均值去噪方法,该方法可以在地震数据中找到相似的构造区域,然后进行分区去噪(Bonar 和 Sacchi,2012)[15]。Canales 提出了一种f-x 反褶积去噪算法,目前应用很广泛,已成为标准的地震信号去噪算法(Canales,1984)[16]。之后,Gulunay (1986)[17]和Hornbostel (1991)[18]进一步优化了f-x 反卷积算法,改进了随机噪声压制效果。Fieire 和Ulrych 根据低秩假设进行奇异值分解,提出了一种基于小波变换的去噪算法 (Freire 和Ulrych,1998)[19]。Cao 和Zhao (2016) 基于低冗余曲波变换提出一种地震数据同时去噪和重建方法,提高了原始曲波变换计算效率(Cao 和Zhao,2016)[20]。近年来,研究人员开发了一种基于字典学习的自适应变换算法(Zhu 等,2017)[21],该算法充分利用了地震数据本身的特点进行去噪,具有更好的去噪效果。空间域去噪算法具有更快的计算效率,变换域去噪算法去噪效果更好,综合类去噪算法结合了空间域去噪算法和变换域去噪算法的特点,以期达到更好的去噪效果。

目前比较流行的基于深度学习方法也被成功的应用于地震数据去噪,学者们提出了各种有效的网络结构,例如VGGNet (Simonyan 和Zisserman,2014)[22],U-Net (Ronneberger,2015)[23],FCN (Shelhamer 等,2017)[24]、CNN (Lecun 等,1998)[25]等,这些网络结构已广泛应用于图像去噪、人脸识别等。Lecun 等人证明CNN 的参数较少,在MNIST 数据集上提供了优异的分类结果(Lecun 等,1998)[25]。Zhang 等人提出了一个17 层的CNN,命名为DnCNN,用于一系列的图像噪声衰减,取得了良好的效果 (Zhang 等,2017)[26]。在地震数据去噪上,Gao 和Zhang 提出了一种地震数据同时重建和去噪的深度学习方法,取得了很好的效果(Gao 和 Zhang,2019)[27]。Song 等通过深卷积自编码神经网络压制地震随机噪声(Song等,2020)[28]。

Tomasi 提出的双边滤波技术是一种有效的空间域非线性滤波技术(Tomasi 和 Manduchi,1998)[29],该方法能够压制图像上的高斯噪声,通过计算局部加权平均,利用被平滑点所在邻域内像素点的几何距离信息和灰度值信息,能够最大限度保护图像的细节信息。但是该方法针对图像上的突变点噪声(椒盐噪声) 压制效果不理想,文章通过对双边滤波方法的深入研究,结合多级中值滤波方法改进双边滤波核函数,提出一种基于多级中值的双边滤波算法,能够同时去除地震数据中的高斯噪声和椒盐噪声。

1 方法原理

1.1 双边滤波器原理

双边滤波方法是一种空间域滤波方法,逐点进行滤波运算,首先以被修复点为中心划定一个邻域,通过分别计算邻域内的点对中心点的贡献,即为双边滤波修复后的点的值。假设u = u(x) 是被噪声污染的地震记录,其中x 为待计算点的位置,令Ω =Ωx(N),是以x 点为中心,N 为半滤波器宽的邻域。令y ∈Ω,这里用uy来表示y 点对应的灰度值,那么双边滤波修复后的点u～

x 可以表示为:

其中ω(x,y) 为双边滤波核函数为

方程(3) 中的ωS我们称之为距离权重,表示在邻域Ω 内的点距离,所计算点x 越远,所对应的距离权重值ωS越小,即对u～x贡献越小;方程(4) 中的ωR表示灰度权重,即邻域Ω 内点的灰度值与所计算点灰度值u(x) 相差越大,所对应的灰度权重值ωR越小,即这个点对u～x 贡献越小,σR、σS为常数,根据噪声水平,取值范围为(0,1)。双边滤波器能够压制图像上的高斯噪声,无法去除椒盐噪声。

1.2 二维多级中值滤波器原理

沿用上一节的邻域,待计算点x 的位置坐标为(a,b) ,令-N ≤k ≤N,可以定义一个基本子窗口

其中W(x) 表示以x 为中心的数值序列。常规二维多级中值滤波器基本子窗口如图1 所示,图1 中的圆点表示图像矩阵u 中某一点x 的一个邻域Ω =Ωx(2) 的像素位置,图1 中(a)、(b)、(c)、(d)中红色圆点分别表示子窗口W1、W2、W3、W4所取的像素位置。

多级中值滤波器输出定义为:

其中:

Ymax(x) = max

1≤i≤4[Zi(x)]max

Ymin(x) = min

1≤i≤4[Zi(x)]min

Zi(x) = median[x(,) ∈Wi(x)],i = 1,2,3,4

图1 二维多级中值滤波器基本子窗口N=2Fig.1 Basic sub window of two-dimensional multilevel median filter N=2

其中median 表示一般中值滤波。多级中值滤波具有方向性,保证了它相比一般中值滤波保护细节信息能力更强。

1.3 基于多级中值的双边滤波器原理

双边滤波器对高斯噪声有很好的压制效果,但是不能去除图像上的椒盐噪声,相反,多级中值滤波器能够压制图像上的椒盐噪声,对高斯噪声压制不明显。实际地震数据中同时存在高斯噪声和椒盐噪声,想要对其进行噪声压制就需要进行多步去噪处理,每一次去噪过程都会造成有效信号的缺失,不利于保留数据中的有效信号。

文章结合多级中值滤波器和双边滤波器提出了一种基于多级中值的双边滤波算法,将双边滤波器核函数灰度权重即方程(4) 中的ux替换为以x 点为中心,N 为半滤波器宽的邻域的多级中值滤波的值。

即多级中值双边滤波器的灰度权重为:

可以看到若所选邻域内存在椒盐噪声,且y 点为椒盐噪声点那么| Ymlm(x) - uy| 取较大值,那么该点灰度权重ωR较小,即这个点对u～x 贡献较小;若所选邻域内不存在椒盐噪声,那么计算领域内点的灰度权重值与双边滤波方法计算的灰度权重值一致,此时该方法等价于双边滤波方法。因此,文章提出的多级中值双边滤波器能够同时有效去除地震记录上的高斯噪声和椒盐噪声。

2 数值实验

2.1 模拟数据实验

首先采用简单的层状模型模拟数据来检验本文方法,模拟的原始不含噪数据如图2a 所示,采用中间激发的方式,101 道检波器接收,时间采样点为1001,采样间隔为4 ms。图2b 为原始数据中加入均值为0 方差为0.01 的高斯噪声,图2c 为在图2b 的基础上加入10%的椒盐噪声。图3 分别为图2 中数据对应的f-k 谱。

下面分别采用多级中值滤波方法、双边滤波方法以及本文提出的方法对含高斯噪声和椒盐噪声数据图2c 进行噪声压制。图4a 为采用多级中值滤波方法的噪声压制结果;图4b 为采用双边滤波方法的噪声压制结果;图4c 为本文方法的噪声压制结果;图5a 为多级中值滤波方法去除的噪声;图5b 为双边滤波方法去除的噪声;图5c 为本文方法去除的噪声。图6a、图6b 和图6c 分别为图4a、图4b 和图4c 对应的f-k谱。可以看到多级中值滤波方法只对椒盐噪声有效果,无法压制高斯噪声,双边滤波方法对高斯噪声有很好的压制效果,无法压制椒盐噪声,本文方法可以同时去除高斯噪声和椒盐噪声。

图2 a-原始不含噪数据;b-含高斯噪声数据;c-含高斯噪声和椒盐噪声数据Fig.2 a-Original data without noise; b-seismic data with Gaussian noise; c-seismic data with Gaussian noise and salt & pepper noise

图3 a-图2a 的f-k 谱;b-图2b 的f-k 谱;c-图2c 的f-k 谱Fig.3 a- f-k spectrum of Fig.2a; b- f-k spectrum of Fig.2b; c- f-k spectrum of Fig.2c

图4 含高斯噪声和椒盐噪声数据压制结果Fig.4 Suppression results of data with Gaussian noise and salt & pepper noise

为了更有效说明问题,进行了单道数据对比,这里选取第50 道数据,如图7 所示,可以很明显看到含噪声数据(蓝线) 中同时存在高斯噪声和椒盐噪声,多级中值滤波结果中残留有高斯噪声,双边滤波结果中残留有椒盐噪声,本文方法结果(红线) 与原始不含噪数据(黑线) 的匹配效果最好。

图5 差剖面a-图2c 与图4a 的差;b-图2c 与图4b 的差;c-图2c 与图4c 的差Fig.5 The difference profiles

图6 a-图4a 的f-k 谱;b-图4b 的f-k 谱;c-图4c 的f-k 谱Fig.6 a- f-k spectrum of Fig.4a; b- f-k spectrum of Fig.4b; c- f-k spectrum of Fig.4c

图7 第50 道数据对比Fig.7 Data comparison of the 50th channel

2.2 实际数据实验

下面采用本文方法对一个二维浅层地震剖面进行去噪,二维地震剖面时间采样点为2 800,采样率为0.1 ms,共240 道,道间距为25 m。二维地震剖面如图8 所示,可以看到该地震剖面上含有明显的噪声,造成部分同相轴不清晰、不连续,图9 为应用本文方法的对该二维地震剖面噪声压制的结果。图10 为本文方法去除的噪声,即图8 和图9 的差,图11 分别为图8 和图9 所对应的f-k 谱,可以看到随机噪声被有效压制。