Bayesian seismic tomography using normalizing flows

An important feature of Bayesian statistics is the opportunity to do sequential inference: the posterior distribution obtained after seeing a dataset can be used as prior for a second inference. However, when Monte Carlo sampling methods are used for inference, we only have a set of samples from the posterior distribution. To do sequential inference, we then either have to evaluate the second posterior at only these locations and reweight the samples accordingly, or we can estimate a functional description of the posterior probability distribution from the samples and use that as prior for the second inference. Here, we investigated to what extent we can obtain an accurate joint posterior from two datasets if the inference is done sequentially rather than jointly, under the condition that each inference step is done using Monte Carlo sampling. To test this, we evaluated the accuracy of kernel density estimates, Gaussian mixtures, mixtures of factor analyzers, vine copulas and Gaussian processes in approximating posterior distributions, and then tested whether these approximations can be used in sequential inference. In low dimensionality, Gaussian processes are more accurate, whereas in higher dimensionality Gaussian mixtures, mixtures of factor analyzers or vine copulas perform better. In our test cases of sequential inference, using posterior approximations gives more accurate results than direct sample reweighting, but joint inference is still preferable over sequential inference whenever possible. Since the performance is case-specific, we provide an R package mvdens with a unified interface for the density approximation methods.

In-line inspection (ILI) data is commonly used in corrosion growth models (CGMs) to predict the corrosion growth in energy pipelines. A hierarchical stochastic corrosion growth model is considered in this paper which considers the variations in the corrosion growth, both spatially and temporally, the inherent measurement error of the ILI tools as well as the model uncertainties. These uncertainties are represented as unknown model variables and are often inferred using a Bayesian method [1], [2] and samples of the unknown parameters’ posterior probability density functions (PDFs) are obtained using Markov Chain Monte Carlo (MCMC) sampling techniques [3]. ILIs can result in massive data sets. In order for MCMC-based inference techniques to yield reasonably accurate results, many samples (approaching infinity) are required. This fact in addition to the massive data sets exponentially increases the scale of the inference problem from an attainable solution to a potentially impossible one that is limited by today’s computing power. For this reason, MCMC-based inference techniques can become inefficient in the cases where ILI datasets are large. The objective is to propose variational inference (VI) as an alternative to MCMC to determine a Bayesian solution for the unknown parameters in complex stochastic CGMs. VI produces approximations of the posterior PDFs by treating the inference as an optimization problem. Variational inference emerged from machine learning for Bayesian inference of large data sets; therefore, it is an appropriate tool to use in the analysis of mass pipeline inspection data[4]–[7]. This paper introduces VI to solve the inference problem and provide a solution for a hierarchical stochastic CGM to describe the defect-specific corrosion growth experienced in pipelines based on excessively large ILI datasets. To gauge the accuracy of the VI implementation in the model, the results are compared to a set of values generated using a stochastic gamma process that represents the corrosion growth process experienced by the pipe.

Geoscientists use observed data to estimate properties of the Earth's interior. This often requires non‐linear inverse problems to be solved and uncertainties to be estimated. Bayesian inference solves inverse problems under a probabilistic framework, in which uncertainty is represented by a so‐called posterior probability distribution. Recently, variational inference has emerged as an efficient method to estimate Bayesian solutions. By seeking the closest approximation to the posterior distribution within any chosen family of distributions, variational inference yields a fully probabilistic solution. It is important to define expressive variational families so that the posterior distribution can be represented accurately. We introduce boosting variational inference (BVI) as a computationally efficient means to construct a flexible approximating family comprising all possible finite mixtures of simpler component distributions. We use Gaussian mixture components due to their fully parametric nature and the ease with which they can be optimized. We apply BVI to seismic travel time tomography and full waveform inversion, comparing its performance with other methods of solution. The results demonstrate that BVI achieves reasonable efficiency and accuracy while enabling the construction of a fully analytic expression for the posterior distribution. Samples that represent major components of uncertainty in the solution can be obtained analytically from each mixture component. We demonstrate that these samples can be used to solve an interrogation problem: to assess the size of a subsurface target structure. To the best of our knowledge, this is the first method in geophysics that provides both analytic and reasonably accurate probabilistic solutions to fully non‐linear, high‐dimensional Bayesian full waveform inversion problems.

