Geophysical Inverse Problems Research Articles

Bayesian Gaussian Mixture Linear Inversion for Geophysical Inverse Problems

A Bayesian linear inversion methodology based on Gaussian mixture models and its application to geophysical inverse problems are presented in this paper. The proposed inverse method is based on a Bayesian approach under the assumptions of a Gaussian mixture random field for the prior model and a Gaussian linear likelihood function. The model for the latent discrete variable is defined to be a stationary first-order Markov chain. In this approach, a recursive exact solution to an approximation of the posterior distribution of the inverse problem is proposed. A Markov chain Monte Carlo algorithm can be used to efficiently simulate realizations from the correct posterior model. Two inversion studies based on real well log data are presented, and the main results are the posterior distributions of the reservoir properties of interest, the corresponding predictions and prediction intervals, and a set of conditional realizations. The first application is a seismic inversion study for the prediction of lithological facies, P- and S-impedance, where an improvement of 30% in the root-mean-square error of the predictions compared to the traditional Gaussian inversion is obtained. The second application is a rock physics inversion study for the prediction of lithological facies, porosity, and clay volume, where predictions slightly improve compared to the Gaussian inversion approach.

Perturbational and nonperturbational inversion of Rayleigh-wave velocities

The inversion of Rayleigh-wave dispersion curves is a classic geophysical inverse problem. We have developed a set of MATLAB codes that performs forward modeling and inversion of Rayleigh-wave phase or group velocity measurements. We describe two different methods of inversion: a perturbational method based on finite elements and a nonperturbational method based on the recently developed Dix-type relation for Rayleigh waves. In practice, the nonperturbational method can be used to provide a good starting model that can be iteratively improved with the perturbational method. Although the perturbational method is well-known, we solve the forward problem using an eigenvalue/eigenvector solver instead of the conventional approach of root finding. Features of the codes include the ability to handle any mix of phase or group velocity measurements, combinations of modes of any order, the presence of a surface water layer, computation of partial derivatives due to changes in material properties and layer boundaries, and the implementation of an automatic grid of layers that is optimally suited for the depth sensitivity of Rayleigh waves.

S-wave velocity profiling from refraction microtremor Rayleigh wave dispersion curves via PSO inversion algorithm

The refraction microtremor method has been increasingly used as an appealing tool for investigating near surface S-wave structure. However, inversion, as a main stage in processing refraction microtremor data, is challenging for most local search methods due to its high nonlinearity. With the development of data optimization approaches, fast and easier techniques can be employed for processing geophysical data. Recently, particle swarm optimization algorithm has been used in many fields of studies. Use of particle swarm optimization in geophysical inverse problems is a relatively recent development which offers many advantages in dealing with the nonlinearity inherent in such applications. In this study, the reliability and efficiency of particle swarm optimization algorithm in the inversion of refraction microtremor data were investigated. A new framework was also proposed for the inversion of refraction microtremor Rayleigh wave dispersion curves. First, particle swarm optimization code in MATLAB was developed; then, in order to evaluate the efficiency and stability of proposed algorithm, two noise-free and two noise-corrupted synthetic datasets were inverted. Finally, particle swarm optimization inversion algorithm in refraction microtremor data was applied for geotechnical assessment in a case study in the area in city of Tabriz in northwest of Iran. The S-wave structure in the study area successfully delineated. Then, for evaluation, the estimated Vs profile was compared with downhole data available around of the considered area. It could be concluded that particle swarm optimization inversion algorithm is a suitable technique for inverting microtremor waves.

Appraisal of geodynamic inversion results: a data mining approach

