Numeric implementation of wave-equation migration velocity analysis operators

Zero-offset migration and velocity analysis

Wavefield reconstruction for zero-offset migration based on the one-way wave-equation is performed by recursive phase-shift in depth starting with data recorded on the surface as boundary conditions. In this configuration, the imaging condition extracts the image as time $t=0$ from the reconstructed wavefield at every location in space. Thus, the surface data need to be extrapolated backward in time which is achieved by selecting the appropriate sign of the phase-shift operation (which depends on the sign convention for temporal Fourier transforms):

$${u}_{z+\Delta z} \left ({\bf m}\right)= e^{- i {k_z}\Delta z}{u}_z \left ({\bf m}\right)\;.$$ (1)

In equation 1, ${u}_z \left ({\bf m}\right)$ represents the acoustic wavefield at depth $z$ for a given frequency $\omega$ at all positions in space ${\bf m}$, and ${u}_{z+\Delta z} \left ({\bf m}\right)$ represents the same wavefield extrapolated to depth $z+\Delta z$. The phase shift operation uses the depth wavenumber ${k_z}$ which is defined by the single square-root (SSR) equation

$${k_z}= \sqrt{ \left [{2 {\omega s} \left ({\bf m}\right)} \right]^2 - \left\vert {{{\bf k}_{\bf m}}} \right\vert^2} \;,$$ (2)

where $s \left ({\bf m}\right)$ represents the spatially-variable slowness at depth level $z$. Equations 1-2 describe wavefield extrapolation using a pseudo-differential operator, which justifies our use of laterally-varying slowness $s \left ({\bf m}\right)$. As indicated earlier, the image is obtained from this extrapolated wavefield by selection of time $t=0$, which is typically implemented as summation of the extrapolated wavefield over frequencies:

$${r}_z \left ({\bf m}\right)= \sum_\omega {u}_z \left ({\bf m}, \omega \right)\;.$$ (3)

Phase-shift extrapolation using wavenumbers computed using equations 1 and 2 is not feasible in media with lateral variation. Instead, implementation of such operators is done using approximations implemented in a mixed space-wavenumber domain (Huang et al., 1999; Stoffa et al., 1990; Ristow and Ruhl, 1994). A brief summary the mixed-domain implementation of the split-step Fourier (SSF) operator is presented in Appendix A.

For velocity analysis, we assume that we can separate the total slowness $s \left ({\bf m}\right)$ into a known background component $s_0 \left ({\bf m}\right)$ and an unknown component $\Delta s \left ({\bf m}\right)$. With this convention, we can linearize the depth wavenumber ${k_z}$ relative to the background slowness $s_0$ using a truncated Taylor series expansion

$${k_z}\approx {k_z}_0 + \left. \frac{d {{k_z}}} {d s} \right\vert _{s_0} \Delta s \left ({\bf m}\right)\;,$$ (4)

where the depth wavenumber in the background medium characterized by slowness $s_0 \left ({\bf m}\right)$ is

$${k_z}_0 = \sqrt{ \left [{2 {\omega s} _0 \left ({\bf m}\righ... ...\right]^2 - \left\vert {{{\bf k}_{\bf m}}} \right\vert^2} \;.$$ (5)

Here, $s_0 \left ({\bf m}\right)$ represents the spatially-variable background slowness at depth level $z$. Using the wavenumber linearization from equation 4, we can reconstruct the acoustic wavefields in the background model using a phase-shift operation

$${u}_{z+\Delta z} \left ({\bf m}\right)= e^{- i {k_z}_0 \Delta z}{u}_z \left ({\bf m}\right)\;.$$ (6)

We can represent wavefield extrapolation using a generic solution to the one-way wave-equation using the notation ${u}_{z+\Delta z} \left ({\bf m}\right)= \mathcal{E}^{-}_{ZOM}\left [2{s_0}_z \left ({\bf m}\right),{u}_z \left ({\bf m}\right) \right]$. This notation indicates that the wavefield ${u}_{z+\Delta z} \left ({\bf m}\right)$ is reconstructed from the wavefield ${u}_z \left ({\bf m}\right)$ using the background slowness $s_0 \left ({\bf m}\right)$. This operation is repeated independently for all frequencies $\omega$. A typical implementation of zero-offset wave-equation migration uses the following algorithm:



A seismic image is produced using migration by wavefield extrapolation as follows: for each frequency, read data at all spatial locations ${\bf m}$; then, for each depth, sum the wavefield into the image at that depth level (i.e. apply the imaging condition) and then reconstruct the wavefield to the next depth level (i.e. perform wavefield extrapolation). The ``-'' sign in this algorithm indicates that extrapolation is anti-causal (backward in time), and the factor ``2'' indicates that we are imaging data recorded in two-way traveltime with an algorithm designed under the exploding reflector model. Wavefield extrapolation is usually implemented in a mixed domain (space-wavenumber), as briefly summarized in Appendix A.

The wavefield perturbation ${\Delta {u}} \left ({\bf m}\right)$ caused at depth $z+\Delta z$ by a slowness perturbation $\Delta s \left ({\bf m}\right)$ at depth $z$ is obtained by subtraction of the wavefields extrapolated from $z$ to $z+\Delta z$ in the true and background models:

$${\Delta {u}}_{z+\Delta z} \left ({\bf m}\right) = e^{- i {k_z}\Delta z}{u}_z \left ({\bf m}\right)- e^{- i {k_z}_0 \Delta z}{u}_z \left ({\bf m}\right)$$

$$= e^{- i {k_z}_0 \Delta z} \left [e^{- i \left. \frac{d {{k_z}}} {d... ...Delta s \left ({\bf m}\right)\Delta z} - 1\right]{u}_z \left ({\bf m}\right)\;.$$ (7)

Here, ${\Delta {u}} \left ({\bf m}\right)$ and ${u} \left ({\bf m}\right)$ correspond to a given depth level $z$ and frequency $\omega$. A similar relation can be applied at all depths and all frequencies.

Equation 7 establishes a non-linear relation between the wavefield perturbation ${\Delta {u}} \left ({\bf m}\right)$ and the slowness perturbation $\Delta s \left ({\bf m}\right)$. Given the complexity and cost of numeric optimization based on non-linear relations between model and wavefield parameters, it is desirable to simplify this relation to a linear relation between model and data parameters. Assuming small slowness perturbations, i.e. small phase perturbations, the exponential function $e^{\pm i \left. \frac{d {{k_z}}} {d s} \right\vert _{s_0} \Delta s \left ({\bf m}\right) \Delta z}$ can be linearized using the approximation $e^{i\phi}\approx 1+i\phi$ which is valid for small values of the phase $\phi$. Therefore the wavefield perturbation ${\Delta {u}} \left ({\bf m}\right)$ at depth $z$ can be written as

$${\Delta {u}} \left ({\bf m}\right) \approx - i \left. \frac{d {{k_z}}} {d s} \right\vert _{s_0} \Delta z\; {u} \left ({\bf m}\right)\Delta s \left ({\bf m}\right)$$

$$\approx - i\Delta z \frac{2\omega {u} \left ({\bf m}\right)\Delta s \left... ..._{\bf m}}} \right\vert}{2 {\omega s} _0 \left ({\bf m}\right)} \right]^2} } \;.$$ (8)