In order to predict future performance of subsurface fluid reservoirs under possible operating scenarios, a dynamic, porous‐medium flow simulation model must be tuned to include representative properties of the reservoir. Estimating subsurface reservoir properties given remotely sensed or borehole‐based observations typically involves finding the solution to a challenging inverse problem. We compare Monte Carlo random sampling to variational inference methods which use optimization to constrain parametrized uncertainties in nonlinear Bayesian inversions. We use them to estimate the posterior probability distribution of reservoir permeability given fluid pressure and seismic measurements. The methods include automatic differentiation variational inference (ADVI), Stein variational gradient descent (SVGD), and a Monte Carlo method called stochastic SVGD (sSVGD), all of which we benchmark against results from Metropolis‐Hastings McMC. We also test an ADVI variant called physically structured variational inference (PSVI): in our implementation this method estimates only spatially‐local correlations between model parameters based on the intuition that such correlations are strong in remote sensing problems in which data only inform about spatial‐averages of local dynamics. We apply the methods to two‐ and three‐dimensional inverse problems of carbon dioxide storage, inspired by the Endurance field, located in the UK North Sea. Results show that PSVI achieves a good balance between mean‐field ADVI and full‐rank ADVI in terms of accuracy of the posterior approximation and computational efficiency. SVGD and sSVGD offer more accurate approximations of the target posterior distribution, but at far higher computational cost. Between them, sSVGD outperforms SVGD, exhibiting better computational efficiency and mitigating the problems of mode collapse and spurious correlations.

In this study, the relationship between nonlinear model properties and inverse problem solutions is analyzed using a numerical technique based on the inverse problem theory formulated by Mosegaard and Tarantola. According to this theory, the inverse problem and solution are defined via convolution and conjunction of probability density functions (PDFs) that represent stochastic information obtained from the model, observations, and prior knowledge in a joint multidimensional space. This theory provides an explicit analysis of the nonlinear model function, together with information about uncertainties in the model, observations, and prior knowledge through construction of the joint probability density, from which marginal solution functions can then be evaluated. The numerical analysis technique derived from the theory computes the component PDFs in discretized form via a combination of function mapping on a discrete grid in the model and observation phase space and Monte Carlo sampling from known parametric distributions. The efficacy of the numerical analysis technique is demonstrated through its application to two well-known simplified models of atmospheric physics: damped oscillations and Lorenz’s three-component model of dry cellular convection. The major findings of this study include the following: (i) Use of a nonmonotonic forward model in the inverse problem gives rise to the potential for a multimodal posterior PDF, the realization of which depends on the information content of the observations and on observation and model uncertainties. (ii) The cumulative effect of observations over time, space, or both could render the final posterior PDF unimodal, even with the nonmonotonic forward model. (iii) A greater number of independent observations are needed to constrain the solution in the case of a nonmonotonic nonlinear model than for a monotonic nonlinear or linear forward model for a given number of degrees of freedom in control parameter space. (iv) A nonlinear monotonic forward model gives rise to a skewed unimodal posterior PDF, implying a well-posed maximum likelihood inverse problem. (v) The presence of model error greatly increases the possibility of capturing multiple modes in the posterior PDF with the nonmonotonic nonlinear model. (vi) In the case of a nonlinear forward model, use of a Gaussian approximation for the prior update has a similar effect to an increase in model error, which indicates there is the potential to produce a biased mean central estimate even when observations and model are unbiased.

