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  

$$k_z= \sqrt{\omega^2\ss^2 - k_{x_s}^2} + \sqrt{\omega^2s({{\bf g},z})^2 - k_{x_g}^2},$$ (1)

where $\omega$ is the temporal frequency, k_xs and k_xg are respectively the wavenumber associated to the source and receiver locations, and $\ss$ and $s({{\bf g},z})$are 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  

$$k_z= \sqrt{\omega^2\ss^2 - \frac{\left(k_{x_m}-k_{x_h}\right... ...^2s({{\bf g},z})^2 - \frac{\left(k_{x_m}+k_{x_h}\right)^2}{4}},$$ (2)

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 (2) twice and rearrange the terms into:  

$$\omega^4\Delta_s + 2\omega^2\left(\Delta_s k_{x_m}k_{x_h}- \... ...k_z^2\left(k_{x_m}^2 + k_{x_h}^2\right) + k_{x_m}^2k_{x_h}^2=0,$$ (3)

where

$$\Delta_s = \ss^2-s({{\bf g},z})^2$$ (4)

$$\Sigma_s = \ss^2+s({{\bf g},z})^2$$ (5)

Equation (3) is a second order equation in $\omega^2$.It has another solution in addition to the desired one. It can be greatly simplified by assuming $\Delta_s \approx 0$.Then equation (3) can be rewritten as  

$$\omega^2=\frac{1}{2\Sigma_s} \left(k_z^2 + k_{x_m}^2 + k_{x_h}^2 + \frac{k_{x_m}^2k_{x_h}^2}{k_z^2}\right).$$ (6)

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

There are few alternatives on how to solve equation (6) 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 $\Sigma_s$ is outside the parentheses in equation (6), 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 (3) to equation (6), not because of the numerical scheme used to solve equation (6).

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

$$\frac{P_{t-\Delta t}-2 P_{t} + P_{t+\Delta t}} {\Delta t ^2}... ..._h}^2 + \frac{k_{x_m}^2k_{x_h}^2}{k_z^2}\right) {\rm FFT}\;\;P.$$ (7)

Using a Fourier method to evaluate the spatial-derivative operators, makes it easy to handle the real limitation of equation (6); that is, the presence of the vertical wavenumber k_z at 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 ($k_{x_h}\neq 0$). In retrospective, the occurrence of problems for waves that overturn at finite offset should not be surprising. Equation (6) 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 $\gamma$and the aperture angle $\alpha$by the relationship

$$\frac{k_{x_m}^2k_{x_h}^2}{k_z^2}=\tan \alpha \tan \gamma,$$ (8)

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

$$\frac{k_{x_m}^2k_{x_h}^2}{k_z^2}a\leq 1.$$ (9)

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