Adaptive SVD Denoising in Time Domain and Frequency Domain

1. Introduction

Determining accurate exploration locations has always been a common demand in the oil and gas industry. People need clear underground information. However, the actual seismic data collected on site often contains a variety of noises [1]. Random noise is one of the most common interferences. Its existence may interfere with data processing steps such as static correction, velocity analysis and superposition [2,3]. Therefore, suppressing the noise is a basic and important work in seismic data processing.

Various methods have been developed to suppress random noise. According to the difference of effective signal and noise spatial distribution, noise can be located and suppressed [4,5]. Methods such as continuous wavelet transform [6,7], curvelet transform, empirical mode decomposition (EMD) and singular value decomposition (SVD) are widely used. The improved wavelet transform can find the local features of the target signal [8] in the time-frequency domain. However, a series of extended methods based on wavelet transform are not ideal for seismic data with complex events. Curvelet transform with multi-directional and multi-scale characteristics can meet this challenge well. It has achieved good results in suppressing random noise in 3D seismic data [9]. The improved two-dimensional discrete curvelet transform is more flexible for the adjustment of scale and angle threshold [10].

Empirical mode decomposition (EMD) has a good denoising ability [11]. The advantage of EMD is that the de-composition is fully automatic. It means that there is no need to set predefined parameters before decomposition [12]. However, EMD is based on experience, so it lacks theoretical foundation based on mathematical physics. It is also prone to aliasing of modal components, so that the accuracy of the calculation results is not high.

Singular value decomposition (SVD) is also used for denoising due to its strong versatility and low matrix requirements [13,14]. The disadvantage of SVD denoising [15,16] is that it is difficult to meet the accuracy and generalization. Local SVD combines the advantages of f-x deconvolution, which can enhance weak events to remove background noise [17]. Due to the statistical characteristics of different domains, SVD can also be used in the frequency domain. Combined with the projection to convex set method, matrix decomposition is performed in the frequency domain [18]. This method can simultaneously realize the denoising and interpolation of seismic data. SVD is introduced into the non-downsampled shear wave transform coefficients [19,20], which can analyze seismic data well [21]. Single-channel SVD and amplitude ratio are introduced to improve the signal-to-noise ratio (SNR) [22], but the calculation efficiency is not ideal. In ordinary singular spectrum analysis methods, adaptive value selection can reduce the time cost of value selection [23]. The difference between effective signal and noise propagation direction can also be used to distinguish noise. Such as polarization filtering, randon transform [24], predictive de-convolution filtering and other methods.

In this paper, we propose an adaptive SVD denoising method in the time domain and frequency domain (ASTF). Figure 1 illustrates the overall structure of the method. First, the Hankel matrices are formed in the time and frequency domain for SVD respectively. The singular values are adaptively selected according to the second-order difference spectrums of singular values to reconstruct the matrices. Then the two matrices are restored to single-channel data. They are assigned initial weights respectively. Weights are updated iteratively through ternary search. When the peak signal-to-noise ratio (PSNR) of the fusion result reaches the maximum, the iteration is stopped. Finally, the processing results of the two domains are added according to the output weight.

The ASTF method has two main advantages. On the one hand, it benefits from the different advantages of distinguishing signal and noise in the time and frequency domain. Faced with different seismic data, this method can effectively suppress random noise. On.

The other hand, ASTF does not require manual selection of singular values. The effectiveness of the processing results is properly guaranteed

Experiments were conducted on the synthetic datasets and field seismic dataset. We compared ASTF with Experiments were conducted on the synthetic datasets and field seismic dataset. We compared ASTF with three effective methods, including damped multi-channel singular spectrum analysis (DMSSA) [25], empirical mode decomposition (EMD) [26], and space-varying median filter (SVMF) [27]. The results show that our method is more effective.

The rest of this paper is organized as follows. In Section 2, we describe the theoretical basis of ASTF in detail. In Section 3, we introduced the experiment and discussed the results. Finally, Section 4 concludes the paper.

2. Theory