Abstract Bayesian sampling based inversions require many thousands or even millions of forward models, depending on how nonlinear or non-unique the inverse problem is, and how many unknowns are involved. The result of such a probabilistic inversion is not a single ‘best-fit’ model, but rather a probability distribution that is represented by the entire model ensemble. Often, a geophysical inverse problem is non-unique, and the corresponding posterior distribution is multimodal, meaning that the distribution consists of clusters with similar models that represent the observations equally well. In these cases, we would like to visualize the characteristic model properties within each of these clusters of models. However, even for a moderate number of inversion parameters, a manual appraisal for a large number of models is not feasible. This poses the question whether it is possible to extract end-member models that represent each of the best-fit regions including their uncertainties. Here, I show how a machine learning tool can be used to characterize end-member models, including their uncertainties, from a complete model ensemble that represents a posterior probability distribution. The model ensemble used here results from a nonlinear geodynamic inverse problem, where rheological properties of the lithosphere are constrained from multiple geophysical observations. It is demonstrated that by taking vertical cross-sections through the effective viscosity structure of each of the models, the entire model ensemble can be classified into four end-member model categories that have a similar effective viscosity structure. These classification results are helpful to explore the non-uniqueness of the inverse problem and can be used to compute representative data fits for each of the end-member models. Conversely, these insights also reveal how new observational constraints could reduce the non-uniqueness. The method is not limited to geodynamic applications and a generalized MATLAB code is provided to perform the appraisal analysis.

Robust total-variation based geophysical inversion using split Bregman and proximity operators

In this paper, we take advantages of split Bregman and proximity operators to formulate geophysical inverse problems in a convex optimization setting. In this type of formulation, we can efficiently consider a general convex data misfit term in order to incorporate more realistic error probability assumptions than the ordinary Gaussian distribution. Furthermore, we can simply impose convex non-quadratic and non-smooth regularization terms as a prior information.Although the proposed formulation can be extended for incorporating any types of convex regularization and data fidelity terms, here we consider misfit term corresponding to the Huber norm and anisotropic total variation regularization as a prior to an unknown model structure. This type of regularization enables us to construct piecewise-constant solutions.Two-dimensional linear and non-linear seismic travel-time tomography are studied to evaluate the efficiency of the proposed formulation for handling a robust measure of the misfit term in recovering blocky solutions. We consider the first arrival travel-times which are contaminated by outliers and also a mixed combination of Gaussian and outliers. Finally, Seinsfeld cross-hole tomography data set is used to investigate the performance of the robust approach on real data sets.

Proper generalized decomposition solution of the parameterized Helmholtz problem: application to inverse geophysical problems

SummaryThe identification of the geological structure from seismic data is formulated as an inverse problem. The properties and the shape of the rock formations in the subsoil are described by material and geometric parameters, which are taken as input data for a predictive model. Here, the model is based on the Helmholtz equation, describing the acoustic response of the system for a given wave length. Thus, the inverse problem consists in identifying the values of these parameters such that the output of the model agrees the best with observations. This optimization algorithm requires multiple queries to the model with different values of the parameters. Reduced order models are especially well suited to significantly reduce the computational overhead of the multiple evaluations of the model.In particular, the proper generalized decomposition produces a solution explicitly stating the parametric dependence, where the parameters play the same role as the physical coordinates. A proper generalized decomposition solver is devised to inexpensively explore the parametric space along the iterative process. This exploration of the parametric space is in fact seen as a post‐process of the generalized solution. The approach adopted demonstrates its viability when tested in two illustrative examples. Copyright © 2016 John Wiley & Sons, Ltd.

Intelligent meshing technique for 2D resistivity inverse problems

In geophysical inverse problems, an a priori structured mesh is often used for inversion and mesh refinement is applied if needed by the user after observation of inversion results. We have developed a new intelligent self-adaptive unstructured finite-element meshing technique for electrical resistivity tomography inverse problems. This new approach uses Harris corner-and-edge detectors that are based on the local autocorrelation function of 2D distribution of pixels. This meshing technique optimizes the size of the inverse problem by refining areas where variations in the physical property structure are sensed to be important. The meshing technique also generates a more appropriate and optimum mesh for the inverse problem that is dependent on the problem itself. Tests on modeled data have demonstrated that the proposed intelligent meshing technique can reduce data misfit, produce a better reconstruction of the true physical properties, and minimize the size of the inverse problem. The synthetic model consists of a conductive dike in a resistive medium. By applying the proposed intelligent meshing technique, the inverse model of the dike is very similar to the inverse model produced using fine meshes, and it is also better reconstructed than the inverse model produced using conventional meshes. We have also applied the intelligent meshing technique to survey data collected for groundwater-saltwater mapping and characterizing the subsurface conductive structure with topography included. Our results indicate that the new meshing technique can produce solutions that are comparable with standard meshing and fine meshing techniques, while optimizing the size of the inverse problem.

