Introduction

With the increasing demand for exploration accuracy, wide azimuth seismic (WAZ) exploration technology has received more and more attention and development. The seismic data obtained by WAZ acquisition contains rich wavefield information with high illumination (Chen et al., 2014; Bucha, 2017; Chen and Song, 2018; Tian et al., 2017; Zu et al., 2017; Zhang et al., 2019; Lu and Feng, 2017; Kim et al., 2017; Wang et al., 2018; Qu et al., 2016; Gan et al., 2016; Zhang et al., 2018; Chen et al., 2017a). However, due to the influence of surface environments, the shot, receiver, azimuth and offset of the original data tend to be unevenly distributed, which is adverse to the subsequent processing and interpretation such as migration, pre-stack attribute analysis, reservoir prediction and fluid identification. Seismic data interpolation is such a processing step to regularize the irregularly sampled seismic traces onto regular grids (Zhang et al., 2016b; Zu et al., 2016b; Zhou and Li, 2018; Fomel, 2003; Chen et al., 2016c; Zhang et al., 2016c; Chen et al., 2016a; Xiang et al., 2016; Zhang et al., 2016a). In the past decades, a number of methods have been developed to reconstruct the missing seismic traces on regular grids. One widely used seismic reconstruction strategy is to transform the noisy seismic data into different domains in order to concisely represent the signal with a few number of selected components that capture the useful information and recover the missing components. These include the methods based on the Fourier transform (Zhou, 2017; Bai and Wu, 2018), Radon transform (Beylkin, 1987), seislet transform (Zhou and Han, 2018b; Xue et al., 2017; Bai and Wu, 2019; Gan et al., 2015a; Chen and Fomel, 2018; Zhou et al., 2018a; Liu et al., 2016b; Gan et al., 2015b,c), curvelet transform (Zu et al., 2016a; Ma and Dimet, 2009; Zu et al., 2015; Candès et al., 2006), dreamlet transform (Wang et al., 2015), wavelet transform (Liu et al., 2016a; Mousavi and Langston, 2016; Xie et al., 2015; Gilles, 2013; Rioul and Vetterli, 1991), dictionary-learning based adaptive transform (Chen et al., 2016b; Anvari et al., 2017; Siahsar et al., 2017b; Wu and Bai, 2018d; Siahsar et al., 2017a; Zu et al., 2018; Wu and Bai, 2018b; Zu et al., 2019). Another type of methods utilize the predictable property of seismic data. In the prediction-based approaches, a prediction error filter is designed such that the predicted data and the existing data have the minimum misfit by solving a least-squares linear inverse problem (Spitz, 1991; Chen et al., 2016c; Fomel, 2002). Considering the aliasing issue in reconstructing regularly sampled seismic data, a prediction error filter is first estimated from the aliasing-free low-frequency components and then applied to aliased high-frequency components (Naghizadeh and Sacchi, 2007). The wave equation based methods can also be used to reconstruct highly incomplete seismic data. These methods connect the seismic record and the subsurface elastic properties via the elastic wave equation. However, these methods depend on a prior information about the subsurface elastic properties and are strictly not data-driven methods. In addition, these methods suffer from losing computational efficiency when applied in the interpolation methods (Canning and Gardner, 1996; Ronen, 1987; Fomel, 2003). A review of the latest methods on reconstructing low-dimensional seismic data is given in Chen et al. (2019).

However, the amplitude and phase variations of the WAZ seismic data are often not accurately recovered by conventional three-dimensional data reconstruction. More accurate reconstruction results can be obtained by simultaneously interpolating in five dimensions of inline, crossline, time, azimuth and offset because sampling along any particular subset of all dimensions is often less than ideal. The pre-stack seismic gathers processed by the five-dimensional (5D) interpolation have higher quality, which can not only improve the imaging accuracy, but also improve the capabilities of fracture prediction and fluid identification. Therefore, a large number of 5D seismic data reconstruction methods and techniques have emerged in the past decade.

The minimum weighted norm interpolation (MWNI) is a constrained inversion algorithm that was successfully applied to 5D seismic data interpolation (Trad, 2009; Liu and Sacchi, 2004; Trad, 2007). Trad (2007) employed this algorithm to regularize the data in the inline-crossline-azimuth-offset frequency domain and obtained pre-stack interpolated results for the migration input. This work proved that the results after 5D interpolation help to improve the fidelity of the migration, which laid the foundation for the subsequent promotion of five-dimensional seismic data reconstruction. Downton et al. (2008) confirmed the validity of the MWNI for 5D interpolation and found that the 5D interpolation can preserve the amplitude and improve the signal-to-noise ratio. The MWNI did successfully interpolate sparse data and reduced migration artifacts (Trad, 2009), but it had difficulty to deal regularly missing data with spatial aliasing. Chiu (2012) proposed an anti-aliasing MWNI method to improve the deficiency of conventional MWNI in processing aliased data. In addition, some other multidimensional interpolation algorithms are also available for the five dimensions (Chopra and Marfurt, 2013). Jin (2010) proposed a 5D interpolation method based on a damped least-norm Fourier inversion (DLNFI). Benefiting from the use of nonuniform discrete Fourier transform, DLNFI breaks the limitation in MWNI that the input data must be binned into the regular grid.

