Finite-difference Eikonal Solver Research Articles

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.

Enhancement in Seismic Imaging using Diffraction Studies and Hybrid Traveltime Technique for PSDM

The accurate migration of seismic data is conditional on the parameters which are nominated. The effective velocity used in residual processing for migration is small compared to the original migration velocity. Considering traveltime computation is a significant part of seismic imaging algorithms. Conventional implementation of Kirchhoff migration is essential for precomputing a traveltime table from the categories involving traditional ray-tracing methods and finite difference eikonal solvers. In this paper, we examine the accuracy using, the eikonal solver and paraxial ray tracing traveltime computation in pre-stack Kirchhoff depth migration. This hybrid traveltime technique can be applied to a variety of problems related to faults, fractures, and complex region. To evaluate the relevance of this identical traveltime technique, we applied on a Marmousi data set.

A comparison on initial-value ray tracing and fast marching eikonal solver for VTI traveltime computing

The Earth's subsurface is an anisotropic medium where the velocity of seismic waves alters in different propagation angles. Omitting anisotropy in seismic imaging not only brings mis-positioning of migrated dipping events but also fails to retain dipping energy during dip-moveout. To account for the efficacy of seismic anisotropy in imaging, an anisotropic wave equation must be engaged. Seismic traveltime computing is fundamental of both Kirchhoff migration and tomography algorithms. Two main categories of traveltime computing involve traditional ray tracing methods and finite difference eikonal solvers. In this study, we present two techniques of initial-value ray tracing and fast marching eikonal solver in isotropic and vertical transverse isotropy (VTI) media, and a comparison between results is demonstrated for more evaluation. Although the ray tracing approach is able to compute multiple arrivals with great precision, the eikonal solver is faster and more robust for traveltime computation. Since the ray tracing result is not a deterministic solution and it depends on the initial circumstance, employing the eikonal solver method are more preferred and suggested.

Factored singularities and high-order Lax–Friedrichs sweeping schemes for point-source traveltimes and amplitudes

In the high frequency regime, the geometrical-optics approximation for the Helmholtz equation with a point source results in an Eikonal equation for traveltime and a transport equation for amplitude. Because the point-source traveltime field has an upwind singularity at the source point, all formally high-order finite-difference Eikonal solvers exhibit first-order convergence and relatively large errors. In this paper, we propose to first factor out the singularities of traveltimes, takeoff angles, and amplitudes, and then we design high-order Lax–Friedrichs sweeping schemes for point-source traveltimes, takeoff angles, and amplitudes. Numerical examples are presented to demonstrate the performance of our new method.

A hybrid approach for tomographic inversion of crosshole seismic first-arrival times

A sequential hybrid approach was presented here to invert crosshole seismic first-arrival times. The proposed tomographic scheme combined a simple simulated annealing algorithm with a linearized smoothness-constrained least-squares inversion. The simulated annealing was implemented to obtain a background velocity distribution used by the linearized inversion for the initial guess. The linearized component was based on the functional description of traveltimes. This indicates a nonlinear function, the eikonal equation, providing traveltimes for a given slowness model. Thus an explicit ray tracing was not required by the linearized scheme. The velocity updates were obtained by a matrix inversion based on an iterative conjugate gradient-like LSQR algorithm. Second-difference regularization was used to stabilize the solutions. The Jacobian matrix giving the partial derivatives with respect to the model parameters was constructed by a finite-difference approximation based on the perturbation of the cell slowness. The traveltimes for both the hybrid and linearized schemes were calculated by a fast finite-difference eikonal solver. The hybrid scheme was tested by using both synthetic and field data sets based on the crosshole geometry. According to the tests studies, the tomograms resulted from the hybrid approach better imaged the subsurface velocity distribution. Also the hybrid optimization was characterized by quicker convergence rate than the conventional optimization based on only the linearized inversion. The tests with the synthetic data set also showed that the hybrid approach yielded a solution having lower rms residual, smaller Euclidean distance and lower relative errors in the cell velocities.

Upper crustal P-wave velocity structure of Kii Peninsula, SW Japan by first-arrival traveltime inversion of the Kawachinagano-Kiwa refraction profile

