An explicit, constrained operator is presented here for wavefield extrapolation in 3D wave equation depth migration. The migration cost and image quality benefit from its reduced number of independent coefficients within the operator, negligible numerical anisotropy, and flexibility that allows for different propagation angles and step sizes in the inline and cross-line directions. In order to further reduce the computational workload we dynamically select operator lengths and extrapolation step sizes based on the wavenumber of the wave components being migrated. The phase-shifted linear interpolation that we propose for interpolating the extrapolated wavefield is suitable for the explicit migration, and significantly improves the accuracy of the result when compared with the linear interpolation typically used in implicit migrations.


Downward extrapolation of the seismic wavefield is the key element of wave equation based migrations because it determines the image quality and computational cost. Extrapolation by finite difference methods is considered to be superior, especially in handling lateral velocity variations. It can be classified into explicit and implicit methods. Accomplished simply by convolving known wavefields with an extrapolation operator, the explicit method is easier to implement, and can be more efficient than the implicit method. In addition, the explicit method can be extended to 3D in a straightforward manner, whereas the implicit method requires inline and crossline splitting of the extrapolation process at each depth step, resulting in phase errors and artifacts.

Even with symmetrical operators (Blacquiere et al., 1989), 3D depth migrations based on explicit extrapolation may still be expensive, due to a full 2D convolution being performed at each depth. The scheme based on the McClellan transform proposed by Hale (1991) is much more efficient, but introduces numerical anisotropy at the same time. Among the many alternatives, or improvements to the scheme (Soubaras, 1992; Sollid and Arntsen, 1994; Biondi and Palacharla, 1995; Mittet, 2002), the constrained operator proposed by Mittet (2002) is attractive. It overcomes the numerical anisotropy problem and still keeps the computational cost comparable to the Hale-McClellan scheme.

Method and Results

Constrained operator

As demonstrated in Figure 1, the explicit finite-difference operator is separated into two areas during the operator design. The operator is a function of the local wavenumber at each output grid location, and is able to handle lateral velocity variations. As usual for explicit finite-difference migration operators, the operator assumes a circular symmetry. What is unique about the operator presented here is that it divides the coefficients into core and constrained areas. In the core area the coefficient layout is the same as that for a full operator. In the constrained area the coefficients in a series of rings are assigned an identical value. Figure 1 shows the layout of a constrained operator with the Δy/Δx ratio of 1.5, x-direction half-length of 12, y-direction half-length of 8, and core half-length of 6.

Fig. 01
Fig. 1. Coefficient layout for one quadrant of a constrained operator. The axes are the inline and cross-line directions for a given depth level. As the operator size increases, the maximum dip that can be migrated increases, however, the computational cost also correspondingly increases. In order to reduce computational cost, the operator design in this example imposes equal coefficient values within each circular band. This design process is pursued prior to the actual migration, and ensures that migration accuracy is still extremely high. From Mittet (2002).

The constrained operator reduces the number of independent coefficients. By factorizing out the coefficients in a ring during the convolution, the number of complex multiplications is reduced. Because the complex multiplication is expensive, using the constrained operator improves the extrapolation efficiency.

The operator half-length, M, in the y direction, is determined from the operator half-length, L, in the x direction, using the relationship M=(Δy/Δx)L. If Δy is greater than Δx, as is usually the case for marine acquisition, the total number of operator coefficients and the computational cost is reduced.

The coefficients of the constrained operator are designed using the exact extrapolator in the frequency-wavenumber domain, that is, the phase shift operator. The design strategy is to optimize the constrained operator such that its wavenumber-domain response approaches that of the exact phase shift operator. For any given local wavenumber, the optimization is performed on the operator’s real and imaginary parts, using an L8 norm for all inline and cross-line wavenumbers (kx and ky) up to the critical wavenumber

Eq. 01

determined by:

Eq. 02

where are the maximum propagation angles considered in the x and y directions.

Eq. 03

are the maximum propagation angles considered in the x and y directions.

We design the constrained operator for a sufficient set of kw and form an operator table. For a given accuracy, the operator length required varies with kw. Our extrapolation process selects the shortest operator that satisfies the accuracy required based on the wavenumber value of the wave component being migrated. This measure avoids using unnecessarily long operators, and can reduce computational cost.

Variable extrapolation step size and interpolation

To reduce the number of extrapolations we use a variable step size that is determined dynamically by the wavenumber value. To interpolate the extrapolated wavefield onto a regular output depth grid we propose a phase-shifted linear interpolation. We do not consider the linear interpolation used in implicit migrations where the interpolation is performed only on the “diffraction” term, and the “thin-lens” term is calculated exactly. In an explicit migration those two terms are not separated, and the use of linear interpolation on the “thin-lens” term will generate a large error. The phase shifts are applied to each extrapolated wavefield frequency prior to the linear interpolation to the output depth grid, and eliminate the influence of the “thin-lens” term.

