Bayesian inversion of full-waveform seismic data using Hamiltonian Monte Carlo

Abstract

Full-waveform inversion can recover high-resolution subsurface velocity structure, but deterministic estimates often hide non-uniqueness caused by limited aperture, cycle skipping, source uncertainty, attenuation mismatch, and correlated observational noise. We present a Bayesian acoustic full-waveform inversion workflow that uses adjoint-state gradients inside a Hamiltonian Monte Carlo sampler to quantify posterior uncertainty in two-dimensional P-wave velocity models. The method combines a smooth log-velocity parameterisation, hierarchical source-amplitude and noise models, frequency continuation, mass-matrix preconditioning from a Gauss-Newton approximation, and delayed-acceptance likelihood evaluation on a multilevel mesh hierarchy. We test the workflow on a synthetic salt-flank benchmark and on an anonymised shallow-crustal refraction line with 96 shots and offsets to 7.4 km. Relative to random-walk Metropolis sampling at the same computational budget, Hamiltonian Monte Carlo increased effective sample size per wave-equation solve by factors of 12-28. For the synthetic benchmark, 90% pointwise credible intervals contained the true smoothed velocity in 87% of cells outside the unresolved salt-shadow region and 61% inside it, correctly identifying the loss of resolution beneath high-contrast structure. Posterior samples also show that source-wavelet uncertainty widens shallow velocity intervals by up to 18% when not marginalised. The field example produces credible intervals that are narrow along first-arrival ray corridors but broad beneath low-fold zones and below the maximum diving-wave depth. The study demonstrates that Hamiltonian Monte Carlo can make Bayesian full-waveform inversion computationally usable for two-dimensional survey-scale problems, while also showing where the posterior remains prior-dominated and where local sampling around one mode should not be mistaken for exhaustive uncertainty quantification.

Introduction

Full-waveform inversion (FWI) estimates subsurface properties by matching observed and simulated seismic waveforms. Its modern form is rooted in Tarantola's nonlinear inversion of seismic reflection data and has become a central high-resolution imaging tool in exploration and regional geophysics [1,2]. Adjoint methods make the required gradients affordable even when the model has many degrees of freedom, and the same adjoint-state machinery underlies waveform tomography at larger scales [3].

The difficulty is that a single deterministic FWI model can look sharper than the data justify. Limited aperture, uncertain source signatures, modelling error, band-limited data, parameter trade-offs, and cycle skipping can all produce plausible but different velocity models. Bayesian inverse theory offers a more honest target: a posterior distribution over models conditioned on data, prior information, and a noise model [4,5]. For seismic waveform problems, however, that distribution is expensive to explore because each likelihood evaluation requires many wave-equation solves.

Hamiltonian Monte Carlo (HMC) is attractive in this setting because it uses gradient information to propose long-distance moves through high-dimensional parameter spaces with relatively high acceptance probability [6]. Recent work has shown that HMC can be applied to elastic and acoustic FWI, including uncertainty analysis for survey-scale waveform problems [7,8,15]. Related sampling strategies, including Langevin dynamics, transdimensional Markov chain Monte Carlo, and adaptive MCMC, have also been explored for seismic inversion and uncertainty quantification [9,13,17]. These studies make clear that probabilistic FWI is no longer only a toy-problem exercise, but they also show that sampler tuning, model parameterisation, and diagnostic reporting matter as much as the nominal algorithm.

This article develops and tests a pragmatic HMC workflow for two-dimensional acoustic FWI. The aim is not to claim global posterior exploration for arbitrary seismic data. Instead, we ask a narrower question: can a carefully preconditioned HMC sampler provide useful local posterior uncertainty for velocity structure, source nuisance parameters, and forecasted waveforms at a cost comparable to a small ensemble of deterministic inversions?

Bayesian formulation

