Tutorials
Juyter Logo

Part 2: Introduction to seismic modelling - 2D acoustic propagator

Dr. Rhodri Nelson
# If running on colab, uncomment the line below and run this cell:
#!pip in stall devito

In this notebook we will highlight various aspects of seismic inversion based on Devito operators. Here we aim to highlight the core ideas behind seismic modelling, where we create a numerical model that captures the processes involved in a seismic survey.

#Goals

  • Introduce the seismic Model class
  • Introduce 'sources' and 'receivers' (sparse operations)
  • Build an acoustic forward propagator

#Modelling workflow

The core process we are aiming to model is a seismic survey, which consists of two main components:

  • Source - A source is positioned at a single or a few physical locations where artificial pressure is injected into the domain we want to model. In the case of land survey, it is usually dynamite blowing up at a given location, or a vibroseis (a vibrating engine generating continuous sound waves). For a marine survey, the source is an air gun sending a bubble of compressed air into the water that will expand and generate a seismic wave.
  • Receiver - A set of microphones or hydrophones are used to measure the resulting wave and create a set of measurements called a Shot Record. These measurements are recorded at multiple locations, and usually at the surface of the domain or at the bottom of the ocean in some marine cases.

In order to create a numerical model of a seismic survey, we need to solve the wave equation and implement source and receiver interpolation to inject the source and record the seismic wave at sparse point locations in the grid.

#The acoustic seismic wave equation

The acoustic wave equation for the square slowness mm, defined as m=1c2m=\frac{1}{c^2}, where cc is the speed of sound in the given physical media, and a source qq is given by:

\begin{cases} &m \frac{\partial^2 u(\mathbf{x},t)}{\partial t^2} - \nabla^2 u(\mathbf{x},t) = q \ \text{in } \Omega \ &u(\mathbf{x},0) = 0 \ &\frac{\partial u(\mathbf{x},t)}{\partial t}|_{t=0} = 0 \end{cases}

with the zero initial conditions to guarantee unicity of the solution. The boundary conditions are Dirichlet conditions:

u(x,t)δΩ=0\begin{equation} u(\mathbf{x},t)|_{\delta\Omega} = 0 \end{equation}
(1)#

where δΩ\delta\Omega is the surface of the boundary of the model Ω\Omega.

#Finite domains

The last piece of the puzzle is the computational limitation. In the field, the seismic wave propagates in every direction to an ''infinite'' distance. However, solving the wave equation in a mathematically/discrete infinite domain is not feasible. In order to compensate, Absorbing Boundary Conditions (ABC) or Perfectly Matched Layers (PML) are required to mimic an infinite domain. These two methods allow us to approximate an unbounded media by damping and absorbing the waves at the limit of the domain to avoid reflections.

The simplest of these methods is the absorbing damping mask. (Note that we will explore the implementation of a more sophisticated boundary condition in the following notebook).The core idea is to extend the physical domain and to add a sponge mask in this extension that will absorb the incident waves. The acoustic wave equation with this damping mask can be rewritten as:

\begin{cases} &m \frac{\partial^2 u(\mathbf{x},t)}{\partial t^2} - \nabla^2 u(\mathbf{x},t) + \eta \frac{\partial u(\mathbf{x},t)}{\partial t}=q \ \text{in } \Omega \ &u(\mathbf{x},0) = 0 \ &\frac{\partial u(\mathbf{x},t)}{\partial t}|_{t=0} = 0 \end{cases}

where η\eta is the damping mask which is equal to 00 inside the physical domain and increasing inside the sponge layer. Multiple choice of profile can be chosen for η\eta from, e.g., linear to exponential.

#Seismic modelling with devito

We describe here a step by step setup of seismic modelling with Devito in a simple 2D case. We will create a physical model of our domain and define a single source along with a set of receivers. But first, we initialize some basic utilities.

import numpy as np
%matplotlib inline

#Define the physical problem

The first step is to define the physical model:

  • What are the physical dimensions of interest
  • What is the velocity profile of this physical domain

We will create a simple velocity model here by hand for demonstration purposes. This model essentially consists of two layers, each with a different velocity: 1.5km/s1.5km/s in the top layer and 2.5km/s2.5km/s in the bottom layer. We will use this simple model a lot in the following tutorials, so we will rely on a utility function to create it again later.

from examples.seismic import Model, plot_velocity

# Define a physical size
shape = (101, 101)  # Number of grid point (nx, nz)
spacing = (10., 10.)  # Grid spacing in m. The domain size is now 1km by 1km
origin = (0., 0.)  # What is the location of the top left corner. This is necessary to define
# the absolute location of the source and receivers

