Introduction

Field obstacles and economic costs in seismic acquisition often result in irregular and incomplete seismic field data (Chiu, 2014). Seismic interpolation is such a processing step to provide 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, thus plays a fundamental role in seismic data processing (Chen et al., 2015; Wang et al., 2010a; Abma and Kabir, 2006; Li et al., 2012; Gan et al., 2016,2015c; Abma and Kabir, 2005; Porsani, 1999).

During the past several decades of continuous research from both academia and industry, there have been many effective methods for seismic data interpolation. In prediction based approaches, a prediction error filter is designed such that the predicted data and the existing data has the minimum misfit by solving a least-squares linear inverse problem (Fomel, 2002; Spitz, 1991). Considering that the irregularly sampled seismic data with large number of missing traces usually suffers from the aliasing problem, the prediction error filter is obtained by solving the inverse problem from low-frequency data that is less aliasing and then applied to high-frequency components (Naghizadeh and Sacchi, 2007). The sparse transform based interpolation assumes that the seismic data is sparse in a certain transformed domain. By iteratively thresholding the transformed domain of the incomplete seismic data, one can gradually attenuate the artifacts in the sparse domain that are caused by the missing traces in the time-space domain, and then gradually recover the missing information (Gan et al., 2015b; Wang et al., 2015a; Gan et al., 2015c; Chen et al., 2014a; Wang et al., 2015b). An important component in this type of methods is the sparse transform. Generally speaking, sparse transforms can be divided into two categories (Chen et al., 2016): analytical transforms, which have fixed basis, and learning-based dictionaries, which iteratively update the basis. A number of fixed basis sparse transforms have been proposed in the literature for restoring seismic data, which include Fourier transform (Sacchi et al., 1998; Chen et al., 2014a; Duijndam et al., 1999), Radon transform (Yu et al., 2007; Wang et al., 2010b; Trad et al., 2002), curvelet transform (Shahidi et al., 2013; Herrmann and Hennenfent, 2008; Liu et al., 2015,2016) and seislet transform (Gan et al., 2015c; Fomel and Liu, 2010; Chen et al., 2014b). The learning-based approach usually utilizes machine learning techniques to construct the dictionary. Another group of methods is based on wave equation. The wave equation theory inherently connects the seismic trace recorded at a specific location from a known shot position to all the other traces assuming the subsurface velocity model is known. Thus, it depends on a prior information that is usually difficult to obtain and it is also limited by the large computational cost aroused from solving the wave equation (Ronen, 1987; Fomel, 2003; Canning and Gardner, 1996).

Recently, rank-reduction based seismic interpolation algorithms become popular (Gao et al., 2015b; Trickett and Burroughs, 2009; Oropeza and Sacchi, 2011). 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 is 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 (Huang et al., 2016,2015; 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., 2015b). Simply speaking, rank-reduction methods are based on the fact that missing traces and random noise will increase the rank of a matrix or tensor. By applying a rank-reduction step, such as truncated singular value decomposition (TSVD), we can intuitively reduce the negative effects caused by the missing traces and random noise. It is worth noting that rank-reduction or low-rank matrix recovery is closely related to sparse transform and compressed sensing (CS) theory because of their utilization of low-dimensional signal structures (Plan, 2011; Candès and Tao, 2010; Recht et al., 2010). Specifically, Recht et al. (2010) discussed the bridge between CS and low-rank matrix recovery. Considering that the multi-dimensional seismic data usually contains significant amount of random noise. A lof of efforts have been put in the area of simultaneous interpolation and denoising of seismic data. Oropeza and Sacchi (2011) introduced a POCS-like iteration into the MSSA algorithm with a variable "alpha" to weight the addition of original data with estimated signal. While this method is effective in attenuating noise while interpolating seismic records, the mathematical relationship between "alpha" and the SNR was not properly understood. Following Oropeza and Sacchi (2011), the weighting strategy was also adopted by Kreimer and Sacchi (2012). Recently the connection between parameter "alpha" and a trade-off parameter "mu" was made in Gao et al. (2015b). Sternfels et al. (2015) solved the simultaneous interpolation and denoising problem via convex optimization, where the problem was formulated as a regularized inverse problem with a sparsity promotion term and solved using the alternating direction method of multipliers (ADMM) algorithm.

In recent years, 5D seismic interpolation has become popular since it takes all physical dimensions of seismic data into consideration for the interpolation, which brings us a stronger constraint when inverting a highly underdetermined inverse problem and thus helps us obtain a much improved reconstruction performance in the case of high missing-traces ratio compared with 2D and 3D interpolation (Kreimer et al., 2013; Liu and Sacchi, 2004; Gao et al., 2015b; Chiu, 2014). The rank-reduction based seismic interpolation has been shown to be an effective method for 5D seismic data interpolation (Gao et al., 2015b). Our recent study shows that when the observed data is extremely noisy, which is often the feature of real seismic data, the reconstructed data using the rank-reduction based algorithm still contains significant residual noise. The residual noise is caused from the fact that the rank-reduced signal space via TSVD is a mixture of signal and noise subspace. We propose an effective method for improving the traditional Cadzow rank-reduction method by modifying the TSVD formula. In order to better decompose the noisy data space into signal and noise subspace, we propose to apply a damping operator to the block Hankel matrix after TSVD. The damping operator is controlled by adjusting an introduced parameter called damping factor. In the special case, when the damping factor is sufficiently large, the improved rank-reduction method reverts the traditional Cadzow rank-reduction method. The proposed method is compared with the traditional method via both synthetic and field 5D seismic dataset and is demonstrated to obtain much better performance.