Automatic noise attenuation based on clustering and empirical wavelet transform

Accepted Manuscript Automatic noise attenuation based on clustering and empirical wavelet transform

Wei Chen, Hui Song PII: DOI: Reference:

S0926-9851(18)30501-9 doi:10.1016/j.jappgeo.2018.09.025 APPGEO 3612

To appear in:

Journal of Applied Geophysics

Received date: Revised date: Accepted date:

3 June 2018 6 August 2018 19 September 2018

Please cite this article as: Wei Chen, Hui Song , Automatic noise attenuation based on clustering and empirical wavelet transform. Appgeo (2018), doi:10.1016/ j.jappgeo.2018.09.025

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT Automatic noise attenuation based on clustering and empirical wavelet transform

Wei Chen* [email protected] and Hui Song [email protected] Key Laboratory of Exploration Technology for Oil and Gas Resources of Ministry of Education

CE

PT

ED

M

AN

US

CR

Corresponding author at: Daxue Road No.111, Caidian District, Wuhan, China, 430100

AC

*

IP

Cooperative Innovation Center of Unconventional Oil and Gas

T

Yangtze University Daxue Road No.111 Caidian District, Wuhan, China, 430100 Hubei

ACCEPTED MANUSCRIPT ABSTRACT Strong noise in seismic data seriously affects many steps in seismic data processing and imaging. While most traditional methods depend on carefully tuned input parameters by humana lot of input parameters that are carefully tuned by human, we are proposingthe method proposed in this paper

T

is an automatic noise attenuation algorithm to faciliteaiming at facilitating a fast preprocessing of

IP

massive prestack seismic data. In the proposed algorithm, the non-stationary seismic data is first

CR

adaptively decomposed into several empirical components adaptively via empirical wavelet transform (EWT) according the frequency contents in the data. Then, the first component is

US

selected to stand for the useful signals. To deal with the residual noise in the roughly estimated

AN

signal infrom the previous step, we propose a clustering based thresholding method. The most dominant signals are detected via a simple clustering steps and other components are damped with

M

an adaptive percentile threshold. The two steps refer to a new automatic algorithm to denoise the

ED

seismic data with high fidelity.We demonstrate the performance of the proposed method via both

AC

CE

PT

synthetic and field data examples.

ACCEPTED MANUSCRIPT INTRODUCTION The attenuationAttenuation of random seismic noise is a long-standing problem in the reflection seismology community (Claerbout, 1976; Yilmaz, 1987; Deighan and Watts, 1997; Welland, 2003; Mallat, 2009; Tian et al., 2013; Zhang et al., 2016b; Bai et al., 2016b; Nazari Siahsar et al.,

T

2016; Shahin et al., 2017; Guo et al., 2017; Ebrahimi et al., 2017; Guo and Mcmechan, 2017; Ryu

IP

et al., 2017; Zhang et al., 2017b; Huang et al., 2018; Wu and Bai, 2018a,b) and passive seismology

CR

community (Han and van der Baan, 2015; Li et al., 2016a,b; Chen, 2018a,c). Because of the crucial benefits of removing random noise to many seismic processing and imaging tasks (Chen et al.,

US

2017; Zhang et al., 2018c,b), a lot of methods to attenuate seismic noise, also known as denoise,

AN

have been investigated in the literature in the past decades (Zu et al., 2017a; Asgedom et al., 2017; Zeng et al., 2017; Jun et al., 2017; Wang et al., 2017a; Zhang et al., 2017a; Zhou and Han, 2018;

M

Zu et al., 2017b, 2018; Zhang et al., 2018c). The classic denoising algorithms include the

ED

prediction based methods (Canales, 1984a; Liu et al., 2012; Liu and Chen, 2013; Chen and Ma, 2014), sparse transform based methods (Liu and Fomel, 2010; Fomel and Liu, 2010; Chen and

PT

Fomel, 2018; Bai and Wu, 2018), and rank-reduction based methods (Chen et al., 2016c; Huang et

CE

al., 2017a). The prediction based methods assume that the seismic data are spatially predictable. The sparse transform based methods are based on the sparsity of the seismic data in a transform

AC

domain, e.g., Fourier transform (Zhong et al., 2016; Zhou, 2017), curvelet domain (Candés et al., 2006; Candés and Plan, 2010; Cao et al., 2014; Liu et al., 2016d; Zu et al., 2017c; Wang and Wang, 2017), wavelet domain (Deighan and Watts, 1997; Welland, 2003; Liu et al., 2016f; Yang et al., 2018), seislet domain (Fomel and Liu, 2010; Liu and Fomel, 2010; Gan et al., 2015b, 2016c; Liu et al., 2016e; Gan et al., 2016a; Wu et al., 2016), Radon domain (Xue et al., 2016b, 2017; Chen, 2018b), shearlet domain (Kong and Peng, 2015; Liu et al., 2016a), polynomial transform domain

ACCEPTED MANUSCRIPT (Chen et al., 2018), and learned dictionary domain (Chen et al., 2016b; Siahsar et al., 2017; Zhou et al., 2017a; Chen, 2017). The rank-reduction based methods assume that the rank of prearranged matrix that is constructed from reflection seismic data is small and thus a denoising step is equal to a rank-reduction step, e.g., using the singular value decomposition (SVD) (Gan et al., 2015a;

T

Zhang et al., 2016a; Xue et al., 2016a; Wang et al., 2017b; Zhou et al., 2017b; Zhou and Zhang,

IP

2017; Zhou et al., 2018). Besides, statistics based methods such as the median/mean filters are also

CR

widely used in the literature (Yang et al., 2015; Gan et al., 2016d; Zu et al., 2016b; Chen et al., 2016a; Bai and Wu, 2017; Xie et al., 2018; Bai et al., 2018; Zhang et al., 2018a; Chen et al., 2019;

US

Zhou and Li, 2018).

AN

Recently, some novel methods have been proposed to suppress the random noise, such as the waveform shaping method (Chen and Jin, 2016), which separates the noise and useful signal by

M

shaping the estimated wavelet and the inverted model to a more admissible model. However, the

ED

waveform shaping method cannot preserve the weak signal. Chen and Fomel (2015) developed a two-step denoising method to avoid signal leakage. Generally speaking, most of the

PT

aforementioned methods need to subdividesub-divide the data into smaller panels where the

CE

events are linear, and cannot analyze well the local characteristics of non- linear and non-stationary signals. Huang et al. (1998) proposed the empirical mode decomposition (EMD) to solve the

AC

non- linear and non-stationary problems in geophysics and other research field (Chen, 2016; Gan et al., 2016b; Liu et al., 2016b, 2017, 2018a,b). EMD can adaptively decompose a signal containing many frequency components into corresponding frequency signals. Thus, the EMD can separate the single signal and noise effectively even for complicated seismic data. The dominant frequency of the intrinsic mode function (IMF) decomposed by the EMD monotonically decreases, and we can remove the high- frequency components to obtain the denoised data when considering the

ACCEPTED MANUSCRIPT dominant frequency of random noise is higher than the useful signals. The EMD and f-x EMD methods (Chen and Ma, 2014) are used to remove noise adaptively, and the residual data contains almost no useful signal. Although the EMD can solve most seismic data denoising problems forwith non- linear and non-stationary signals, Wu and Huang (2009) found that EMD cannot

T

