The imaging condition used in reverse-time migration requires that the source wavefield (computed via a forward recursion) and the receiver wavefield (computed via a backwards recursion) must be made available at the same time in an implementation of the algorithm. Several strategies to organize the calculation can be employed, differing in balance between memory and computation. This paper describes and compares these different approaches, and argues that strategies favoring computational complexity over memory (to the point where disk i/o can be avoided) are attractive for 3D prestack migrations. An example of 3D reverse-time migration applied to wide-azimuth data from the Gulf of Mexico is presented to support the claim.

Reverse-Time Migration (RTM) was introduced in the late 1970''s (Hemon, 1978) but despite showing promising imaging capabilities (Baysal et al., 1983; Whitmore, 1983; McMechan, 1983; Loewenthal and Mufti, 1983), it was not used in practice due to its stringent requirements, both in terms of computation and memory. Until recently, RTM was therefore largely confined to 2D and/or post-stack imaging but computer technology has now reached the point where 3D prestack RTM isfeasible (Yoon et al., 2003, 2004; Bednar and Bednar, 2006; Farmer et al., 2006; Guitton et al., 2007; Jones et al., 2007). The core of the algorithm is the crosscorrelation of two wavefields at the same time level, one computed by stepping forward in time, the other computed by stepping backwards in time. The forward recursion is usually carried out first, in which case the entire time history must be made available during the backwards recursion in order to compute the imaging condition. Several options can be used to arrange the calculation, differing in the amount of memory and computation requirements. The next section describes the realization of a particular approach to RTM based on leapfrog time stepping for a scalar field (pressure). The following section compares the various strategies which can be used to organize the computations. Finally, an example of 3D RTM applied to wide-azimuth data from the Gulf of Mexico is presented.

The forward problem involves marching this scheme for n = 0,1, . . . ,N -1. The simulation of synthetic seismic data dn is related to pn by a sampling operator Sn which extracts time samples of the wavefield at receiver positions, at the time sample rate of the output data traces. The implementation used for the examples shown below does not assume any relation between the computation grid and the acquisition geometry, or between the simulation time step and the sample rate of the seismic traces. It uses bilinear interpolation for the spatial variables and cubic interpolation in time. Similarly, the source field is adjoint-interpolated onto the computational grid. The Laplacian is approximated using an eighthorder centered difference scheme with optimized coefficients (Ye and Chu, 2005; Etgen, 2007) and Perfectly Matched Layers (PML) absorbing boundary conditions (Cohen, 2001) are used to simulate unbounded domain wave propagation.