3D seismic data denoising using two-dimensional sparse coding scheme

TL;DR

This paper introduces a seismic data denoising method based on two-dimensional sparse coding with learned overcomplete dictionaries, improving noise reduction and computational efficiency in seismic signal processing.

Contribution

It proposes a novel 2D sparse coding scheme using tensor modeling and learning-based dictionary training for seismic data denoising.

Findings

Reduces computational costs in seismic denoising.

Enhances signal-to-noise ratio in synthetic and land survey data.

Demonstrates effectiveness of tensor-based 2D sparse coding.

Abstract

Seismic data denoising is vital to geophysical applications and the transform-based function method is one of the most widely used techniques. However, it is challenging to design a suit- able sparse representation to express a transform-based func- tion group due to the complexity of seismic data. In this paper, we apply a seismic data denoising method based on learning- type overcomplete dictionaries which uses two-dimensional sparse coding (2DSC). First, we model the input seismic data and dictionaries as third-order tensors and introduce tensor- linear combinations for data approximation. Second, we ap- ply learning-type overcomplete dictionary, i.e., optimal sparse data representation is achieved through learning and training. Third, we exploit the alternating minimization algorithm to solve the optimization problem of seismic denoising. Finally we evaluate its denoising performance on synthetic seismic data and land data survey. Experiment results show that the two-dimensional sparse coding scheme reduces computational costs and enhances the signal-to-noise ratio.

Tables

Table 1. Table 1: Comparison of SNR using different approaches.

Models	Original	2DSC	K-SVD	MSSA
Synthetic	0.1403	13.8268	10.2818	9.5052
Field	13.6809	26.2310	24.0547	21.7438

Equations

$${\mathcal{D}}_1\ast {\mathcal{X}}_1+\cdots +{\mathcal{D}}_r\ast {\mathcal{X}}_r=\mathcal{D}\ast \mathcal{X},$$

$$\mathcal{Y}=\mathcal{M}+\epsilon ,$$

$$\begin{array}{rl}f(\hat{\mathcal{D}},\hat{\mathcal{X}},\hat{\mathcal{M}})=& \lambda {\Vert \mathcal{M}-\mathcal{Y}\Vert }_2^2+\sum _{i,j}{\mu }_{ij}{\Vert {\mathcal{X}}_{ij}\Vert }_0+\\ & \sum _{i,j}{\Vert {\mathcal{D}}^H{\mathcal{X}}_{ij}-{\mathcal{M}}_{ij}\Vert }_2^2,\end{array}$$

$$\begin{array}{rl}\underset{\mathcal{D},\mathcal{X}}{\min }& \frac{1}{2}{\Vert \mathcal{Y}-\mathcal{D}\ast \mathcal{X}\Vert }_F^2+\beta {\Vert \mathcal{X}\Vert }_1,\\ \text{s.t. }& {\Vert \mathcal{D}(:,j:))\Vert }_F^2\le 1,j\in \left\{1,2,\cdots ,k\right\},\end{array}$$

$$\underset{\mathcal{X}\in {\mathbb{R}}^{(r\times n\times k)}}{\min }\frac{1}{2}{\Vert \mathcal{Y}-\mathcal{D}\ast \mathcal{X}\Vert }_F^2+\beta {\Vert \mathcal{X}\Vert }_1.$$

$$\underset{\mathcal{X}\in {\mathbb{R}}^{(r\times n\times k)}}{\min }f(\mathcal{X})+\beta g(\mathcal{X}),$$

$${\mathcal{X}}_{p+1}=\text{arg}\underset{\mathcal{X}}{\min }$$

$$\frac{L_{p+1}}{2}{\Vert \mathcal{X}-{\mathcal{X}}_p\Vert }_F^2+\beta g(\mathcal{X}).$$