overcome the mode mixing problem due to signal interruption. When mode mixing problem

IP

occurs, different frequency components are mixed in one or more IMFs, which causes difficulties

CR

in interpreting the final time frequency distribution. For this reason, Wu and Huang (2009) superposed signals with white noise to effectively reduce mode mixing, which is the so-called

US

Ensemble Empirical Mode Decomposition (EEMD). Nevertheless, the conventional EEMD

AN

denoising method cannot effectively reduce the noise in the high frequencyhigh- frequency IMFs. Especially for noise suppression in multichannel seismic data, selecting the same high- frequency

M

components for different traces will more or less cause the loss of useful signals. Moreover, the

ED

EMD-based and EEMD-based methods lack mathematical basis and thus are not easy to control. On the other hand, some tests show that the EMD-based or EEMD-based method behaves

PT

just like an adaptive filter bank, resulting in the existence of some indescribable modes. In order to

CE

overcome the aforementioned negative features of the EMD-based and EEMD-based methods, a novel method called empirical wavelet transform (EWT) is recently proposed, which has better

AC

denoising performance than EEMD and has mature mathematical theory basis based on wavelet and amplitude modulated-frequency modulated (AM-FM) framework (Gilles, 2013; Liu et al., 2016c). EWT-based method can also adaptively decompose a signal containing many frequency components into corresponding different frequency signals as EMD or EEMD does, but the modes modeling from raw signal can be reasonably interpreted. It is worth mentioning that the EWT has been applied in time-frequency analysis in many fields, including the seismic data analysis field

ACCEPTED MANUSCRIPT (Liu et al., 2016c). However, it has never been applied to denoise seismic data so far. While denoising viausing the EWT decomposition enjoys the benefit of being completely automatic, the denoised data still contains some residual noise. To conquer the issue of remaining noise, we develop a novel clustering guided EWT thresholding method. The waveforms are detected via a

T

parameter- free unsupervised clustering algorithm. Other “non-waveform” components are

IP

damped with a percentile thresholding operation. In this way, the proposed method can

CR

automatically attenuate the random noise in seismic data without damaging the useful signals. We present several examples demonstrating the detailed implementation steps of the proposed

US

algorithm, including a single-channel seismic trace, multi-channel synthetic data, and complicated

AN

real seismic data. All examples demonstrate the successful performance of the proposed method.

M

THEORY

ED

Empirical wavelet transform

PT

We denote the analyzed signal as s  t  . In the time domain, wavelet sets  i , j  are defined by mother wavelet function  u ,v with mean value zero, scale factor v  v > 0  and shift factor

CE

u  u > 0  . We can get the wavelet transform Ws  u, v  of signal s  t  using inner product

AC

operation between signal s  t  and wavelet function u  u > 0  . We consider a normalized Fourier axis with period 2 , and make discussions in consideration of Fourier frequency domain as   0,   in order to satisfy the Shannon theorem. Next we disperse the Fourier frequency domain  0,   into N continuous subintervals, so we get N  1 discontinuity points and denote the nth one as n 1 , the first one as 0 = 0 and the last one as N =  . Then we can get the nth subinterval as dn = n1 , n  . It is clear that

ACCEPTED MANUSCRIPT nN=1 dn = 0,   . If we center on discontinuity point n , a transitional zone denoted as Tn whose width is 2 n can be defined. Empirical wavelet is defined as a band-pass filter on the domain d n . In addition, for the

T

half of the width of the transition zone  n , we generally have  n = n , and 0 <  < 1 . Similar to

IP

the structure strategy of Littlewood-Paley wavelet and Meyer wavelet, for n > 0 , we can define

CR

the empirical scale function as: 1,  | 1   )n  cos    1   |  1   )n ))], ˆn  )   2 2n  1   )n   | 1   )n  0, else

(1)

AN

US



and empirical wavelet as:

M

1, 1   )n  w  1    n 1  cos  2  21n 1    1    n 1   ,   1    n 1   | 1   )n 1 ˆ n     sin  2  21 n    1    n   ,   1    n    1    n 0, else 







(2)

CE

PT

ED



AC

In the above two equations,   x  is 0 or 1, and generally we have

  x  = x 4  35  84 x  70 x 2  20 x3  .

(3)

How to divide the Fourier spectrum is of vital importance in EWT, which is directly related to the adaptivity of the decomposition of the original signal. Here we imitate the classical wavelet transform theory to define the EWT, so the detailed coefficients are defined by the inner product of empirical wavelet and original signal:

ACCEPTED MANUSCRIPT Ws  n, t  =  s, n  = s   n   t d .

(4)

Noise attenuation via EWT

T

Ws  n, t  in equation 4 can also be expressed as cn (t ) and the summation of cn (t ) reconstructs

IP

the original signal s(t ) . In an EMD-like format, we have

i =1

CR

N

s(t ) = cn (t ),

(5)

US

where cn (t ) denotes the n th component and is also called intrinsic mode function (IMF). N denotes the number of components. Different components after EWT decomposition usually

AN

correspond to different frequency bands. The first component, i.e., c1 (t ) is usually corresponding

M

to the most distinct/dominant frequency band, which indicates the frequency band of the reflection

ED

signals. The beauty of the EWT method is that the decomposition process is fully automatic. Although EMD has a similar advantage of being adaptive, it is remained as empirical. Because of

PT

the lack of mathematical support, the modes in EMD are not flexibly controlled.

CE

Considering that the first component after EWT decomposition mainly represents the useful signals, a straightforward application of the EWT in denoising is to only select the first

AC

component as the denoised data. As a matter of fact, this method works quite well and due to its single-channel nature, it preserves the edges in seismic data very effectively. Here, we will first use a single-channel example to illustrate the denoising performance using a simple EWT decomposition step. Figure 1 shows a single-trace denoising example using EWT. The example is used to show how the adaptive EWT decomposition can help separate the noise from the signal. Figure 1(a) shows the clean seismic trace. Figure 1(b) shows the noisy data, with a SNR equal to 0.29 dB. Figure 1(c) shows the denoised data, with a improved SNR equal to

ACCEPTED MANUSCRIPT 4.56 dB. 1(d) shows the removed noise using the EWT decomposition. The EWT decomposition methods simply choose the first component after EWT decomposition as the denoised data. Note that the method is completely automatic. To see the difference in the spectra, we plot a comparison of different data in the frequency domain in Figure 2. To see different IMFs after EWT

T

decomposition clearly, we show the five IMFs in Figure 3. Note that the number of IMFs using the

IP

EWT method is adaptively defined by the algorithm itself. It is clear that the second to fift h

CR

components almost contain no distinct signal energy. It is worth mentioning that such EWT decomposition based denoising method has not been reported publicly so farused before.

US

However, as we can see from Figure 1(c), although improved a lot using the simple

AN

adaptive EWT decomposition, there is still a lot of residual noise. A thresholding strategy is preferable to obtain a cleaner data. In the next subsection, we will introduce a clustering guided

ED

M

thresholding algorithm.

Clustering guided EWT thresholding

PT

It can be observed fromin Figures 1 and 3 that after EWT decomposition, the shape of the Ricker

CE

waveform is reconstructed very well reconstructed. The residual noise is mainly spreading over the non-waveform areas. Thus, we propose to first detect those waveform positions with high

AC

