因各种干扰的存在, 地震资料中会存在随机噪声, 有些衰减噪声的滤波处理方法会造成地震反射中的不连续性信息被平滑而变得模糊^[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(α)直接取决于其形状参数α:

$ \kappa_{x}(\alpha)=\frac{\mathit{\Gamma}\left(\frac{1}{\alpha}\right) \mathit{\Gamma}\left(\frac{5}{\alpha}\right)}{\mathit{\Gamma}\left(\frac{3}{\alpha}\right)^{2}} $

(1)

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

$ \kappa=\frac{\mu_{4}}{\sigma^{4}} $

(2)

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

$ \sigma_{y}^{2}=\sigma_{x}^{2}\left(1+\frac{\sigma_{n}^{2}}{\sigma_{x}^{2}}\right) $

(3)

y的四阶中心矩满足:

$ \mu_{4}(y)=3 \sigma_{x}^{4}\left(1+\frac{\sigma_{n}^{2}}{\sigma_{x}^{2}}\right)^{2}+\sigma_{x}^{4}\left(\kappa_{x} \alpha-3\right) $

(4)

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

$ \kappa_{y}=\frac{\kappa_{x} \alpha-3}{\left(1+\frac{\sigma_{n}^{2}}{\sigma_{x}^{2}}\right)^{2}}+3 $

(5)

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

$ \boldsymbol{b}(m, n)=\cos \left[\frac{(2 m+1) p \pi}{2 N}\right] \cos \left[\frac{(2 n+1) q \pi}{2 N}\right] $

(6)

式中:b(m, n)为生成的基函数矩阵; m, n为数据在空间域中的坐标; p, q为基函数在频率域中的坐标; N为数据分割的块大小, 根据经验一般取值为8。m, n, p, q的取值范围为0~N－1。用N×N个以b(m, n)为基的滤波器对数据进行卷积, 生成响应数据, y_i表示第i个响应数据, 估计该响应数据的方差和峰度获得$\hat \sigma _{{y_i}}^2 $和$ {{\hat \kappa }_{{y_i}}}$, 对在范围(2, 3, …, N²)内的每个分量i执行此过程, 通过计算(7)式的极值来估计噪声方差和峰度值:

$ \mathop {{\mathop{\rm argmin}\nolimits} }\limits_{{\kappa_x}, \sigma _n^2} \sum\limits_{i = 2}^{{N^2}} {\left| {\frac{{{\kappa _x} - 3}}{{{{\left( {1 + \frac{{\sigma _n^2}}{{\hat \sigma _{yi}^2 - \sigma _n^2}}} \right)}^2}}} + 3 - {{\hat \kappa }_{{y_i}}}} \right|} $

(7)

具体步骤如下:

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

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

3) 估计每个响应数据y_i的方差$\hat \sigma _{{y_i}}^2 $和峰度$ {{\hat \kappa }_{{y_i}}}$;

4) 将求取的$\hat \sigma _{{y_i}}^2 $和$ {{\hat \kappa }_{{y_i}}}$代入公式(7), 计算极值来估计噪声峰度和方差。

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

1.2 多道加权核范数最小化

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

$ d\left(y_{i}, y_{j}\right)=\frac{\left\|y_{i}-y_{j}\right\|_{2}^{2}}{N^{3}} $

(8)

式中:“‖·‖₂²”表示两个数据块之间的欧式距离; N为分块的大小, 根据实际数据进行选定, 为充分保留数据的细节特征, 其取值不宜过大。如果两个数据块的相似性小于或等于预设阈值, 则认为它们是相似的。

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

将WNNM算法应用于数据降噪处理, 由Y_j来估计X_j。使用基于尺度不变性的噪声水平估计方法计算的噪声方差$\sigma _n^2 $来归一化F范数的保真项$\parallel {\mathit{\boldsymbol{Y}}_j} - {\mathit{\boldsymbol{X}}_j}\parallel _F^2 $ ^[24], 可以得到:

$ {{\mathit{\boldsymbol{\hat X}}}_j} = \mathop {{\mathop{\rm argmin}\nolimits} }\limits_{{\mathit{\boldsymbol{X}}_j}} \frac{1}{{\sigma _n^2}}\left\| {{\mathit{\boldsymbol{Y}}_j} - {\mathit{\boldsymbol{X}}_j}} \right\|_F^2 + {\left\| {{\mathit{\boldsymbol{X}}_j}} \right\|_{\mathit{\boldsymbol{w}}, *}} $

(9)

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

$ w_{i}=c \sqrt{m} /\left[\lambda_{i}\left(\boldsymbol{X}_{j}\right)+\varepsilon\right] $

(10)

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

$ \hat{\lambda}_{i}\left(\boldsymbol{X}_{j}\right)=\sqrt{\max \left[\lambda_{i}^{2}\left(\boldsymbol{Y}_{j}\right)-m \sigma_{n}^{2}, 0\right]} $

(11)

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

将上述步骤应用于每个分块并将所有分块组合, 可以获得去噪后的数据$ {\hat x}$, 具体步骤如下:

1) 对于含噪数据, 初始化$ \hat{x}^{0}=y, y^{0}=y$, 给定迭代次数L, 迭代权重参数δ;

2) 迭代执行$y(k)=\hat{x}^{k-1}+\delta\left(y-\hat{y}^{k-1}\right) $;

3) 对y(k)中的每一小块y_j寻找相似分块集合Y_j, 对Y_j进行奇异值分解[U, Σ, V]=SVD(Y_j), 根据(10)式计算权重向量w, 得到$ \hat{\boldsymbol{X}}_{j}=\boldsymbol{U S}_{w} \boldsymbol{\mathit{\Sigma} V}^{\mathrm{T}}$, 其中, $\boldsymbol{S}_{w} \boldsymbol{\mathit{\Sigma}}_{i i}=\max \left(\boldsymbol{\mathit{\Sigma}}_{i i}-w_{i}, 0\right) $;

4) 聚集${{\mathit{\boldsymbol{\hat X}}}_j} $来重构去噪数据${{\hat x}^k} $, 重复步骤2), 直至循环结束。经过上述步骤得到的${\hat x} $便是去噪后的数据。

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

2 理论模型试算 2.1 模型数据去噪效果分析

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

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

2.2 断面保持效果分析

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

3 实际数据处理

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

4 结论

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

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

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