Born modelling with the scalar wave equation

Derivation

We will start with a short but slightly hand-wavy derivation. Using a Taylor expansion,

(1)\[\frac{1}{\left(c_0+h_c\right)^2} \approx \frac{1}{c_0^2}-\frac{2h_c}{c_0^3}.\]

Inserting this into the scalar wave equation, with \(c_0\) as the background velocity model and \(h_c\) as a perturbation to it, and assuming that the wavefield can be approximated by a summation of unscattered waves, \(u_0\), and singly scattered waves, \(u_1\), yields

(2)\[\Delta \left(u_0+u_1\right) - \left(\frac{1}{c_0^2}-\frac{2h_c}{c_0^3}\right)\left(\ddot u_0+\ddot u_1\right) = f,\]

where \(\Delta\) is the Laplacian, \(\ddot u_0\) is the second time derivative of wavefield \(u_0\), and \(f\) is the source time series. Disregarding terms above first order scattering, we can separate this into two terms that match our intuition for scattering,

(3)\[\Delta u_0 - \frac{1}{c_0^2}\ddot u_0 = f,\]

and

(4)\[\Delta u_1 - \frac{1}{c_0^2}\ddot u_1 = - \frac{2h_c}{c_0^3}\ddot u_0.\]

This shows that Born modelling consists of two wavefields propagating using the background velocity model: one for unscattered waves (\(u_0\)) that uses the regular source term (\(f\)), and one for singly scattered waves (\(u_1\)) that is generated by interactions of \(u_0\) with the scattering model \(h_c\).

Now for a slightly different derivation. We apply the Taylor expansion again, but this time to the wavefield,

(5)\[u^t|_{c_0+h_c} \approx u^t|_{c_0} + \left.\frac{d u^t}{d c}\right|_{c_0}h_c,\]

where \(u^t\) is the wavefield at time step \(t\). Using the scalar wave equation, with the second order accurate finite difference approximation of the second order time derivative,

(6)\[u^t = c^2h_t^2\left(\Delta u^{t-1}-f^{t-1}\right)+2u^{t-1}-u^{t-2},\]

where \(h_t\) is the finite difference step size, we can derive the scattered term in Equation (5) that is produced by the model perturbation \(h_c\). Differentiating (6),

(7)\[\frac{d u^t}{d c} = 2ch_t^2\left(\Delta u^{t-1}-f^{t-1}\right)+c^2h_t^2\Delta\frac{d u^{t-1}}{d c}+2\frac{d u^{t-1}}{d c}-\frac{d u^{t-2}}{d c}.\]

Using, from Equation (6),

\[c^2h_t^2\left(\Delta u^{t-1}-f^{t-1}\right) = u^t -2u^{t-1}+u^{t-2},\]

rearranged into

\[2ch_t^2\left(\Delta u^{t-1}-f^{t-1}\right) = \frac{2}{c}\left(u^t -2u^{t-1}+u^{t-2}\right),\]

we get, on substituting into Equation (7),

(8)\[\frac{d u^t}{d c} = \frac{2}{c}\left(u^t-2u^{t-1}+u^{t-2}\right)+c^2h_t^2\Delta\frac{d u^{t-1}}{d c}+2\frac{d u^{t-1}}{d c}-\frac{d u^{t-2}}{d c}.\]

This implies (by comparison with the scalar wave equation, such as Equation (6)) that we can calculate \(\left.\frac{d u^t}{d c}\right|_{c_0}h_c\) by solving the wave equation with a source term of

\[-\frac{2h_c}{c_0^3h_t^2}\left(\left.u^{t+1}\right|_{c_0}-\left.2u^{t}\right|_{c_0}+\left.u^{t-1}\right|_{c_0}\right),\]

using the linearity of the wave equation to scale by \(h_c\). This is the same as before (in Equation (7)).

Implementation

Deepwave’s Born propagator simultaneously propagates two wavefields: the background wavefield (\(u_0\) and \(u^t\) in the two derivations above, respectively), and the scattered wavefield (\(u_1\) and \(\left.\frac{d u^t}{d c}\right|_{c_0}h_c\)). The details of how these wavefields are propagated, such as the PML, are the same as for the regular scalar propagator. The only difference is that at each time step the interaction of the background wavefield with the scattering model is calculated and used as a source in the scattered wavefield.

The output state from the Born propagator contains twice as many elements as the regular scalar propagator because it is propagating two wavefields. Receiver data can be recorded from both of these wavefields by specifying the bg_receiver_locations (record from the background wavefield) and receiver_locations (record from the scattered wavefield) input parameters.

The velocity model that is provided as input to the Born propagator is used as the background velocity model \(c_0\), and the scattering model that is provided as input is used as \(h_c\). For a small \(h_c\) (so the linear approximation in Equation (1) is reasonable), the output of the regular scalar propagator with the velocity model \(c_0+h_c\) should be approximately the same as the sum of the background and scattered outputs of the Born propagator with velocity model \(c_0\) and scattering model \(h_c\).

One propagation time step can be performed by applying Equation (17) to update the background wavefields. To update the scattered wavefield, the Deepwave implementation uses Equation (7) rather than (8), yielding the matrix

\[\begin{split}\begin{pmatrix} u^{t+1} \\ u^{t} \\ z^{t} \\ p^{t} \\ u_s^{t+1} \\ u_s^{t} \\ z_s^{t} \\ p_s^{t} \\ r^{t} \\ r_s^{t} \\ \end{pmatrix} = \begin{pmatrix} c^2h_t^2(1+b)\left((1+b)\partial_x^2 +\partial_x(b)\partial_x\right) + 2 & -1 & c^2h_t^2a & c^2h_t^2(1+b)\left(\partial_x a\right) & 0 & 0 & 0 & 0 & -c^2h_t^2 \\ 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ b\left((1+b)\partial_{x}^2+\partial_x(b)\partial_x\right) & 0 & a & b\left(\partial_x a\right) & 0 & 0 & 0 & 0 & 0\\ b\partial_x & 0 & 0 & a & 0 & 0 & 0 & 0 & 0 \\ 2ch_ch_t^2(1+b)\left((1+b)\partial_x^2 +\partial_x(b)\partial_x\right) & 0 & 2ch_ch_t^2a & 2ch_ch_t^2(1+b)\left(\partial_x a\right) & c^2h_t^2(1+b)\left((1+b)\partial_x^2 +\partial_x(b)\partial_x\right) + 2 & -1 & c^2h_t^2a & c^2h_t^2(1+b)\left(\partial_x a\right) & -2ch_ch_t^2 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & b\left((1+b)\partial_{x}^2+\partial_x(b)\partial_x\right) & 0 & a & b\left(\partial_x a\right) & 0\\ 0 & 0 & 0 & 0 & b\partial_x & 0 & 0 & a & 0\\ \delta_r & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & \delta_r^s & 0 & 0 & 0 & 0\\ \end{pmatrix} \begin{pmatrix} u^{t} \\ u^{t-1} \\ z^{t-1} \\ p^{t-1} \\ u_s^{t} \\ u_s^{t-1} \\ z_s^{t-1} \\ p_s^{t-1} \\ f^t \end{pmatrix}\end{split}\]