Joint Seismic Data Denoising and Interpolation with Double-Sparsity Dictionary Learning

TL;DR

This paper introduces a novel double-sparsity dictionary learning method for simultaneous seismic data denoising and interpolation, improving data quality by adaptively learning from corrupted data and preserving subtle features.

Contribution

It extends previous denoising techniques to handle both interpolation and denoising simultaneously using an adapted double-sparsity dictionary learning approach.

Findings

Preserves subtle seismic features better than traditional methods

Avoids pseudo-Gibbs artifacts in reconstructed data

Outperforms directional multiscale transform methods like curvelets

Abstract

Seismic data quality is vital to geophysical applications, so methods of data recovery, including denoising and interpolation, are common initial steps in the seismic data processing flow. We present a method to perform simultaneous interpolation and denoising, which is based on double-sparsity dictionary learning. This extends previous work that was for denoising only. The original double sparsity dictionary learning algorithm is modified to track the traces with missing data by defining a masking operator that is integrated into the sparse representation of the dictionary. A weighted low-rank approximation algorithm is adopted to handle the dictionary updating as a sparse recovery optimization problem constrained by the masking operator. Compared to traditional sparse transforms with fixed dictionaries that lack the ability to adapt to complex data structures, the double-sparsity dictionary learning method learns the signal adaptively from selected patches of the corrupted seismic data while preserving compact forward and inverse transform operators. Numerical experiments on synthetic seismic data indicate that this new method preserves more subtle features in the dataset without introducing pseudo-Gibbs artifacts when compared to other directional multiscale transform methods such as curvelets.

Figures

Equations

