Stochastic Eikonal Tomography

• The paper introduces stochastic eikonal tomography as a travel-time inversion method that estimates uncertainty by modeling the velocity field as a random variable.
• It leverages Bayesian physics-informed neural networks, Gaussian processes, and variational inference to capture posterior distributions and quantify uncertainty in subsurface imaging.
• Advanced computational solvers, including fast marching and adjoint-state methods with factorization techniques, ensure accurate first-arrival computations in complex media.

Searching arXiv for recent and relevant work on stochastic eikonal tomography and related inverse problems. Stochastic eikonal tomography denotes a family of inverse problems in which first-arrival travel times, eikonal fields, or closely related high-frequency singular observables are used to recover not only a single deterministic medium but also statistical structure, posterior uncertainty, or even a full probability law. In the literature, the governing forward model is most often the eikonal equation for traveltime or phase delay, while adjacent work on random scattering shows that high-frequency singular amplitudes can encode line integrals whose correlations are tomographic projections of stochastic moments. Taken together, these works define stochastic eikonal tomography as an extension of classical traveltime imaging from coefficient recovery to uncertainty-aware or distributional inference (Gou et al., 2022, Muir, 2023, Gao et al., 31 Jul 2025, Caro et al., 2018).

1. Conceptual scope

At its narrowest, stochastic eikonal tomography is travel-time tomography with an explicit probabilistic model: the unknown slowness or velocity field is treated as random, the data are noisy, and the output is a posterior distribution or uncertainty map rather than a single point estimate. This is the role played by Bayesian physics-informed neural networks for subsurface tomography based on the eikonal equation, by Gaussian-process formulations of surface-wave eikonal tomography, and by variational-inference-based 3D eikonal imaging with per-voxel uncertainty (Gou et al., 2022, Muir, 2023, Gao et al., 31 Jul 2025).

A broader usage is also suggested by the literature. Some works are deterministic inverse methods for the eikonal equation, but they supply the forward models, linearizations, adjoint operators, and direct reconstruction formulas that stochastic methods would require. Examples include direct probing for the point-source eikonal equation, PINN-based joint $V_P$–$V_S$ inversion, adjoint-state traveltime tomography, convexification for 3D travel-time tomography, and fast marching or fast sweeping solvers for factored or anisotropic eikonal equations (Ito et al., 2023, Song et al., 2024, Chen et al., 2024, Klibanov et al., 2024, Treister et al., 2016, Waheed et al., 2013).

A third strand is only indirectly eikonal but is strongly tomographic in the same high-frequency sense. In time-domain inverse scattering from a random potential, the measured jumps of leading singularities become ray integrals of the medium, and correlations of those jumps become affine-subspace integrals of stochastic moments. Under stronger assumptions, the same data determine the full law of the random field as an $H^2(\mathbb R^n)$-valued random variable (Caro et al., 2018).

A recurring distinction in this literature is between inferential stochasticity and merely stochastic optimization. Random collocation in a PINN or stochastic gradient descent in variational inference is not, by itself, stochastic tomography in the sense of recovering uncertainty in the medium. By contrast, Bayesian PINNs, Gaussian-process formulations, and variational posterior models explicitly return uncertainty-aware reconstructions (Song et al., 2024, Gou et al., 2022, Muir, 2023, Gao et al., 31 Jul 2025).

2. Forward models, first arrivals, and computational solvers

The standard isotropic point-source formulation is the eikonal equation

$$|\nabla T(x_s,x)|^2=\frac{1}{v^2(x)}, \qquad T(x_s,x_s)=0,$$

or equivalently, in slowness form,

$|\nabla u(x)|=f(x), \qquad u(x_0)=0.$

In these formulations, the viscosity solution is interpreted as the least travel time, or first-arrival travel time, from source to receiver in the high-frequency limit (Ito et al., 2023, Gou et al., 2022). The same high-frequency interpretation underlies surface-wave eikonal tomography, where phase velocity satisfies $C_p=1/|\nabla \tau|$ after neglecting amplitude-curvature terms in the Helmholtz-based relation (Muir, 2023).

