From downward continuation to time stepping

We would like to propagate the recorded wavefield backward in time instead of downward into the Earth, but we also would like to preserve the computational advantages of propagating the recorded wavefield in midpoint-offset coordinates. The advantage of the midpoint-offset coordinates derives from the focusing of the reflected wavefield towards zero offset as it approaches the reflector. The wavefield focuses towards zero-offset during downward continuation because we are essentially datuming the whole data set to an increasingly deeper level in the Earth. It is thus reasonable to start our derivation from the double square root (DSR) equation, that is the main tool for datuming prestack data. As we will see later, this choice of a starting point limits the usefulness of the final result.

The DSR equation in the frequency-wavenumber domain is  

 292#292 (119)

where 2#2 is the temporal frequency, k_xs and k_xg are respectively the wavenumber associated to the source and receiver locations, and 293#293 and 294#294are the slowness at the source and receiver locations. We first start by rewriting the DSR in terms of midpoint x_m and half offset x_h as  

 295#295 (120)

where k_xm and k_xh are respectively the wavenumber associated to the midpoint x_m and the half-offset x_h.

Then, to obtain a time marching equation, we first square equation () twice and rearrange the terms into:  

 296#296 (121)

where

297#297 (122)

298#298 (123)

Equation () is a second order equation in 299#299.It has another solution in addition to the desired one. It can be greatly simplified by assuming 300#300.Then equation () can be rewritten as  

 301#301 (124)

This is the basic equation solved for the numerical examples shown in this paper. Notice that when k_xh is equal to zero, equation () degenerates to the well-known equation used for reverse-time migration of zero-offset data ().

There are few alternatives on how to solve equation () numerically. The simplest one is to use finite-differences for approximating the time derivative, and Fourier transforms for evaluating the spatial-derivative operators. Because the slowness term 302#302 is outside the parentheses in equation (), using Fourier transforms does not preclude the use of a spatially variable slowness field. Strong lateral velocity variations would cause problems because of the approximations needed to go from equation () to equation (), not because of the numerical scheme used to solve equation ().

The time marching scheme that I used can be summarized as;  

 303#303 (125)

Using a Fourier method to evaluate the spatial-derivative operators, makes it easy to handle the real limitation of equation (); that is, the presence of the vertical wavenumber 132#132at the denominator. Waves propagating horizontally have an effective infinite velocity, making a finite-difference solution unstable, no matter how small the extrapolation time step. Unfortunately, this is a major obstacle for migrating overturned events, which is one of the main goals for developing a reverse time migration in midpoint-offset coordinates. The problem exists only for finite offset data (304#304). In retrospective, the occurrence of problems for waves that overturn at finite offset should not be surprising. Equation () was derived from the DSR that cannot model data for which the source leg overturns at different depth than the receiver leg.

For non-overturning events the problem can be sidestepped. The spatial wavenumbers are related to the reflector geological dip angle 66#66and the aperture angle 121#121by the relationship

305#305 (126)

By simple trigonometry is also possible to show that for non-overturned events  

 306#306 (127)

In the Fourier domain it is straightforward to include condition () in the time-marching algorithm and thus to avoid instability without suppressing reflected energy.