$$\begin{array}{rl}\underset{\mathcal{D}\in {\mathbb{R}}^{(m\times r\times k)}}{\min }& \frac{1}{2}{\Vert \mathcal{Y}-\mathcal{D}\ast \mathcal{X}\Vert }_F^2,\\ \text{s.t. }& {\Vert \mathcal{D}(:,j:))\Vert }_F^2\le 1,j\in \left\{1,2,\cdots ,k\right\}.\end{array}$$

$$\begin{array}{rl}\underset{{\hat{\mathcal{D}}}^{(l)}}{\min }& \sum _{l=1}^k{\Vert {\hat{\mathcal{Y}}}^{(l)}-{\hat{\mathcal{D}}}^{(l)}\ast {\hat{\mathcal{X}}}^{(l)}\Vert }_F^2,l\in \left\{1,2,\cdots ,k\right\},\\ \text{s.t.}& \sum _{l=1}^k{\Vert {\hat{\mathcal{D}}}^{(l)}(:,j)\Vert }_F^2\le k,j\in \left\{1,2,\cdots ,k\right\}.\end{array}$$

$$\begin{array}{rl}L(\hat{\mathcal{D}},\Lambda )=& \sum _{l=1}^r{\Vert {\hat{\mathcal{Y}}}^{(l)}-{\hat{\mathcal{D}}}^{(l)}\ast {\hat{\mathcal{X}}}^{(l)}\Vert }_F^2+\\ & \sum _{j=1}^r{\lambda }_j(\sum _{l=1}^k{\Vert {\hat{\mathcal{D}}}^{(l)}(:,j)\Vert }_F^2-k),\end{array}$$

$$\begin{array}{r}SNR=10{\text{log}}_{10}\frac{{\Vert \mathcal{M}\Vert }_2^2}{{\Vert \mathcal{M}-\mathcal{Y}\Vert }_2^2}.\end{array}$$

Full text

1 Introduction

Underground image construction and rock property estimation are important tasks for seismic exploration. However, seismic data may be subject to too serious noise contamination to carry out further processing and interpretation. In many cases, therefore, noise mitigation is essential for extracting useful information from the raw measurements (Zhai,, 2014).

Transform-based function method is one of the most commonly used techniques for seismic denoising. This kind of methods commonly require the conversion of the data into different domains by using methods such as the Fourier transform, the Radon transform, and the curvelet transform (Dashtian and Sahimi,, 2015; Wang et al.,, 2013). However, the bases of the above transforms are pre-chosen and fixed during the process, so they lack of adaptability to dynamic data structures. Tang et al., (2012) introduce learning-type overcomplete dictionaries in which sparse data representation is achieved through learning and training the raw seismic data. As an effective dictionary learning algorithm, K-Singular Value Decomposition (K-SVD) has been applied to the problem of seismic denoising (Beckouche and Ma,, 2014). However, the main drawback is its high computational complexity. In addition, as these algorithms mainly deal with vector-valued input data, the vectorization preprocess would break apart the local proximity and destruct the object structures of seismic data. Oropeza and Sacchi, (2011) present a rank reduction algorithm for noise attenuation of seismic records which adopts the multichannel singular spectrum analysis (MSSA). But this will cause strong residual noise and thus affect the following processing and interpretation tasks (Chen et al.,, 2016).

The basic idea that motivates us to address these challenges lies in two aspects: 1) tensor representation is able to preserve the local proximity and to capture the elementary object structures of seismic data; and 2) we exploit tensor-linear combinations to approximate the seismic data. Then the number of required bases can be significantly reduced, which can reduce the computational complexity.

In this paper, we apply a seismic data denoising method based on learning-type overcomplete dictionaries which adopts the two dimensional sparse coding (2DSC) (Jiang et al.,, 2017; Liu et al.,, 2015). In the method, the input seismic data and dictionary are represented as third-order tensors, and the tensor-linear combinations (Kilmer et al.,, 2013) are used for data approximation. Figure 1 shows the tensor-linear combination. Alternating minimization algorithm which consists of a sparse coding learning step and a dictionary learning step can be used to solve the optimization problem. For sparse coding, we apply an iterative shrinkage thresholding algorithm (Xu and Yin,, 2014) based on the tensor-product, which is directly implemented in the tensor space. For dictionary learning, we show that it can be solved efficiently by transforming to a Lagrange dual problem in the frequency domain.

