Five-dimensional interpolation in the OVT domain for wide-azimuth seismic data


 To improve the imaging quality of wide-azimuth seismic data and enhance the uniformity of the attributes between adjacent bins, we developed a novel interpolation method in the offset-vector tiles (OVT) domain for wide-azimuth data. The orthogonal matching pursuit (OMP) interpolation method based on the Fourier transform is a frequency-domain processing technique based on discrete Fourier interpolation that achieves the goal of anti-aliasing by extracting the weight factor in the effective band from low-frequency data without aliasing. For data reconstruction, the OMP-based data interpolation technique in the OVT domain comprehensively uses the seismic data in five dimensions: the vertical and horizontal coordinates, time, offset and azimuth. Compared with conventional three-dimensional data interpolation, five-dimensional interpolation in the OVT domain is more accurate and achieves better results in practical applications.


Introduction
The demand for high-resolution oil and gas exploration has prompted continuous progress in seismic technology, and thus, seismic exploration has developed gradually, not only from two-dimensional (2D) methods to three-dimensional (3D) techniques but also from high-precision 3D methods to wide-azimuth (WAZ) methods. The WAZ approach can resolve the problems of complex structures and concealed reservoirs.
Accordingly, from exploration to development, the ability to characterise reservoirs and predict fractures and the accuracy of fluid identification have gradually improved. As a result, because WAZ seismic exploration has more advantages in detecting structures, lithologies and fluids than narrowazimuth (NAZ) seismic exploration, the former has become the mainstream direction of research in the ongoing development of seismic technology (Brian & Tim 2007). WAZ seismic exploration refers to seismic acquisition with horizontal and vertical ratios greater than 0.5. When the ratio is equal to 1.0, the acquisition mode is called fullazimuth seismic acquisition; that is, the data are uniformly sampled in every azimuth. Currently, WAZ acquisition is generally deployed for seismic exploration. Compared with NAZ exploration, WAZ exploration has immeasurable advantages. For example, WAZ seismic acquisition can increase the illumination in all azimuth directions and obtain a more complete seismic wavefield, suppress near-surface scattering noise, thereby enhancing the signal-to-noise ratio of the seismic data, and improve the imaging quality of steeply dipping strata. In the same context, it can be beneficial to study the amplitude variation with the azimuth, which boosts the ability to identify fractures and fluids. With the increasingly extensive application of WAZ seismic exploration both domestically and internationally, the corresponding data processing techniques have also developed rapidly. In particular, the offset-vector tile (OVT) technique is an effective method aimed at WAZ seismic data. Different from conventional processing methods, the OVT technique is relatively new. This approach seeks mainly to convert the conventional processing domain into a five-dimensional (5D) domain (one of time and four of space) that features both offset and azimuth information in addition to the conventional 3D coordinates (one of time and two space variables). The number of folds during WAZ seismic exploration is relatively high, and thus, the densities of shots and traces are relatively large; as a result, the volume of seismic data collected in this way is relatively large, generally reaching the level of TGB. Fortunately, the OVT technique can extract large amounts of information from such high-quality data volumes, especially azimuth and offset information. In the OVT domain, through a series of 5D data processing steps, prestack gathers featuring offset and azimuth information can be ultimately obtained and then used to interpret 5D seismic data, enabling the prediction of reservoirs and the identification of fluids (Heloise 2016).
Many scholars have done a lot of research about interpolation techniques. Generally, they can be grouped into three categories: (i) interpolation methods based on wavefield continuation, such as seismic data reconstruction using a combination of DMO and inverse DMO methods. This kind of method requires the previous knowledge of velocity information of the subsurface medium. If the subsurface information is unknown or the velocity accuracy is low, it will negatively affect the reconstruction results (Ronen 1987). (ii) Interpolation methods based on predictive filtering in the f-x-y domain (Wang 2002;Naghizadeh & Sacchi 2007). Spitz (1991) proposed an anti-aliasing 2D f-x interpolation method with a prediction filter for equally spaced traces exhibiting spatial aliasing. The disadvantage of this kind of method is that it can only reconstruct regularly sampled data and cannot handle irregularly sampled data. (iii) Reconstruction methods based on mathematical transformations, etc. (Trad et al. 2002;Liu & Sacchi 2004;Hennenfent et al. 2010). Four reconstruction formulas for non-uniform sampling of bandwidth-limited signals based on the traditional sampling theorem was proposed (Wang 2003;Yen 2003). The minimum weighted norm interpolation algorithm in the Fourier transform domain was proposed by Liu & Sacchi (2004), but it is only suitable for regular observation systems and cannot effectively suppress spatial aliasing. Compact Fourier interpolation is a minimum mean-square-error interpolation algorithm that can interpolate irregular data to arbitrary locations. Trad (2009) extended this algorithm to interpolate in four spatial dimensions at once to overcome the sampling limitations in 3D land seismic surveys. Xu et al. (2010) proposed the anti-leakage Fourier transform to solve the problem of spectrum leakage caused by irregular sampling. More recently, Kreimer et al. (2013) presented a new algorithm that uses nuclear norm minimisation to solve the tensorcompletion problem in exploration seismology. A focus on matrix and tensor-completion techniques has boosted the development of methods that simultaneously denoise and reconstruct seismic volumes (Kreimer et al. 2013;Gao et al. 2015Gao et al. , 2017. The addition of a robust misfit function into the parallel matrix factorisation algorithm (Xu et al. 2015), a singular value decomposition free tensor-completion technique for seismic data reconstruction is proposed by Carozzi & Sacchi (2019).

Data transformation and characteristics in OVT domain
The concept of OVT was first proposed by Vermeer (1998), after which many researchers (Williams & Jenner 2002) conducted in-depth studies on azimuth gathers and attributes. In this way, the basic theory for the OVT domain processing of WAZ data was formed. When using the OVT technique to process WAZ seismic data, the data need to be first transformed into the OVT domain, where the unique advantages of WAZ acquisition can be used. Four steps are taken to form OVT domain data, as follows. (i) The traces are extracted from all cross-spreads in the geometry, thereby generating a cross-spread gather. Each cross-spread is indicated by its inline and crossline numbers. (ii) The cross-spread gather is divided into small OVTs by every source line and receiver line. Each OVT is composed of a limited range of shots and receivers; therefore, the offset and azimuth are roughly the same (figure 1a). (iii) A coordinate system is established based on each cross-spread. The origin is the intersection between the source line and the receiver line, the horizontal axis is the receiver line and the source line is the vertical axis. The coordinate of an OVT is defined by its projection in the system (figure 1b), and the offset of an OVT is approximately the distance between it and the coordinate origin, while the azimuth angle of an OVT is the angle between the vertical axis and the line connecting the OVT center to the origin. (iv) According to the previous steps, after extracting the cross-spread gathers from the geometry for the whole survey area and collecting all the OVTs with the same coordinates, an OVT gather for the whole area can be produced by lining up the same OVTs with their corresponding inline and crossline numbers, as shown in figure 1c. An OVT gather is a single-fold data set containing both offset and azimuth information, that is, 5D data.
Because the offset and azimuth of each OVT gather are roughly the same, each gather has good energy consistency. Regardless of whether an OVT gather is near-offset, medium-offset or far-offset, there are no problems such as an energy imbalance caused by offset grouping in conventional processing.

Anti-aliasing Fourier interpolation: orthogonal matching pursuit (OMP)
The matching pursuit (MP) algorithm is a transform method based on the theory of compressed sensing and sparse representation. The MP method has been widely used in the fields of time-frequency analysis and the reconstruction of seismic data. One application of this algorithm is MP Fourier interpolation, an effective data interpolation method developed in recent years (Gilles et al. 2010;Jin et al. 2019), which uses Fourier basis functions to model the input data, after which the maximum coefficients are selected and placed in the sparse spectrum at each iteration until the residual value between the sparse representation and the input can be neglected. After implementing an inverse Fourier transform toward the resulting sparse spectrum, the interpolation results at any point can finally be obtained. The key of this method is the selection of the maximum sparse spectrum.
In the real case of data interpolation, we often find that the aliased signal has an energy level comparable to that of the actual signal; therefore, the aliased signal is likely to be selected, which would lead to terrible results. Thus, MP interpolation cannot effectively avoid the risk of aliasing, especially for the high-frequency components of real data. Considering the limitations and solutions of this approach, this paper adopts anti-aliasing Fourier interpolation based on OMP. Similar to MP, OMP also uses Fourier basis functions to represent the data and functions in the frequency domain. Through multiple iterations, one by one, OMP selects a Fourier weight coefficient from the redundant space. To remove the aliased signals, in each iteration, a filter operator is introduced to recalculate the coefficient of the selected Fourier component so that the goal of anti-aliasing is achieved Stephen 2014).
Assuming the input data are represented as , N x is the number of total traces, 0 ≤ x j ≤ x max is the coordinate of the ith trace and the Fourier expansion of the input data f (x j ) is given by The inverse Fourier transform is: where k l is the wavenumber, l is the spatial wavenumber index, k l = 2 l ∕ x max (l = − N k∕ 2 , − N k∕ 2 + 1, ⋯ , N k∕ 2 − 1) and N k is the number of wavenumbers. Increasing x max means increasing the number of traces involved in the calculation. Moreover, it is beneficial to sample the wavenumber domain to improve the interpolation accuracy. In field data processing x max general value is approximately equal to four times the maximum offset of the data.
On the basis of equation (1), the calculation is carried out in the frequency-wavenumber (f-k) domain, where one component can be added each time. After step m, the residual can be expressed as where P l is the index of the coefficient selected at the lth step.
, and the iteration is performed bỹ The index of the maximum coefficient P m can be found by calculating equation (4).
Next, the least-square function is established (equation (5)), and all the coefficients are calculated: (5) The iteration can be stopped when the residual R m+1 (x j ) calculated via equation (5) is small enough.
The specific implementation steps of the anti-aliasing Fourier interpolation algorithm based on OMP are listed here: (i) The discrete Fourier transform is applied in both time and space to the input seismic data to obtain the Fourier spectrum in the f-k domain. (ii) The energy curve of the Fourier spectrum is calculated along different slopes, and the energy curve is used as a weight factor to act on the whole Fourier spectrum. (iii) The maximum energy component of the weighted Fourier spectrum (the effective signal) is selected and then added to the unweighted original Fourier spectrum. After weighting, the energy of the aliased signal decreases because of the small weights, whereas the energy of the effective signal increases because of the larger weights so that the aliasing signal will not be selected. (iv) An inverse Fourier transform is performed on the Fourier spectral component obtained in step (iii), and the result is output to the original input position. (v) The iteration result of step (iv) is subtracted from the input data, and the next iteration is performed. Figure 2 shows the implementation process of antialiasing based on OMP. Figure 2a shows the input data and figure 2b shows the corresponding spectrum. For the spectrum of lower-frequency data without aliasing, the energy curve is calculated along different slopes (figure 2be). The energy spectrum is 'stretched' to a higher-frequency band, where the energy curve is calculated ( figure 2f). This energy spectrum is used as the weight factor and applied to the spectrum of the higher-frequency data before selecting the maximum energy component (figure 2g). Due to the large weights for the real signals and the small weights for the aliasing signals, the risk of selecting the aliasing energy is avoided, thereby achieving the aim of anti-aliasing. For the complicated field data, this necessitates adoption of the windowing strategies. The windowing function method is mainly to select the effective signal, the size of the window generally contains the energy of the low-frequency effective signal as much as possible. The frequency range of the window in figure 2 is 0-30 Hz, and the wave number range is −45-45. The window function is relatively suitable for linear data, but for the actual complex data, because the local linearisation law is present, windowing function strategies are also feasible.
The OMP-based anti-aliasing Fourier interpolation algorithm commonly uses three dimensions: the vertical coordinate, the horizontal coordinate and time. This algorithm can be easily extended to higher dimensions, and the stability of the transformation estimate increases with increasing dimension.

