1996, 1999, 2006
J.A. Sethian

Applications to Seismology
You are currently in the
topic outlined in red.
Overview of and references for papers on theory Overview of and references for papers on link to 
hyperbolic equations Overview of and references for level_set methods Overview of and references for on stationary 
formulation Overview of and references for Narrow Band formulation Overview of and references for papers on Fast Marching Methods Work on unstructured mesh versions of level set and fast marching methods Coupling interface methods to complex physics Adaptive mesh refinement Applications to semiconductor modeling Applications to geometry Applications to medical imaging Applications to constructing geodesics on surfaces Applications to seismology and travel times Applications to combustion Applications to fluid mechanics Applications to materials sciences Applications to robotics Applications to computer graphics Applications to CAD/CAM Applications to mesh generation

Click on navigable flow chart to go to new topic

click on any text
to go to a new topic.

First Arrivals

The key idea is to use Fast Marching Methods to compute the first arrivals in seismology. Because the technique so efficiently solves the Eikonal equation, it gives a very fast way of computing these answers. The ability to compute the first arrivals is an important part of what is known as three-dimensional prestack migration, which is the process by which one builds a model for the Earth's interior.

Three-dimensional (3D) prestack migration of surface seismic data is a tool for imaging the earth's subsurface when complex geological structures and velocity fields are present. The most commonly used imaging techniques applied to 3-D prestack surveys are methods based on the Kirchhoff integral, because of its flexibility in imaging irregularly sampled data and its relative computational efficiency. In order to perform this Kirchhoff migration, one approximately solves the wave equation with a boundary integral method. The reflectivity at every point of the earth's interior is computed by summing the recorded data on multidimensional surfaces; the shapes of the summation surfaces and the summation weights are computed from the Green's functions of the single scattering wave-propagation experiment

In 3-D seismic surveying, seismic waves are generated by surface sources (shots), and the reflected waves are recorded at surface receivers (geophones). The Green's function describes the energy of the wavefield back-scattered from the reflector point at all possible source and receiver combinations.

For 3-D prestack Kirchhoff depth migration, the Green's functions are represented by five-dimensional (5D) tables; these tables are functions of the source/receiver surface locations (x,y) and of the reflector position (x,y,z) in the earth's interior.

The key element of 3-D prestack Kirchhoff depth migration is the calculation of traveltime tables used to parameterize the asymptotic Green's functions. This is where the Fast Marching Method comes in; it allows one to compute these travel time tables *extremely* fast.


The figures show slices through the three-dimensional velocity and corresponding structural images obtained from migration on prestack data obtained from a given data set. On the left, the figure shows a depth slice through a velocity cube at a depth of $1220$ meters; on the right, the corresponding migrated image slice is shown. The salt/sediment interface and the semicircular fault cutting through the salt body are imaged with high resolution. The right figure compares the velocity model on the left with the corresponding migrated line on the right for a different slice. The sediment images are imaged at the correct locations, together with the salt body borders. The areas with lesser quality are under the salt, most probably due to the multiple reflected arrivals at this spot from the water bottom and intra-salt reflections, and also close to the left side of the top of the salt, most probably due to the use of first arrivals in the Fast Marching Method.

Multiple Arrivals

Computing first arrivals can be done efficiently and accurately using Fast Marching Methods. However, there are cases in which later arrivals contain important information. A good example of the important of these later arrivals is in geophysical imaging, in which one tries to predict what lies beneath the earth's surface by sending waves into the ground and recording their reflection. In this case, the first returning wave might not contain all the information, and later arrivals might contain more energy which can be used to more accurately predict what lies beneath.
Energy contained in many arrivals through a point
The standard approach is to use a Lagrangian ray tracing method to track the additional arrivals: this has some problems when rays diverge, and can also become a bit tricky in three dimensions. Instead, the problem can be transformed into a boundary value problem by adding an extra dimension to the problem. Let's start by thinking about a two dimensional problem. At each red spot in the domain, we can imagine three variables: the two x and y coordinates, and an additional coordinate theta corresponding to the takeoff angle from that point. Then we can imagine an arrow, leaving that point, passing and bending through the medium (the same way that light bends as passes through different media) and landing somewhere on the boundary: the place and direction it lands are called escape coordinates.
We can write down a set of equations for the escape coordinates: what's great is that they don't contain time anymore, and that means we need only solve a boundary value problem, similar to the Eikonal

equation (though with an extra dimension). We can do this by solving using a variant of fast marching methods and ordered upwind methods. Start a surface at the boundary (in three dimensional phase space): for every point on this boundary, we know the escape position and angle, since it is the same as the starting point (it's already on the exit!). Then, we can systematically march inwards, reaching back to the known escape values, and eventually cover the entire phase space cube: this is what is shown below.

Below is one result from this technique: Waves propagate from the top point through a region that contains a slowness disk in the center: you can easily see the waves propagate around the slow part, and double back on themselves, creating multiple arrivals.

New Book and Resource on Level Set and Fast Marching Methods


  1. Three dimensional traveltimes computation using the Fast Marching Method : Sethian, J.A. and Popovici, M., Geophysics, 64, 2, 1999.

    We present a fast algorithm for solving the eikonal equation in three dimensions, based on the Fast Marching Method (FMM). The algorithm is of order $O(N \log N)$, where $N$ is the total number of grid points in the computational domain. The algorithm can be used in any orthogonal coordinate system, and global constructs the solution to the Eikonal equation for each point in the coordinate domain. The method is unconditionally stable, and constructs solutions consistent with the exact solution for arbitrarily large gradient jumps in velocity. In addition, the method resolves any overturning propagation wavefronts.

    We begin with the mathematical foundation for solving the Eikonal equation using the based on the Fast Marching Method (FMM), and follow with the numerical details. We then examples of traveltime propagation through the SEG/EAGE salt dome and analyze the errors in several velocity media.

    The algorithm allows for any shape of the initial wavefront. While a point source is the most commonly used initial condition, initial plane waves can be used for controlled illumination or for downward continuation of the traveltime field from one depth to another, or from a topographic depth surface to another. The algorithm presented here is designed for computing first arrival traveltimes. Nonetheless, since it exploits the Fast Marching Method for solving the Eikonal equation, we believe that it is the fastest of all possible consistent schemes to compute first arrivals. We suggest several modifications to the basic algorithm that can lead to a most energetic arrivals traveltime computation module.

    Download publications

  2. Fast Phase Space Computation of Multiple Arrivals : Fomel, S, and Sethian, J.A., of the National Academy of Sciences, 99, 11, pp. 7329-7334, 2002

    We present a fast, general computational technique for computing the phase-space solution of static Hamilton-Jacobi equations. Starting with the Liouville formulation of the characteristic equations, we derive ``Escape Equations'' which are static, time-independent Eulerian PDEs. They represent all arrivals to the given boundary from all possible starting configurations. The solution is numerically constructed through a `one-pass' formulation, building on ideas from semi-Lagrangian methods, Dijkstra-like methods for the Eikonal equation, and Ordered Upwind Methods. To compute all possible trajectories corresponding to all possible boundary conditions, the technique is of computational order O(N \log N), where N is the total number of points in the computational phase-space domain; any particular set of boundary conditions is then extracted through rapid post-processing. Suggestions are made for speeding up the algorithm in the case when the particular distribution of sources is provided in advance. As an application, we apply the technique to the problem of computing first, multiple, and most energetic arrivals to the Eikonal equation.

    Download publications