# Scalar wave equation seismic amplitude derivations

## Green's functions
The wavefield at $(\mathbf{x}, t)$ due to source term $f(\mathbf{x}, t)$ is given by
$$ u(\mathbf{x}, t) = \int \int G(\mathbf{x}, t; \mathbf{x'}, t')f(\mathbf{x'}, t') d\mathbf{x'}dt', $$
where $G(\mathbf{x}, t; \mathbf{x'}, t')$ is the Green's function for the wave equation between $(\mathbf{x'}, t')$ and $(\mathbf{x}, t)$.

The scalar wave equation is
$$ \frac{1}{c^2(\mathbf{x})}\partial_{tt}u(\mathbf{x}, t) - \nabla u(\mathbf{x}, t) = f(\mathbf{x}, t) $$

For the scalar wave equation when the wave speed is constant, the Green's functions are given by

Dimension | Green's function
--- | ---
**1D** | $\frac{c}{2}\Theta(\tau - r/c)$
**2D** | $\frac{c}{2\pi \sqrt{c^2\tau^2 - r^2}}\Theta(\tau - r/c)$
**3D** | $\frac{\delta(\tau-\frac{r}{c})}{4 \pi r}$

where $r = ||\mathbf{x}-\mathbf{x'}||_2$, $\tau=t-t'$ for the causal Green's function or $\tau=t'-t$ for the anticausal Green's function, and $\Theta(x)$ is a step function.

## Direct wave

The expected amplitude of the direct, unscattered wave in a constant wave speed model can thus be computed using the equations above.

For a discretized simulation, with a source that is only active in a single cell of the model (at position $\mathbf{x_s}$), so $f(\mathbf{x}, t) = f(t)\delta(\mathbf{x} - \mathbf{x_s})v$ (where $v$ is the volume of the cell) and $r = ||\mathbf{x}-\mathbf{x_s}||_2$, the amplitudes are thus:

\begin{align}
u_{1D}(\mathbf{x}, t) &= \frac{\Delta x\Delta tc}{2}\sum_{t'=0}^{t - r/c} f(t') \\
u_{2D}(\mathbf{x}, t) &= \frac{\Delta x^2\Delta tc}{2\pi}\sum_{t'=0}^{t - r/c}\frac{f(t')}{\sqrt{c^2(t-t')^2 - r^2}} \approx \mathcal{F^{-1}} \left( \frac{i}{4}H^{(1)}_0\left(\frac{-2\pi \omega r}{c}\right)f(\omega)\Delta x^2 \right) \\
u_{3D}(\mathbf{x}, t) &= \frac{\Delta x^3f(t-\frac{r}{c})}{4 \pi r}
\end{align}

## Scattered wave

To calculate the expected amplitude of a wave that has been scattered by a wave speed anomaly in a single cell (a "point scatterer"), we need to derive the scattering amplitude.

Separating the wave speed into a smooth background ($c_0$, which we will assume varies slowly enough to be considered locally constant) and short wavelength changes ($\Delta c$), and performing a linear expansion around $\Delta c(\mathbf{x})=0$,
$$\frac{1}{c^2(\mathbf{x})} = \frac{1}{(c_0+\Delta c)^2(\mathbf{x})} \approx \frac{1}{c_0^2(\mathbf{x})} - 2 \frac{\Delta c(\mathbf{x})}{c_0^3(\mathbf{x})}.$$
We also assume that the wavefield can be split into a term that only senses the smooth model, and a scattered term,
$$u(\mathbf{x}, t) = u_0(\mathbf{x}, t) + u_{sc}(\mathbf{x}, t).$$
For the scalar wave equation, this gives
$$ \left(\frac{1}{c_0^2(\mathbf{x})}-2\frac{\Delta c(\mathbf{x})}{c_0^3(\mathbf{x})}\right)\partial_{tt}(u_0(\mathbf{x}, t)+u_{sc}(\mathbf{x},t)) - \nabla (u_0(\mathbf{x}, t)+u_{sc}(\mathbf{x},t)) = f(\mathbf{x}, t), $$
and
$$ \frac{1}{c_0^2(\mathbf{x})}\partial_{tt}u_0(\mathbf{x}, t) - \nabla u_0(\mathbf{x}, t) = f(\mathbf{x}, t). $$
Subtracting the second of these from the first gives
$$ \frac{1}{c_0^2(\mathbf{x})}\partial_{tt}u_{sc}(\mathbf{x}, t) - \nabla u_{sc}(\mathbf{x}, t) = 2\frac{\Delta c(\mathbf{x})}{c_0^3(\mathbf{x})}\partial_{tt}u(\mathbf{x}, t). $$

Making the Born approximation, we assume that multiply scattered waves have negligible amplitude. That is, the scattered wavefield $u_{sc}$ is approximated as being only generated by interactions of $u_0$ with scatterers, and further interactions of $u_{sc}$ with the scatterers are ignored. This is expressed by changing $u$ to $u_0$ in the previous equation,
$$ \frac{1}{c_0^2(\mathbf{x})}\partial_{tt}u_{sc}(\mathbf{x}, t) - \nabla u_{sc}(\mathbf{x}, t) = 2\frac{\Delta c(\mathbf{x})}{c_0^3(\mathbf{x})}\partial_{tt}u_0(\mathbf{x}, t). $$
Note that this is in the form of the wave equation, with the right-hand side as the source term. We can use this to calculate the scattered wavefield. In terms of Green's functions,
$$ u_0(\mathbf{x}, t) = \int \int G(\mathbf{x}, t; \mathbf{x'}, t') f(\mathbf{x'}, t') d\mathbf{x'} dt', $$
$$ u_{sc}(\mathbf{x}, t) = 2\int \int G(\mathbf{x}, t; \mathbf{x'}, t') \frac{\Delta c(\mathbf{x'})}{c_0^3(\mathbf{x'})}\partial_{t't'}u_0(\mathbf{x'}, t') d\mathbf{x'} dt'. $$

For a discretized simulation with a single cell source at $\mathbf{x_s}$ as above, and a wave speed model that is the constant value $c$ everywhere except at a single cell ($\mathbf{x_p}$) where it is $c+\Delta c$ (a point scatterer),  with $r=||\mathbf{x}-\mathbf{x_p}||_2$, and with $u_{1D}$, $u_{2D}$, and $u_{3D}$ from the direct wave calculation above, and using the direct wave amplitudes derived above, with $u_{p}(\mathbf{x}, t)$ being the direct wave that originates from the scatterer location with an amplitude that is the second time derivative of the direct wave that arrives at that location from the real source, the expected amplitudes of the scattered wavefields are:
\begin{equation}
u^{sc}(\mathbf{x}, t) = \frac{2 \Delta c}{c^3} u_{p}(\mathbf{x}, t).
\end{equation}

## Scatterer image/Gradient

Let's assume we only know the constant background wave speed $c$, and we want to find the location and amplitude of a point scatterer. The only information we have to help us find it is the wavefield recorded on the surface of the model that contains the scatterer, due to sources on the surface, which I will call $d(t)$. This is the typical situation when applying RTM, LSRTM, or FWI.

We'll define a mean squared error cost function that compares the wavefield at the surface when using the constant model $c$ with the data that was provided to us, $d(t)$,
$$ J = \frac{1}{T}||u_{sc}(\mathbf{x_{s}}, t) - d(t)||_2^2, $$
where $\mathbf{x_s}$ are the coordinates of the observation surface.

The derivative of this cost function with respect to changes in the model, where $T$ is the maximum recorded time, is
$$ \frac{\partial J}{\partial \Delta c(\mathbf{x})} = \frac{1}{T}\int_0^T2(u_{sc}(\mathbf{x_{s}}, t) - d(t))\frac{\partial u_{sc}}{\partial \Delta c(\mathbf{x})}dt. $$

Using our result from earlier, we get
$$ \frac{\partial u_{sc}(\mathbf{x_{s}}, t)}{\partial \Delta c(\mathbf{x})} = \frac{2v}{c_0^3(\mathbf{x})}\int_0^t G(\mathbf{x_s}, t; \mathbf{x}, t') \partial_{t't'}u_0(\mathbf{x}, t') dt', $$
where $u_0$ is, as before, the wavefield propagated in the background model, and $v$ is the volume of one cell. So,
$$ \frac{\partial J}{\partial \Delta c(\mathbf{x})} = \frac{1}{T}\int_0^T\frac{4v}{c_0^3(\mathbf{x})}(u_{sc}(\mathbf{x_{s}}, t) - d(t)) \int_0^t G(\mathbf{x_s}, t; \mathbf{x}, t') \partial_{t't'}u_0(\mathbf{x}, t') dt'dt. $$

Since $G(\mathbf{x_s}, t; \mathbf{x}, t') = 0 \forall t'>t$ when using the causal Green's function, the upper limit of the inner integral can be increased to $\infty$ (so it no longer depends on $t$ and thus the order of integration can be switched), and since the model is time invariant and due to source-receiver reciprocity, $G^+(\mathbf{x_s}, t; \mathbf{x}, t') = G^-(\mathbf{x}, t'; \mathbf{x_s}, t)$, we can write this as
$$ \frac{\partial J}{\partial \Delta c(\mathbf{x})} = \frac{1}{T}\int_0^\infty\int_T^0\frac{4v}{c_0^3(\mathbf{x})}(u_{sc}(\mathbf{x_{s}}, t) - d(t)) G^-(\mathbf{x}, t'; \mathbf{x_s}, t) \partial_{t't'}u_0(\mathbf{x}, t') dtdt'. $$

This is a product involving $u_0$ and a new wavefield I will call $u_r$, obtained by propagating the source $u_{sc}(\mathbf{x_{s}}, t) - d(t)$ backward in time:
$$ \frac{\partial J}{\partial\Delta c(\mathbf{x})} = \frac{4}{Tc_0^3(\mathbf{x})}\int u_r(\mathbf{x}, t)\partial_{tt}u_0(\mathbf{x}, t) dt. $$

This is known as an adjoint-state method means of calculating the gradient.

db6bcbb37f8bad1bfa8d2a1b94d8f6e68d684663
ipynb
Jupyter Notebook
test/amplitude_derivation.ipynb
vkazei/deepwave
032bb06328673f4f824fbca20f09ba7bb277c8d1
MIT
2018-07-16T13:57:09.000Z
2022-03-24T04:08:27.000Z
test/amplitude_derivation.ipynb
vkazei/deepwave
032bb06328673f4f824fbca20f09ba7bb277c8d1
MIT
2018-07-14T15:44:13.000Z
2022-03-25T09:35:08.000Z
test/amplitude_derivation.ipynb
vkazei/deepwave
032bb06328673f4f824fbca20f09ba7bb277c8d1
MIT
2018-12-02T14:42:59.000Z
2022-03-21T15:52:52.000Z
Qwen/Qwen-72B
1. YES
2. YES