Knowledge of the Earth's interior relies on indirect information collected at or near the surface. Typically, data do not uniquely constrain the subsurface properties and are contaminated by noise, and therefore the solution of ill-posed inverse problems is required. Geophysicists have traditionally addressed such problems through deterministic approaches, seeking a single 'best-fitting' model defined by optimality criteria that reflect our understanding of the problem at hand. In recent decades, Bayesian approaches have become increasingly common, as they characterise the posterior probability distribution of the model conditioned on the observations, thereby quantifying uncertainty.In this context, Markov chain Monte Carlo (MCMC) methods have emerged as a crucial tool as they allow sampling from posterior distributions of arbitrary complexity. At the core of many MCMC algorithms lies the Metropolis-Hastings scheme. This combines a proposal distribution with a probabilistic acceptance criterion to construct a Markov chain that has the desired target distribution as its stationary distribution. The algorithm is versatile as it rests on mild technical conditions on the proposal, and is thus widely adopted across a broad range of geoscientific inference problems. Yet when the parameter space is large or the forward models are computationally expensive---both common scenarios in geophysical applications---it can become inefficient, resulting in poor chain mixing and slow convergence to the target (posterior) distribution. These challenges underscore the importance of effective proposal mechanisms.In this presentation, we introduce a novel approach to designing Metropolis-Hastings proposals based on adaptive transport maps. The framework is inspired by recent developments from the field of Applied Mathematics linking Bayesian inference with Optimal Transport theory. The idea is to find a monotone, nonlinear transformation to recast a (complex) target probability distribution into a (simpler) reference distribution that is more amenable to standard MCMC steps. Our key contribution is to parameterise these transformations using invertible neural networks, ensuring monotonicity while gaining the flexibility and expressiveness that neural architectures afford.The proposed method progresses iteratively. We begin with a standard sampling strategy (e.g., a random-walk Metropolis) to obtain initial draws from the target distribution. These samples inform the training of an invertible neural network that learns to map from the target to a simpler reference distribution, specifically a standard (multivariate) Gaussian. Subsequent proposals are then generated in the reference space, either as global independence moves or local perturbations, and are accepted or rejected following a suitably modified Metropolis-Hastings criterion. As more samples accumulate, the network's parameters are updated, improving overall sampling efficiency. Ultimately, the approach yields not only an ensemble of samples representative of the desired target distribution---just as in standard MCMC---but also a compact, learned representation of it in the network's weights.We illustrate the proposed paradigm both theoretically and through examples. Preliminary results indicate that transport-map-enhanced MCMC has the potential to significantly accelerate Bayesian sampling across a range of applications.

BackgroundTranslating a known metabolic network into a dynamic model requires reasonable guesses of all enzyme parameters. In Bayesian parameter estimation, model parameters are described by a posterior probability distribution, which scores the potential parameter sets, showing how well each of them agrees with the data and with the prior assumptions made.ResultsWe compute posterior distributions of kinetic parameters within a Bayesian framework, based on integration of kinetic, thermodynamic, metabolic, and proteomic data. The structure of the metabolic system (i.e., stoichiometries and enzyme regulation) needs to be known, and the reactions are modelled by convenience kinetics with thermodynamically independent parameters. The parameter posterior is computed in two separate steps: a first posterior summarises the available data on enzyme kinetic parameters; an improved second posterior is obtained by integrating metabolic fluxes, concentrations, and enzyme concentrations for one or more steady states. The data can be heterogenous, incomplete, and uncertain, and the posterior is approximated by a multivariate log-normal distribution. We apply the method to a model of the threonine synthesis pathway: the integration of metabolic data has little effect on the marginal posterior distributions of individual model parameters. Nevertheless, it leads to strong correlations between the parameters in the joint posterior distribution, which greatly improve the model predictions by the following Monte-Carlo simulations.ConclusionWe present a standardised method to translate metabolic networks into dynamic models. To determine the model parameters, evidence from various experimental data is combined and weighted using Bayesian parameter estimation. The resulting posterior parameter distribution describes a statistical ensemble of parameter sets; the parameter variances and correlations can account for missing knowledge, measurement uncertainties, or biological variability. The posterior distribution can be used to sample model instances and to obtain probabilistic statements about the model's dynamic behaviour.