PDE-based geophysical modelling using finite elements: examples from 3D resistivity and 2D magnetotellurics

We present a general finite-element solver, escript, tailored to solve geophysical forward and inverse modeling problems in terms of partial differential equations (PDEs) with suitable boundary conditions. Escript's abstract interface allows geoscientists to focus on solving the actual problem without being experts in numerical modeling. General-purpose finite element solvers have found wide use especially in engineering fields and find increasing application in the geophysical disciplines as these offer a single interface to tackle different geophysical problems. These solvers are useful for data interpretation and for research, but can also be a useful tool in educational settings. This paper serves as an introduction into PDE-based modeling with escript where we demonstrate in detail how escript is used to solve two different forward modeling problems from applied geophysics (3D DC resistivity and 2D magnetotellurics). Based on these two different cases, other geophysical modeling work can easily be realized. The escript package is implemented as a Python library and allows the solution of coupled, linear or non-linear, time-dependent PDEs. Parallel execution for both shared and distributed memory architectures is supported and can be used without modifications to the scripts.

Solving probabilistic inverse problems rapidly with prior samples

Owing to the increasing availability of computational resources, in recent years the probabilistic solution of non-linear, geophysical inverse problems by means of sampling methods has become increasingly feasible. Nevertheless, we still face situations in which a Monte Carlo approach is not practical. This is particularly true in cases where the evaluation of the forward problem is computationally intensive or where inversions have to be carried out repeatedly or in a timely manner, as in natural hazards monitoring tasks such as earthquake early warning. Here,we present an alternative to Monte Carlo sampling, inwhich inferences are entirely based on a set of prior samples-that is, samples that have been obtained independent of a particular observed datum. This has the advantage that the computationally expensive sampling stage becomes separated from the inversion stage, and the set of prior samples-once obtained-can be reused for repeated evaluations of the inverse mapping without additional computational effort. This property is useful if the problem is such that repeated inversions of independent data have to be carried out.We formulate the inverse problem in a Bayesian framework and present a practical way to make posterior inferences based on a set of prior samples. We compare the prior sampling based approach to a Markov Chain Monte Carlo approach that samples from the posterior probability distribution. We show results for both a toy example, and a realistic seismological source parameter estimation problem. We find that the posterior uncertainty estimates obtained based on prior sampling can be considered conservative estimates of the uncertainties obtained by directly sampling from the posterior distribution. Published by Oxford University Press on behalf of The Royal Astronomical Society.

Minimax approach to inverse problems of geophysics

A new approach is suggested for solving the inverse problems that arise in the different fields of applied geophysics (gravity, magnetic, and electrical prospecting, geothermy) and require assessing the spatial region occupied by the anomaly-generating masses in the presence of different types of a priori information. The interpretation which provides the maximum guaranteed proximity of the model field sources to the real perturbing object is treated as the best interpretation. In some fields of science (game theory, economics, operations research), the decision-making principle that lies in minimizing the probable losses which cannot be prevented if the situation develops by the worst-case scenario is referred to as minimax. The minimax criterion of choice is interesting as, instead of being confined to the indirect (and sometimes doubtful) signs of the “optimal” solution, it relies on the actual properties of the information in the results of a particular interpretation. In the hierarchy of the approaches to the solution of the inverse problems of geophysics ordered by the volume and quality of the retrieved information about the sources of the field, the minimax approach should take special place.

Geodetic inversion for spatial distribution of slip under smoothness, discontinuity, and sparsity constraints

In geodetic data inversion, insufficient observational data and smoothness constraints for model parameters make it difficult to clearly resolve small-scale heterogeneous structures with discontinuous boundaries. We therefore developed a novel regularization scheme for the inversion problem that uses discontinuity, sparsity, and smoothness constraints. In order to assess its usefulness and applicability, the proposed method was applied to synthetic displacements calculated by a ring-shaped and sharply varying afterslip distribution on a plate interface. The afterslip was obtained from reasonable numerical simulation of earthquake generation cycle with a rate- and state- dependent friction law and realistic three-dimensional plate geometry. The obtained afterslip distribution was heterogeneous, and the discontinuous boundary was sharper than that obtained by using smoothness constraint only. The same inversion test was conducted with a smoothly varying circular slip distribution with large slips inside the ring-shaped distribution. The method accurately reproduces the smooth distribution of the slip area as well as the ring-shaped distribution. Therefore, the method could be applied to any slip distribution, with both discontinuous and continuous boundaries. Adopting this method for measured data will make it possible to obtain detailed heterogeneous distributions of physical structures on fault planes. The proposed method is therefore applicable to various geophysical inversion problems that exhibit discontinuous heterogeneity.