SUMMARY We investigated the upper crustal P-wave velocity structure down to 5 km in the Kii Peninsula, SW Japan by a tomographic inversion of the first-arrival traveltimes from the seismic refraction experiment performed by the Research Group for Explosion Seismology (RGES) in 1988. The observations were carried out on a profile running in the N–S direction from Kawachinagano, Osaka Prefecture to Kiwa, Mie Prefecture. The profile extends for about 65 km across the major geological zones, which characterize the geological features of SW Japan. Six shots were fired and the generated seismic waves were recorded at 86 temporary observation sites. The inversion scheme was applied to 359 first-arrival times to delineate the velocity structure along the profile. In the tomographic scheme, velocity estimation was achieved by an iterative, linearized least-squares inversion. The Jacobian matrix was constructed via a finite-difference approximation by perturbing the slownesses of the cells instead of performing ray tracing. Traveltime calculations were carried out by using a fast finite-difference eikonal solver. Velocity updates were obtained by a matrix inversion algorithm using a conjugate gradient least-squares scheme. In addition, model covariance and model resolution matrices were obtained to assess the velocity image. Two low-velocity zones observed in the northern and southern parts of the profile are the most prominent features of the tomogram. The P-wave velocity structure is generally consistent with the surface geology, and the fault zones associated with the Median (MTL) and Gobo-Hagi (GHTL) tectonic lines across the peninsula.

Traveltime tomography of crosshole radar data without ray tracing

Abstract We presented a least-squares traveltime inversion algorithm for crosshole ground-penetrating radar (GPR) direct-arrival data. The proposed scheme used the eikonal equation as traveltime functional and thus avoided tracing rays during inversion. The Jacobian matrix is constructed by a finite-difference approximation via the perturbation of slowness. The solutions were obtained by an iteratively linearized inversion approach. A smoothness-type regularization was implemented to stabilize the solutions. Traveltime calculations in forward modeling were performed by a finite-difference eikonal solver that allows modeling wavefronts. Matrix inversions were achieved by using conjugate gradient least-squares (CGLS) and LSQR algorithms. Broyden's method was used to accelerate the calculation of the Jacobian matrix when the number of model parameters was large. We tested the proposed method on three synthetic data sets and on a field data set from the Boise Hydrogeophysical Research Site (BHRS), Idaho; and we compared our model for the field data with the one obtained by a ray-tracing-based algorithm. This comparison indicated that the suggested inversion scheme was able to generate a solution as good as the one resulting from a conventional ray-based scheme. The synthetic data were obtained from simple to complex subsurface velocity distributions, including low- and high-velocity anomalies. Additionally, an image quality analysis was performed by calculating model covariance and model resolution matrices for one of the synthetic models having a complex subsurface structure and for the model resulting from the field data. All inversions were characterized by fast and stable convergence. Tests with noisy data sets indicated that the tomograms were relatively insensitive to noise in the data. It was also observed that the LSQR algorithm produced better results than the CGLS did in the tests with the synthetic models having complex subsurface structures. We considered the proposed technique to be an efficient traveltime inversion scheme for crosshole radar data.

Seismic first-arrival tomography with functional description of traveltimes

A two-dimensional smoothness-constrained least-squares inversion scheme was applied to seismic first-arrival time data. The inversion scheme was based on the functional description of traveltimes; thus, it did not require a step of ray tracing, and traveltimes were obtained by a finite-difference eikonal solver. The Laplacian difference of cell slownesses was used for the smoothness constraint. Model velocities were obtained by matrix inversion including the QR decomposition and iterative LSQR method. The Jacobian matrix of the partial derivatives was constructed by a finite-difference approximation based on the perturbation of the cell slowness. Since the construction of the Jacobian matrix was the most time-consuming step of the inversion scheme, Broyden's update was used for this matrix, and it was replaced by its numerical approximation obtained from Broyden's method after the third iteration and for all subsequent iterations to expedite the inversion process. This significantly improved the computational performance of the scheme by reducing the computer time for the calculation of the Jacobian matrix. The algorithm was tested by using a number of synthetic and field data sets. The test studies included both surface seismic refraction and crosshole seismic data acquisition configurations. Also image appraisal analyses were performed for the solutions obtained from both surface and crosshole field data sets by calculating model covariance and model resolution matrices. The presented algorithm yielded satisfactory results during the test studies. The stability and fast convergence rate were the main characteristics of the algorithm.

