Seismic random noise attenuation using modified wavelet thresholding

In seismic exploration, random noise deteriorates the quality of acquired data. This study analyzed existing denoising methods used in seismic exploration from the perspective of random noise. Wavelet thresholding offers a new approach to reducing random noise in simulation results, synthetic data, and real data. A modified wavelet threshold function was developed by considering the merits and demerits of conventional soft and hard thresholding schemes. A MATLAB (matrix laboratory) simulation model was used to compare the signal-to-noise ratios (SNRs) and mean square errors (MSEs) of the soft, hard, and modified threshold functions. The results demonstrated that the modified threshold function can avoid the pseudo-Gibbs phenomenon and produce a higher SNR than the soft and hard threshold functions. A seismic convolution model was built using seismic wavelets to verify the effectiveness of different denoising methods. The model was used to demonstrate that the modified thresholding scheme can effectively reduce random noise in seismic data and retain the desired signal. The application of the proposed tool to a real raw seismogram recorded during a land seismic exploration experiment located in north China clearly demonstrated its efficiency for random noise


Introduction
Seismic exploration is the only way for petroleum companies to obtain firsthand geological data on oil and gas reservoirs.The process of identifying the geological structure of a reservoir involves three steps: generating seismic waves with artificial seismic sources, acquiring seismic data with wave detectors, and analyzing and processing the collected seismic data with a computer.In actual seismic exploration, the signals received by wave detectors are usually mixed with noise from the environment; this interference inevitably influences the signal analysis and processing, and the useful signal may be submerged in severe cases.Therefore, denoising is a key step to acquire reliable seismic exploratory information.
Noise can generally be divided into two categories: regular and irregular.Regular noise, which is also called coherent noise, possesses a certain dominant frequency and apparent velocity.Examples include acoustic waves, shallow refracted waves, multiple waves, and infrasonic waves.This type of noise can be suppressed by utilizing the differences in the frequency spectra and apparent velocity between the noise and desired useful signal.A representative method involves one-dimensional frequency-based filtering and two-dimensional apparent velocity domain filtering via the Fourier transform [Galiana-Merino et al. 2013].
Irregular noise, which is also called random noise, has no definite apparent velocity or dominant frequency; in seismic records, it appears disordered.Therefore, it is impossible to remove random noise from the signal by targeting specific differences in the spectra and apparent velocity along the propagating direction.The random noise seems irregular but conforms to its own statistical regularity.The main methods for reducing random noise include f − x domain predictive filtering, x − p transform [Kappus et al. 1990, Turner 1990, Zhou and Greenhalgh 1994], and K − L transform [Hunt and Kubler 1984].The second-order autoregressive (AR) model is used for the f − x predictive filtering technique, which uses the energy minimization principle to predict the error operator in the two-dimensional frequencyspatial domain.The x − p transform can distinguish the primary reflection from multiple reflections.The K − L transform is an optimal orthogonal transformation based on statistical signal processing; it is mostly used in clustering analysis, pattern recognition, optimization, and denoising.In the context of seismic exploration,
the K − L transform is used to reduce seismic noise after a suitable eigenvalue is selected; it can maintain signal coherence while denoising [Golestani et al. 2013, Xue et al. 2014].
During seismic exploration, the reflection seismogram is based on varying wavelets.Waves reflected from different interfaces arrive at the wave detectors at different times; they have different waveforms and frequency spectra, and the seismic record is a superposition of these seismic wavelets.Existing denoising techniques such as f − x filtering, x − p transform, and K − L transform cannot distinguish such wavelet shapes in the time or frequency domain.The difference between an interference wave and the useful signal in the time-scale domain can only be utilized for the discrimination if the seismic record is decomposed in the time-frequency domain.
Different frequency components have different signal-to-noise ratios (SNRs) and must be treated differently.Furthermore, different frequency components with different SNRs also have different effects on the resolution; a higher SNR helps improve the resolution.Therefore, seismic signal denoising should be done at different scales specific to different frequency components, or the denoising effect will likely be impaired.
The classical Fourier transform cannot perform analyses in the time and frequency domains simultaneously [Chakraborty and Okaya 1995].It is only capable of analysis in the frequency domain because the time characteristic of the signal disappears after the Fourier transform.The short-time Fourier transform can overcome the drawback of the classical Fourier transform, but it is limited to a single resolution.The modified wavelet transform overcomes this single resolution problem of the short-time Fourier transform.In addition, it provides two-dimensional resolution within a single frequency.Thus, it is an effective tool for analyzing the change in the frequency distribution of a signal over time [Bartosch and Seidl 1999].
The wavelet transform can provide good time and frequency localization as well as decompose signal information into extremely small parts.When a wavelet transform is applied to two-dimensional decomposition of a seismic record, the desired signal can be accurately separated from the noise based on differences in details in the wavelet domain.Moreover, the wavelet transform's multiscale resolution properties can help with identifying the contributions of different frequency components having different SNRs to the resolution [Singh and Khare 2014].Wide time windows are used to improve the frequency resolution of low-frequency components, and narrow time windows are used to increase the time resolution of high-frequency components.Because it can achieve both time and frequency resolutions, the wavelet transform has the potential to become a useful and powerful technique for seismic signal denoising.
In signal analysis, a wavelet transform is used for filtering, wave compression, and transmission [Zhu et al. 2009, Yang andRen 2011].In image processing, it is usually used for image compression, classification, identification, and diagnosis [Liu et al. 2012].However, research on applying the wavelet transform to seismic signal processing remains in the early stages.
One of the key steps of wavelet transform denoising is the choice of a suitable threshold function.The soft and hard threshold functions are commonly used.The hard threshold function is discontinuous at the threshold limit; a relatively violent oscillation occurs at the point of abrupt change when a wavelet coefficient is reconstructed.This is known as the pseudo-Gibbs phenomenon.Researchers tend to utilize the soft threshold function when applying a wavelet transform.Cao et al. constructed a second-generation wavelet by using a lifting scheme.By quantizing the threshold with a soft threshold function, they were able to separate the signal into two isolated sub-signals.They then obtained the approximate signal and high-frequency effective signal by performing this operation repeatedly [Cao and Chen 2005].
The undecimated discrete wavelet transform (UDWT) method uses soft thresholding for noise reduction and does not suffer from the limitations of the conventional discrete wavelet transform (DWT) [Goudarzi and Riahi 2012].Quadfeul and Aliouane combined a denoising algorithm based on the discrete wavelet transform decomposition with the continuous wavelet transform and used a level 5 Haar wavelet to attenuate random noise without using classical bandpass filters [Ouadfeul and Aliouane 2014].
It is worth noting that, when the absolute value of the wavelet coefficient is greater than the preset threshold value, the soft threshold function will entirely subtract this threshold value without thorough consideration.This results in the loss of useful data in the reconstructed signal.Consequently, the signal shows an extremely smooth waveform.Therefore, determining how to modify the threshold function by incorporating the merits of both soft and hard threshold functions is a very important issue to improving the SNR.
Based on a theoretical analysis of wavelet denoising, we propose a new method for wavelet thresholding.The numerical computing software MATLAB (matrix laboratory) was used to produce simulations of the Heavy sine signal, and a comparative analysis was performed to assess the feasibility of the new denoising method.We constructed a seismic convolution model to examine the effectiveness of this method for seis-mogram denoising.Finally, the proposed tool was applied to a real raw seismogram record to show its efficiency for random noise attenuation.

