Search

We present a framework to solve geophysical inverse problems using the Hamiltonian Monte Carlo (HMC) method, with a focus on Bayesian tomography. Recent work in the geophysical community has shown the potential for gradient-based Monte Carlo sampling for a wide range of inverse problems across several fields. Many high-dimensional (non-linear) problems in geophysics have readily accessible gradient information which is unused in classical probabilistic inversions. Using HMC is a way to help improve traditional Monte Carlo sampling while increasing the scalability of inference problems, allowing access to uncertainty quantification for problems with many free parameters (>10'000). The result of HMC sampling is a collection of models representing the posterior probability density function, from which not only "best" models can be inferred, but also uncertainties and potentially different plausible scenarios, all compatible with the observed data. However, the amount of tuning parameters required by HMC, as well as the complexity of existing statistical modeling software, has limited the geophysical community in widely adopting a specific tool for performing efficient large-scale Bayesian inference. This work attempts to make a step towards filling that gap by providing an HMC sampler tailored for geophysical inverse problems (by e.g. supplying relevant priors and visualizations) combined with a set of different forward models, ranging from elastic and acoustic wave propagation to magnetic anomaly modeling, traveltimes, etc.. The framework is coded in the didactic but performant languages Julia and Python, with the possibility for the user to combine their own forward models, which are linked to the sampler routines by proper interfaces. In this way, we hope to illustrate the usefulness and potential of HMC in Bayesian inference. Tutorials featuring an array of physical experiments are written with the aim of both showcasing Bayesian inference and successful HMC usage. It additionally includes examples on how to speed up HMC e.g. with automated tuning techniques and GPU computations.

Methods that rely on Gaussian statistics require a choice of a mean and covariance to describe a Gaussian probability distribution. This is the case using for example kriging, sequential Gaussian simulation, least-squares collocation, and least-squares-based inversion, to name a few examples. Here, an approach is presented that provides a general description of a likelihood function that describes the probability that a set of, possibly noisy, data of both point and/or volume support is a realization from a Gaussian probability distribution with a specific set of Gaussian model parameters. Using this likelihood function, the problem of inferring the parameters of a Gaussian model is posed as a non-linear inverse problem using a general probabilistic formulation. The solution to the inverse problem is then the a posteriori probability distribution over the parameters describing a Gaussian model, from which a sample can be obtained using, e.g., the extended Metropolis algorithm. This approach allows detailed uncertainty and resolution analysis of the Gaussian model parameters. The method is tested on noisy data of both point and volume support, mimicking data from remote sensing and cross-hole tomography.

Some multiple-point-based sampling algorithms, such as the snesim algorithm, rely on sequential simulation. The conditional probability distributions that are used for the simulation are based on statistics of multiple-point data events obtained from a training image. During the simulation, data events with zero probability in the training image statistics may occur. This is handled by pruning the set of conditioning data until an event with non-zero probability is found. The resulting probability distribution sampled by such algorithms is a pruned mixture model. The pruning strategy leads to a probability distribution that lacks some of the information provided by the multiple-point statistics from the training image, which reduces the reproducibility of the training image patterns in the outcome realizations. When pruned mixture models are used as prior models for inverse problems, local re-simulations are performed to obtain perturbed realizations. Consequently, these local re-simulations lead to additional pruning in the set of conditioning data, which further deteriorates the pattern reproduction. To mitigate this problem, it is here suggested to combine the pruned mixture model with a frequency matching model. The multiple-point statistics of outcome realizations from this combined model has improved degree of match with the statistics from the training image. An efficient algorithm that samples this combined model is suggested. Finally, a tomographic cross-borehole inverse problem with prior information expressed by the combined (prior) model is used to demonstrate the effect of pattern reproducibility on the resolution of an inverse problem.

From a probabilistic point-of-view, the solution to an inverse problem can be seen as a combination of independent states of information quantified by probability density functions. Typically, these states of information are provided by a set of observed data and some a priori information on the solution. The combined states of information (i.e. the solution to the inverse problem) is a probability density function typically referred to as the a posteriori probability density function. We present a generic toolbox for Matlab and Gnu Octave called SIPPI that implements a number of methods for solving such probabilistically formulated inverse problems by sampling the a posteriori probability density function. In order to describe the a priori probability density function, we consider both simple Gaussian models and more complex (and realistic) a priori models based on higher order statistics. These a priori models can be used with both linear and non-linear inverse problems. For linear inverse Gaussian problems we make use of least-squares and kriging-based methods to describe the a posteriori probability density function directly. For general non-linear (i.e. non-Gaussian) inverse problems, we make use of the extended Metropolis algorithm to sample the a posteriori probability density function. Together with the extended Metropolis algorithm, we use sequential Gibbs sampling that allow computationally efficient sampling of complex a priori models. The toolbox can be applied to any inverse problem as long as a way of solving the forward problem is provided. Here we demonstrate the methods and algorithms available in SIPPI. An application of SIPPI, to a tomographic cross borehole inverse problems, is presented in a second part of this paper.

The frequency matching method defines a closed form expression for a complex prior that quantifies the higher order statistics of a proposed solution model to an inverse problem. While existing solution methods to inverse problems are capable of sampling the solution space while taking into account arbitrarily complex a priori information defined by sample algorithms, it is not possible to directly compute the maximum a posteriori model, as the prior probability of a solution model cannot be expressed. We demonstrate how the frequency matching method enables us to compute the maximum a posteriori solution model to an inverse problem by using a priori information based on multiple point statistics learned from training images. We demonstrate the applicability of the suggested method on a synthetic tomographic crosshole inverse problem.

Markov chain Monte Carlo methods such as the Gibbs sampler and the Metropolis algorithm can be used to sample solutions to non-linear inverse problems. In principle, these methods allow incorporation of prior information of arbitrary complexity. If an analytical closed form description of the prior is available, which is the case when the prior can be described by a multidimensional Gaussian distribution, such prior information can easily be considered. In reality, prior information is often more complex than can be described by the Gaussian model, and no closed form expression of the prior can be given. We propose an algorithm, called sequential Gibbs sampling, allowing the Metropolis algorithm to efficiently incorporate complex priors into the solution of an inverse problem, also for the case where no closed form description of the prior exists. First, we lay out the theoretical background for applying the sequential Gibbs sampler and illustrate how it works. Through two case studies, we demonstrate the application of the method to a linear image restoration problem and to a non-linear cross-borehole inversion problem. We demonstrate how prior information can reduce the complexity of an inverse problem and that a prior with little information leads to a hard inverse problem, practically unsolvable except when the number of model parameters is very small. Considering more complex and realistic prior information thus not only makes realizations from the posterior look more realistic but it can also reduce the computation time for the inversion dramatically. The method works for any statistical model for which sequential simulation can be used to generate realizations. This applies to most algorithms developed in the geostatistical community.

Linear inverse Gaussian problems are traditionally solved using least squares-based inversion. The center of the posterior Gaussian probability distribution is often chosen as the solution to such problems, while the solution is in fact the posterior Gaussian probability distribution itself. We present an algorithm, based on direct sequential simulation, which can be used to efficiently draw samples of the posterior probability distribution for linear inverse problems. There is no Gaussian restriction on the distribution in the model parameter space, as inherent in traditional least squares-based algorithms. As data for linear inverse problems can be seen as weighed linear averages over some volume, block kriging can be used to perform both estimation (i.e. finding the center of the posterior Gaussian pdf) and simulation (drawing samples of the posterior Gaussian pdf). We present the kriging system which we use to implement a flexible GSLIB-based algorithm for solving linear inverse problems. We show how we implement such a simulation program conditioned to linear average data. The program is called VISIM as an acronym for Volume average Integration SIMulation. An effort has been made to make the program efficient, even for larger scale problems, and the computational efficiency and accuracy of the code is investigated. Using a synthetic cross-borehole tomography case study, we show how the program can be used to generate realizations of the a posteriori distributions (i.e. solutions) from a linear tomography problem. Both Gaussian and non-Gaussian a priori model parameter distributions are considered.

We propose a smooth formulation of multiple-point statistics that enables us to solve inverse problems using gradient-based optimization techniques. We introduce a differentiable function that quantifies the mismatch between multiple-point statistics of a training image and of a given model. We show that, by minimizing this function, any continuous image can be gradually transformed into an image that honors the multiple-point statistics of the discrete training image. The solution to an inverse problem is then found by minimizing the sum of two mismatches: the mismatch with data and the mismatch with multiple-point statistics. As a result, in the framework of the Bayesian approach, such a solution belongs to a high posterior region. The methodology, while applicable to any inverse problem with a training-image-based prior, is especially beneficial for problems which require expensive forward simulations, as, for instance, history matching. We demonstrate the applicability of the method on a two-dimensional history matching problem. Starting from different initial models we obtain an ensemble of solutions fitting the data and prior information defined by the training image. At the end we propose a closed form expression for calculating the prior probabilities using the theory of multinomial distributions, that allows us to rank the history-matched models in accordance with their relative posterior probabilities., Mathematical Geosciences, 47 (4), ISSN:1874-8961, ISSN:0882-8121, ISSN:1874-8953, ISSN:1573-8868

We present a new method for solving the history matching problem by gradient-based optimization within a probabilistic framework. The focus is on minimizing the number of forward simulations and conserving geological realism of the solutions. Geological a priori information is taken into account by means of multipoint statistics borrowed from training images. Then production data and prior information are integrated into a single differentiable objective function, minimizer of which has a high posterior value. Solving the proposed optimization problem for an ensemble of different starting models, we obtain a set of solutions honouring both data and prior information.

Many inverse problems are most naturally formulated as implicit problems where data cannot be expressed in closed form as a function of the unknown model parameters. Notable examples are cases where the data satisfy a partial differential equation with model parameters as constants. When the problem can be recast as an ordinary explicit inverse problem the price to be paid is often that the functional relationship between data and model cannot be formulated in closed form. The consequence is typically that a computer time intensive numerical algorithm is needed to perform the forward calculation, thereby making an iterative solution of the problem unreasonably time consuming or impossible. In this paper we present a probabilistic procedure to solve rather general, implicit inverse problems through Monte Carlo sampling of feasible solutions. Our starting point is a probabilistic formulation of inverse problems, and our goal is to produce near-independent samples from the posterior distribution in model space. From these samples important information on the model and its resolution can be obtained. The proposed algorithm is applied to an implicit real-data problem involving analysis of the fluid motions at the Earth's core-mantle boundary from geomagnetic data.

In certain inverse problems it is useful to be able to compute solutions which are, in some sense, as simple as possible. For example,k one may wish to compute solutions which are piecewise constant and with as few discontinuities as possible. Such solutions are suited to describe models, e.g., geological layers, where the coarse structure is more important than the fine structure. A natural generalization of piecewise constant functions is piecewise polynomial solutions. In this paper we present a new algorithm which is capable of computing solutions that are piecewise polynomials, without having to specify a priori the positions of the break points between the polynomial pieces.

We present a general Monte Carlo full-waveform inversion strategy that integrates a priori information described by geostatistical algorithms with Bayesian inverse problem theory. The extended Metropolis algorithm can be used to sample the a posteriori probability density of highly nonlinear inverse problems, such as full-waveform inversion. Sequential Gibbs sampling is a method that allows efficient sampling of a priori probability densities described by geostatistical algorithms based on either two-point (e.g., Gaussian) or multiple-point statistics. We outline the theoretical framework for a full-waveform inversion strategy that integrates the extended Metropolis algorithm with sequential Gibbs sampling such that arbitrary complex geostatistically defined a priori information can be included. At the same time we show how temporally and/or spatiallycorrelated data uncertainties can be taken into account during the inversion. The suggested inversion strategy is tested on synthetic tomographic crosshole ground-penetrating radar full-waveform data using multiple-point-based a priori information. This is, to our knowledge, the first example of obtaining a posteriori realizations of a full-waveform inverse problem. Benefits of the proposed methodology compared with deterministic inversion approaches include: (1) The a posteriori model variability reflects the states of information provided by the data uncertainties and a priori information, which provides a means of obtaining resolution analysis. (2) Based on a posteriori realizations, complicated statistical questions can be answered, such as the probability of connectivity across a layer. (3) Complex a priori information can be included through geostatistical algorithms. These benefits, however, require more computing resources than traditional methods do. Moreover, an adequate knowledge of data uncertainties and a priori information is required to obtain meaningful uncertainty estimates. The latter may be a key challenge when considering field experiments, which will not be addressed here.

In order to move beyond simplified covariance based a priori models, which are typically used for inverse problems, more complex multiple-point-based a priori models have to be considered. By means of marginal probability distributions ‘learned’ from a training image, sequential simulation has proven to be an efficient way of obtaining multiple realizations that honor the same multiple-point statistics as the training image. The frequency matching method provides an alternative way of formulating multiple-point-based a priori models. In this strategy the pattern frequency distributions (i.e. marginals) of the training image and a subsurface model are matched in order to obtain a solution with the same multiple-point statistics as the training image. Sequential Gibbs sampling is a simulation strategy that provides an efficient way of applying sequential simulation based algorithms as a priori information in probabilistic inverse problems. Unfortunately, when this strategy is applied with the multiple-point-based simulation algorithm SNESIM the reproducibility of training image patterns is violated. In this study we suggest to combine sequential simulation with the frequency matching method in order to improve the pattern reproducibility while maintaining the efficiency of the sequential Gibbs sampling strategy. We compare realizations of three types of a priori models. Finally, the results are exemplified through crosshole travel time tomography.

Geophysical inference is characterized by non-linear relationships between model and data parameters, large model spaces describing spatial distributions of media properties, and intensive computations related to the numerical resolution of the forward problem. Although sampling approaches are convenient to solve such inverse problems, sometimes the involved computations are demanding. We consider here sampling techniques directed to improve the efficiency of sampling procedures in large real geophysical applications. We propose a sampling algorithm incorporating classical importance sampling within a two-step (or multistep) Markov chain sampler set-up. The first step of the algorithm is a Metropolis sampler ergodic to an importance density function and the second step is a Metropolis sampler correcting from the bias introduced by the importance density function; the combined algorithm samples the posterior probability density function asymptotically. The importance density is a combination of the prior density and an importance likelihood function, obtained as an approximation of the data likelihood function. Hence, the importance density is an approximation of the posterior density. Although rejection rates of the combined algorithm are larger than the rejection rates of the simple step Metropolis algorithm, the computational sampling efficiency is improved when the calculation of the importance density is much easier than the calculation of the actual posterior density. The algorithm is characterized by high acceptance rates in its second step because the first step serves as a barrier rejecting unlikely samples. Data likelihood approximations can be obtained in several ways : using simplified geophysical simulation, using the likelihood with partial (smaller) sets of the observed data, or using information obtained by preliminary analysis of the data. The technique has been successfully used for the inversion of phase arrival times from an offset vertical seismic profile (OVSP) data.

We have presented a new algorithm PP-TSVD that computes piecewise polynomial solutions to ill-posed problems, without a priori knowledge about the positions of the break points. In particular, we can compute piecewise constant functions that describe layered models. Such solutions are useful, e.g., in seismological problems, and the algorithm can also be used as a preprocessor for other methods where break points/discontinuities must be incorporated explicitly.

In ‘inverse problems’ data from indirect measurements are used to estimate unknown parameters of physical systems. Uncertain data, (possibly vague) prior information on model parameters, and a physical theory relating the model parameters to the observations are the fundamental elements of any inverse problem. Using concepts from probability theory, a consistent formulation of inverse problems can be made, and, while the most general solution of the inverse problem requires extensive use of Monte Carlo methods, special hypotheses (e.g., Gaussian uncertainties) allow, in some cases, an analytical solution to part of the problem (e.g., using the method of least squares).

Hansen and Mosegaard (2008), discuss the use of sequentialGaussian simulation as a tool to generate realizations of the aposteriori probability density function of linear inverseGaussian problems in the formulation of Tarantola and Valette(1982). They compare sequential Gaussian simulation to themethod of kriging through error simulation as proposed byJournel and Huijbregts (1978, p. 405) and applied by Deutschand Journel (1998). This method was later applied to linearaverage data (Go´mez-Herna´ndez et al., 2004; Gloaguen et al.,2004, 2005).In Section 5.4.3 and Fig. 13 Hansen and Mosegaard (2008)suggest that kriging through error simulation will not accuratelyhonornoisy dataobservations. This is however nottrueand due totheir implementation of kriging with error simulation. Marcotteet al. (2008) demonstrate how to account for data with uncorre-lated Gaussian noise using conditional simulation through errorsimulation. This can be generalized to using a noise modeldescribed by a correlated data covariance model as shown below.For noise free data a conditional simulation, m

Static Light Scattering.- Thermodynamic Geometry Determines Optimal Temperature Profile in Distillation Column.- Inversion of Seismic AVO Data.- Generalized Fourier Transforms and Inverse Scattering.- Multi-Step Samplers for Improving Efficiency in Probabilistic Geophysical Inference.- Full Waveform Inversion of Seismic Data for a Viscoelastic Medium.- Integrating Surface Reflection Data Inversion and Offest-VSP Modeling for Imaging of Converted S-Waves in Gulf Coast Hydrocarbon Reservoirs.- Inversion for the Moment Tensor and Source Location Using Love- and Rayleigh-Type Waves.- Resolution of Cole-Cole Parameters Based on Induced Polarization Data.- The Resolution Function in Linearized Born and Kirchoff Inversion.- SVD Analysis of a 3D Inverse Thermal Model.- Optimization Tools for Tikhonov Regularization of Nonlinear Equations Using the L-Curve and its Dual.- The PP-TSVD Algorithm for Image Restoration Problems.- Smooth versus Sharp Frechet Kernels in Time-Distance Helioseismology.- Covariance Modelling for Merging of Multi-Sensor Ocean Surface Data.- Inversion of 2D Geoelectrical Data by Iterated Multichannel Deconvolution.- Image Restoration as an Inverse Problem.- Resolution Studies of Fluid Flow Models Near the Core-Mantle Boundary Using Bayesian Inversion.- Inversion of GPS Occultation Data for Atomspheric Profiling.- A Limited Memory BFGS Method for an Inverse Problem in Atmospheric Imaging.

Summary In a probabilistic formulation of geophysical data inversion, prior information on subsurface parameters is combined with new information obtained from geophysical data, and an a posteriori probability distribution is defined in the model space. This distribution contains all accessible information about the considered inverse problem. Most often, realistic a priori information is highly complex, and the theory linking data with model parameters is nonlinear. All this leads to an a posteriori probability in the model space which may be difficult to describe (it may be multimodal, some moments may not be defined, etc.). When optimization techniques are used to obtain the maximum likelihood model, it is important that they are not trapped in secondary maxima. Except for very simple problems, this excludes all deterministic optimization methods, like gradient methods: only well designed Monte Carlo methods have a high probability of converging to the globally optimal model within practical timelimits. Obtaining the maximum likelihood model is, however, usually not sufficient. We often wish to have information on resolution power of the data, i.e., on the posterior uncertainties in the model space. Hence, we need to obtain information on the posterior probability as a whole. In the typical case we have a large number of parameters - for instance, when imaging the Earth’s interior using seismic reflection data - and limiting ourselves to inspection of the marginal probability densities of interest may be insufficient. But it is possible to pseudorandomly generate a large collection of models according to the posterior probability distribution, and display them in such a way that information on their relative likelihoods is conveyed to the spectator. One way is to display them sequentially, as the images of a movie. By using crude Monte Carlo methods it may take too long to sample the model space densely enough for visiting narrow maxima of the posterior probability distribution. However, this problem can be overcome by sampling the model space with a density proportional to the considered probability density function. The importance sampling method (Metropolis et al. [1953]) was designed to approximately perform such a sampling. Importance sampling can be generalized in such a way that it allows direct combination of a priori information over the model space with information obtained by theoretical analysis of noise contaminated data. Input to the main algorithm is Monte Carlo generated a priori models, whose frequencies of ocurrence are proportional to the prior probabilities. Another piece of input information is contained in a likelihood function describing the theoretical relationship between model and data, and the statistical properties of the noise. The output of the algorithm is a posteriori models whose density in the model space is proportional to the a posteriori probabilities. These output models can be used for direct, sequential inspection (a “movie”) and for extraction of information on model parameters (e.g. marginal and joint distributions, averages, correlations, approximate upper and lower bounds for single parameters or functions of several parameters).

The general inverse problem is characterized by at least one of the following two complications: (1) data can only be computed from the model by means of a numerical algorithm, and (2) the a priori model constraints can only be expressed via numerical algorithms. For linear problems and the so-called `weakly nonlinear problems', which can be locally approximated by a linear problem, analytical methods can provide estimates of the best fitting model and measures of resolution (nonuniqueness and uncertainty of solutions). This is, however, not possible for general problems. The only way to proceed is to use sampling methods that collect information on the posterior probability density in the model space. One such method is the inverse Monte Carlo strategy for resolution analysis suggested by Mosegaard and Tarantola. This method allows sampling of the posterior probability density even in cases where prior information is only available as an algorithm that samples the prior probability density. Once a collection of models sampled according to the posterior is available, it is possible to estimate, not only posterior model parameter covariances, but also resolution measures that are more useful in many applications. For example, posterior probabilities of the existence of interesting Earth structures like discontinuities and flow patterns can be estimated. These extended possibilities for resolution analysis may also provide new insight into problems that are usually treated by means of analytical methods.

Monte Carlo methods have become important in analysis of nonlinear inverse problems where no analytical expression for the forward relation between data and model parameters is available, and where linearization is unsuccessful. In such cases a direct mathematical treatment is impossible, but the forward relation materializes itself as an algorithm allowing data to be calculated for any given model. Monte Carlo methods can be divided into two categories: the sampling methods and the optimization methods. Monte Carlo sampling is useful when the space of feasible solutions is to be explored, and measures of resolution and uncertainty of solution are needed. The Metropolis algorithm and the Gibbs sampler are the most widely used Monte Carlo samplers for this purpose, but these methods can be refined and supplemented in various ways of which the neighbourhood algorithm is a notable example. Monte Carlo optimization methods are powerful tools when searching for globally optimal solutions amongst numerous local optima. Simulated annealing and genetic algorithms have shown their strength in this respect, but they suffer from the same fundamental problem as the Monte Carlo sampling methods: no provably optimal strategy for tuning these methods to a given problem has been found, only a number of approximate methods.