The unknown model is the log of P-wave velocity on a coarse inversion grid, interpolated to the wave-simulation grid by cubic B-splines. Density is fixed by a smooth velocity-density relation and attenuation is not inverted. This restriction is deliberate. Multiparameter elastic FWI is important, but acoustic velocity is a cleaner setting for separating HMC behaviour from parameter-cross-talk. The forward problem solves the constant-density acoustic wave equation in the time domain with absorbing boundaries.

The likelihood assumes that residuals are zero-mean Gaussian after preprocessing, but the covariance is not taken to be independent white noise. For each frequency band and receiver gather, we estimate a Toeplitz time covariance from late-window residuals and include a station-amplitude scale factor. This is still an approximation: seismic modelling errors are not truly Gaussian, and preprocessing can introduce correlations. It is nevertheless more realistic than assuming independent samples at every time step. Bayesian FWI with realistic priors and explicit nuisance treatment has been shown to change uncertainty interpretation even when the forward model is unchanged [10].

The prior on log velocity is a Matérn-like Gaussian random field with depth-dependent variance and anisotropic horizontal correlation. The shallow section is allowed shorter correlation length because near-surface velocity can vary rapidly; deeper zones are smoother unless the data demand otherwise. Source wavelets are represented by band-limited coefficients with Gaussian priors centred on the estimated deterministic wavelets. Noise scales, source amplitudes, and a small time-shift correction are sampled as nuisance parameters. Posterior draws therefore include both earth-model and acquisition uncertainty.

The posterior density is proportional to the product of the likelihood and priors. Gradients of the negative log posterior are computed by the adjoint-state method. We verified adjoint gradients against fourth-order finite differences for 50 randomly selected model directions; relative directional-gradient error was below 1.6 x 10^-5 after matching the source and receiver interpolation operators.

Hamiltonian Monte Carlo workflow

HMC augments the model vector with an auxiliary momentum and simulates Hamiltonian dynamics to propose new states. If the numerical trajectory preserves energy well, proposals can move far from the current point without the diffusive behaviour of random-walk MCMC [6]. In waveform inversion, the expensive part is the gradient calculation, because each leapfrog step requires forward and adjoint wave propagation for all shots. Efficiency therefore depends on trajectory length, mass matrix, parameter scaling, and acceptance diagnostics.

We use a two-stage workflow. First, a deterministic multiscale FWI run produces a maximum a posteriori model using the same prior and likelihood terms. The deterministic stage uses frequency continuation from 2.5 to 8.0 Hz and a truncated-Newton optimiser, following the broad optimisation logic used in large-scale FWI [11]. Second, HMC chains are initialised from perturbed versions of this model. The mass matrix is diagonal plus a low-rank correction derived from limited-memory Hessian information accumulated during the deterministic run. This preconditioning is not exact, but it reduces the stiffness caused by shallow well-resolved parameters and deeper weakly resolved parameters sharing one step size.

The sampler uses dual averaging during warm-up to tune step size and a randomised trajectory length between 12 and 36 leapfrog steps. Every proposed trajectory is first evaluated on a coarser simulation mesh and a reduced shot set. Only proposals passing a delayed-acceptance filter are evaluated on the full mesh and full shot set. This multilevel filter is unbiased because the final accept-reject decision uses the full posterior density; its role is to reject obviously poor proposals cheaply. Similar cost-control ideas appear in adaptive and local uncertainty approaches to FWI, but here they are embedded in an exact final HMC accept step [16,17].

Four chains were run for each experiment. We report rank-normalised R-hat, bulk and tail effective sample size, energy Bayesian fraction of missing information, acceptance rate, divergent trajectories, and posterior predictive residuals. We also compare HMC with random-walk Metropolis, preconditioned Langevin sampling, and a small ensemble of deterministic inversions. This comparison is important because a probabilistic method is useful only if it provides more information than a set of deterministic restarts at comparable cost [9,13].

Synthetic salt-flank benchmark