Fig. 02
Fig. 2. Source wavefields obtained with phase-shifted linear and linear interpolation. As observed, the application of an appropriate phase shift correct after wavefield extrapolation (left) significantly reduces the artifacts that might arise when interpolating to the final output depth grid.

Figure 2 shows a part of the source wavefield snapshot obtained with the constrained operator and phase-shifted linear interpolation. The extrapolation interval is 25 m and the interpolation interval is 5 m. The source function used contains 4 wavelets with a flat frequency spectrum between 1 – 39 Hz. The model has a constant velocity of 2000 m/s. The wavefield obtained with the phase-shifted linear interpolation is clean, and does not have obvious amplitude distortions and phase errors. For comparison, also shown in Figure 2, is the result obtained with linear interpolation. We note that the linear interpolation generates strong noise cutting through the wavefronts, and obvious amplitude distortions and phase errors can be observed, especially around the low dip portions of the wavefront.


To illustrate the performance of the explicit 3D wave equation migration, we have applied it to the SEG/EAGE salt model dataset. Figure 3 compares Kirchhoff and explicit finite-difference migrations. In comparison to Kirchhoff migration, the image from the wave equation migration has a lower noise level, does a better job imaging the salt boundaries, and, most importantly, it produces a much more focused and continuous image.

Fig. 03
Fig. 3. Images from the SEG/EAGE model dataset, obtained using Kirchhoff and explicit wave equation migrations. The image on the left is obtained using energetic raytracing and Kirchhoff migration. The image in the middle is obtained using the constrained operator and extrapolation with a variable step size, plus the phase-shifted linear interpolation. The image on the right is the velocity model used in both cases. Note that the wave equation migration yields much better focussing, and less artifacts.

We have also applied the explicit wave equation migration to a real dataset from the southern North Sea (Figure 4). Not only does the explicit finite-difference migration yield significantly better steep dip imaging in the centre area of the example shown, but the coherency and resolution of low dip events are also much better than the Kirchhoff result. A thin, high-velocity chalk layer confronts imaging quality for deeper events in this region, as evidenced in the Kirchhoff result.

Fig. 04
Fig. 4. Explicit finite-difference vs. Kirchhoff migration; example from the southern North Sea. Note the better steep dip imaging, and the better imaging below the fast flat layer on the left. The wave equation migration is able to significantly improve steep dip imaging and general resolution, whilst being pursued with an efficient implementation.


The constrained operator is accurate and reliable for wavefield extrapolation in explicit depth migrations. The smaller number of independent coefficients offered by the operator leads to a reduced migration computational cost. Its flexibility to allow for a larger grid interval and shorter operator length in the cross-line direction makes it a useful and efficient method for migrating marine data. The measures we take, including the dynamic selection of the operator length and extrapolation step size based on the wavenumber, also contribute to the efficiency of the method. The phase-shifted linear interpolation method proposed for interpolating the extrapolated wavefield is suitable for explicit migration, and is more accurate than the linear interpolation.



The authors would like to thank Rune Mittet and John Anderson for their useful suggestions and discussions. We also thank Sandy Carroll, Boris Tsimelzon and Peter Wijnen for their assistance with processing.

This article is based on the Best Petroleum Paper which was presented by Andrew Long at the 17th ASEG Conference and Exhibition, held in Sydney during August 2004. It was first published in the October 2004 issue of the ASEG's Preview magazine and is re-published here with permission of the ASEG.


About the Author(s)


Biondi, B. and Palacharla, G., 1995, 3D depth migration by rotated McClellan filters: Geophys. Prosp., 43, 1005-1020.

Blacquiere, G., Debeye, H. W. J., Wapenaar, C. P. A., and Berkhout, A. J., 1989, 3-D table driven migration: Geophys. Prosp., 37, 925-958.

Hale, D., 1991, 3-D depth migration via McClellan transformations: Geophysics, 56, 1778-1785.

Mittet, R., 2002, Explicit 3D depth migration with a constrained operator: 72nd Annual International Meeting, SEG, Salt Lake City, Expanded Abstracts, 1148-1151.

Ren, J., Gerrard, C., McClean, J., Orlovich, M., and Long, A., 2004, An efficient 3D explicit depth migration: 17th Annual ASEG Meeting, Sydney, Extended Abstracts.

Sollid, A., and Arntsen, B., 1994, Cost effective 3D one-pass depth migration: Geophys. Prosp., 42, 755-776.

Soubaras, R., 1992, Explicit 3-D migration using equiripple polynomial expansion and Laplacian synthesis: 62nd Annual International Meeting, SEG, New Orleans, Expanded Abstracts, 905-908.


Join the Conversation

Interested in starting, or contributing to a conversation about an article or issue of the RECORDER? Join our CSEG LinkedIn Group.

Share This Article