OMP 5D data interpolation in the OVT domain
The OVT domain processing technology is effective for wide-azimuth seismic processing. Since the OVT domain data is 5D data, data interpolation or reconstruction in the OVT domain is 5D. The conventional interpolation and regularisation methods for seismic data are generally implemented in the common offset domain with 3D interpolation. Because the seismic traces in the common offset domain arrive from different directions within the Earth, the data involved in the interpolation calculation also come from different directions. Thus, the azimuth information is lost during interpolation, and the self-similarity of the interpolated data is not good, so spatial aliasing is more likely to occur. However, in actual data processing, it is very difficult to uniformly 532 extract single-fold data in the common offset domain. Artificial offset grouping causes the middle-offset data to be dense, and these data correspond to strong energy in the gather, whereas the near-offset and far-offset data are sparse, so these data correspond to the weak energy in the gather.
The implementation of the OMP anti-aliasing 5D interpolation is similar to that of the OMP 3D algorithm, which is operated in the frequency domain. By establishing the Fourier coefficient, the corresponding algorithm is used to reconstruct the output data. On the basis of the original three dimensions, the 5D interpolation algorithm adds two dimensions: azimuth and offset. Since the sampling interval of the time dimension is regular, 5D interpolation is aimed mainly at the four spatial dimensions that describe seismic data (Daniel 2009;Massimiliano et al. 2010;Satinder & Kurt 2013;Fernanda & Mauricio 2019). The 5D interpolation method comprehensively considers five dimensions, namely the vertical and horizontal coordinates, time, offset and azimuth. By defining the output geometry and rebuilding the output gather, this method can regularise the offset and azimuth data within different tiles. The 5D interpolation algorithm in the OVT domain has high accuracy, and the fidelity of the interpolated data is higher than that of the conventional 3D interpolated data in terms of the amplitude, frequency and phase characteristics.
3D prestack seismic data can be represented by a 5D function f (t, I m , X n , O p , A q ) with independent variables of time t, inline I, crossline X, offset O and azimuth A, where m, n, p, q are the indexes of I, X, O, A, respectively. In the frequency space domain, the function can be expressed as F(I m , X n , O p , A q ), and 5D interpolation interpolates the four dimensions of F(I m , X n , O p , A q ).
Based on the expansion of the Fourier orthogonal basis, the discrete Fourier transform of F(I m , X n , O p , A q ) is expanded to: where the discrete Fourier coefficient is: By performing discrete Fourier transform with respect to X n for the function k I m , we obtained its expansion formula as: where the discrete Fourier coefficient is: ] .
(9) By performing discrete Fourier transform with respect to O p for the function k I m k X n , we obtained its expansion formula as: where the discrete Fourier coefficient is: In the same way, by performing discrete Fourier transform with respect to A q for the function k I m k X n k O p , we obtained its expansion formula as: where the discrete Fourier coefficient is: In these equations, K I m , K X n , K O p , K A q denote the corresponding sampling number in the I m , X n , O p , A q dimensions, respectively, and D I m , D X n , D O p , D A q denote the corresponding sampling intervals.
534 Similarly, substituting equation (14) into equations (11) and (13) subsequently, we obtain: where, k I m k X n k O p k A q is the discrete Fourier transform coefficient. Equation (15) is the discrete Fourier transform of 5D seismic data. Finally, based on the expansion of the Fourier orthogonal basis function, the inverse of the discrete Fourier transform of F(I m , X n , O p , A q ) is: where Therefore, by successively calculating the coefficients of the Fourier expansion function for each of the four dimensions (inline, crossline, offset, azimuth angle), we can do the iterative process of steps (i) to (v) in section 3.1 to obtain the corresponding expansion function with minimised the error function. When obtaining the expansion function, all the expansion coefficients are recalculated after adding a new expansion term. In contrast, the conventional method keeps previously calculated expansion coefficient unchanged and calculates an additional expansion function and its coefficient.
5D interpolation in the OVT domain uses data with the same azimuth angle, thus its accuracy is higher than that of the common offset domain interpolation without considering the azimuth angle. The two dimensions of the offset and the azimuth angle are added, which further improves the stability and the accuracy of the interpolation algorithm. It is convenient and efficient to apply the proposed 5D interpolation algorithm in the OVT domain due to the advantage of preserving the azimuth angle of gathers in the OVT domain. For each input OVT domain gathers, three new OVT domain gathers will be generated, and the total data volume becomes four times the original.