probability and damp the remaining areas using a simple thresholding strategy. In order to maintain the automatic benefit of EWT, both waveform detection and thresholding steps should also be automatic. First, in aan easier way, the thresholding step can use a percentile strategy that is previously used in Chen et al. (2014). The percentile thresholding methods adaptively define the threshold by preserving a certain percentage of the maximal values and damping the preserved values with the adaptively calculated cut-off threshold. In our method, we fixed the percentage as

ACCEPTED MANUSCRIPT 50%. A reason for not choosing a small percentage is that too small percentage has the danger of clearing off all the low-amplitude but useful signals. The next more difficult task is to automatically detect the waveforms. A possible solution may be using the clustering method, which is a parameter-free unsupervised machine learning technique.

T

An important component of a clustering algorithm is the distance measure between data

dimensional data p

,

(6)

US

d (xi , x j ) = xi  x j

CR

IP

points. The following Minkowski Metric is a common way for measuring distance for an N

where d (xi , x j ) denotes the distance between x i and x j . x i denotes the i th N dimensional 

denotes p -norm of the input vector. When p = 2

AN

data point. x i and x j are both vectors.

p

M

, it refers to the Euclidean distance, when p = 1 , it refers to the Manhattan metric.

ED

Here, we use a similar clustering method that was proposed by Chen (2018a) However,

1. Mean M

PT

different from Chen (2018a), we only define two features as

CE

M (i) =

1 N

iw

d (i),

(7)

iw

AC

2. Power E

iw

E (i ) = d 2 (i ),

(8)

iw

w denotes half of window length. The window length is fixed to be the length of the wave length. d (i) denotes the input seismic. The Fuzzy c- means (FCM) is used in our method to cluster the

waveforms. The FCM clustering method allows one piece of data to belong to two or more clusters (dunn, 1973; Bezdek, 1981). It is based on the minimization of the following objective function

ACCEPTED MANUSCRIPT N

C

J = uim, j xi  c j

2

,1  m  ,

(9)

i =1 j =1

where m is any real number greater than 1, ui , j is the degree of membership of x i in the cluster

j , x i is the i th of d -dimensional measured data ( d is the dimension of x i ), c j is the d

CR

1. Decompose the input seismic trace into multiple IMFs.

IP

We summarize our proposed algorithm as the following steps

T

-dimension center of the cluster.

US

2. Select the first IMF as the initial guess of the useful signals.

3. Apply FCM clustering algorithm to roughly estimate the signal waveform locations.

AN

4. Damp the remaining “non-waveform” components using a percentile thresholding method. For the synthetic example shown in Figure 1, we show the clustered waveform locations in

M

Figure 4. The red line indicates the label for waveform clustering, value 1 means waveforms and

ED

value 0 means non-waveforms. Figure 5 shows a comparison between denoising methods. It is

PT

clear that the proposed method obtains a very excitingencouraging denoising performance. The SNRs of the data in Figure 5(c) and Figure 5(d) are 4.56 dB and 7.88 dB, respectively.

CE

In the proposed method, the clustering process is an automatic process and the EWT

AC

method is also an automatic method, meaning that the proposed method is a fully automatic method. With other well-tuned methods, we may get better performance, but will lose the adaptivity enjoyed from the fully automatic method.

EXAMPLES In this section, we will test the proposed algorithm on one synthetic example and one field data example. We will compare the performance from the proposed method with thatthose from the

ACCEPTED MANUSCRIPT classic f  x deconvolution method, which is denoted by FX, the singular spectrum analysis (SSA) methodthe EMD method, and the traditional curvelet thresholding. The appendix introduces the basics of the three comparing methods. To quantitatively compare the performance for the synthetic example, we use the signal- to-noise ratio (SNR) (Liu et al., 2009; Qu et al., 2015,

T

2016; Bai et al., 2016a; Zu et al., 2016a; Gan et al., 2016e; Huang et al., 2017b; Wang et al., 2018)

dclean 22 d denoised  dclean

CR

SNR = 10log10

IP

as the metric, which can be defined as follows:

2 2

.

(10)

US

where dclean denotes the exact solution, i.e., the clean data, and d denoised denotes the denoised

AN

data. Both dclean and d denoised denote vectors constructed from the 2D matrix.

M

The synthetic data example is shown in Figure 6, where (a) 6a shows the clean data and (b) shows the noisy data. Figures 6c and 6d show the spectra of the clean and noisy synthetic da ta. The

ED

synthetic data is very realistic, because it contains several typical morphological structures that are

PT

common in seismic data, e.g., the crossing and curving events, faults, and discontinuities. Figure 7 shows a comparison of denoised data using four different methods. Figure 7a shows the result

CE

from the FX method. For the FX method, we use a prediction filter of 6 points length. Figure 7b

AC

shows the result from the EMD method, the first intrinsic mode function is removed to attenuate the noise. The EMD method is also a single-trace method. Figure 7c shows the result from the curvelet method, we preserves the 5% largest coefficients in the curvelet domain to attenuate the noise. For the proposed method, as mentioned previously, we do not specify any parameters. From the denoised results, it is straightforward to conclude that both FX and SSA methods damagethe FX method damages a lot of energy around the edges. It is clear that the discontinuities in Figure 6 are smeared a lot. Since the FX and SSA methods are bothmethod is based on the spatial coherency

ACCEPTED MANUSCRIPT of the seismic events, they it will inevitably make the denoised data spatially smoothed. However, both curvelet method and the proposed method seem to obtain a successful preservations of the edges. For computational cost comparison, the FX method takes 0.04s, the EMD method takes 0.55s, the curvelet method takes 0.15s, and the proposed algorithm takes 1.18s. This example

T

contains 501 time samples and 64 spatial traces. The computation is done on a PC station equipped

IP

with an Intel Core i7 CPU clocked at 3.1 GHz and 16 GB of RAM. It is obvious that although the

CR

proposed obtains the best performance, it takes the longest computing time. The comparison of denoised data in the spectra is shown in Figure 8. It is obvious that the multichannel methods

US

remove more noise than the single-channel methods. But the proposed method removes more

AN

noise than its competitive method, i.e., the EMD method.

To compare the denoising performance, it is useful to compare the removed noise. The

M

removed noise is calculated by subtracting the denoised data (each subfigure in Figure 7) from the

ED

noisy data (in Figure 6b). Ideally, the removed noise should only contain random noise. If there is spatially coherent energy in the removed noise, it indicates that the denoising method damages

PT

some useful energy. We show the removed noise corresponding to four different methods in

CE

Figure 9. Figure 9a corresponds to the FX method, which contains significant coherent energy. Figure 9b shows the noise from the SSAEMD method, which also contains a lot of coherent

AC

energy, especially around those edge positions. The damage of EMD method may come from the fact that around the edge positions in this example, the waveform is more complicated and the EMD method is difficult to separate the signal and noise. Figure 9c shows the noise from the curvelet method. It is clear that the coherent noise is relatively weak compared with the FX and SSAEMD methods, but the low- frequency spatially coherent events are still visible, which indicates the damages to the useful signals using the curvelet method. Figure 9d, however, does

ACCEPTED MANUSCRIPT not contain any coherent energy. From this comparison, we conclude that only the proposed method among the four aforementioned methods can preserve all the useful signals. To make a quantitative comparison, we calculate the SNRs corresponding to the FX, SSAEMD, curvelet, and the proposed methods as 14.87 dB, 13.369.37 dB, 14.98 dB, and 16.32 dB, respectively. Note that

