基于曲波噪声估计的三维块匹配地震资料去噪

孙成禹 刁俊才 李文静

(①中国石油大学(华东)地球科学与技术学院，山东青岛 266580； ②青岛海洋国家实验室海洋矿产资源评价与探测技术功能实验室，山东青岛 266071； ③东方地球物理公司物探技术研究中心，河北涿州 072751)

0 引言

地震数据中噪声严重影响资料处理和解释结果[1]，压制噪声对后期地震资料的处理、反演和解释等具重要作用[2]。地震资料常用的去噪方法一般分为空间域方法和变换域方法两大类。空间域去噪方法包括均值滤波、中值滤波[3]和维纳滤波[4]等，变换域去噪方法包括傅里叶变换、小波变换、曲波变换和Randon变换等阈值方法[5-6]。这些去噪算法本质都是利用了图像本身的局部相关性，没有充分挖掘、利用图像的非局部相关性，因此在压制噪声的同时也会损害部分有效信号，不能精确保持断层、裂缝等地质体的边缘特征，甚至会出现伪影等现象[7]。

针对这一问题，Buades等[8]提出利用原始数据的结构相似性进行非局部去噪，与传统的局部去噪算法相比，该算法能有效地保护图像的边缘细节信息。在非局部去噪算法中，Kostadin等[9]提出的三维块匹配(Block-Matching 3D，BM3D)去噪算法利用块匹配和三维变换域滤波技术进行串联去噪，该方法综合了非局部算法和变换域算法的优势[10-11]，对随机噪声具有很好的去噪效果。BM3D去噪是将图像分割成许多个不同的小块，根据不同小块之间的相似性进行块匹配，然后在三维变换域中去除噪声，充分利用了数据的自相似性和冗余性信息，能较好地保留信号细节[12]。韩玉兰等[13]针对BM3D算法运算量大、运行时间过长的问题，采用积分图计算块的相似性，代替相应的滤波过程，提高了计算效率；李文静等[14]将改进后的BM3D算法引入地震数据处理领域，得到了较为理想的去噪效果。但该方法依然存在一些不足，滤波阈值参数需要根据经验人为设定[15]，而阈值选取的不确定性严重影响基础估计部分的降噪效果。由于在实际计算中，噪声方差对阈值大小和块匹配相似度的影响最大，而实际数据噪声的方差却是未知的[16]，极大地限制了该算法在地震资料去噪中的应用。

在变换域地震资料去噪方法中，Candès等[17]在脊波变换的基础上提出了第一代曲波变换。张恒磊等[18]利用曲波尺度分解得到较为准确的噪声方差估计值，使后续的曲波变换阈值选取更加精确。但是由于实际地震数据的复杂性，单一曲波变换方法的去噪精度不足，处理结果常常存在伪影和过度平滑等现象。针对这些问题，薛诗桂等[19]将曲波变换和循环平移技术结合，消除了曲波变换产生的伪吉布斯效应。姚振岸等[20]将各向异性扩散滤波和曲波变换结合保护地震数据的边界特征；杨会等[21]将曲波变换和二维经验模态分解结合，在去除噪声的同时较好地保护了弱信号。曹静杰等[22]将曲波变换作为稀疏变换自适应地进行稀疏反演去噪。

为了弥补曲波变换和常规BM3D去噪方法的缺陷，本文结合两种方法的优势，提出了一种基于曲波噪声估计的BM3D地震资料去噪方法。该方法无需人为设定相关参数，根据曲波变换得到的地震资料噪声估计作为先验信息，准确计算块匹配相似度，自适应选取合适的滤波阈值进行去噪，可以得到更好的效果。数值模型和实际资料测试表明，与常规BM3D去噪方法、曲波变换去噪法相比，该方法去噪效果明显、实用性较高。

1 基本原理

1.1 三维块匹配去噪原理

