基于噪声水平估计的加权核范数最小化噪声压制方法研究

唐 杰,张文征,戚瑞轩,李 聪

(中国石油大学(华东)地球科学与技术学院,山东青岛266580)

因各种干扰的存在,地震资料中会存在随机噪声,有些衰减噪声的滤波处理方法会造成地震反射中的不连续性信息被平滑而变得模糊[1-2]。为了克服过度平滑的问题,可以根据地震数据中的噪声水平采用自适应的滤波处理方法。噪声水平估计的常见方法有:基于小波分解的方法、基于主成分分析的方法和基于尺度不变性的方法[3]。基于小波分解的噪声估计方法根据小波变换的特点来估计噪声的标准方差,小波变换后数据的能量主要集中在尺度大的子带,而尺度小的高频子带其系数较小、能量较低。当噪声比例较高时,可将最高频率子带的系数全部看成噪声,由此来估计噪声的标准方差。基于主成分分析的噪声估计方法利用主成分分析(PCA)提取数据中特征值最小的区域,估计噪声水平[4-5]。这两种方法得到的估计结果低于真实噪声水平。而基于尺度不变性的估计方法利用峰度值不随尺度改变、只由噪声引起其系统变化的特性进行噪声估计,效果要明显优于前两种方法[6],尤其对边缘信息丰富的低噪声地震数据进行估计时效果更好[7-8]。

传统压制随机噪声的方法可分为空间域和变换域两种。空间域方法包括:均值滤波、中值滤波及各向异性扩散滤波等[9-10]。变换域方法包括:傅里叶变换滤波、小波变换滤波、曲波变换阈值去噪等[11-12]。均值滤波是滤除高斯噪声的常规算法,均值滤波算法将窗口内所有数据点的均值作为中心数据的值,对所有数据都采用无差别滤波,这样会使得非噪声数据点的值也被改变,容易破坏边缘细节信息,造成数据模糊。独立分量分析技术(ICA)是实现盲源分离的有效手段,它考虑数据的高阶统计特性,从而更加全面地表达数据的本质特征[13]。K次迭代奇异值分解(K-SVD)数据降噪算法是一种基于K-means算法和奇异值分解(SVD)的方法[14],该方法降噪计算时间较长,且对高比例的噪声数据降噪效果略显不足[15]。为了解决这些问题,DONG等[16]提出了非局部集中稀疏表示(NCSR)数据修复算法,对信噪比较低的地震数据的降噪效果有了较大提升,但是它构建用于稀疏表示所需字典时的计算量比较大,此外还可能过度地平滑数据,造成数据失真。为了克服这种方法的缺点,GU等[17]提出了加权核范数最小化(WNNM)降噪算法,该算法能根据矩阵奇异值刻画数据差异,给定不同的权值,突显数据中重要的信息,具有良好的去噪效果。在此基础上发展起来的多道加权核范数最小化噪声压制算法考虑多道之间的噪声差异,利用各道之间的冗余信息,同时考虑区分不同噪声的统计性信息,去噪效果较为理想[18-19]。

本文研究了基于噪声水平估计的加权核范数最小化噪声压制方法,在基于尺度不变性的噪声方差估计方法的基础上,采用加权核范数最小化方法去除地震数据中的随机噪声,采用理论模型数据和实际地震资料验证了方法的有效性。

1 方法原理

基于噪声水平估计的加权核范数最小化噪声压制方法采用分块处理的思路。首先利用块匹配的思想对地震数据进行分块处理,然后利用基于尺度不变性的噪声水平估计方法来计算分块数据中的随机噪声方差,利用所得到的方差来归一化加权核范数最小化算法中F范数的保真项,保证加权核范数最小化算法能在压制随机噪声的同时保护有效信号,从而实现根据地震数据本身噪声水平自适应压制随机噪声的目的。

1.1 基于尺度不变性的噪声水平估计

为了实现自适应的去噪处理,需要知道地震数据中随机噪声的分布特征,本文采用噪声方差估计来衡量数据中的噪声水平[20-21]。

