Klaus Mosegaard

Monte Carlo sampling of solutions to inverse problems

Probabilistic formulation of inverse problems leads to the definition of a probability distribution in the model space. This probability distribution combines a priori information with new information obtained by measuring some observable parameters (data). As, in the general case, the theory linking data with model parameters is nonlinear, the a posteriori probability in the model space may not be easy to describe (it may be multimodal, some moments may not be defined, etc.). When analyzing an inverse problem, obtaining a maximum likelihood model is usually not sucient, as we normally also wish to have infor

Monte Carlo analysis of inverse problems

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.

Probabilistic Approach to Inverse Problems

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).

A New Seismic Velocity Model for the Moon from a Monte Carlo Inversion of the Apollo Lunar Seismic Data

A reanalysis of the Apollo lunar seismic data and the subsequent application of an inverse Monte Carlo method to P and S-wave arrival times has resulted in a more detailed lunar velocity structure than previously obtainable. The velocity is seen to increase from the surface down to the base of the crust at 45±5 km depth. The results furthermore indicate a constant velocity upper mantle extending to 560±15 km km depth, separated from a more complex high velocity middle mantle by an increase in velocity of 1.0 km/s. In addition, the moonquake locations have been improved. The shallow moonquakes are found to be located in the depth range 50–220 km. The majority of deep moonquakes are concentrated in the depth range 850–1000 km with an apparently rather sharp lower boundary.

Linear inverse Gaussian theory and geostatistics

Inverse problems in geophysics require the introduction of complex a priori information and are solved using computationally expensive Monte Carlo techniques (where large portions of the model space are explored). The geostatistical method allows for fast integration of complex a priori information in the form of covariance functions and training images. We combine geostatistical methods and inverse problem theory to generate realizations of the posterior probability density function of any Gaussian linear inverse problem, honoring a priori information in the form of a covariance function describing the spatial connectivity of the model space parameters. This is achieved using sequential Gaussian simulation, a well-known, noniterative geostatisticalmethod for generating samples of a Gaussian random field with a given covariance function. This work is a contribution to both linear inverse problem theory and geostatistics. Our main result is an efficient method to generate realizations, actual solutions rat...

Depth to Moho in Greenland: receiver-function analysis suggests two Proterozoic blocks in Greenland

Abstract The GLATIS project (Greenland Lithosphere Analysed Teleseismically on the Ice Sheet) with collaborators has operated a total of 16 temporary broadband seismographs for periods from 3 months to 2 years distributed over much of Greenland from late 1999 to the present. The very first results are presented in this paper, where receiver-function analysis has been used to map the depth to Moho in a large region where crustal thicknesses were previously completely unknown. The results suggest that the Proterozoic part of central Greenland consists of two distinct blocks with different depths to Moho. North of the Archean core in southern Greenland is a zone of very thick Proterozoic crust with an average depth to Moho close to 48 km. Further to the north the Proterozoic crust thins to 37–42 km. We suggest that the boundary between thick and thin crust forms the boundary between the geologically defined Nagssugtoqidian and Rinkian mobile belts, which thus can be viewed as two blocks, based on the large difference in depth to Moho (over 6 km). Depth to Moho on the Archean crust is around 40 km. Four of the stations are placed in the interior of Greenland on the ice sheet, where we find the data quality excellent, but receiver-function analyses are complicated by strong converted phases generated at the base of the ice sheet, which in some places is more than 3 km thick.

Inverse problems with non-trivial priors: efficient solution through sequential Gibbs sampling

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.

Monte Carlo estimation and resolution analysis of seismic background velocities

The complete solution to an inverse problem, including information on accuracy and resolution, is given by the a posteriori probability density in the model space. By running a modified simulated annealing, samples from the model space can be drawn in such a way that their frequencies of occurrence approximate their a posteriori likelihoods. Using this method, maximum likelihood estimation and uncertainty analysis of seismic background velocity models are performed on multioffset seismic data. The misfit between observed and synthetic waveforms within the time windows along computed multioffset travel times, is used as an objective function for the simulated annealing approach. The real medium is modeled as a series of layers separated by curved interfaces. Lateral velocity variations within the layers are determined by interpolation from specified values at a number of sampling points. The input data consists of multioffset seismic data. Additionally, zero-offset times are used to migrate the reflectors in time to the depth domain. The multioffset times are calculated by an efficient ray-tracing algorithm which allows inversion of a large number of seismograms. The a posteriori probability density for this problem is highly multidimensional and highly multimodal. Therefore, the information contained in this distribution cannot be adequately represented by standard deviations and covariances. However, by sequentially displaying a large number of images, computed from the a posteriori background velocity samples and the data, it is possible to convey to the spectator a better understanding of what information we really have on the subsurface.

VISIM: Sequential simulation for linear inverse problems

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.

New information on the deep lunar interior from an inversion of lunar free oscillation periods

We have obtained crucial information on the deep lunar interior through a Monte Carlo inversion of a number of fundamental lunar spheroidal free oscillations. The results indicate a homogeneous upper mantle with a velocity and density of 3.751.1 km/s and 3.30.5 g/cm 3 . A transition between the upper and middle mantle in the range 520-580 km depth also seems to be implied. Mid- dle mantle velocities and densities are 5.41.2 km/s and 3.50.2 g/cm 3 . Moreover, a velocity decrease is seen in the lower mantle (1150-1450 km depth) consistent with earlier inferences of partial melt in this region. The density seems to increase gradually from 1250 km depth to the center of the moon, where it reaches values indicative of an FeS or silicate composition.