Therefore, there are three major advantages of our scheme: 1) the size of dictionary can be significantly reduced since the tensor-linear combination are used for seismic data approximation; 2) shifting invariant of tensor-product means that the seismic data can be generated by the shifted versions of bases without explicitly storing them; 3) the bases of dictionary can be trained and adaptively driven by target data, then it is able to develop dictionary transform-based functions suitable for special seismic signals.

NOTATION AND PRELIMINARY

A third-order tensor is denoted as $\mathcal{Y}\in\mathbb{R}^{m\times n\times k}$. The discrete Fourier transform (DFT) along the third dimension of $\mathcal{Y}$ is denoted as $\widehat{\mathcal{Y}}$. The $\ell_{1}$-norm and Frobenius norm are denoted as $\left\|\mathcal{Y}\right\|_{1}=\sum_{i,\,j,\,k}\left|\mathcal{Y}(i,\,j,\,k)\right|$, and $\left\|\mathcal{Y}\right\|_{F}=(\,\sum_{i,\,j,\,k}\sqrt{\mathcal{Y}(i,\,j,\,k)^{2}}\,)$. The superscript $H$ of $\mathcal{D}^{(l)^{H}}$ represents the conjugate transpose. Furthermore, we need two definitions.

Definition1: The tensor-product $\mathcal{Y}=\mathcal{D}*\mathcal{X}$ of $\mathcal{D}\in\mathbb{R}^{m\times r\times k}$ and $\mathcal{X}\in\mathbb{R}^{r\times n\times k}$ is a tensor $\mathcal{Y}\in\mathbb{R}^{(m\times n\times k)}$, where $\mathcal{Y}(i\,,j,\,:)=\sum_{q=1}^{r}\mathcal{D}(i\,,q,\,:)*\mathcal{X}(q\,,j,\,:)$ and the $*$ denotes the circular convolution operation.

Definition2: The tensor-linear combination of the tensor bases $\left\{\mathcal{D}_{j}\right\}_{j=1}^{r}\in\mathbb{R}^{m\times 1\times k}$ with the corresponding tensor coefficient $\left\{\mathcal{X}_{j}\right\}_{j=1}^{r}\in\mathbb{R}^{1\times 1\times k}$ are defined as:

[TABLE]

where $\mathcal{D}\in\mathbb{R}^{(m\times r\times k)}$ with $\mathcal{D}(:\,,j,\,:)=\mathcal{D}_{j}$, and $\mathcal{X}\in\mathbb{R}^{(m\times r\times k)}$.

Problem Formulation

We use the following model to describe the noisy seismic data:

[TABLE]

where $\mathcal{M}$ denotes the clean seismic data to be estimated, $\mathcal{Y}$ is the noisy seismic data and $\varepsilon$ is the random noise. To reduce noise $\varepsilon$ is to estimate $\mathcal{M}$ as accurate as possible from the acquired noisy data $\mathcal{Y}$.

Since seismic data sets are usually so complex with different geometrical features, we use the learning-based sparse dictionary to represent them. This dictionary is overcomplete and its elements are called atoms. In this overcomplete space domain, we adaptively search for an optimal composition of base atoms to sparsely represent the target signal through learning and training (Zhu et al.,, 2015).

Divide the seismic data $\mathcal{M}$ into several sub-blocks $\mathcal{M}_{ij}$ $(i,\,j\in\Omega)$ with fixed size. Assume each sub-block $\mathcal{M}_{ij}$ can be represented by a linear combination of several atoms extracted from the fixed dictionary $\mathcal{D}$ and the corresponding coefficients are denoted by $\mathcal{X}_{ij}$. Then we can formulate the problem with a general form, to minimize the following function with three penalties (Tang et al.,, 2012):

