SIAM News Blog
SIAM News

# Discrete Seismic Tomography

A seismic “camera” images Earth’s interior by probing it with elastic waves. The camera consists of transmitters and receivers placed near or on the surface of Earth. The transmitter generates elastic waves with an explosive source, which scatter from Earth’s discontinuities. Some of the scattered waves reflect back to the surface and are measured by receivers or geophones. Figure 1. Schematic of a seismic experiment. The red star represents the transmitter that generates waves (blue curves), and the green triangles denote the receivers. The transmitter moves along the x-axis and a seismic experiment is conducted for each position.
Figure 1 illustrates a typical seismic acquisition. The task of forming an image of Earth from the scattered waves is known as seismic tomography. Here, an image constitutes a map of Earth’s spatially-varying properties — for example, sound speed (referred to as velocity from this point on), density, anisotropy, and attenuation. A few challenges concerning seismic tomography are as follows:

1. Partial information: scattered waves are measured at only a few locations
2. Nonlinearity: the relation between the measurements and properties of Earth is highly nonlinear
3. Non-uniqueness: different images may produce similar measurements
4. Sharp contrasts: large discontinuities in Earth give rise to complex scattering effects
5. Noise: an improper physics model introduces a non-random, coherent noise.

These challenges make it difficult to retrieve an image of Earth from seismic measurements.

We address the challenge caused by large discontinuities in Earth, for instance, by salt bodies. Subsurface salt formations have much larger velocities than their surrounding sediments. These formations are of particular interest because hydrocarbon reservoirs are commonly found underneath or next to them. Figure 2 depicts a typical salt model. Accurate imaging of these bodies is essential for better hydrocarbon exploration. We are ultimately interested in the shape of these bodies, and imaging them is the topic of discrete seismic tomography. Figure 2. Two-dimensional salt model provided by BP. The model consists of two salt bodies (denoted by solid yellow colors). Dark blue color signifies the water.

### The Problem

One can mathematically pose the seismic tomography as an optimization problem of finding Earth’s medium parameters so as to minimize the least-squares misfit between the simulated and measured wavefield. The simulated wavefield is the solution of the wave equation for given medium properties. Since this problem is ill-posed, we introduce a regularization function alongside the least-squares misfit term. A general procedure for optimizing such cost-functions is to generate a sequence of iterates, such that the iterates converge to a minimizer. This sequence follows from the function’s gradient. However, we must note that the sequence can converge to a global minimum only if the cost function is convex. This is generally not the case, meaning that we must choose a good initial guess to converge to a local solution.

Many popular regularization approaches fail in the presence of salt bodies, usually because it is difficult to capture the sharp contrast between the layers and salt with these methods. Moreover, one cannot incorporate the salt’s almost constant velocity, which is assumed to be known in advance, with these techniques.

### The Solution Figure 3. Demonstration of the level-set approach. The top row denotes the true image. The model is 3 km deep and 10 km wide. The velocities of the layers vary from 1.5 km/s to 3.5 km/s, and the salt has a velocity of 4.5 km/s. The middle row represents the classical reconstruction while the bottom row depicts the reconstruction with a level-set approach. The red dotted line denotes true boundary of the salt.
To address this issue, discrete tomography imposes a regularization function that constrains the solution space. This constraint is a semi-discrete set consisting of a range of values for the velocity of Earth’s layers and a constant value for the velocity of the salt. For example, if Earth’s sediment layers have velocities between 1.5 and 3 km/s and the salt has a 4.5 km/s velocity, the set is {[1.5, 3], 4.5} km/s. This set is non-convex, meaning that the corresponding regularization function—which is an indicator function to this set—is non-convex and non-differentiable. To tackle this, we utilize a level-set method that separates the layers from the salt. In this approach we define an auxiliary function whose positive value denotes the salt and negative values indicate the sediment layers. A level-set of this function yields the boundary between the layers and the salt (see Figure 3). For convenience, we represent the auxiliary function with a convex combination of radial basis functions (RBFs). If we fix their centers and width, the auxiliary function’s boundary is uniquely determined by the coefficients of these basis functions.

With this formulation, we now must find the velocity profile of the layers and the coefficients of the RBFs so that the least-squares misfit is smallest. An alternating minimization technique, which minimizes the velocity of layers and the RBF coefficients alternatively, can solve the problem. The salt’s shape is then retrieved from the level-set of the auxiliary function. This approach is demonstrated in Figure 3.

Ultimately, the level-set based approach incorporates prior information about the salt velocity while simultaneously introducing nonlinearity to separate the layers from the salt. This nonlinearity can potentially create a problem during the inversion process. We will address this issue in future work.

The authors presented this research during a minisymposium at the 2017 SIAM Annual Meeting, which took place in Pittsburgh, Pa.

Acknowledgments: This work is part of the "Computational Sciences for Energy Research" initiative, a research program of the Netherlands' Organization for Scientific Research and Shell.