Fast Sweeping Method Research Articles

Topography-dependent eikonal tomography based on the fast-sweeping scheme and the adjoint-state technique

First-arrival traveltime tomography has been widely used for upper crustal velocity modeling, but it usually suffers from the problem of complex surface topography. To overcome this problem, we have developed a new topography-dependent eikonal tomography scheme that combines a developed accurate and efficient traveltime modeling method and introduces a flexible and robust adjoint inversion scheme in the presence of irregular topography. A surface-flattening scheme is used to handle the irregular surface, where the real model is discretized by curvilinear grids and the irregular free surface is mathematically flattened through the transformation from Cartesian to curvilinear coordinates. Based on this parameterization, the forward traveltime modeling is conducted by a monotone fast-sweeping method that discretizes the factored topography-dependent eikonal equation with a point-source condition. This algorithm can circumvent the source-singularity problem and decrease the numerical error in the vicinity of a point source in the curvilinear system. Then, the gradient-based inversion is used to minimize the misfit function, which is achieved by a matrix-free adjoint-state method without cumbersome ray tracing and explicit estimation of the Fréchet derivative matrix in the curvilinear coordinate system. The new tomographic scheme is evaluated through numerical examples with different seismic structures with complex topography, and then applied to a wide-angle profile acquired in the northeastern Tibetan Plateau. The results validate the effectiveness and efficiency of our tomography scheme in constructing shallow crustal velocity models with irregular topography.

Spatial convergence study of Lax-Friedrichs WENO fast sweeping method on the SN transport equation with nonsmoothness

Spatial convergence study of the Lax-Friedrichs fast sweeping method based on the classical third order WENO scheme (WENO3) without a priori smoothness indicators is performed on the variants of Larsen’s benchmark problem with nonsmooth solutions. For comparison, third order upwind (UPWD3), weighted diamond difference, step method, bilinear discontinuous finite element method are also included. In optically thick regimes of the cases with continuous angular fluxes, WENO3 detects steep gradients of angular fluxes and handles them as discontinuities to suppress the unphysical oscillation, which improves the accuracy significantly compared with UPWD3. In the cases with discontinuous angular fluxes, WENO3 achieves non-oscillatory solutions with limited numerical diffusion. Nonlinear Gauss-Seidel iteration is necessary to obtain convergent solution for the nonlinearity of WENO3, which leads to low computational efficiency. The convergence behavior of nonlinear Gauss-Seidel iteration is affected by nonsmoothness of the exact solution, mesh number and the relaxation parameter.

Single Pass Computation of First Seismic Wave Travel Time in Three Dimensional Heterogeneous Media With General Anisotropy

We present a numerical method for computing the first arrival travel-times of seismic waves in media defined by a general Hooke tensor, in contrast with previous methods which are limited to a specific subclass of anisotropic media, such as “tilted transversally isotropic” (TTI) media or “tilted orthorhombic” (TOR) media in Waheed et al. (G 80: 49-58, 2015), Bouteiller et al. (GJlI 212: 1498-1522, 2017). Our method proceeds in a single pass over the discretized domain, similar to the fast marching method, whereas existing methods for these types of anisotropy require multiple iterations, similar to the fast sweeping method. We introduce a new source factorization model, making it possible to achieve third-order accuracy in smooth media. We also validate our solver by comparing it with the solution of the elastic wave equation in a 3D medium with general anisotropy.

An iterative factored topography-dependent eikonal solver for anisotropic media

Accurate and efficient eikonal solvers for heterogeneous media play an important role in many areas of seismology, such as seismic tomography, migration, and earthquake localization. Incorporating seismic anisotropy and complex topography remain a computational challenge for finite-difference eikonal solvers. In recent years, the topography-dependent eikonal equation (TDEE) has been proposed as an effective way to calculate seismic traveltimes for isotropic and anisotropic media with irregular topography. However, the Lax-Friedrichs sweeping method used in previous studies to approximate the viscosity solution of TDEE for anisotropic media is more dissipative and needs a much higher number of iterations to converge. In addition, the TDEE solution for the initial point source has an upwind source singularity, which makes all TDEE solvers, even the high-order ones, exhibit polluted convergence and relatively large errors that propagate from the point source to the entire computational domain. To solve these problems, we have formulated the factored topography-dependent anisotropic eikonal (FTDAE) equation in tilted transversely isotropic (TTI) media using the factorization principle. Then, the resulting quartic equation can be numerically solved by using a fixed-point iteration technique based on the simpler elliptical anisotropic eikonal (EAE) equation with a high-order source term due to anelliptical anisotropy introduced by TTI media. At each iteration, the unknown traveltime in the EAE equation is factored into two functions: One of the functions is specified analytically to capture the source singularity, such that the unknown factor is differentiable in the source neighborhood and could be solved by the fast sweeping method. Numerical examples indicate that our FTDAE solver can treat source singularity successfully and achieve high accuracy after just a few iterations, independently of the mesh size, which could provide a more efficient and robust tool for traveltime calculation in the presence of seismic anisotropy and complex surfaces.

