Travel-time tomography

Part 7 — Travel-time, surface-wave, and joint inversion

Learning objectives

Set up the first-arrival travel-time tomography inverse problem
Run a SIRT (Simultaneous Iterative Reconstruction Technique) inversion end-to-end on a crosshole geometry
Recognise straight-ray resolution limits: along-ray smearing of recovered anomalies
Understand the role of Tikhonov regularisation in suppressing noise
Map this back to the modern PINN-FATT framework

The first three sections of Part 7 solved the FORWARD eikonal problem: given a velocity model $c(x)$ , compute travel times. We now confront the INVERSE problem: given observed travel times at receivers, recover $c(x)$ . This is travel-time tomography — one of the oldest and most widely-deployed seismic-inverse methods, dating to Aki-Lee 1976 for earthquake-source studies and Dines-Lytle 1979 for industry surveys.

The forward and inverse problems

For each source $x_s$ , the FORWARD problem is the eikonal we already know:

|\nabla T_s(x)|^2 = s(x)^2 = 1/c(x)^2, \quad T_s(x_s) = 0 .

The data is a matrix of picked first-break times $t_{\mathrm{obs}}(x_r; x_s)$ at receiver-source pairs. The INVERSE problem minimises

J(s) = \sum_{s, r} \bigl( T_{\mathrm{pred}}(x_r; x_s) - t_{\mathrm{obs}}(x_r; x_s) \bigr)^2 + \lambda \, \|\nabla s\|_2^2 ,

where the second term is the Tikhonov smoothness regulariser preventing checkerboard artefacts. The challenge is non-uniqueness: many slowness fields can fit the same finite set of travel times, especially with limited source-receiver coverage.

SIRT — the iterative classical baseline

SIRT (Gilbert 1972) is the workhorse of seismic tomography codes. It iteratively backprojects timing residuals along ray paths:

Trace straight rays from each source to each receiver through the current slowness field $s_k(x)$ .
Predict travel times $t_{\mathrm{pred}} = \int_{\mathrm{ray}} s_k(x),dl$ .
Backproject residual $\Delta t = t_{\mathrm{obs}} - t_{\mathrm{pred}}$ uniformly along the ray:

\Delta s_{i} = \frac{\sum_{\mathrm{rays} \, j \, \mathrm{through} \, i} \Delta t_j \cdot w_{ij} / L_j}{\sum_{\mathrm{rays} \, j \, \mathrm{through} \, i} w_{ij}}

where $w_{ij}$ is the path length of ray $j$ in cell $i$ and $L_j$ is the total length of ray $j$ .

This is Kaczmarz-style projection: each ray contributes its mismatch evenly to all cells it visits. The key SIRT difference from algebraic ART is that all rays update simultaneously (averaged) rather than one at a time — this gives much better noise robustness at the cost of slower convergence.

After each backprojection step we apply a 5-point Laplacian smoothing as Tikhonov regularization:

s^{(\mathrm{smoothed})}_i = (1 - \alpha) s_i + \alpha \, \overline{s}_i ,

where $\overline{s}_i$ is the average of the four neighbour cells. We use $\alpha = 0.18$ in the widget below — light smoothing that suppresses checkerboard noise without over-blurring real structure.

Try it: SIRT in 30 iterations

The widget runs:

True model: 2.5 km/s background with a Gaussian low-velocity anomaly ( $-0.7$ km/s amplitude, σ_x = 0.30, σ_z = 0.20) centred at $(1.0, 0.5)$ km. Crosshole geometry: 6 sources in the LEFT borehole (yellow dots at x ≈ 0), 11 receivers in the RIGHT borehole (cyan squares at x ≈ 2 km). Rays travel horizontally across the domain at all depths — this is the canonical configuration that lets surface-distant sensors illuminate buried structure. Surface-only geometries with straight rays cannot resolve depth because all rays stay near the surface and never visit the buried anomaly. Real-world surface tomography works because Fermat's principle bends rays into faster layers, providing depth sensitivity that straight-ray algorithms cannot exploit.
Forward modelling: FSM solver from §7.1 generates $t_{\mathrm{obs}}$ for all source-receiver pairs.
Inversion: SIRT starts from a homogeneous 2.5 km/s prior (no anomaly knowledge). 30 iterations with $\alpha = 0.18$ . Updates the c-field inline every 3 iterations so you can watch the anomaly emerge.

Four panels: true model with source-receiver geometry, recovered model evolving over iterations, the difference map (c_inv − c_true on a diverging colormap), and RMS travel-time residual on a log-y scale.

Expected behaviour:

RMS residual drops 1-2 orders of magnitude (from ~30-50 ms initial to ~1-3 ms after 30 iterations).
Recovered anomaly amplitude reaches 60-80% of the true -0.7 km/s — never 100% because SIRT regularisation always smooths slightly. Increase the number of sources and receivers, or relax $\alpha$ , to push higher.
HORIZONTAL SMEARING along the ray-cross direction: the anomaly appears more elongated horizontally than the true Gaussian. This is because crosshole rays cross the anomaly at low angle and the SIRT backprojection spreads the residual along the entire ray, not just at the anomaly location. Adding source-receiver pairs with more diverse incidence angles (e.g., a third borehole, or surface-to-borehole "VSP" geometry) sharpens the recovered anomaly.