This reference is for an abstract only. A full paper was not submitted for this conference. Abstract AVO (amplitude-varying-with-offset) is a seismic prestack inversion technique for estimating elastic parameters of the subsurface. The seismic AVO inversion problem can be formulated as a linearized or nonlinear inverse problem. It is a multidimensional and highly ill-posed inverse problem. It is affected by strong noise and measurement uncertainty. Therefore, in seismic AVO analysis, the goal is not only to find the model that best-fits the data, but also to characterize the uncertainty of the analysis. Uncertainty characterization of seismic analysis makes the geophysical interpretation more reliable. The Bayesian formalism generally constitutes a powerful approach for solving many seismic inverse problems. It allows combining available prior knowledge with the information contained in the seismic data. The solution of a Bayesian inverse problem is given by the joint posterior distribution of model parameters. Stochastic approximations of the true posterior are exemplified by the MCMC method and more recently by the particle filtering (PF) method. For most seismic inverse problems, stochastic approximation techniques are computationally unappealing. The maximum a posteriori (MAP) is the simplest Bayesian method for estimating model parameters. It provides the explanation that maximizes the posterior distribution. The maximum likelihood (ML) is a popular non-Bayesian method for estimating model parameters. It provides the explanation that maximizes the likelihood of the data. It assumes that the model parameters are deterministic and omits the pertinent prior knowledge on them. In practice, both MAP and ML point estimators have severe problems with over-fitting the data and model order estimation. In this paper, the recently-developed variational Bayesian learning is proposed as a new method for solving the linearized seismic AVO inversion problem. Variational Bayesian learning allows the true joint posterior distribution of model parameters to be approximated by a simpler approximating ensemble for which the required inferences are tractable. The main advantage in resorting to variational Bayesian learning is its robustness to the over-fitting problem, which is major in practice. This probabilistic machine learning method allows finding an optimum balance between the representation power and the model complexity of the data.

<p>Seismic travel time tomography is used widely to image the Earth's interior structure and to infer subsurface properties. Tomography is an inverse problem, and computationally expensive nonlinear inverse methods are often deployed in order to understand uncertainties in the tomographic results. Monte Carlo sampling methods estimate the posterior probability distribution which describes the solution to Bayesian tomographic problems, but they are computationally expensive and often intractable for high dimensional model spaces and large data sets due to the curse of dimensionality. We therefore introduce a new method of variational inference to solve Bayesian seismic tomography problems using optimization methods, while still providing fully nonlinear, probabilistic results. The new method, known as normalizing flows, warps a simple and known distribution (for example a Uniform or Gaussian distribution) into an optimal approximation to the posterior distribution through a chain of invertible transforms. These transforms are selected from a library of suitable functions, some of which invoke neural networks internally. We test the method using both synthetic and field data. The results show that normalizing flows can produce similar mean and uncertainty maps to those obtained from both Monte Carlo and another variational method (Stein varational gradient descent), at significantly decreased computational cost. In our tomographic tests, normalizing flows improves both accuracy and efficiency, producing maps of UK surface wave speeds and their uncertainties at the finest resolution and the lowest computational cost to-date, allowing results to be interrogated efficiently and quantitatively for subsurface structure.</p>

Existing studies on sparse Bayesian learning for structural damage detection usually assume that the posterior probability density functions follow standard distributions which facilitate to circumvent the intractable integration problem of the evidence by means of numerical sampling or analytical derivation. Moreover, the uncertainties of each mode are usually quantified as a common parameter to simplify the calculation. These assumptions may not be realistic in practice. This study proposes a sparse Bayesian method for structural damage detection suitable for standard and nonstandard probability distributions. The uncertainty corresponding to each mode is assumed as different. Variational Bayesian inference is developed and the posterior probability density functions of each unknown are individually derived. The parameters are found to follow the gamma distribution, whereas the distribution of the damage index cannot be directly obtained because of the nonlinear relationship in its posterior probability density function. The delayed rejection adaptive Metropolis algorithm is then adopted to generate numerical samples of the damage index. The coupled damage index and parameters in the variational Bayesian inference are successively calculated via an iterative process. A laboratory-tested frame is utilised to verify the effectiveness of the proposed method. The results indicate that the sparse damage can be accurately detected. The proposed method has the advantage of high accuracy and broad applicability.