T

the SNR of the noisy data (shown in Figure 6b) is 6.49 dB. We find that the proposed method

IP

obtains the largest SNR improvement. It is also worth mentioning that all three methods other than

CR

the proposed method are tuned carefully to ensure the best performance to be presented. To further examine the performance, our proposed method is applied to a post-stack field

US

seismic data set, which consists of 219 traces and 1251 samples with a sample interval of 4 ms, as

AN

shown in Figure 10a. The data set mainly includes discontinuous events, nonstationary events, and faults. In other words, the non- linear and non-stationary signals in this data set are widespread.

M

From Figure 10a, it is obvious that the useful signal is blurred by noise. The signal is no longer

ED

mapped to a superposition of simple harmonics but rather a superposition of non- linear and non-stationary ones. We apply the four methods, namely, the FX, SSA,EMD, curvelet, and the

PT

proposed methods to the field data, and obtain the four denoised data, which are shown in Figure

CE

11. All four methods seem to obtain much improved results in terms of noise level, but. However, both SSAFX method and the curvelet method smooth too much the events and some details are no

AC

longer existing in the denoised data. The comparison in spectra is shown in Figure 12, where we can see clearly that the proposed method removes a little less noise than the multi-channel methods, e.g., the FX and curvelet methods, but removes more noise than the other single-channel method, i.e., the EMD method. When compared in the removed noise sections, as shown in Figure 13, the difference among all four methods becomes much more apparent. The SSAFX method and the curvelet

ACCEPTED MANUSCRIPT method, damage a significant amount of useful signals, as highlighted by the red frame boxes. The FXEMD method and the proposed method seem to obtain a good compromise between signal preservation and noise removal. Nevertheless, when observed carefully, there are still some signals in the noise from the FXEMD method. To compare the performance in details, we zoomed

T

a part from the noise sections and show the zoom- in comparison among different methods in

IP

Figure 14. The zooming areas correspond to the red frame boxes are shown in Figure 13. From

CR

Figure14, it is very obvious that only the proposed method preserves well the useful reflection signals well. Note that all the plots use the same amplitude clip value. As a complement, we also

US

show a comparison of the zoomed denoised data in the same area in Figure 15.

AN

The proposed method is only applicable to post-stack data, but also can obtain satisfactory results for pre-stack data. Figure 16(a) shows a realistic pre-stack data. Figure 16(b) shows the

M

denoised result using the proposed algorithm. It is salient that the reflection events are clearer and

PT

ED

more spatially coherent after the proposed denoising step.

CONCLUSIONS

CE

We presenthave proposed a new denoising algorithm for enhancing the reflection signals of seismic data. The advantage of the proposed method is that it can be automatically implemented,

AC

with no need for human interference. The adaptivity of the proposed method is partially due to the adaptive property of the empirical wavelet transform (EWT). To further reduceattenuate the residual noise, a clustering guided percentile thresholding step is also proposed. We found that the clustering method can effectively detect the positions of most signals and leave other components damped withby a percentile thresholding. The proposed method is a single-trace method and thus it preserves the spatial discontinuity well. It also has the potential to be used to denoise

ACCEPTED MANUSCRIPT microseismic data. The synthetic and field data examples show that the proposed method can successfully remove a significant amount of noise without damaging the useful signals.

ACKNOWLEDGEMENTS

T

This work was supported by the National Natural Science Foundation of China (Grant No.

AC

CE

PT

ED

M

AN

US

CR

IP

41804140).

ACCEPTED MANUSCRIPT References Asgedom, E. G., O. C. Orji, and W. Soellner, 2017, Rough-sea deghosting of single-sensor seismic data using the knowledge of the sea surface shape: Journal of Seismic Exploration, 26, 105–123.

T

Bai, M., X. Chen, J. Wu, Y. Chen, G. Liu, and E. Wang, 2016a, Multiple-component gaussian

CR

Geophysics-Chinese Edition, 59, no. 9, 3379–3393.

IP

beam reverse-time migration based on attenuation compensation: Chinese Journal of

Bai, M., X. Chen, J. Wu, G. Liu, Y. Chen, H. Chen, and Q. Li, 2016b, Q-compensated migration

US

by gaussian beam summation method: Journal of Geophysics and Engineering, 13, no. 1,

AN

35–48.

Bai, M., and J.Wu, 2017, Efficient deblending using median filtering without correct normal

M

moveout - with comparison on migrated images: Journal of Seismic Exploration, 26,

ED

455–479.

——, 2018, Seismic deconvolution using iteartive transform-domain sparse inversion: Journal of

PT

Seismic Exploration, 27, no. 2, 103–116.

CE

Bai, M., J. Wu, J. Xie, and D. Zhang, 2018, Least-squares reverse time migration of blended data with low-rank constraint along structural direction: Journal of Seismic Exploration, 27, no.

AC

1, 29–48.

Bezdek, J. C., 1981, Pattern recognition with fuzzy objective function algoritms: Plenum Press, New York. Canales, L., 1984a, Random noise reduction: 54th Annual International Meeting, SEG, Expanded Abstracts, 525–527.

ACCEPTED MANUSCRIPT ——, 1984b, Random noise reduction: 54th Annual International Meeting, SEG, Expanded Abstracts, 525–527. Candès, E. J., and Y. Plan, 2010, A probabilistic and ripless theory of compressed sensing: IEEE Transactions on Information Theory, 57, 7235–7254.

T

Candès, E. J., J. Romberg, and T. Tao, 2006, Robust uncertainty principles: Exact signal

IP

reconstruction from highly incomplete frequency information: IEEE Transactions on

CR

Information Theory, 52, 489–509.

Cao, J., Y. Wang, and B. Wang, 2014, Accelerating seismic interpolation with a gradient

US

projection method based on tight frame property of curvelet: Exploration Geophysics, 46,

AN

253–260.

Chen, Y., 2016, Dip-separated structural filtering using seislet thresholding and adaptive empirical

M

mode decomposition based dip filter: Geophysical Journal International, 206, no. 1,

ED

457–469.

——, 2017, Fast dictionary learning for noise attenuation of multidimensional seismic data:

PT

Geophysical Journal International, 209, no. 1, 21–31.

CE

——, 2018a, Automatic microseismic event picking via unsupervised machine learning: Geophysical Journal International, 212, no. 1, 88–102.

AC

——, 2018b, Automatic velocity analysis using high-resolution hyperbolic radon transform: Geophysis, 83, no. 4, A53–A57. ——, 2018c, Non-stationary least-squares complex decomposition for microseismic noise attenuation: Geophysical Journal International, 213, no. 3, 1572–1585.

ACCEPTED MANUSCRIPT Chen, Y., H. Chen, K. Xiang, and X. Chen, 2016a, Geological structure guided well log interpolation for high- fidelity full waveform inversion: Geophysical Journal International, 207, no. 2, 1313–1331. ——, 2017, Preserving the discontinuities in least-squares reverse time migration of

T

simultaneous-source data: Geophysics, 82, no. 3, S185–S196.

IP

Chen, Y., X. Chen, Y. Wang, and S. Zu, 2019, The interpolation of sparse geophysical data:

CR

Surveys in Geophysics, doi: 10.1007/s10712–018–9501–3.

Chen, Y., and S. Fomel, 2015, Random noise attenuation using local signal-and-noise