Analysis of practical application and effects
To verify the effectiveness of the proposed OMP interpolation algorithm, we perform 5D interpolation and regularisation in the OVT domain on the WAZ field data of a survey area. The data are 3D wide-azimuth seismic data from the Y88 area of the Qijia-Gulong sag in the northern Songliao Basin of China. The bin size of the data is 5 × 5 m, the folds are 256 and the acquisition aspect ratio is 0.85. Before performing 5D regularisation in the OVT domain, many pre-processing techniques have been adopted to improve the signal-to-noise ratio and resolution of the data, including amplitude-preserving denoising and high-precision static correction, near-surface absorption attenuation compensation, resolution-enhancing deconvolution and the other five processing techniques (shown in figure 3). Then, the preprocessed data are extracted into the OVT domain, where we can perform 5D data regularisation, denoising, migration imaging and other processing. Figure 4a shows the prestack gather before regularisation, which is sorted from prestack 5D seismic data, the other three dimensions are fixed during sorting and the horizontal axis represents the offset (LINE number 800, CDP number 660) from which we can find that the energy of the gather without regularisation is obviously weak in the near-and far-offset regions, while the energy is strong in the middle-offset part, particularly from 1200 to 1400 ms (the red box). Figure 4b displays the regularised gather obtained via the conventional 3D method, which is improved to some extent, and the energy of the events is enhanced in comparison with figure 4a. figures is 1000-1600 ms. The results are clearly superior to those in figure 4a and b. OVT domain regularisation obviously leads to an improved gather with a more balanced distribution of energy. Moreover, the continuity of the events is improved; in particular, the abrupt energy changes between each of the middle-, far-and near-offset parts are invisible. Figure 5 shows a comparison of the stacked profiles before (figure 5a) and after (figure 5b) OVT domain regularisation (LINE number 565, CDP number changes horizontally, as shown in figure 5). Figure 5a reveals numerous missing traces caused by the varying geometry, which would affect the interpretation of follow-up work using the profile. After 5D OVT domain interpolation, the missing traces are recovered well and the events exhibit improved continuity. The amplitudes of the recovered traces also display good consistency with those of the adjacent traces and no obvious distortion is observed. Figure 6 shows a comparison of the amplitude slices before and after OVT domain regularisation, The time corresponding to the amplitude slice is 1200 ms. The amplitude slices after regularisation are clearer, and the details are more  obvious (black ellipse). Since OVT domain 5D interpolation processing needs to consider the orientation information, the algorithm is more complicated and the calculation efficiency will also be reduced. For these actual data, the 3D interpolation processing of the same line data takes 30 minutes, while the OVT domain 5D interpolation processing takes about 45 minutes.

Conclusions
Due to the limitation imposed by variations in the field geometry, the spatial sampling of seismic data is not uniform. To overcome this challenge, this paper proposes a WAZ 5D data interpolation method in the OVT domain based on OMP. The proposed method can be summarised in the following aspects: (1) OVT domain data processing is a novel processing concept for the processing of WAZ data that presents unique advantages. OVT domain data processing retains azimuth information, which is conducive to 5D data interpolation, regularisation and azimuth anisotropy analysis.
(2) The OMP data interpolation algorithm uses the Fourier transform as the basis function. This method is an iterative algorithm that creates frequency slices 537 in the frequency domain, and it has an anti-aliasing capability.
(3) The OMP 5D interpolation algorithm in the OVT domain takes five dimensions into consideration, namely, the vertical and horizontal coordinates, time, offset and azimuth. By defining the output geometry and rebuilding the output gathers, the proposed algorithm can regularise the offset and azimuth within different tiles.
In addition, the algorithm achieves high accuracy, and the amplitude, frequency and phase characteristics of the interpolated data are of high fidelity in comparison with the original data.