A parallel evolution strategy for an earth imaging problem in geophysics

In this paper we propose a new way to compute a rough approximation solution, to be later used as a warm starting point in a more refined optimization process, for a challenging global optimization problem related to earth imaging in geophysics. The warm start consists of a velocity model that approximately solves a full-waveform inverse problem at low frequency. Our motivation arises from the availability of massively parallel computing platforms and the natural parallelization of evolution strategies as global optimization methods for continuous variables. Our first contribution consists of developing a new and efficient parametrization of the velocity models to significantly reduce the dimension of the original optimization space. Our second contribution is to adapt a class of evolution strategies to the specificity of the physical problem at hands where the objective function evaluation is known to be the most expensive computational part. A third contribution is the development of a parallel evolution strategy solver, taking advantage of a recently proposed modification of these class of evolutionary methods that ensures convergence and promotes better performance under moderate budgets. The numerical results presented demonstrate the effectiveness of the algorithm on a realistic 3D full-waveform inverse problem in geophysics. The developed numerical approach allows us to successfully solve an acoustic full-waveform inversion problem at low frequencies on a reasonable number of cores of a distributed memory computer.

2D joint inversion of geophysical data using petrophysical clustering and facies deformation

Geologic expertise and petrophysical relationships can be brought together to provide prior information while inverting multiple geophysical data sets. The merging of such information can result in more realistic solution in the distribution of the model parameters. We have evaluated the geophysical inverse problem in terms of Gaussian random fields with mean functions controlled by petrophysical relationships and covariance functions controlled by a prior geologic cross section, including the definition of spatial boundaries for the geologic facies. The petrophysical relationship problem is formulated as a regression problem upon each facies. The inversion of the geophysical data is performed in a Bayesian framework. We have developed the usefulness of this strategy using a first synthetic case for which we have performed a joint inversion of gravity and galvanometric resistivity data with the stations located at the ground surface. The joint inversion was used to recover the density and resistivity distributions of the subsurface. In a second step, we considered the possibility that the facies boundaries were deformable and their shapes were inverted as well. We used the level-set approach to perform such deformation preserving a priori topological properties of the facies throughout the inversion. With the help of a priori facies petrophysical relationships and the topological characteristics of each facies, we made a posteriori inference about multiple geophysical tomograms based on their corresponding geophysical data misfits. We have applied this method to a second synthetic case showing that we can recover the heterogeneities inside the facies, the mean values for the petrophysical properties, and, to some extent, the facies boundaries using the 2D joint inversion of gravity and galvanometric resistivity data.

Iterative estimation of reflectivity and image texture: Least-squares migration with an empirical Bayes approach

In many geophysical inverse problems, smoothness assumptions on the underlying geology are used to mitigate the effects of nonuniqueness, poor data coverage, and noise in the data and to improve the quality of the inferred model parameters. Within a Bayesian inference framework, a priori assumptions about the probabilistic structure of the model parameters can impose such a smoothness constraint, analogous to regularization in a deterministic inverse problem. We have considered an empirical Bayes generalization of the Kirchhoff-based least-squares migration (LSM) problem. We have developed a novel methodology for estimation of the reflectivity model and regularization parameters, using a Bayesian statistical framework that treats both of these as random variables to be inferred from the data. Hence, rather than fixing the regularization parameters prior to inverting for the image, we allow the data to dictate where to regularize. Estimating these regularization parameters gives us information about the degree of conditional correlation (or lack thereof) between neighboring image parameters, and, subsequently, incorporating this information in the final model produces more clearly visible discontinuities in the estimated image. The inference framework is verified on 2D synthetic data sets, in which the empirical Bayes imaging results significantly outperform standard LSM images. We note that although we evaluated this method within the context of seismic imaging, it is in fact a general methodology that can be applied to any linear inverse problem in which there are spatially varying correlations in the model parameter space.

Regularized sparse-grid geometric sampling for uncertainty analysis in non-linear inverse problems