BM3D去噪是图像处理中的一种方法。其基本算法是：将含噪数据看作一个图像，将其划分为若干块，并根据各块与其他块的相似度进行匹配分组；对组内数据进行三维线性变换，在变换域内对噪声进行滤波后反变换回原来的域内，并将每个组内的块像素返回到图像原来的位置；最后对重叠的块进行加权平均，得到去噪后的结果。为了提高匹配分组的正确性，通常对含噪数据进行前置硬约束阈值滤波作为地震信号的基础估计，获取维纳滤波的收缩系数。因此，整个过程就包含基础估计和重新估计两个阶段。每个阶段也都包含块匹配、变换域去噪和聚集三个环节。

在基础估计的块匹配阶段，将含噪地震数据视为二维图像，划分成尺寸为N×N的多个块Zi(i=1，2，…，n，n为块总数)，依次选取各块作为参考块Zr(r=1，2，…，n)，定义块与参考块之间的距离为

(1)

(2)

(3)

(4)

(5)

(6)

将维纳滤波应用于噪声数据矩阵Yr，即将噪声数据矩阵三维变换后的系数与维纳滤波收缩系数相乘，再进行反变换后得到新的估计值为

(7)

最后将不同群组中的重叠块进行加权平均得到重新估计值

(8)

1.2 基于曲波噪声估计的BM3D地震资料去噪方法

常规的BM3D去噪方法需要根据噪声方差σ2作为先验信息求取滤波阈值和块匹配参数。在实际地震资料的处理中，常用的噪声方差估计方法是中值噪声方差估计和基于小波分析的细节估计。当原始地震资料信噪比较高、噪声方差较小时，这两种方法的估计值与真值较为接近。但是，当原始地震资料的噪声水平较高时，前者的估计值与真值偏差较大，而后者则不能完整地估计噪声，导致估计值比实际噪声方差偏小。

(9)

根据σ可以进一步确定BM3D去噪算法中的滤波阈值γ1及收缩系数φr。这样滤波参数成为噪声强度σ的自适应函数，弱反射和边界细节信息可以得到有效的保护。

BM3D去噪算法结合了非局部算法和变换域算法的优势，可以有效去除高斯随机噪声。由于事先无法得到噪声方差，需要人为设定滤波参数，因而去噪效果难以保证。本文将曲波变换估计出的噪声方差作为先验信息，求取适用于BM3D去噪算法的输入噪声强度σ，进而自适应地调整滤波阈值和块匹配参数，完成去噪处理。当σ较小时，可直接应用式(1)计算块匹配相似度；当σ为中等或较大时，直接块匹配可能会导致错误的块匹配分组，需要先对图像块进行预滤波处理，即先对两个图像块进行正则二维线性变换，对其进行阈值滤波，然后再计算两者之间的距离，即

(10)

式中：γ0=λ2Dσ为二维预滤波处理的阈值(λ2D为二维系数)；T2D为正则二维线性变换算子。去噪阈值参数通常都设为常量，为了有效地保护有效信号，将其设定为噪声强度σ的自适应函数。当σ较小时，无需进行预滤波处理，图像块可以取值较小(如N取为8)； 当σ中等时，需要进行预滤波处理，滤波阈值γ0取为1.0，图像块应适当增大(如N取为9～10)；当原始资料噪声强度σ较大时，二维预滤波阈值γ0可取为2.0， 图像块应取值较大(如N取为11～12)。

2 模型试算

为了测试自适应选取阈值的改进方法的优越性，建立楔状体模型，对其正演数据加入随机噪声，得到含噪记录如图1a所示。分别利用常规BM3D方法和改进的BM3D法进行去噪处理，结果如图1b和图1c所示。

图1 楔状体模型正演数据两种方法去噪结果对比(a)原始含噪数据； (b)常规BM3D去噪结果； (c)本文方法去噪结果

选用峰值信噪比和信噪比两个参数对去噪效果进行定量比较。二者分别定义为

(11)

(12)

对图1a所示的模型正演数据，变化所加入噪声的强度，分别使用常规BM3D和改进BM3D方法进行去噪处理，并计算其PSNR和SNR，结果如图2

所示。从图中曲线可以看出，改进BM3D去噪方法相对常规BM3D方法去噪后SNR和PSNR均更高，尤其是低信噪比情况下去噪效果更显著。

