Numeric implementation of wave-equation migration velocity analysis operators

Survey-sinking migration and velocity analysis

Wavefield reconstruction for multi-offset migration based on the one-way wave-equation under the survey-sinking framework (Claerbout, 1985) is implemented similarly to the zero-offset case by recursive phase-shift of prestack wavefields

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

followed by extraction of the image at time $t=0$. Here, ${\bf m}$ and ${\bf h}$ stand for midpoint and half-offset coordinates, respectively, defined according to the relations

$${\bf m} = \frac{{\bf r}+{\bf s}}{2}$$ (13)

$${\bf h} = \frac{{\bf r}-{\bf s}}{2} \;,$$ (14)

where ${\bf s}$ and ${\bf r}$ are coordinates of sources and receivers on the acquisition surface. In equation 12, ${u}_z\left ({\bf m},{\bf h}\right)$ represents the acoustic wavefield for a given frequency $\omega$ at all midpoint positions ${\bf m}$ and half-offsets ${\bf h}$ at depth $z$, and ${u}_{z+\Delta z}\left ({\bf m},{\bf h}\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 double square-root (DSR) equation

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

$$+ \sqrt{ \left [{ {\omega s} \left ({\bf m}+{\bf h}\right)} \right]^2 - \left\vert {\frac{{{\bf k}_{\bf m}}+{{\bf k}_{\bf h}}}{2}} \right\vert^2} \;.$$ (15)

The image is obtained from this extrapolated wavefield by selection of time $t=0$, which is usually implemented as summation over frequencies:

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

Similarly to the derivation done in the zero-offset case, we can assume the separation of the extrapolation slowness $s \left ({\bf m}\right)$ into a background component $s_0 \left ({\bf m}\right)$ and an unknown perturbation component $\Delta s \left ({\bf m}\right)$. Then we can construct a wavefield perturbation ${\Delta {u}}\left ({\bf m},{\bf h}\right)$ at depth $z$ and frequency $\omega$ related linearly to the slowness perturbation $\Delta s \left ({\bf m}\right)$. Linearizing the depth wavenumber given by the DSR equation 15 relative to the background slowness $s_0 \left ({\bf m}\right)$, we obtain

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

where the depth wavenumber in the background medium is

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

$$+ \sqrt{ \left [{ {\omega s} _0\left ({\bf m}+{\bf h}\right)} \righ... ...- \left\vert {\frac{{{\bf k}_{\bf m}}+{{\bf k}_{\bf h}}}{2}} \right\vert^2} \;.$$ (18)

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

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

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},{\bf h}\right)= \mathcal{E}^{-}_{SSM}\left [{s_0}_z \left ({\bf m}\right),{u}_z\left ({\bf m},{\bf h}\right) \right]$. This notation indicates that the wavefield ${u}_{z+\Delta z}\left ({\bf m},{\bf h}\right)$ is reconstructed from the wavefield ${u}_z\left ({\bf m},{\bf h}\right)$ using the background slowness $s_0 \left ({\bf m}\right)$. This operation is repeated independently for all frequencies $\omega$. A typical implementation of survey-sinking wave-equation migration uses the following algorithm:



This algorithm is similar to the one described in the preceding section for zero-offset migration, except that the wavefield and image are parametrized by midpoint and half-offset coordinates and that the depth wavenumber used in the extrapolation operator is given by the DSR equation using the background slowness $s_0 \left ({\bf m}\right)$. Wavefield extrapolation is usually implemented in a mixed domain (space-wavenumber), as briefly summarized in Appendix A.

Similarly to the derivation of the wavefield perturbation in the zero-offset migration case, we can write the linearized wavefield perturbation for survey-sinking migration as

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

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

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

$$- i\Delta z \frac{\omega {u}\left ({\bf m},{\bf h}\right)\Delta s\l... ...}}} \right\vert}{2 {\omega s} _0\left ({\bf m}+{\bf h}\right)} \right]^2} } \;.$$ (20)

Equation 20 defines the forward scattering operator $\mathcal{F}^{-}_{SSM}\left [{u}\left ({\bf m},{\bf h}\right),s_0 \left ({\bf m}\right),\Delta s\left ({\bf m},{\bf h}\right) \right]$, producing the scattered wavefield ${\Delta {u}}\left ({\bf m},{\bf h}\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},{\bf h}\right)$. 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},{\bf h}\right)= \sum_\omega {\Delta {u}}\left ({\bf m},{\bf h},\omega \right)\;.$$ (21)

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

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

for every frequency $\omega$. Then, similarly to the procedure used in the zero-offset case, the slowness perturbation at depth $z$ 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}-{\bf h}\right) \approx + i \left. \frac{d {{k_z}}_s}{d s} \right\vert _{s_0} \Delta z\; \o... ...ine{{u}\left ({\bf m},{\bf h}\right)} {\Delta {u}}\left ({\bf m},{\bf h}\right)$$

$$\approx + i\Delta z \frac{\omega \overline{{u}\left ({\bf m},{\bf h}\right)... ...}}} \right\vert}{2 {\omega s} _0\left ({\bf m}-{\bf h}\right)} \right]^2} } \;,$$ (23)

and

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

$$\approx + i\Delta z \frac{\omega \overline{{u}\left ({\bf m},{\bf h}\right)... ...}}} \right\vert}{2 {\omega s} _0\left ({\bf m}+{\bf h}\right)} \right]^2} } \;.$$ (24)

Equations 23-24 define the adjoint scattering operator $\mathcal{A}^{+}_{SSM}\left [{u}\left ({\bf m},{\bf h}\right),s_0 \left ({\bf m}\right),{\Delta {u}}\left ({\bf m},{\bf h}\right) \right]$, producing the slowness perturbation $\Delta s \left ({\bf m}\right)$ from the scattered wavefield ${\Delta {u}}\left ({\bf m},{\bf h}\right)$, based on the background slowness $s_0 \left ({\bf m}\right)$ and background wavefield ${u}\left ({\bf m},{\bf h}\right)$. A typical implementation of survey-sinking forward and adjoint scattering follows the algorithms:







These algorithms are similar to the ones described in the preceding section for zero-offset migration, except that the wavefield and image are parametrized by midpoint and half-offset coordinates. Furthermore, the two square-roots corresponding to the DSR equation update the slowness model separately, thus characterizing the source and receiver propagation paths to the image positions. Both forward and adjoint scattering algorithms require the background wavefield, ${u}\left ({\bf m},{\bf h}\right)$, to be precomputed at all depth levels. 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