Resolution limits: why straight-ray geometries always smear

The sensitivity of a travel-time observation $t(x_r; x_s)$ to a slowness perturbation $\delta s(x)$ is concentrated along the ray path connecting $x_s$ and $x_r$ . Backprojecting residuals along that ray spreads the correction over its ENTIRE LENGTH, not just at the spot where the anomaly actually sits. So SIRT cannot localise an anomaly better than the FAN of ray angles passing through it. With crosshole geometry, that fan is roughly the angle subtended by the source array as seen from the anomaly — typically 30-60°.

The fix in production tomography: use diverse geometry (multi-borehole + surface), use FREQUENCY-DEPENDENT sensitivity kernels (banana-doughnut wavefronts of finite-frequency tomography, Dahlen et al 2000), or invert SURFACE-WAVE DISPERSION which has frequency-dependent depth sensitivity (§7.6). Single-band first-arrival tomography is fundamentally limited by ray-path geometry.

The PINN-FATT extension

Modern travel-time tomography uses neural networks for both the forward and inverse pieces. The PINN-FATT formulation (Smith et al 2022, bin Waheed-Alkhalifah-Haghighat 2022) parameterises:

The slowness field $s_{\mathrm{NN}}(x; \xi)$ as a small MLP with output passed through a softplus or sigmoid to ensure positivity.
The travel time $T_{\mathrm{NN}}(x; x_s; \theta)$ as a factored EikoNet (§7.3) conditioned on source position via input concatenation.

Both networks train JOINTLY against

\mathcal{L} = \underbrace{\sum_{s,r} (T_{\mathrm{NN}}(x_r; x_s) - t_{\mathrm{obs}})^2}_{\text{data}} + \underbrace{\frac{1}{N_c} \sum_k (|\nabla T_{\mathrm{NN}}|^2 - s_{\mathrm{NN}}^2)^2}_{\text{eikonal residual}} + \underbrace{\lambda \|\nabla s_{\mathrm{NN}}\|^2}_{\text{smoothness}} .

Three advantages over classical SIRT:

Continuous, smooth $s(x)$ via the implicit network parameterisation — no grid discretisation artefacts.
Bent-ray forward modelling via EikoNet — straight-ray SIRT systematically underestimates travel-time anomalies in regions with strong velocity variation.
Source-conditioning — train ONCE on a base velocity model, then evaluate at any source for free. Useful in microseismic monitoring where new sources arrive constantly (§7.5).

The downside: PINN-FATT training is slower (~minutes on GPU vs SIRT's seconds) and requires careful loss-weight tuning. The factored eikonal of §7.3 is essential — without it, the eikonal-residual loss couples too strongly with the data term and the joint training stalls. We do not implement PINN-FATT in the widget here — at the in-browser scale, the SIRT widget already demonstrates the inversion principle, and the additional cost of joint network training does not pay off pedagogically. Production codes (HypoSVI, Smith 2022) implement the full pipeline.

Operational notes

Three lessons that apply to ANY travel-time tomography (classical or PINN):

Initial model matters more than algorithm. SIRT, conjugate gradient, LSQR, PINN-FATT — all fail if the homogeneous initial $c$ is too far from truth. In production, start from a 1-D depth-trend model derived from sonic logs or refraction arrivals.
Source-receiver geometry sets resolution. No algorithm can recover what the rays don't illuminate. Run a CHECKERBOARD TEST first: invert synthetic data with checkerboard anomalies of known size to see what wavelengths are recoverable in the actual survey geometry.
Pick quality dominates everything. Travel-time picks have ±5-10 ms uncertainty in good data, ±50+ ms in noisy data. Tomography RMS residuals below pick-error reflect overfitting to noise. Cross-validate on held-out picks.

References

Gilbert, P. (1972). Iterative methods for the three-dimensional reconstruction of an object from projections. J. Theor. Biol. 36, 105–117. The original SIRT paper.
Aki, K., Lee, W.H.K. (1976). Determination of three-dimensional velocity anomalies under a seismic array using first P arrival times from local earthquakes. J. Geophys. Res. 81, 4381–4399. First crustal-scale travel-time tomography.
Dines, K.A., Lytle, R.J. (1979). Computerized geophysical tomography. Proc. IEEE 67(7), 1065–1073. Industry tomography review.
Smith, J.D., Ross, Z.E., Azizzadenesheli, K., Muir, J.B. (2022). HypoSVI: Hypocentre inversion with Stein variational inference and physics informed neural networks. Geophys. J. Int. 228, 698–710. The PINN-FATT framework with uncertainty quantification.
bin Waheed, U., Alkhalifah, T., Haghighat, E., Song, C., Virieux, J. (2022). PINNTomo: Seismic tomography using physics-informed neural networks. arXiv:2104.01588. Joint slowness + travel-time PINN tomography.