This paper introduces an efficiency improvement to the sparse-grid geometric sampling methodology for assessing uncertainty in non-linear geophysical inverse problems. Traditional sparse-grid geometric sampling works by sampling in a reduced-dimension parameter space bounded by a feasible polytope, e.g., a generalization of a polygon to dimension above two. The feasible polytope is approximated by a hypercube. When the polytope is very irregular, the hypercube can be a poor approximation leading to computational inefficiency in sampling. We show how the polytope can be regularized using a rotation and scaling based on principal component analysis. This simple regularization helps to increase the efficiency of the sampling and by extension the computational complexity of the uncertainty solution. We demonstrate this on two synthetic 1D examples related to controlled-source electromagnetic and amplitude versus offset inversion. The results show an improvement of about 50% in the performance of the proposed methodology when compared with the traditional one. However, as the amplitude versus offset example shows, the differences in the efficiency of the proposed methodology are very likely to be dependent on the shape and complexity of the original polytope. However, it is necessary to pursue further investigations on the regularization of the original polytope in order to fully understand when a simple regularization step based on rotation and scaling is enough.

Systematic evaluation of sequential geostatistical resampling within MCMC for posterior sampling of near-surface geophysical inverse problems

SUMMARY We critically examine the performance of sequential geostatistical resampling (SGR) as a model proposal mechanism for Bayesian Markov-chain-Monte-Carlo (MCMC) solutions to near-surfacegeophysicalinverseproblems.Focusingonaseriesofsimpleyetrealisticsynthetic crosshole georadar tomographic examples characterized by different numbers of data, levels of data error and degrees of model parameter spatial correlation, we investigate the efficiency of three different resampling strategies with regard to their ability to generate statistically independentrealizationsfromtheBayesianposteriordistribution.Quiteimportantly,ourresults show that, no matter what resampling strategy is employed, many of the examined test cases require an unreasonably high number of forward model runs to produce independent posterior samples,meaningthattheSGRapproachascurrentlyimplementedwillnotbecomputationally feasible for a wide range of problems. Although use of a novel gradual-deformation-based proposal method can help to alleviate these issues, it does not offer a full solution. Further, we find that the nature of the SGR is found to strongly influence MCMC performance; however no clear rule exists as to what set of inversion parameters and/or overall proposal acceptance rate will allow for the most efficient implementation. We conclude that although the SGR methodology is highly attractive as it allows for the consideration of complex geostatistical priors as well as conditioning to hard and soft data, further developments are necessary in the context of novel or hybrid MCMC approaches for it to be considered generally suitable for near-surface geophysical inversions.

Improved analysis‐error covariance matrix for high‐dimensional variational inversions: application to source estimation using a 3D atmospheric transport model

Variational methods are widely used to solve geophysical inverse problems. Although gradient‐based minimization algorithms are available for high‐dimensional problems (dimension >106), they do not provide an estimate of the errors in the optimal solution. In this study, we assess the performance of several numerical methods to approximate the analysis‐error covariance matrix, assuming reasonably linear models. The evaluation is performed for a CO2 flux estimation problem using synthetic remote‐sensing observations of CO2 columns. A low‐dimensional experiment is considered in order to compare the analysis error approximations to a full‐rank finite‐difference inverse Hessian estimate, followed by a realistic high‐dimensional application. Two stochastic approaches, a Monte‐Carlo simulation and a method based on random gradients of the cost function, produced analysis error variances with a relative error <10%. The long‐distance error correlations due to sampling noise are significantly less pronounced for the gradient‐based randomization, which is also particularly attractive when implemented in parallel. Deterministic evaluations of the inverse Hessian using the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm are also tested. While existing BFGS preconditioning techniques yield poor approximations of the error variances (relative error >120%), a new preconditioner that efficiently accumulates information on the diagonal of the inverse Hessian dramatically improves the results (relative error <50%). Furthermore, performing several cycles of the BFGS algorithm using the same gradient and vector pairs enhances its performance (relative error <30%) and is necessary to obtain convergence. Leveraging those findings, we proposed a BFGS hybrid approach which combines the new preconditioner with several BFGS cycles using information from a few (3–5) Monte‐Carlo simulations. Its performance is comparable to the stochastic approximations for the low‐dimensional case, while good scalability is obtained for the high‐dimensional experiment. Potential applications of these new BFGS methods range from characterizing the information content of high‐dimensional inverse problems to improving the convergence rate of current minimization algorithms.

