地球物理学报  2015, Vol. 58 Issue (12): 4568-4575   PDF    
GNMF小波谱分离在地震勘探噪声压制中的应用
田雅男1,2, 李月2, 林红波2, 吴宁2    
1. 吉林大学地球科学学院, 长春 130061;
2. 吉林大学通信工程学院, 长春 130012
摘要: 地震勘探资料噪声压制及信噪比提高是整个地震勘探信号处理过程中的重要任务,随着地震勘探深度的增加及其复杂性,人们对地震数据质量的要求越来越高.勘探环境的复杂化使得采集到的地震资料中有效信号被大量噪声淹没,无法清晰辨识,严重影响后续的数据处理与解释.小波去噪是地震勘探中常用且发展较成熟的一种方法,但是其涉及到的阈值函数选取问题一直令人困扰,虽然已有多种阈值函数被提出,但仍存在各自的缺陷.本文利用小波分解在时域及频域良好的信号细节体现特性,引入模式识别中的非负矩阵分解(NMF)谱分离思想,针对小波系数阈值优化问题,提出了一种小波域图非负矩阵分解(GNMF)消噪算法.该方法首先在小波分解基础上,利用GNMF算法实现小波分解系数谱中信号分量与噪声分量的谱分离,然后通过反变换重构各分离子谱对应的子信号,最后利用K均值聚类算法将得到的多个子信号划分为信号类及噪声类,最终得到重构信号及分离噪声.合成记录和实际地震资料的消噪结果验证了新方法在提高信号与噪声分离准确性和精度方面的有效性,同时新方法避免了阈值选取造成的噪声压制不理想或有效成分损失问题.与小波消噪结果的对比及数值分析也说明了新方法在噪声压制及有效成分保持方面的优势.
关键词: 地震勘探     噪声压制     小波分解     谱分离     信噪比    
Application of GNMF wavelet spectral unmixing in seismic noise suppression
TIAN Ya-Nan1,2, LI Yue2, LIN Hong-Bo2, WU Ning2    
1. College of Earth Sciences, Jilin University, Changchun 130061, China;
2. College of Communication and Engineering, Jilin University, Changchun 130012, China
Abstract: Seismic noise suppression and signal-to-noise ratio (SNR) improvement is an important task in the process of seismic signal processing. With the development of seismic exploration in depth and complexity, the requirement of the quality of seismic data is becoming higher. The complexity of acquisition environment makes seismic data mixed with a lot of noise, which makes the effective signal difficult to identify. It directly affects the follow-up data processing and interpretation process. Wavelet denoising is a common and mature method, which is used in seismic exploration. However, the selection of threshold function is always a troubling problem, which hindering its performance. Although many improved methods have been proposed, they still had some shortcomings, respectively.

The proposed method applies the popular graph nonnegative matrix factorization (GNMF) spectral separation theory to seismic random noise suppression. In the constraint function of GNMF, an additional item is added to the conventional constraint function of NMF, which plays an important role in the performance and accuracy improvement of GNMF for the wavelet coefficients spectrums unmixing. It makes full use of the wavelet decomposition characteristics that embodiment the details of the signals in time domain and frequency domain. The novel method first uses GNMF to separate the wavelet coefficient spectrums into some sub-spectrums, and then reconstructs the corresponding sub-signals from these sub-spectrums through inverse transform. Then, it classifies the sub signals into signal class and noise class by K-mean clustering algorithm. The sum of the sub-signals in signal class is the effective signal and the sum of the sub-signals in noise class is the separated noise.

The novel method effectively improves the accuracy and precision of separation signal and noise in the spectral space. Meanwhile, it avoids the worse noise suppression or serious energy loss problems, which is caused by threshold selection. The experimental results on synthetic records and actual seismic data both show the effectiveness and advantages of the proposed algorithm in the aspect of noise suppression and effective composition maintain. Firstly, based on wavelet decomposition, GNMF algorithm is used to decompose the signal components and the noise components. Then, the sub signals are decomposed into signal and noise by using K means clustering algorithm. A new method avoids the problem of noise suppression caused by the threshold selection and the loss of the effective components. Meanwhile, the effectiveness of the method in improving the accuracy and precision of the signal and noise separation is verified. Compared with the results of wavelet de-noising, the proposed method is more effective and has great advantages in noise suppression and effective maintenance.

In order to solve the problem of threshold selection in wavelet denoising method, the wavelet method is applied to the optimization problem of wavelet coefficient threshold in time domain and frequency domain. A novel method based upon GNMF is proposed. The novel method resolves the problem of the threshold selection in the traditional wavelet de-noising method and it can suppress the noise and protects the effective signal components effectively. The simulation and real data results both show that the novel wavelet denoising method has obvious advantages over the traditional wavelet de-noising method.