A fast sweeping method for calculating qP-wave traveltimes in 3-D vertical transversely isotropic media using a quadratic equation

SUMMARY Traveltime computations are an important aspect of seismic data processing applications such as traveltime tomography, migration and seismic source localization. Seismic anisotropy is a widespread feature of the Earth. Solutions to the eikonal equation that account for such anisotropy are needed for high-resolution seismic imaging and inversion. The fast sweeping method (FSM) has been widely used in computing the first-arrival traveltimes for anisotropic media because it does not need to expand the wave front from the point of the smallest traveltime. To apply FSM on strong anisotropic media, one has to solve the slowness equation derived from the Christoffel equation. All the previous developed FSM methods transform the quartic coupled slowness surface equation of quasi-P (qP) and quasi-SV (qSV) waves to the quartic equation in terms of the unknown traveltime, then numerically solve this quartic equation to compute the first-arrival traveltimes of the qP waves. However, the computational cost is significantly increased due the numerically solving the quartic equation, especially for the 3-D problems. In this study, we find a way to transform the quartic slowness equation into a quadratic one if a specific triangular-pyramid stencil around a target point is used. As the quadratic equation has the analytical solution and does not need a numerical solver, the computational efficiency of the scheme is greatly improved. We apply this methodology to develop an efficient 3-D FSM to compute the first-arrival traveltimes for qP waves in 3-D vertical transversely isotropic (VTI) media. We use both layered VTI model and complex VTI model to demonstrate the efficiency of the proposed method to obtain accurate traveltimes in 3-D VTI media involving strong anisotropy.

Absolutely convergent fixed-point fast sweeping WENO methods for steady state of hyperbolic conservation laws

Fixed-point iterative sweeping methods were developed in the literature to efficiently solve steady state solutions of Hamilton-Jacobi equations and hyperbolic conservation laws. Similar as other fast sweeping schemes, the key components of this class of methods are the Gauss-Seidel iterations and alternating sweeping strategy to achieve fast convergence rate. A family of characteristics of the corresponding hyperbolic partial differential equations (PDEs) are covered in a certain direction simultaneously in each sweeping order. Furthermore, good properties of fixed-point iterative sweeping methods which are different from other fast sweeping methods include that they have explicit forms and do not involve inverse operation of nonlinear local systems, and they can be applied to general hyperbolic equations using any monotone numerical fluxes and high order approximations easily. In a recent article (Wu et al. 2016) [31], a fifth order fixed-point sweeping WENO scheme was designed for solving steady state of hyperbolic conservation laws, and it was shown that the scheme converges to steady state solution much faster than the regular total variation diminishing (TVD) Runge-Kutta time-marching approach by stability improvement of high order schemes with a forward Euler time-marching. An open problem is that for some benchmark numerical examples, the iteration residue of the fixed-point sweeping WENO scheme hangs at a truncation error level, or even higher error levels, instead of settling down to machine zero. This issue makes it difficult to determine the convergence criterion for the iteration and challenging to apply the method to complex problems. To solve this issue, in this paper we apply the multi-resolution WENO scheme developed in Zhu and Shu (2018) [45] to the fifth order fixed-point sweeping WENO scheme and obtain an absolutely convergent fixed-point fast sweeping method for steady state of hyperbolic conservation laws, i.e., the residue of the fast sweeping iterations converges to machine zero / round off errors for all benchmark problems tested. Extensive numerical experiments, including solving difficult problems such as the shock reflection, supersonic flow past plates, and supersonic, transonic and subsonic flows past an airfoil, etc., are performed to show the accuracy, computational efficiency, and absolute convergence of the presented fifth order sweeping scheme.

A Rotating-Grid Upwind Fast Sweeping Scheme for a Class of Hamilton-Jacobi Equations

We present a fast sweeping method for a class of Hamilton-Jacobi equations that arise from time-independent problems in optimal control theory. The basic method in two dimensions uses a four point stencil and is extremely simple to implement. We test our basic method against Eikonal equations in different norms, and then suggest a general method for rotating the grid and using additional approximations to the derivatives in different directions in order to more accurately capture characteristic flow. We display the utility of our method by applying it to relevant problems from engineering.

