Introduction

Reflected seismic signals describing the underlying geological structure are usually contaminated by random noise. Irregular interference noise could be brought into the recorded data during the field data collection. In addition, some processing procedures, such as velocity analysis, deconvolution and migration, also lead to erratic noise (Yilmaz 1987). Noise interference may yield unrealistic artifacts in inversion or imaging results, hindering the extraction of geological information. Therefore, denoising is an important processing step. However, it is challenging to effectively remove noise from noisy seismic data.

Various denoising methods have been developed and applied to eliminate seismic random noise, such as median filtering (Bednar 1983), predictive filtering (Canales 1984; Gulunay 1986; Spitz and Deschizeaux 1994; Abma and Claerbout 1995; Guo et al. 1995; Kang et al. 2003), K-L transform (Hemon and Mace 1978; Jones and Levy 1987; Al-Yahya 1991; Liu 1999; Kasina 2010), SVD (Freire and Ulrych 1988; Bekara and van der Baan 2007), wavelet transform (Zhang and Ulrych 2003; Fantine et al. 2019), shearlet transform (Zhang and van der Baan 2019), spectra analysis (Ebrahim et al. 2018) and time–frequency peak filtering (Li et al. 2018), etc. Among them, applying the coherence of multitrace seismic signals and the randomness of noise is a promising scheme. Both K-L transform and SVD belong to this category. The K-L transform was first introduced to seismic data processing for enhancing the signal-to-noise ratio of seismic records by Hemon and Mace (1978). The K-L filtering can be used to remove noise effectively, but it has no fast algorithm. To overcome expensive computational cost, Liu (2016) applied a self-organization learning algorithm to perform orthogonal decomposition of input seismic data. According to the actual needs of filtering, only a few eigenimages are extracted, thus reducing computation costs. Freire and Ulrych (1988) first applied SVD technique in seismic data processing. They derived the relationship between the K-L reconstruction approach and SVD and separated VSP wave fields from the viewpoint of eigenimage filtering. Bekara and van der Baan (2007) made use of SVD technique to suppress noise in seismic data, enhancing signal-to-noise ratio. Existing SVD techniques work in time domain.

In the research presented in this paper, we propose a novel denoising technique which works by SVD filtering in frequency domain. The frequency spectrum of random interference noise is not localized and shows poor coherence. But in the same area, the dominant frequency and bandwidth of effective signals in seismic traces are nearly the same, i.e., the frequency spectra of effective signals show a high degree of trace-to-trace correlation. Based on this point of view, firstly, the input seismic data are transformed to the frequency domain via the Fourier transform. Then, the frequency spectra are decomposed into eigenimages by means of SVD. In terms of selected singular values corresponding to the target signal, eigenimages are extracted to reconstruct the filtered spectra. The last step is to take the inverse Fourier transform, back to the time domain. The outline of this paper is as follows. First, we present the procedure of SVD for noise attenuation in the frequency domain. Then, we show results for synthetic and real data.

Methodology

Singular value decomposition (SVD)

Let \( {\mathbf{X}} \) represent a seismic data matrix which contains \( m \) traces each with \( n \) data points, i.e.,

$$ {\mathbf{X}} = (x_{ij} )_{m \times n} ,\quad i = 1,2, \ldots ,m;\;j = 1,2, \ldots ,n $$

generally \( m < n \). The SVD of \( {\mathbf{X}} \) is given by Freire and Ulrych (1988)

$$ {\mathbf{X}} = {\mathbf{USV}}^{T} = \sum\limits_{k = 1}^{r} {\sigma_{k} {\mathbf{u}}_{k} {\mathbf{v}}_{k}^{T} } , $$
(1)

