You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
Copy file name to clipboardExpand all lines: FWI_forward_adjoint/FWIpart1.md
+17-19Lines changed: 17 additions & 19 deletions
Original file line number
Diff line number
Diff line change
@@ -3,7 +3,7 @@
3
3
4
4
# Full-Waveform Inversion - Part 1: forward modeling
5
5
6
-
Mathias Louboutin<sup>1</sup>\*, Philipp Witte<sup>1</sup>, Michael Lange<sup>2</sup>, Navjot Kurjeka<sup>2</sup>, Fabio Luporini<sup>2</sup>, Gerard Gorman<sup>2</sup>, and Felix J. Herrmann<sup>1,3</sup>
6
+
Mathias Louboutin<sup>1</sup>\*, Philipp Witte<sup>1</sup>, Michael Lange<sup>2</sup>, Navjot Kukreja<sup>2</sup>, Fabio Luporini<sup>2</sup>, Gerard Gorman<sup>2</sup>, and Felix J. Herrmann<sup>1,3</sup>
7
7
8
8
<sup>1</sup> Seismic Laboratory for Imaging and Modeling (SLIM), The University of British Columbia
Since its re-introduction by Pratt (1999), Full-waveform 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 this topic, publications focus mostly on technical aspects, making this topic inaccessible for a broader audience due to the lack of simple introductory resources for geophysical newcomers. This is part one of three tutorials that attempt to provide an introduction and software to help people getting started. We hope to accomplish this by providing a hands-on walkthrough of FWI using Devito [2], a system based on domain specific languages that automatically generates code for time-domain finite-differences. In this capacity, 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 adjoint (in part two) wave equations. Devito releases the user from recurrent and time-consuming coding of performant time-stepping codes and allows the user to concentrate on the geophysics of the problem rather than on low-level implementation details of wave-equation simulators. This tutorial covers the conventional adjoint-state formulation of full-waveform tomography [5] that underlies most of the current methods referred to as full-waveform inversion [Matthias add ref]. While other formulations have been developed to improve the convergence of FWI for poort starting models, we will in these tutorials concentrate on the standard formulation, which relies on the combination of a forward/adjoint pair of propagators and a correlation-based 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.
18
+
Since its re-introduction by Pratt (1999), Full-waveform 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 this 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. This is part one of three tutorials that together attempt to provide an introduction to the technique and software to help people get started. We hope to accomplish this by providing a hands-on walkthrough of FWI using Devito [2], a system based on domain-specific languages that automatically generates code for time-domain finite-differences. In this capacity, 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 adjoint (in part two) wave equations. Devito frees the user from the recurrent and time-consuming development of performant time-stepping codes and allows the user to concentrate on the geophysics of the problem rather than on low-level implementation details of wave-equation simulators. This tutorial covers the conventional adjoint-state formulation of full-waveform tomography [5] that underlies most of the current methods referred to as full-waveform inversion [Matthias add ref]. 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 correlation-based 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.
19
19
20
-
Full-waveform inversion 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. The 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. This tutorial is accompanied by a Jupyter notebook - **`forward_modeling.ipynb`** -, in which 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 two of this series, we will address how to calculate model updates, i.e. gradients of the FWI objective function, via adjoint modeling. Finallt, in part three we will demonstrate how to use this gradient as part of an optimization framework for inverting an unknown velocity model.
21
20
22
-
Installation instructions for Devito are detailed at the end of the paper and required to execute the notebook.
21
+
Full-waveform inversion 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. The 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. This tutorial is accompanied by a Jupyter notebook - **`forward_modeling.ipynb`** -, in which 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 two 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 three we will demonstrate how to use this gradient as part of an optimization framework for inverting an unknown velocity model.
22
+
23
+
Devito is required to execute the notebook and can be installed following the instructions at the end of the paper.
23
24
24
25
## Wave simulations for inversion
25
26
@@ -54,32 +55,29 @@ At the core of Devito's symbolic API are two symbolic types that behave like Sym
54
55
55
56
* `Function` objects represent a spatially varying function
56
57
discretized on a regular Cartesian grid. For example, a function
57
-
symbol `f = Function(name='f', grid=model.grid space_order=2)`
58
-
is denoted symbolically as `f(x, y)`. Auto-generated symbolic
59
-
expressions for finite-difference derivatives are provided by these
60
-
objects via shorthand expressions, where the `space_order` parameter
58
+
symbol `f = Function(name='f', grid=model.grid, space_order=2)`
59
+
is denoted symbolically as `f(x, y)`. The objects provide auto-generated symbolic expressions for finite-difference derivatives through shorthand expressions like `f.dx` and `f.dx2` for the first and second derivative in `x`, where the `space_order` parameter
61
60
defines the order of the created finite-difference stencil.
62
61
63
-
* `TimeFunction` objects represent a time-dependent function
64
-
that includes leading dimension $\text{time}$, for example `g(time, x, y)`. In
62
+
* `TimeFunction` objects represent a time-dependent function that has $\text{time}$ as the leading dimension, for example `g(time, x, y)`. In
65
63
addition to spatial derivatives `TimeFunction` symbols also provide
66
64
time derivatives `g.dt` and `g.dt2`, as well as options to store
67
65
the entire data along the time axis.
68
66
69
-
To demonstrate Devito's symbolic capabilities, let us consider a time-dependent function $\mathbf{u}(\text{time}, x, y)$ representing the discrete forward wavefield. We can define this as a `TimeData` object in Devito:
67
+
To demonstrate Devito's symbolic capabilities, let us consider a time-dependent function $\mathbf{u}(\text{time}, x, y)$ representing the discrete forward wavefield. We can define this as a `TimeFunction` object in Devito:
70
68
71
69
```python
72
-
u = TimeData(name="u", shape=model.shape_domain, time_order=2,
70
+
u = TimeFunction(name="u", shape=model.shape_domain, time_order=2,
73
71
space_order=2, save=True, time_dim=nt)
74
72
```
75
73
76
74
where the `grid` object provided by the `model` defines the size
77
75
of the allocated memory region, `time_order` and `space_order` define
78
-
the default discretization order of derived derivative expressions.
76
+
the default discretization order of the derived derivative expressions.
79
77
80
78
We can now use this symbolic representation of our wavefield to
81
79
generate simple discretized stencil expressions for finite-difference
82
-
derivatives approximations using shorthand expressions, such as `u.dx` and `u.dx2` to
80
+
derivative approximations using shorthand expressions, such as `u.dx` and `u.dx2` to
83
81
denote $\frac{\partial u}{\partial x}$ and $\frac{\partial^2
84
82
u}{\partial x^2}$ respectively.
85
83
@@ -122,7 +120,7 @@ In Python, we can rearrange our `pde` expression automatically using the SymPy u
122
120
stencil = Eq(u.forward, solve(pde, u.forward)[0])
123
121
```
124
122
125
-
This `stencil` expression now represents the finite-difference approximation derived from Equation @WEstencil. This expression includes the finite-difference approximation of the Laplacian and the damping term. The `stencil` expression defines the update for a single time step only, but since the wavefield `u` is a `TimeData` object, Devito knows that we will be solving a time-dependent problem over a number of time steps.
123
+
This `stencil` expression now represents the finite-difference approximation derived from Equation @WEstencil. This expression includes the finite-difference approximation of the Laplacian and the damping term. The `stencil` expression defines the update for a single time step only, but since the wavefield `u` is a `TimeFunction` object, Devito knows that we will be solving a time-dependent problem over a number of time steps.
126
124
127
125
128
126
### Setting up the acquisition geometry
@@ -147,7 +145,7 @@ Devito also provides a special function for setting up a Ricker wavelet called `
147
145
148
146
The `src.inject` function now injects the current time sample of the Ricker wavelet (weighted with $\frac{s^2}{\mathbf{m}}$ as shown in equation @WEdisa) into the updated wavefield `u.forward` at the specified coordinates. Again, the parameter `offset` represents the size of the absorbing layer as shown in Figure #model (i.e. the source position is shifted by `offset`).
149
147
150
-
To extract the wavefield at a predetermined set of receiver locations, there is an according 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` (with `ndim=2`, since we have a two-dimensional example).
148
+
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` (with `ndim=2`, since we have a two-dimensional example).
151
149
152
150
```python
153
151
# create receiver array from receiver coordinates
@@ -208,11 +206,11 @@ In Figure 3, we show the resulting shot record. A movie of snapshots of the forw
208
206
209
207
## Conclusions
210
208
211
-
In this first part of our tutorial, we have demonstrated how to set up the discretized forward acoustic wave equations and its associated wave propagator with at runtime code generation. While we limited our discussion to the constant density acoustice 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 two of our tutorial, we will show how to calculate a valid gradient of the FWI objective using the adjoint state method. In part three, we will demonstrate how to set up a complete matrix-free and scalable optimization framework for acoustic FWI.
209
+
In this first part of the tutorial, we have demonstrated how to set up the discretized forward acoustic wave equations and its associated wave propagator with runtime code generation. While we limited our discussion to the constant density acoustice 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 two of our tutorial, we will show how to calculate a valid gradient of the FWI objective using the adjoint state method. In part three, we will demonstrate how to set up a complete matrix-free and scalable optimization framework for acoustic FWI.
212
210
213
211
### Installation
214
212
215
-
This tutorial and the coming second part are based on Devito version 3.0.3. It also require to install the full software with examples, not only the code generation API. To install devito
213
+
This tutorial and the coming second part are based on Devito version 3.0.3. It requires the installation of the full software with examples, not only the code generation API. To install devito
0 commit comments