$$\{\hat{\mathbf{D}},$$

$$\mathbf{D}=\mathbf{\Phi A},$$

$$\{\hat{\mathbf{A}},\hat{\mathbf{X}}\}=\underset{\mathbf{A},\mathbf{X}}{\mathrm{argmin}}\Vert \mathbf{Y}-\mathbf{\Phi AX}{\Vert }_F^2$$

$$\text{mbox}s.t.\{\begin{array}{ll}\Vert {\mathbf{x}}_i{\Vert }_1\le t,& \forall i=1,\dots ,R\\ \Vert {\mathbf{a}}_j{\Vert }_1\le p,\Vert {\mathbf{\Phi a}}_j{\Vert }_2=1& \forall j=1,\dots ,L.\end{array}$$

$${\mathcal{I}}_j=\{r|1\le r\le R,x_{jr}\ne 0\}.$$

$${\Vert {\mathbf{Y}}_{{\mathcal{I}}_j}-{\mathbf{\Phi AX}}_{{\mathcal{I}}_j}\Vert }_F^2$$

$$\left\{{\mathbf{a}}_j,{\mathbf{X}}_{j,{\mathcal{I}}_j}^T\right\}=\underset{\mathbf{a},\mathbf{x}}{\mathrm{argmin}}{\Vert {\mathbf{E}}_j-{\mathbf{\Phi ax}}^T\Vert }_F^2\text{mbox}s.t.\{\begin{array}{l}\Vert \mathbf{a}{\Vert }_1\le p\\ \Vert \mathbf{\Phi a}{\Vert }_2=1.\end{array}$$

$$\left\{{\mathbf{a}}_j,{\mathbf{X}}_{j,{\mathcal{I}}_j}^T\right\}=$$

$$\{\hat{\mathbf{s}},\hat{\mathbf{x}}\}=\underset{\mathbf{s},\mathbf{x}}{\mathrm{argmin}}\Vert \mathbf{s}-\mathbf{Dx}{\Vert }_2^2+\mu \Vert \mathbf{x}{\Vert }_1+\lambda \Vert \mathbf{s}-\mathbf{p}{\Vert }_2^2,$$

$$\{\hat{\mathbf{s}},{\hat{\mathbf{x}}}_i\}=\underset{\mathbf{s},{\mathbf{x}}_i}{\mathrm{argmin}}\sum _i\Vert {\mathcal{R}}_i\mathbf{s}-\mathbf{D}{\mathbf{x}}_i{\Vert }_2^2+\sum _i{\mu }_i\Vert {\mathbf{x}}_i{\Vert }_1+\lambda \Vert \mathbf{s}-\mathbf{w}{\Vert }_2^2.$$

$$m_i=\{\begin{array}{ll}1,& i^{\text{mbox}th}\text{element corresponds to available data}\\ 0,& i^{\text{mbox}th}\text{element corresponds to a missing trace}.\end{array}$$

$$\left\{\hat{\mathbf{s}},\hat{\mathbf{A}},{\hat{\mathbf{x}}}_i\right\}=$$

$${\mathbf{x}}_i=\underset{\mathbf{x}}{\mathrm{argmin}}\Vert \mathbf{x}{\Vert }_1$$

$$\text{mbox}s.t.\Vert ({\mathcal{R}}_i\mathbf{m})\odot ({\mathcal{R}}_i\hat{\mathbf{s}}-\mathbf{\Phi Ax}){\Vert }_2^2\le \Vert {\mathcal{R}}_i\mathbf{m}{\Vert }_0^{}{\sigma }^2,\forall i=1,\dots ,R,$$

$$\left\{{\mathbf{a}}_k,{\mathbf{X}}_{k,{\mathcal{I}}_k}^T\right\}=\underset{\mathbf{a},\mathbf{x}}{\mathrm{argmin}}\Vert {\mathbf{M}}_k\odot ({\mathbf{E}}_k-{\mathbf{\Phi ax}}^T){\Vert }_F^2\text{mbox}s.t.\{\begin{array}{l}\Vert \mathbf{a}{\Vert }_1\le p\\ \Vert \mathbf{\Phi a}{\Vert }_2=1,\end{array}$$

$${\mathbf{E}}_k^{\prime }\triangleq {\mathbf{M}}_k\odot {\mathbf{E}}_k+(1-{\mathbf{M}}_k)\odot (\mathbf{\Phi }{\mathbf{a}}_{\text{mbox}old}{\mathbf{x}}_{\text{mbox}old}^T)$$

$$\{{\mathbf{a}}_{\text{mbox}new},{\mathbf{x}}_{\text{mbox}new}\}=\underset{\mathbf{a},\mathbf{x}}{\mathrm{argmin}}{\Vert {\mathbf{E}}_k^{\prime }-{\mathbf{\Phi ax}}^T\Vert }_F^2\text{mbox}s.t.\{\begin{array}{l}\Vert \mathbf{a}{\Vert }_1\le p\\ \Vert \mathbf{\Phi a}{\Vert }_2=1.\end{array}$$

$$\hat{\mathbf{s}}=\underset{\mathbf{s}}{\mathrm{argmin}}\sum _i\Vert {\mathcal{R}}_i\mathbf{s}-{\mathbf{\Phi Ax}}_i{\Vert }_2^2+\lambda \Vert \mathbf{m}\odot (\mathbf{s}-\mathbf{w}){\Vert }_2^2.$$

$$\hat{\mathbf{s}}={\left(\lambda \text{mbox}diag(\mathbf{m})+\sum _i{\mathcal{R}}_i^{†}{\mathcal{R}}_i\right)}^{-1}\left(\lambda (\mathbf{m}\odot \mathbf{w})+\sum _i{\mathcal{R}}_i^{†}\mathbf{\Phi A}{\mathbf{x}}_i\right).$$

$${\mathbf{E}}_k^{\prime }\triangleq {\mathbf{M}}_k\odot {\mathbf{E}}_k+(1-{\mathbf{M}}_k)\odot (\mathbf{\Phi }{\mathbf{a}}_{\text{mbox}old}{\mathbf{x}}_{\text{mbox}old}^T)$$

$$\{{\mathbf{a}}_{\text{mbox}new},{\mathbf{x}}_{\text{mbox}new}\}=\underset{\mathbf{a},\mathbf{x}}{\mathrm{argmin}}{\Vert {\mathbf{E}}_k^{\prime }-{\mathbf{\Phi ax}}^T\Vert }_F^2\text{mbox}s.t.\{\begin{array}{l}\Vert \mathbf{a}{\Vert }_1\le p\\ \Vert \mathbf{\Phi a}{\Vert }_2=1.\end{array}$$

$$\hat{\mathbf{s}}=\underset{\mathbf{s}}{\mathrm{argmin}}\sum _i\Vert {\mathcal{R}}_i\mathbf{s}-{\mathbf{\Phi Ax}}_i{\Vert }_2^2+\lambda \Vert \mathbf{m}\odot (\mathbf{s}-\mathbf{w}){\Vert }_2^2.$$

$$\hat{\mathbf{s}}={\left(\lambda \text{mbox}diag(\mathbf{m})+\sum _i{\mathcal{R}}_i^{†}{\mathcal{R}}_i\right)}^{-1}\left(\lambda (\mathbf{m}\odot \mathbf{w})+\sum _i{\mathcal{R}}_i^{†}\mathbf{\Phi A}{\mathbf{x}}_i\right).$$

$$\{\begin{array}{l}\hat{\mathbf{x}}=\underset{\mathbf{x}}{\mathrm{argmin}}\Vert \mathbf{x}{\Vert }_1\text{mbox}s.t.\Vert \mathbf{m}\odot (\mathbf{w}-\mathbf{\Phi x}){\Vert }_2^2\le \Vert \mathbf{m}{\Vert }_0^{}{\sigma }^2\\ \hat{\mathbf{s}}=\mathbf{\Phi }\hat{\mathbf{x}}\end{array}$$

$$\text{mbox}PSNR=20\log \frac{{\mathbf{s}}_{\text{mbox}max}\sqrt{N_x{N}_z}}{\Vert \mathbf{s}-\hat{\mathbf{s}}\Vert },$$

Full text

Keywords— denoising, interpolation, double-sparsity dictionary learning

A seismic dataset is an ensemble of time-domain wiggle traces collected from an array of receivers. In exploration geophysics a seismic wavefield recording is processed to produce estimates of various properties of the Earth’s subsurface. However, the recorded data may suffer not only from correlated and uncorrelated noise, but also from missing traces due to various constraints (see Figure 1 for examples), such as malfunctioning sensors, limited budget, lack of permission to access the field, etc. Many seismic applications, such as event detection, migration, and inversion have strict requirements on the quality of data. Preconditioning of the recorded data is typically needed for those applications to yield satisfactory results. In this work we concentrate on two major tasks: denoising that attenuates the noise and interpolation that reconstructs the missing traces, and propose a scheme based on dictionary learning that fulfills both of these goals. For 2D seismic data, the interpolation task is actually equivalent to a special case of inpainting in imaging processing.

In recent years, sparse representation of seismic signals via transform-domain methods has attracted considerable attention in seismic data recovery. This model suggests that natural seismic signal are compressible, or well approximated, by a linear combination of only a few atoms from a dictionary. Imposing sparsity constraints on the coefficients of the linear combination could efficiently eliminate the anomalies in the signal recovered from the contaminated seismic data. Good results in both denoising and interpolation have been reported using different transforms/dictionaries, such as wavelets [1], curvelets [2, 3], contourlets [4], seislets [5], etc.

The transform-domain methods above assume specific underlying regularity of the data described by analytic models, resulting in transforms according to implicit but fixed dictionaries that process the data section as a whole. Alternatively, data driven sparse dictionaries [6] or tight frames [7] can be learned directly from the dataset. This ad hoc learned dictionary, which is typically in the form of explicit matrices for small patches instead of full scale dataset, can better adapt the complex data characteristics. For instance, better denoising results obtained using double sparse dictionary learning rather than generic transforms such as curvelet, contourlet, or seislet have been reported in [6, 8]. As an effective dictionary learning algorithm, K-Singular Value Decomposition (K-SVD) [9] has been adopted in seismic denoising problems [10, 11]. However, the major drawback of K-SVD is its high computational complexity. To this end, in this paper we propose to apply the efficient double sparsity dictionary learning approach [12], which further squeezes out the redundancy in atoms of the learned dictionary.

This paper is organized as follows. In Section 1 we revisit the double sparsity dictionary learning. Then we provide the scheme to solve the joint denoising and interpolation problem in Section 2. In Section 3 we give numerical simulations of the proposed scheme and Section 4 concludes this paper.

1 Double Sparsity Dictionary Learning

Given a training set $\mathbf{Y}=[\mathbf{y}_{1},\mathbf{y}_{2},\dots,\mathbf{y}_{R}]\in\mathbb{R}^{N\times R}$ that contains $R$ training samples, the dictionary learning process looks for $\mathbf{D}=[\mathbf{d}_{1},\ldots,\mathbf{d}_{L}]\in\mathbb{R}^{N\times L}$ such that $\mathbf{Y}\approx\mathbf{D}\mathbf{X}$, where $\mathbf{X}=[\mathbf{x}_{1},\mathbf{x}_{2},\dots,\mathbf{x}_{R}]\in\mathbb{R}^{L\times R}$ is a coefficient matrix that is sparse. Here we assume atoms in $\mathbf{D}$ are normalized without loss of generality, i.e., $\forall j\;\|\mathbf{d}_{j}\|_{2}=1$, and $N\leq L$. This dictionary learning process could be accomplished by solving the following tractable convex $\ell_{1}$ minimization problem:

[TABLE]

where $\|\cdot\|_{F}$ is the Frobenius norm. This $\ell_{1}$-norm minimization problem is relaxed from the $\ell_{0}$-norm minimization, which is well known to be NP-hard and can not be solved directly. Under certain conditions [13], problem (1) yields an exact solution to the $\ell_{0}$-norm problem. Hereafter, we use the $\ell_{1}$-norm to measure the sparsity level. Here we follow the convention of K-SVD and divide the seismic data into many small overlapping patches to mitigate the computation burden, the advantages of this approach over the training using whole data is discussed in [6].

Because of the physical properties of seismic wave propagation and reflectivity, atoms with similar geometric features are commonly observed in learned dictionaries [6]. Therefore, it is efficient to employ an off-the-shelf base dictionary to represent the atoms in the learned dictionary. The form we choose is

[TABLE]

where $\mathbf{\Phi}\in\mathbb{R}^{N\times L}$ is a fixed base dictionary and $\mathbf{A}=[\mathbf{a}_{1},\mathbf{a}_{2},\dots,\mathbf{a}_{L}]\in\mathbb{R}^{L\times L}$ is a sparse matrix to be learned in which each column satisfies $\|\mathbf{a}_{i}\|_{1}\leq p$ for some sparsity level $p$. Intuitively, a base dictionary $\mathbf{\Phi}$ that fits well to the expected data regularities should give better final outputs. The advantage of the model in (2) is that the output of dictionary learning is reduced from finding all the entries in a full matrix $\mathbf{D}$ to finding the sparse matrix $\mathbf{A}$, which significantly impacts the efficiency of computation, storage, and transmission. More importantly, with fewer degrees of freedom, such a dictionary model reduces the chance of overfitting the noise in the training set and produces robust results even with a limited number of training examples. For a fixed generic base dictionary $\mathbf{\Phi}$, the dictionary learning process for $\mathbf{D}$ in (2) becomes

[TABLE]

[TABLE]

Because both the actual learned dictionary $\mathbf{A}$ and the representation coefficients $\mathbf{X}$ are sparse matrices, this model is called double sparsity dictionary learning [12] and (1) can be solved by the sparse K-SVD algorithm. To be honest, the name sparse K-SVD is a bit misleading, since the columns in $\mathbf{A}$ are updated one at a time using sparse coding rather than by rank-1 approximation as in the original K-SVD method [9]. To deal with the optimization involving two arguments, the authors of [12] developed a scheme that alternates the updating of $\mathbf{A}$ or $\mathbf{X}$ while keeping the other one fixed.

We now adapt their approach to the missing traces scenario. When updating $\mathbf{a}_{j}$, one column of $\mathbf{A}$, we determine the column index set $\mathcal{I}_{j}$ of the training samples in $\mathbf{Y}$ whose representations use $\mathbf{d}_{j}=\mathbf{\Phi}\mathbf{a}_{j}$ which are found from nonzero entries in $\mathbf{X}$

[TABLE]

Then the objective functional for $\mathbf{a}_{j}$ in (1) can be written as

[TABLE]

where $\mathbf{E}_{j}=\mathbf{Y}_{\mathcal{I}_{j}}-\sum\limits_{i\neq j}\mathbf{\Phi a}_{i}\mathbf{X}_{i,\mathcal{I}_{j}}$ is the residual matrix without the contribution of $\mathbf{d}_{j}$. Therefore, the resulting problem to update $\mathbf{a}_{j}$ and $\mathbf{X}_{j,\mathcal{I}_{j}}$ is given by

[TABLE]

By virtue of Lemma 1 in [12], (6) is equivalent to the following sparse coding problem:

[TABLE]

which updates only those columns within the index set $\mathcal{I}_{j}$ at the $j^{\mbox{th}}$ row of $\mathbf{X}$.

2 Simultaneous Seismic Data Denoising and Interpolation

Besides random noise, the dictionary learning method is applicable to data distortion as well. For instance, distortions caused by missing traces and additive noise can be handled jointly using dictionary learning. The interpolation of missing traces is crucial for many seismic applications, because inadequate or irregularly spaced traces in the acquired seismic dataset could produce strong artifacts in the following seismic processing stages. In practice, trace interpolation, along with denoising, has become one essential step in industrial seismic data preprocessing workflows.

Previously, a variety of methods have been developed for seismic dataset interpolation. Early work proposed a trace interpolation method by wave-equation methods based on the principles of wave physics [14]. Later, methods based on the Fourier transform [15, 16, 17] were adopted to reconstruct irregularly sampled seismic signals. In the past decade, multi-scale transform methods have been widely used to fill gaps among traces based on the sparsity of seismic wave fronts in the transform domain [1, 2, 18, 19, 20, 5]. These methods process the dataset as a whole.

To start with we consider a small $n_{z}\times n_{x}$ patch in the 2D seismic data section which is then vectorized into $\mathbf{p}\in{}^{N}$, where $N=n_{z}n_{x}$. The noisy patch $\mathbf{p}$ is composed of the signal part $\mathbf{s}$ and noise $\mathbf{n}$, i.e., $\mathbf{p}=\mathbf{s}+\mathbf{n}$. For a given dictionary $\mathbf{D}=\mathbf{\Phi}\mathbf{A}$, denoising can be achieved by:

[TABLE]

where $\mu$ and $\lambda$ control the balance among the above three terms: fidelity of the denoised result to the sparse model, sparsity level, and close fit to the original data. For the general case, where $\mathbf{w}\in{}^{N_{z}N_{x}}$ is the vector representation of the entire noisy seismic section, we define a patching operator as $\bm{\mathcal{R}}_{i}\mathbf{w}=\mathbf{p}_{i}\in{}^{N}$ which takes values from the $i^{\mbox{th}}$ noisy patch and then reshapes them into a vector. In terms of the patches, the denoising problem for the whole seismic data section is written as:

[TABLE]

Again, a local optimal solution can be obtained by the alternating optimization over $\mathbf{s}$ and $\mathbf{x}_{i}$. When $\mathbf{s}$ is fixed (9) is a sparse coding problem, when $\mathbf{s}$ is fixed, equation (9) is a sparse coding problem, and when $\mathbf{x}_{i}$ is fixed it simplifies into a Tikhonov-regularized least squares problem which has a closed-form solution.

In denoising problems, dictionary learning approaches typically employ all the available data to train the dictionary. However, when doing interpolation we have to keep track of missing data during the learning process to avoid artifacts, so the basic assumption is that the locations of all missing data are known. Thus, we use a mask operator in the learning process to mute those missing traces. The mask vector is denoted by $\mathbf{m}\in\{0,1\}^{N_{z}N_{x}}$ whose elements are

[TABLE]

With $\odot$ denoting element-wise multiplication between two matrices or two vectors, the joint denoising and interpolation optimization problem becomes

[TABLE]

where the minimization over $\mathbf{A}$ incorporates the dictionary learning step as well.

There are three alternating steps to solve this optimization problem. After initializing $\hat{\mathbf{s}}=\mathbf{w}$ and using a fixed $\mathbf{A}$, the sparse representation basis pursuit denoising (BPDN) [21] problem for each patch $\bm{\mathcal{R}}_{i}\hat{\mathbf{s}}$ becomes

[TABLE]

where the mask $\mathbf{m}$ guarantees that the missing traces are not taken into account and $\sigma^{2}$ here is the variance of the noise assuming additive white Gaussian noise.

Then, in the process of updating each column $\mathbf{a}_{k}$ of the matrix $\mathbf{A}$ using the fixed $\hat{\mathbf{s}}$ and calculated $\mathbf{x}_{i}$, the following problem, which replaces (6), needs to be solved:

[TABLE]

where the matrix $\mathbf{M}_{k}$ collects $\bm{\mathcal{R}}_{i}\mathbf{m}$ in columns for those $i\in\mathcal{I}_{k}$ and it has the same size as $\mathbf{E}_{k}$. Different from (6), this problem is a weighted low-rank approximation problem. Unfortunately, due to the element-wise mask matrix $\mathbf{M}_{k}$, we cannot explicitly find the simple form by Lemma 1 in [12] as we did for (6) and (7). Alternatively, Nati and Jaakkola [22] put forward a simple but effective iterative algorithm that converges to the local minima of the objective function in (13). The algorithm is based on the expectation-maximization (EM) procedure in which the expectation step fills in the current estimate of $\mathbf{\Phi ax}^{T}$ for all missing elements in $\mathbf{M}_{k}\odot\mathbf{E}_{k}$ and the maximization step updates $\mathbf{\Phi ax}^{T}$ from the filled-in version of $\mathbf{M}_{k}\odot\mathbf{E}_{k}$.

Concretely, Algorithm 1 presents the iterative EM-based algorithm that solves (13). Every time $\mathbf{a}$ and $\mathbf{d}$ are estimated, with the names $\mathbf{a}_{\mbox{old}}$ and $\mathbf{x}_{\mbox{old}}$, they are used to fill in $\mathbf{M}_{k}\odot\mathbf{E}_{k}$ by generating a new observation matrix

[TABLE]

in the expectation step. Then, in the maximization step, $\mathbf{a}$ and $\mathbf{d}$ are updated by the filled-in observation matrix $\mathbf{E}_{k}^{\prime}$

[TABLE]

The problem in the form of (19) can be solved with using Lemma 1 in [12] as in sparse K-SVD. The EM procedure converges to a local minimum very quickly, within only a few ($K\approx 5$) iterations.

Finally, with $\mathbf{A}$ and all $\mathbf{x}_{i}$ obtained, the last remaining problem of (11) for the interpolation result $\hat{\mathbf{s}}$ is the least-squares problem

[TABLE]

which has the closed-form solution

[TABLE]

Note that the mask $\bm{\mathcal{R}}_{i}\mathbf{m}$ has been removed in front of the reconstruction misfit $\bm{\mathcal{R}}_{i}\mathbf{s}-\mathbf{\Phi Ax}_{i}$ in (20) since at this point the entire $\mathbf{s}$ is now being restored including the missing traces.

The detailed implementation of the proposed seismic data recovery method can be found in Algorithm 3, where the atom replacing is a trick borrowed from [6] which replaces those duplicated and rarely used atoms in the learned dictionary and in turn improves the efficiency of the dictionary learning.

Concretely, Algorithm 1 presents the iterative EM-based algorithm that solves (13). Every time $\mathbf{a}$ and $\mathbf{d}$ are estimated, with the names $\mathbf{a}_{\mbox{old}}$ and $\mathbf{x}_{\mbox{old}}$, they are used to fill in $\mathbf{M}_{k}\odot\mathbf{E}_{k}$ by generating a new observation matrix

[TABLE]

in the expectation step. Then, in the maximization step, $\mathbf{a}$ and $\mathbf{d}$ are updated by the filled-in observation matrix $\mathbf{E}_{k}^{\prime}$

[TABLE]

The problem in the form of (19) can be solved with using Lemma 1 in [12] as in sparse K-SVD. The EM procedure converges to a local minimum very quickly, within only a few ($K\approx 5$) iterations.

Finally, with $\mathbf{A}$ and all $\mathbf{x}_{i}$ obtained, the last remaining problem of (11) for the interpolation result $\hat{\mathbf{s}}$ is the least-squares problem

[TABLE]

which has the closed-form solution

[TABLE]

Note that the mask $\bm{\mathcal{R}}_{i}\mathbf{m}$ has been removed in front of the reconstruction misfit $\bm{\mathcal{R}}_{i}\mathbf{s}-\mathbf{\Phi Ax}_{i}$ in (20) since at this point the entire $\mathbf{s}$ is now being restored including the missing traces.

The detailed implementation of the proposed seismic data recovery method can be found in Algorithm 3, where the atom replacing is a trick borrowed from [6] which replaces those duplicated and rarely used atoms in the learned dictionary and in turn improves the efficiency of the dictionary learning.

3 Numerical Simulation

The following experiments provide recovery performance results for the double-sparsity dictionary learning method when the seismic dataset has missing traces and is corrupted with additive random noise. Figure 1(a) shows the original dataset provided by BP [23, 24] where the number of receivers is $N_{x}=240$ and each trace has $N_{z}=384$ time samples. Noisy seismic datasets are shown in Figure 1(b,c) with 33% and 50% missing traces whose indices are randomly selected between 1 and 240. Note that all the missing traces have Not-a-Number (NaN) values and their corresponding values in the mask vector $\mathbf{m}$ are set to zeros. For the valid (non-missing) traces, white Gaussian noise with $\sigma=0.1$ is added. The value of $\sigma$ is chosen after we remove the mean value of each trace and then normalize its range to one, i.e., divide by the difference between max and min values. After this normalization the largest absolute value on a valid trace is approximately 0.5.

First, as baseline experiments, the fixed multi-scale contourlet and curvelet transforms are used for seismic dataset recovery (denoising and interpolation jointly). The BPDN method (implemented by the package SPGL1 [21]) is used to find the sparse representation of all the valid traces and then the missing traces are inferred via inverse transform operations as follows

[TABLE]

where $\mathbf{\Phi}$ refers to the dictionary of the contourlet/curvelet synthesis operator. Figure 2 presents the restoration results based on the BPDN method using the contourlet and curvelet transforms for the 33% missing traces case. The performance using contourlets can achieve PSNR = 27.50 dB while using the curvelets can achieve PSNR = 28.12 dB. Still, just like in the pure denoising scenario [6], pseudo-Gibbs artifacts are quite obvious in the recovery results.

Next, recovery experiments are carried out following the procedure in Algorithm 3. In the patch-based interpolation framework, one can fill “holes” whose sizes are smaller than that of the atoms [25]. Therefore, to guard against clusters of missing traces, the patch size is set to a slightly larger size $n_{z}\times n_{x}=24\times 24$, and a non-redundant DCT dictionary $\mathbf{\Phi}$ of size $N\times N=576\times 576$ is selected as the base dictionary. The DCT basis elements are purely real and orthogonal so its computation is very efficient. Currently, the DCT is a widely adopted transform in many well known patch-based image processing schemes, e.g., JPEG and MPEG. Similarly, a total number of 10,000 overlapping patches are randomly selected from the corrupted seismic dataset for dictionary learning and the sparse matrix $\mathbf{A}$ is initialized to $\mathbf{I}_{576\times 576}$. The atom sparsity level $p$ is set to 50. For the case with 33% missing traces, Figure 3(a) shows the non-redundant DCT base dictionary, while the learned sparse matrix $\mathbf{A}$ after $K=20$ training iterations is visualized in Figure 3(b). Based on our simulations, we notice that the essence of the dictionary learning is the “learning” process. The choice of initial dictionary $\mathbf{\Phi}$ does not make a big difference. For example, when we initialize $\mathbf{\Phi}$ with a curvelet dictionary, the improvement in the reconstruction results is negligible. However, the computing time increases significantly, because of the redundancy of the curvelet dictionary. The overall dictionary, $\mathbf{D}=\mathbf{\Phi A}$ of size $576\times 576$, is visualized in Figure 3(c).

Based on this double-sparsity learned dictionary, the recovery result can be obtained by (21). Throughout this paper we define the PSNR as

[TABLE]

where $\mathbf{s}_{\mbox{max}}$ is the maximum possible value of the seismic data after normalization, and $N_{x}$, $N_{z}$ are the number of traces and time samples per trace, respectively. The measured performance has been improved to PSNR = 32.11 dB and PSNR = 30.31 dB, as shown in Figure 4(a) and 4(c) for 33% and 50% missing traces, respectively. The corresponding error panels are shown in Figures 4(b) and 4(d). When compared to the contourlet and curvelet transforms, the double-sparsity result exhibits no pseudo-Gibbs artifacts around the wave fronts. More experiments were performed in which the percentage of missing traces ranges from 10% to 60% and the PSNR performance curves are provided in Figure 5(a). The result with the double-sparsity dictionary learning method based on Algorithm 3, which is a modified version of the sparse K-SVD algorithm, yields significantly better PSNR values than the recovery with fixed transforms. In order to test the performance of the proposed scheme against noise, we also run the recovery simulation for different $\sigma$ values to obtain Figure 5(b). The double-sparsity recovery method provides correct and robust interpolation results for $\sigma$ up to 0.25, considering that the dynamic range of each trace is only 1 after normalization.

4 Conclusions

For seismic datasets contaminated by random noise and missing traces, we presented a double-sparsity dictionary learning scheme to recover the data from these two types of distortions simultaneously. The main contribution of this work lies in the extension of the sparse K-SVD algorithm with a masking operator that tracks the missing data locations during the dictionary learning process. In addition, in order to solve the optimization involving the introduced masking operator, we adopt a weighted low-rank approximation algorithm to handle the dictionary updating. Numerical simulations on a benchmark dataset illustrate the validity of this new approach and its advantages over fixed transform approaches in the sense of yielding restoration with better PSNR and greatly reduced pseudo-Gibbs artifacts.

Reference