where \( r = rank ({\mathbf{X}} ) \), \( {\mathbf{u}}_{k} \) is the kth eigenvector of the covariance matrix \( {\mathbf{XX}}^{T} \), also known as the propagation vectors, \( {\mathbf{v}}_{k} \) is the kth eigenvector of \( {\mathbf{X}}^{T} {\mathbf{X}} \), also known as the eigenwavelets, and \( \sigma_{k} \) is the kth singular value of \( {\mathbf{X}} \). The matrices \( {\mathbf{U}} = ({\mathbf{u}}_{1} ,{\mathbf{u}}_{2} ,\ldots ,{\mathbf{u}}_{r} ) \) and \( {\mathbf{V}} = ({\mathbf{v}}_{1} ,{\mathbf{v}}_{2} ,\ldots ,{\mathbf{v}}_{r} ) \) are both orthogonal. \( {\mathbf{S}} = {\text{diag(}}\sigma_{1} ,\sigma_{2} , \ldots ,\sigma_{r} ) \) is a diagonal matrix in which the singular values are ordered in decreasing magnitude as \( \sigma_{1} \ge \sigma_{2} \ge \cdots \ge \sigma_{r} \). According the matrix theory, \( \sigma_{k} \) can be shown to be the positive square roots of the eigenvalues of \( {\mathbf{XX}}^{T} \). The term \( {\mathbf{u}}_{k} {\mathbf{v}}_{k}^{T} \) is an (\( m \times n \)) matrix, which is referred to as the kth eigenimage of \( {\mathbf{X}} \). Eigenimages form an orthogonal basis for the representation of \( {\mathbf{X}} \).

The SVD filtering depends on the coherence among seismic traces. The SVD-filtered data \( {\tilde{\mathbf{X}}} \) is obtained by taking only the contribution of the partial eigenimages into account. This can perform by selection of singular values. For instance, the enhanced signal is

$$ {\tilde{\mathbf{X}}} = \sum\limits_{k = 1}^{p} {\sigma_{k} {\mathbf{u}}_{k} {\mathbf{v}}_{k}^{T} } , $$
(2)

and removed noise is given by

$$ {\tilde{\mathbf{N}}} = \sum\limits_{k = p + 1}^{r} {\sigma_{k} {\mathbf{u}}_{k} {\mathbf{v}}_{k}^{T} } . $$
(3)

SVD denoising in frequency domain

The basis of this work lies in: the frequency spectrum of random interference noise has poor coherence, but in the same area, the dominant frequency band of effective signals shows a high degree of trace-to-trace correlation.

The steps of SVD denoising technique in frequency domain are as follows:

  1. (1)

    Input seismic data \( {\mathbf{X}} = (x_{ij} )_{m \times n} \), \( i = 1,2, \ldots ,m \); \( j = 1,2, \ldots ,n \).

  2. (2)

    Take the split radix FFT of each trace.

    \( {\hat{\mathbf{X}}} = {\text{FFT[}}{\mathbf{X}} ]= {\hat{\mathbf{X}}}_{R} + j{\hat{\mathbf{X}}}_{I} . \)

  3. (3)

    Perform the SVD transform for \( {\hat{\mathbf{X}}}_{R} \) and \( {\hat{\mathbf{X}}}_{I} \), respectively.

    Owing to the conjugate symmetry of discrete Fourier transform (Lyons 2001), the SVD transforms of both \( {\hat{\mathbf{X}}}_{R} \) and \( {\hat{\mathbf{X}}}_{I} \) are done only for the first half length of frequency sample points (within folding frequency).

  4. (4)

    Analyze the distribution of singular values associated with \( {\hat{\mathbf{X}}}_{R} \) and \( {\hat{\mathbf{X}}}_{I} \), determining filter parameter \( p \).

  5. (5)

    Perform the filtering.

    Compute \( {\mathbf{\tilde{\hat{X}}}}_{R} \) and \( {\mathbf{\tilde{\hat{X}}}}_{I} \) by Eq. 2.

  6. (6)

    Obtain \( {\mathbf{\tilde{\hat{X}}}} = {\mathbf{\tilde{\hat{X}}}}_{R} + j{\mathbf{\tilde{\hat{X}}}}_{I} \).

  7. (7)

    Take the split radix IFFT of each trace.

    \( {\tilde{\mathbf{X}}} = {\text{IFFT[}}{\mathbf{\tilde{\hat{X}}}} ]. \)

The conventional filtering method in frequency domain designs a band-pass filter in terms of the frequency range of effective signals, enhancing signals and rejecting noise. Its drawback is that it fails to remove the noise component within frequency pass band and preserve effective signals in low and high frequency. The presented SVD method in frequency domain is different from the conventional frequency filtering. Instead of defining a frequency filter, we design a SVD filter in eigensubspace. It is capable of attenuating the noise component within frequency pass band and protects effective signals in low and high frequency.

Experiment with synthetic data