Another alternative 5D interpolation method utilizes wavefront attributes such as wavefront curvatures and propagation angles. Xie and Gajewski (2017) proposed a wavefront-attribute-based 5D interpolation (5D WABI) via a global optimization strategy instead of pragmatic search approach, which can take advantage of the wide, rich, or full azimuth acquisitions. Application on a synthetic seismic data have shown that the 5D WABI method is better at preserving diffractions than the damped rank-reduction method but at the expense of significantly lower computational efficiency. In recent years, dictionary learning and machine learning are applied to the reconstruction of 5D simple data (Jia et al., 2018; Yu et al., 2015; Jia and Ma, 2017). Data driven tight frame (DDTF) is a kind of dictionary-learning method, which can simultaneously denoise and interpolate 5D seismic data (Yu et al., 2015). In DDTF, a sparsity-promoting algorithm is used to build the dictionary which can represent the observed data and estimate the complete data. Jia and Ma (2017) combined the DDTF with a classic machine learning method named support vector regression (SVR) to optimize the learning, which obtained better performance than the Gauss SVR method. With the continuous improvement of intelligent methods, learning-based 5D data reconstruction will be a hot research topic in the future.

Rank-reduction based algorithms have great potentials for 5D seismic data interpolation. The basic assumption of these methods is that the fully sampled noise-free seismic data can be characterized as a low-rank matrix or tensor and the rank of the matrix or tensor increases when there are missing traces or noise in the data (Zhou et al., 2018b). The missing traces denote the zero-value seismic traces when the seismic data is binned from irregular station positions to regular grids. By solving a low-rank tensor completion problem via convex optimization, the missing traces can be accurately recovered with the information of all dimensions (Kreimer et al., 2013). In addition, this method has stronger denoising ability than the previous methods. From other perspectives of rank minimization, Ely et al. (2015) estimated the complete data tensor via tensor singular value decomposition (SVD) and parallel matrix factorization, respectively. Subsequently, Gao et al. (2017) further developed a new and fast low-rank tensor completion method based on parallel square matrix factorization and advocated that reshaping the complete data tensor into almost square or square matrices can improve the reconstruction quality. However, when the signal-to-noise ratio of the observed seismic data is very low, the Cadzow rank-reduction method via truncated singular value decomposition (TSVD) does not achieve reasonable reconstruction result, i.e., the result may still contains a large amount of residual noise. Chen et al. (2016c) proposed a damped rank-reduction method to further suppress the residual noise by introducing a damping operator to the block Hankel matrix after TSVD. In the 5D reconstruction of field seismic data, the damped rank-reduction method achieved better performance than the Cadzow rank-reduction method (Trickett and Burroughs, 2009; Oropeza and Sacchi, 2011). Throughout the paper, the Cadzow rank-reduction method is referred to as the traditional rank-reduction method.

The basic assumption of the rank-reduction methods is that the Hankel matrix formulated from the useful seismic signal is of low rank and its rank is equal to the number of linear/planar events (or dipping components) in the seismic data (Siahsar et al., 2017c; Wu and Bai, 2019; Oropeza and Sacchi, 2011; Huang et al., 2016; Wu and Bai, 2018c,a; Bai et al., 2019; Wu and Bai, 2018e; Bai et al., 2018b,a; Zhou and Han, 2018a; Chen et al., 2017b). However, in general the real seismic data is composed of nonlinear events, which are more complicated than the linear events. A common strategy for addressing this issue is to implement the algorithm in local patches since the events can be approximately viewed as linear in small patches (Zhang et al., 2017). However, this strategy poses another difficulty of choosing an appropriate rank for each local processing window. A practical implementation of the rank-reduction method may require the predefined rank to be relatively large in order to avoid the damage of useful signal due to inappropriate assumption of the structural complexity. This conservative selection of rank would make the resulted data contain a significant amount of residual noise. In this paper, we develop a relatively adaptive rank-reduction method, by which we can suppress the strong residual noise even when using a large rank. We first calculate a set of optimal weighting coefficients to weight the singular values as a first step by solving an optimization problem that minimizes the Frobenius-norm difference between the approximated signal components and the exact signal components. The rank-reduction method based on the optimal weighting strategy can be further improved when connected with the damped rank-reduction method. The resulted method is named as the optimally damped rank-reduction method, and can potentially provide the best solution to the rank-reduction based high-dimensional seismic reconstruction problem in an adaptive way. We will use comprehensive numerical tests and detailed analysis to demonstrate the superiority of the presented algorithm based on several synthetic examples and a real seismic data.