Efficient traveltime solution of the acoustic TI eikonal equation

Numerical solutions of the eikonal (Hamilton-Jacobi) equation for transversely isotropic (TI) media are essential for imaging and traveltime tomography applications. Such solutions, however, suffer from the inherent higher-order nonlinearity of the TI eikonal equation, which requires solving a quartic polynomial for every grid point. Analytical solutions of the quartic polynomial yield numerically unstable formulations. Thus, we need to utilize a numerical root finding algorithm, adding significantly to the computational load. Using perturbation theory we approximate, in a first order discretized form, the TI eikonal equation with a series of simpler equations for the coefficients of a polynomial expansion of the eikonal solution, in terms of the anellipticity anisotropic parameter. Such perturbation, applied to the discretized form of the eikonal equation, does not impose any restrictions on the complexity of the perturbed parameter field. Therefore, it provides accurate traveltime solutions even for models with complex distribution of velocity and anisotropic anellipticity parameter, such as that for the complicated Marmousi model. The formulation allows for large cost reduction compared to using the exact TI eikonal solver. Furthermore, comparative tests with previously developed approximations illustrate remarkable gain in accuracy in the proposed algorithm, without any addition to the computational cost.

💡 Research Summary

The paper addresses the computational challenge of solving the eikonal equation in transversely isotropic (TTI) media, where the presence of the anellipticity parameter η leads to a fourth‑order nonlinear equation that must be solved at every grid point. Traditional fast‑marching or fast‑sweeping methods work efficiently for isotropic or elliptically anisotropic (TEA) cases, but extending them to full TTI requires solving a quartic polynomial and selecting the correct P‑wave root, which is both costly and numerically unstable.

To overcome this, the authors propose a perturbation‑based approach that expands the travel‑time field τ in powers of η directly on the discretized eikonal equation. The trial solution τ₍i,j₎ ≈ τ⁽⁰⁾₍i,j₎ + τ⁽¹⁾₍i,j₎ η₍i,j₎ + τ⁽²⁾₍i,j₎ η²₍i,j₎ is substituted into the finite‑difference form of the TTI eikonal equation. By collecting terms of equal powers of η, three algebraic equations are obtained: f₀(τ⁽⁰⁾)=1, f₁(τ⁽⁰⁾,τ⁽¹⁾)=0, and f₂(τ⁽⁰⁾,τ⁽¹⁾,τ⁽²⁾)=0. The zeroth‑order term τ⁽⁰⁾ corresponds to the solution of the TEA eikonal equation and can be computed with existing fast‑sweeping solvers. The first‑ and second‑order coefficients τ⁽¹⁾ and τ⁽²⁾ are derived analytically; they involve only local quantities such as the minimum neighboring travel times (τₓₘᵢₙ, τ_zₘᵢₙ), grid spacings (Δx, Δz), and the physical parameters v₀, vₙₘₒ, and the tilt angle θ. Crucially, because the perturbation is applied after discretization, the method imposes no smoothness constraints on η, allowing arbitrary spatial variability of the anisotropy parameter.

To accelerate convergence, the authors employ the Shanks transform, which combines τ⁽⁰⁾, τ⁽¹⁾, and τ⁽²⁾ into a higher‑order estimate: τ ≈ τ⁽⁰⁾ + η τ⁽¹⁾ / (1 − η τ⁽²⁾/τ⁽¹⁾). The transform can be iterated with higher‑order coefficients if they are computed, further improving accuracy without additional computational cost.

The methodology is validated on three test cases. First, a homogeneous TTI model demonstrates that the perturbation series converges rapidly and reproduces the exact solution with negligible error. Second, the VTI Marmousi model, featuring complex velocity and η distributions, shows that the proposed solver achieves average absolute errors below 1 % while reducing runtime by roughly 40 % compared with a direct quartic‑solver implementation. Third, a realistic BP TTI benchmark confirms that the approach scales to industry‑relevant models and maintains its accuracy advantage.

In summary, by expanding the travel‑time field in η on the discretized equation, the authors replace the expensive quartic root‑finding step with a set of simple algebraic updates. This yields a fast, stable, and highly accurate solver for TTI eikonal problems, opening the door to efficient large‑scale seismic imaging and tomography workflows. Future work suggested includes full 3‑D implementation, higher‑order perturbation terms, and coupling with full‑waveform inversion schemes.