Introduction

Seismic data obtained from complex geological environments may suffer from severe noise contamination (Simon and Jianwei 2013). Compared with data collected from simple geological conditions, both the resolution and signal-to-noise ratio (SNR) of the seismic record are greatly reduced.

In terms of seismic noise suppression, there are several classic methods that can be applied to field seismic data processing, such as adaptive filtering (Han et al. 2015), f-k filtering (Mostafa 2012), polynomial fitting (Sanyi and Shangxu 2013), independent component analysis (Chuntian and Wei 2010), time frequency peak filtering (Chao and Yue 2015), f-x domain filtering (Yangkang and Jitao 2014) and τ-p domain filtering (Famoush et al. 2012). F-k filtering is efficacious at eliminating surface wave interference, whereas τ-p domain filtering mainly focuses on linear interference. However, these mature technologies cannot operate effectively in the presence of strong and complex random noise, even results in false seismic events (Chao and Yue 2015).

In recent years, multi-resolution analysis in the transform domain has opened up new methods for seismic data denoising (Shucong et al. 2014; Mei et al. 2012; Andrzej et al. 2014; Chengming et al. 2014). The multi-scale wavelet transform has been widely used in seismic processing (Rongfeng et al. 2003; Shucong et al. 2014). However, with regard to reflection events, the wavelet transform encounters a directional limitation and shift variance. These shortcomings result in serious distortion, especially around peaks and troughs in the signal (Chengming et al. 2014; Ivan et al. 2005). More recently, some new methods based on sparse representation have been developed, including Ridgelets (Mei et al. 2012), Contourlets (Yang et al. 2009), Curvlets (Andrzej et al. 2014) and Shearlets (Chengming et al. 2014). These offer a promising framework for seismic data denoising, although the non-subsampled property of the Curvlets and Shearlets algorithms increases their time complexity.

To solve the problems described above, this paper presents an improved method based on the dual-tree complex wavelet transform (DTCWT). Our approach is based on variations in different directions and scales between the seismic signal and the noise. An adaptive threshold is applied in low-frequency scales where the signal dominates the noisy record. In high-frequency scales, where the noise and the signal are of similar strength, a block-matching (BM) algorithm and singular value decomposition are applied. This approach takes temporal continuity and spatial correlation of reflected events into account.

The paper is organized as follows. The fundamentals of the DTCWT are summarized in Section “Fundamentals of the dual-tree complex wavelet transform,” and the proposed method to mitigate noise is described in Section “Combined noise suppression method in the DTCWT domain.” In Section “Simulations and application,” the performance of the proposed scheme is evaluated using simulated seismic signals and measurement data. Section “Conclusion” presents our conclusions.

Fundamentals of the dual-tree complex wavelet transform

Unlike the classical discrete wavelet transform (DWT), the DTCWT is constructed via a complex-valued scaling function and complex-valued wavelet function. The complex-valued wavelet function ψ c (t) is constructed as follows; the complex-valued scaling function is defined similarly (Hongye et al. 2014; Ivan et al. 2005):

$${\psi_c}\left( t \right) = {\psi_r}\left( t \right) + j{\psi_i}\left( t \right)$$
(1)

where ψ r (t) is even and real, jψ i (t) is odd and imaginary, but ψ i (t) is real. In addition, ψ r (t) and ψ i (t) form a Hilbert transform pair.

Consider the 2-D DTCWT associated with the row–column implementation of the 1D-DTCWT:

$$\psi \left( {x,y} \right) = \psi \left( x \right)\psi \left( y \right)$$
(2)

where ψ(x) and ψ(y) are given by formula (1). To represent an integrated real two-dimensional signal completely, the row or column of the complex conjugate filter is required. Three sub-bands are produced in both the first and second quadrants, corresponding to six directions in space: ±15°, ±45°, ±75°.

Compared with 2-D DWT, 2-D DTCWT has a number of excellent properties for seismic signal processing, such as shift invariance, adequate directions, anti-aliasing and limited data redundancy (Ivan et al. 2005). In particular, because of the shift invariance property, the total energy at scale j is nearly constant. In contrast to the DWT, the coefficients do not oscillate between positive and negative values around singularities. That is, the signal obtains a sufficiently sparse representation in the DTCW domain.

Combined noise suppression method in the DTCWT domain

Noise suppression can also be achieved using a fixed hard threshold. However, the signal and the noise have different features in different sub-bands. Clearly, a fixed threshold cannot take this into account (Hongye et al. 2014). Suppose that the noisy seismic data model is:

$$Y = S + N$$
(3)

where Y indicates the noisy record, S is the pure seismic data, and N represents the noise model. To recover the pure data S from Y as effectively as possible, we propose a method that combines an adaptive threshold, BM and singular value decomposition (SVD) in the DTCWT domain.

Denoising at low-frequency scales

The frequency of seismic reflected events is generally <50 Hz, whereas that of random noise is much higher. As a result, there is a large amount of signal and only a small amount of background noise at low-frequency scales, as shown in Fig. 1a. This makes it relatively easy to suppress the noise.

