Wavefield scattering

Imaging by wavefield extrapolation (WE) is based on recursive continuation of the wavefields ($\u$) from a given depth level to the next, by means of an extrapolation operator (${\bf E}$):  

$$\u_{z+\Delta z}= {\bf E}_z\left[\u_{z } \right] \;.$$ (1)

Here and hereafter, we use the following notation conventions: ${\bf A}\left[x \right]$ stands for the linear operator ${\bf A}$ applied to x, and $f\left(x \right)$ stands for function f of argument x.

At any depth z, the wavefield ($\widetilde{\u}$), extrapolated through the background medium characterized by the background velocity ($\tilde{s}$), interacts with medium perturbations ($\Delta s$) and creates wavefield perturbations ($\Delta \mathcal V$):

$$\Delta \mathcal V_{z }= {\bf S}_z\left({\widetilde{\u}}_{z }\right)\left[{\Delta s}_{z}\right] \;.$$ (2)

${\bf S}$ is a scattering operator relating slowness perturbations to wavefield perturbations. The total wavefield perturbation at depth $z+\Delta z$ is the sum of the perturbation accumulated up to depth z from all depths above ($\Delta \u_{z }$), plus the scattered wavefield from depth ($\Delta \mathcal V_{z }$)extrapolated one depth step ($\Delta z$):  

$$\Delta \u_{z+\Delta z}= {\bf E}_z\left[\Delta \u_{z } \right... ...detilde{\u}}_{z }\right)\left[{\Delta s}_{z}\right] \right] \;.$$ (3)

We can use the recursive equation (3) to compute a wavefield perturbation, given a precomputed background wavefield and a slowness perturbation. In a more compact notation, we can write equation (3) as follows:  

$$\Delta \u= \left({\bf 1}- {\bf E}\right)^{-1} {\bf E}{\bf S}\Delta s\;,$$ (4)

where $\Delta \u$ and $\Delta s$ stand respectively for the wavefield and slowness perturbations at all depth levels, and ${\bf E}$, ${\bf S}$ and ${\bf 1}$ are respectively the wavefield extrapolation operator, the scattering operator and the identity operator. In our current implementation, ${\bf S}$ refers to a first-order Born scattering operator.

From the wavefield perturbation ($\Delta \u$), we can compute an image perturbation ($\Delta \r$) by applying an imaging condition, $\Delta \r= {\bf I}\Delta \u$. For example, the imaging operator, (${\bf I}$) can be a simple summation over frequencies. If we accumulate all scattering, extrapolation and imaging into a single operator ${\bf L}= {\bf I}\left({\bf 1}- {\bf E}\right)^{-1} {\bf E}{\bf S}$,we can write a simple linear expression relating an image perturbation ($\Delta \r$) to a slowness perturbation ($\Delta s$):  

$$\Delta \r= {\bf L}\Delta s\;.$$ (5)

For wave-equation migration velocity analysis, we use equation (5) to estimate a perturbation of the slowness model from a perturbation of the migrated image by minimizing the objective function  

$$J \left(\Delta s\right)= \left\vert\left\vert{\bf W}\left(\Delta \r- {\bf L}\Delta s\right)\right\vert\right\vert^2 \;,$$ (6)

where ${\bf W}$ is a weighting operator related to the inverse of the data covariance, indicating the reliability of the data residuals. Since, in most practical cases, the inversion problem is not well conditioned, we need to add constraints on the slowness model via a regularization operator. In these situations, we use the modified objective function  

$$J \left(\Delta s\right) = \left\vert\left\vert{\bf W}\left(\... ...^2 \left\vert\left\vert\AA \Delta s\right\vert\right\vert^2 \;.$$ (7)

$\AA$ can be a regularization operator which penalizes rough features of the model, and $\epsilon$ is a scalar parameter which balances the relative importance of the data residual, ${\bf W}\left(\Delta \r- {\bf L}\Delta s\right)$, and the model residual, $\left(\AA \Delta s\right)$.

An essential element of our velocity analysis method is the image perturbation, $\Delta \r$. For the purposes of the optimization problem in equation (7), this is object is known and has to be precomputed, together with the background wavefield used by the operator ${\bf L}$.In the next section, we discuss how $\Delta \r$is estimated in practice.