Eigenrays in 3D heterogeneous anisotropic media: Part III -- Kinematics, Finite-element implementation

TL;DR

This paper develops a finite-element numerical method for ray tracing in complex 3D anisotropic media, enabling accurate computation of stationary ray paths including caustics, based on the theory from Part I.

Contribution

It introduces a finite-element approach using Hermite interpolation for solving the nonlinear kinematic ray tracing equations in anisotropic media, advancing numerical accuracy and efficiency.

Findings

Demonstrates the method's accuracy on canonical examples

Shows improved handling of caustics in ray paths

Validates the approach's efficiency and robustness

Abstract

Following the theory presented in Part I, where we derived the Euler-Lagrange, nonlinear, second-order kinematic ray tracing equation for smooth heterogeneous general anisotropic media, this part is devoted for its numerical finite-element solution. For a given initial-guess trajectory between two fixed endpoints, discretized with a set of nodes, we update the location and direction of the ray trajectories at the nodes, to obtain the nearest stationary ray path. Starting with the Euler-Lagrange equation derived in Part I, we apply the weak formulation and the Galerkin method to reduce this second-order, ordinary differential equation into a nonlinear, local, first-order, weighted residual algebraic equation set. The solution is based on a finite element approach with the Hermite polynomial interpolation, for computing the ray characteristics between the nodes. The Hermite interpolation is a natural choice, since, in addition to the nodal locations, also the ray directions at the nodes are independent degrees of freedom of the finite element method. This is in particular important in anisotropic media. We then construct the target function to be minimized, where we distinguish between two functions triggered by the two possible types of stationary rays: minimum traveltime and saddle-point solutions, mainly in cases where the stationary paths include caustics. The target functions include two penalty terms related to two essential constraints. The first is related to the distribution of the nodes along the ray and the second enforces the normalization of the ray direction to a unit length. The minimization process involves the computation of the traveltime gradient vector and the Hessian matrix. The latter is used also to analyze the type of the stationary ray. Finally, we demonstrate the efficiency and accuracy of the proposed method along three canonical examples.