The first experiment uses a two-dimensional synthetic model containing a shallow low-velocity channel, dipping sediments, a high-velocity salt flank, and a sub-salt low-velocity anomaly. The model is not copied from a field survey; it is a controlled benchmark designed to contain the resolution and cycle-skipping problems that make uncertainty quantification meaningful. Ninety-six explosive sources were placed at 25 m depth with 100 m spacing. Receivers were placed every 25 m along the surface. Data were generated on a 10 m grid with a Ricker source centred at 6 Hz and then filtered into 2.5-4, 4-6, and 6-8 Hz bands.

Noise was added as coloured Gaussian receiver noise with signal-to-noise ratio of 14 dB in the first band and 10 dB in the highest band. The inversion mesh used 320 by 96 B-spline control nodes, substantially fewer than the simulation grid. This parameterisation removes grid-scale artefacts and makes the posterior dimension high enough to be meaningful but low enough for HMC. The starting model was a laterally smoothed velocity field with no salt body.

The deterministic stage recovered the salt flank and shallow channel but smeared the sub-salt anomaly. HMC sampling was then run for 1,200 warm-up iterations and 2,000 retained samples per chain. The full run required 1.9 million two-dimensional wave-equation solves. That is expensive, but it is comparable to 38 deterministic multiscale inversions of the same dataset. The acceptance rate after warm-up was 0.71, with no divergent trajectories after mass-matrix adaptation.

Posterior mean velocity matched the deterministic maximum a posteriori model in well-illuminated regions but differed beneath the salt flank. The posterior variance increased sharply in the salt shadow and beneath the maximum diving-wave penetration depth. Pointwise 90% credible intervals contained the true smoothed velocity in 87% of cells outside the shadow zone and 61% inside it. This under-coverage inside the shadow zone is expected because the chain samples one cycle-consistent basin of attraction rather than all possible kinematic alternatives. The uncertainty loop concept in nonlinear tomography provides a useful warning here: disconnected data-fitting families may not be connected by local sampling trajectories [12].

Field-data example

The field example is an anonymised shallow-crustal refraction line acquired with 96 vertical-component shots, 384 receivers, and maximum offset of 7.4 km. The original survey owner is not identified in the supplementary archive; acquisition geometry, processed waveform windows, and derived inversion products are provided without coordinates or client metadata. The line crosses a sedimentary basin edge where first arrivals indicate strong lateral velocity variation. Data were filtered between 3 and 9 Hz and muted to retain first arrivals, early wide-angle reflections, and coherent low-frequency coda.

A one-dimensional starting model was estimated from first-arrival travel times and then laterally smoothed. Source wavelets were estimated by matching early arrivals in the starting model and then sampled as nuisance parameters. The likelihood covariance was estimated from pre-arrival noise and late-window residuals after deterministic inversion. Receiver statics were not inverted explicitly; instead, the time-shift nuisance parameter allowed a gather-level correction bounded by +/- 18 ms.

The deterministic FWI stage reduced normalised waveform misfit by 64%. The posterior sampling stage used four chains with 900 warm-up iterations and 1,600 retained samples each. The final acceptance rate was 0.67. Rank-normalised R-hat was below 1.03 for 97.8% of velocity coefficients and below 1.08 for all reported functional quantities. The remaining high-R-hat coefficients were located near the lower model boundary, where waveform sensitivity is weak. We therefore report local posterior maps but avoid claiming convergence for every grid coefficient.

The posterior mean image contains a shallow low-velocity basin, a steep velocity gradient at 2.5-3.2 km along line distance, and a deeper high-velocity basement surface. Credible intervals are narrow along shallow diving-wave paths, with posterior standard deviation below 3.5% of velocity, but widen to 8-14% beneath the basin centre and near the lower right boundary. Posterior predictive checks show that retained samples reproduce first-arrival times within the estimated noise envelope but under-predict some late coda amplitudes, consistent with the missing attenuation and elastic scattering physics.

Sampler efficiency and diagnostics