Near a point source, the eikonal field is singular, so many methods use factorizations. A common form is $T=T_0\tau$, with $T_0(x_s,x)=|x-x_s|/v(x_s)$, which regularizes the source singularity and lets the network or discrete solver approximate the smoother factor $\tau$ (Gou et al., 2022). A later joint $V_P$–$V_S$0 PINN formulation replaces that with

$V_S$1

so that the factorization depends only on source distance and no longer on source-point background velocity. The paper argues that this eliminates a nuisance dependence of the original multiplicative factorization and improves convergence speed and multiparameter inversion accuracy (Song et al., 2024).

Accurate first-arrival forward solves remain central because stochastic inversion amplifies any systematic forward-model bias. For point sources, a fast marching algorithm for the factored eikonal equation solves

$V_S$2

and derives sensitivities with respect to squared slowness $V_S$3 for use in travel-time tomography. The factored formulation removes the source singularity while retaining the computational structure of fast marching (Treister et al., 2016). In anisotropic media, a perturbative fast sweeping solver for the acoustic TI and TTI eikonal equation approximates the local quartic update by a polynomial expansion in the anellipticity parameter $V_S$4, followed by a Shanks transform. In the reported homogeneous TTI example, the preferred Shanks version had maximum error $V_S$5 ms and computational cost $V_S$6 of the direct quartic-based solver cost (Waheed et al., 2013).

Neural surrogates have also been proposed for repeated first-arrival computation. The Neural Eikonal Solver solves one-point and two-point eikonal problems with improved factorization, a Hamiltonian-based non-symmetric loss, Gaussian activation, and reciprocity-preserving symmetrization. In the reported tests it achieved relative-mean-absolute error of about $V_S$7-$V_S$8 from the second-order factored Fast Marching Method, with $V_S$9-$H^2(\mathbb R^n)$0 times lower errors and $H^2(\mathbb R^n)$1-$H^2(\mathbb R^n)$2 times faster training than prior neural solvers; the two-point formulation was proposed as a compact surrogate for applications involving millions of source-receiver pairs (Grubas et al., 2022).

For large-scale seismic tomography, anisotropic adjoint-state implementations generalize the forward model beyond isotropic media. TomoATT solves an anisotropic eikonal equation in spherical coordinates,

$H^2(\mathbb R^n)$3

for slowness $H^2(\mathbb R^n)$4 and azimuthal anisotropy parameters $H^2(\mathbb R^n)$5 and $H^2(\mathbb R^n)$6, using the fast sweeping method, multiplicative factorization to eliminate source singularity, and adjoint-state gradients (Chen et al., 2024).

3. Deterministic inversion backbones

The deterministic core of eikonal tomography is strongly ill posed. A direct probing study of the point-source eikonal equation explicitly states that the inverse problem is “highly” or “severely” ill-posed, and gives a concrete nonuniqueness example in which boundary travel times on an upper segment are independent of the inclusion contrast $H^2(\mathbb R^n)$7 once $H^2(\mathbb R^n)$8 is sufficiently large (Ito et al., 2023). This nonuniqueness is intrinsic to first-arrival data: once rays bypass a sufficiently slow inclusion, increasing the slowness further need not change the observed first arrivals.

In low-contrast media, the same paper linearizes the eikonal equation around a homogeneous background and shows that the travel-time perturbation behaves like a fanbeam transform of the slowness perturbation. This yields a direct probing reconstruction by inverse fanbeam tomography and filtered back projection,

$H^2(\mathbb R^n)$9

followed by a nonlinear refinement step (Ito et al., 2023). In high-contrast media, the method instead assumes a background $$|\nabla T(x_s,x)|^2=\frac{1}{v^2(x)}, \qquad T(x_s,x_s)=0,$$0, linearizes around its traveltime $$|\nabla T(x_s,x)|^2=\frac{1}{v^2(x)}, \qquad T(x_s,x_s)=0,$$1, and uses an adjoint transport equation to back-project boundary residuals. This produces a direct imaging estimate of $$|\nabla T(x_s,x)|^2=\frac{1}{v^2(x)}, \qquad T(x_s,x_s)=0,$$2 without posing a full least-squares optimization (Ito et al., 2023).