Key words: Seismic exploration     Noise suppression     Wavelet decomposing     Spectral unmixing     Signal-to-noise ratio (SNR)    
 1 引言

随着地震勘探向着深度和复杂性方向发展,人们对地震数据质量的要求变得越来越高.数据采集环境的复杂化使得地震资料中混入大量噪声,致使有效信号淹没其中难以辨识,直接影响后续的数据处理和解释工作.目前国内外已有很多不同的去噪方法,其中比较经典的方法包括中值滤波(Liu,2009)、维纳滤波(Chen,2006)、F-X滤波(Abma,1995)、Curvelet变换(Neelamani,2008)、小波去噪(Zhang,2001)、经验模态分解(李雪英等,2011王云专等,2010)和隐马尔科夫模型(王金芳等,2009)等等.这些方法虽然可以在一定程度上压制噪声,但仍存在一些不足.

小波变换是一种经典的时频分析方法并已成功应用于地震勘探噪声压制.与短时傅里叶变换不同,它可以根据信号局部特征调节窗长,在时域和频域都具有很好的信号细节体现特性.基于小波的消噪方法通常是利用信号和噪声具有不同的尺度系数,对相应系数施加某种方式的非线性处理(如选取某种阈值函数),从而实现噪声消减.但小波阈值函数的选取一直是令人困扰的问题,阈值选取不当要么造成有效信号成分的严重损失,要么造成大量噪声残留.在地震勘探噪声压制中最常用的阈值估计方法包括三类:小波系数最大值法(Wu,2006),相邻尺度小波变换的空间相关性法(张三宗等,2006Xu,1994)和小波阈值收缩法(Donoho,1994Donoho,1995).然而这些阈值估计方法在叠前地震数据信噪比较低时,很难获得准确的阈值来划分信号和噪声分量,故而存在分量判别不准确问题.

本文为优化小波方法中阈值选取问题,利用小波分析法在时域及频域良好的信号细节体现特性,引入图非负矩阵分解(GNMF)算法,提出了一种地震勘探记录噪声压制的GNMF小波谱分离方法.由于GNMF算法在传统非负矩阵分解(NMF)基础上添加了散度广义准则对解加以约束,因此可以在谱空间内提高噪声和信号分离的准确性和精度.新算法成功解决了阈值选取不当造成的有效能量损失或噪声压制不理想现象.

2 GNMF小波谱分解法 2.1 小波去噪及其缺点

在过去几年里,小波变换发展迅速并成功应用于地震勘探资料消噪处理.它的主要优点是可以根据信号的局部特性调节窗长,在时域及频域都具有良好的局部信号细节体现特性.基于小波的去噪方法通常是利用信号和噪声具有不同的尺度系数,通过对这些小波系数施加某种非线性操作实现消噪过程(Alex and er,2006).图 1给出了小波去噪的过程示意图.在小波分解过程中,以三次小波分解为例,首先含噪信号X分别经过低通和高通滤波器并进行降采样处理得到低频小波系数cA1和高频小波系数cD1,实现一次分解;再对低频小波系数cA1进行相同操作得到低频小波系数cA2和高频小波系数cD2,实现二次分解;继续对低频小波系数cA2进行分解得到低频小波系数cA3和高频小波系数cD3,实现三次分解.低频系数cA1,cA2和cA3中主要包含信号的低频信息,而高频系数cD1,cD2和cD3中主要包含信号的高频信息.由图 1的结构图可知,经过三次分解后含噪信号X可以表示为

图 1 小波去噪过程 Fig. 1 Process of wavelet denoising

由于噪声通常都包含在高频系数中,因此利用某种阈值函数对分解结果进行处理,将小波系数与设定阈值进行数值对比,如果小波系数大于设定阈值,就保留这些系数;反之,将其置零,实现噪声成分消减.接着通过对保留的小波系数进行重构与重组,得到最终的去噪结果.但由于地震勘探信号的特点及其采集环境的复杂性,噪声谱与信号谱在中低频部分是彼此重叠的,利用现有的阈值设定方法去噪后仍存在较多的低频噪声.为了充分去除残留低频噪声,本文引入模式识别中的非负矩阵分解(NMF)思想,解决小波消噪中的阈值优化问题.

2.2 GNMF算法实现