A Continuum model for pedestrian flow with explicit consideration of crowd force and panic effects

This paper proposes a second-order pedestrian model that comprises two types of equations: continuity equation and a set of transport equations. To complete the model, we develop Eikonal equations to explicitly consider the effects of the collective decisions of individuals and crowd pressure on pedestrian dynamics. Then, the crowd movement is simulated using a set of partial differential equations under appropriate initial and boundary conditions. Based on the stability requirements derived by performing a standard linear stability analysis, suitable parameters are selected to test the model in a numerical example. The proposed second-order system is then solved using the characteristic-wise third-order weighted essentially non-oscillatory (WENO3) scheme, and the Eikonal equations are solved using the fast sweeping method. The numerical results indicate the effectiveness of the model because the derived local flow-density relationship produces a second peak in the high-density region, which is consistent with previous empirical studies. Besides, the applicability of the model to an unstable condition is verified through the simulation of complex phenomena such as stop-and-go waves. Furthermore, the estimate of crowd pressure in the simulation results can be used as a risk-level indicator for crowd management and control.

An iterative fast sweeping method for the eikonal equation in 2D anisotropic media on unstructured triangular meshes

Computation of traveltimes and raypaths is important for anisotropic tomography inversions. The eikonal-equation-based method outperforms traditional ray methods by producing more accurate results. However, most existing eikonal solvers are formulated on structured regular meshes, which are no longer accurate for models with the presence of irregular topography and subsurface interfaces. To solve the eikonal equation in vertically transversely isotropic (VTI) or tilted transversely isotropic (TTI) models with irregular geometry, we have formulated a new iterative fast sweeping method (FSM) on unstructured triangular meshes. The fixed-point iteration is implemented to capture the high-order nonlinear terms therein, and an FSM on unstructured triangular meshes is implemented to solve the resulting elliptically anisotropic eikonal equation at every iteration. We test the new algorithm for direct arrivals and reflected arrivals, and then we use the calculated traveltimes to track the raypath in VTI/TTI media. Numerical tests demonstrate the validity and accuracy of the new method for models with rough topography and subsurface interface.

Image-guided ray tracing and its applications

Eikonal solvers have important applications in seismic data processing and inversion, the so-called image-guided methods. To this day, in image-guided applications, the solution of the eikonal equation is implemented using partial-differential-equation solvers, such as fast-marching or fast-sweeping methods. We have found that alternatively, one can numerically integrate the dynamic Hamiltonian system defined by the image-guided eikonal equation and reconstruct the solution with image-guided rays. We evaluate interesting applications of image-guided ray tracing to seismic data processing, demonstrating the use of the resulting rays in image-guided interpolation and smoothing, well-log interpolation, image flattening, and residual-moveout picking. Some of these applications make use of properties of the ray-tracing system that are not directly obtained by eikonal solvers, such as ray position, ray density, wavefront curvature, and ray curvature. These ray properties open space for a different set of applications of the image-guided eikonal equation, beyond the original motivation of accelerating the construction of minimum distance tables. We stress that image-guided ray tracing is an embarrassingly parallel problem that makes its implementation highly efficient on massively parallel platforms. Image-guided ray tracing is advantageous for most applications involving the tracking of seismic events and imaging-guided interpolation. Our numerical experiments using synthetic and real data sets indicate the efficiency and robustness of image-guided rays for the selected applications.

A hybrid fast sweeping method for the isotropic eikonal equation

Eikonal solvers in the Cartesian coordinates often suffer from source singularity due to the plane-wave assumption for the wavefront near the source point. Traveltime errors induced by the source singularity near the source point will spread to the whole domain and reduce the accuracy of subsequent traveltimes. The source singularity can be avoided if the eikonal equation is formulated and solved in the polar/spherical coordinates. However, the grid in the polar/spherical coordinates will cause sparse non-uniform resampling in the region far away from the source point after the traveltime is transformed back to the Cartesian mesh, hence reducing the uniform accuracy in this region. To deal with the source singularity and maintain uniform accuracy, we introduce a hybrid fast sweeping method that computes the traveltime near the source in the polar/spherical coordinates and computes the traveltime far away from the source in the Cartesian coordinates. The source singularity near the source point is resolved in the polar/spherical coordinates, and the uniform accuracy is achieved by switching to the Cartesian coordinates away from the source point. Numerical examples are presented to demonstrate the method.

A Level-Set Adjoint-State Method for Transmission Traveltime Tomography in Irregular Domains

