Seismic data interpolation plays a fundamental role in seismic data processing, which provides the regularly sampled seismic data for the following workflows like high-resolution processing, wave-equation migration, multiple suppression, amplitude-versus-offset (AVO) or amplitude-versus-azimuth (AVAZ) analysis, and time-lapse studies (Zhong et al., 2016; Chen et al., 2015; Gan et al., 2015a; Wang et al., 2015; Zhang et al., 2016; Chiu, 2014). Three main types of interpolation approaches have been proposed in the literature for handling with the interpolation problem. The first type of approach is based on prediction (Spitz, 1991; Porsani, 1999). A prediction operator is designed from the low-frequency components of seismic data and is applied to high-frequency components. However, the predictive filtering method can only be applied to regularly sampled seismic data. The second type is a transformed domain method (Liu et al., 2015; Candès et al., 2006a; Gan et al., 2015b; Chen et al., 2014b), which is based on compressive sensing theory (Donoho, 2006; Candès et al., 2006b) to achieve a successful recovery using highly incomplete available data (Chen et al., 2014a; Sacchi et al., 1998; Wang, 2003). Compressive sensing (CS) is a relatively new paradigm (Candès et al., 2006b; Liu et al., 2016; Gan et al., 2016; Donoho et al., 2006; Donoho, 2006) in signal processing that has recently received a lot of attention. The theory indicates that the signal which is sparse under some basis may still be recovered even though the number of measurements is insufficient as required by the Shannon's criterion. Naghizadeh and Sacchi (2007) propose a multistep autoregressive strategy which combines the first two types of methods to reconstruct irregular seismic data. The third type is based on the wave equation. This type of method utilizes the inherent constraint of the seismic data from wave equation to interpolate seismic data, thus it depends on the known velocity model, which also becomes its limitation (Fomel, 2003; Canning and Gardner, 1996). Recently, an increasing number of researchers have developed algorithms connecting the interpolation and deblending (Berkhout, 2008) problems for the irregular sampled simultaneous-source data (Li et al., 2013), which provides new recipes for conventional seismic interpolation problem.

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.