A different deterministic strategy is convexification. For the 3D problem of travel-time tomography in a circular cylinder with sources on the axis and measurements on the whole boundary, the coefficient inverse problem for the eikonal equation is transformed by introducing $$|\nabla T(x_s,x)|^2=\frac{1}{v^2(x)}, \qquad T(x_s,x_s)=0,$$3, differentiating with respect to the source coordinate $$|\nabla T(x_s,x)|^2=\frac{1}{v^2(x)}, \qquad T(x_s,x_s)=0,$$4, expanding in a specially constructed basis, and minimizing a Carleman-weighted residual functional. The resulting method is proved globally strongly convex for sufficiently large Carleman parameter and comes with global-convergence and noise-stability statements in the semi-discrete setting (Klibanov et al., 2024).

Adjoint-state traveltime tomography provides another backbone. TomoATT minimizes weighted combinations of absolute traveltimes, common-source differential arrival times, and common-receiver differential arrival times, using forward eikonal solves and adjoint equations to construct sensitivity kernels for slowness and azimuthal anisotropy (Chen et al., 2024). The implementation employs multiple-grid parameterization, kernel-density normalization, and multi-level parallelization. The paper reports that each inversion starts with only three simple input files and completes within 2 hours using 64 processors; in a Parkfield application, 80 iterations on 64 processors took 57 minutes, and in a Thailand teleseismic application 80 iterations took 124 minutes (Chen et al., 2024).

Fast marching also supports inversion directly. The factored fast marching paper formulates travel-time tomography as minimizing a data misfit plus regularization over squared slowness $$|\nabla T(x_s,x)|^2=\frac{1}{v^2(x)}, \qquad T(x_s,x_s)=0,$$5, and derives a sparse Jacobian for the factored eikonal solver. Because the linearized operator is lower triangular in fast-marching acceptance order, Jacobian applications reduce to forward substitution, which is especially valuable in Gauss–Newton inversion (Treister et al., 2016).

4. Probabilistic formulations and uncertainty quantification

The most explicit stochastic formulations treat the eikonal inverse problem as Bayesian inference. In Bayesian physics-informed neural networks for subsurface tomography, two neural networks represent the factorized travel-time field $$|\nabla T(x_s,x)|^2=\frac{1}{v^2(x)}, \qquad T(x_s,x_s)=0,$$6 and the normalized velocity field $$|\nabla T(x_s,x)|^2=\frac{1}{v^2(x)}, \qquad T(x_s,x_s)=0,$$7, while the eikonal residual enters as a physics likelihood. The posterior is

$$|\nabla T(x_s,x)|^2=\frac{1}{v^2(x)}, \qquad T(x_s,x_s)=0,$$8

with Gaussian priors and Gaussian likelihoods for data misfit and residual misfit; posterior predictive means and variances are estimated by Monte Carlo over posterior samples (Gou et al., 2022). The paper compares mean-field Gaussian variational inference with Stein Variational Gradient Descent and reports that SVGD was closer to the analytic posterior in a 1D benchmark and consistently better calibrated than VI in 2D examples, while VI was overconfident (Gou et al., 2022).

Surface-wave Bayesian eikonal tomography using Gaussian processes shifts the stochastic object from the slowness field to the phase-delay field. The observed phase delays are modeled as

$$|\nabla T(x_s,x)|^2=\frac{1}{v^2(x)}, \qquad T(x_s,x_s)=0,$$9

with $|\nabla u(x)|=f(x), \qquad u(x_0)=0.$0 and $|\nabla u(x)|=f(x), \qquad u(x_0)=0.$1 a Gaussian process. Because differentiation is linear, the posterior for $|\nabla u(x)|=f(x), \qquad u(x_0)=0.$2 is available analytically, and the nonlinear map from gradient to phase velocity,

$|\nabla u(x)|=f(x), \qquad u(x_0)=0.$3

is handled by a saddlepoint approximation to the induced non-Gaussian posterior (Muir, 2023). A central implication is that standard eikonal tomography often underestimates uncertainty because it treats the interpolated phase-delay field as known exactly rather than uncertain (Muir, 2023).