# Define a velocity profile. The velocity is in km/s
v = np.empty(shape, dtype=np.float32)
v[:, :51] = 1.5
v[:, 51:] = 2.5

# With the velocity and model size defined, we can create the seismic model that
# encapsulates this properties. We also define the size of the absorbing layer as 10 grid points
model = Model(vp=v, origin=origin, shape=shape, spacing=spacing,
              space_order=2, nbl=10, bcs="damp")

plot_velocity(model)
Operator `initdamp` ran in 0.01 s
Operator `pad_vp` ran in 0.01 s
<Figure size 576x432 with 2 Axes>

#Acquisition geometry

To fully define our problem setup we also need to define the source that injects the wave to the model and the set of receiver locations at which to the wavefield is sampled. The source time signature will be modelled using a 'Ricker' wavelet whose form is defined as

q(t)=(12π2f02(t1f0)2)eπ2f02(t1f0).\begin{equation} q(t) = (1-2\pi^2 f_0^2 (t - \frac{1}{f_0})^2 )e^{- \pi^2 f_0^2 (t - \frac{1}{f_0})}. \end{equation}
(2)#

To include such a source signature in our model we first need to define the time duration of our model and the timestep size, which is dictated by the CFL condition and our grid spacing. Luckily, our Model utility provides us with the critical timestep size, so we can fully discretize our model time axis as an array:

from examples.seismic import TimeAxis

t0 = 0.  # Simulation starts a t=0
tn = 1000.  # Simulation last 1 second (1000 ms)
dt = model.critical_dt  # Time step from model grid spacing

time_range = TimeAxis(start=t0, stop=tn, step=dt)

Recall that we can view the documentation for the TimeAxis object (or any other object via):

print(time_range.__doc__)

    Data object to store the TimeAxis. Exactly three of the four key arguments
    must be prescribed. Because of remainder values it is not possible to create
    a TimeAxis that exactly adhears to the inputs therefore start, stop, step
    and num values should be taken from the TimeAxis object rather than relying
    upon the input values.

    The four possible cases are:
    start is None: start = step*(1 - num) + stop
    step is None: step = (stop - start)/(num - 1)
    num is None: num = ceil((stop - start + step)/step);
                 because of remainder stop = step*(num - 1) + start
    stop is None: stop = step*(num - 1) + start

    Parameters
    ----------
    start : float, optional
        Start of time axis.
    step : float, optional
        Time interval.
    num : int, optional
        Number of values (Note: this is the number of intervals + 1).
        Stop value is reset to correct for remainder.
    stop : float, optional
        End time.
    

(Note: Values of the time steps are stored in time_range.time_values.)

We position the source at a depth of 20m20m and in the center of the xx axis (xsrc=(500m,20m)\mathbf{x}_{src}=(500m, 20m)), with a peak wavelet frequency of 10Hz10Hz.

from examples.seismic import RickerSource

f0 = 0.010  # Source peak frequency is 10Hz (0.010 kHz)
src = RickerSource(name='src', grid=model.grid, f0=f0,
                   npoint=1, time_range=time_range)

# Set source coordinates
src.coordinates.data[0, 0] = np.array(model.domain_size[0]) * .5
src.coordinates.data[0, 1] = 20.  # Depth is 20m

# We can plot the time signature to see the wavelet via:
src.show()
<Figure size 576x432 with 1 Axes>

Similarly to our source object, we can now define our receiver geometry as a symbol of type Receiver. It is worth noting here that both utility classes, RickerSource and Receiver are thin wrappers around Devito's SparseTimeFunction type, which encapsulates sparse point data and allows us to inject and interpolate values into and out of the computational grid. As we have already seen, both types provide a .coordinates property to define the position within the domain of all points encapsulated by that symbol.

In this example we will position receivers at the same depth as the source, every 10m10m along the x axis. The rec.data property will be initialized, but left empty, as we will compute the receiver readings during the simulation.

from examples.seismic import Receiver

# Create symbol for 101 receivers
rec = Receiver(name='rec', grid=model.grid, npoint=101, time_range=time_range)

# Prescribe even spacing for receivers along the x-axis
rec.coordinates.data[:, 0] = np.linspace(0, model.domain_size[0], num=101)
rec.coordinates.data[:, 1] = 20.  # Depth is 20m

# We can now show the source and receivers within our domain:
# Red dot: Source location
# Green dots: Receiver locations (every 4th point)
plot_velocity(model, source=src.coordinates.data,
              receiver=rec.coordinates.data[::4, :])
<Figure size 576x432 with 2 Axes>

#Finite-difference discretization