用y表示某个含噪数据,x表示与y对应的不含噪声的原始数据,η表示噪声,则:y=x+η。假定噪声高斯分布,x可以用广义高斯分布拟合,x的广义高斯分布的峰度值κx(α)直接取决于其形状参数α:

式中:Γ是标准伽马函数;κx(α)表示数据x的峰度值。可以看出,峰度值κ与形状参数α成反比。在实际地震剖面的噪声估计中,峰度值κ为用方差平方归一化处理的四阶中心矩:

式中:μ4表示四阶中心距;σ2为方差。由于噪声的独立性,含噪数据y的方差满足:

y的四阶中心矩满足:

利用前面计算的方差的平方进行归一化处理:

根据上述结果,可以在给定无噪数据的情况下,从含噪数据中获得边缘滤波响应分布的峰度值。但是,在实际处理中,无噪数据是未知的,因此假设原始无噪数据具有尺度不变特性,其边缘滤波器响应数据分布的峰度值是未知常数,而在数据中加入噪声会导致整个数据的峰度值发生变化。用以离散余弦变换(DCT)为基的滤波器对图像作卷积处理,基函数为:

具体步骤如下:

1) 根据输入的含噪数据y和规定的分块大小N,由公式(6)得到N×N个以离散余弦变换为基的滤波器矩阵b;

2) 用每个滤波器矩阵b对图像进行卷积,生成响应数据yi;

图1给出了采用不同方法得到的随机噪声水平估计结果,噪声为白噪。从图1可以看出,基于尺度不变性的噪声估计方法的估计水平要优于小波分解估计和主成分分析方法,后二者估计结果与实际结果相比偏低。可以明显看出,在添加中、低比例噪声时,基于尺度不变性的噪声估计方法具有较好且稳定的估计结果,虽然添加高比例噪声时其估计效果不如低比例噪声时理想,但是估计结果仍然优于其它两种算法。因此基于尺度不变性的噪声估计方法在估计复杂低比例噪声时,效果较好。

图1 不同方法的噪声水平估计结果

1.2 多道加权核范数最小化

对于叠后含噪数据y的局部分块yj,可以采用块匹配的方法在数据中找到它的非局部分块[22]。将地震数据分成指定大小的数据块,选取yj作为参考块,其它为候选块,从候选块中找出与参考块相似的块,将这些非局部相似分块组成矩阵Yj,其中两个数据块之间的相似性采用(8)式表示:

假设Xj和Nj分别是无噪数据分块矩阵和噪声数据分块矩阵,则Yj=Xj+Nj,Xj为低秩矩阵,低秩矩阵逼近方式可以用于由Yj估计Xj。聚集所有经过降噪的分块,则可以估计整个数据[23]。

式中:‖·‖w,*表示加权矩阵奇异值的和,该方法的主要问题是确定权重向量w。权重向量与对应的奇异值成反比,权重公式为:

式中:λi(Xj)是数据块组Xj对应的第i个奇异值;c为常数;m是Yj中相似块的数量。为了避免分母为0,令ε=10-16。λi(Xj)的奇异值不能直接获得,可以近似为:

其中,λi(Yj)是Yj的第i个奇异值。

根据实际数据选择迭代次数L和迭代权重参数δ,L的设置原则为处理结果趋于稳定,根据试验结果分析,当L=6,δ=0.1时,去噪效果较好。

2 理论模型试算

2.1 模型数据去噪效果分析

建立如图2a所示的合成记录剖面,图2b为含噪数据,图2c为采用本文方法得到的去噪结果,图2d为滤除的噪声剖面。图3为图2黑色框区域的放大图像。由图2和图3可以看出,去除的噪声中基本不包含有效的结构信息,说明本文去噪方法能够保持剖面的边界特征;去噪剖面中基本没有随机噪声的残留,说明本文方法在未伤害有效信号的前提下,能够有效地压制加入的噪声。

图2 模型数据采用本文方法去噪前、后的结果a 合成数据; b 含噪数据; c 去噪结果; d 去除的噪声