Figure 1 represents a noisy synthetic section that consists of 48 traces with 250 samples per trace. The sampling interval of the data is 2 ms. The Nyquist frequency is 250 Hz. The wavelet on each trace is a 40 Hz Ricker wavelet with zero phase. The section has been corrupted with additive Gaussian noise (mean μ = 0, variance σ = 30.0). There are three horizontal reflectors, the lowest one of which contains a fault. By means of split radix FFT, the dataset is transformed into the frequency domain. Figure 2a, b shows the real and imaginary components of the frequency spectra, respectively. SVD is successively applied to Fig. 2a, b. The distribution curves of singular values corresponding to real and imaginary parts are separately plotted in Fig. 2c, d. It can be seen that there is little difference between distribution trends of the two curves in Fig. 2c, d. Therefore, in the presented SVD denoising scheme on frequency spectra, the determination of filter parameter \( p \) needs only the singular values curve of real part (or imaginary, either of them). Figure 2e–g shows the first three eigenimages of the real component in Fig. 2a. We take \( p = 3 \) to perform SVD filtering for Fig. 2a, b and then go back to time domain via inverse FFT. The denoised section is displayed in Fig. 3.

Fig. 1
figure 1

Synthetic seismic data with additive Gaussian noise, the lowest reflector shows a fault

Full size image
Fig. 2
figure 2

Singular values decomposition of the frequency spectra. a The real part of the frequency spectra of the seismic data in Fig. 1. b The imaginary part of the frequency spectra of the seismic data in Fig. 1c Singular values of the real part, with the blue line representing whole singular values and the red line representing the selected singular values to be used in reconstruction. d Singular values of the imaginary part, with the blue line representing whole singular values and the green line representing the selected singular values to be used in reconstruction. e The first eigenimage of the real part. f The second eigenimage of the real part. g The third eigenimage of the real part

Full size image
Fig. 3
figure 3

Denoised result. The first three singular values in Fig. 2c, d are selected for reconstruction

Full size image

Real data examples

Two real data examples considered here are taken from land surveys in a basin in northwestern China using a dynamite source. The ground weathered layer is thick and heterogeneous. As a result, explosion and receiving conditions are less good. The recorded data contain strong noise. The stacked sections are often contaminated by stacking noise, largely due to static correction, velocity analysis and deconvolution.

Figure 4 shows a stacked section consisting of 130 traces with 1000 samples per trace. The existence of random noise reduces the quality of seismic section. The deep reflection is concealed by the strong background noise. To understand the characteristics of the seismic data, we perform the frequency spectrum analysis. The magnitude spectrum of the sixth seismic trace in the section is given in Fig. 5. It can be seen that the noises occupy a broad frequency band and show a strong energy level. The dominant frequency of effective signal is about in the range of (10, 70) Hz. The conventional frequency filtering is less capable of attenuating the noise within the range of frequency pass band. Meanwhile, the effective information within the rejecting range may be cut off. Instead of the frequency filter, we apply a SVD filter to the real and imaginary components of the frequency spectra, which is the Fourier transform of input seismic data. The noise reduction is performed by selecting singular values used in the reconstruction of eigenimages. Figure 6 shows the distribution curve of singular values corresponding to the frequency spectra of Fig. 4, plotted in blue line. Figures 7, 8, 9 and 10 are, respectively, the results of reconstruction associated with singular values in Fig. 6a–d. In Fig. 7, only the contribution of the first ten eigenimages is taken into account. As a result, the random noise is eliminated completely, but a little loss of local horizon details occurred (e.g., the area marked by the ellipse). Figure 8, the reconstruction of Fig. 6b (from \( \sigma_{11} \) to \( \sigma_{15} \)), shows missed details. Finally, we make use of the first 15 singular values (Fig. 6c) to extract the effective signals. The enhanced signals are shown in Fig. 9. Figure 10 shows the removed noises (the reconstruction of Fig. 6d with Eq. 3). Clearly, there is no loss of geological information. Applying a band-pass filter with a parameter of (0, 10, 70, 80) Hz to Fig. 4, the denoised data are displayed in Fig. 11. The background noise remains strong, and the continuity of events is less good. Figure 12 is the difference plot between Figs. 4 and 11. The magnitude spectrum of the sixth seismic trace in Fig. 9 is also plotted in Fig. 5. Figure 13a, b shows, respectively, the frequency spectra of the raw data (Fig. 4) and the denoised section (Fig. 9). Inspection of Figs. 5 and 13 shows the presented method is capable of attenuating the noise within the range of frequency pass band. In addition, the denoised section has a broad frequency band, which is beneficial to preventing decrease in the resolution of seismic data (Yilmaz 1987; Fons ten et al. 2013).