Wavelet thresholding
The wavelet transform decomposes the original wavelet into a series of wavelets that are shifted and scaled.After this wavelet transform, the wavelet coefficient corresponding to signals possessing important information has a relatively large amplitude, but such signals are few in number.On the other hand, the wavelet coefficient corresponding to noise has a small amplitude.Basic procedures for wavelet transform denoising include selecting the proper thresholds for different scale, threshold quantization (i.e., set the wavelet coefficients smaller than the threshold to zero and retaining those that are bigger than the threshold), reverse transformation on the remaining wavelets, and obtaining the reconstructed signals.
Wavelet thresholding is an important denoising method for the wavelet transform.Research has shown that wavelet thresholding provides a nearly optimal effect in terms of the least mean squares error [Donoho andJohnston 1994, Borcea et al. 2010].
Selecting a suitable threshold and threshold function are the most important steps in wavelet thresholding.Choosing a proper wavelet function for the decomposition and determining the number of decomposition layers are also important.Thus, a complete algorithm for wavelet thresholding can only be completely designed to realize the optimum denoising effect if these factors are perfectly coordinated with each other.
The wavelet transform can be expressed as where, W f (a,b) is the coefficient of the transformed signal f (t), a is the scaling factor used to generate fundamental waves with different frequencies, and b is the shift factor used to localize frequency spectra.Together, a and } determine both the position of the center of f (t) and the temporal width in the analysis.We choose the proper wavelet function and number of decomposition levels to decompose the noisy signal X for denoising.Figure 1 presents the schematic for wavelet decomposition at three levels.Ultimately, the high-frequency components of the signal dominated by noise are mostly contained in CD 1 , CD 2 , and CD 3 , while CA 1 , CA 2 , and CA 3 comprise the low-frequency components dominated by the desired signal and should be free from noise.The high-frequency coefficients CD 1 , CD 2 , and CD 3 are quantized by using the chosen threshold quantization rule.The inverse wavelet transform is then performed on the high-frequency coefficients from the first level to the third level and the low-frequency coefficients at the third level to yield the reconstructed signals.