图3 图2的局部放大结果a 合成数据; b 含噪数据; c 去噪结果

对图2b所示含噪数据分别采用SVD方法、正则化非平稳自回归(RNA)方法和本文方法进行去噪处理,得到的结果如图4所示。采用SVD方法得到的去噪结果(图4a)噪声残留严重,去噪效果较差,去噪剖面上存在伪影现象;RNA方法的去噪结果较SVD方法要好(图4b),但是去噪剖面上仍然残留噪声;与前两种方法相比,本文方法去噪效果最好(图4c),噪声得到了压制。对比图4中箭头所指区域,本文方法对剖面中不连续性信息保留效果较好,说明本文方法能够在保护有效信号的同时,压制数据中的随机噪声。

2.2 断面保持效果分析

为了验证本文方法的边界保持效果,建立含断层的模型,如图5a所示;对图5a所示模型加入噪声得到的结果如图5b所示;采用本文方法对图5b 进行去噪处理后的结果如图5c所示;图5d为去除的噪声。由图5d可以明显看出,本文方法能够有效压制原始剖面中的随机噪声,而且能较好地保护边缘处的不连续性信息,噪声剖面中无有效信号。

图4 对模型数据采用不同方法得到的去噪结果a SVD方法; b RNA方法; c 本文方法

图5 含断层模型采用本文方法去噪前、后的结果a 合成数据; b 含噪数据; c 去噪结果; d 去除的噪声

3 实际数据处理

对某工区实际地震资料采用本文方法进行去噪处理。图6a为过该工区内的一条叠加地震剖面,该剖面含有较高比例的随机噪声,同相轴模糊、边界及断点不够清晰、地质接触关系模糊,给后续处理解释工作带来极大困难。采用本文方法处理后的地震剖面如图6b所示。对比图6a和图6b可以看出,去噪后的地震剖面中同相轴连续性得到加强,断点清晰、断层展布特征清楚。图6c和图6d分别为典型区域(图6a和图6b中黑框所示区域)放大显示结果。从图6c和图6d可以看出,原始数据中的噪声得到有效衰减,处理后的地震数据中反射波同相轴的连续性得到了增强,且断层特征保持良好,说明本文方法能为层位追踪、断层识别和构造解释等后续处理提供信噪比较高的地震数据。图7a 和图7b分别给出了对实际地震资料采用本文方法进行去噪处理前、后的瞬时振幅属性。对比图7a 和图7b可以看出,去噪前瞬时振幅属性较为杂乱,难以识别断点和岩性突变点;采用本文方法去噪处理后改善了地震资料的品质,瞬时振幅属性变得清晰和容易识别,不连续性信息得到很好的突出,增强了地质体边缘的清晰度(如图中箭头所指区域)。

图6 实际地震资料采用本文方法去噪处理前、后的地震剖面a 去噪前的数据; b 去噪后的数据; c 图6a局部放大结果; d 图6b局部放大结果

图7 对实际地震资料采用本文方法进行去噪处理前(a)、后(b)的瞬时振幅属性

4 结论

地震资料去噪处理需要既能够衰减噪声,提高信噪比,同时又能保留和增强地震反射信息中有效的不连续性信息。本文研究了基于噪声水平估计的加权核范数最小化噪声压制方法,首先采用基于尺度不变性噪声估计方法,对块匹配后的地震数据进行噪声估计,根据噪声方差估计归一化加权核范数最小化噪声压制算法中的保真项,实现对地震数据的自适应去噪处理。主要结论包括:

1) 用本文方法实现了对地震数据的去噪处理,并保护有效信号,尤其在同相轴不连续区域,避免了采用简单阈值处理方法产生的伪吉布斯震荡现象;

2) 数值试验和实际数据的应用结果表明,本文方法能够在不破坏有效信息的同时,自适应地衰减地震数据中的随机噪声,处理后的地震数据信噪比得到提高,有利于后期的构造解释、断层和断点识别、层位追踪、几何属性提取等。