A large-scale voxelwise formulation appears in the 3D Kīlauea study. There, the unknown $|\nabla u(x)|=f(x), \qquad u(x_0)=0.$4 is a 3D voxel grid with $|\nabla u(x)|=f(x), \qquad u(x_0)=0.$5 voxels, and the approximate posterior is chosen as a multivariate lognormal distribution with diagonal covariance in log-space. Variational inference minimizes

$|\nabla u(x)|=f(x), \qquad u(x_0)=0.$6

with $|\nabla u(x)|=f(x), \qquad u(x_0)=0.$7 set to $|\nabla u(x)|=f(x), \qquad u(x_0)=0.$8 and regularization supplied implicitly by annealed gradient smoothing (Gao et al., 31 Jul 2025). The real-data application uses 18 receivers, 45,310 earthquakes, 378,074 P arrivals, and 280,152 S arrivals; uncertainty is computed by sampling the variational posterior and taking per-voxel standard deviations (Gao et al., 31 Jul 2025). The paper reports $|\nabla u(x)|=f(x), \qquad u(x_0)=0.$9 km resolution in the upper half of the domain in the best-covered areas and at least $C_p=1/|\nabla \tau|$0 km in much of the ray-covered region, and it estimates a deep magma-column volume of $C_p=1/|\nabla \tau|$1, melt fraction $C_p=1/|\nabla \tau|$2 at 11 km depth, and total melt volume $C_p=1/|\nabla \tau|$3 (Gao et al., 31 Jul 2025).

A representative cross-section of these probabilistic formulations is as follows.

Formulation	Stochastic quantity	Main mechanism
BPINN tomography (Gou et al., 2022)	Posterior over network parameters; predictive mean and variance of velocity and traveltime	Physics-informed likelihood with VI or SVGD
GP eikonal tomography (Muir, 2023)	Posterior over phase delay, phase-delay gradient, squared slowness, and phase velocity	GP interpolation with analytic gradient posterior and saddlepoint approximation
Kīlauea VI tomography (Gao et al., 31 Jul 2025)	Voxelwise lognormal posterior over 3D velocity	Variational inference with adjoint eikonal gradients and stochastic optimization

Deterministic PINN-based joint $C_p=1/|\nabla \tau|$4–$C_p=1/|\nabla \tau|$5 inversion is closely related but should be distinguished from these stochastic formulations. The PINNPStomo framework jointly fits P- and S-wave first-arrival traveltimes while enforcing two eikonal equations and using a source-conditioned traveltime network plus a source-independent velocity network. Its relevance to stochastic tomography is structural rather than inferential: it provides a differentiable inverse parameterization, but it does not itself return posterior uncertainty (Song et al., 2024).

5. Correlation imaging and tomography of random-medium laws

An important extension of stochastic eikonal thinking arises from random scattering rather than direct eikonal inversion. In the wave-equation model

$C_p=1/|\nabla \tau|$6

with $C_p=1/|\nabla \tau|$7 a compactly supported $C_p=1/|\nabla \tau|$8-valued random field, the measurable quantity is the jump of the leading singularity of the scattered wave across characteristic hyperplanes $C_p=1/|\nabla \tau|$9. For each incident direction, the distant jump amplitude is exactly a ray transform of the random potential:

$T=T_0\tau$0

Thus the first singular arrival carries a line integral of the medium, in direct analogy with high-frequency geometric tomography (Caro et al., 2018).

The stochastic object of interest is the hierarchy of moment maps

$T=T_0\tau$1

Substituting the jump formula into the $T=T_0\tau$2-point exterior correlation data yields

$T=T_0\tau$3

where $T=T_0\tau$4 is a $T=T_0\tau$5-dimensional affine subspace (Caro et al., 2018). The measured singularity correlations are therefore Radon-type projections of the $T=T_0\tau$6-th moment field on the higher-dimensional configuration space $T=T_0\tau$7.

The inverse results are correspondingly distributional. Under the assumptions that $T=T_0\tau$8 almost surely and is almost surely compactly supported, every moment map $T=T_0\tau$9 is uniquely determined by the exterior correlation data $T_0(x_s,x)=|x-x_s|/v(x_s)$0. For Gaussian random fields, the first two datasets $T_0(x_s,x)=|x-x_s|/v(x_s)$1 and $T_0(x_s,x)=|x-x_s|/v(x_s)$2 determine the full probability distribution. Under the additional exponential moment assumption