2.1. SVD Denoising

Assuming that the matrix ${\mathit{N}}_{\mathit{m}\times \mathit{n}}$ contains seismic signal and noise, it is expressed as:

N = X + D,

(1)

where ${\mathit{X}}_{\mathit{m}\times \mathit{n}}$ represents pure signal and ${\mathit{D}}_{\mathit{m}\times \mathit{n}}$ represents noise. We perform SVD on ${\mathit{N}}_{\mathit{m}\times \mathit{n}}$ as:

N = UΣV^T,

(2)

where ${\mathit{U}}_{\mathit{m}\times \mathit{m}}$ is an orthogonal matrix composed of eigenvectors NN^T, ${\mathit{V}}_{\mathit{n}\times \mathit{n}}$ is an orthogonal matrix composed of eigenvectors N^TN, Σ is a diagonal matrix composed of singular values in descending order. The first 10% or even 1% of the singular values of Σ can contain more than 90% of the energy information of N [28]. In other words, small singular values mean that their contribution to the effective signal is small.

Assuming that r is the rank of N, the effective rank of the signal matrix X is r′(r′< r). The best estimate of X can be obtained when the r′ largest singular values of N (σ₁ > σ₂ > ⋯> σ_r′) are selected to reconstruct the matrix. The discarded singular values represent the information in the noise matrix D.

2.2. Frequency Domain SVD Denoising

The seismic data N is Fourier transformed to obtain frequency domain data Nf. A certain characteristic frequency of m seismic traces is extracted to form an array (f₁, f₂, …, f_m) to construct a Hankel matrix Hf, which can be expressed as:

$$\mathbf{H}\mathbf{f}=\left[\begin{array}{ccc}{\mathrm{f}}_1& \dots & {\mathrm{f}}_{\mathrm{n}-\mathrm{l}+1}\\ ⋮& & ⋮\\ {\mathrm{f}}_{\mathrm{l}}& \dots & {\mathrm{f}}_{\mathrm{n}}\end{array}\right],$$

(3)

where l = ⌊${\displaystyle \frac{\mathrm{n}}{2}}$ ⌋ + 1. The elements on the anti-diagonal line of the matrix are the same.

If the seismic data contains k linear events without noise, the rank of the corresponding Hankel matrix is k. When the random noise in the seismic data increases, k also increases. SVD can be used to reduce the rank of the matrix to suppress random noise. Similarly, if there are missing seismic events, the corresponding column value in the matrix is 0. The corresponding frequency domain sample value is still 0, when converted to the f-x domain. k will increase if the number of zero-valued samples in the Hankel matrix increases. A low-rank matrix similar to the original matrix without missing traces can be obtained, when the Hankel matrix is reduced by SVD. It is transformed into frequency domain data. Then it performs the inverse Fourier transform to obtain the denoising result.

2.3. Adaptive Selection of Singular Values

In actual reconstruction, the selection of the singular values depend on the correlation of the effective signal. Singular values are usually artificially selected in traditional SVD methods. If these subjective methods are applied, it is difficult to guarantee the effectiveness of the experimental results. In order to meet this challenge, we propose to adaptively select singular values.

In the time domain, we convert the single-channel seismic data into a Hankel matrix Ht. Then the singular value matrix Σ is obtained by performing SVD on Ht. Σ is a diagonal matrix composed of singular values σ_x in descending order:

$$\mathsf{\Sigma }={\left[\begin{array}{ccccc}{\sigma }_1& 0& \dots & 0& 0\\ 0& {\sigma }_2& & ⋮& ⋮\\ ⋮& & \ddots & 0& 0\\ 0& \dots & 0& {\sigma }_r& 0\end{array}\right] }_{mxn}$$

(4)

where r is the rank of matrix Σ, the singular value of the matrix satisfies σ₁ ≥σ₂ ≥⋯≥σ_r ≥0. The effective signal can be regarded as a series of coherent linear events. Random noise is messy because of its complexity. The decay speed of the singular value is from fast to slow in the Σ. The values of the fast decay part are highly correlated with the coherent information of the Hankel matrix. They represent effective signals with great energy. The values of the slow decay part mainly correspond to the weakly correlated random noise. According to the singular value second-order difference spectrum, we select the singular value corresponding to the effective signal. The first difference spectrum B of Σ is:

$${\mathbf{B}}_{\mathrm{i}}={\mathsf{\sigma }}_{\dot{\mathrm{i}}}-{\mathsf{\sigma }}_{\dot{\mathrm{i}}+1},i = 1, 2, \dots , r -1.$$

(5)

Second-order difference spectrum S of Σ is:

$${\mathbf{S}}_{\mathbf{j}}={\mathbf{B}}_{\mathbf{j}}-{\mathbf{B}}_{\mathbf{j}+1}, j=1,2,\dots , r-2.$$

(6)

Because the decay speed of singular value changes from fast to slow. The critical value of fast and slow decay speed is the key to selecting singular value. S does not absolutely change from large to small, how to find the critical value is very important. We propose to apply a sliding time window method to find it. The threshold for judging whether it is a critical value is set as a constant percentage of the average value in the initial time window. With such a parameterization, it will stop searching as soon as the average value of the sliding window is less than the threshold. An advantage of this method is that it is actually a data-driven way to tune the parameter. It is worth noting that the reconstructed matrix is often not a Hankel matrix, which does not meet our expectations. This problem can be solved by a diagonal averaging method [14].

2.4. Adaptive Weight Fusion

Time domain SVD methods have a good denoising effect when they are directly applied on the horizontal events. But if the events are tilted or bent, they need to add cumbersome and complicated corrections. Frequency domain SVD methods can usually overcome these limitations. They avoid assuming that the events of signals are horizontal [29]. When the frequency bandwidth of the signal is wide, however, frequency domain SVD methods will also attenuate part of the signal when suppressing random noise. We propose to adaptively fusing the denoised results in the time domain and frequency domain. It combines the advantages of time domain and frequency domain, which is described in this subsection.

The denoised results in the time domain and frequency domain are assigned initial weights for fusion. The fused matrix TFN can be expressed as:

TFN = ωTN + (1 −ω)FN,

(7)

where ω ∈ [0, 1] is the weight of the time domain denoised matrix, TN is the denoised matrix in the time domain, FN is the denoised matrix in the frequency domain. In the experiment, we discovered that the influence of ω on the peak signal-to-noise ratio (PSNR) of the final result is a convex function, as shown in Figure 2.

We use ternary search to iteratively update ω and stop the iteration when PSNR reaches its maximum value:

$${\omega }^{\ast }=\underset{\omega }{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{a}\mathrm{x}}10{\mathrm{l}\mathrm{o}\mathrm{g}}_{10}\left({\displaystyle \frac{{∥\mathrm{m}\mathrm{a}\mathrm{x}\left(\mathbf{X}\right)-\mathrm{m}\mathrm{i}\mathrm{n}\left(\mathbf{X}\right)∥}_2}{\mathrm{M}\mathrm{S}\mathrm{E}}}\right),$$

(8)

where

$$\min \mathrm{MSE}=\underset{\omega }{\min }{\displaystyle \frac{1}{n}}\sum _{c=1}^n{({\mathrm{TFN}}_{\mathbf{c}}-{\mathrm{X}}_{\mathit{c}})}^2,$$

(9)

where n represents the total number of samples of the seismic data, ${\mathrm{TFN}}_{\mathbf{c}}$ represents the cth sample of the denoised data, ${\mathrm{X}}_{\mathit{c}}$ represents the cth sample of the clean data.

3. Experiments

In this section, synthetic datasets with simulated Gaussian-type noise is used to evaluate how effective ASTF is. We also use a field seismic dataset to illustrate the performance of the proposed method in practice. On the synthetic datasets, we utilize the widely used PSNR metric.

3.1. Adaptive Weight Fusion

