Abstract
In seismic exploration, effective seismic signals can be seriously distorted by and interfered with noise, and the performance of traditional seismic denoising approaches can hardly meet the requirements of high-precision seismic exploration. To remarkably enhance signal-to-noise ratios (SNR) and adapt to high-precision seismic exploration, this work exploits the non-subsampled contourlet transform (NSCT) and threshold shrink method to design a new approach for suppressing seismic random noise. NSCT is an excellent multiscale, multidirectional, and shift-invariant image decomposition scheme, which can not only calculate exact contourlet transform coefficients through multiresolution analysis but also give an almost optimized approximation. It has better high-frequency response and stronger ability to describe curves and surfaces. Specifically, we propose to utilize the superior performance NSCT to decomposing the noisy seismic data into various frequency sub-bands and orientation response sub-bands, obtaining fine enough transform high frequencies to effectively achieve the separation of signals and noises. Besides, we use the adaptive Bayesian threshold shrink method instead of traditional handcraft threshold scheme for denoising the high-frequency sub-bands of NSCT coefficients, which pays more attention to the internal characteristics of the signals/data itself and improve the robustness of method, which can work better for preserving richer structure details of effective signals. The proposed method can achieve seismic random noise attenuation while retaining effective signals to the maximum degree. Experimental results reveal that the proposed method is superior to wavelet-based and curvelet-based threshold denoising methods, which increases synthetic seismic data with lower SNR from −8.2293 dB to 8.6838 dB, and 11.8084 dB and 9.1072 dB higher than two classic sparse transform based methods, respectively. Furthermore, we also apply the proposed method to process field data, which achieves satisfactory results.
1. Introduction
In recent years, high-precision seismic exploration has been a key subject in modern seismic exploration. This technique will be hindered if the noise in acquired seismic signals cannot be removed perfectly. The traditional seismic data denoising approaches can hardly meet the requirements of high-precision seismic exploration because the level and complexity of the accompanying noise in seismic signals have significantly increased due to the increasingly complex exploration environment and the increase in exploration depth with the extension of field of seismic exploration. So, it is crucial to design new effective techniques to remarkably enhance the signal-to-noise ratio (SNR).
At present, many seismic denoising approaches have been proposed including the initial seismic data denoising method [1], traditional transform domain based denoising methods [2–4], sparse transform based methods [5–8] for solving multitasks, learning-based methods [9, 10], and other methods [11–14]. Actually, as the most common seismic noise, random noise can penetrate the whole time domain and severely distort and interfere with effective seismic data. Thus, since Canales [1] first developed a random noise reduction approach, a lot of random noise attenuation methods have been presented on this basis, such as the sparse transform based approaches, the empirical mode decomposition (EMD) based approaches, and fast dictionary learning-based approaches [10]. Chen and Ma [4] removed random noise with predictive filtering of f-x empirical mode decomposition. Chen and Fomel [6] developed an EMD-Seislet transform based method to remove the seismic random noise. Liu et al. [7] presented variational mode decomposition for suppressing the random noise of seismic data. In reality, one type of the intensively used and most efficient seismic random noise attenuation approaches are based on the sparse transform of multiscale geometric analysis. Zhang and Lu [2] removed noise and improved the resolution of seismic data by using the applied wavelet transform. Neelamani et al. [3] attenuated random noise with the curvelet transform, and subsequently, several variants [8, 12] with good results have been reported. Lin [14] proposed a three-dimensional (3-D) steerable pyramid decomposition-based suppression method of seismic random noise. Sang et al. [15] presented an unconventional technique on the basis of a proximal classifier with consistency (PCC) in transform domain for attenuating seismic random noise, and they also proposed another seismic denoising approach [16] via the deep neural network and simultaneously suppressed seismic coherent and incoherent noises [17] based on the deep neural network.
High SNR data are the important guarantee of high-precision seismic exploration. But, the existing transform domain based methods are difficult to obtain higher SNR data due to not fine enough transform high frequencies such as wavelet transform or curvelet transform. Compared with the existing sparse transform based methods, that is, wavelet-based transform and curvelet-based transform methods, the NSCT presents multiscale, multidirectional, and shift-invariant decomposition scheme, which has better high-frequency response and stronger ability to describe curves and surfaces. Besides, they often conduct rough threshold operation by using manual threshold processing methods such as hard thresholding or soft thresholding. There is often a loss of effective signals. Therefore, to remarkably enhance SNRs and adapt to high precision seismic exploration, we exploit an effective seismic data denoising method in this paper. The contributions are as follows:(i)We propose to utilize the new sparse transform technique, non-subsampled contourlet transform (NSCT), to decomposing the noisy seismic data into various frequency sub-bands and orientation response sub-bands, obtaining fine enough transform high frequencies to effectively achieve the separation of signals and noises.(ii)We use the adaptive Bayesian threshold shrink method instead of the traditional handcraft threshold scheme for denoising the high-frequency sub-bands of NSCT coefficients, which pays more attention to the internal characteristics of the signals/data itself and improve the robustness of method, which can work better for preserving richer structure details of effective signals.(iii)We conduct the experiments on synthetic and field data, which reveals that our approach is superior to the wavelet and the curvelet transform based classical ones, achieving higher signal-to-noise ratio (SNR) values.
The remainder of this paper is organized as follows. We present our method in Section 2. Experiments and performance evaluation are presented in Sections 3. Conclusion is drawn in Section 4.
2. Method
In this paper, we focus on transform domain based thresholding methods due to their good performance. The wavelet-based thresholding scheme is the most classic method for seismic data denoising. Wavelets can sparsely represent one-dimensional (1-D) digital data with smoothed point discontinuities and have been successfully used for representing digital signals [18]. However, wavelets cannot efficiently handle higher dimensional data because of the usual presence of other kinds of singularities. As a matter of fact, curvelets [19], contourlets [20], bandelets [21], and some other image/signal representations can take the advantages of the anisotropic regularity of a surface along edges, but these representations all have their own disadvantages, such as lack of a multiresolution geometry representation for curvelets, extremely limited clear directional features for contourlets, and computationally expensive geometry optimization for bandelets. The non-subsampled contourlet transform (NSCT) [22] is an excellent multiscale, multidirectional, and shift-invariant image decomposition scheme, which can not only calculate exact contourlet coefficients through multiresolution analysis but also give an almost optimized approximation. It has better high-frequency response and stronger ability to describe curves and surfaces. Therefore, we attempt to utilize NSCT to denoise seismic data in this paper.
2.1. NSCT for Seismic Data
The NSCT [22] primarily consists of a cascade of non-subsampled pyramid filter bank (NSPFB) and non-subsampled direction filter bank (NSDFB). First, the NSPFB is utilized to decompose an image and the sub-bands obtained are used as inputs of the NSDFB to generate decomposition results of the initial image in multiple directions and dimensions. The NSCT conducts K-level decomposition on an image to produce one low-frequency (LF) and several high-frequency (HF) sub-bands, and the size of all these sub-bands is identical with that of the original image. So, the full reconstruction of NSCT is possible since NSPFB and NSDFB can both be completely rebuilt.
The shift-invariant filtering structure of the NSCT results in its multiscale feature. By using a bank of non-subsampled 2-D two-channel filters, we have sub-band decomposition like the Laplacian pyramid. Figure 1 shows a 3-stage non-subsampled pyramid (NSP) decomposition, whose expansion has an alike concept with the 1-D non-subsampled wavelet transform (NSWT) using the àtrous algorithm [22]. For a J-stage decomposition, the redundancy will be J + 1. The region and its complement are the ideal passband support of the low- and high-pass filter at the j-th stage, respectively. Upsampling the filter for the first stage can yield the filters for subsequent stages; thus, no additional filter is needed to give the multiscale property. Our structure differs from that of the separable NSWT. Particularly, our structure produces one bandpass image in each stage leading to J + 1 redundancy. However, the NSWT generates three directional sub-bands in each stage, leading to 3J + 1 redundancy. The advantage of NSP lies in its ability to generate better filters due to its generality.
The NSCT has following steps. First, the NSPFB is used to multidimensionally decompose an image into an HF and an LF sub-band. In multilevel decomposition, an image is finally decomposed to an LF sub-band and a set of HF sub-bands if only its LF sub-band is further iteratively filtered. The redundancy of an X-level NSPFB decomposition will be X + 1. With regard to the X-level low-pass and the bandpass filter, the ranges of the ideal support of frequency domain are and , respectively. The decomposition does not extra filter during acquiring multidimension properties. So, the generated redundancy of X + 1 will be generated in each stage with a bandpass image, and the structure is significantly superior to that of the wavelet transform. Then, these sub-bands can be decomposed along singular points and multiple directions in various dimensions, and the directions are integrated. The NSDFB also belongs to a two-channel filter bank which comprises decomposition filters and synthesis filters satisfying Bézout’s identity:
To adopt ideal support of the frequency domain, two channels are decomposed by and . Then, and are upsampled instead of subsampled by all sampling matrices in each level to get direction filters in the subsequent levels. This completes the image decomposition and an NSCT transform is schematically presented in Figures 2 and 3. Figure 2 displays an overview of NSCT which has a filter bank for dividing the 2-D frequency plane into sub-bands plotted in the bottom left quarter of Figure 1. By using NSCT, we decompose the 2-D seismic signal data into two shift-invariant components: an NSP structure to ensure multiscale properties (Part 1 of Figure 2) and an NSDFB structure to give directionality (Part 2 of Figure 2). The obtained idealized frequency partitioning diagram is presented in Figure 3. The structure consists in a bank of filters that splits the 2-D frequency plane into several sub-bands. In this paper, we use this mode of non-downsampling to reduce the sampling distortion in the filters and obtain translation invariance, in which the size of the directional sub-band at each scale is the same as that of the original 2-D seismic signal matrix. The NSCT has more details to be preserved, and the decomposition can better maintain the edge information and contour structure of the seismic signals.
Figure 4 shows processing results of implementing the two-level NSCT on the synthesized seismic signal data with noise (Figure 4(a)) to yield a low-pass sub-band (Figure 4(b)) and a set of high-pass sub-bands (Figures 4(c) and 4(d)). Here, two and four shearing directions are used for the coarser and the finer scale, respectively. We can see from Figure 4 that the LF record (Figure 4(b)) decomposed by the NSCT basically contains the effective synthesized seismic signals. For HF records (Figure 4(d)) of scale 1, all they contain is noise, while HF records (Figure 4(c)) of scale 2 contain partially effective signals and noise, which needs to be further processed via signal-noise separation.
(a)
(b)
(c)
(d)
2.2. Denoising Seismic Data Using the Threshold Shrink in the NSCT Domain
To remarkably suppress seismic random noise while not damaging the effective signal in the process of denoising, this paper proposes a novel NSCT-based scheme with an adaptive threshold value setting for suppressing seismic random noise. The NSCT with the properties of multiscale, multidirection, and relative optimized sparsity can not only calculate the exact contourlet coefficients through multiresolution analysis but also give an almost optimized approximation. As the NSCT is multidirectional, large coefficients can be obtained when the direction of the NSCT basic function is approximately the same as the direction of the seismic signals, while small coefficients can be obtained when they have large difference. Thus, the random noises are distributed on small coefficients, so we can remove smaller coefficients to achieve random noise attenuate using an appropriate thresholding operator.
We first analyze the sparsity of NSCT before giving the steps of denoising. It is known that the degree of approximation of the decomposed effective data determines the effect of noise suppression [23]. That is to say, the denoising effect depends on the sparsity of the approach. Figure 5 presents the reconstruction error on synthesized data (Figure 6(a)) in the wavelet transform domain, curvelet transform domain, and NSCT domain. Clearly, the construction error of NSCT is the smallest at the same percentage of coefficients, and it is approximate to zero at 6% coefficient, showing its optimal sparsity. In Figure 6, we compare the high-frequency coefficients of NSCT, curvelet transform, and wavelet transform, where we can clearly see that the NSCT represents the curvature more accurately.
(a)
(b)
(c)
(d)
Generally, one threshold is used for the whole image/signals (or sub-bands) in signal denoising techniques based on the threshold shrink. Obviously, the threshold value should be smaller if the signals contain more effective information, and it should be larger if the signals have more smooth regions. For the two cases, a larger threshold value should be correspondingly used for a higher noise level. Evidently, detailed information with an optimal threshold value does not function adequately for smooth regions and vice versa. Therefore, setting for threshold value can be further optimized by introducing adaptive threshold for different regions in seismic signals to exploit the fact that most signals consist of smooth regions and effective seismic signal information.
Specifically, the two-dimensional noisy seismic data can be calculated bywhere denotes time, represents the trace number, and , , and represent the effective seismic data, the additive random noise, and the noisy observed seismic signals, respectively. Signal will be recovered from .
In NSCT-based denoising approaches, a threshold is properly set for the NSCT coefficients so that seismic signals can be retrieved from the acquired noisy seismic data. The proposed seismic random noise attenuation has the following main steps.
Step 1. Decomposing noisy seismic data with a K-level NSCT to yield one low-pass sub-band and one set of high-pass sub-bands , with the current scale k, the decomposition orientation j, and the total number J of decomposition directions.
Step 2. Calculating denoising threshold values of all sub-bands . The level adaptive Bayesian threshold [24] is used and calculated as below:(i)Using the robust median estimator to calculate noise variance from sub-bands:(ii)Using the maximum likelihood estimator (MLE) [24] to estimate signal variance for the noisy coefficients of each detail sub-band : where , and m and n denote the size of seismic signals(iii)Calculating discriminating threshold with the near exponential prior of NSCT coefficients across scales: where k denotes the current scale(iv)Calculating denoising threshold of each sub-band for :where denotes the standard deviation of sub-band
Step 3. Processing the noise-related NSCT coefficients in high-frequency sub-bands with the well-known soft-thresholding method [25]:where and are the matrices of the coefficients after and before denoising in the NSCT domain, respectively. denotes the adaptive Bayesian threshold
Step 4. Reconstructing the denoised seismic data by conducting an inverse NSCT on these denoised NSCT sub-bands.
The workflow of our method is figuratively presented in Figure 7, where a two-level NSCT is performed on synthetic data to yield a low-pass sub-band and a set of high-pass sub-bands, and two and four shearing directions are used for coarser and finer scale, respectively.
3. Experiments
In this section, we evaluate our proposed method with two classic sparse transform based methods (wavelet- and the curvelet-based thresholding schemes) on synthetic seismic data and field data. The hardware used in this experiment is Intel 6226R CPU @2.90 GHz processor, 93 GB memory and NVIDIA RTX 3090 24G graphics card. The software used is Matlab R2016b.
3.1. Synthetic Seismic Example
To demonstrate the performance, an example of hyperbolic-events synthetic data is used. Figure 8(a) presents the synthesized signals of 150 traces with 1 ms time sampling interval. The Ricker wavelet is expressed by
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
Figure 8(b) is the corresponding noisy synthetic data, which is denoised by using two existing methods, namely, the wavelet- and the curvelet-based threshold denoising approach and the proposed approach. These three methods use the same threshold method. We set up K = 2 and J = 2 for the proposed approach. So, we can acquire the denoising threshold by computing formulas (3-6) step by step; thus, we can further obtain the denoised high-frequency results by the soft-thresholding method. Figures 8(c)–8(e) show the obtained results by three methods, respectively. Obviously, the result of our approach is much superior to the ones by other two approaches. Concretely, the results are evaluated by using SNRs [26]:where and represent the noise-free data and the noisy or denoised data, respectively. The resulted SNR values for Figures 8(b)–8(e)) are −2.2063 dB, 5.8541 dB, 8.2496 dB, and 13.2078 dB, respectively. Obviously, wavelet- and the curvelet-based approaches implement insufficient noise removal, while our approach conducts good performance in attenuating most random noise and the SNR value has been significantly improved. Figures 8(f)–8(h) show the removed noise sections by these three approaches, respectively. The wavelet- and the curvelet-based approaches lose part of useful signals (red arrow). Obviously, the proposed approach does not harm any useful signals. Besides, Tables 1 and 2 present the summary of the results with various SNRs before and after denoising. Our approach shows the better denoising performance, especially for the low SNR seismic data.
In addition, to validate the processing result of our method, real noisy seismic data are measured with same excitation and reception in the identical data area. Figure 9 presents the acquired noisy signals (Figure 9(a)) and the denoising result with the proposed approach (Figure 9(b)). We can see that several highlighted effective signals, clearer interlayer structure, and improved event continuity can be observed from the patterns of the denoised data, which significantly improves the SNR value. Figure 10(a) shows the real stacked profile. Similarly, after processing using our proposed approach, effective signals are highlighted; information between layers is richer and noises are effectively suppressed, significantly improving the SNR (Figure 10(b)); it can be seen from the removed noise that there is basically no effective signal loss.
(a)
(b)
(a)
(b)
(c)
4. Conclusion
This article presents a novel NSCT-based seismic random noise denoising method. The superior performance NSCT with an appropriate thresholding operator brings excellent denoising results for seismic signals. The proposed method can achieve seismic random noise attenuation while retaining effective signals to the maximum degree. The experiments are performed with both synthesized and real seismic signals and the results demonstrate effectiveness of our approach compared with existing ones. In the future, we will consider deep learning based techniques to denoise seismic data with low SNR in view of the powerful learning ability and feature recognition ability, which aim to highlight effective signals and suppress false signals.
Data Availability
The data included in this paper are available without any restriction.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.