Seismic data reconstruction using a multi-step auto-regressive method

The main purpose of exploration seismology is data gathering, data processing, and finally obtaining an interpretable image of subsurface layers. Sometimes, because of problems such as undesirable area topography, instrument defects, and environmental constraints, we have data with missing spatial samples. Reconstruction and recovery of the missing data can be carried out using interpolation and reconstruction methods. There are many reconstruction and interpolation methods. One of the most useful methods to reconstruct missing data is the auto-regressive model. This method refers to the techniques that model the evolution of a signal as a function of its past/future samples (Lau et al., 2002; Takalo et al., 2005). Also, it has a wide range of applications in signal processing including noise suppression (Canales, 1984), parametric spectral analysis (Marple, 1987), and signal interpolation and reconstruction (Sacchi and Ulrych, 1996; Porssani, 1999; Spits, 1991; Naghizade and Sacchi, 2007). The autoregressive reconstruction methods were introduced by Spitz (1991). Spitz (1991) proposed computing prediction filters (autoregressive operators) from low frequencies to predict interpolated traces at high frequencies. This methodology is applicable only if the original seismic section is regularly sampled in space. Conversely, irregularly sampled data can be reconstructed using Fourier methods. In this case, the Fourier coefficients of the irregularly sampled data are retrieved by inverting the inverse Fourier operator with a band limiting and/or a sparsity constraint (Sacchi et al., 1998; Zwartjes and Gisolf, 2006). In this paper, a reconstruction method has been introduced that combines a Fourier-based method and an auto-regressive model to reconstruct the missing data. The method includes a two-stage algorithm. The first step of the proposed algorithm involves the reconstruction of the irregularly missing spatial data on a regular grid at low frequencies using a Fourier-based algorithm called the minimum-weighted norm (Liu and Sacchi, 2004) method. Fourier reconstruction methods are well suited to reconstruct seismic data in the low-frequency (non-aliased) portion of the Fourier spectrum. The reconstruction problem is well-conditioned at low frequencies where only a few wavenumbers are required to honor the data. This makes the problem well-posed; therefore, it is quite easy to obtain a low frequency spatial reconstruction of the data. Seismic data at low frequencies are band-limited in the wavenumber domain. Due to the band-limited nature of the wavenumber spectra at low frequencies, this portion of the data can be reconstructed with high accuracy (Duijndam et al., 1999). Then, prediction filter components are computed for all frequency bands from the low-frequency portion of the reconstructed data using the auto-regressive method. Finally, these prediction filters are used to reconstruct the missing data.  The aforementioned equations show that one can predict the data samples using past/future samples (forward/backward equations). It is important to stress that the technique presented in this paper can only be used to reconstruct data that live on a regular grid with missing observations. The results of the application of the algorithm on both synthetic and real seismic data showed and confirmed the performance of the method.