Infrasonic signals propagate from an atmospheric source via media with stochastic and fast space-varying conditions. Hence, their travel time, the amplitude at sensor recordings and even manifestation in the so-called “shadow zones” are random. Therefore, the traditional least-squares technique for locating infrasonic sources is often not effective, and the problem for the best solution must be formulated in probabilistic terms. Recently, a series of papers has been published about Bayesian Infrasonic Source Localization (BISL) method based on the computation of the posterior probability density function (PPDF) of the source location, as a convolution of a priori probability distribution function (APDF) of the propagation model parameters with likelihood function (LF) of observations. The present study is devoted to the further development of BISL for higher accuracy and stability of the source location results and decreasing of computational load. We critically analyse previous algorithms and propose several new ones. First of all, we describe the general PPDF formulation and demonstrate that this relatively slow algorithm might be among the most accurate algorithms, provided the adequate APDF and LF are used. Then, we suggest using summation instead of integration in a general PPDF calculation for increased robustness, but this leads us to the 3D space-time optimization problem. Two different forms of APDF approximation are considered and applied for the PPDF calculation in our study. One of them is previously suggested, but not yet properly used is the so-called “celerity-range histograms” (CRHs). Another is the outcome from previous findings of linear mean travel time for the four first infrasonic phases in the overlapping consecutive distance ranges. This stochastic model is extended here to the regional distance of 1000 km, and the APDF introduced is the probabilistic form of the junction between this travel time model and range-dependent probability distributions of the phase arrival time picks. To illustrate the improvements in both computation time and location accuracy achieved, we compare location results for the new algorithms, previously published BISL-type algorithms and the least-squares location technique. This comparison is provided via a case study of different typical spatial data distributions and statistical experiment using the database of 36 ground-truth explosions from the Utah Test and Training Range (UTTR) recorded during the US summer season at USArray transportable seismic stations when they were near the site between 2006 and 2008.

Bayesian model updating provides a rigorous framework to account for uncertainty induced by lack of knowledge about engineering systems in their respective mathematical models through updates of the joint probability density function (PDF), the so-called posterior PDF, of the unknown model parameters. The Markov chain Monte Carlo (MCMC) methods are currently the most popular approaches for generating samples from the posterior PDF. However, these methods often found wanting when sampling from difficult distributions (e.g., high-dimensional PDFs, PDFs with flat manifolds, multimodal PDFs, and very peaked PDFs). This paper introduces a new multi-level sampling approach for Bayesian model updating, called Sequential Gauss-Newton algorithm, which is inspired by the Transitional Markov chain Monte Carlo (TMCMC) algorithm. The Sequential Gauss-Newton algorithm improves two aspects of TMCMC to make an efficient and effective MCMC algorithm for drawing samples from difficult posterior PDFs. First, the statistical efficiency of the algorithm is enhanced by use of the systematic resampling scheme. Second, a new MCMC algorithm, called Gauss-Newton MCMC algorithm, is proposed which is essentially an M-H algorithm with a Gaussian proposal PDF tailored to the posterior PDF using the gradient and Hessian information of the negative log posterior. The effectiveness of the proposed algorithm for solving the Bayesian model updating problem is illustrated using three examples with irregularly shaped posterior PDFs.

Generalized mixed linear, threshold, and logistic sire models and Markov chain, Monte Carlo simulation procedures were used to estimate genetic parameters for calving rate and calf survival in a multibreed beef cattle population. Data were obtained from a 5-generation rotational crossbreeding study involving Angus, Brahman, Charolais, and Hereford (1969 to 1995). Gelbvieh and Simmental bulls sired terminal-cross calves from a sample of generation 5 cows. A total of 1,458 cows sired by 158 bulls had a mean calving rate of 78% based on 4,808 calving records. Ninety-one percent of 5,015 calves sired by 260 bulls survived to weaning. Mean heritability estimates and standard deviations for daughter calving rate from posterior distributions were 0.063 +/- 0.024, 0.150 +/- 0.049, and 0.130 +/- 0.047 for linear, threshold, and logistic models, respectively. For calf survival, mean heritability estimates and standard deviations from posterior distributions were 0.049 +/- 0.022, 0.160 +/- 0.058, and 0.190 +/- 0.078 from linear, threshold, and logistic models, respectively. When transformed to an underlying normal scale, linear sire, mixed model, heritability estimates were similar to threshold and logistic sire mixed model estimates. Posterior density distributions of estimated heritabilities from all models were normal. Spearman rank correlations between sire EPD across statistical models were greater than 0.97 for daughter calving rate and for calf survival. Sire EPD had similar ranges across statistical models for daughter calving rate and for calf survival.