US

orthogonalization: Geophysics, 80, WD1–Wd9.

AN

——, 2018, EMD-seislet transform: Geophysics, 83, no. 1, A27–A32. Chen, Y., S. Fomel, and J. Hu, 2014, Iterative deblending of simultaneous-source seismic data

M

using seislet-domain shaping regularization: Geophysics, 79, no. 5, V179–V189.

ED

Chen, Y., W. Huang, Y. Zhou, W. Liu, and D. Zhang, 2018, Plane-wave orthogonal polynomial

214, 2207–2223.

PT

transform for amplitude-preserving noise attenuation: Geophysical Journal International,

CE

Chen, Y., and Z. Jin, 2016, Simultaneously removing noise and increasing resolution of seismic data using waveform shaping: IEEE Geoscience and Remote Sensing Letters, 13,

AC

102–104.

Chen, Y., and J. Ma, 2014, Random noise attenuation by f -x empirical mode decomposition predictive filtering: Geophysics, 79, V81–V91. Chen, Y., J. Ma, and S. Fomel, 2016b, Double-sparsity dictionary for seismic noise attenuation: Geophysics, 81, V17–V30.

ACCEPTED MANUSCRIPT Chen, Y., D. Zhang, Z. Jin, X. Chen, S. Zu, W. Huang, and S. Gan, 2016c, Simultaneous denoising and reconstruction of 5D seismic data via damped rank-reduction method: Geophysical Journal International, 206, 1695–1717. Claerbout, J. F., 1976, Fundamentals of geophysical data processing: Blackwell Scientific

T

Publications.

IP

Deighan, A. J., and D. R.Watts, 1997, Ground-roll suppression using the wavelet transform:

CR

Geophysics, 62, 1896–1903.

Dunn, J. C., 1973, A fuzzy relative of the isodata process and its use in detecting compact

US

well-separated clusters: Journal of Cybernetics, 3, 32–57.

AN

Ebrahimi, S., A. R. Kahoo, Y. Chen, and M. J. Porsani, 2017, A high-resolution weighted AB semblance for dealing with amplitude-variation-with-offset phenomenon: Geophysics, 82,

M

no. 2, V85–V93.

ED

Fomel, S., and Y. Liu, 2010, Seislet transform and seislet frame: Geophysics, 75, V25–V38. Gan, S., Y. Chen, S. Wang, X. Chen, W. Huang, and H. Chen, 2016a, Compressive sensing for

PT

seismic data reconstruction using a fast projection onto convex sets algorithm based on the

CE

seislet transform: Journal of Applied Geophysics, 130, 194–208. Gan, S., Y. Chen, S. Zu, S. Qu, and W. Zhong, 2015a, Structure-oriented singular value

AC

decomposition for signal enhancement of seismic data: Journal of Geophysics and Engineering, 12, 262–272. Gan, S., S. Wang, Y. Chen, J. Chen, W. Zhong, and C. Zhang, 2016b, Improved random-noise attenuation using f -x empirical mode decomposition and local similarity: Applied Geophysics, 13, 127–134.

ACCEPTED MANUSCRIPT Gan, S., S. Wang, Y. Chen, and X. Chen, 2016c, Simultaneous-source separation using iterative seislet-frame thresholding: IEEE Geoscience and Remote Sensing Letters, 13, no. 2, 197–201. Gan, S., S. Wang, Y. Chen, X. Chen, and K. Xiang, 2016d, Separation of simultaneous sources

T

using a structural-oriented median filter in the attened dimension: Computers &

IP

Geosciences, 86, 46–54.

CR

Gan, S., S. Wang, Y. Chen, S. Qu, and S. Zu, 2016e, Velocity a nalysis of simultaneous-source data using high-resolution semblance-coping with the strong noise: Geophysical Journal

US

International, 204, 768–779.

AN

Gan, S., S. Wang, Y. Chen, Y. Zhang, and Z. Jin, 2015b, Dealiased seismic data interpolation

Sensing Letters, 12, 2150–2154.

M

using seislet transform with low-frequency constraint: IEEE Geoscience and Remote

ED

Gilles, J., 2013, Empirical wavelet transform: IEEE Transactions on Signal Processing, 61, 3999–4010.

PT

Guo, P., and G. A. Mcmechan, 2017, Evaluation of three first-order isotropic viscoelastic

CE

formulations based on the generalized standard linear solid: Journal of Seismic Exploration, 26, 199–226.

AC

Guo, T., H. Chen, G.Wang, H. Wang, M. Cheng, and P. Chen, 2017, Seismic characteristics of complex halo-anhydrites and their effects on the underlying carbonatites: A case study of the right bank of the amu darya river: Journal of Seismic Exploration, 26, 381–397. Han, J., and M. van der Baan, 2015, Microseismic and seismic denoising via ensemble empirical mode decomposition and adaptive thresholding: Geophysics, 80, KS69–KS80.

ACCEPTED MANUSCRIPT Huang, N. E., Z. Shen, S. R. Long, M. C. Wu, H. H. Shih, Q. Zheng, N.-C. Yen, C. C. Tung, and H. H. Liu, 1998, The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis: Proceeding of the Royal Society of London Series A, 454, 903–995.

T

Huang, W., R. Wang, and Y. Chen, 2018, Regularized non-stationary morphological

IP

reconstruction algorithm for weak signal detection in micro-seismic monitoring:

CR

Methodology: Geophysical Journal International, 213, no. 2, 1189–1211. Huang, W., R. Wang, Y. Yuan, S. Gan, and Y. Chen, 2017a, Signal extraction using

US

randomized-order multichannel singular spectrum analysis: Geophysics, 82, no. 2,

AN

V59–V74.

Huang, W., R. Wang, S. Zu, and Y. Chen, 2017b, Low- frequency noise attenuation in seismic and

ED

International, 211, 1318–1340.

M

microseismic data using mathematical morphological filtering: Geophysical Journal

Jun, H., H. Jin, and C. Shin, 2017, Application of efficient frequency-domain full waveform

PT

inversion using time-domain encoded simultaneous sources: Journal of Seismic

CE

Exploration, 26, 141–169.

Kong, D., and Z. Peng, 2015, Seismic random noise attenuation using shearlet and total

AC

generalized variation: Journal of Geophysics and Engineering, 12, no. 12, 1024–1035. Li, H., R. Wang, S. Cao, Y. Chen, and W. Huang, 2016a, A method for low-frequency noise suppression based on mathematical morphology in microseismic monitoring: Geophysics, 81, no. 3, V159–V167.

ACCEPTED MANUSCRIPT Li, H., R. Wang, S. Cao, Y. Chen, N. Tian, and X. Chen, 2016b, Weak signal detection using multiscale morphology in microseismic monitoring: Journal of Applied Geophysics, 133, 39–49. Liu, C., D. Wang, B. Hu, and T. Wang, 2016a, Seismic deconvolution with shearlet sparsity

T

constrained inversion: Journal of Seismic Exploration, 25, no. 5, 433–445.

IP

Liu, G., and X. Chen, 2013, Noncausal f-x-y regularized nonstationary prediction filtering for

CR

random noise attenuation on 3D seismic data: Journal of Applied Geophysics, 93, 60–66.

US

Liu, G., X. Chen, J. Du, and K. Wu, 2012, Random noise attenuation using f -x regularized nonstationary autoregression: Geophysics, 77, V61–V69.