$T_0(x_s,x)=|x-x_s|/v(x_s)$3

equality of all exterior singularity-correlation data implies equality in law as $T_0(x_s,x)=|x-x_s|/v(x_s)$4-valued random variables (Caro et al., 2018).

This work is not formulated as an eikonal equation, but its mechanism is closely allied to stochastic eikonal tomography. The relevant observables are first singularities or first arrivals, the measured quantity at infinity is a line integral along a geometric propagation line, and the inversion reconstructs moments and laws of a random medium by Radon inversion on higher-dimensional spaces (Caro et al., 2018). A plausible implication is that stochastic eikonal tomography need not be restricted to posterior estimation for a single slowness field; it can also mean tomography of a hierarchy of random-medium statistics.

6. Limitations, misconceptions, and open directions

The dominant physical limitation is shared by nearly all formulations: the eikonal model is a high-frequency, first-arrival approximation. Later arrivals, amplitude information, finite-frequency sensitivity, and full waveform effects are excluded in the direct probing method, in Bayesian PINNs, in the GP surface-wave formulation, in PINNPStomo, and in the Kīlauea variational-inference workflow (Ito et al., 2023, Gou et al., 2022, Muir, 2023, Song et al., 2024, Gao et al., 31 Jul 2025). Even when multipathing is partially represented by adjoint kernels rather than a single ray, the forward model remains fundamentally first-arrival (Chen et al., 2024).

A second limitation is severe ill-posedness and nonuniqueness. The direct probing paper makes this explicit by exhibiting high-contrast inclusions that yield identical first-arrival boundary data once the contrast is sufficiently large (Ito et al., 2023). This implies that posterior distributions in genuinely stochastic formulations can be broad, multimodal, or insensitive in certain directions, especially when only first arrivals are assimilated. The Kīlauea study similarly emphasizes nonlinearity, ill-posedness, and dependence on initialization and regularization, even though it mitigates these with variational entropy and annealed gradient smoothing (Gao et al., 31 Jul 2025).

A third issue is that uncertainty quantification is often approximate. Mean-field Gaussian VI in BPINNs can be overconfident (Gou et al., 2022). The Kīlauea posterior is a diagonal lognormal variational family, so it cannot represent strong spatial correlations or multimodality (Gao et al., 31 Jul 2025). Gaussian-process eikonal tomography is fully Bayesian at the field-interpolation stage, but hyperparameters are estimated by type-II maximum likelihood rather than fully marginalized (Muir, 2023). In adjacent random-scattering work, full-law recovery requires strong support and exponential-integrability assumptions and idealized measurements of singularity jumps at infinity (Caro et al., 2018).

A frequent misconception is that any randomized eikonal workflow is stochastic tomography. The literature distinguishes more sharply. Random collocation points in a PINN and Adam-based Monte Carlo optimization over a variational family are stochastic optimization devices; they become stochastic tomography only when the output is a posterior or uncertainty-aware model of the medium (Song et al., 2024, Gao et al., 31 Jul 2025). Conversely, stochastic tomography is not exclusively Bayesian: moment recovery and law recovery from singularity correlations constitute tomography on probability distributions without introducing a posterior over parameters (Caro et al., 2018).

Several deterministic advances point toward likely extensions. A plausible extension is to combine source-conditioned traveltime networks, adjoint-state gradients, and factored source regularizations with Bayesian likelihoods and priors, especially for joint $T_0(x_s,x)=|x-x_s|/v(x_s)$5–$T_0(x_s,x)=|x-x_s|/v(x_s)$6 inversion and uncertain acquisition geometries, although the deterministic papers stop short of that step (Song et al., 2024, Grubas et al., 2022, Chen et al., 2024). Another plausible extension is to place explicit priors on geometric level-set variables and phasewise slowness fields in discontinuous-media tomography, turning multilayer level-set inversion into a hierarchical stochastic eikonal model (Li et al., 18 Oct 2025). What the current literature already establishes is that stochastic eikonal tomography is not a single algorithmic recipe, but a spectrum of mathematically linked approaches in which first-arrival geometry, probabilistic modeling, and integral geometry intersect.