TY - JOUR
AU - Huang, Yiming
AB - Abstract 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. wide azimuth, OVT domain, five-dimensional, interpolation, orthogonal matching pursuit 1. 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 narrow-azimuth (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 full-azimuth 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 tensor-completion 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.2015, 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). 2. 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. Figure 1. Open in new tabDownload slide OVT partitioning in a cross-spread and an OVT gather in plan view. (a) OVT partitioning in a cross-spread. (b) OVT partitioning method. (c) OVT gather in plan view. Figure 1. Open in new tabDownload slide OVT partitioning in a cross-spread and an OVT gather in plan view. (a) OVT partitioning in a cross-spread. (b) OVT partitioning method. (c) OVT gather in plan view. 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. 3. 5D data interpolation in the OVT domain 3.1. 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 (Ali et al.2010; Stephen 2014). Assuming the input data are represented as |$f( {{x_j}} )( {j = 0, \cdots ,{N_x} - 1} )$|⁠, |$\ {N_x}$| is the number of total traces, |$0 \le {x_j} \le {x_{\rm {max}}}$| is the coordinate of the ith trace and the Fourier expansion of the input data |$f( {{x_j}} )$| is given by $$\begin{equation}f( {{x_j}}) \approx \sum\limits_l {\tilde{f}} (l){\rm{exp}}\left( {\frac{{i2\pi l{x_j}}}{{{x_{\rm {max}}}}}} \right).\end{equation}$$(1) The inverse Fourier transform is: $$\begin{equation}\tilde{f}(l) \approx \sum\limits_{j = 0}^{{N_x} - 1} {f( {{x_j}})} {\rm{exp}}\left( {\frac{{ - i2\pi l{x_j}}}{{{x_{\rm{max}}}}}} \right),\end{equation}$$(2) where |$\ {k_l}\ $| is the wavenumber, |$l\ $|is the spatial wavenumber index, |$\ {k_l} = {\hbox{${2\pi l}$} \!\mathord{/ {\vphantom {{2\pi l} {{x_{\rm{max}}}}}}} \!\hbox{${{x_{\rm{max}}}}$}}$||$( {l = - {\hbox{${{N_k}}$} \!\mathord{/ {\vphantom {{{N_k}} 2}}} \!\hbox{$2$}}, - {\hbox{${{N_k}}$} \!\mathord{/ {\vphantom {{{N_k}} 2}}} \!\hbox{$2$}} + 1, \cdots ,{\hbox{${{N_k}}$} \!\mathord{/ {\vphantom {{{N_k}} 2}}} \!\hbox{$2$}} - 1} )$| and |$\ {N_k}$| is the number of wavenumbers. Increasing |${x_{\rm{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_{\rm{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 $$\begin{equation}{R^m}( {{x_j}}) = f( {{x_j}}) - \sum\limits_{l = 0}^{m - 1} {\tilde{f}( {{P_l}})} {\rm{exp}}\left( {\frac{{i2\pi {P_l}{x_j}}}{{{x_{\rm{max}}}}}} \right),\end{equation}$$(3) where |${P_l}\ $|is the index of the coefficient selected at the lth step. After initialising |$m = 0$|⁠, |${R^0}( {{x_j}} ) = f({x_j})$|⁠, and the iteration is performed by $$\begin{equation}\tilde{R}( l) = \sum\limits_{j = 0}^{{N_x} - 1} {{R^m}( {{x_j}})} {\rm{exp}}\left( {\frac{{ - i2\pi l{x_j}}}{{{x_{\rm{max}}}}}} \right).\end{equation}$$(4) 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: $$\begin{equation}\phi( {m + 1}) = {\sum\limits_j {\left[ {f( {{x_j}}) - \sum\limits_{l = 0}^m {\tilde{f}\left( {{P_l}} \right)} {\rm{exp}}\left( {\frac{{i2\pi {P_l}{x_j}}}{{{x_{\rm{max}}}}}} \right)} \right]} ^2}.\end{equation}$$(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: 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. 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. 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. 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. 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 anti-aliasing 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 2b–e). 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. Figure 2. Open in new tabDownload slide The f-k spectra and energy curves for anti-aliasing Fourier interpolation based on OMP. (b) f-k spectrum of the input seismic data. (c) Low-frequency spectrum without aliasing. (d) Energy distribution along different slopes in the f-k spectrum. (e) Calculating the energy curve along different slopes. (f) Energy curve for high-frequency data after stretching the spectrum. (g) Applying weights to the high-frequency spectrum. (h) The seismic data after interpolation. Figure 2. Open in new tabDownload slide The f-k spectra and energy curves for anti-aliasing Fourier interpolation based on OMP. (b) f-k spectrum of the input seismic data. (c) Low-frequency spectrum without aliasing. (d) Energy distribution along different slopes in the f-k spectrum. (e) Calculating the energy curve along different slopes. (f) Energy curve for high-frequency data after stretching the spectrum. (g) Applying weights to the high-frequency spectrum. (h) The seismic data after interpolation. 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. 3.2. 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 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: $$\begin{equation}F({I_m},{X_n},{O_p},{A_q}) = \sum\limits_{{k_{{I_m}}} = 0}^{{K_{{I_m}}}} {{\beta _{{k_{{I_m}}}}}\centerdot \exp \big[ {{\rm{i}}2\pi {k_{{I_m}}}{I_m}/({K_{{I_m}}}{D_{{I_m}}})} \big]} ,\end{equation}$$(6) where the discrete Fourier coefficient is: $$\begin{eqnarray}{\beta _{{k_{{I_m}}}}} &&= \frac{1}{{{K_{{I_m}}}{D_{{I_m}}}}}\sum\limits_{{k_{{I_m}}} = 0}^{{K_{{I_m}}}} {F({I_m},{X_n},{O_p},{A_q})}\nonumber\\ &&\quad\centerdot \exp \big[ { - {\rm{i}}2\pi {k_{{I_m}}}{I_m}/({K_{{I_m}}}{D_{{I_m}}})} \big] .\end{eqnarray}$$(7) By performing discrete Fourier transform with respect to |${X_n}$| for the function |${\beta _{{k_{{I_m}}}}}$|⁠, we obtained its expansion formula as: $$\begin{equation}{\beta _{{k_{{I_m}}}}}({X_n},{O_p},{A_q}) = \sum\limits_{{k_{{X_n}}} = 0}^{{K_{{X_n}}}} {{\chi _{{k_{{I_m}}}{k_{{X_n}}}}}\centerdot \exp \big[ {{\rm{i}}2\pi {k_{{X_n}}}{X_n}/({K_{{X_n}}}{D_{{X_n}}})} \big]} ,\end{equation}$$(8) where the discrete Fourier coefficient is: $$\begin{equation}{\chi _{{k_{{I_m}}}{k_{{X_n}}}}} = \frac{1}{{{K_{{X_n}}}{D_{{X_n}}}}}\sum\limits_{{k_{{X_n}}} = 0}^{{K_{{X_n}}}} {{\beta _{{k_{{I_m}}}}}\centerdot \exp \big[ { - {\rm{i}}2\pi {k_{{X_n}}}{X_n}/({K_{{X_n}}}{D_{{X_n}}})} \big]} .\end{equation}$$(9) By performing discrete Fourier transform with respect to |${O_p}$| for the function |${\chi _{{k_{{I_m}}}{k_{{X_n}}}}}$|⁠, we obtained its expansion formula as: $$\begin{eqnarray}{\chi _{{k_{{I_m}}}{k_{{X_n}}}}}({O_p},{A_q}) &&= \sum\limits_{{k_{{O_p}}} = 0}^{{K_{{O_p}}}} {{\xi _{{k_{{I_m}}}{k_{{X_n}}}{k_{{O_p}}}}}}\nonumber\\ &&\quad\centerdot \exp \big[ {{\rm{i}}2\pi {k_{{O_p}}}{O_p}/({K_{{O_p}}}{D_{{O_p}}})} \big] ,\end{eqnarray}$$(10) where the discrete Fourier coefficient is: $$\begin{eqnarray}{\xi _{{k_{{I_m}}}{k_{{X_n}}}{k_{{O_p}}}}}\!\! &&= \frac{1}{{{K_{{O_p}}}{D_{{O_p}}}}}\sum\limits_{{k_{{O_p}}} = 0}^{{K_{{O_p}}}} {{\chi _{{k_{{I_m}}}{k_{{X_n}}}}}}\nonumber\\ &&\quad\centerdot \exp \big[ {{\rm{ - i}}2\pi {k_{{O_p}}}{O_p}/({K_{{O_p}}}{D_{{O_p}}})} \big] .\end{eqnarray}$$(11) In the same way, by performing discrete Fourier transform with respect to |${A_q}$| for the function |${\xi _{{k_{{I_m}}}{k_{{X_n}}}{k_{{O_p}}}}}$|⁠, we obtained its expansion formula as: $$\begin{eqnarray}{\xi _{{k_{{I_m}}}{k_{{X_n}}}{k_{{O_p}}}}}({A_q}) &&= \sum\limits_{{k_{{A_q}}} = 0}^{{K_{{A_q}}}} {{\alpha _{{k_{{I_m}}}{k_{{X_n}}}{k_{{O_p}}}{k_{{A_q}}}}}}\nonumber\\ &&\centerdot \exp \big[ {{\rm{i}}2\pi {k_{{A_q}}}{A_q}/({K_{{A_q}}}{D_{{A_q}}})} \big] ,\end{eqnarray}$$(12) where the discrete Fourier coefficient is: $$\begin{eqnarray}{\alpha _{{k_{{I_m}}}{k_{{X_n}}}{k_{{O_p}}}{k_{{A_q}}}}}\!\! &&= \frac{1}{{{K_{{A_q}}}{D_{{A_q}}}}}\sum\limits_{{k_{{A_q}}} = 0}^{{K_{{A_q}}}} {{\xi _{{k_{{I_m}}}{k_{{X_n}}}{k_{{O_p}}}{k_{{A_q}}}}}}\nonumber\\ &&\centerdot \exp \big[ {{\rm{ - i}}2\pi {k_{{A_q}}}{A_q}/({K_{{A_q}}}{D_{{A_q}}})} \big] .\end{eqnarray}$$(13) 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. Substituting equation (7) into equation (9), we obtain: $$\begin{eqnarray}{\chi _{{k_{{I_m}}}{k_{{X_n}}}}} \!\!&&= \frac{1}{{{K_{{X_n}}}{D_{{X_n}}}\centerdot {K_{{I_m}}}{D_{{I_m}}}}}\sum\limits_{{k_{{X_n}}} = 0}^{{K_{{X_n}}}} {\sum\limits_{{k_{{I_m}}} = 0}^{{K_{{I_m}}}} {F({I_m},{X_n},{O_p},{A_q})}}\nonumber\\ &&\quad\centerdot \exp \left[ { - {\rm{i}}2\pi \left(\frac{{{k_{{I_m}}}{I_m}}}{{{K_{{I_m}}}{D_{{I_m}}}}} + \frac{{{k_{{X_n}}}{X_n}}}{{{K_{{X_n}}}{D_{{X_n}}}}}\right)} \right] .\end{eqnarray}$$(14) Similarly, substituting equation (14) into equations (11) and (13) subsequently, we obtain: $$\begin{eqnarray}{\alpha _{{k_{{I_m}}}{k_{{X_n}}}{k_{{O_p}}}{k_{{A_q}}}}} &=& \frac{1}{{{K_{{I_m}}}{D_{{I_m}}}\centerdot {K_{{X_n}}}{D_{{X_n}}}\centerdot {K_{{O_p}}}{D_{{O_p}}}\centerdot {K_{{A_q}}}{D_{{A_q}}}}}\nonumber\\ &&\times \,\sum\limits_{{k_{{A_q}}} = 0}^{{K_{{A_q}}}} {\sum\limits_{{k_{{O_p}}} = 0}^{{K_{{O_p}}}} {\sum\limits_{{k_{{X_n}}} = 0}^{{K_{{X_n}}}} {\sum\limits_{{k_{{I_m}}} = 0}^{{K_{{I_m}}}} {F({I_m},{X_n},{O_p},{A_q})}}}}\nonumber\\ &&\centerdot \exp \left[ { - {\rm{i}}\Phi } \right]\! ,\end{eqnarray}$$(15) where, |${\alpha _{{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: $$\begin{eqnarray}F({I_m},{X_n},{O_p},{A_q}) &=& \sum\limits_{{k_{{A_q}}} = 0}^{{K_{{A_q}}}} {\sum\limits_{{k_{{O_p}}} = 0}^{{K_{{O_p}}}} {\sum\limits_{{k_{{X_n}}} = 0}^{{K_{{X_n}}}} {\sum\limits_{{k_{{I_m}}} = 0}^{{K_{{I_m}}}} {{\alpha _{{k_{{I_m}}}{k_{{X_n}}}{k_{{O_p}}}{k_{{A_q}}}}}}}}}\nonumber\\ &&\times \exp \left[ {{\rm{i}}\Phi } \right]\!,\end{eqnarray}$$(16) where $$\begin{eqnarray}\Phi &=& \frac{{2\pi }}{{{K_{{I_m}}}{D_{{I_m}}}}}{k_{{I_m}}}{I_m} + \frac{{2\pi }}{{{K_{{X_n}}}{D_{{X_n}}}}}{k_{{I_m}}}{X_n} + \frac{{2\pi }}{{{K_{{O_p}}}{D_{{O_p}}}}}{k_{{O_p}}}{O_p} \nonumber\\ &&\,+\, \frac{{2\pi }}{{{K_{{A_q}}}{D_{{A_q}}}}}{k_{{A_q}}}{A_q}.\end{eqnarray}$$(17) 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. 4. 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 pre-processed data are extracted into the OVT domain, where we can perform 5D data regularisation, denoising, migration imaging and other processing. Figure 3. Open in new tabDownload slide 3D seismic data processing technology before OVT domain. Figure 3. Open in new tabDownload slide 3D seismic data processing technology before OVT domain. 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. Figure 4c displays the regularised gather obtained via the OMP algorithm in the OVT domain; Since the processing is mainly for the target layer, the time window size in the 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. Figure 4. Open in new tabDownload slide Comparison of gathers before and after OVT domain regularisation. (a) Before regularisation. (b) After conventional 3D regularisation. (c) After 5D regularisation in the OVT domain. Figure 4. Open in new tabDownload slide Comparison of gathers before and after OVT domain regularisation. (a) Before regularisation. (b) After conventional 3D regularisation. (c) After 5D regularisation in the OVT domain. Figure 5. Open in new tabDownload slide Comparison of stacked profiles before and after OVT domain regularisation. (a) Before regularisation. (b) After 5D regularisation. Figure 5. Open in new tabDownload slide Comparison of stacked profiles before and after OVT domain regularisation. (a) Before regularisation. (b) After 5D regularisation. Figure 6. Open in new tabDownload slide Comparison of amplitude slices before and after OVT domain regularisation. (a) Before regularisation. (b) After 5D regularisation. Figure 6. Open in new tabDownload slide Comparison of amplitude slices before and after OVT domain regularisation. (a) Before regularisation. (b) After 5D regularisation. 5. 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: 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. The OMP data interpolation algorithm uses the Fourier transform as the basis function. This method is an iterative algorithm that creates frequency slices in the frequency domain, and it has an anti-aliasing capability. 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. Acknowledgements Thanks to WesternGeco Geophysical Company and Daqing Oilfield Co., Ltd for providing technical support. Funding The support of both the National Natural Science Foundation of China Youth Fund Project (grant no. 41804138) and the China Geological Survey Project ‘Strategic Survey of Shale Oil in the Northern Songliao Basin’ (grant no. DD20190114) is appreciated. Conflict of interest statement. None declared. References Ali Ö. , Massimiliano V., Kemal Ö., Dirk-Jan van M., Kurt E., 2010 . Crossline wavefield reconstruction from multicomponent streamer data: part 2–joint interpolation and 3D up/down separation by generalized matching pursuit , Geophysics , 75 , WB69 – WB85 . Google Scholar OpenURL Placeholder Text WorldCat Brian B. , Tim S., 2007 . Multi-azimuth and wide-azimuth seismic: shallow to deep water, exploration to production , The Leading Edge , 26 , 450 – 458 . Google Scholar OpenURL Placeholder Text WorldCat Carozzi F. , Sacchi M.D., 2019 . Robust tensor-completion algorithm for 5D seismic-data reconstruction , Geophysics , 84 , V97 – V109 . Google Scholar Crossref Search ADS WorldCat Daniel T. , 2009 . Five-dimensional interpolation: recovering from acquisition constraints , Geophysics , 74 , V123 – V132 . Google Scholar OpenURL Placeholder Text WorldCat Fernanda C. , Mauricio D. S., 2019 . Robust tensor-completion algorithm for 5D seismic-data reconstruction , Geophysics , 84 , V97 – V109 . Google Scholar OpenURL Placeholder Text WorldCat Gao J.J. , Chen X., Sacchi M.D., 2017 . Five-dimensional seismic reconstruction using parallel square matrix factorization , IEEE Transactions on Geoscience and Remote Sensing , 55 , 2124 – 2135 . Google Scholar Crossref Search ADS WorldCat Gao J.J. , Stanton A., Sacchi M.D., 2015 . Parallel matrix factorization algorithm and its application to 5D seismic reconstruction and denoising , Geophysics , 80 , V173 – V187 . Google Scholar Crossref Search ADS WorldCat Gilles H. , Lloyd F., Felix J.H., 2010 . Nonequispaced curvelet transform for seismic data reconstruction: a sparsity-promoting approach , Geophysics , 75 , WB203 – WB210 . Google Scholar OpenURL Placeholder Text WorldCat Heloise B.L. , 2016 . Image-log calibration of fracture-azimuth and fracture-density attributes from OVT prestack depth-migrated data , The Leading Edge , 35 , 164 – 169 . Google Scholar OpenURL Placeholder Text WorldCat Hennenfent G.L. , Fenelon L., Herrmann F.J., 2010 . Nonequispaced curvelet transform for seismic data reconstruction: a sparsity-promoting approach , Geophysics , 75 , WB203 – WB210 . Google Scholar Crossref Search ADS WorldCat Jin K.C. , Mauricio S., Jianjun G., 2019 . Computational efficient multidimensional singular spectrum analysis for prestack seismic data reconstruction , Geophysics , 84 , V111 – V119 . Google Scholar OpenURL Placeholder Text WorldCat Kreimer N. , Stanton A., Sacchi M.D., 2013 . Tensor completion based on nuclear norm minimization for 5D seismic data reconstruction , Geophysics , 78 , V273 – V284 . Google Scholar Crossref Search ADS WorldCat Liu B. , Sacchi M.D., 2004 . Minimum weighted norm interpolation of seismic records , Geophysics , 69 , 1560 – 1568 . Google Scholar Crossref Search ADS WorldCat Massimiliano V. , Ali Ö., Kemal Ö., Kurt E., 2010 . Crossline wavefield reconstruction from multicomponent streamer data: Part 1–multichannel interpolation by matching pursuit (MIMAP) using pressure and its crossline gradient , Geophysics , 75 , WB53 – WB67 . Google Scholar OpenURL Placeholder Text WorldCat Naghizadeh M. , Sacchi M.D., 2007 . Multistep autoregressive reconstruction of seismic records , Geophysics , 72 , V111 – V118 . Google Scholar Crossref Search ADS WorldCat Ronen J. , 1987 . Wave equation trace interpolation , Geophysics , 52 , 973 – 984 . Google Scholar Crossref Search ADS WorldCat Satinder C.H. , Kurt J.M., 2013 . Preconditioning seismic data with 5D interpolation for computing geometric attributes , The Leading Edge , 32 , 1456 – 1460 . Google Scholar OpenURL Placeholder Text WorldCat Spitz S. , 1991 . Seismic trace interpolation in the f-x domain , Geophysics , 56 , 85 – 794 . Google Scholar Crossref Search ADS WorldCat Stephen K.C. , 2014 . Multidimensional interpolation using a model-constrained minimum weighted norm interpolation , Geophysics , 79 , V191 – V199 . Google Scholar OpenURL Placeholder Text WorldCat Trad D. , 2009 . Five-dimensional interpolation: recovering from acquisition constraints , Geophysics , 74 , V123 – V132 . Google Scholar Crossref Search ADS WorldCat Trad D. , Ulrych T.J., Sacchi M.D., 2002 . Accurate interpolation with high-resolution time-variant radon transforms , Geophysics , 67 , 644 – 656 . Google Scholar Crossref Search ADS WorldCat Vermeer G.J.O. , 1998 . Creating image gathers in the absence of proper common-offset gather exploration , Geophysics , 29 , 636 – 642 . Google Scholar OpenURL Placeholder Text WorldCat Wang Y.H. , 2002 . Seismic trace interpolation in the f-x-y domain , Geophysics , 67 , 1232 – 1239 . Google Scholar Crossref Search ADS WorldCat Wang Y.H. , 2003 . Sparseness-constrained least-squares inversion: application to seismic wave reconstruction , Geophysics , 68 , 1633 – 1638 . Google Scholar Crossref Search ADS WorldCat Williams M. , Jenner E., 2002 . Interpreting seismic data in the presence of azimuthal anisotropy or azimuthal anisotropy in the presence of the seismic interpretation , The Leading Edge , 21 , 771 – 774 . Google Scholar Crossref Search ADS WorldCat Xu S. , Zhang Y., Lambare G., 2010 . Anti-leakage Fourier transform for seismic data regularization in higher dimensions , Geophysics , 75 , WB113 – WB120 . Google Scholar Crossref Search ADS WorldCat Xu Y.R. , Hao W.Y., Su Z., 2015 . Parallel matrix factorization for low-rank tensor completion , Inverse Problems and Imaging , 9 , 601 – 624 . Google Scholar Crossref Search ADS WorldCat Yen J.L. , 2003 . On nonuniform sampling of bandlimited signals , IRE Transactions on Circuit Theory , 3 , 251 – 257 . Google Scholar Crossref Search ADS WorldCat © The Author(s) 2021. Published by Oxford University Press on behalf of the Sinopec Geophysical Research Institute. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.
TI - Five-dimensional interpolation in the OVT domain for wide-azimuth seismic data
JF - Journal of Geophysics and Engineering
DO - 10.1093/jge/gxab032
DA - 2021-08-09
UR - https://www.deepdyve.com/lp/oxford-university-press/five-dimensional-interpolation-in-the-ovt-domain-for-wide-azimuth-ln2TRy2mIQ
SP - 529
EP - 538
VL - 18
IS - 4
DP - DeepDyve
ER -