梁晨曦,张固澜,2,李 磊,罗 帆,罗一梁,李 勇,段 景,杜 皓

(1.西南石油大学地球科学与技术学院,四川成都610500;2.油气藏地质与开发国家重点实验室,四川成都610500;3.中国石油集团东方地球物理勘探有限责任公司物探技术研究中心,河北涿州072751)

在地震勘探过程中,随机噪声的存在影响地震数据的后续处理精度,因此压制噪声对后期地震资料的处理、反演和解释等具重要作用。地震资料去噪方法一般分为原始数据域方法[1-2]和变换域方法[3-6]。原始数据域方法主要包括中值滤波[7-8]、维纳滤波[9]、各向异性扩散滤波等[10];变换域方法主要包括离散余弦变换[11-12]、小波变换[13-14]、曲波变换等[15]。这些算法由于缺少对非局部信息的利用,因此在去除噪声的同时会损害部分有效信号。

DABOV等[16]提出的块匹配三维协同滤波(Block-matching and 3D collaborative Filtering,BM3D)方法将数据域与变换域结合,充分利用数据自相似性和冗余性信息对随机噪声进行压制,其核心步骤如下:①将原始数据分解成多个相同大小的分块数据[17-18];②对各分块数据进行二维变换,即先在列(行)方向进行一维变换,直至所有列(行)都处理结束,在此基础上,再在行(列)方向进行一维变换,直至所有行(列)都处理结束,得到二维变换域分块数据;③在二维变换域,按目标块和各分块数据的相似程度和给定的阈值进行块匹配[19-20](按照匹配后相似块的块顺序号进行),得到分组后的三维数据体;④对分组后的三维数据体在块序号方向进行一维变换得到三维变换域数据,并进行三维协同滤波(硬阈值滤波和经验维纳滤波)去噪处理;⑤对分块数据块聚集重构图像,得到最终的去噪结果。

BM3D方法在去噪[21-22]的同时可较好地保留原始数据的有效信息[23],因此被广泛用于图像处理和地震数据处理等[24-25]领域。刘向乐等[26]提出的变换域方法为小波变换的BM3D方法,在图像处理中取得了较好的应用效果。任婷婷等[27]提出的变换域方法为离散余弦变换(discrete cosine transform,DCT)的BM3D方法,可用于地震数据随机噪声压制。孙成禹等[28]将曲波变换与BM3D相结合,提出了基于曲波估计噪声的BM3D方法,并用于实际地震数据的随机噪声压制,该方法利用曲波变换估计地震数据噪声强度,并基于估计的噪声强度自适应调整BM3D中的滤波阈值和块匹配参数。张欢等[29]提出结合主成分分析的BM3D方法,该方法通过对地震数据进行降维,得到地震数据噪声估计值,并基于估计噪声值规避传统噪声预估值的局限性,优化了BM3D对地震数据的噪声估计。

BM3D方法的块匹配及分组结果受块间能量差异影响较大,去噪精度有待提升;另外,BM3D方法中的参数选取(如滤波阈值等)需经过多次试验才能确定,参数选取效率较低。为提高BM3D方法的去噪效果和参数选取效率,本文提出了基于分块能量归一化的BM3D,介绍了该方法的基本原理及具体实现步骤,并将其与传统方法(中值滤波、DCT滤波、BM3D)进行对比。最后通过模型数据和实际数据测试与应用,验证本文方法的有效性。

1 块匹配三维协同滤波基本原理

BM3D由基础估计和最终估计组成,其核心步骤如图1所示,主要分为以下3个方面:①硬阈值块匹配及分组;②三维协同滤波,即基础估计阶段的硬阈值滤波和最终估计阶段的经验维纳滤波;③块聚集。

图1 传统BM3D算法基本流程

1.1 硬阈值块匹配及分组

(1)

且

(2)

(3)

其中,{ }表示集合。

假设Di[m][n],Dj[m][n]和Gj[m][n]分别为分块数据Di,Dj和相似块匹配结果Gj的离散表示,则:

(4)

且

(5)

(6)

且

(7)