Fig. 4
figure 4

A stacked section contaminated with random interference noise. The deep reflection is concealed by the strong background noise

Full size image
Fig. 5
figure 5

Frequency spectra of the raw and filtered seismogram (trace 6)

Full size image
Fig. 6
figure 6

Singular values of the frequency Spectra of Fig. 4. The blue line represents whole singular values. The red and green lines represent the selected singular values to be used in reconstruction

Full size image
Fig. 7
figure 7

Denoised section via SVD filtering. The first ten singular values are used. The background noise is eliminated completely, but a little loss of geological details occurred

Full size image
Fig. 8
figure 8

Reconstruction of Fig. 6b (from 11th to 15th singular values) shows missed details

Full size image
Fig. 9
figure 9

Reflection signals separated from Fig. 4 using SVD filtering method. The first 15 singular values are used. The weak reflection can be recognized more clearly, and it is convenient to interpret

Full size image
Fig. 10
figure 10

Removed random noise. It is the difference plot between Figs. 4 and 9. No loss of geological information

Full size image
Fig. 11
figure 11

Filtered section by means of a band-pass filter with a parameter of (0, 10, 70, 80) Hz. Background noise remains strong, and events are distorted

Full size image
Fig. 12
figure 12

Difference plot between Figs. 4 and 11. High-frequency components of the coherent signal are attenuated

Full size image
Fig. 13
figure 13

Comparison of the frequency spectra. a The raw data. b SVD-filtered section

Full size image

Figure 14 shows another stacked section consisting of 216 traces with 3000 samples per trace. There is a plenty of random interference noise in the section. Random noise is the main factor that affects the quality of the section. In this case study example, we compare the performance of median filtering, curvelet domain filtering and the presented denoising method for noise removal.

Fig. 14
figure 14

A stacked section contaminated with random interference noise. Random noise is the main factor that affects the quality of the section

Full size image

Figures 15, 16 and 17 show the results after noise removal using, respectively, median filtering, curvelet domain filtering and the presented denoising method. The main drawback of median filtering is its tendency to retain some background noise. Median filtering is least effective in removing background noise (see Fig. 15). The curvelet domain filtering succeeds in suppressing background noise, but it did less well in terms of protection of coherent signals (see Fig. 16). The performance of the presented denoising method is better than that of the curvelet domain filtering in terms of suppression of background noise. Inspection of Figs. 15, 16 and 17 shows that the presented denoising method performs best in terms of protecting coherent signals. The continuity of seismic events in Fig. 17 is better than that in Fig. 16. Figures 18, 19 and 20 show the removed noises (i.e., difference plots, section before minus section after filtering) for median filtering, curvelet domain filtering and the presented denoising method.

Fig. 15
figure 15

Reflection signals separated from Fig. 14 using median filtering. Some background noise remains

Full size image
Fig. 16
figure 16

Reflection signals separated from Fig. 14 using curvelet domain filtering. The curvelet domain filtering succeeds in suppressing background noise, but it did less well in terms of protection of coherent signals

Full size image
Fig. 17
figure 17

Reflection signals separated from Fig. 14 using the presented denoising method. The uncorrelated background noise level has been reduced and trace-to-trace coherency thus enhanced

Full size image
Fig. 18
figure 18

Removed noises (the difference plot between Figs. 14, 15)

Full size image
Fig. 19
figure 19

Removed noises (the difference plot between Figs. 14, 16)

Full size image
Fig. 20
figure 20

Removed noises (the difference plot between Figs. 14, 17)

Full size image

Conclusion

In this paper, we have developed a novel denoising scheme which works by SVD filtering in frequency domain. This approach is better than the traditional band-pass filtering in removing the background noise. Instead of the frequency filtering, we perform the eigenimage filtering of the frequency spectra by selecting singular values to be used in the reconstruction, suppressing the random noise. Compared with the traditional band-pass filtering, the presented method is capable of attenuating the interference noise components within the range of frequency pass band and protects effective signals in high frequency. Through comparison to the median filtering and the curvelet domain filtering, the presented denoising method performs better in removing background noise and protecting reflection events.