A practical grid-based method for tracking multiple refraction and reflection phases in three-dimensional heterogeneous media

We present a practical grid-based method in 3-D spherical coordinates for computing multiple phases comprising any number of reflection and transmission branches in heterogeneous layered media. The new scheme is based on a multistage approach which treats each layer that the wave front enters as a separate computational domain. A finite-difference eikonal solver known as the fast-marching method (FMM) is reinitialized at each interface to track the evolving wave front as either a reflection back into the incident layer or a transmission through to the adjacent layer. Unlike the standard FMM, which only finds first arrivals, this multistage approach can track those later arriving phases explicitly caused by the presence of discontinuities. Notably, the method does not require an irregular mesh to be constructed in order to connect interface nodes to neighbouring velocity nodes which lie on a regular grid. To improve accuracy, local grid refinement is used in the neighbourhood of a source point where wave front curvature is high. The method also provides a way to trace reflections from an interface that are not the first arrival (e.g. the global PP phase). These are computed by initializing the multistage FMM from both the source and receiver, propagating the two wave fronts to the reflecting interface, and finding stationary points of the sum of the two traveltime fields on the reflecting interface. A series of examples are presented to test the efficiency, accuracy and robustness of the new scheme. As well as efficiently computing various global phases to an acceptable accuracy through the ak135 model, we also demonstrate the ability of the scheme to track complex crustal phases that may be encountered in coincident reflection, wide-angle reflection/refraction or local earthquake surveys. In one example, a variety of phases are computed in the presence of a realistic subduction zone, which includes several layer pinch-outs and a subducting slab. Our numerical tests show that the new scheme is a practical and robust alternative to conventional ray tracing for finding various phases in layered media at a variety of scales.

First-Arrival Traveltime and Amplitude Calculation From Monochromatic Two-Way Wave Equation in Frequency Domain

By adding one complex frequency (the real part for frequency and the imaginary part for damping factor) into the frequency-domain wave equation, the energy or amplitude after first-arrival wave can be damped. Therefore, the calculation of first-arrival traveltime and amplitude can be transformed into the picking of the maximum energy and its phase in the monochromatic wavefield. Then, the first-arrival traveltime and amplitude can be extracted from the phase and amplitude of monochromatic frequency-domain wavefield. In order to obtain the optimum combination of the frequency and damping factor, we propose the parameter analysis scheme to scan all possible combinations. The numerical examples are compared with finite-difference eikonal solvers and mostenergetic amplitude. The results show that our algorithm has good applicability to multi-shot and multi-receiver computation in complex media.

Traveltime computation for 3D anisotropic media by a finite‐difference perturbation method

ABSTRACTThe first‐order perturbation theory is used for fast 3D computation of quasi‐compressional (qP)‐wave traveltimes in arbitrarily anisotropic media. For efficiency we implement the perturbation approach using a finite‐difference (FD) eikonal solver. Traveltimes in the unperturbed reference medium are computed with an FD eikonal solver, while perturbed traveltimes are obtained by adding a traveltime correction to the traveltimes of the reference medium. The traveltime correction must be computed along the raypath in the reference medium. Since the raypath is not determined in FD eikonal solvers, we approximate rays by linear segments corresponding to the direction of the phase normal of plane wavefronts in each cell. An isotropic medium as a reference medium works well for weak anisotropy. Using a medium with ellipsoidal anisotropy as a background medium in the perturbation approach allows us to consider stronger anisotropy without losing computational speed. The traveltime computation in media with ellipsoidal anisotropy using an FD eikonal solver is fast and accurate. The relative error is below 0.5% for the models investigated in this study. Numerical examples show that the reference model with ellipsoidal anisotropy allows us to compute the traveltime for models with strong anisotropy with an improved accuracy compared with the isotropic reference medium.

