Fullwaveform inversion, Part 1: Forward modeling
This tutorial originally appeared as a featured article in The Leading Edge in December 2017 — see issue. 
Since its reintroduction by Pratt (1999)^{[1]}, fullwaveform inversion (FWI) has gained a lot of attention in geophysical exploration because of its ability to build high resolution velocity models more or less automatically in areas of complex geology. While there is an extensive and growing literature on the topic, publications focus mostly on technical aspects, making this topic inaccessible for a broader audience due to the lack of simple introductory resources for newcomers to geophysics. We will accomplish this by providing a handson walkthrough of FWI using Devito, a system based on domainspecific languages that automatically generates code for timedomain finitedifferences.^{[2]}
As usual, this tutorial is accompanied by all the code you need to reproduce the figures. Go to http://github.com/seg/tutorials2017 and follow the links. In the Notebook, we describe how to simulate synthetic data for a specified source and receiver setup and how to save the corresponding wavefields and shot records. In Part 2 of this series, we will address how to calculate model updates, i.e. gradients of the FWI objective function, via adjoint modeling. Finally, in Part 3 we will demonstrate how to use this gradient as part of an optimization framework for inverting an unknown velocity model.
Contents
Introduction
Devito provides a concise and straightforward computational framework for discretizing wave equations, which underlie all FWI frameworks. We will show that it generates verifiable executable code at run time for wave propagators associated with forward and (in Part 2) adjoint wave equations. Devito frees the user from the recurrent and timeconsuming development of performant timestepping codes and allows the user to concentrate on the geophysics of the problem rather than on lowlevel implementation details of waveequation simulators. This tutorial covers the conventional adjointstate formulation of fullwaveform tomography^{[3]} that underlies most of the current methods referred to as fullwaveform inversion.^{[4]} While other formulations have been developed to improve the convergence of FWI for poor starting models, in these tutorials we will concentrate on the standard formulation that relies on the combination of a forward/adjoint pair of propagators and a correlationbased gradient. In part one of this tutorial, we discuss how to set up wave simulations for inversion, including how to express the wave equation in Devito symbolically and how to deal with the acquisition geometry.
What is FWI?
FWI tries to iteratively minimize the difference between data that was acquired in a seismic survey and synthetic data that is generated from a wave simulator with an estimated (velocity) model of the subsurface. As such, each FWI framework essentially consists of a wave simulator for forward modeling the predicted data and an adjoint simulator for calculating a model update from the data misfit. This first part of this tutorial is dedicated to the forward modeling part and demonstrates how to discretize and implement the acoustic wave equation using Devito.
Wave simulations for inversion
The acoustic wave equation with the squared slowness , defined as with being the unknown spatially varying wavespeed, is given by:
( )
where is the Laplace operator, is the seismic source, located at and is a spacedependent dampening parameter for the absorbing boundary layer.^{[5]} As shown in Figure 1, the physical model is extended in every direction by nbpml
grid points to mimic an infinite domain. The dampening term attenuates the waves in the dampening layer and prevents waves from reflecting at the model boundaries. In Devito, the discrete representations of and are contained in a model
object that contains a grid object with all relevant information such as the origin of the coordinate system, grid spacing, size of the model and dimensions time, x, y
:
model = Model(vp=v, # A velocity model. origin=(0, 0), # Top left corner. shape=(101, 101), # Number of grid points. spacing=(10, 10), # Grid spacing in m. nbpml=40) # boundary layer.
In the Model
instantiation, vp
is the velocity in km/s, origin
is the origin of the physical model in meters, spacing
is the discrete grid spacing in meters, shape
is the number of grid points in each dimension and nbpml is the number of grid points in the absorbing boundary layer. Is is important to note that shape is the size of the physical domain only, while the total number of grid points, including the absorbing boundary layer, will be automatically derived from shape
and nbpml
.
Symbolic definition of the wave propagator
To model seismic data by solving the acoustic wave equation, the first necessary step is to discretize this partial differential equation (PDE), which includes discrete representations of the velocity model and wavefields, as well as approximations of the spatial and temporal derivatives using finitedifferences (FD). Unfortunately, implementing these finitedifference schemes in lowlevel code by hand is error prone, especially when we want performant and reliable code. The primary design objective of Devito is to allow users to define complex matrixfree finitedifference approximations from highlevel symbolic definitions, while employing automated code generation to create highly optimized lowlevel C code. Using the symbolic algebra package SymPy to facilitate the automatic creation of derivative expressions, Devito generates computationally efficient wave propagators.^{[6]}
At the core of Devito's symbolic API are symbolic types that behave like SymPy function objects, while also managing data:

Function
objects represent a spatially varying function discretized on a regular Cartesian grid. For example, a function symbolf = Function(name='f', grid=model.grid, space_order=2)
is denoted symbolically asf(x, y)
. The objects provide autogenerated symbolic expressions for finitedifference derivatives through shorthand expressions likef.dx
andf.dx2
for the first and second derivative inx
. 
TimeFunction
objects represent a timedependent function that has time as the leading dimension, for exampleg(time, x, y)
. In addition to spatial derivativesTimeFunction
symbols also provide time derivativesg.dt
andg.dt2
. 
SparseFunction
objects represent sparse components, such as sources and receivers, which are usually distributed sparsely and often located off the computational grid — these objects also therefore handle interpolation onto the model grid.
To demonstrate Devito's symbolic capabilities, let us consider a timedependent function u(time,x,y) representing the discrete forward wavefield:
u = TimeFunction(name="u", grid=model.grid, time_order=2, space_order=2, save=True, time_dim=nt)
where the grid object provided by the model defines the size of the allocated memory region, time_order
and space_order
define the default discretization order of the derived derivative expressions.
We can now use this symbolic representation of our wavefield to generate simple discretized expressions for finitedifference derivative approximations using shorthand expressions, such as u.dt
and u.dt2
to denote du/dt and (d^2 u)/(dt^2 ) respectively:
>>> u.dt u(time  dt, x, y)/(2*dt) + u(time + dt, x, y)/(2*dt) >>> u.dt2 2*u(time, x, y)/dt**2 + u(time  dt, x, y)/dt**2 + u(time + dt, x, y)/dt**2
Using the automatic derivation of derivative expressions, we can now implement a discretized expression for Equation 1 without the source term q(x,y,t;x_s,y_s). The
model object, which we created earlier, already contains the squared discrete slowness model.m
and damping term model.damp
as Function
objects:
pde = model.m * u.dt2  u.laplace + model.damp * u.dt
If we write out the (second order) second time derivative u.dt2
as shown earlier and ignore the damping term for the moment, our pde
expression translates to the following discrete the wave equation:
( )
with time being the current time step and dt being the time stepping interval. To propagate the wavefield, we rearrange to obtain an expression for the wavefield u(time+dt) at the next time step. Ignoring the damping term once again, this yields:
( )
We can rearrange our pde
expression automatically using the SymPy utility function solve
, then create an expression which defines the update of the wavefield for the new time step u(time+dt), with the command u.forward
:
stencil = Eq(u.forward, solve(pde, u.forward)[0])
stencil
represents the finitedifference approximation derived from Equation 3, including the finitedifference approximation of the Laplacian and the damping term.
Although it defines the update for a single time step only, Devito knows that we will be solving a timedependent problem over a number of time steps because the wavefield u is a TimeFunction
object.
Setting up the acquisition geometry
The expression for time stepping we derived in the previous section does not contain a seismic source function yet, so the update for the wavefield at a new time step is solely defined by the two previous wavefields. However as indicated in Equation 1, wavefields for seismic experiments are often excited by an active (impulsive) source q(x,y,t;x_s), which is a function of space and time (just like the wavefield u). To include such a source term in our modeling scheme, we simply add the the source wavefield as an additional term to Equation 3:
( )
Because the source appears on the righthand side in the original equation (Equation 1), the term also needs to be multiplied with dt^2/m (this follows from rearranging Equation 2, with the source on the righthand side in place of 0). Unlike the discrete wavefield u however, the source q is typically localized in space and only a function of time, which means the timedependent source wavelet is injected into the propagating wavefield at a specified source location. The same applies when we sample the wavefield at receiver locations to simulate a shot record, i.e. the simulated wavefield needs to be sampled at specified receiver locations only. Source and receiver both do not necessarily coincide with the modeling grid.
Here, RickerSource
acts as a wrapper around SparseFunction
and models a Ricker wavelet with a peak frequency f0
and source coordinates src_coords
:
f0 = 0.010 # kHz, peak frequency. src = RickerSource(name='src', grid=model.grid, f0=f0, time=time, coordinates=src_coords)
The src.inject
function now injects the current time sample of the Ricker wavelet (weighted with dt^2/m as shown in Equation 4) into the updated wavefield u.forward
at the specified coordinates.
src_term = src.inject(field=u.forward, expr=src * dt**2 / model.m, offset=model.nbpml)
To extract the wavefield at a predetermined set of receiver locations, there is a corresponding wrapper function for receivers as well, which creates a SparseFunction
object for a given number npoint
of receivers, number nt
of time samples, and specified receiver coordinates rec_coords
:
rec = Receiver(name='rec', npoint=101, ntime=nt, grid=model.grid, coordinates=rec_coords)
Rather than injecting a function into the model as we did for the source, we now simply save the wavefield at the grid points that correspond to receiver positions and interpolate the data to their exact possibly of the computational grid location:
rec_term = rec.interpolate(u, offset=model.nbpml)
Forward simulation
We can now define our forward propagator by adding the source and receiver terms to our stencil object:
op_fwd = Operator([stencil] + src_term + rec_term)
The symbolic expressions used to create Operator
contain sufficient metainformation for Devito to create a fully functional computational kernel. The dimension symbols contained in the symbolic function object (time, x, y)
define the loop structure of the created code,while allowing Devito to automatically optimize the underlying loop structure to increase execution speed.
The size of the loops and spacing between grid points is inferred from the symbolic Function
objects and associated model
. grid
object at runtime. As a result, we can invoke the generated kernel through a simple Python function call by supplying the number of time steps time
and the timestep size dt
. The user data associated with each Function
is updated inplace during operator execution, allowing us to extract the final wavefield and shot record directly from the symbolic function objects without unwanted memory duplication:
op_fwd(time=nt, dt=model.critical_dt)
When this has finished running, the resulting wavefield is stored in u.data
and the shot record is in rec.data
. We can easily plot this 2D array as an image, as shown in Figure 2.
Conclusions
In this first part of the tutorial, we have demonstrated how to set up the discretized forward acoustic wave equations and associated wave propagator with runtime code generation. While we limited our discussion to the constant density acoustic wave equation, Devito is capable of handling more general wave equations, but this is a topic beyond this tutorial on simulating waves for inversion. In Part 2 of our tutorial, we will show how to calculate a valid gradient of the FWI objective using the adjointstate method. In Part 3, we will demonstrate how to set up a complete matrixfree and scalable optimization framework for acoustic FWI.
Acknowledgments
This research was carried out as part of the SINBAD II project with the support of the member organizations of the SINBAD Consortium. This work was financially supported in part by EPSRC grant EP/L000407/1 and the Imperial College London Intel Parallel Computing Centre.
Related links
References
 ↑ Pratt, R. G., 1999, Seismic waveform inversion in the frequency domain, part 1: Theory and verification in a physical scale model: GEOPHYSICS, 64, 888–901. http://dx.doi.org/10.1190/1.1444597
 ↑ Lange, M., Kukreja, N., Louboutin, M., Luporini, F., Zacarias, F. V., Pandolfo, V., Gorman, G., 2016, Devito: Towards a generic finite difference DSL using symbolic python: 6th workshop on python for highperformance and scientific computing. http://dx.doi.org/10.1109/PyHPC.2016.9
 ↑ Tarantola, A., 1984, Inversion of seismic reflection data in the acoustic approximation: GEOPHYSICS, 49, 1259–1266. http://dx.doi.org/10.1190/1.1441754
 ↑ Virieux, J., and Operto, S., 2009, An overview of fullwaveform inversion in exploration geophysics: GEOPHYSICS, 74, WCC1–WCC26. http://dx.doi.org/10.1190/1.3238367
 ↑ Cerjan, C., Kosloff, D., Kosloff, R., and Reshef, M., 1985, A nonreflecting boundary condition for discrete acoustic and elastic wave equations: GEOPHYSICS, 50, 705–708. http://dx.doi.org/10.1190/1.1441945
 ↑ Meurer A, Smith CP, Paprocki M, et al., 2017, SymPy: symbolic computing in Python. PeerJ Computer Science 3:e103 http://dx.doi.org/10.7717/peerjcs.103
Corresponding authors
 Corresponding author: Mathias Louboutin, Seismic Laboratory for Imaging and Modeling (SLIM), The University of British Columbia,
mloubouteoas.ubc.ca
 Philipp Witte, Seismic Laboratory for Imaging and Modeling (SLIM), The University of British Columbia
 Michael Lange, Imperial College London, London, UK
 Navjot Kukreja, Imperial College London, London, UK
 Fabio Luporini, Imperial College London, London, UK
 Gerard Gorman, Imperial College London, London, UK
 Felix J. Herrmann, Seismic Laboratory for Imaging and Modeling (SLIM), The University of British Columbia, now at Georgia Institute of Technology, USA