近几年来,非负矩阵分解(NMF)成为模式识别领域中一种流行的矩阵分解方法,它可以有效地获 取低维空间内数据潜在的结构信息(Ren,2014Liu,2012Virtanen,2007Buciu,2004).假设一个M×N维的非负矩阵 X,NMF试图寻找一个M×K维的非负矩阵 U和一个N×K维的非负矩阵 V,使得

利用经典的欧氏距离目标函数 F0=‖ XUV T2,可以得到关于 UV 的更新法则如下:

在此基础上利用局部几何结构实现数据样本点建模,并在目标函数中引入假设条件:如果两个数据点在原始数据空间离得很近,那么在新坐标基的表 达空间依然离得很近(Jeong,2009Cai,2011Yang,2011),进而得出GNMF算法.其目标函数就是在原目标函数基础增加了正则项 R:

其中 z 为样本表达空间的点,Tr(·)为矩阵的迹,W 是加权矩阵,用来控制数据点间的相近程度,D 是对角阵,D jj=∑l W jl. 通过最小化式(5)中的 R 就可以实现前面的假设条件,并增强数据表达空间内数据点的归类精度.于是,GNMF的目标函数可表示为 F0 与正则项 R 之和:
其中 λ1 是正则因子.通过推导可以得到GNMF中关于 UV 的更新法则如下:
与式(3)和(4)的更新法则相比,约束条件的增加提高了数据特征提取的精度和准确性,为后续GNMF算法应用于小波谱分离提供重要的前提条件.

2.3 GNMF实现小波谱分离

首先,我们利用小波去噪方法对人工合成的单道含噪地震信号进行去噪处理.图 2a是一道由雷克子波仿真生成的纯净地震信号,加入高斯白噪声后得到信噪比为-3 dB的含噪信号如图 2b.图 2c到2f分别给出了基于db4小波对该含噪信号进行三阶分解并系数重构后得到的子信号.其中图 2c为低频系数cA1重构出的子信号,图 2d,2e,2f分别为高频系数cD1,cD2和cD3重构出的子信号.图 2g到 2j为上述子信号对应的时频谱.由这些子信号的时域波形和时频谱可以看到高频子信号中包含的主要是随机噪声,有效的信号分量主要集中在图 2c中的低频子信号里.

图 2 db4小波分解 (a)和(b)分别为纯净信号和含噪信号;(c)—(f)为小波系数cA1,cD1,cD2和cD3对应的重构信号;(g)—(j)为(c)—(f)对应的时频谱. Fig. 2 ‘db4’ wavelet decomposition (a)Pure signal;(b)Noisy signal;(c)—(f)Sub-signals of cA1,cD1,cD2 and cD3;(g)—(j)Spectrums of(c)—(f).

但是图 2c中的低频子信号中仍然存在一些低频噪声未被分离.由图 2g同样可以看出低频系数谱中仍有部分残留噪声谱与信号谱混叠在一起.对于这些噪声,简单的阈值函数法无法有效将其消除.于是为优化小波去噪方法中阈值选取问题,进一步去除残留噪声,本文引入GNMF算法实现小波分解系数谱中信号分量与噪声分量的有效分离.

首先在上述小波分解前提下,利用GNMF算法将图 2g到2j中的每个系数谱分别分解为两个子谱,再通过短时傅里叶反变换重构出与每个子谱相对应的子信号,如图 3所示.

图 3 谱分解后重构子信号 图 2的小波分解得到的两个子信号;(a、e)对应图 2c;(b、f)对应图 2d;(c、g)对应图 2e;(d、h)对应图 2f. Fig. 3 Reconstructed sub-signals after wavelet spectrum unmixing (a) and (e)are the two subsignals of Fig. 2c;(b) and (f)are the two subsignals of Fig. 2d; (c) and (g)are the two subsignals of Fig. 2e;(d) and (h)are the two subsignals of Fig. 2f.

图 3a和3e是图 2g中系数谱的两个子谱,图 3b和3f是图 2h中系数谱的两个子谱,图 3c和3g是图 2i中系数谱的两个子谱,图 3d和3h是图 2j中系数谱的两个子谱.从图中的子信号波形可以看出,图 3b到3h均为重构出的噪声分量,将这些分量求和便得到分离出的噪声,图 3a为重构出的信号成分.将其与小波去噪结果图 2c进行波形对比可以看出,利用GNMF算法对小波分解系数谱进行分解后,低频噪声得到了进一步的有效分离,同时有效信号幅度没有出现严重衰减现象,依然保持完整.可见,新方法在小波阈值优化方面具有明显效果.

3 实验结果 3.1 合成地震记录

