Exploding-reflector migration by downward continuation

Recursive downward continuation migration algorithms start with an operator, R, that marches the wavefield down into the earth,

$$q(x,y,z + \Delta z,\omega) = R(z,\omega) \; q(x,y,z,\omega),$$ (61)

with initial conditions determined by the reflection data recorded on the surface,

$$q(x,y,z=0,\omega)=d(x,y,\omega).$$ (62)

Taking $\Delta z=1$, these equations can be rewritten in matrix form as

$$\begin{eqnarray} {\bf q}_{z + 1,\omega} & = & {\bf R}_{z,\omega} \; {\bf q}_{z,\omega} \\ {\bf q}_{1,\omega} & = &{\bf d}_{\omega}, \end{eqnarray}$$ (63) (64)

where ${\bf d}_{\omega}$ is a single frequency component of data recorded at the surface, ${\bf q}_{z,\omega}$ is a vector containing a single frequency component of the wavefield at depth-step, z, and ${\bf R}_{z,\omega}$ is the corresponding extrapolation operator.

If the downward-continuation operator, ${\bf R}_{z,\omega}$, contains purely a phase-shift, then its adjoint, ${\bf R}_{z,\omega}'$ will fully describe the inverse process of upward-continuation. However, for the amplitudes to be treated accurately, ${\bf R}_{z,\omega}$ must respect the physics of wave propagation. Stolt and Benson (1986) show that v(z) extrapolators based on WKBJ Green's functions contain an amplitude term as well as a phase term, and for v(x,y,z) earth models this effect is even more pronounced. So while kinematically-correct extrapolators are pseudo-unitary, true-amplitude depth extrapolators are not. In this introductory section, I follow conventional seismic processing methodology, and treat the depth extrapolator as a unitary operator; however, in section , I discuss how to model amplitudes correctly.

The recursion in equations () and () can be rewritten as:  

$$\left( \matrix{ 1 & 0 & 0 &...& 0 \cr -{\bf R}_{1,\omega} & ... ...atrix{ {\bf d}_{\omega} \cr 0 \cr 0 \cr ...\cr 0 \cr } \right),$$ (65)

or more simply,  

$${\bf D}_{\omega} \; {\bf q}_{\omega} = {\bf Z} \; {\bf d}_{\omega},$$ (66)

where ${\bf q}_{\omega}$ now contains ${\bf q}_{z,\omega}$ for all depth-steps, ${\bf Z}$ is a zero-padding operator, and ${\bf D}_{\omega}$ is the extrapolation matrix in equation () that can be inverted rapidly by recursion.

Equation () encapsulates the idea that we can reconstruct the wavefield at every depth-step in the earth from the wavefield at the surface by inverting matrix, ${\bf D}_{\omega}$.

As a final step, to produce a migrated image we need to invoke an imaging condition. For the case of exploding-reflector (zero-offset) migration [e.g. Claerbout (1995)], we need to extract the image, $\hat{\bf m}$, corresponding to t=0. In the temporal frequency domain, we can do this by summing over frequency. This stacking process is described by the matrix equation,

$$\left( \matrix{ \hat{\bf m}_{1} \cr \hat{\bf m}_{2} \cr \hat... ...\omega_2} \cr ...\cr {\bf q}_{\omega_{N_\omega}} \cr } \right),$$

or again more simply,

$$\hat{\bf m} = \left( \matrix{{\bf I}_{N_{xy}N_{\omega}} &{\... ...\omega}}} \right) \; {\bf q} = {\bf \Sigma}_{\omega}\; {\bf q},$$ (67)

where ${\bf m}$ is the reflectivity model, ${\bf u}$ is the total wavefield, ${\bf I}_N$ is the identity matrix of rank N, and ${\bf \Sigma}_{\omega}$ is the (exploding-reflector) imaging operator that sums over frequencies.

The process of imaging by exploding-reflector migration can then be summarized as the chain of composite operators:

$$\hat{\bf m}= {\bf \Sigma}_{\omega} \; {\bf D}^{-1} \; {\bf Z}_{N_\omega} \; {\bf d},$$ (68)

where

$${\bf D} = \left(\matrix{ {\bf D}_{\omega_1} & 0 &... & 0 \cr... ... d}_{\omega_2} \cr ... \cr {\bf d}_{\omega_{N_\omega}}}\right).$$