Markov Chain Monte Carlo Non-linear Geophysical Inversion with an Improved Proposal Distribution: Application to Geo-electrical Data

نوع مقاله : مقاله پژوهشی


1 موسسه ژئوفیزیک دانشگاه تهران

2 استادیار موسسه ژئوفیزیک دانشگاه تهران


Geophysical inverse problems seek to provide quantitative information about geophysical characteristics of the Earth’s subsurface for indirectly related data and measurements. It is generally formulated as an ill-posed non-linear optimization problem commonly solved through deterministic gradient-based approaches. Using these methods, despite fast convergence properties, may lead to local minima as well as impend accurate uncertainty analysis. On the contrary, formulating a geophysical inverse problem in a probabilistic framework and solving it by constructing the multi-dimensional posterior probability density (PPD) allow for complete sampling of the parameter space and uncertainty quantification. The PPD is numerically characterized using Markov Chain Monte Carlo (MCMC) approaches. However, the convergence of the MCMC algorithm (i.e. sampling efficiency) toward the target stationary distribution highly depends upon the choice of the proposal distribution. In this paper, we develop an efficient proposal distribution based on perturbing the model parameters through an eigenvalue decomposition of the model covariance matrix in a principal component space. The covariance matrix is retrieved from an initial burn-in sampling, which is itself initiated using a linearized covariance estimate. The proposed strategy is first illustrated for inversion of hydrogeological parameters and then applied to synthetic and real geo-electrical data sets. The numerical experiments demonstrate that the presented proposal distribution takes advantage of the benefits from an accelerated convergence and mixing rate compared to the conventional Gaussian proposal distribution.



عنوان مقاله [English]

Markov Chain Monte Carlo Non-linear Geophysical Inversion with an Improved Proposal Distribution: Application to Geo-electrical Data

نویسندگان [English]

  • Zahra Tafaghod Khabaz 1
  • Reza Ghanati 2
1 University of Tehran Institute of Geophysics
2 Assistant professor of Institute of Geophysics, University of Tehran
چکیده [English]

Observations of the altered data consist of information regarding the Earth’s interior physical properties which are not directly available to the surface or borehole measurements but are inferred by solving an inverse problem or an inductive reasoning process in a sense of logic. The data-earth interaction is described by a model which states the physical theory, an appropriate subsurface parameterization, and a statistical representation of the data error which may also be characterized by parameters in the model. Geophysical inverse problems are generally formulated as an ill-posed non-linear optimization problem commonly solved through Newton-based methods. The significant property of the gradient-based approaches (e.g., steepest descent, conjugate gradient, and Landweber iterative scheme) is their fast convergence toward the final solution, but the local linearization of the inverse solution hinders reliable uncertainty appraisal (i.e., underestimated or overestimated uncertainty). In addition, these algorithms may get trapped in local minima if the initial model is far from the convergence region of a global minimum of the cost function so that slight variations in starting model can lead to a notably different subsurface model. A remedy to the exiting limitations is to employ derivative-free global direct search techniques, e.g., sampling of the posterior probability distribution using Markov Chain Monte Carlo (MCMC) algorithms in the Bayesian framework. The MCMC sampling is essentially a guided-random walk through the probable parts of the posterior model space (Tarantola, 1987; Sambridge and Mosengaard, 2002). With MCMC algorithms, the influence of the initial model diminishes as the model space sampling progresses. The implementation of the MCMC algorithms requires the sampling of a large number of models, and consequently, the forward calculation of these models can be computationally expensive. Despite the advantages of the Bayesian inversion over the Newton-based optimization approaches, computational efficiency remains a critical factor for high-dimensional model spaces. However, with increasing computational power, the application of the Bayesian inversion methods toward large-scale problems is highly growing. In recent years, several variants of the MCMC algorithms have been proposed in the mathematical and geophysical literature. Beginning from the Metropolis-Hastings algorithm and Gibbs method (Metropolis and Ulam 1949; Hastings 1970; Geman and Geman 1984; Gelfand and Smith 1990; Gelfand et al. 1990), different investigations on enhancing the efficiency of the MCMC algorithms have been implemented for better performance of the classical approaches. For instance, Haario et al (2001) proposed a self-tuning algorithm namely the adaptive metropolis algorithm in which proposal values are sampled from a multivariate normal distribution with covariance matrix generated using the accepted samples of the chain. However, it was shown by Cui et al (2011) that the proposal distribution used by the adaptive Metropolis method can be suboptimal in cases where the posterior distribution is non-Gaussian. Later, they suggested the delayed rejection adaptive Metropolis method (Haario et al, 2006). Ter Braak (2006) introduced the differential evolution Monte Carlo algorithm. Vrugt et al (2008) applied the idea of combining the differential evolution technique and the adaptive Metropolis to hydrological data. An accelerated variant of the MCMC technique is the parallel tempering algorithm (Swendsen and Wang, 1986) in which multiple chains are simulated but are allowed to swap information. This method can be much more efficient than the classic Metropolis-Hastings algorithm for problems involving multi-modal posterior distributions (Higdon et al. 2002; Dettmer and Dosso, 2012; Sambridge, 2013). Recent examples of parallel tempering as applied to geophysical inversion can be found in Dosso et al. (2012), Ray et al (2013), Blatter et al (2018), and Blatter et al (2021). Despite significant improvements in MCMC sampling methods, choosing an appropriate proposal distribution, used to generate trial moves in the Markov chain, is of crucial importance for an efficient inversion process which creates a well-mixed Markov-Chain, that is, chain samples widely over the model parameters space avoiding both high rejection rates and small steps. Note that the mixing of a Markov chain is the number of steps the Markov chain must take before its probability distribution reaches the stationary distribution (Andrieu and Thoms, 2008). The common choice for proposal density is a Gaussian distribution centered on the current model and variance tuned by the user

کلیدواژه‌ها [English]

  • Test