为了进一步说明本文方法相对于常规方法的优越型，使用图3a所示的模型数据进行测试。对模型数据加入随机噪声，得到信噪比为3dB(即信号与噪声的能量之比为2)的含噪记录，如图3b所示。分别使用常规BM3D方法、曲波变换滤波方法和本文方法进行去噪处理，结果如图3c、图3e和图3g所示，去噪前后的差值剖面即去除的噪声如图3d、图3f和图3h所示。由图可以看出，加入随机噪声后模型断点处的反射基本淹没在噪声中，断层边缘位置模糊不清。 常规BM3D方法能够在一定程度 上去除部分噪声，但去噪后同相轴连续性降低(图3c)；曲波变换滤波后的剖面(图3e)同相轴弯曲部分连续性降低，断层边缘有效信号在残差剖面中明显(图3f中红框处)，对断层细节信息保持较差。本文方法不仅能明显去除随机噪声，且去噪后反射同相轴清晰，断层边界特征也得到了良好的保持(图3g)。从图3h所示的去除噪声剖面上可以看出，本方法未将有效信号作为噪声去除，很好地保留了有效信息，去噪效果相对最佳。图4为三种方法去噪后的SNR和PSNR对比，可以看出本文方法的去噪效果最好。

图2 两种方法去噪效果的定量对比(a)峰值信噪比； (b)信噪比

图3 不同方法去噪效果对比(a)原始模型数据； (b)加入随机噪声的模型数据； (c)常规BM3D方法去噪结果； (d)常规BM3D方法去除的噪声； (e)曲波变换方法去噪结果； (f)曲波变换方法去除的噪声； (g)本文方法去噪结果； (h)本文方法去除的噪声

图4 不同方法去噪结果的定量对比(a)SNR； (b)PSNR A：常规BM3D方法； B：曲波变换法； C：本文方法

3 实际资料试算

图5a为实际叠后地震剖面，其中含有随机噪声，影响了资料品质。分别使用常规BM3D法、曲波变换法和本文方法进行去噪，结果如图5b、图5c和图5d所示。

图5 不同方法实际资料去噪结果对比

(a)原始剖面； (b)常规BM3D方法； (c)曲波变换法； (d)本文方法

为了说明本方法对断层边界信息的保持能力，从图5中截取存在多个断层的红框区域进行放大显示，如图6a～图6d所示。可见：常规BM3D去噪结果边界信息不够清楚，弱反射层同相轴连续性较差(图6b)；曲波变换方法保持边界能力较差，断层边缘模糊(图6c)；本文去噪方法在去噪后断点清晰，断层解释更准确，较弱的反射同相轴连续性强(图6d)。图6e～图6g为三种方法去除的噪声剖面，可见常规BM3D法和本文方法去除的噪声剖面均匀，表明未去除有效信号；而曲波变换去除的噪声剖面中部分断点处存在有效信号(红框所示)。综合去噪后剖面和残差剖面可以看出，本文的方法去噪效果最好。

图6 断层区局部去噪效果对比放大显示(a)实际资料去噪前； (b)常规BM3D法去噪结果； (c)曲波变换法去噪结果； (d)本文方法去 噪结果； (e)常规BM3D法去除的噪声； (f)曲波变换法去除的噪声； (g)本文方法去除的噪声

4 结论

BM3D方法在地震资料去噪中具有良好的效果，但因无法准确估计噪声方差并选取合适的处理参数，在实际应用中受到一定的限制。本文联合曲波噪声估计和BM3D去噪方法，通过曲波分析，实现了滤波阈值和块匹配参数的自适应选取和去噪处理，提高了去噪精度和计算效率，具有较强的实际意义。理论模型及实际资料处理结果表明：

(1)基于曲波噪声估计的BM3D去噪方法能准确选取滤波参数，去噪后剖面的信噪比和峰值信噪比更高；

(2)本文方法避免了参数选取的多次测试，提高了计算效率；

(3)由于实现了滤波阈值的自适应选取，本文方法较好地保护弱反射信号和边界细节信息，为后续地震资料的精细解释提供保障。