[TABLE]

where $\lambda$ and $\mu_{ij}$ denote penalty functions weights. If given $\widehat{\mathcal{M}}$ and $\widehat{\mathcal{D}}$ we can deduce $\widehat{\mathcal{X}}$ by minimizing equation (3). Once $\widehat{\mathcal{X}}$ is deduced, then $\widehat{\mathcal{M}}=\mathcal{D}^{H}\widehat{\mathcal{X}}$ can be obtained and dictionary $\widehat{\mathcal{D}}$ can be also updated by minimizing function (3).

Algorithm 1 Algorithm for 2DSC

  Input: seismic data: $\mathcal{Y}\in\mathbb{R}^{m\times n\times k}$, the number of atoms: $r$, sparsity regularizer: $\beta>0$, maximum iterative steps: num,

01: Initialization: Randomly initialize $\mathcal{D}\in\mathbb{R}^{(m\times r\times k)}$ , $\mathcal{X}:=0\in\mathbb{R}^{(m\times r\times k)}$, and Lagrange dual variables $\lambda$, set $\mathcal{C}_{1}=\mathcal{X}_{0}\in\mathbb{R}^{(r\times 1\times k)}$, $t_{1}=1$,

02: for iter = 1 to num **do

**03:  //**Tensor Coefficients Learning

**04:  Set $L^{P}=\eta^{p}(\sum_{l=1}^{k}\left\|\widehat{\mathcal{D}}^{(l)^{H}}\widehat{\mathcal{D}}^{(l)}\right\|_{F})$,

05:  Compute $\bigtriangledown f(\mathcal{C}_{p})$,

06:  Compute $\mathcal{X}_{p}$ via $\mathbf{Prox}_{\beta/L_{p}}(\mathcal{C}_{p}-\frac{1}{L_{p}}\bigtriangledown f(\mathcal{C}_{p}))$,

07:  $t_{p+1}=\frac{1+\sqrt{1+4t_{p}^{2}}}{2}$,

08:  $\mathcal{C}_{p+1}=\mathcal{X}_{p}+\frac{t_{p-1}}{t_{p+1}}(\mathcal{X}_{p}-\mathcal{X}_{p-1})$,

09:  //**Tensor Dictionary Learning

**10:  $\widehat{\mathcal{Y}}=\text{fft}(\mathcal{Y},[\,],3)$, $\widehat{\mathcal{X}}=\text{fft}(\mathcal{X},[\,],3),$

11:  for $l$ = 1 to $k$ do

12:    Solving $\Lambda$ by Newton’s method,

13:    $\widehat{\mathcal{D}}^{(l)}=(\widehat{\mathcal{Y}}^{(l)}\widehat{\mathcal{X}}^{(l)H})(\widehat{\mathcal{X}}^{(l)}\widehat{\mathcal{X}}^{(l)H}+\Lambda)^{-1}$,

14:  **end for

**15:  $\mathcal{D}=\text{ifft}(\widehat{\mathcal{D}},[\,],3)$,

16: **end for

Output:** $\mathcal{D},\mathcal{X}.$

2DSC model

In this section, we apply the two-dimensional sparse coding to dictionary learning algorithm. We represent the training sample of 3-D seismic data as a third-order tensor $\mathcal{Y}\in\mathbb{R}^{m\times n\times k}$, then the adaptive dictionary learning process can be described as the following problems:

[TABLE]

where $\mathcal{D}\in\mathbb{R}^{(m\times r\times k)}$ is the tensor dictionary where each lateral slice $\mathcal{D}(:,j,:)=\mathcal{D}_{j}$ is a basis, $\mathcal{X}\in\mathbb{R}^{(r\times n\times k)}$ is the tensor coefficient. The parameter $\beta$ balances the approximation error and the sparsity of the tensor coefficients, and $r$ is the number of atoms.