Choosing the wavelet function
The wavelet functions provided by MATLAB include Morlet, Mexican hat, Meyer, Haar, Daubechies, Symlets, Coiflets, and Spline biorthogonal.If the wavelet functions are orthogonal and have scaling functions, the transformed signal have the same number of points as the original signal.This allows a fast algorithm to be used to reduce the amount of calculation.Wavelet functions that satisfy these requirements include Haar, Daubechies, Coiflets, and Symlets.
The Haar wavelet is discontinuous in the time domain and thus does not provide superior performance as a basic wavelet.
The Daubechies wavelet of order N (dbN) waveletis characterized by the fact that the order of vanishing moments increases with the order (sequence N).Larger vanishing moments indicate greater smoothness along with stronger localization of the frequency domain and better distinction among frequency bands.However, this weakens the compacted support in the time domain, which significantly increases the amount of calculation and degrades the real-time capability.Furthermore, the dbN wavelet is not symmetric (non-zero phase) except when N = 1, which results in phase distortion during the decomposition and reconstruction of the signal [Yi et al. 2012].
The wavelet function W(t) of Coiflet has a 2N moment equal to zero, just like the 2N − 1 moment of the scaling function {(t).Both have a support length of 6N − 1. W(t) and {(t) of Coiflet have better symmetry than dbN.
Symlet is an approximate symmetric wavelet function proposed by Daubechies as an improvement of the dbN wavelet function.The support length of a symN wavelet is 2N − 1, and the vanishing moment is N.The wavelet also has good regularity.The symN wavelet has Figure 1.Schematic of wavelet decomposition.
the same properties as the dbN wavelet in terms of the continuity, support length, and filter length but has better symmetry, which helps with attenuating the phase distortion during signal analysis and reconstruction [Garg 2016].Based on the above factors, the sym6 wavelet was selected in this study.

Choosing the threshold
A threshold is generally expressed as m = v√2ln(N), where v denotes the standard deviation of the noise and N denotes the number of sampling points [Azzalini et al. 2005].
If the standard deviation is unavailable in practical computation, it can be estimated by using the median absolute deviation (MAD): where x is the wavelet coefficient of the high-frequency components at the corresponding level of the wavelet decomposition and median is the median value of the wavelet coefficient series.
To avoid excess elimination of wavelet coefficients, the thresholding technique is modified by making use of the attenuation properties of noise in the wavelet transform domain.As shown in Equation (3), the modified threshold tends to gradually decrease with increasing j, which is consistent with the attenuation properties of noise: (3) Here, v is the standard deviation of noise, N is the number of sampling points, and j denotes the decomposition level.

Choosing the threshold function
Donoho proposed two threshold functions: hard and soft [To et al. 2009, Paul et al. 2013].In the soft threshold function, as shown in Equation ( 4), the threshold m is subtracted from a wavelet coefficient with an absolute value that is greater than m, or the wavelet coefficient is set to zero if its absolute value is less than m (Figure 2a): (4) In the hard threshold function, as shown in Equation (5), the wavelet coefficient is unaltered if its absolute value is greater than m or set to zero if its absolute value is less than m (Figure 2b): (5) where x is the wavelet coefficient, m is the threshold, and y is the quantized coefficient.
Wavelet coefficients are unbiased when filtered with the hard threshold function, but the function is discontinuous at the threshold points m and −m.During signal reconstruction, strong oscillations occur at the signal's inflection points; this is called the pseudo-Gibbs phenomenon.In contrast, the soft threshold function is continuous and therefore avoids the pseudo-Gibbs phenomenon during signal reconstruction.However, the soft threshold function is bound to cause a deviation in the estimate because it indiscriminately subtracts m from all wavelet coefficients whose absolute values are higher than the threshold.This is reflected by the fact that the signal reconstructed from wavelet coefficients using the soft threshold function is too smooth, which reduces the effectiveness of denoising.
For this reason, it is necessary to create a new threshold function that overcomes the drawbacks of the soft and hard threshold functions.The new threshold function needs to meet two conditions: the estimated coefficient value should approach zero as the wavelet coefficient's absolute value approaches m, and .as the wavelet coefficient's absolute value diverges from m, the coefficient should approach the value estimated by the hard threshold function.Only under these conditions can the pseudo-Gibbs phenomenon be avoided for better denoising effectiveness.Thus, we propose the following modified threshold function: (6 where m is the adjusting factor, m is the threshold, and x is the wavelet coefficient.
The above equation illustrates the following findings: -when m→ +∞, the new function is equivalent to the soft threshold function; when m→ 0, the new function is equivalent to the hard threshold function; -when |x|→ ∞, lim |x|→∞ y = sign(|x|)(x), and the function is equivalent to the hard threshold function; this means that the function y = x is an asymptote of the modified threshold function; -when |x|=m, solving the equation yields y = sign(x)(|x|− m) = 0, which indicates that the modified threshold function has discontinuity.In this manner, the pseudo-Gibbs phenomenon is eliminated.

Evaluating the denoising effectiveness
Common methods for evaluating the denoising effectiveness are divided into two categories: objective and subjective.One of the subjective evaluations for evaluating the denoising effectiveness is to observe the resolution of seismic section.The clearer effective reflection on the seismic section means the better denoising effect.Objective evaluation criteria are based on the deviation and mainly consist of the signal-tonoise ratio (SNR) and the mean square error (MSE): where f (i) represents the original signal, f (t) represents the signal reconstructed at a particular scale, and n is the length of the signal.

Simulation experiment
To verify the performance of the algorithm, we used MATLAB to construct a simulation model.The model included the Heavy sine signal, which is widely used in algorithm simulations, a total of 1024 sampling points, and white Gaussian noise with a standard deviation of 5.The sym6 wavelet was used for five levels of decomposition.The wavelet coefficients at each of the five levels were quantized by using the soft, hard, and modified threshold (with adjusting factors of 1, 5, and 10) functions.
Figure 3a shows the original signals, Figure 3b shows the same signals with additive random noise, and Figures 3c-g show the reconstructed signals after denoising.
Figure 3c shows the reconstructed signal after denoising with the soft threshold function.It has relatively fine smoothness as a whole.However, the waveforms at points A and B are too smooth and fail to fully retain the waveforms of the original signal.Figure 3d shows the reconstructed signal after denoising with the hard threshold function.It has glitches at several points and relatively poor smoothness.Figure 3e shows the reconstructed signal after denoising with the modified threshold function and the adjusting factor m = 1, the result is similar to that of the hard threshold function with a serious glitch in the waveform.Figure 3f shows the signal result of the modified threshold function when m = 5.It has relatively fine smoothness and general agreement between the waveforms of the reconstructed and original signals.Figure 3g shows the signal when m = 10.The result is close to that of the soft threshold function and also displays excessive smoothness at points A and B.
These results demonstrate the following benefits of the modified threshold function: (i) it restrains local burrs compared with the hard threshold function and (ii) improves the overall smoothness of the signal compared with the soft threshold function.The optimal adjusting factor appeared to be 5 rather than 1 or 10.
Table 1 indicates that the modified thresholding scheme performed better than the other two schemes in terms of SNR, especially when the adjusting factor was set to 5. The modified threshold function with an adjusting factor of 5 produced the highest SNR values, increased the original signal from 14.2292 to 19.9929, and had the lowest MSE.

Application to synthetic data
In seismic exploration, a reflection seismogram is the convolution of seismic wavelets and the reflection coefficient series.The mathematic model can be expressed as where b(t) is the seismic wavelet, p(t) represents the reflection coefficient series, and n(t) represents noise.
In practice, seismic wavelets are too complex to extract from the data.For this reason, the commonly used Ricker wavelet was introduced as a substitute for seismic wavelet in forward modeling, as expressed in Equation (10) [Kim et al.2013, Zheng et al. 2013]: (10) where f p is the main frequency of the wavelet and t is the propagation time.
Figure 4 shows the waveforms of the Ricker wavelet with main frequencies of 25 Hz.
The reflection coefficient p(t) reflects the redistribution of energy during the reflection and transmission processes and is given by (11) where t 1 and t 2 denote the densities of the upper and lower layers of a stratum, respectively, and o 1 and o  represent the propagation velocities of the seismic wave across the stratum's upper and lower layers, respectively.Assume horizontal strata with a uniform medium and a Ricker wavelet with a dominant frequency of 25 Hz as the simulated seismic signal.For the sake of clarity, only one reflecting layer and a single reflection are considered.The reflection depth is 100 m.The top and bottom velocities are 2000 and 3000 m/s, respectively.The minimum offset distance (distance between the seismic source and first detector) is 10 m.The channel number is 50.The recording duration is 600 s.Then, Figure 5 shows the approximate geometric graph that represents the propagation path of the reflected seismic wave.
Figure 6 shows the denoising result for the artificially synthesized seismic record. Figure 6a shows the seismic convolution model, and Figure 6b shows the seismic record incorporating Gaussian white noise.Figure 6c shows the seismic record of the modified wavelet thresholding when the adjusting factor is 5. Figure 6d shows the differences in the waveforms before and after denoising for the data of a single seismic channel.
Figures 6a and 6b show that the seismic section undergoes severe compromise after noise is added.Thus, the original seismic waveform is extremely hard to distinguish, which adversely affects understanding of the seismic data.Figure 6c shows the processed seismic section; the modified algorithm eliminates most of the random noise and effectively retains the desired reflections.Figure 6d

Application to real data
To check the efficiency of this technique, we tested the algorithm on a raw seismogram recorded during a land seismic exploration experiment located in north China.Figure 7 shows the denoised results of the real seismic records.
Figure 7a shows the group of seismic records extracted from the raw seismogram records.The useful signals were submerged by noise in the original seismic section, especially in the red box area.Figure 7b shows the seismic section after the modified algorithm was utilized for denoising.
Figures 7a and 7b show that the random noise can be effectively eliminated with the modified denoising algorithm; and the effective reflections become clear.This result further proves that the modified algorithm is capable of retaining useful signals while effectively removing random noises.
To offer a clear picture, two typical channels from the actual seismic data were selected to evaluate the amplitude spectra before and after denoising.Figure 8 shows the comparison of amplitude spectra before and after denoising.Figure 8a-b show the amplitude spectra of the two channels data respectively.Figure 8 allows us to know that the wavelet thredholding can reduce the high frequency noises and reserve the useful signal in evidence.

Conclusions
Seismic field data contain information about underground geological structures, but the data often cannot be directly used for geological interpretation because they are susceptible to interference from environmental noise.Thus, the original seismic data must be denoised in order to acquire high-quality seismic data to form a reliable basis for subsequent geological interpretation.Because seismic signal processing is a critical task in seismic exploration, the processing techniques have a direct bearing on progress in the field of seismic exploration.
This study analyzed existing denoising methods used in seismic exploration from the perspective of random noise.The simulations and convolution modeling results provide a basis for the application of wavelet thresholding to reducing the random noise in seismic signals.
(1) A new wavelet-based denoising method was proposed based on an analysis of conventional wavelet thresholding and threshold functions for denoising.Compared to existing denoising methods, this new method fully considers the time-frequency domain and multiscale resolution and is able to remove random noise.
(2) To test the new denoising method, a simulation experiment was conducted on the Heavy sine signal using MATLAB.The results showed that signals denoised by hard thresholding had a higher SNR but lower smoothness than those denoised by soft thresholding.Moreover, the pseudo-Gibbs phenomenon occurred in signals produced by hard thresholding.Compared to these two methods, the modified threshold function performed well in terms of SNR and overall smoothness, especially when the adjusting factor was 5.
(3) A seismic convolution model was built.The modified wavelet thresholding with an adjusting factor of 5 was used to denoise the model with Gaussian white noise.The results showed that the random noise can be suppressed and that the resolution and SNR can be improved.
(4) Denoising was performed with the modified wavelet thresholding on seismic data extracted from a seismic experiment in north China.The results showed that the modified algorithm effectively eliminated the random noise and retained the useful signals, which further proves the feasibility of the modified wavelet thresholding for practical seismic data processing.Seismograms used in Section 5 (Application to real data) were taken from a raw seismogram record in a land seismic exploration experiment located in north China.

Figure 2 .
Figure 2. Schematic of hard and soft threshold functions.(a) Soft threshold function; (b) hard threshold function.

Figure 3 .
Figure 3. Signal denoising process.(a) Original signal; (b) noisy signal containing white Gaussian noise with a standard deviation of 5; (c) the reconstructed signal after denoising with the soft threshold function; (d) the reconstructed signal after denoising with the hard threshold function; (e) the reconstructed signal after denoising with the modified threshold function, m = 1; (f ) the reconstructed signal after denoising with the modified threshold function, m = 5; (g) the reconstructed signal after denoising with the modified threshold function, m = 10.
Figure5shows the approximate geometric graph that represents the propagation path of the reflected seismic wave.Figure6shows the denoising result for the artificially synthesized seismic record.Figure6ashows the seismic convolution model, and Figure6bshows the seismic record incorporating Gaussian white noise.Figure6cshows the seismic record of the modified wavelet thresholding when the adjusting factor is 5. Figure6dshows the differences in the waveforms before and after denoising for the data of a single seismic channel.Figures6a and 6bshow that the seismic section undergoes severe compromise after noise is added.Thus, the original seismic waveform is extremely hard to distinguish, which adversely affects understanding of the seismic data.Figure6cshows the processed seismic section; the modified algorithm eliminates most of the random noise and effectively retains the desired reflections.Figure6dcompares the waveforms of the original signal, noisy signal, and denoised signal of a single seismic channel.The results show that the modified wavelet thresholding eliminates the random noise in the signals.In addition, the waveforms of the original and denoised signals are almost matching at the crest and trough, which confirms that the denoised signal has relatively fine quality.

Figure 4 .
Figure 4. Ricker wavelet with a main frequency of 25 Hz.

Figure 5 .
Figure 5.The propagation path of the reflected seismic wave.

Figure 6 .
Figure 6.Denoising for artificially synthesized seismic record.(a) Seismic convolution model; (b) seismic record incorporating Gaussian white noise; (c) seismic record of the modified wavelet thresholding, m = 5; (d) comparison of the waveforms before and after denoising for the data of a single seismic channel.

Figure 7 .
Figure 7.The denoised results of the real seismic records.(a) Original seismic records; (b) denoised seismic records.

Figure 8 .
Figure 8.The comparison of amplitude spectra before and after denoising.(a) Amplitude spectra of the 131th channel data; (b) amplitude spectra of the 195th channel data.

Table 1 .
Compare the SNRs and MSEs of the different wavelet thresholdings.