由(1)式、(2)式、(4)式和(5)式可知,τ1和τ2越小,搜索到的相似块越少;τ1和τ2越大,搜索到的相似块越多。由(3)式和(5)式可知,欧氏距离d(Bi,Bj)和d(Di,Dj)受分块数据总能量差异影响较大,并不能完全反映分块数据间的相似度。因此,在分块数据总能量差异较大的情况下,τ1和τ2的选择会影响BM3D的计算效率,需多次试验;且搜索到的相似块不能真正反映搜索块与目标块间的相似度,影响了三维协同滤波结果,最终影响BM3D去噪效果。

1.2 三维协同滤波

(8)

且

(9)

其中,|·|表示取模。

(10)

且

(11)

其中,δ∈[0,+∞)。

由(8)式至(11)式可知,λ1越大,噪声压制效果越明显;δ越大,噪声压制效果越明显。在分块数据总能量差异较大情况下,λ1和δ会影响BM3D计算效率和三维协同滤波结果,最终影响BM3D去噪效果,其值的选择需多次试验。

1.3 块聚集

H0[x][y]=

(12)

H1[x][y]=

(13)

2 块能量归一匹配三维协同滤波基本原理

BM3D方法易受块间能量差异影响且参数选取效率不高,相似块匹配精度、处理参数选取效率及滤波效果都有待提升。为消除上述影响,本文提出块能量归一BM3D方法,主要包括以下3个方面:①能量归一软/硬阈值块匹配及分组;②能量归一软/硬阈值三维协同滤波;③能量归一块聚集。其核心步骤如图2 所示(改进部分为图2中红色模块),与BM3D方法的不同之处主要体现在①和②两方面。

图2 改进的BM3D算法基本流程

2.1 能量归一软/硬阈值块匹配及分组

(14)

(15)

式中:max(·)表示取最大值;μ1为实数,且μ1∈[0,1];Ei表示分块数据Fi的能量。

变换域软/硬阈值相似块匹配可表示为:

(16)

且

(17)

式中:软阈值τ2为常数,且τ2∈[0,1];硬阈值τ3为常数,且τ3∈[0,1]。

(18)

2.2 能量归一软/硬阈值三维协同滤波

(19)

式中:μ2为实数,且μ2∈[0,1]。

2.3 能量归一块聚集

最终估计块聚集过程如下:①滤波数据分块;②原始图像重构。实现方式与传统方法一致。

3 模型数据试算

通过模型试算来对比分析传统BM3D方法与本文方法的相似块匹配精度和滤波效果。本文方法中,变换域方法均为二维DCT,图中色标数值均代表振幅值。

3.1 模型试算一

图3a为12×12个像素点的无噪理论模型数据的数字显示,其中,图3a右上角、左下角和右下角蓝框区域像素值分别为红框区域像素值的2倍、5倍和10倍;图3b为图3a的图像显示;图3c为对图3a加入随机噪声后的图像显示。

图3c中黄色方块为目标块,利用传统的块匹配方法和本文改进的块匹配方法,进行相似块匹配处理。图4为不同τ1的传统块能量归一硬阈值块匹配结果。图5为不同τ2的块能量归一软阈值块匹配结果。图6为不同τ3的块能量归一硬阈值块匹配结果。对比图3至图6可知,传统块匹配结果受块间能量差异影响,精度不高;且随着硬阈值的不断增大,其相似块匹配精度提升效果不明显;基于块能量归一的软/硬阈值块匹配均消除了块间能量差异影响,其相似块匹配精度较高。

图3 无噪理论模型的数字显示(a)、无噪理论模型的图像显示(b)和含噪理论模型的图像显示(c)

图4 不同τ1的传统块能量归一硬阈值块匹配结果a τ1=16;b τ1=200;c τ1=1000

图5 不同τ2的块能量归一软阈值块匹配结果a τ2=0.08;b τ2=0.20;c τ2=0.50

图6 不同τ3的块能量归一硬阈值块匹配结果a τ3=0.02;b τ3=0.05;c τ3=0.20

3.2 模型试算二

