Diffusion Posterior Sampling for Seismic Inversion

Seismic inversion is ill-posed. Many subsurface models fit the data, and the choice of regularisation decides which one you are shown. Smoothness and damping are priors wearing an engineering disguise: they are doing exactly the work a prior does, without being written down as one or being available for inspection.

Diffusion posterior sampling replaces the disguise with the thing itself. A score-based diffusion model is trained as a generative prior over plausible velocity models [1,2]. At inference, the gradient of the data likelihood is injected into each reverse-diffusion step, so the trajectory targets the posterior rather than the prior alone [3]. Run many independent chains and you get an ensemble of data-consistent models whose spread is the uncertainty.

For imaging the lithosphere-asthenosphere boundary from S-receiver functions, where the S-to-P conversion arrives as a precursor to the direct S wave in a window relatively free of crustal multiples [4], this is an appealing fit. The target is spatially structured and high-dimensional, which is where MCMC mixes badly and where a learned prior has the most to offer.

What the method is, in three lines

A score network is trained on an ensemble of geologically plausible velocity models, learning the gradient of the log-density of noised versions of them. Sampling starts from noise and integrates the reverse stochastic differential equation toward a clean model, moving along the learned score at each step [1].

Data consistency is what makes it an inversion rather than a generator. At each reverse step, the gradient of the data likelihood, evaluated through a forward operator, is added to the update [3]. The trajectory is pulled toward models that both look like the training ensemble and explain the recorded waveforms.

Relative to a deterministic inversion this trades an opaque smoothness penalty for an explicit, learned, sample-able prior. Relative to MCMC it produces independent samples without random-walk mixing pathologies in a multimodal, high-dimensional space.

The ordering that costs money

stop the reverse trajectory after this many of 100 steps20 steps, 80% of the cost saved40 steps, 60% of the cost saved60 steps, 40% of the cost saved80 steps, 20% of the cost saved100 steps, 0% of the cost saved

Stopping at step 60 saves 40% of the sampling cost. The picture has looked finished since step 58, and the boundary depth can still move 4.03 km. About 1.95 km of that is genuine posterior spread; the rest is sampler error, and nothing in the ensemble distinguishes the two.

How much of each feature is decided, by stopping point. A feature counts as settled at 90 percent.

feature	20 steps	40 steps	60 steps	80 steps	100 steps	still moving here
gross crustal structure	81%	96%	99%	100%	100%	0.01
layer count and ordering	55%	80%	91%	96%	98%	0.09
boundary depth	30%	52%	66%	77%	84%	4.03 km
boundary sharpness	23%	40%	54%	64%	72%	0.46

One property of the sampler decides whether a run is worth trusting, and it is not the one the cost objection points at.

The reverse trajectory does not resolve the model uniformly. Coarse structure is decided early and fine structure late, which is the ordering our note describes and which the exhibit puts numbers on. By step 60 of 100, gross crustal structure is 99% decided and layer count and ordering 91%. The picture has looked finished since step 58.

The boundary depth at that point is 66% decided, and can still move about 4 km. Boundary sharpness is barely past half.

Every reverse step costs a network evaluation, and sampling cost is the method's main practical objection: many evaluations per sample, many samples per posterior. So there is a standing incentive to stop once the image looks right, and the image looks right early. Stopping at step 60 saves 40% of the sampling cost and leaves the quantity anyone commissioned the study for still in motion.

What makes this worse than an ordinary speed-accuracy trade is that the damage is invisible in the output. Run the full trajectory and roughly 2 km of movement remains across independent chains; that residual is the posterior spread, and it is the honest uncertainty you were trying to measure. Stop early and the extra movement is sampler error sitting on top of it. In the ensemble the two look identical: a spread of boundary depths, distributed plausibly, that you will read as uncertainty and report as uncertainty. The under-converged run does not announce itself by looking rough. It announces itself by looking uncertain, which is the one appearance you were expecting anyway.

The practical discipline follows from that. Convergence has to be checked on the fine-scale quantity that matters, not on the appearance of the field, and a sampler budget has to be justified against the feature you are measuring rather than against the point at which the image stops changing visibly.

What the prior does not guarantee

Three limits deserve to be stated where a reader will see them rather than in an appendix.

The learned prior can produce structure that is plausible under the training ensemble and unsupported by the data, particularly where the S-receiver function is uninformative. The data-consistency term reduces this and does not remove it, and a biased training ensemble biases every sample drawn from it.

Posterior sampling with diffusion priors is approximate. The likelihood-gradient injection is a guidance heuristic, and exact posterior correctness is not guaranteed [3]. The ensemble is a good estimate of the posterior, not a draw from it.

And the forward operator caps everything. Ray-theoretical approximations bound the accuracy of the likelihood no matter how good the prior is, which is a limit this method shares with every other one in the series.

None of that argues against the approach. It argues for reporting the prior ensemble, the sampler budget and the convergence check alongside the result, because all three are load-bearing and none of them is visible in the picture.

Key takeaways

1. 1A deterministic inversion's smoothness penalty is a prior in disguise; a score-based diffusion model makes it explicit, learned and sample-able.
2. 2Injecting the data-likelihood gradient at each reverse step turns generation into posterior sampling, and the spread of independent chains is the uncertainty.
3. 3The trajectory resolves coarse structure early and fine structure late: by step 60 of 100 the layering is 91% decided while the boundary depth is 66% decided.
4. 4The picture looks finished from step 58, so the standing incentive to stop early removes precisely the fine-scale part that carries the answer.
5. 5An under-converged run does not look rough, it looks uncertain: the extra movement is indistinguishable in the ensemble from genuine posterior spread.
6. 6Report the prior ensemble, the sampler budget and the convergence check with the result, because all three are load-bearing and none is visible in the output.

References

1. Y. Song, J. Sohl-Dickstein, D. P. Kingma, A. Kumar, S. Ermon, B. Poole. Score-Based Generative Modeling through Stochastic Differential Equations. ICLR 2021. arXiv:2011.13456
2. J. Ho, A. Jain, P. Abbeel. Denoising Diffusion Probabilistic Models. NeurIPS 2020. arXiv:2006.11239
3. H. Chung, J. Kim, M. T. Mccann, M. L. Klasky, J. C. Ye. Diffusion Posterior Sampling for General Noisy Inverse Problems. ICLR 2023. arXiv:2209.14687
4. X. Yuan, R. Kind, X. Li, R. Wang. The S receiver functions: synthetics and data example. Geophysical Journal International, 165, 555 to 564, 2006. doi:10.1111/j.1365-246X.2006.02885.x
5. The low-dimensional, amortized complement to this approach: Amortized Bayesian inversion of crustal structure.