Damped rank-reduction method for 5D seismic data restoration

Huang et al. (2016) introduced a damping factor into TSVD to effectively attenuate the residual random noise for 3D dataset. Here, we further extend this damped algorithm to rank-reduction method for 5D seismic data interpolation in the presence of random noise. Since the missing data at each frequency slice can be regarded as the random noise, the derivation is similar to Huang et al. (2016) except for the extra reconstruction process based on the weighted projection onto convex set (POCS) like method in a level-four block Hankel matrix. Appendix B gives a detailed introduction on how to derive a damped TSVD in order to reduce the residual noise existing in the rank-reduced matrix $\tilde{\mathbf{M}}^{(4)}$ . Here, we just conclude the approximation of $\mathbf{S}$ as follows:

$$\mathbf{S} = \mathbf{U}_1^{M^{(4)}} \Sigma_1^{M^{(4)}}\mathbf{T}(\mathbf{V}_1^{M^{(4)}})^H,$$ (12)

$$\mathbf{T} =\mathbf{I}-(\Sigma_1^{M^{(4)}})^{-N}\hat{\delta}^N,$$ (13)

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

The damped TSVD of rank-reduction algorithm for 5D seismic data interpolation can be represented in operator notation as follows:

$$\mathbf{S}=\mathcal{R}_d\mathbf{M}^{(4)}.$$ (14)

where we use $\mathcal{R}_d$ as the rank-reduction operator using the damped TSVD (equations 12 and 13) while $\mathcal{R}$ is used for conventional TSVD. According to our experience, the damped TSVD is powerful for attenuating random noise while interpolating missing data. However, it is still SVD-based algorithm that takes the largest portion of computational time. The randomized singular value decomposition (RSVD) strategy (Oropeza and Sacchi, 2011; Rokhlin et al., 2009; Halko et al., 2011) can be used to accelerate the proposed method for better efficiency. Specifically according to our numerical tests, RSVD can save at least 45% computational cost compared with the classic TSVD, and cost saving is different with data dimensions, the bigger the dimension is, the faster the RSVD is. Thus RSVD can increase the efficiency significantly.

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 the damped TSVD:

$$\hat{\mathbf{D}}=\mathcal{A}\mathbf{S}=\mathcal{A}\mathcal{R}_d\m... ...M}^{(4)}=\mathcal{A}\mathcal{R}_d\mathcal{H}\mathbf{D}=\mathcal{F}_d\mathbf{D}.$$ (15)

$\mathcal{A}$ denotes the averaging operator. The operator $\mathcal{F}_d$ represents improved rank-reduction filter. If $\mathcal{R}_d$ is substituted by the traditional TSVD operator $\mathcal{R}$, $\mathcal{F}$ can represent conventional rank-reduction filter similarly.

In this paper, we pay our attention to the reconstruction of noisy data. Our final goal is to recover the useful signal from the observed noisy and incomplete data. It is worth mentioning that this type of problem often arises in the field of collaborative filtering ( $\mathbf{CF}$). Similarly to those arising in the field of $\mathbf{CF}$, seismic data reconstruction can be interpreted as a matrix completion problem. Missing traces in seismic data will introduce missing value in the observed data $\mathbf{D}_{k1,k2,k3,k4}$. The corresponding level-four block Hankel matrix $\mathbf{M}^{(4)}$ is expected to have rank $K$. Random noise and missing data will increase the rank of $\mathbf{M}^{(4)}$. Consequently, rank-reduction via the damped TSVD with the damping factor can be implemented to recover the useful signal and to reconstruct the missing records for 5D seismic volumes.

The modified weighted POCS-like algorithm based on conventional rank-reduction method for 5D seismic data interpolation can be shown as $\mathbf{D}_n=a_n \mathbf{D}_{obs} + (1-a_n\mathcal{S})\circ \mathcal{F}\mathbf{D}_{n-1},\quad n=1,2,3,\cdots,n_{max}$. Revised by the improved rank-reduction filter $\mathcal{F}_d$, therefore our proposed improved rank-reduction algorithm for simultaneous 5D seismic data reconstruction and denoising can be formulated as follows:

$$\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}$$ (16)

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. Symbol $\circ$ denotes the Hadamard product (entry-wise product) of two matrices with the same size.

A complete and detailed algorithm workflow of the proposed damped rank-reduction algorithm for 5D seismic data reconstruction and denoising is provided as follows:

The iteration terminates after all $F$ frequencies are finished. For each frequency slice, 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$, where $\parallel \cdot \parallel_F$ denotes the Frobenius norm and $\epsilon$ denotes a small tolerance value.