3D seismic data reconstruction via DMSSA

Nevertheless, $\tilde{\mathbf{M}}$ is actually still contaminated with residual random noise. Huang et al. (2016) suggest using damped MSSA (DMSSA) algorithm to attenuate the residual random noise that the conventional MSSA fails to. Here, we further extend the DMSSA algorithm to simultaneous reconstruction and denoising of 3D seismic data. Since the missing data at each frequency slice can be regarded as the random noise, the derivation are similar with Huang et al. (2016) except for the extra reconstruction process based on the weighted POCS-like method. Therefore, we conclude the approximation of $\mathbf{S}$ as:

$\displaystyle \mathbf{S}$ $\displaystyle = \mathbf{U}_1^M \Sigma_1^M\mathbf{T}(\mathbf{V}_1^M)^H,$ (8)
$\displaystyle \mathbf{T}$ $\displaystyle =\mathbf{I}-(\Sigma_1^M)^{-K}\hat{\delta}^K.$ (9)

where $\hat{\delta}$ denotes the maximum element of $\Sigma_2^M$ and $K$ denotes the damping factor. It is worth mentioning that the greater the $K$ is, the weaker the damping is, and equation 8 degrades to equation 7 when $K\rightarrow +\infty$.

The revised TSVD for DMSSA algorithm can be represented in operator notation as follows:

$\displaystyle \mathbf{S}=\mathcal{R}_d\mathbf{M}$ (10)

where we use $\mathcal{R}_d$ as the rank reduction operator for DMSSA while $\mathcal{R}$ is used for conventional MSSA.

The filtered data is recovered with random noise attenuated and missing traces reconstructed by properly averaging along the anti-diagonals of the low-rank matrix $\mathbf{S}$ obtained via revised TSVD in DMSSA:

$\displaystyle \hat{\mathbf{D}}=\mathcal{A}\mathbf{S}=\mathcal{A}\mathcal{R}_d\mathbf{M}=\mathcal{A}\mathcal{R}_d\mathcal{H}\mathbf{D}=\mathcal{F}_d\mathbf{D}$ (11)

$\mathcal{A}$ denotes the averaging operator. And the operator $\mathcal{F}_d$ represents DMSSA filter. If $\mathcal{R}_d$ is changed by the traditional TSVD rank reduction operator $\mathcal{R}$, $\mathcal{F}$ can represent conventional MSSA filter similarly:

$\displaystyle \mathcal{F}=\mathcal{A}\mathcal{R}\mathcal{H}.$ (12)

In this paper, we pay our attention to the reconstruction of noisy data, more specifically, those noisy data with very low SNR. Here, we recall the traditional weighted POCS-like algorithm (Oropeza and Sacchi, 2011) as follows:

$\displaystyle \mathbf{D}_n=a_n \mathbf{D}_{obs} + (1-a_n\mathcal{S})\circ \mathcal{F}\mathbf{D}_{n-1},\qquad n=1,2,3,\cdots,n_{max}$ (13)

Substituting the MSSA filter $\mathcal{F}$ with the DMSSA filter $\mathcal{F}_d$, the improved DMSSA algorithm for simultaneous seismic data reconstruction and denoising can be formulated as follows:

$\displaystyle \mathbf{D}_n=a_n \mathbf{D}_{obs} + (1-a_n\mathcal{S})\circ \mathcal{F}_d\mathbf{D}_{n-1},\qquad n=1,2,3,\cdots,n_{max}$ (14)

where $\mathbf{D}_0=\mathbf{D}_{obs}$. $a_n$ is an iteration-dependent scalar that linearly decreases from $a_1=1$ to $a_{n_{max}}=0$. $\mathcal{S}$ denotes the sampling factor. $\mathcal{S}$ equals 1 at the observation point while 0 at the missing traces. $\circ$ denotes the Hadamard product (entry-wise product) of two same size matrices. The algorithm stops when either a maximum number of iterations $n_{max}$ is reached or $\Arrowvert \mathbf{D}_n - \mathbf{D}_{n-1}\Arrowvert_F^2 \leq \epsilon$.

The MSSA based approach is based on the assumption that seismic data is of low rank at each frequency slice in the frequency-space domain. Since seismic data are locally composed of plane-wave components, it intuitively satisfies such assumption. According to our experience, when the SNR is relatively high, in other words with less noise, the performance of MSSA and DMSSA based approaches is similar to each other. However, when the SNR is very low, the proposed DMSSA based approach will outperform the MSSA based approach a lot.