1.1 Solution

Problem (4) can be solved by an iterative method that alternates between sparse tensor coefficients of the examples based on the current dictionary and a process of updating the dictionary atoms to better fit the data (Algorithm 1). The update of the dictionary is combined with an update of the sparse representations, thereby accelerating convergence. The sparse coding learning step and the dictionary learning step are showd as follows:

1.1.1 Learning Tensor Coefficients

Given the dictionary $\mathcal{D}\in\mathbb{R}^{(m\times r\times k)}$, solving the tensor sparse representations of seismic data $\mathcal{X}\in\mathbb{R}^{(r\times n\times k)}$ are converted to the following problem as:

[TABLE]

Rewrite (5) as:

[TABLE]

then (6) can be sloved by Iterative Shrinkage Thresholding algorithm based on Tensor-product (ISTA-T), which is similar to FISTA (Beck and Teboulle,, 2009). Then at the $p+1$-th iteration, $B_{p}+1$ can be updated by:

[TABLE]

1.1.2 Tensor Dictionary Learning

For learning the dictionary $\mathcal{D}$ while fixed $\mathcal{X}$, the optimization problem is:

[TABLE]

Decompose (8) into $k$ nearly-independent problems by DFT as follows:

[TABLE]

Then, we adopt the Lagrange dual for solving (9) in frequency domain and the Lagrangian of (9) is:

[TABLE]

where $\lambda_{j}\geq 0$, $j\in\left\{1,\,2,\,\cdots,\,k\right\}$ is a dual variable, and $\Lambda=diag(\lambda)$. The advantage of Lagrange dual is that the number of optimization variables is much smaller than which in the primal problem (8).

EVALUATION

In this section, we use a synthetic example and a field data example to compare the denoising performance of the 2DSC scheme whih other two methods. To numerically test the performance and to quantify the comparison, we simulate noisy data by adding Gaussian white noise and define the criterion for comparison as signal-to-noise ratio (SNR):

[TABLE]

Here, $\mathcal{M}$ is noise-free data, $\mathcal{Y}$ is the denoised data.

1.2 Synthetic data example

To evaluate the method effectiveness, we use a Ricker wavelet with central frequency of 60 Hz to generate a simple 3-D model in which there are two dipping planes located around 300 and 900 m on the vertical axis, respectively. The clean and noisy data are shown in Figure 2(a) and 2(b), respectively. After applying the 2DSC method, the noise suppression section is shown in Figure 2(c). These three figures show that this method can effectively suppress noise. For comparison, we also used the K-SVD and MSSA methods to process the model. The comparison of SNR is shown in the first row of Table 1. The 2DSC can get better results considering the preservation of useful signals and removal of noise. And its calculated SNR also is the highest among all the three approaches.

1.3 Filed data example

To verify the practicality of this method, it was applied to F3 land seismic data. The original data is shown in Figure 3(a). The simulated noisy data with Gaussian white noise is shown in Figure 3(b). The denoised results from the three different approaches are shown in Figure 4. The comparison of SNR is shown in the second row of Table 1. From the comparison of denoised results, we conclude that, while all of the three methods can obtain acceptable results, the proposed 2DSC method produces the cleanest image (Figure 4(a)) with the least amount of noisy energy left. The SNR comparison also confirms the previous observation as 2DSC-based method can obtain the highest SNR.

CONCLUSION

In this paper, we applied a dictionary learning algorithm to seismic denoising which adopts 2DSC scheme. Alternating minimization algorithm which consists of a sparse coding learning step and a dictionary learning step was used to solve the problem. We evaluated its denoising performance on synthetic seismic data and land data survey. The results show that this method is a good choice for suppressing noise as it can reduce computational costs and enhance the signal-to-noise ratio.

ACKNOWLEDGMENTS

We would like to acknowledge financial support from the National Natural Science Foundation of China (Grant No.U1562218)