为了测试新方法的消噪性能,我们将其应用于一张80道的合成地震记录,每道1024点,采样频率为1000 Hz,主频为30 Hz,道间距为50 m,层速度 分别为1800 m·s-1和2500 m·s-1,信噪比为-3 dB,如图 4a所示.在含噪记录中存在大量随机噪声,局部同相轴变得模糊.图 4b和4c分别为采用小波去噪和本文提出的GNMF小波系数谱分离法对该含噪记录进行噪声压制后的多道效果图.

图 4 合成记录消噪结果 (a)含噪信号;(b)小波去噪结果;(c)GNMF小波谱分离去噪结果;(d)两种方法去噪后残留噪声MSE对比;(e)两种方法去噪后时域波形对比. Fig. 4 Results of synthetic records (a)Noisy signal;(b)Result of wavelet denoising;(c)Result of GNMF wavelet spectal unmixing;(d)MSEs of the denoising results after the two methods;(e)Waveforms comparison of the denoising results after the two methods.

通过对比可以看出,两种方法去噪后大部分随机噪声都得到了有效压制,同相轴变得清晰连贯;但是相比之下,图 4c中GNMF小波系数谱分解法处理后残留的背景噪声明显少于图 4b中小波去噪处理结果;同时有效成分保持完整,没有发生衰减和畸变.为了对比两种方法去噪后信号中残留噪声情况,图 4d绘制了含噪记录中噪声及两种方法去噪后残留噪声的多道均方差平均值曲线.从图中可以看出,含噪记录中噪声均方差水平较高,GNMF小波谱分解法去噪后的残留噪声均方差较小波变换去噪后残留噪声均方差降低了约0.03,可见新方法对噪声的压制效果更好.为了更清晰地对比两种方法的效果 差异,我们分别从图 4b图 4c中选取一道数据进行波形对比,如图 4e所示.从图中两个波峰放大图可以看出,本文提出的新方法在有效信号幅度保持 方面与传统小波去噪相比更完整,更重要的是两种方法在噪声压制方面新方法较小波去噪具有明显优势.

为了对两种方法效果进行量化对比,表 1给出了两种方法处理后的信噪比、幅度衰减及均方误差.

表 1 信噪比、幅度衰减及RMSE对比 Table 1 Comparisons of SNRs,amplitude loss and RMSE

从这些数据可以判断出小波去噪后信噪比提高了9.8019 dB,GNMF小波系数谱分解法去噪后信 噪比提高了15.0763 dB,比小波去噪提高了5.2744 dB;同理,新方法去噪后信号幅度衰减比小波去噪降低 了7.4114%,均方差比小波去噪降低了0.0331,可 见GNMF小波系数谱分解法充分压制了随机噪声的同时,完整地保持了有效信号.

3.2 实际地震资料

最后,我们将GNMF方法应用于实际地震记录噪声压制.图 5a是某地区的一张168道,每道6145点的实际共炮点地震记录,其采样频率为1000 Hz.样本中不仅存在随机噪声,还有直流干扰噪声,且左 下角A区域随机噪声较强,信噪比较低.分别采用 GNMF小波谱分离及小波去噪方法对其进行噪声压制,如图 5b图 5c.

图 5 实际地震记录消噪结果 (a)含噪记录;(b)小波去噪结果;(c)GNMF小波谱分离去噪结果. Fig. 5 Results of real seismic records (a)Noisy signal;(b)Result of wavelet denoising;(c)Result of GNMF wavelet spectal unmixing.

对比实验结果可以看出,小波去噪后区域A颜色明显变浅,大部分随机噪声得到有效压制,区域B内同相轴也变得清晰;但是两区域内仍有噪声残留.观察图 5c可以看出,GNMF小波谱分离去噪后,A和B两区域内颜色进一步变浅,说明残留的随机噪声得到进一步压制,同时同相轴清晰、连贯.

4 结论

本文为解决小波去噪方法中阈值选取问题,利用小波方法在时域和频域良好的信号细节体现特性,将模式识别中的GNMF思想应用到小波系数阈值优化问题,提出了一种GNMF小波谱分离消噪方法,并应用于地震勘探资料噪声压制.新方法解决了传统小波消噪方法中的阈值选取造成的噪声压制不理想及有效信号成分衰减或损失问题,有效压制噪声的同时能够完好保持地震勘探信号.仿真实验结果及实际资料处理结果都表明新方法与小波去噪相比可实现噪声的进一步充分消减,同时对有效成分保持完整,具有明显优势.