Computed tomography (CT) angiography (CTA) is a non-invasive imaging method used to detect arteries and examine various brain diseases. When CTA is performed for follow-up or postoperative evaluation, reproducibility of vessel delineation is required. A reproducible and stable contrast enhancement can be achieved by manipulating the factors affecting it. Previous studies have investigated several factors that alter the contrast enhancement of arteries. However, no reports establishing the effect of different operators on contrast enhancement exist. To assess the differences between inter-operator arterial contrast enhancement in cerebral CTA using Bayesian statistical modeling. Image data were obtained using a multistage sampling method from the cerebral CTA scans of patients who underwent the process between January 2015 and December 2018. Several Bayesian statistical models were developed, and the objective variable was the mean CT number of the bilateral internal carotid arteries after contrast enhancement. The explanatory variables were sex, age, fractional dose (FD), and the operator's information. The posterior distributions of the parameters were computed via Bayesian inference using the Markov chain Monte Carlo (MCMC) method, with the Hamiltonian Monte Carlo method employed as the algorithm. The posterior predictive distributions were computed using the posterior distributions of the parameters. Finally, the differences between inter-operator arterial contrast enhancement on the CT number in cerebral CTA were estimated. The posterior distributions showed that all parameters representing the difference between operators included zero at the 95% credible intervals (CIs). The maximum mean difference between inter-operator CT number in the posterior predictive distribution was only 12.59 Hounsfield units (HUs). The Bayesian statistical modeling results suggest that contrast enhancement of cerebral CTA examination between operator-to-operator differences in postcontrast CT number was small compared to those within-operator differences resulting from factors not considered in the model.

The Bayesian approach to inverse problems, in which the posterior probability distribution on an unknown field is sampled for the purposes of computing posterior expectations of quantities of interest, is starting to become computationally feasible for partial differential equation (PDE) inverse problems. Balancing the sources of error arising from finite-dimensional approximation of the unknown field, the PDE forward solution map and the sampling of the probability space under the posterior distribution are essential for the design of efficient computational Bayesian methods for PDE inverse problems. We study Bayesian inversion for a model elliptic PDE with an unknown diffusion coefficient. We provide complexity analyses of several Markov chain Monte Carlo (MCMC) methods for the efficient numerical evaluation of expectations under the Bayesian posterior distribution, given data δ. Particular attention is given to bounds on the overall work required to achieve a prescribed error level ε. Specifically, we first bound the computational complexity of ‘plain’ MCMC, based on combining MCMC sampling with linear complexity multi-level solvers for elliptic PDE. Our (new) work versus accuracy bounds show that the complexity of this approach can be quite prohibitive. Two strategies for reducing the computational complexity are then proposed and analyzed: first, a sparse, parametric and deterministic generalized polynomial chaos (gpc) ‘surrogate’ representation of the forward response map of the PDE over the entire parameter space, and, second, a novel multi-level Markov chain Monte Carlo strategy which utilizes sampling from a multi-level discretization of the posterior and the forward PDE. For both of these strategies, we derive asymptotic bounds on work versus accuracy, and hence asymptotic bounds on the computational complexity of the algorithms. In particular, we provide sufficient conditions on the regularity of the unknown coefficients of the PDE and on the approximation methods used, in order for the accelerations of MCMC resulting from these strategies to lead to complexity reductions over ‘plain’ MCMC algorithms for the Bayesian inversion of PDEs.