Equation 8 defines the zero-offset forward scattering operator $\mathcal{F}^{-}_{ZOM}\left [{u} \left ({\bf m}\right),2s_0 \left ({\bf m}\right),\Delta s \left ({\bf m}\right) \right]$, producing the scattered wavefield ${\Delta {u}} \left ({\bf m}\right)$ from the slowness perturbation $\Delta s \left ({\bf m}\right)$, based on the background slowness $s_0 \left ({\bf m}\right)$ and background wavefield ${u} \left ({\bf m}\right)$ at a given frequency $\omega$. The image perturbation at depth $z$ is obtained from the scattered wavefield using the time $t=0$ imaging condition, similar to the procedure used for imaging in the background medium:

$$\Delta {r} \left ({\bf m}\right)= \sum_\omega {\Delta {u}} \left ({\bf m}, \omega \right)\;.$$ (9)

Given an image perturbation $\Delta {r} \left ({\bf m}\right)$, we can construct the scattered wavefield by the adjoint of the imaging condition

$${\Delta {u}} \left ({\bf m}, \omega \right)= \Delta {r} \left ({\bf m}\right)\;,$$ (10)

for every frequency $\omega$. Then, the slowness perturbation at depth $z$ and frequency $\omega$ caused by a wavefield perturbation at depth $z$ under the influence of the background wavefield at the same depth $z$ can be written as

$$\Delta s \left ({\bf m}\right) \approx + i \left. \frac{d {{k_z}}} {d s} \right\vert _{s_0} \Delta z\; \overline{{u} \left ({\bf m}\right)} {\Delta {u}} \left ({\bf m}\right)$$

$$\approx + i\Delta z \frac{2\omega \overline{{u} \left ({\bf m}\right)} {\... ... k}_{\bf m}}} \right\vert}{2 {\omega s} \left ({\bf m}\right)} \right]^2} } \;.$$ (11)

Equation 11 defines the adjoint scattering operator $\mathcal{A}^{+}_{ZOM}\left [{u} \left ({\bf m}\right),2s_0 \left ({\bf m}\right),{\Delta {u}} \left ({\bf m}\right) \right]$, producing the slowness perturbation $\Delta s \left ({\bf m}\right)$ from the scattered wavefield ${\Delta {u}} \left ({\bf m}\right)$, based on the background slowness $s_0 \left ({\bf m}\right)$ and background wavefield ${u} \left ({\bf m}\right)$. A typical implementation of zero-offset forward and adjoint scattering uses the following algorithms:







The forward zero-offset wave-equation MVA operator follows a similar pattern to the implementation of the downward continuation operator: for each frequency and for each depth, read the slowness perturbation $\Delta s$ at all spatial locations ${\bf m}$, then apply the scattering operator ( ${\textbf{\sc w.s.}}$) given equation 11 to the slowness perturbation and accumulate the scattered wavefield for downward continuation; then, apply the imaging condition ( ${\textbf{\sc i.c.}}$) producing the image perturbation $\Delta {r}$ at depth $z$ and then reconstruct the scattered wavefield backward in time using the wavefield extrapolation operator ( ${\textbf{\sc w.r.}}$) to the next depth level. The adjoint zero-offset wave-equation MVA operator also follows a similar pattern to the implementation of the downward continuation operator: for each frequency and for each depth, reconstruct the scattered wavefield forward in time using the wavefield extrapolation operator ( ${\textbf{\sc w.r.}}$) to the next depth level, then apply the adjoint of the imaging condition ( ${\textbf{\sc i.c.}}$) by adding the image to the scattered wavefield and then apply the adjoint wavefield scattering operator ( ${\textbf{\sc w.s.}}$) to obtain the slowness perturbation $\Delta s$. Both forward and adjoint scattering algorithms require the background wavefield, ${u}$, to be precomputed at all depth levels, although more efficient implementations using optimal checkpointing are possible (Symes, 2007). Scattering and wavefield extrapolation are implemented in the mixed space-wavenumber domain, as briefly explained in Appendix A.

Numeric implementation of wave-equation migration velocity analysis operators