参考文献
[1] Abma R, J Claerbout. 1995. Lateral prediction for noise attenuation by t-x and f-x techniques. Geophysics, 60:1887-1896.
[2] Buciu I, Pitas I. 2004. Application of nonnegative and local non negative matrix factorization to facial expression recognition. Proceedings of the 17th IEEE International Conference on Pattern Recognition. Cambridge, UK, 288-291.
[3] Cai D, He X, Han J. 2011. Graph regularized nonnegative matrix factorization for data representation. IEEE Trans. Pattern Anal. Mach. Intell., 331548-1560.
[4] Chen J D, Benesty Huang. 2006. New insights into the noise reduction wiener filter. IEEE Ransactions on Udio, Peech and Anguage Rocessing, 14(4):1218-123.
[5] Donoho D L, Johnstone. 1994. Ideal time-frequency denoising Technical Report, Dept. of Statistics, Stanford University.
[6] Donoho D L. 1995. De-noising by soft-thresholding. IEEE Trans. Information Theory, 41(3): 613-627.
[7] Alexander D. 2006. Theory and seismic applications of the eigen-image discrete wavelet transform, European association of geoscientists & engineers. Geophysical Prospecting, 54:441-461.
[8] Jeong S Y, Kim K H, Jeong J H. 2009. Semi-blind disjoint non-negative matrix factorization for extracting source from single channel noisy mixture. IEEE Workshop on Applications of Signal Processing to Audio and Acoustics, 73-76.
[9] Li X Y, Sun D, Hou X H, et al. 2011. Comparison of empirical mode decomposition and generalized S transform for prestack noise suppression. Progress in Geophysics, 26(6):2039-2045.
[10] Liu H, Wu Z, Li X. 2012. Constrained nonnegative matrix factorization for image representation. IEEE Trans. Pattern Anal. Mach. Intell., 34, 1299-1311.
[11] Liu Y, Liu C, Wang D. 2009. A 1D time-varying median filter for seismic random, spike-like noise elimination. Geophysics, 74(1), V17.
[12] Neelamani R, Baumstein A I, Gillard D G. 2008. Coherent and Random noise attenuation using the curvelet transform. The Leading Edge, 27:240-248.
[13] Ren W Y, Li G H, Tu D. 2014. nonnegative matrix factorization with regularizations. IEEE Journal on Emerging and Selected Topics in Ircuits and Ystems, 4(1):153-164.
[14] Virtanen T. 2007. Monaural sound source separation by nonnegative matrix factorization with temporal continuity and sparseness criteria. IEEE Ransactions on Audio, Speech and Language Processing,15:1066-1074.
[15] Wang J F, Li Y, Wang J B, et al. 2009. Random noise suppression based on Hidden Markov model smoothing estimation. Progress in Geophysics, 24(5):1861-1867.
[16] Wang Y Z H, Lan J T, Long Y S H. 2010. Random noise suppression based on S transform. Progress in Geophysics, 25(2):562-567.
[17] Wu Z C, Liu T Y, Hua C H. 2006. Wavelet threshold de-noising based on higher-order statistics in attenuating random noise. 76th SEG Technical Program Expanded Abstracts, 2857-2861.
[18] Xu Y, Weaver J B, Healy D M, et al. 1994. Wavelet transform domain filters: A spatially selective noise filtration technique. IEEE Trans. Image Processing, 747-758.
[19] Yang Z Y, Zhou G X, Xie S L. 2011. Blind spectral unmixing based on sparse nonnegative matrix factorization. IEEE Transactions on Image Processing, 20:1112-1125.
[20] Zhang L, Bao P, Pan Q. 2001. Threshold analysis in wavelet-based denoising. IEEE Electronics Letters, 37(24): 1485-1486.
[21] Zhang S Z, Xu Y X. 2006. Higher-order correlative stacking for seismic data in the wavelet domain. Chinese J. Geophys. (in Chinese), 49(2):554-560.
[22] 李雪英, 孙丹, 侯相辉等. 2011. 基于广义S变换、经验模态分解叠前去噪方法的比较. 地球物理学进展, 26(6):2039-2045, doi:10.3969/j.issn.1004-2903.2011.06.019.
[23] 王金芳, 李月, 王金宝等. 2009. 基于隐马尔可夫模型平滑估计的随机噪声压制方法. 地球物理学进展, 24(5): 1861-1867, doi:10.3969/j.issn.1004-2903.2009.05.042.
[24] 王云专, 兰金涛, 龙玉沙. 2010. 基于S变换的随机噪声压制方法. 地球物理学进展, 252:562-567, doi:10.3969/j.issn.1004-2903.2010.02.026.
[25] 张三宗, 徐义贤. 2006. 地震记录小波域高阶相关叠加技术. 地球物理学报, 492:554-560.