In this paper a fast sweeping method for computing the numerical solution of Eikonal equations on a rectangular grid is presented. The method is an iterative method which uses upwind difference for discretization and uses Gauss-Seidel iterations with alternating sweeping ordering to solve the discretized system. The crucial idea is that each sweeping ordering follows a family of characteristics of the corresponding Eikonal equation in a certain direction simultaneously. The method has an optimal complexity of O(N) for N grid points and is extremely simple to implement in any number of dimensions. Monotonicity and stability properties of the fast sweeping algorithm are proven. Convergence and error estimates of the algorithm for computing the distance function is studied in detail. It is shown that 2 n Gauss-Seidel iterations is enough for the distance function in n dimensions. An estimation of the number of iterations for general Eikonal equations is also studied. Numerical examples are used to verify the analysis.

https://www.ams.org/mcom/2005-74-250/S0025-5718-04-01678-3/S0025-5718-04-01678-3.pdf

A fast sweeping method for Eikonal equations

SUMMARY A method has been developed to compute seismic reflection traveltimes in complex 3-D velocity models with complex 3-D reflector geometry. An existing finitedifference algorithm for calculating first-arrival traveltimes was modified to handle large, sharp velocity contrasts properly. The modified algorithm is faster and more accurate than several alternative schemes, and was incorporated in a procedure to compute reflection traveltimes. Snell's law for reflections is used in the vicinity of the reflecting interface. The reflector model is allowed to vary smoothly in depth, increasing the accuracy compared with a discretized reflector model. Reflection traveltimes are computed simultaneously for all revivers, requiring only two applications of the finite-difference algorithm for each shot. This results in a significant saving in computation time in comparison with other algorithms which require one pass of the finite-difference algorithm for each shot and each receiver. The reflection traveltime procedure is well suited for incorporation in inversion schemes for 3-D velocity and reflector structure.

3-D finite-difference reflection travel times

Summary

In the summer of 1998, project INDEPTH recorded a 400 km long NNW–SSE wide-angle seismic profile in central Tibet, from the Lhasa terrane across the Banggong-Nujiang suture (BNS) at about 89.5°E and into the Qiangtang terrane. Analysis of the P-wave data reveals that (1) the crustal thickness is 65 ± 5 km beneath the line; (2) there is no 20 km step in the Moho in the vicinity of the BNS, as has been suggested to exist along-strike to the east based on prior fan profiling; (3) a thick high-velocity lower crustal layer is evident along the length of the profile (20–35 km thick, 6.5–7.3 km s−1); and (4) in contrast to the southern Lhasa terrane, there is no obvious evidence of a mid-crustal low-velocity layer in the P-wave data, although the data do not negate the possibility of such a layer of modest proportions.



Combining the results from the INDEPTH III wide-angle profile with other seismic results allows a cross-section of Moho depths to be constructed across Tibet. This cross-section shows that crustal thickness tends to decrease from south to north, with values of 70–80 km south of the middle of the Lhasa terrane, 60–70 km in the northern part of the Lhasa terrane and the Qiangtang terrane, and less than 60 km in the Qaidam basin. The overall northward thinning of the crust evident in the combined seismic observations, coupled with the essentially uniform surface elevation of the plateau south of the Qaidam basin, is supportive of the inference that northern Tibet until the Qaidam basin is underlain by somewhat thinner crust, which is isostatically supported by relatively low-density, hot upper mantle with respect to southern Tibet.

Crustal structure of central Tibet as derived from project INDEPTH wide-angle seismic data

Determining the depth and geometry of a landfill’s lower boundary is a difficult task. Potential field methods generally lack the necessary depth resolution, and seismic reflection data are usually contaminated by source‐generated noise in the time range of interest (<50 ms). To address this problem, we have developed a surface 2-D tomographic refraction scheme that is based on a fast finite‐difference eikonal solver and an inversion method that incorporates appropriate damping and smoothing constraints. This new scheme has been applied to a first‐arrival traveltime data set collected across adjacent landfills in northern Switzerland. High‐quality seismic data were collected along five profiles that crossed the landfills and two that sampled undisturbed natural sediments. Seismic waves generated from multiple shots were recorded on large numbers of closely spaced receivers during quiet evening periods. Reliability of the resultant velocity tomograms was estimated on the basis of (1) ray diagrams, (2) plot...

Refraction tomography over a buried waste disposal site

The original fast sweeping method, which is an efficient iterative method for stationary Hamilton-Jacobi equations, relies on natural ordering provided by a rectangular mesh. We propose novel ordering strategies so that the fast sweeping method can be extended efficiently and easily to any unstructured mesh. To that end we introduce multiple reference points and order all the nodes according to their $l^p$-metrics to those reference points. We show that these orderings satisfy the two most important properties underlying the fast sweeping method: (1) these orderings can cover all directions of information propagating efficiently; (2) any characteristic can be decomposed into a finite number of pieces and each piece can be covered by one of the orderings. We prove the convergence of the new algorithm. The computational complexity of the algorithm is nearly optimal in the sense that the total computational cost consists of $O(M)$ flops for iteration steps and $O(M{\rm log}M)$ flops for sorting at the predetermined initialization step which can be efficiently optimized by adopting a linear time sorting method, where $M$ is the total number of mesh points. Extensive numerical examples demonstrate that the new algorithm converges in a finite number of iterations independent of mesh size.

/pdf/fast-sweeping-methods-for-eikonal-equations-on-triangular-1d3kqs7vlm.pdf

Fast Sweeping Methods for Eikonal Equations on Triangular Meshes

Curved‐ray tomographic traveltime inversion, reverse‐time migration and various other seismic modeling applications require the calculation of traveltime and raypath information throughout a two‐ or three‐dimensional medium. When arbitrary velocity distributions and curved rays are involved, traditional ray shooting or bending procedures can be time consuming and error prone. A two‐dimensional dynamic programming traveltime computation technique, based upon Fermat’s principle, uses simple calculus techniques and a systematic mapping scheme to determine first‐arrival times on a uniform grid, given an arbitrary, discrete velocity distribution. It accurately handles large contrast, discontinuous velocity distributions, including those that generate caustics. First arrival seismic energy can travel either as transmitted waves, diffracted waves, or headwaves, and this technique models all types. The traveltime computations begin with starting values computed near the source location. Then, the mapping systemat...

Curtis Kruse

Papers

A dynamic programming approach to first arrival traveltime computation in media with arbitrarily distributed velocities