AN

Liu, G., S. Fomel, L. Jin, and X. Chen, 2009, Stacking seismic data using local correlation: Geophysics, 74, V43–V48.

M

Liu, W., S. Cao, and Y. Chen, 2016b, Applications of variational mode decomposition in seismic

ED

time-frequency analysis: Geophysics, 81, no. 5, V365–V378.

PT

——, 2016c, Seismic time-frequency analysis via empirical wavelet transform: IEEE Geo-science and Remote Sensing Letters, 13, 28–32.

CE

Liu, W., S. Cao, Y. Chen, and S. Zu, 2016d, An effective approach to attenuate random noise based on compressive sensing and curvelet transform: Journal of Geophysics and

AC

Engineering, 13, no. 2, 135–145. Liu, W., S. Cao, S. Gan, Y. Chen, S. Zu, and Z. Jin, 2016e, One-step slope estimation for dealiased seismic data reconstruction via iterative seislet thresholding: IEEE Geoscience and Remote Sensing Letters, 13, 1462–1466. Liu, W., S. Cao, Z. Jin, Z. Wang, and Y. Chen, 2018a, A novel hydrocarbon detection approach via high-resolution frequency-dependent avo inversion based on variational mode

ACCEPTED MANUSCRIPT decomposition: IEEE Transactions on Geoscience and Remote Sensing, 56, no. 4, 2007–2024. Liu, W., S. Cao, Y. Liu, and Y. Chen, 2016f, Synchrosqueezing transform and its applications in seismic data analysis: Journal of Seismic Exploration, 25, 27.

T

Liu, W., S. Cao, Z. Wang, K. Jiang, Q. Zhang, and Y. Chen, 2018b, A novel approach for seismic

IP

time-frequency analysis based on high-order synchrosqueezing transform: IEEE

CR

Geoscience and Remote Sensing Letters, 15, no. 8, 1159–1163.

Liu, W., S. Cao, Z. Wang, X. Kong, and Y. Chen, 2017, Spectral decomposition for hydrocarbon

US

detection based on vmd and teager-kaiser energy: IEEE Geoscience and Remote Sensing

AN

Letters, 14, no. 4, 539–543.

Liu, Y., and S. Fomel, 2010, Oc-seislet: Seislet transform construction with differential offset

M

continuation: Geophysics, 75, WB235–WB245.

ED

Mallat, S. G., 2009, A wavelet tour of signal processing: The sparse way: Academic Press. Nazari Siahsar, M. A., S. Gholtashi, A. R. Kahoo, H. Marvi, and A. Ahmadifard, 2016, Sparse

PT

time-frequency representation for seismic noise reduction using low-rank and sparse

CE

decomposition: Geophysics, 81, V117–V124. Qu, S., H. Zhou, Y. Chen, S. Yu, H. Zhang, J. Yuan, Y. Yang, and M. Qin, 2015, An effective

AC

method for reducing harmonic distortion in correlated vibroseis data: Journal of Applied Geophysics, 115, 120–128. Qu, S., H. Zhou, R. Liu, Y. Chen, S. Zu, S. Yu, J. Yuan, and Y. Yang, 2016, Deblending of simultaneous-source seismic data using fast iterative shrinkage-thresholding algorithm with firm-thresholding: Acta Geophysica, 64, no. 4, 1064–1092.

ACCEPTED MANUSCRIPT Ryu, D., A. Kim, W. Ha, and C. Shin, 2017, Robustness of laplace domain waveform inversions to cycle skipping: Journal of Seismic Exploration, 26, 251–266. Shahin, A., M. Myers, and L. Hathon, 2017, Carbonates' dual-physics modeling aimed at seismic reservoir characterization: Journal of Seismic Exploration, 26, 331–349.

T

Siahsar, M. A. N., V. Abolghasemi, and Y. Chen, 2017, Simultaneous denoising and interpolation

IP

of 2d seismic data using data-driven non-negative dictionary learning: Signal Processing,

CR

141, 309–321.

Tian, Y., Y. Li, H. Lin, and H. Ma, 2013, A sparse nmf-su for seismic random noise attenuation:

US

IEEE Geoscience & Remote Sensing Letters, 10, 607–611.

AN

Tufts, D., and R. Kumaresan, 1982, Estimation of frequencies of multiple sinusoids: Making linear prediction perform like maximum likelihood: Proceedings of the IEEE, 70, 975–989.

M

Wang, H., Z. Jin, D. Tian, W. Xu, and Y. Tan, 2017a, An improved q-tomography method for

ED

estimation of seismic attenuation through non-linear inversion: Journal of Seismic Exploration, 26, 171–182.

PT

Wang, R., and Y. Wang, 2017, Seismic reectivity inversion by curvelet deconvolution–a

331–349.

CE

comparative study and further improvements: Journal of Seismic Exploration, 26,

AC

Wang, Y., X. Ma, H. Zhou, and Y. Chen, 2018, L1-2 minimization for exact and stable seismic attenuation compensation: Geophysical Journal International, 213, no. 3, 1629–1646. Wang, Y., H. Zhou, S. Zu, W. Mao, and Y. Chen, 2017b, Three-operator proximal splitting scheme for 3D seismic data reconstruction: IEEE Geoscience and Remote Sensing Letters, 14, no. 10, 1830–1834. Welland, G., 2003, Beyond wavelets: Elsevier, 10.

ACCEPTED MANUSCRIPT Wu, J., and M. Bai, 2018a, Attenuating seismic noise via incoherent dictionary learning: Journal of Geophysics and Engineering, 15, doi: 10.1088/1742–2140/aaaf57. ——, 2018b, Incoherent dictionary learning for reducing crosstalk noise in least-squares reverse time migration: Computers and Geosciences, 114, 11–21.

T

Wu, J., R. Wang, Y. Chen, Y. Zhang, S. Gan, and C. Zhou, 2016, Multiples attenuation using

IP

shaping regularization with seislet domain sparsity constraint: Journal of Seismic

CR

Exploration, 25, 1–9.

Wu, Z., and N. Huang, 2009, Ensemble empirical mode decomposition: A noise-assisted data

US

analysis method: Advances in Adaptive Data Analysis, 1, 1–41.

AN

Xie, J., B. Di, D. Schmitt, J. Wei, and Y. Chen, 2018, Estimation of δ and c13 of organic-rich shale from laser ultrasonic technique (LUT) measurement: Geophysis, 83, no. 4, C137–C152.

M

Xue, Y., F. Chang, D. Zhang, and Y. Chen, 2016a, Simultaneous sources separation via an

ED

iterative rank-increasing method: IEEE Geoscience and Remote Sensing Letters, 13, no.

PT

12, 1915 – 1919.

Xue, Y., M. Man, S. Zu, F. Chang, and Y. Chen, 2017, Amplitude-preserving iterative deblending

CE

of simultaneous source seismic data using high-order radon transform: Journal of Applied

AC

Geophysics, 139, 79–90. Xue, Y., J. Yang, J. Ma, and Y. Chen, 2016b, Amplitude-preserving nonlinear adaptive multiple attenuation using the high-order sparse radon transform: Journal of Geophysics and Engineering, 13, 207–219. Yang, W., R. Wang, J. Wu, Y. Chen, S. Gan, and W. Zhong, 2015, An efficient and effective common reection surface stacking approach using local similarity and plane-wave attening: Journal of Applied Geophysics, 117, 67–72.