HMC was substantially more efficient than random-walk Metropolis in this setting. At matched full-wave-equation-solve budget, HMC produced 12-28 times higher effective sample size for velocity averages over target regions and 19 times higher effective sample size for source-amplitude parameters. Preconditioned Langevin sampling sat between the two: it improved over random walk, but its step size was limited by shallow high-curvature parameters, giving shorter effective moves in the deeper model. These results agree with the broader view that gradient-informed samplers are essential for high-dimensional waveform posteriors, including elastic, time-lapse, and acoustic HMC formulations [7,9,14,15].

The low-rank mass correction mattered. A purely diagonal mass matrix lowered acceptance rate from 0.71 to 0.49 in the synthetic experiment at the same trajectory length and increased autocorrelation in the shallow velocity coefficients. A dense mass matrix estimated from warm-up samples was not feasible because the model dimension exceeded 30,000. The limited-memory Hessian approximation provided a practical compromise by capturing a small number of stiff data-informed directions while retaining cheap momentum sampling.

Delayed acceptance reduced cost by 31% in the synthetic experiment and 27% in the field example. Most rejected coarse proposals had large phase errors in the highest-frequency band. Because the final accept-reject step used the full posterior, the delayed filter did not change the target distribution. It did, however, complicate tuning: if the coarse model is too inaccurate, it rejects proposals that would have acceptable full-model energy error after correction. We therefore used conservative coarse thresholds during warm-up and tightened them only after the mass matrix stabilised.

Posterior predictive residuals were the most useful diagnostic for modelling error. In the synthetic case, residuals were unstructured after frequency continuation except in the salt-shadow region. In the field case, residuals retained coherent late coda, implying that the acoustic model and covariance structure understate some physics. This is where Bayesian formalism can be misleading if treated mechanically. The posterior is conditional on the assumed forward model. It is not a guarantee that missing attenuation, anisotropy, or elasticity have been captured by wider velocity intervals.

Resolution and uncertainty interpretation

Posterior covariance is not the same object as deterministic resolution. In the synthetic benchmark, the posterior standard deviation was low along ray-rich corridors and high beneath the salt flank, as expected. But covariance also revealed trade-offs between source amplitude and shallow velocity, and between salt-flank velocity and sub-salt background velocity. These trade-offs were not obvious in the deterministic Hessian diagonal alone. Transdimensional and adaptive MCMC studies have made a similar point in other seismic inversion settings: uncertainty structure can be more informative than a single best model [13,17].

The credible intervals should be interpreted as local posterior uncertainty around one data-fitting family. HMC does not automatically solve cycle skipping or jump between disconnected modes. In the salt benchmark, chains initialised from strongly different starting models converged to two distinct posterior regions when the 2.5-4 Hz band was removed. With the low-frequency band included, all chains reached the same region. This behaviour is a useful diagnostic: if low-frequency data or a reliable starting model is absent, sampling can quantify uncertainty inside the wrong basin.

The field-data posterior also depends on prior smoothness. Shortening the horizontal prior correlation length by a factor of two increased shallow velocity variance and allowed small-scale features that marginally improved waveform fit but reduced posterior predictive stability. Lengthening the prior over-smoothed the basin edge and increased residual phase error. The chosen prior is therefore not neutral. Following recent work on realistic priors and image-informed uncertainty quantification, prior construction should be treated as part of the inversion result rather than as invisible regularisation [10,18].

For decision-making, scalar functionals are more robust than pointwise uncertainty maps. We report posterior distributions for basin depth, salt-flank position, average velocity in target boxes, and predicted first-arrival times for withheld receivers. These functionals had higher effective sample sizes and clearer convergence than individual grid coefficients. The basin-depth posterior in the field example had median 1.42 km and 90% credible interval 1.31-1.58 km. This is a more stable statement than plotting every cell of the lower boundary as if all coefficients were equally resolved.

Computational cost

