Ignoring density in waveform inversion

Application to seismic inversion

Application of the Born approximation as expressed in (14) requires knowledge of the residual wavefield everywhere in the image space. Unfortunately, the full wavefield, and thus the residual wavefield, is only known at the receivers. For simplicity we assume that receivers are located at all x-locations on the surface, or that the receiver wavefield is unaliased and can be perfectly recovered. We define $w_n(\vec{r},\omega,s)$ as the background wavefield for shot $s$ at the $n^{th}$ iteration. This wavefield is computed by forward modeling the shot field through the $n^{th}$ slowness model. The data residual $\Delta d_n(x,\omega,s)$ is computed by selecting the background wavefield at ${z=0}$ and subtracting from the recorded data $d(x,\omega,s)$, or

$$\Delta d_n(x,\omega,s) = d(x,\omega,s)-\left\vert w_n(\vec{r},\omega,s) \right\vert _{z=0}.$$ (15)

The objective of the inversion is to minimize the $l^2$ norm of $\Delta d$.

Our implementation uses the linear forward operator to compute a step length at each iteration. Substituting $w_n$ for the background wavefield and selecting only the scattered field at the receivers, the frequency domain expression for the linear forward operator becomes

$$\Delta d_n(x,\omega,s) = \left\vert -\int d\vec{r}'2\omega^... ...r},\omega;\vec{r}')w_n(\vec{r}',\omega,s) \right\vert _{z=0}.$$ (16)

Since we use time-domain finite-difference modeling, it is useful to express the operator in the time domain. The $-\omega^2$ factor is applied as a second time derivative to the $w_n$ wavefield and the multiplication of $w_n$, and $G_n$ becomes a convolution along the time dimension, yielding

$$\Delta d_n(x,t,s) = \left\vert\int d\vec{r}'2\sigma_0(\vec{... ...c{r},t;\vec{r}')*\ddot{w}_n(\vec{r}',t,s) \right\vert _{z=0}.$$ (17)

The forward operator is implemented in two steps. First, the background wavefield $w_n$ is computed by propagating the source field forward in time. Next, the background wavefield is scaled by $-2\sigma_n\Delta\sigma$ and used as a new source field that is also propagated forward in time.

The gradient direction $\Delta\sigma$ for each step of the inversion is computed using the adjoint of the forward operator. The independent variables used in the forward operator are $\vec{r}$, $\vec{r}'$, $\omega$, and $s$. The forward operator integrates over $\vec{r}'$ and selects data at ${z=0}$, so the adjoint is expressed by integrating over the remaining variables and injecting data, expressed here as multiplying with a delta function, at ${z=0}$:

$$\Delta\sigma_n(\vec{r}') = -\int\!\!\!\int\!\!\!\int dsd\ome... ...G^*_n(\vec{r},\omega;\vec{r}')\delta(z)\Delta d_n(x,\omega,s)$$ (18)

This integral represents reverse-time migration of the data residual. We show a simplified expression by defining a new wavefield $res_n$ that represents the propagation of the data residual. The time axis of the Green's function is reversed due to the complex conjugate in the frequency domain:

$$res_n(\vec{r}',t,s) = \int d\vec{r}G_n(\vec{r},-t;\vec{r}')*\delta(z)\Delta d_n(x,t,s).$$ (19)

In practice, the integral is computed by forward propagating the time-reversed data residual. Due to Green's function reciprocity, integration over $\vec{r}$ is equivalent to the integration over $\vec{r}'$ in (17). The wavefield $res_n$ is then substituted into the time-domain expression for the adjoint operator where integration over frequencies is exchanged for integration over time, and the time axis of the background wavefield is reversed:

$$\Delta\sigma_n(\vec{r}') = 2\sigma_n(\vec{r}') \int\!\!\!\int dsdt \ \ddot{w}_n(\vec{r}',-t,s)\cdot res_n(\vec{r}',t,s).$$ (20)

With both the forward and adjoint linear seismic modeling expressions defined, we have all of the building blocks needed to invert for $\sigma$. We use a non-linear variation of conjugate gradients following Claerbout (2004). The method differs from linear conjugate gradients in that for each iteration the operators, which depend on $w_n$, change and the data residual $\Delta d_n$ is recomputed.

Ignoring density in waveform inversion