ACCEPTED MANUSCRIPT Yang, Y., D. Li, T. Tong, D. Zhang, Y. Zhou, and Y. Chen, 2018, Denoising controlled-source electromagnetic data using least-squares inversion: Geophysis, 83, no. 4, E229–E244. Yilmaz, O., 1987, Seismic data processing: Society of Exploration Geophysicists. Zeng, Q., Y. Guo, R. Jiang, J. Ba, H. Ma, and J. Liu, 2017, Fluid sensitivity of rock physics

T

parameters in reservoirs: Quantitative analysis: Journal of Seismic Exploration, 26,

IP

125–140.

CR

Zhang, D., Y. Chen, W. Huang, and S. Gan, 2016a, Multi-step damped multichannel singular spectrum analysis for simultaneous reconstruction and denoising of 3D seismic data:

US

Journal of Geophysics and Engineering, 13, 704–720.

AN

Zhang, G., Z. Wang, and Y. Chen, 2018a, Deep learning for seismic lithology prediction: Geophysical Journal International, 215, 1368–1387.

M

Zhang, H., Q. Liu, and J. Hao, 2017a, Least-squares reverse-time migration toward “true"

ED

reectivity: Journal of Seismic Exploration, 26, 183–198. Zhang, P., Y. Dai, R. Wang, and Y. Tan, 2017b, A quantitative evaluation method based on emd

PT

for determining the accuracy of time-varying seismic wavelet extraction: Journal of

CE

Seismic Exploration, 26, 267–292. Zhang, Q., Y. Chen, H. Guan, and J. Wen, 2016b, Well-log constrained inversion for lithology

AC

characterization: a case study at the jz25-1 oil field, china: Journal of Seismic Exploration, 25, 121–129.

Zhang, Q., W. Mao, and Y. Chen, 2018b, Attenuating crosstalk noise of simultaneous-source least-squares reverse time migration with gpu-based excitation-amplitude imaging condition: IEEE Transaction on Geosciences and Remote Sensing, in press.

ACCEPTED MANUSCRIPT Zhang, Q., W. Mao, H. Zhou, H. Zhang, and Y. Chen, 2018c, Hybrid-domain simultaneous-source full waveform inversion without crosstalk noise: Geophysical Journal International, 215, no. 3, 1659–1681. Zhong, W., Y. Chen, S. Gan, and J. Yuan, 2016, L1/2 norm regularization for 3D seismic data

T

interpolation: Journal of Seismic Exploration, 257–268.

IP

Zhou, Y., 2017, A POCS method for iterative deblending constrained by a blending mask: Journal

CR

of Applied Geophysics, 138, 245–254.

Seismic Exploration, 27, no. 1, 69–88.

US

Zhou, Y., and W. Han, 2018, Multiples attenuation in the presence of blending noise: Journal of

AN

Zhou, Y., and S. Li, 2018, Simultaneous deblending and interpolation using structure-oriented filters: Journal of Applied Geophysics, 150, 230–243.

M

Zhou, Y., S. Li, J. Xie, D. Zhang, and Y. Chen, 2017a, Sparse dictionary learning for seismic noise

ED

attenuation using a fast orthogonal matching pursuit algorithm: Journal of Seismic Exploration, 26, no. 5, 433–454.

PT

Zhou, Y., S. Li, D. Zhang, and Y. Chen, 2018, Seismic noise attenuation using an online subspace

CE

tracking algorithm: Geophysical Journal International, 212, no. 2, 1072–1097. Zhou, Y., C. Shi, H. Chen, J. Xie, G. Wu, and Y. Chen, 2017b, Spike-like blending noise

AC

attenuation using structural low-rank decomposition: IEEE Geoscience and Remote Sensing Letters, 14, no. 9, 1633–1637. Zhou, Y., and S. Zhang, 2017, Robust noise attenuation based on nuclear norm minimization and a trace prediction strategy: Journal of Applied Geophysics, 147, 52–67.

ACCEPTED MANUSCRIPT Zu, S., H. Zhou, H. Chen, H. Zheng, and Y. Chen, 2017a, Two field trials for deblending of simultaneous source surveys: why we failed and why we succeeded?: Journal of Applied Geophysics, 143, 182–194. Zu, S., H. Zhou, Y. Chen, , S. Qu, X. Zou, H. Chen, and R. Liu, 2016a, A periodically varying code

T

for improving deblending of simultaneous sources in marine acquisition: Geophysics, 81,

IP

V213–V225.

CR

Zu, S., H. Zhou, Y. Chen, X. Pan, S. Gan, and D. Zhang, 2016b, Interpolating big gaps using inversion with slope constraint: IEEE Geoscience and Remote Sensing Letters, 13,

US

1369–1373.

AN

Zu, S., H. Zhou, Q. Li, H. Chen, Q. Zhang, W. Mao, and Y. Chen, 2017b, Shot-domain deblending using least-squares inversion: Geophysics, 82, no. 4, V241–V256.

M

Zu, S., H. Zhou, W. Mao, F. Gong, and W. Huang, 2018, 3D deblending of simultaneous source

ED

data based on 3D multi-scale shaping operator: Journal of Applied Geophysics, 151, 274–289.

PT

Zu, S., H. Zhou, W. Mao, D. Zhang, C. Li, X. Pan, and Y. Chen, 2017c, Iterative deblending of

CE

simultaneous-source data using a coherency-pass shaping operator: Geophysical Journal

AC

International, 211, no. 1, 541–557.

ACCEPTED MANUSCRIPT APPENDIX A INTRODUCTION OF THE COMPARING METHODS FX method f  x deconvolution is one of the most widely used approaches for random noise attenuation. Let

T

st ( x) be the seismic signal located at trace x and time t . If the slope of a linear event with

IP

constant amplitude in a seismic section is  , then:

CR

st ( x  1) = st  xh (1),

(A-1)

US

where h denotes the spatial interval. Equation A-1 can be directly transformed into the frequency domain utilizing the time-shift property of the Fourier transform:

AN

S f ( x  1) = S f (1)ei 2 fxh .

(A-2)

M

For a specific frequency f 0 , from equation A-2 we can obtain a linear recursion, which is given

ED

by:

S f ( x  1) = a f (1) S f ( x),

where a f (1) = e

 i 2 f0h

0

0

(A-3)

0

PT

0

. This recursion is also known as an auto-regressive (AR) model of order

CE

1 (Canales, 1984b). Similarly, superposition of p linear events in the t  x domain can be

AC

represented by an AR model of order p (Tufs and Kumaresan, 1982) as the following equation: S f ( x  1) = a f (1)S f ( x)  a f (2)S f ( x  1)  0

0

0

0

0

 a f ( p)S f ( x  1  p), 0

0

(A-4)

Equation A-4 can be formulated as a convolutional form: (A-5)

d = f * a,

where d denotes the vector composed of S f ( x  1)( x = 1, 2, 0

composed of S f ( x)( x = 1, 2, 0

, X ) , f denotes the vector

, X ) , a denotes the vector composed of a f ( x)( x = 1, 2, 0

, X) ,

ACCEPTED MANUSCRIPT and X denotes the number of traces. Equation A-5 can be formulated as a matrix vector form: (A-6)

d = Fm,

where F is the convolution matrix composed by f .