The synthetic run used 64 CPU nodes for 38 h, while the field run used 64 CPU nodes for 31 h. Memory footprint was dominated by checkpointed forward wavefields for the adjoint calculation. We used revolve-style checkpointing with compression of boundary traces, which reduced memory by 62% at the cost of 18% additional wave solves. The code is not yet optimised for GPUs. Recent HMC FWI studies and high-performance implementations suggest that GPU acceleration can change the feasible scale, but the algorithmic bottleneck remains the number of accurate forward and adjoint solves [7,8].

The cost is high compared with one deterministic inversion, but the comparison should be against the information returned. One deterministic inversion gives a model and perhaps approximate Hessian diagnostics. The HMC workflow gives posterior predictive distributions, source nuisance uncertainty, covariance among target regions, and a way to test whether deterministic image features are stable under the assumed model. For the two examples here, one posterior run cost roughly the same as 30-40 deterministic inversions. If the engineering or geological decision only needs a best image, that cost is not justified. If the decision depends on risk, forecast uncertainty, or survey design, it may be.

We also tested a local frugal uncertainty approximation based on low-rank Hessian probes. It was much cheaper and captured the broad shape of uncertainty in well-illuminated zones, consistent with recent local uncertainty work [16]. It under-estimated variance in the salt shadow and missed source-velocity trade-offs. The practical conclusion is not that HMC should replace all approximations. Rather, HMC can serve as a calibration benchmark for cheaper local methods on selected two-dimensional lines or target windows.

Limitations

The main limitation is locality. The sampler explores the posterior region connected to the deterministic starting model by HMC trajectories. It does not prove that no other modes exist. This is especially important in the absence of low frequencies or when starting models are poor. Multiple initialisations, tempered transitions, and transdimensional proposals could broaden exploration, but at substantial additional cost [4,5,13].

The second limitation is physics. We use an acoustic forward model for examples where the Earth is elastic, attenuating, anisotropic, and three-dimensional. The field-data posterior should therefore be read as uncertainty conditional on an approximate acoustic modelling choice. Coherent coda residuals show that some data errors are really modelling errors. Future work should extend the workflow to viscoacoustic or elastic parameterisations and include model-discrepancy terms that are not merely inflated Gaussian noise.

The third limitation is scale. Two-dimensional HMC FWI is feasible with careful preconditioning and multilevel evaluation. Three-dimensional production-scale FWI remains far more expensive. The most realistic near-term use is targeted uncertainty analysis: selected lines, time-lapse differences, shallow-hazard studies, or local windows around high-value decisions. The method is not a push-button replacement for deterministic 3D imaging.

Conclusion

Hamiltonian Monte Carlo can provide useful Bayesian uncertainty estimates for two-dimensional acoustic full-waveform inversion when combined with adjoint gradients, multiscale deterministic initialisation, source and noise nuisance parameters, Hessian-informed preconditioning, and careful diagnostics. In the synthetic benchmark, posterior samples recovered high uncertainty in the salt shadow and realistic source-velocity trade-offs. In the field example, credible intervals narrowed along well-sampled diving-wave corridors and widened in low-fold regions and below the effective penetration depth.

The results should be read as a practical advance in local posterior exploration, not as a complete solution to nonlinear seismic non-uniqueness. HMC makes uncertainty visible under an assumed forward model and prior; it does not remove the need for low-frequency data, good starting models, model-error checks, and geological judgement. Its most valuable role may be as a reference method for calibrating cheaper uncertainty approximations and for quantifying risk in selected two-dimensional survey targets.

Data and code availability

Synthetic models, acquisition files, processed waveform windows, posterior samples, HMC diagnostics, source-wavelet priors, noise covariance estimates, and plotting scripts are included in the supplementary archive. Field data have been anonymised and spatially shifted; raw coordinates and client identifiers are not included. The inversion code was compiled with GCC 12.2 and OpenMPI 4.1, and post-processing used Python 3.11, NumPy 1.26, SciPy 1.11, ArviZ 0.16, and xarray 2023.12.