Deterministic models of interpretation for optimizing the locations and depths of the boreholes for verifying the anomalies in gravity

We present the algorithm and results of the first experience in synthesizing the updated deterministic models of interpretation and simplest statistical techniques for determining the coordinates and depths of the verification boreholes based on the results of gravity prospecting. The suggested algorithm implements the ideas of guaranteed approach and construction of generalized assembly algorithms for finding the separate variants of interpretation if the a priori information is available for the sources of the field. In this paper, we present the first results of joint application of the deterministic techniques for constructing the representative subset of the admissible solutions of the inverse problem of gravity prospecting and describe the simplest statistical methods for their processing for elaborating the substantiated recommendations on determining the surface coordinates and depths of the verification boreholes. The suggested algorithms rely on the concept of a guaranteed approach to the inverse problems of geophysics and on the methodological idea of reconciling the different admissible interpretations of gravity data intended aimed at increasing the reliability of the conclusions. For finding the individual admissible variants of solution of the inverse problem, we use the finite-element descriptions of density heterogeneities and modern modifications of the assembly method suggested by V.V. Strakhov, which are flexible for the formalizations and allowance for heterogeneous a priori information.

Comparative analysis of the solution of linear continuous inverse problems using different basis expansions

Geophysical inverse problems try to infer the value of a physical property of the earth from data measured at the boundary of the domain. Model parameterization is a key concept to make the inverse problem less ill conditioned. In this contribution we compare the performance of four different basis functions, Fourier, Pixel, Haar and Daubechies Wavelets, for a 1D linear continuous inverse problem. Adopting the right parameterization also reduces the number of dimensions in which the inverse problem is going to be solved, allowing easy posterior analysis. To compare the different basis expansions we have studied a simple 1D linear inverse problem in gravimetric inversion for a density anomaly with Gaussian shape. For this simple toy-problem we show that the Fourier basis gives a better reconstruction using Filon's quadrature to avoid numerical instabilities caused by highly oscillating functions, both, in the noise-free and noisy cases. Besides, the Fourier base is the one that provides the lowest reconstruction error for the number of basis terms and the highest condition number. The Haar and pixel (piecewise-continuous functions) bases provide similar results, although the pixel base needs more terms to achieve lower reconstruction errors. Also, the pixel basis is the one that has the lowest condition number that is obviously related to the way the energy in the system matrix is distributed. Besides, the Daubechies Db2 basis expansion is the one that has the highest system rank, it is more difficult to apply to project the integral kernel, and provides the worst results, compared to the other basis expansions. Finally we propose how to generalize this methodology to linear inverse problems in several dimensions and to non-linear problems. A separate paper will be devoted to this subject.

Efficient trans-dimensional Bayesian inversion for geoacoustic profile estimation

This paper considers the efficiency of trans-dimensional (trans-D) Bayesian inversion based on reversible-jump Markov-chain Monte Carlo (rjMCMC) sampling, with application to geophysical inverse problems for a depth-dependent earth or seabed model of an unknown number of layers (seabed acoustic reflectivity inversion is the specific example). Trans-D inversion is applied to sample the posterior probability density over geoacoustic/geophysical parameters for a variable number of layers, providing profile estimates with uncertainties that include the uncertainty in the model parameterization. However, the approach is computationally intensive. The efficiency of rjMCMC sampling is largely determined by the proposal schemes which are applied to generate perturbed values for existing parameters and new values for parameters assigned to layers added to the model. Several proposal schemes are considered here, some of which appear new for trans-D geophysical inversion. Perturbations of existing parameters are considered in a principal-component space based on an eigen-decomposition of the unit-lag parameter covariance matrix (computed from successive models along the Markov chain, a diminishing adaptation). The relative efficiency of proposing new parameters from the prior versus a Gaussian distribution focused near existing values is examined. Parallel tempering, which employs a sequence of interacting Markov chains in which the likelihood function is successively relaxed, is also considered as a means to increase the acceptance rate of new layers. The relative efficiency of various proposal schemes is compared through repeated inversions with a pragmatic convergence criterion.