图7a为含有复杂断层的无噪合成地震数据,图7b 为对图7a加入随机噪声后的合成地震数据。图8a、图8b、图8c和图8d分别为对图7b使用二维中值滤波、二维DCT滤波、传统BM3D方法和本文方法后的滤波结果,滤除的噪声分别如图9a、图9b、图9c和图9d所示。图10为4种方法去噪后的峰值信噪比(peak signal to noise ratio,PSNR)[30]的对比结果。

由图7至图10可知:二维中值滤波和二维DCT滤波后的剖面断点不清晰,边缘保持能力不强且对有效信息伤害较大(如图8a、图8b、图9a和图9b中黑色箭头所示);传统的BM3D方法滤后的剖面断点较清晰,边缘保持能力较强,且对有效信号伤害较小(如图8c 中黑色箭头所示);本文方法较传统的BM3D方法的滤波效果有所提升,且滤除的噪声剖面中几乎不含有效信号(如图9c、图9d中黑色箭头所示)。对比4种方法的PSNR值可知,本文方法峰值信噪比最大,较传统方法去噪效果有明显提升。

图7 原始无噪合成地震数据(a)和含噪合成地震数据(b)

图8 采用不同滤波方法的滤波结果a 二维中值滤波;b 二维DCT滤波;c 传统BM3D方法;d 本文方法

图9 采用不同滤波方法滤除的噪声a 二维中值滤波;b 二维DCT滤波;c 传统的BM3D方法;d 本文方法

图10 4种方法去噪后的峰值信噪比对比

4 实际地震数据应用

采用传统BM3D滤波方法和本文改进的BM3D滤波方法对实际三维地震数据进行去噪处理,具体步骤如下:①对三维地震数据体按Inline(或Crossline)线读取数据,并对读取的Inline线剖面数据进行去噪处理;②重复步骤①,直至所有的Inline线(或Crossline线)地震数据去噪处理都结束,得到去噪后的三维地震数据体。

图11a、图11b分别为某三维地震数据中Inline100和Crossline100原始数据。图12a和图12b分别为对图11a采用传统BM3D方法和能量归一块匹配软阈值BM3D方法进行滤波处理后得到的剖面,滤除的噪声分别如图12c和图12d所示。图13a和图13b 分别为对图11b采用传统BM3D方法和能量归一块匹配软阈值BM3D方法进行滤波处理后得到的剖面,滤除的噪声分别如图13c和图13d所示。图14a 是某三维地震数据中1s等时切片,图14b和图14c 分别为对图14a采用传统BM3D方法和能量归一块匹配软阈值BM3D方法进行滤波处理后得到结果,滤除的噪声分别如图14d和图14e所示。

图11 Inline100(a)和Crossline100(b)原始数据

由图12、图13和图14可见:两种方法滤波处理后的剖面上同相轴连续性和信噪比均得到有效提升;改进的BM3D方法的滤波结果对构造刻画精度较传统BM3D方法的滤波结果有所提升;另外,传统BM3D方法滤除的噪声中含有较多有效信号,而改进的BM3D方法滤除的噪声中几乎不含有效信号,去噪后反射同相轴更清晰且参数选取效率较高。

图12 对Inline100原始数据采用传统BM3D方法和本文方法进行滤波处理后得到的剖面a 传统BM3D方法;b 本文方法;c 传统BM3D方法滤除的噪声;d 本文方法滤除的噪声

图13 对Crossline100原始数据采用传统BM3D方法和本文方法进行滤波处理后得到的剖面a 传统BM3D方法;b 本文方法;c 传统BM3D方法滤除的噪声;d 本文方法滤除的噪声

图14 1s处等时切片a 原始数据;b 传统BM3D方法滤波处理后的结果;c 本文方法滤波处理后的结果;d 传统BM3D方法滤除的噪声;e 本文方法滤除的噪声

5 结论

1)本文提出的基于分块数据总能量归一的块匹配及分组算法,可有效提升相似块匹配精度,从而提升BM3D的应用效果;

2)本文提出的基于分块数据总能量归一的软/硬阈值滤波方法,通过软阈值滤波减少参数试验次数,从而提高了参数选取效率;

3)本文方法提升了去噪效果和参数选取效率,可被广泛用于地震资料的去噪处理。