Bayesian seismic waveform inversion: Parameter estimation and uncertainty analysis

The goal of geophysical inversion is to make quantitative inferences about the Earth from remote observations. Because the observations are finite in number and subject to uncertainty, these inferences are inherently probabilistic. A key step is to define what it means for an Earth model to fit the data. This requires estimation of the uncertainties in the data, both those due to random noise and those due to theoretical errors. But the set of models that fit the data usually contains unrealistic models; i.e., models that violate our a priori prejudices, other data, or theoretical considerations. One strategy for eliminating such unreasonable models is to define an a priori probability density on the space of models, then use Bayes theorem to combine this probability with the data misfit function into a final a posteriori probability density reflecting both data fit and model reasonableness. We show here a case study of the application of the Bayesian strategy to inversion of surface seismic field data. Assuming that all uncertainties can be described by multidimensional Gaussian probability densities, we incorporate into the calculation information about ambient noise, discretization errors, theoretical errors, and a priori information about the set of layered Earth models derived from in situ petrophysical measurements. The result is a probability density on the space of models that takes into account all of this information. Inferences on model parameters can be derived by integration of this function. We begin by estimating the parameters of the Gaussian probability densities assumed to describe the data and model uncertainties. These are combined via Bayes theorem. The a posteriori probability is then optimized via a nonlinear conjugate gradient procedure to find the maximum a posteriori model. Uncertainty analysis is performed by making a Gaussian approximation of the a posteriori distribution about this peak model. We present the results of this analysis in three different forms: the maximum a posteriori model bracketed by one standard deviation error bars, pseudo-random simulations of the a posteriori probability (showing the range of typical subsurface models), and marginals of this probability at selected depths in the subsurface. The models we compute are consistent both with the surface seismic data and the borehole measurements, even though the latter are well below the resolution of the former. We also contrast the Bayesian maximum a posteriori model with the Occam model, which is the smoothest model that fits the surface seismic data alone.