Born Modelling with the Scalar Wave Equation

This section describes Born modelling, a technique used to simulate wave propagation in a perturbed medium, within the Deepwave framework. It is for readers who already know the scalar propagator and want to model, or invert for, a small perturbation to the wavespeed.

Minimal example

deepwave.scalar_born() takes the same arguments as deepwave.scalar() plus a scattering perturbation dv, and returns background and scattered receiver data together:

import torch
import deepwave

v = 1500 * torch.ones(40, 40)
dv = torch.zeros(40, 40)
dv[20, 20] = 50.0  # a small wavespeed perturbation
out = deepwave.scalar_born(
    v, grid_spacing=5, dt=0.001,
    dv=dv,
    source_amplitudes=deepwave.wavelets.ricker(25, 200, 0.001, 0.05).reshape(1, 1, -1),
    source_locations=torch.tensor([[[5, 20]]]),
    sc_receiver_locations=torch.tensor([[[35, 20]]]),
    pml_width=10,
    pml_freq=25,
)
out.sc_receiver_amplitudes  # the scattered (Born) data, shape [1, 1, 200]

Derivation

Born modelling assumes that the wavefield in a heterogeneous medium can be approximated as a sum of an unscattered (background) wavefield and a scattered wavefield. The scattered wavefield is generated by small perturbations in the medium. We present two derivations.

Derivation 1: Intuitive Approach via Taylor Expansion

Starting with a Taylor expansion of the inverse squared velocity:

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

where \(c_0\) is the background velocity and \(h_c\) is a small perturbation. If we substitute this into the scalar wave equation, and assume the total wavefield \(u\) can be approximated as the sum of an unscattered wave \(u_0\) and a singly scattered wave \(u_1\) (i.e., \(u = u_0 + u_1\)), we get:

(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 operator, \(\ddot u_0\) is the second time derivative of \(u_0\), and \(f\) is the source term. By disregarding terms beyond first-order scattering, we can separate this into two equations that intuitively describe the unscattered and scattered wavefields:

Unscattered Wavefield:

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

This equation describes the propagation of the unscattered wavefield \(u_0\) in the background medium \(c_0\), driven by the original source \(f\).

Scattered Wavefield:

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

This equation shows that the singly scattered wavefield \(u_1\) propagates in the background medium \(c_0\) and is generated by an effective source term proportional to the perturbation \(h_c\) and the second time derivative of the unscattered wavefield \(\ddot u_0\).

Derivation 2: Perturbation Theory Approach

Alternatively, we can apply a Taylor expansion directly to the wavefield with respect to the velocity perturbation:

(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 a second-order accurate finite difference approximation for the second 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 time step size. By differentiating (6) with respect to \(c\) and performing algebraic rearrangements, we can derive the source term for the scattered wavefield. This reaches the same conclusion as the first approach: the scattered wavefield can be modelled by solving the wave equation with a source term that depends on the background wavefield and the velocity perturbation.

Implementation

Deepwave’s Born propagator simultaneously simulates two wavefields:

  1. Background Wavefield: This corresponds to \(u_0\) (or \(u^t\)) from the derivations, representing the wavefield propagating in the unperturbed background medium.

  2. Scattered Wavefield: This corresponds to \(u_1\) (or \(\left.\frac{d u^t}{d c}\right|_{c_0}h_c\)), representing waves generated by interactions with the scattering model.

The propagation details, including the Perfectly Matched Layer (PML), are identical to those used for the regular scalar propagator. The key difference is that at each time step, the interaction between the background wavefield and the scattering model is computed and then used as a source term for the scattered wavefield. At the edge of the model both models are treated alike: inside the absorbing layer (and the finite-difference halo around it) the background velocity is the edge value extended outwards, and so is the scattering model. The scattered wavefield is therefore the exact linearisation of scalar with respect to the velocity for any perturbation, including one that is non-zero in the outermost cells of the model, whose extension into the layer both the linearisation of scalar and the Born propagator include. acoustic_born and elastic_born extend their perturbations the same way. Deepwave 0.0.27 and earlier instead filled the layer with zeros on scalar_born, so a perturbation of the outermost cells now produces different numbers; a call whose dv is non-zero there says so, once, under deepwave.MigrationWarning.

Inputs and Outputs

The Born propagator takes two primary model inputs:

  • Background Velocity Model (\(c_0\)): Provided as the v input.

  • Scattering Model (\(h_c\)): Provided as the scatter input.

For a small \(h_c\) (ensuring the linear approximation in (1) is valid), the output of the regular scalar propagator with a velocity model of \(c_0+h_c\) should be approximately equal to the sum of the background and scattered outputs from the Born propagator (with inputs \(c_0\) and \(h_c\)).

The output of the Born propagator carries the final wavefields of both the background and scattered fields, reflecting the propagation of two distinct wavefields, plus one opaque auxiliary_state field holding both fields’ absorbing-boundary memory variables together, and receiver data. Passing the whole output back in as initial_state continues the propagation in a later call (see Usage). Receiver data can be recorded from both wavefields:

  • bg_receiver_locations: Records data from the background wavefield.

  • sc_receiver_locations: Records data from the scattered wavefield.

Internal Time Stepping (Matrix Form)

Internally, one propagation time step involves updating both wavefields. The Deepwave implementation uses a matrix form to efficiently compute these updates. This matrix, while complex, essentially combines the update rules for both the background and scattered wavefields, including their interactions and PML effects. This detailed matrix representation is primarily an internal implementation detail for optimised performance:

\[\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}\]

Specifying sources and receivers

The channel-keyed mapping form is the canonical way to specify a survey:

out = deepwave.scalar_born(
    ...,
    sources={"u": (amplitudes, source_locations)},
    receivers={"u": receiver_locations},
    bg_receivers={"u": receiver_locations},
)

This propagator’s source channels are u and its receiver channels are u (the z channel is absent in 2-D). receivers= means the scattered receivers and bg_receivers= the background ones, matching the order of the returned tuple (background first, scattered second) and the sc_/bg_ flat spellings. The recorded data comes back as out.receivers[channel] and out.bg_receivers[channel].

The equivalent flat source_amplitudes_*/source_locations_*/ sc_receiver_locations keywords are equally permanent and produce bit-identical results; giving one channel in both forms raises. See Usage for the full vocabulary.