Fig. 1
figure 1

Two sub-bands. a Low-frequency sub-band. b High-frequency sub-band

Full size image

To prevent the signal loss caused by the fixed threshold, we propose an adaptive method. The adaptive threshold T is chosen as follows:

$$T = c \cdot \sigma$$
(4)

where \(\sigma\) is the standard deviation of each sub-band. Anisotropy causes the content of the useful signal to vary with different directions of the scale j. The energy of each sub-band is calculated by:

$${E_k} = \mathop \sum \limits_r \mathop \sum \limits_c {\left( {{d_{r,c}}} \right)^2}\;,\quad k = 1,2 \ldots 6,\;$$
(5)

where d r,c is the coefficient on column c and row r.

The direction of a reflection event is simple and symmetrical. This simplicity means that the signal energy is usually concentrated in two direction sub-bands. The energy in these two sub-bands is equal and greater than that in the other four sub-bands because it contains energy from both the signal and the noise. For the two sub-bands with maximal energy, c is set to 1.5; c is set to 2.5–3 for the other sub-bands in the case of synthetic data.

Denoising at high-frequency scales

However, at high-frequency scales, the situation is not conducive to noise suppression. A small amount of signal and most of the noise are mixed together with trivial differences, as shown in Fig. 1b. This makes it difficult to remove the noise through ordinary threshold processing. We propose a novel noise suppression method for high-frequency sub-bands that combines BM and SVD. The noise is attenuated using low-rank approximation based on signal characteristics that can be distinguished from the noise.

Block matching

The seismic events exhibit temporal continuity in each channel and spatial correlation between adjacent channels. The additive random noise does not contain these characteristics (Zhipeng et al. 2012). Therefore, the waveforms on each channel are similar. Using this similarity, we can roughly distinguish between the signal and the noise. To identify the similarity between reflection events, BM is employed as a crucial filtering step.

The matching of similar blocks is illustrated in Fig. 2. The original data are divided into blocks Y x of size N 1 × N 1. The blocks overlap by N s . The search window is set to N h  × N h . Select a reference block Y R and calculate the Euclidean distance between Y R and candidate blocks Y x in the search window as:

$$d\left( {{Y_R},{Y_x}} \right) = \frac{{{|| r^\prime }\left( {Y_R} \right) - {r^\prime }\left( {Y_x} \right)||_2^2}}{{{{\left( {N_1} \right)}^2}}}$$
(6)

where r′ is a pre-threshold processing to prevent the erroneous grouping caused by large noise variance.

Fig. 2
figure 2

Schematic of matching

Full size image

When the distance is less than a given threshold τ, two blocks are judged to be similar. All similar sub-blocks are stored in a 3-D array Z R , along with the current reference block. Obviously, if the reference block contains the useful signal, the blocks in Z R should all contain reflected event.

Local- and non-local denoising

For seismic signals, each block containing the pure reflected event gives a low-rank matrix because of its similarity. When the signal is contaminated by noise, the matrix becomes full rank. To achieve the low-rank equivalent, each block in Z R is subjected to SVD. For any sub-block Y x , there must exist two orthogonal matrices \({U_{1D}} \in {R^{{N_1} \times {N_1}}}\) and \({V_{1D}} \in {R^{{N_1} \times {N_1}}}\) that satisfy the following formula:

$${Y_x} = {U_{1D}}{\varLambda_{1D}}V_{1D}^T$$
(7)

where Λ 1D is an N 1 × N 1 diagonal matrix whose diagonal elements are the singular values of Y x in non-increasing order. The larger singular values mainly represent the signal, and the smaller singular values mainly represent the noise. This view is demonstrated in Fig. 3. We retain the larger singular values and replace the smaller singular values with zeros. The matrix can then be restructured as low-rank estimation. In this step, the local similarity has been taken into consideration to suppress the noise in each block.

Fig. 3
figure 3

Singular value of synthetic seismic signal blocks. The vertical coordinate is the singular value, and the horizontal coordinate is the index of the singular value

Full size image

The BM process results in the blocks within a group being similar. To further suppress noise, the discrete cosine transform is applied for the corresponding position points in different sub-blocks of the group. The transform coefficients are then subjected to a hard threshold. In this step, the non-local similarity is considered to suppress the noise in neighboring blocks.

Aggregation

The inverse discrete cosine transform is applied to obtain a basic estimation \({\hat Y_x}\) of all blocks in the grouping. Because the blocks are overlapped with each other, the final estimation \({\hat Y^{\text{est}}}\) is reconstructed by a weighted summation. The weights ω ht R are calculated as follows:

$$\omega_R^{ht} = \frac{1}{{{\sigma^2} \times {N_R}}}$$
(8)

where N R denotes the number of nonzero coefficients and σ 2 represents the noise variance. Overlapping blocks are aggregated in the following manner:

$${\hat Y^{\text{est}}} = \frac{{\mathop \sum \nolimits_{Y_R} \mathop \sum \nolimits_{Z_R} \omega_R^{ht} \times {{\hat Y}_x}}}{{\mathop \sum \nolimits_{Y_R} \mathop \sum \nolimits_{Z_R} \omega_R^{ht}}}$$
(9)

The results of processing the high-frequency scale using different methods are shown in Fig. 4. Figure 4a shows the unprocessed high-frequency sub-band containing events submerged in strong noise. Figure 4b shows the denoising results using a fixed threshold. Although the events are effectively retained, an extensive amount of noise also remains. Compared with thresholding, BM and SVD obtain clearer events and less noise, as shown in Fig. 4c.

Fig. 4
figure 4

Denoising results at the high-frequency scale. a Noisy coefficients, b result using fixed hard threshold, c result using BM and SVD

Full size image

Complexity

The time complexity of the algorithm depends on the size of the input seismic record and the parameters described above. Assuming that the size of the noisy sub-band is N × N, the time complexity of the adaptive threshold is N 2. For the BM and SVD processes, assume a block size of 1 × 1. The step between two blocks is N s , and the search window is N h  × N h . The time complexity is

$$\frac{{{N^2} \times N_h^2}}{N_S^2} + \mathcal{O}\left( {N^3} \right)$$
(10)

where the first term is the complexity of BM and the second term is the complexity of SVD.

It is clear that BM and SVD have much greater complexities than the threshold. Therefore, considering the denoising result of high-frequency scales, BM and SVD are applied instead of the adaptive threshold at fine scales. To ensure a reasonable runtime, BM and SVD are not used at coarse scales.

Simulations and application

Experimental results and analysis

To verify the performance of the proposed algorithm, we tested our method using the 100-channel synthetic seismic record. For clarity, only 50 traces are plotted as shown in Fig. 5a. The record contains three reflected events with dominant frequencies of 45, 35 and 25 Hz. The sample frequency is 1000 Hz. The apparent velocities of these events are 1000, 1400 and 1500 m/s, respectively. And the geophone interval is 10 m. White Gaussian noise of −5 dB is embedded. The noisy signal is shown in Fig. 5b. In Fig. 5c, d, we can clearly see that the method proposed in this paper retains reflected events well and suppresses noise effectively, compared with the fixed threshold in DTCWT domain. The SNR increased to 13.88 dB instead of the 11.27 dB achieved with a fixed threshold.

Fig. 5
figure 5

Random noise attenuation result (synthetic seismic signal). a Synthetic signal. b Noisy synthetic signal with SNR = −5 dB, c result based on threshold, SNR = 11.27 dB. d Result based on our method, SNR = 13.88 dB

Full size image

The waveform and spectrum of a single synthetic seismic channel (the 30th) are plotted in Figs. 6 and 7, respectively. The results in Fig. 6b demonstrate that our method retains the signal amplitude satisfactorily while effectively suppressing the noise. In particular, the noise is attenuated from 200 to 700 Hz, as shown in Fig. 7b.

Fig. 6
figure 6

Part of noise attenuation result of synthetic seismic signal

Full size image
Fig. 7
figure 7

Amplitude spectrum of part of synthetic seismic signal

Full size image

Real seismic record processing

To test the practicality of the proposed method, we examined field seismic data acquired from a certain zone of China. These field data consist of common-shot-point records with 168 traces and record length of 6 s. The sampling frequency is 1000 Hz, and the geophone interval is 30 m. In this seismic record, the reflection events are buried in strong random noise. As shown in Fig. 8a, there are some areas in which the reflection events are too disordered to be identified. Moreover, there is some other type of noise, which may be owing to the collection device.

Fig. 8
figure 8

Noise attenuation result of real-field data. a Real-field data. b The result of random noise attenuation method based on threshold. c The result of our method

Full size image

The processing results of thresholding in the DTCWT domain and the proposed method are shown in Fig. 8b, c, respectively. For a detailed comparison, the denoising results of the two methods are shown in Fig. 9b, c. In Figs. 8b and 9a, we can clearly see that the random noise is suppressed effectively, but some residual noise remains. By comparison, Figs. 8c and 9b show that our method has suppressed the noise more effectively, especially in the marked areas. Moreover, the events appear more continuous and smoother.

Fig. 9
figure 9

Part of noise attenuation result of real-field data. a The result of random noise attenuation method based on threshold. b The result of our method

Full size image

Conclusion

In this paper, we have proposed an improved DTCWT denoising method for seismic random noise attenuation. The difference in various directions and frequencies between seismic signals and random noise in the DTCWT domain is the foundation of denoising. According to the similarity of seismic signals, BM was introduced to handle high-frequency scales. Despite of improved noise suppression, the proposed method raises the computation time. The local similarity in each sub-block allows the SVD to be obtained from a low-rank equivalent to achieve noise suppression. And then, noise is further reduced based on non-local similarity. Experimental results using synthetic records and field data prove that the proposed algorithm recovers events with better continuity and suppresses noise more effectively than threshold processing in the DTCWT domain.

However, using a single trace record, SVD leads to signal loss. Therefore, further research is needed to improve this issue.