We propose an efficient PDE based approach for solving the first-arrival traveltime tomography problem in irregular domains. We consider the mismatch functional based on the $L_1$ or $L_2$ fidelity term, respectively, and we compute both functional gradients with respect to the sound speed by the adjoint-state method. The novelty of the proposed method consists of two aspects. First, since the tomography problem is formulated in an irregular domain, we develop new efficient level-set based fast sweeping methods for solving the eikonal and adjoint-state equations in the irregular domain. Second, since the computed adjoint state is associated with a source singularity originated from the gradient of a point source eikonal solution, we propose a novel decomposition approach to remove this singularity, which just requires solving the same adjoint-state equation again but with a different boundary condition. Finally, we show the performance of the proposed algorithm on several numerical examples including a synthetic data set from ultrasound computerized tomography.

An adaptive finite-difference method for accurate simulation of first-arrival traveltimes in heterogeneous media

Accurate simulation of first-arrival traveltimes is the foundation of seismic depth imaging and macro velocity model building for full waveform inversion. In this paper, an adaptive finite-difference method based on Gauss-Seidel iterations has been developed to achieve a more accurate calculation of first-arrival traveltimes in heterogeneous media. The plane wave, spherical wave and refracted wave approximations of local eikonal solver are used to simulate the seismic traveltimes, respectively. Then, a novel updating strategy based on Fermat principle is proposed to make the accurate choice of factor and traveltime, adaptively. The new algorithm can guarantee the high accuracy of traveltime simulation at the grid point not only near the source point, but also at a great distance from the source point. With the combination of global fast sweeping method, the traveltime simulation is unconditionally stable and all possible propagating directions of different seismic waves are considered to describe the complex propagating process due to the velocity model with sharp variations. Several numerical examples were solved to demonstrate that the adaptive finite-difference method can fulfill the accurate, efficient and robust simulation of first-arrival traveltimes.

Implicit fast sweeping method for hyperbolic systems of conservation laws

Implicit time-accurate methods are often used to integrate stiff problems where explicit schemes impose severe time step restrictions. This paper presents an efficient numerical framework based on the fast sweeping method (FSM) for solving linear and nonlinear hyperbolic systems of conservation laws. The solution at each discrete location is computed by sweeping the numerical domain in several predetermined directions that follow the causality of the characteristic families. The use of a fractional step strategy eliminates the need for a solution selection criterion while one-sided stencils limit the number of sweeps to at most 2d for d space dimensions. This work focuses on the first-order implicit upwind method since it constitutes the building block for high-order conservative schemes. For problems where the degree of stiffness evolves over time, implicit-explicit hybridization can be accomplished with the same algorithm by simply switching the stencil at each time level. As opposed to traditional implicit solvers, the sweeping method does not require a local time linearization of the fluxes thereby preserving the nonlinear stability properties of the original implicit scheme. It also avoids the large computational and memory requirements associated with solving large block-diagonal systems of equations. A series of one- and two-dimensional test cases are presented for the inviscid Burgers' equation and the reactive Euler equations. The results indicate that the implicit FSM can allow for a major reduction in the number of time steps even in the presence of discontinuous solution profiles.

A Locking Sweeping Method Based Path Planning for Unmanned Surface Vehicles in Dynamic Maritime Environments

Unmanned surface vehicles (USVs) are new marine intelligent platforms that can autonomously operate in various ocean environments with intelligent decision-making capability. As one of key technologies enabling such a capability, path planning algorithms underpin the navigation and motion control of USVs by providing optimized navigational trajectories. To accommodate complex maritime environments that include various static/moving obstacles, it is important to develop a computational efficient path planning algorithm for USVs so that real-time operation can be effectively carried out. This paper therefore proposes a new algorithm based on the fast sweeping method, named the locking sweeping method (LSM). Compared with other conventional path planning algorithms, the proposed LSM has an improved computational efficiency and can be well applied in dynamic environments that have multiple moving obstacles. When generating an optimal collision-free path, moving obstacles are modelled with ship domains that are calculated based upon ships’ velocities. To evaluate the effectiveness of the algorithm, particularly the capacity in dealing with practical environments, three different sets of simulations were undertaken in environments built using electronic nautical charts (ENCs). Results show that the proposed algorithm can effectively cope with complex maritime traffic scenarios by generating smooth and safe trajectories.

A fast-marching eikonal solver for tilted transversely isotropic media

