Minimum weighted norm wavefield reconstruction for AVA imaging

Seismic wavefield reconstruction is posed as an inversion problem where, from inadequate and incomplete data, we attempt to recover the data we would have acquired with a denser distribution of sources and receivers. A minimum weighted norm interpolation method is proposed to interpolate prestack volumes before wave-equation amplitude versus angle imaging. Synthetic and real data were used to investigate the effectiveness of our wavefield reconstruction scheme when preconditioning seismic data for wave-equation amplitude versus angle imaging. I N T R O D U C T I O N Robust inversion of prestack seismic data is an important step towards the estimation of rock properties and fluid indicators. Recent developments have demonstrated that waveequation amplitude versus angle (AVA) imaging is an emerging and powerful methodology for the accurate estimation of AVA gathers (Mosher, Foster and Hassanzadeh 1997; Prucha, Biondi and Symes 1999). It is well known that waveequation imaging methods require regularly acquired wavefields. In order to maximize the benefits of wave-equation AVA imaging when working with irregularly sampled data, Kuehl and Sacchi (2003) and Wang, Kuehl and Sacchi (2003) proposed least-squares migration methods to account for missing observations in the prestack volume. A more economical alternative entails the reconstruction of the acquired wavefield before wave-equation AVA imaging. In this case, the interpolation/resampling problem can be posed as an inversion problem where, from inadequate and incomplete data, an attempt is made to recover a properly sampled version of the original seismic wavefield. The problem, however, is underdetermined and, as is well known, the solution is not unique. Missing observations lie in the null space of the sampling operator and therefore a regularization strategy is required to retrieve a unique and stable data reconstruction. Regularization methods are used not only to retrieve stable ∗E-mail: [email protected] and unique data reconstructions but also to impose desirable features on the spectrum of the reconstructed wavefield. For example, minimum norm spectral regularization can be used when we assume that seismic data are band-limited in the spatial wavenumber domain (Cary 1997; Duijndam, Schonewille and Hindriks 1999; Hindriks and Duijndam 2000). Similarly, a regularization derived using the Cauchy criterion can be used to obtain a high-resolution (sparse) discrete Fourier transform that can be used to perform the synthesis of the data at new spatial positions (Sacchi, Ulrych and Walker 1998; Zwartjes and Duijndam 2000). In the minimum weighted norm interpolation (MWNI) method (Liu and Sacchi 2001, 2003), we have used a spectral weighted norm regularization term that incorporates a priori knowledge of the energy distribution of the signal to be interpolated. The technique can be used to interpolate large portions of data simultaneously, along any number of spatial dimensions. It is important to stress that the MWNI algorithm is quite efficient; the computational cost of the interpolation relies on fast Fourier transforms (FFTs) in conjunction with a preconditioned conjugate-gradient scheme to accelerate convergence. We present examples that illustrate the application of the MWNI algorithm to 2D/3D prestack seismic data regularization. We also test the effectiveness of our interpolation strategy at the time of reconstructing data before 2D/3D wave-equation AVA imaging (Mosher et al. 1997; Prucha et al. 1999). It is important to stress that rather than performing the classical comparison of data before and after interpolation, we have preferred to compare migrated images C © 2005 European Association of Geoscientists & Engineers 787 788 M.D. Sacchi and B. Liu and to extract AVA curves obtained from common-image gathers before and after interpolation. Interpolation strategies based on prediction error filters have been shown to yield excellent results in problems where there is enough information to predict the missing data properly (Claerbout 1991; Spitz 1991). However, prediction error filters are unable to handle situations where there are large gaps in the data. In this case, methods that exploit data multidimensionality should be preferred. The type of algorithm proposed in this paper can handle large segments of missing information by using spatial data from as many coordinates as possible. In the first part of this paper, we introduce the interpolation problem as an inverse problem where a wavenumber-domain norm is used to regularize the inverse problem. Then we briefly review the problem of estimating AVA gathers using waveequation migration methods. The final part of the paper is devoted to numerical and real data examples. We first examine the 2D problem with a simple synthetic model and with the Marmousi model. Then the 3D reconstruction problem is studied with a simple 3D synthetic data set and a field data set from the Western Canadian Sedimentary Basin. These examples are used to highlight the strengths of our algorithm when preconditioning data for wave-equation AVA imaging. P R O B L E M F O R M U L AT I O N The multidimensional or N-dimensional (ND) interpolation is carried out along the vector of spatial coordinates u for each temporal frequency ω. In other words, we denote the seismic data at one monochromatic temporal frequency component as D(u, ω), where, for instance, u = [xs, xr] denotes a 2D data volume defined in terms of source and receiver positions; similarly, u = [mx, my, hx] defines a constant azimuth 3D data volume in terms of the two midpoint positions mx, my and the in-line offset hx. It should be remembered that ND refers to the number of spatial dimensions of the reconstruction problem. The methodology presented in this paper can be easily adapted to handle the following cases: N = 1 : interpolation of one single gather; N = 2 : prestack 2D interpolation or 3D post-stack interpolation; N = 3 : interpolation of 3D prestack common-azimuth volumes; N = 4 : full spatial interpolation (multi-azimuth interpolation). The multidimensional wavefield D(u, ω) at a monochromatic frequency component ω can be organized via lexicographic ordering in a length-M vector x = [x1, x2, . . . xM]. The vector x defines the desired observations on the regular grid. The actual observations (recorded data) are given by the elements of the vector y = [xn(1), xn(2), xn(3), . . . xn(N)], where the set N = {n(1), n(2), n(3), . . . , n(N)} is used to indicate the position of the known samples or observations in the regular grid. It is clear that with this type of sampling, one is assigning recorded traces to the nearest spatial position in the desired grid (binning). We now define the sampling matrix T with elements Ti,j = δn(i),j, where δ denotes the Kronecker operator. It is quite simple to show that the complete data and the observations are connected by the following linear system of equations: