Waveform Inversion Algorithm

Waveform inversion seeks to reconstruct the earth model parameters, e.g. velocity and density, from the waveform data. In most cases the waveform inversion algorithm will consist of minimizing a waveform misfit function by a gradient method. This is explained below using the example of inverting acoustic seismic data in the frequency domain.

1.

Misfit Function: Assume harmonic energy emanating from a source at $\bf {r}_s$ and recorded by a geophone at $\bf {r}_g$. The observed pressure spectrums are denoted by $\tilde{P}(\bf {r}_g,\bf {r}_s)_{obs}$, where the unknown source wavelet w(t) has the spectrum denoted by $\tilde{W}(\omega)$. The pressure misfit function $\epsilon$ is given by:

 

$$\epsilon 1/2 \sum_{\omega} \sum_g \sum_s \vert \Delta \tilde{ P}(\bf {r}_g,\bf {r}_s)\vert^2 ,$$ = (3.1)

where $\tilde{P}(\bf {r}_g,\bf {r}_s)$ is the predicted pressure field, the data residual $\Delta \tilde P(\bf {r}_g,\bf {r}_s)$ is defined by $\Delta \tilde P(\bf {r}_g,\bf {r}_s) = \tilde{ P} (\bf {r}_g,\bf {r}_s) - \tilde{P}(\bf {r}_g,\bf {r}_s)_{obs}$, and the summation is over source frequencies and the source+geophone coordinates. The predicted pressure is usually obtained by a finite-difference solution to the Helmholtz equations.

2.

Slowness Gradient and Frechet Derivative: Using a steepest descent method (see Chapter 1), the slowness model $s(\bf {r})^{i}$ at the i^th iterate is updated by

 

$$s( \bf {r})^{i+1} s(\bf {r})^{i} - \kappa \gamma (\bf {r}) ,$$ = (3.2)

where $\kappa$ is the step length and $\gamma(\bf {r})$ is the misfit gradient defined as:

$$\gamma (\bf {r}) = \frac{\delta \epsilon }{\delta s(\bf {r})} ,$$

 

$$\sum_{\omega} \sum_g \sum_s Real [ \Delta \tilde{P} (\bf {r}_g,\b... ...}_s)^* \frac{\partial \tilde{P}(\bf {r}_g,\bf {r}_s) }{ \partial s(\bf {r})} ].$$ = (3.3)

From here on out the Real notation will be dropped because the sum over frequencies is symmetrical and the imaginary part of the summand is anti-symmetric in $\omega$.

To find an explicit expression for the gradient, recall the perturbed Lippmann-Scwhinger equation (see appendix 1 and Stolt and Benson (1986)) which says:

 

$$\delta \tilde{P}(\bf {r}_g,\bf {r}_s) 2 \omega^2 \tilde{W}(\omega) \int_{V'} s(\bf {r}') \delta s(\bf {r}') \tilde G(\bf {r}' \vert \bf {r}_g) \tilde G(\bf {r}'\vert\bf {r}_s) dx'^3,$$ =

where $\tilde G(\bf {r}\vert \bf {r}_s)$ is the Green's function for a source at $\bf {r}_s$ and an observer at $\bf {r}$. Here, $\delta \tilde{P}(\bf {r}_g,\bf {r}_s)$ is the linearized pressure change with respect to a perturbation $\delta s(\bf {r})$ in the slowness model. The linearized field is equivalent to the solution using the Born approximation, where the slowness perturbations are too weak to generate significant multiples.

In equation 3.4, set the perturbation in slowness $\delta s(\bf {r})^{i}$ equal to $\delta (\bf {r}'-\bf {r}) \delta s(\bf {r})$ and rearrange to get the Frechet derivative:

 

$$\frac{ \delta \tilde{P}(\bf {r}_g,\bf {r}_s) }{ \delta s(\bf {r})} 2 s(\bf {r}) \omega^2 \tilde{W}(\omega) \tilde G(\bf {r}\vert \bf {r}_g) \tilde G(\bf {r}\vert\bf {r}_s).$$ = (3.4)

Substituting equation 3.5 into 3.3 yields the expression for the gradient:

$$\gamma (\bf {r}) = 2 s(\bf {r}) \sum_{\omega} \sum_g \sum_s... ...lde G(\bf {r}\vert \bf {r}_g) \tilde G(\bf {r}\vert\bf {r}_s).$$ (3.5)

It will be shown in the next section that this equation is equivalent to migration of the data residuals.

3.

Slowness Update: Substituting equation 3.6 into 3.2 yields the slowness update equation:

$$s(\bf {r})^{i+1} = s(\bf {r})^{i} - \kappa 2 s(\bf {r})^{i} ... ...lde G(\bf {r}\vert \bf {r}_g) \tilde G(\bf {r}\vert\bf {r}_s),$$ (3.6)

and steps 1-3 are repeated until convergence to $s(\bf {r})$.

Equation 3.7 can be physically interpreted by using the high frequency Green's function for the special case of a single source. We will see that the first iterate of the steepest descent method is equivalent to prestack migration.