Fast and accurate traveltime computation for quasi-P waves in anisotropic media is an essential ingredient of many seismic processing and interpretation applications such as Kirchhoff modeling and migration, microseismic source localization, and traveltime tomography. Fast-sweeping methods are widely used for solving the anisotropic eikonal equation due to their flexibility in solving general equations compared to the fast-marching method. However, it has been observed that fast sweeping can be much less efficient than fast marching for models with curved characteristics and practical grid sizes. By representing a tilted transversely isotropic (TTI) equation as a sequence of elliptically isotropic (EI) eikonal equations, we determine that the fast-marching algorithm can be used to compute fast and accurate traveltimes for TTI media. The tilt angle is absorbed into the description of the effective EI model; therefore, the adopted approach does not compromise on the solution accuracy. Through tests on benchmark synthetic models, we test our fast-marching algorithm and discover considerable improvement in accuracy by using factorization and a second-order finite-difference stencil. The adopted methodology opens the door to the possibility of using the fast-marching algorithm for a wider class of anisotropic eikonal equations.

A hybrid-order strategy to accelerate a high-order fast sweep method

A hybrid-order strategy to accelerate a high-order fast sweep method

Traveltime Calculations for qP, qSV, and qSH Waves in Two‐Dimensional Tilted Transversely Isotropic Media

AbstractThis paper presents a fast sweeping method (FSM) to calculate the first‐arrival traveltimes of the qP, qSV, and qSH waves in two‐dimensional (2D) transversely isotropic media, whose symmetry axis may have an arbitrary orientation (tilted transverse isotropy [TTI]). The method discretizes the anisotropic eikonal equation with finite difference approximations on a rectangular mesh and solves the discretized system iteratively with the Gauss‐Seidel iterations along alternating sweeping orderings. At each mesh point, a highly nonlinear equation is solved to update the numerical solution until its convergence. For solving the nonlinear equation, an interval that contains the solutions is first determined and partitioned into few subintervals such that each subinterval contains one solution; then, the false position method is applied on these subintervals to compute the solutions; after that, among all possible solutions for the discretized equation, a causality condition is imposed, and the minimum solution satisfying the causality condition is chosen to update the solution. For problems with a point‐source condition, the FSM is extended for solving the anisotropic eikonal equation after a factorization technique is applied to resolve the source singularities, which yields clean first‐order accuracy. When dealing with the triplication of the qSV wave, solutions corresponding to the minimal group velocity are chosen such that continuous solutions are computed. The accuracy, efficiency, and capability of the proposed method are demonstrated with numerical experiments.

Hybrid Fast Sweeping Methods for Anisotropic Eikonal Equation in Two-Dimensional Tilted Transversely Isotropic Media

We present hybrid fast sweeping methods for computing first-arrival traveltime of the qP, qSV and qSH waves in two-dimensional tilted transversely isotropic media, based on solving the anisotropic eikonal equation. A factorization approach is applied to resolve the source singularity near the point source, which leads to a factored anisotropic eikonal equation whose solutions can be computed with high accuracy. The proposed methods solve the factored equation in a neighborhood of the point source with the size of the neighborhood independent of the mesh, and solve the original equation outside the neighborhood. The methods enjoy all the appealing features, such as efficiency, accuracy and convergence, of the usual fast sweeping method. Furthermore, the “super-convergence” property of the first-order fast sweeping method, i.e., both its numerical solution and gradient are first-order accurate, allows us to design a second-order fast sweeping method based on a linear discontinuous Galerkin formulation. As a post-processing procedure of the first-order method, the second-order method reduces the local degrees of freedom from three to one in the linear discontinuous Galerkin formulation, which implies a simple local updating formula, hence an efficient second-order scheme. Numerical experiments are presented to demonstrate the proposed methods.

A Unified Shape-From-Shading Approach for 3D Surface Reconstruction Using Fast Eikonal Solvers

Object shape reconstruction from images has been an active topic in computer vision. Shape-from-shading (SFS) is an important approach for inferring 3D surface from a single shading image. In this paper, we present a unified SFS approach for surfaces of various reflectance properties using fast eikonal solvers. The whole approach consists of three main components: a unified SFS model, a unified eikonal-type partial differential image irradiance (PDII) equation, and fast eikonal solvers for the PDII equation. The first component is designed to address different reflectance properties including diffuse, specular, and hybrid reflections in the imaging process of the camera. The second component is meant to derive the PDII equation under an orthographic camera projection and a single distant point light source whose direction is the same as the camera. Finally, the last component is targeted at solving the resultant PDII equation by using fast eikonal solvers. It comprises two Godunov-based schemes with fast sweeping method that can handle the eikonal-type PDII equation. Experiments on several synthetic and real images demonstrate that each type of the surfaces can be effectively reconstructed with more accurate results and less CPU running time.