T

However, equation A-6 is based on clean signal model. In reality, the seismic data is

IP

composed of random noise. Thus, we have to solve a from the noise corrupted observation d

CR

based on some optimization schemes. Based on equation A-4, we can formulate an optimization

US

problem based on the minimum prediction error energy assumption. The predictive error filter can be solved by minimizing the following objective function:



2 2

denotes the squares of L2 norm.

M

where

(A-7)

AN

J = Fm  d 22 ,

Taking derivatives of the cost function A-7 with respect to m , and setting the result to

ED

zero, we can obtain the following equation:

(A-8)

PT

FT d = FT Fm,

CE

where []T denotes transpose. Note that FT F is a Toeplitz form and thus can be efficiently solved using Levinson’s recursion. In order to stabilize the recursion for solving aˆ , we need to

AC

add a small perturbation to the diagonal of the Toeplitz matrix:

aˆ = (FT F  I)1 FT d.

(A-9)

Finally, the estimated clean data (denoised data) can be expressed as:

dˆ = Faˆ .

EMD method

(A-10)

ACCEPTED MANUSCRIPT The aim of empirical mode decomposition (EMD) is to empirically decompose a non-stationalnon-stationary signal into a finite set of subsignals, which are calledtermed intrinsic mode functions (IMF) and are considered to be stable. The IMFs satisfy two cond itions: (1) in the whole data set, the number of extrema and the number of zero crossings must either equal or differ

T

at most by one; and (2) at any point, the mean value of the envelope defined by the local maxima

IP

and the envelope defined by the local minima is zero (Huang et al., 1998).

CR

Provided that s(t ) , cn (t ) , r (t ) , and N denote the original non-stationalnon-stationary

US

signal, the separated Intrisic Mode Function (IMF)IMFs, the residual, and the number of IMFs, respectively., Tthe mathematical principle of EMD can be expressed as: N

n =1

AN

s(t ) = cn (t )  r (t ).

(A-11)

, N ).

ED

subsignals cn (t ) ,( n = 1, 2,

M

For a non-stationalnon-stationary signal s(t ) , using equation A-11, we get a finite set of

PT

A special property of IMFEMD is that the IMFs represents different oscialla tions embedded in the data, andwhere the oscillating frequency for each subsignal cn (t ) decreases as

CE

the sequence number of the IMF becomes larger (we call it a frequency decreasing property in the

AC

following context). This property results from the sifting algorithm used to implement the decomposition. For non-stationary seismic data, we can assume the noise lays in the first component after EMD decomposition. Since EMD is an adaptive method, the resulting EMD method is also automatic.

Curvelet method Under a tight frame, the forward curvelet transform can be simply described as the inner product

ACCEPTED MANUSCRIPT between the signal f and curvelets  j ,l ,k

c( j, l , k ) =< f ,  j ,l ,k >=  2f (x) j ,l ,k (x)dx

(A-12)

R

where, f denotes the signal,  j ,l ,k is curvelets, c( j, l , k ) is the curvelet coefficient, j, l , k

T

represent scale, orientation and position, respectively.

IP

Thus, the curvelets after dilating, translation and rotation can produce the tight frame in

A-13:

US

f'(x) = c( j, l , k )   j ,l ,k

CR

square- integrable space, which implies the reconstruction equation can be written as equation

j ,l , k

(A-13)

AN

where  j ,l ,k denotes inverse curvelets, f' represents inverse curvelet transform.

M

In the curvelet method, a forward curvelet transform is first applied to the raw seismic data, a thresholding operation is then applied to the coefficients, finally an inverse curvelet transform is

ED

applied to the transform domain to obtain the denoised data. In this paper, we use a intelligent

PT

threshold value selection strategy called percentile thresholding for the curvelet method. We preserve a constant percentage of the largest coefficients to achieve the thresholding. Thus, the

CE

threshold value is intelligently adjusted by the percentage we use. For example, we can select a

AC

percentage less than 5% to reject the small-amplitude coefficients.

ACCEPTED MANUSCRIPT Fig. 1 Single-trace denoising example. (a) Clean data. (b) Noisy data. (c) Denoised data using EWT decomposition. (d) Removed noise using EWT decomposition. Fig. 2 Spectra of the single-trace denoising example. (a) Spectrum of the clean data. (b) Spectrum of the noisy data. (c) Spectrum of the denoised data using EWT decomposition. (d)

T

Spectrum of the removed noise using EWT decomposition.

IP

Fig. 3 Decomposed results using EWT. (a)-(e): First to fifth components.

CR

Fig. 4 Clustered waveforms. The red line indicates the label for waveform clustering. Fig. 5 A comparison between denoising methods. (a) Clean data. (b) Noisy data. (c)

US

Denoised data using EWT decomposition. (d) Denoised data using the clustering guided EWT

AN

thresholding method.

Fig. 6 Synthetic example. (a) Clean data. (b) Noisy data. (c) Spectra of the clean data. (d)

M

Spectra of the noisy data.

ED

Fig. 7 Comparison of the denoised data for the synthetic example. (a) Denoised data using FX method. (b) Denoised data using SSA method.Denoised data using EMD method. (c) Denoised

PT

data using the curvelet method. (d) Denoised data using the proposed method.

CE

Fig. 8 Comparison of the denoised data in spectra. (a) Spectra of denoised data using FX method. (b) Spectra of denoised data using EMD method. (c) Spectra of denoised data using the

AC

curvelet method. (d) Spectra of denoised data using the proposed method. Fig. 9 Comparison of the removed noise for the synthetic example. (a) Removed noise using FX method. (b) Removed noise using SSA method.Removed noise using EMD method. (c) Removed noise using the curvelet method. (d) Removed noise using the proposed method. Fig. 10 Field data example. Fig. 11 Comparison of the denoised data. (a) Denoised data using FX method. (b)

ACCEPTED MANUSCRIPT Denoised data using SSA method.Denoised data using EMD method. (c) Denoised data using the curvelet method. (d) Denoised data using the proposed method. Fig. 12 Comparison of the denoised field data in spectra. (a) Spectra of denoised data using FX method. (b) Spectra of denoised data using EMD method. (c) Spectra of denoised data using

T

the curvelet method. (d) Spectra of denoised data using the proposed method.

IP

Fig. 13 Comparison of the removed noise. (a) Removed noise using FX method. (b)

CR

Removed noise using SSA method.Removed noise using EMD method. (c) Removed noise using the curvelet method. (d) Removed noise using the proposed method.

US

Fig. 14 Zoom- in comparison of the removed noise. (a) Removed noise using FX method.

AN

(b) Removed noise using SSA method.Removed noise using the EMD method. (c) Removed noise using the curvelet method. (d) Removed noise using the proposed method.

M

Fig. 15 Zoom- in comparison of the denoised data in the same area as in Figure 14. (a)

ED

Denoised data using FX method. (b) Denoised data using the EMD method. (c) Denoised data using the curvelet method. (d) Denoised data using the proposed method.

PT

Fig. 16 Denoising test for a pre-stack data. (a) Noisy data. (b) Denoised data using the

AC

CE

proposed method.

Figure 1

Figure 2

Figure 3

Figure 4

Figure 5

Figure 6

Figure 7

Figure 8

Figure 9

Figure 10

Figure 11

Figure 12

Figure 13

Figure 14

Figure 15

Figure 16