The finite-difference approximation is derived from Taylor expansions of the continuous field after removing the error term.

#Time discretization

We only consider the second order time discretization for now. From the Taylor expansion, the second order discrete approximation of the second order time derivative is:

2u(x,t)t2=u(x,t+Δt)2u(x,t)+u(x,tΔt)Δt2+O(Δt2).\begin{equation} \begin{aligned} \frac{\partial^2 u(\mathbf{x},t)}{\partial t^2} = \frac{\mathbf{u}(\mathbf{x},\mathbf{t+\Delta t}) - 2 \mathbf{u}(\mathbf{x},\mathbf{t}) + \mathbf{u}(\mathbf{x},\mathbf{t-\Delta t})}{\mathbf{\Delta t}^2} + O(\mathbf{\Delta t}^2). \end{aligned} \end{equation}
(3)#

where u\mathbf{u} is the discrete wavefield, Δt\mathbf{\Delta t} is the discrete time-step (distance between two consecutive discrete time points) and O(Δt2)O(\mathbf{\Delta t}^2) is the temporal discretization error term. The discretized approximation of the second order time derivative is then given by dropping the error term. This derivative is represented in Devito by u.dt2 where u is a TimeFunction object.

#Spatial discretization

We define the discrete Laplacian as the sum of the second order spatial derivatives in the three dimensions:

Δu(x,y,z,t)=j=1j=k2[αj(u(x+jdx,y,z,t)+u(xjdx,y,z,t)+u(x,y+jdy,z,t)+u(x,yjdy,zt)+u(x,y,z+jdz,t)+u(x,y,zjdz,t))]+3α0u(x,y,z,t).\begin{equation} \begin{aligned} \Delta \mathbf{u}(\mathbf{x},\mathbf{y},\mathbf{z},\mathbf{t})= \sum_{j=1}^{j=\frac{k}{2}} \Bigg[\alpha_j \Bigg(& \mathbf{u}(\mathbf{x+jdx},\mathbf{y},\mathbf{z},\mathbf{t})+\mathbf{u}(\mathbf{x-jdx},\mathbf{y},\mathbf{z},\mathbf{t}) + \\ &\mathbf{u}(\mathbf{x},\mathbf{y+jdy},\mathbf{z},\mathbf{t})+\mathbf{u}(\mathbf{x},\mathbf{y-jdy},\mathbf{z}\mathbf{t}) + \\ &\mathbf{u}(\mathbf{x},\mathbf{y},\mathbf{z+jdz},\mathbf{t})+\mathbf{u}(\mathbf{x},\mathbf{y},\mathbf{z-jdz},\mathbf{t})\Bigg) \Bigg] + \\ &3\alpha_0 \mathbf{u}(\mathbf{x},\mathbf{y},\mathbf{z},\mathbf{t}). \end{aligned} \end{equation}
(4)#

where αj\alpha_j are the appropriate finite difference weights for the approximation. This derivative is represented in Devito by u.laplace.

#Wave equation

With the space and time discretization defined, we can fully discretize the wave-equation with the combination of time and space discretizations and obtain the following second order in time and kthk^{th} order in space discrete stencil to update one grid point at position x,y,z\mathbf{x}, \mathbf{y},\mathbf{z} at time t\mathbf{t}, i.e.

u(x,y,z,t+Δt)=2u(x,y,z,t)u(x,y,z,tΔt)+Δt2m(x,y,z)(Δu(x,y,z,t)+q(x,y,z,t)).\begin{equation} \begin{aligned} \mathbf{u}(\mathbf{x},\mathbf{y},\mathbf{z},\mathbf{t+\Delta t}) = &2\mathbf{u}(\mathbf{x},\mathbf{y},\mathbf{z},\mathbf{t}) - \mathbf{u}(\mathbf{x},\mathbf{y}, \mathbf{z},\mathbf{t-\Delta t}) +\\ & \frac{\mathbf{\Delta t}^2}{\mathbf{m(\mathbf{x},\mathbf{y},\mathbf{z})}} \Big(\Delta \mathbf{u}(\mathbf{x},\mathbf{y},\mathbf{z},\mathbf{t}) + \mathbf{q}(\mathbf{x},\mathbf{y},\mathbf{z},\mathbf{t}) \Big). \end{aligned} \end{equation}
(5)#
# In order to represent the wavefield u and the square slowness we need symbolic objects 
# corresponding to time-space-varying field (u, TimeFunction) and 
# space-varying field (m, Function)
from devito import TimeFunction

# Define the wavefield with the size of the model and the time dimension
u = TimeFunction(name="u", grid=model.grid, time_order=2, space_order=2)

