Recently, rank reduction based seismic interpolation algorithms become popular (Oropeza and Sacchi, 2011; Trickett and Burroughs, 2009). The rank reduction methods for seismic data reconstruction can be divided into two main categories. The first category applies rank reduction to multilevel block Hankel/Toeplitz matrices formed from the entries of the tensor. In other words, multidimensional seismic data are rearranged into a block Hankel or Toeplitz matrix, and a rank reduction algorithm is used to improve the signal-to-noise ratio (SNR) and to reconstruct the data. Such algorithm is usually named as the Cadzow method (Trickett and Burroughs, 2009) or multichannel singular spectrum analysis (MSSA) method (Oropeza and Sacchi, 2011). The other category of rank reduction methods encompasses techniques that are based on dimensionality reduction of multi-linear arrays or tensors. In this case, the multichannel seismic data is viewed as a multi-linear array, and dimensionality reduction techniques are directly applied to the multi-linear array (Gao et al., 2015). For example, Kreimer and Sacchi (2012) adopt the high-order singular value decomposition (HOSVD) to solve the 5D seismic data reconstruction problem in the frequency-space domain.
In this paper, we focus on simultaneous reconstruction and denoising of 3D seismic data using the MSSA algorithm. MSSA is a data-driven algorithm developed from research on alternative tools for the analysis of multichannel time series, which is based on the truncated singular value decomposition (TSVD) (Golub and Loan, 1996) of the Hankel matrix. MSSA is also an extension of singular spectrum analysis (SSA) (Vautard and Ghil, 1989), which is used to analyze 1D time series. Missing traces and random noise increase the rank of the appropriately constructed Hankel matrix that is composed of seismic data at a given frequency slice. To some extent, missing data in each frequency slice performs like random noise at the first iteration of weighted projection onto convex sets (POCS) like framework. Similar to other simultaneous seismic data reconstruction and denoising approaches, MSSA transforms the noisy data with missing traces into a domain where signal and noise are mapped onto separate subspaces and then removes the noise. Weighted POCS-like method, which is widely used for seismic data reconstruction and first introduced to the community of seismic exploration by Abma and Kabir (2006), is in charge of the reconstruction procedure. Many numerical experiments, however, indicate that the random noise cannot be completely removed using the conventional MSSA algorithm and there still exist some reconstruction errors. The main reasons is that the traditional TSVD can only decompose the data into a noise subspace and a signal-plus-noise subspace. Huang et al. (2016) suggest using damped MSSA (DMSSA) algorithm to better decompose the data into the signal subspace and noise subspace for random noise attenuation. In order to overcome the defect mentioned above, we extend the DMSSA algorithm further to simultaneous reconstruction and denoising of 3D seismic data. We first review the theory of traditional MSSA algorithm, then we introduce the proposed DMSSA based reconstruction and denoising framework. Next, we introduce the main components in the Matlab package and correlate the interior functions with the mathematical symbols in the theory sections. Finally, one synthetic and one field data examples are used to demonstrate the performance using the proposed algorithm and the presented open-source Matlab package. The first synthetic example is reproducible from the provided Matlab package. The field data example is, however, not in the public domain, therefore we do not provide the segy file of the field data in the package. But the field data example is just a straightforward application of the Matlab package.