A Slowness Matching Eulerian Method for Multivalued Solutions of Eikonal Equations

Traveltime, or geodesic distance, is locally the solution of the eikonal equation of geometric optics. However traveltime between sufficiently distant points is generically multivalued. Finite difference eikonal solvers approximate only the viscosity solution, which is the smallest value of the (multivalued) traveltime (“first arrival time”). The slowness matching method stitches together local single-valued eikonal solutions, approximated by a finite difference eikonal solver, to approximate all values of the traveltime. In some applications, it is reasonable to assume that geodesics (rays) have a consistent orientation, so that the eikonal equation may be viewed as an evolution equation in one of the spatial directions. This paraxial assumption simplifies both the efficient computation of local traveltime fields and their combination into global multivalued traveltime fields via the slowness matching algorithm. The cost of slowness matching is on the same order as that of a finite difference solver used to compute the viscosity solution, when traveltimes from many point sources are required as is typical in seismic applications. Adaptive gridding near the source point and a formally third order scheme for the paraxial eikonal combine to give second order convergence of the traveltime branches.

Paraxial geometrical optics for quasi-P waves: theories and numerical methods

The quasi-P wave in anisotropic solids is of practical importance in obtaining maximal imaging resolution in seismic exploration. The geometrical optics term in the asymptotic expansion for the wave characterizes the high frequency part of the quasi-P wave by using two functions: a phase (traveltime) function satisfying an eikonal equation and an amplitude function satisfying a transport equation. Based on a paraxial eikonal equation satisfied by the traveltime corresponding to the quasi-P downgoing waves, two new advection equations for take-off angles provide essential ingredients for computing amplitude functions on uniform Cartesian grids. However, the radiation problem of the eikonal equation has an upwind singularity at the point source which renders all finite-difference eikonal solvers to be first-order accurate near the source. Extending an isotropic adaptive eikonal solver to the paraxial quasi-P eikonal equations can treat this singularity efficiently and yield highly accurate traveltimes and amplitudes. Numerical experiments for quasi-P traveltimes and amplitudes in transversely isotropic media with vertical symmetry axes verify that the numerical methods are efficient and accurate.

Wavefronts of linear elastic waves: local convexity and modeling

Seismic techniques incorporating high frequency asymptotic representation of the 3D elastic Green’s function require efficient solution methods for the computation of traveltimes. For finite difference eikonal solvers, upwind differences are requisite to sharply resolve discontinuities in the traveltime derivatives. In anisotropic media, the direction of energy propagation is not in general tangent to the wavefront normal, while finite difference eikonal solvers compute the solution based on the traveltime gradients and wavefront normal. Local convexity of the wavefronts in transverse isotropic (TI) media is proved to show that wavefront normal determines the upwind direction of the energy propagation. The eikonal equations for the traveltimes in TI media of a generally inclined symmetry axis (ITI) are derived in a way that the eikonal solvers fit conveniently. A stable, second-order, shock-capturing, upwind finite difference scheme is suggested for solving ITI eikonal equations in regular grids in 3D. Numerical experiments are presented to demonstrate the efficiency of the algorithm.

Comparison of six different methods for calculating traveltimes

We study six different methods for the calculation of seismic traveltimes. All methods yield traveltimes at all points of a regular grid.The methods examined comprise three different variants of finite‐difference (FD) eikonal solvers, the graph method, wavefront construction and a combined FD and Runge–Kutta method.The main points of investigation are computational time, accuracy and memory requirements. We took measures to obtain a high level of both general validity and clear understanding of the results. We used a profiling program to be able to measure the time that the actual core algorithm needs, thus avoiding any overhead of highly system‐dependent in‐/output operations.The comparison shows that no single method is the most appropriate but that the choice depends on the task to be fulfilled. The FD eikonal solver that uses expanding squares proves to be best suited for models which are not too complicated because it offers the best compromise between speed and accuracy, whereas wavefront construction should be applied to complex media because of its superior reliability which then justifies the much higher computational times.