# We can now write the PDE
pde = model.m * u.dt2 - u.laplace + model.damp * u.dt

# The PDE representation is as on paper
pde
# This discrete PDE can be solved in a time-marching way updating u(t+dt) from the previous time step
# Devito as a shortcut for u(t+dt) which is u.forward. We can then rewrite the PDE as 
# a time marching updating equation known as a stencil using customized SymPy functions
from devito import Eq, solve

stencil = Eq(u.forward, solve(pde, u.forward))

#Source injection and receiver interpolation

With a numerical scheme to solve the homogenous wave equation, we need to add the source to introduce seismic waves and to implement the measurement operator, and interpolation operator. This operation is linked to the discrete scheme and needs to be done at the proper time step. The semi-discretized in time wave equation with a source reads:

u(x,y,z,t+Δt)=2u(x,y,z,t)u(x,y,z,tΔt)+Δt2m(x,y,z)(Δu(x,y,z,t)+q(x,y,z,t)).\begin{equation} \begin{aligned} \mathbf{u}(\mathbf{x},\mathbf{y},\mathbf{z},\mathbf{t+\Delta t}) = &2\mathbf{u}(\mathbf{x},\mathbf{y},\mathbf{z},\mathbf{t}) - \mathbf{u}(\mathbf{x},\mathbf{y}, \mathbf{z},\mathbf{t-\Delta t}) +\\ & \frac{\mathbf{\Delta t}^2}{\mathbf{m(\mathbf{x},\mathbf{y},\mathbf{z})}} \Big(\Delta \mathbf{u}(\mathbf{x},\mathbf{y},\mathbf{z},\mathbf{t}) + \mathbf{q}(\mathbf{x},\mathbf{y},\mathbf{z},\mathbf{t}) \Big). \end{aligned} \end{equation}
(6)#

It shows that in order to update u\mathbf{u} to the time t+Δt\mathbf{t+\Delta t} we must inject the value of the source term q\mathbf{q} at time t\mathbf{t}. In Devito, this corresponds to the update of uu at index t+1t+1 (t = time implicitly) with the source of time tt.

# Finally we define the source injection and receiver read function to generate the corresponding code
src_term = src.inject(field=u.forward, expr=src * dt**2 / model.m)

# Create interpolation expression for receivers
rec_term = rec.interpolate(expr=u.forward)

#Devito operator and solve

After constructing all the necessary expressions for updating the wavefield, injecting the source term and interpolating onto the receiver points, we can now create the Devito operator that will generate the C code at runtime. When creating the operator, Devito's two optimization engines will log which performance optimizations have been performed:

  • DSE: The Devito Symbolics Engine will attempt to reduce the number of operations required by the kernel.
  • DLE: The Devito Loop Engine will perform various loop-level optimizations to improve runtime performance.

Note: The argument subs=model.spacing_map causes the operator to substitute values for our current grid spacing into the expressions before code generation. This reduces the number of floating point operations executed by the kernel by pre-evaluating certain coefficients.

from devito import Operator

op = Operator([stencil] + src_term + rec_term, subs=model.spacing_map)

Now we can execute the create operator for a number of timesteps. We specify the number of timesteps to compute with the keyword time and the timestep size with dt. First, lets run for a few timesteps and visualize the wave-field:

op(time=100, dt=model.critical_dt)
Operator `Kernel` ran in 0.01 s
PerformanceSummary([(PerfKey(name='section0', rank=None), PerfEntry(time=0.0013609999999999991, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])), (PerfKey(name='section1', rank=None), PerfEntry(time=2e-06, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])), (PerfKey(name='section2', rank=None), PerfEntry(time=0.0001079999999999998, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])
from examples.seismic import plot_image
plot_image(u.data[0, :, :], cmap="seismic")
<Figure size 576x432 with 2 Axes>

Then, lets reset the wavefield and run the simulation for the full duration:

u.data[:] = 0.0
op(time=time_range.num-1, dt=model.critical_dt)
Operator `Kernel` ran in 0.01 s
PerformanceSummary([(PerfKey(name='section0', rank=None), PerfEntry(time=0.00613100000000002, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])), (PerfKey(name='section1', rank=None), PerfEntry(time=2.0000000000000005e-05, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])), (PerfKey(name='section2', rank=None), PerfEntry(time=0.0004960000000000034, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])

After running our operator kernel, the data associated with the receiver symbol rec.data has now been populated due to the interpolation expression we inserted into the operator. This allows us the visualize the shot record:

from examples.seismic import plot_shotrecord

plot_shotrecord(rec.data, model, t0, tn)
<Figure size 576x432 with 2 Axes>