This subsection introduces four synthetic datasets and a field seismic dataset used in the experiment. The four synthetic seismic datasets are generated by changing seismic wave speed, stratum and other conditions. Table 1 lists shot points, traces, samples, sampling intervals and source information of these datasets.

3.2. Synthetic Data

In this subsection, we show the denoising performance of our method on synthetic datasets. In order to verify the feasibility and effectiveness of our method, we compared the other three methods on four synthetic datasets with three noise levels (SNR = 1dB, 3dB and 5dB). The results show that our method has obvious denoising effect and can retain most of the effective signals. We first use single-channel seismic data to illustrate the method. Figure 3 shows a clean seismic trace and a noisy data by adding some random noise, the SNR of which is 1dB.

Figure 4 shows the denoising effects of DMSSA, EMD, SVMF and ASTF. Although four methods have achieved obvious results on the noisy data in Figure 3, EMD and SVMF seem to have strong residual noise. ASTF suppresses more noise and its result is relatively smooth. Figure 5 plots the added Gaussian noise and the noise of four different noise suppression methods. EMD and SVMF are less effective in suppressing noise than the other two methods. DMSSA damages a lot of effective signals while suppressing more noise. When ASTF suppresses more noise, it can retain more effective signals.

Table 2 shows the comparison results of the four methods on four complete synthetic datasets with three noise levels. The highest PSNR is indicated in bold. The results can show that ASTF performs well, especially when the signal-to-noise ratio is low, it is obviously better than other methods. Note that for synthetic data, we define the best performance by obtaining the highest PSNR, whereas for field data, the best performance is defined visually.

Figure 6 shows the denoised results in the time domain and frequency domain before fusion. Figure 7 shows the noise suppressed in the time domain and frequency domain respectively. It can be seen from Figure 6 and Figure 7 that both time domain and frequency domain processing can suppress noise. The effective signal processed in the time domain is not smooth enough. Part of the effective signal is lost when processed in the frequency domain. The denoised results of the two domains are fused according to a certain weight, which can well take into account the preservation of the signal and the suppression of the noise.

3.3. Field Data

This subsection introduces the comparison between ASTF and three effective methods (including DMSSA, EMD and SVMF) on a field dataset. Table 1 and Figure 8a show the information of the dataset. The seismic events are damaged to a certain extent below the first arrival. Figure 8b shows the denoised result of ASTF. It can be seen from this that ASFT makes the seismic event clear, and it effectively suppresses the noise.

In order to observe more clearly, we zoom in on the area framed by the black rectangle in Figure 8a,b, and compared the denoising effect with the other three methods. Figure 9a–e show part of the field data and the denoised results of DMSSA, EMD, SVMF and ASTF. The waveform of the signal is smooth if there is no noise. Figure 9a shows that the events are interfered by noise. Compared with the other three results, our method has the smoothest events after denoising. Compared with the other three results, our method has the smoothest event after denoising. The effect can be particularly shown from the lower right corner area of these figures.

Figure 10 plots the noise comparison of four different noise suppression methods. It can be seen that the EMD and SVMF methods are less effective in suppressing noise than the other two methods. DMSSA and EMD lose some effective signals. ASTF achieved a good denoising effect while not obviously attenuating the effective signals.

4. Conclusions

Field seismic data often contains a lot of random noise, which sets high demands for signal enhancement. We propose a denoising method in which the time domain and frequency domain adaptive SVD denoising results are fused. It uses a data-driven model to adaptively select the singular values of the reconstruction matrix. This reduces the time cost of selection and provides an appropriate guarantee for a good denoising effect. Due to the limitations of single-domain processing, the suppressed noise may not be ideal. This fusion method can combine their advantages well to restore a smooth effective signal. This method is compared with DMSSA, EMD and SVMF. Experimental results show that it can obtain the highest peak signal-to-noise ratio. In addition, the experimental results also show that it can effectively restore smooth effective signals without obviously attenuating the signal energy. In the field of image denoising, scholars also use other matrix decomposition methods to suppress noise. Exploring how image denoising migrates to seismic data denoising will be our future work.