Joint geophysical, petrophysical and geologic inversion using a dynamic Gaussian mixture model

by   Thibaut Astic, et al.

We present a framework for petrophysically and geologically guided inversion to perform multi-physics joint inversions. Petrophysical and geological information is included in a multi-dimensional Gaussian mixture model that regularizes the inverse problem. The inverse problem we construct consists of a suite of three cyclic optimizations over the geophysical, petrophysical and geological information. The two additional problems over the petrophysical and geological data are used as a coupling term. They correspond to updating the geophysical reference model and regularization weights. This guides the inverse problem towards reproducing the desired petrophysical and geological characteristics. The objective function that we define for the inverse problem is comprised of multiple data misfit terms: one for each geophysical survey and one for the petrophysical properties and geological information. Each of these misfit terms has its target misfit value which we seek to fit in the inversion. We detail our reweighting strategies to handle multiple data misfits at once. Our framework is modular and extensible, and this allows us to combine multiple geophysical methods in a joint inversion and to distribute open-source code and reproducible examples. To illustrate the gains made by multi-physics inversions, we apply our framework to jointly invert, in 3D, synthetic potential fields data based on the DO-27 kimberlite pipe case study (Northwest Territories, Canada). The pipe contains two distinct kimberlite facies embedded in a host rock. We show that inverting the datasets individually, even with petrophysical information, leads to a binary geologic model consisting of background or kimberlite. A joint inversion, with petrophysical information, can differentiate the two main kimberlite facies of the pipe.


page 7

page 18

page 19

page 21

page 22

page 24

page 32


Petrophysically and geologically guided multi-physics inversion using a dynamic Gaussian mixture model

In a previous paper, we introduced a framework for carrying out petrophy...

An Inversion Tool for Conditional Term Rewriting Systems – A Case Study of Ackermann Inversion

We report on an inversion tool for a class of oriented conditional const...

A comparative study of structural similarity and regularization for joint inverse problems governed by PDEs

Joint inversion refers to the simultaneous inference of multiple paramet...

A probabilistic approach to emission-line galaxy classification

We invoke a Gaussian mixture model (GMM) to jointly analyse two traditio...

Efficient learning methods for large-scale optimal inversion design

In this work, we investigate various approaches that use learning from t...

Multi-Regularization Reconstruction of One-Dimensional T_2 Distributions in Magnetic Resonance Relaxometry with a Gaussian Basis

We consider the inverse problem of recovering the probability distributi...

Simultaneous shot inversion for nonuniform geometries using fast data interpolation

Stochastic optimization is key to efficient inversion in PDE-constrained...

1 Introduction

Mineral deposits, or other geologic features, are characterized by different physical properties, and hence multiple geophysical surveys are used to delineate them. For example, kimberlites have signatures that depend upon density, magnetic permeability, and electrical conductivity. They are often discovered through data collected in airborne surveys thanks to their signature of a circular low gravity anomaly, high magnetic response, and sometimes negative electromagnetic transient response (Macnae, 1995; Keating and Sailhac, 2004; Bournas et al., 2018). Although a joint interpretation of several individually inverted datasets can significantly improve our understanding of the subsurface (Oldenburg et al., 1997; Devriese et al., 2017; Kang et al., 2017; Melo et al., 2017), multiple cases studies have shown that joint inversions can reveal information that was not accessible through individual geophysical dataset inversions (Doetsch et al., 2010; Jegen et al., 2009; Kamm et al., 2015). However, those inversions require a coupling term that uses a priori information about the physical properties or geologic structures. Coupling methods for joint inversion schemes generally use one or a combination of those two types of information.

The firsts frameworks for joint inversion focused on linking geophysical models through their structural similarities. Haber and Oldenburg (1997) defined the structure of a model in terms of the absolute value of its spatial curvature and compared different models to see if variations occurred at the same locations. This idea was further developed by Gallardo and Meju (2003) with the introduction of the concept of cross-gradient. This approach has become commonly used, and Gallardo and Meju (2011) provides an in-depth review of the method. However, this strategy assumes that physical properties contrasts are co-located, and that is not valid for many problems. This drawback can be overcome by using other coupling methods.

The second coupling approach uses physical property relationships to link geophysical models. The framework proposed in this paper belongs to this family. Some of the earliest works used an experimental constitutive formula as their physical properties constraint (Afnimar et al., 2002; Onizawa et al., 2002). Stefano et al. (2011) combined this approach with the above mentioned cross-gradient method for sub-salt imaging. Some stochastic frameworks have also been proposed. Grana and Della Rossa (2010); Grana et al. (2017) proposed a Bayesian approach for the prediction of lithological facies and fluid properties from seismic attributes. Shamsipour et al. (2012) used geostatistical techniques of cokriging and conditional simulation for jointly inverting gravity and magnetic data assuming the density and magnetic susceptibility auto- and cross-covariances follow a linear model of coregionalization. On the deterministic side, of which the framework we present belongs, recent frameworks used clustering techniques such as the fuzzy C-means algorithm, which was first used in Paasche and Tronicke (2007) and further expanded in Lelièvre et al. (2012). This approach adds a clustering term to the objective function. It allows to formulate less strict relationships between physical properties. An important contribution in the last few years is the work by Sun and Li (2015, 2016, 2017) which further developed the method. Besides the addition of the fuzzy C-means term, they developed an iterative update to the cluster centres throughout the algorithm; this starts introducing the notion of uncertainty for the petrophysical data. Giraud et al. (2017) focused on reducing uncertainties in stochastic geological modelling by linking geophysics and geological modelling through petrophysical information.

This study extends the framework for petrophysically and geologically guided inversion (PGI) developed in Astic and Oldenburg (2019) as a coupling term to perform joint inversions. This framework (Fig. 1) generalizes concepts presented in Grana et al. (2017), Giraud et al. (2017) and Sun and Li (2017)

. We show how geological and petrophysical information represented as a Gaussian mixture model (GMM) can be used to couple multiple voxel-based geophysical inversions. Each contrasting geological unit is represented by a multi-dimensional Gaussian distribution which summarizes its characteristics and petrophysical signature. The resulting objective function incorporates both petrophysical and geological information into a single smallness term in the regularization, and the iteration steps are decomposed into a suite of cyclic optimization problems across the geophysical, petrophysical and geological data.

In this paper, we start by reviewing the key concepts of PGI. We then generalize those concepts for the case of multiple physical properties and governing equations for performing joint inversions. We delineate the numerical implementation considerations that are necessary when dealing with multiple datasets and our strategy for reaching multiple target misfits. We also highlight our flexible and modular implementation within the SimPEG framework (Cockett et al., 2015; Heagy et al., 2017). Finally, we demonstrate the advantages of performing joint inversion by using a synthetic model of the DO- Tli Kwi Cho kimberlite, Northwest Territories, Canada (Jansen and Witherly, 2004).

Figure 1: A graphical representation for our PGI framework. Each diamond is an optimization process that requires data (shown in rectangular boxes) as well as prior information provided by the other processes.

2 Key concepts for PGI joint inversion

In this section, we present our representation of petrophysical and geological information as a Gaussian mixture model (GMM) and its inclusion into a geophysical inverse problem.

2.1 Notation conventions

In terms of notation for parameters, we use the following conventions:

  • Lowercase italic symbols are used for scalar values, such as the trade-off parameter .

  • Bold lowercase symbols designate vectors, such as the geophysical model


  • Bold uppercase symbols designate matrices, such as the multidimensional Gaussian distribution covariance matrix .

For a joint inversion, the geophysical model is likely to contain multiple physical properties. Several surveys might point to the same physical property (e.g. gravity and gravity gradiometry) or one survey might take several physical properties (e.g. electromagnetic surveys). We thus adopt the following notations for the geophysical model indices:


A row of the matrix represents all the physical properties that live in the same location. A column represents a single physical property at all the mesh cells. For clarity, we stay consistent throughout this study with the index notation. The index always refers to the number of model’s parameters, from 1 to . We will then denote a single cell’s physical properties as:


Likewise, we denote the vector model for a single physical property on the whole mesh with the index from 1 to with an exponent:


Lastly, the geophysical model at iteration of an inversion will be denoted with parentheses .

2.2 Background: the Tikhonov inverse problem

Our inverse problem can be considered as an augmentation of the well-established Tikhonov inversion, first introduced in (Tikhonov and Arsenin, 1977). This is labeled process in our PGI framework (Fig. 1). In the Tikhonov approach, the geophysical inverse problem is cast as an optimization problem in which an objective function , such as shown in equation 5, is minimized. Using the same notation convention as Oldenburg and Li (2005), the goal is to find a solution that minimizes:


In equation 5, the vector is the geophysical model, which represents physical properties on a mesh. The term is the data misfit, is the model regularization function and is a positive scalar that adjusts the relative weighting between the two terms. A value of is sought so that the data misfit is below an acceptable target misfit (Parker, 1977). The regularization term contains a priori information about the expected features of the geophysical model. It is usually composed of a smallness term that regulates the values the model can take and smoothness terms in all directions which control the variations between adjacent locations, weighted by scalar weights :


The smallness term is essential for our framework. Next, we present our representation of petrophysical and geological information and how we include it into the regularization term.

2.3 Gaussian Mixture Model for petrophysical and geological information

To include physical properties information into the inversion, we need to define a way to model the distribution of the physical properties mathematically. Consider a petrophysical dataset with rock samples. For each sample, we have measured physical properties. We will denote those measurements for the various samples (indexed by ) by the vectors . We assume that for each rock sample of the dataset we know its current geological classification, denoted by where is the number of distinct geological units. A sample belonging to unit is noted .

Figure 2: Example of a two-dimensional Gaussian mixture model with three rock units. The background is coloured according to the geological classification evaluated by equation 11. Thicker contour lines indicate a higher probability density. On the left and bottom, we provide the 1D projections of the total and individual probability distributions for each physical property.

We represent the petrophysical signature of each geological unit as a multidimensional Gaussian probability distribution. If the physical properties of the unit do not follow a Gaussian distribution, we can often project the physical properties into a transformed space where their distribution appears approximately Gaussian, such as a log-space for electrical conductivity or magnetic susceptibility (see the notion of mapping in the SimPEG framework as defined in Cockett et al. (2015); Kang et al. (2015)). The Gaussian distribution representing the physical properties distribution for each unit is defined by its mean (vector of size ), and its covariance matrix (matrix of size ), plus its proportion . Their respective definitions are given in equations 7 to 9, with labeled samples for each unit (total of samples):


The expected probability function to observe an unlabeled physical property data point can be written as a Gaussian Mixture Model (GMM), which simply sums the probability of each known unit:


The variable holds the GMM global variables . With some modifications, it can also represent nonlinear relationships (such as polynomial as presented in Onizawa et al. (2002)); see Appendix B for such an example.

We denote

as the categorical variable for the labels, i.e. the geological classification (or membership). Given an unlabeled rock sample with physical properties

, we define its membership as the most probable geological unit:


The categorical variable is key for building a quasi-geological model (Li et al., 2019) from the physical properties model obtained by inversion. This corresponds to the process in our framework (Fig. 1).

In Fig. 2, we present an example of a GMM with three distinct rock units characterized by two physical properties (two-dimensional distributions). The background is coloured according to the geological identification that would be made at each location using equation 11. On the bottom and the left of the Fig. 2 we represent the GMM probability distribution for each property individually. We notice that, while all three units are distinct in the two-dimensional space, they can overlap significantly when only considering one property at the time. For property , the rock units and are distinct while rock unit is indistinguishable from rock unit . For property , rock units and are now distinct while unit is indistinguishable from either rock units or . This highlights the gains that can be made by jointly inverting for several physical properties. Units that might not be distinguishable in one survey may be in another one, and this distinction will guide the inversions towards separating them in both physical properties. In the following section, we show how to use this probability distribution that links the various physical properties as a priori information to regularize the multi-physics inverse problem.

2.4 Gaussian mixture model definition

Suppose we are given petrophysical and geological information about the geophysical model that we want to recover in our inversion. Petrophysical information can include an expected number of distinct units, the mean physical property value for each geological unit, or their covariance. Geological information can consist of an expectation of encountering certain rock units in particular locations. We need to design a probability distribution that offers the maximum flexibility for representing that petrophysical and geological a priori knowledge.

In Astic and Oldenburg (2019), we developed a GMM prior to include petrophysical and geological information in inversions involving only one type of geophysical survey with only one physical property to recover. In the equation 14

below, we propose an extended version of that GMM prior that is designed to couple multiple physical properties as well as geological information. Note that now the means are vectors the size of the number of different physical properties and the scalar variance becomes a full positive-definite matrix of the same size. The parameters are both spatially (index

) and lithologically (index ) dependent.



  • is the number of distinct rock units.

  • is the number of active cells in the mesh.

  • represents the physical property values at location .

  • describes the membership to a rock unit.

  • is the a priori probability of observing rock unit at location . It can be either constant over the whole area, then denoted by , or locally determined by a priori geological knowledge.

  • contains the means of the physical properties of rock unit .

  • is the covariance matrix of the physical properties of rock unit .

  • is a weighting term at location , used for example to include depth or sensitivity weighting. We define it from weights ; this is a scalar value for each cell and physical property (with different physical properties). is defined as a diagonal matrix made of the combination of all the weights at the specific location : , with the same notation convention as for the model . This allows the weighting to be different for each physical property. We can thus weigh each physical property accordingly to the survey on which it depends.

  • holds the GMM global variables .

This GMM probability distribution (equation 14), representing the current geological and petrophysical knowledge, is used to define the smallness term in the regularization. In the next subsection, we develop the modified objective function for the inverse problem.

2.5 PGI Joint Inversion

In Astic and Oldenburg (2019), we demonstrated how using the GMM defined in equation 14 modifies the smallness term of the Tikhonov inverse problem. The smallness term is still being defined as a least-squares misfit between the current model, but the reference model is now updated at each iteration, as is the smallness matrix . This dynamic reference model and smallness matrix are determined based on the current geophysical model and the geological and petrophysical prior information. Generalizing the result obtained in Astic and Oldenburg (2019) to multiple physical properties, we obtain a single smallness term that encompasses and couple all of the model parameters:


When the petrophysical signature of the rock units are not known, it is possible to learn the parameters of the GMM (Astic and Oldenburg, 2019) (Fig. 1, process ). The extension to multiple physical properties of that learning process can be found in Appendix A.

2.6 Petrophysical and geological target misfit

The PGI smallness expresses a misfit between the petrophysical and geological information and the geophysical model. In Astic and Oldenburg (2019) we defined a measure , similar to the PGI smallness term but without the weights , for which we defined a target misfit . We can generalize the same approach here:


Each term in equation 19 is a multivariate Gaussian variable of dimension (the number of physical properties) of mean with an identity covariance. Thus,

follows a chi-squared distribution with a sample size equal to the number of model parameters. Thus its target misfit value is defined as its expectation, such as defined in

Parker (1977) for a geophysical data misfit, and is given by:


with the number of active cells in the mesh and the number of physical properties. This generalizes the result obtained in Astic and Oldenburg (2019) ().

Our algorithm stops when all the target misfits, geophysical and petrophysical, are fitted. In the next section, we present our approach for handling multiple geophysical data misfits as well as an additional petrophysical misfit.

3 Numerical considerations for reaching multiple target misfits

We have multiple geophysical data misfits that we wish to fit. In addition, the inclusion of petrophysical and geological data with PGI adds another data misfit term; this also needs to reach its target misfit. In this section, we provide our strategies to choose and dynamically adjust the various parameters of the objective function to find a solution that fits all the data. An algorithm that summarizes the whole framework is provided in Appendix D.

3.1 Notations and objective function for multiple geophysical data misfits

We present in the equation below the objective function we are minimizing during the geophysical inversion process (Fig. 1, process 1). The terms contains multiple geophysical data misfits, each defined as a least-squares norm. The smallness term is our coupling defined in equation 15. The smoothness terms, one for each direction and physical property, are contained within :


where is the data misfit of the survey and the scalars are weighting terms for each data misfit. These become important to balance the various terms so that a solution is sought such that each is below its target misfit . This is developed later in this section. The forward operator generates the predicted data for the survey from . The notation denotes the subset of model parameters associated with the survey. Note that multiple surveys can be associated with the same physical property, for example, gravity and gravity gradiometry both depend on density contrasts. Similarly one survey can be sensitive to several physical properties; for example electromagnetic surveys are sensitive to electrical conductivity and magnetic-susceptibility. The data measured by the survey is symbolized by and the uncertainty on those measurements by the matrix . In the smoothness terms (equation 26), the smoothness operator (usually first or second order difference) is represented by the matrix .

3.2 Review of previous strategies for balancing various geophysical data misfits and coupling terms

In this section, we present strategies that have been previously published with regards to balancing multiple geophysical data misfits along with coupling terms.

In the context of a single geophysical data misfit, we developed, in Astic and Oldenburg (2019), our empirical approach to fit a geophysical data misfit, as well as a petrophysical and geological data misfit. Our approach was based on updating the and values at each iteration. Similar work for linear problems has also been detailed in Sun and Li (2018) where the authors explore several parameters to optimize the clustering of their model, they do not define a target misfit for the petrophysical and geological information.

When several geophysical methods are involved, we need a strategy to scale the multiple geophysical data misfits so that each is properly fitted. This is where the weights in equation 25 become important. Several approaches have been published in the recent literature. Sun and Li (2016) and Sosa et al. (2013) normalize each geophysical survey data misfit by its number of data and keep those weights constant. In our experience, keeping the weights constant has led us to overfit some surveys while underfitting others. On the other hand, Lelièvre et al. (2012) devised a rigorously defined, but computationally expensive, strategy for balancing two geophysical data misfits. Their approach relies on adjusting the trade-off parameter until the two data misfits are Pareto-optimal. They then adjust the relative weights of the two surveys to fit both geophysical surveys. They next reinforce the importance of the coupling term before going into another round of adjustments of the trade-off and surveys weights parameters. Their approach is somewhat related to what we developed in Astic and Oldenburg (2019), in the sense that they focus first on fitting the geophysical data misfit terms and then adjust the coupling term. On the contrary, the strategy presented in Gallardo and Meju (2004) favoured the coupling over the geophysical misfits.

Here our goal is to define a practical, computationally inexpensive, strategy for balancing the geophysical data misfits as well as the coupling term. We design this strategy to work for any number of surveys and to not require additional simulations at each parameter update. For the full algorithm, the reader can refer to Appendix D.

3.3 Scaling the regularization terms

Two types of weighting act on the regularization terms; they are the trade-off parameter and the parameters. In Astic and Oldenburg (2019), we developed a strategy for cooling and warming to find a solution to the inverse problem that fits the geophysical and the petrophysical and geological information. This philosophy is still appropriate in a joint inversion framework and is what we use in this study.

For the smoothness terms, we use the weights to scale each physical property representation by its expected maximum amplitude (available through the GMM means if provided). This helps to equalize their contribution to the objective function value when parameters have widely different scales (like density contrast and magnetic susceptibility contrast). Note that the scaling of the physical properties in our extended smallness term (equation 15) is taken care of by the covariance matrices of the GMM.

3.4 Balancing the geophysical data misfits

We now define our strategy to scale the multiple geophysical data misfits to reach all the target values. We use the weights defined in equation 25. We dynamically update each weight based on its current misfit value as compared to the respective targets for the other surveys.

For initialization, a common normalization is to use the number of data as is done in previous strategies presented above. Here, we also propose a novel initialization. Our motivation is to not only account for the number of data but also for the relative noise levels, the scale of the physical properties, and the numerical operators, which all play major roles in the data misfit metrics. For example, density contrasts are often (not always) within while magnetic-susceptibility is most often below SI. As for the associated operators, a given body with a density contrast of can create a gravity response of a few mGal while the same body with a magnetic-susceptibility of

SI will produce a magnetic response in the hundreds of nT. There are two steps to our initialization. First, we account for the range of predicted data and noise levels of each geophysical survey, and second, we address the variability in scales of the physical properties. In the first part, for each geophysical data misfit term, we estimate the highest eigenvalue, denoted

, of , with the linearization (or sensitivity matrix, see McGillivray and Oldenburg (1990)) of the physical operator (equation 27). We choose to use a Rayleigh quotient with a power method (Ascher and Greif, 2011), which allows us to get a quick approximation of the . To account for the scale of physical properties, we multiply the highest eigenvalue estimate by the maximum expected value of the associated model parameters (available through the GMM means). The weight is defined as the inverse of the resulting product (equation 28). Because we will update those relative weights later, we prefer to always normalize the sum of the weights to unity (equation 30). This is to keep the importance of the total term relatively similar before and after adjusting the weights.

then we normalize the sum: (29)

This approach worked well on synthetic examples with linear operators (such as gravity and magnetic). This is what is used in the synthetic example in the present paper.

Good initial weights can help the inversion converge faster, but they have not appeared crucial in our experiments. What ensures the progress of all data misfits, while limiting the possibilities of overfitting, is to update the weights during the inversion. Our approach is philosophically similar to what we use for and to balance the geophysical data misfit and the petrophysical and geological misfit at each iteration in the inversion. If one geophysical data misfit reaches its target value before the others, we use the ratio of its current value with its target to warm the weights of the other geophysical data misfit terms and then normalize the sum of the weights to be equal to unity again. If several surveys are below their respective targets, we simply use the median of the ratios to warm the weights of the still unfitted surveys. Thus at any iteration of the geophysical inverse problem, if an ensemble of surveys has reached their respective targets, we warm the weights of the remaining surveys that are not yet fitted as:


This strategy is what we used in the case study in Astic et al. (2020) where we invert together three field datasets, including complex petrophysical and geological relationships. The influence of the initial weights has proven minimal compared to the reweighting strategy for the fitting of all datasets.

3.5 A flexible interface for joint inversions in SimPEG

We implemented our framework as part of the open-source software SimPEG (Cockett et al., 2015; Heagy et al., 2017). SimPEG is designed to be a modular, extensible, object-oriented framework for simulations and inversions of geophysical data. In particular, two important features enabled us to focus our implementation efforts:

  1. the composability of objective functions in the data misfit and the regularization terms,

  2. the use of directives which orchestrate updates to components of the inversion at each iteration.

The first point streamlines the implementation of joint inversions. In the code, each misfit term is a Python object that has properties, such as the weights used to construct , and methods, including functions to evaluate the misfit given a model as well as derivatives for use in the optimization routines. To construct a joint inversion with gravity and magnetic data, we first define each misfit term independently and then simply sum them. We use operator-overloading in Python so that when we express the addition of two objective functions in code, the evaluation of this creates a combo-objective function. This is an object that has the same evaluation and derivative methods as the individual data misfits, and thus readily inter-operates with the rest of the simulation and inversion machinery in SimPEG. In pseudo-code, the construction of the data misfit as expressed in equation 25 is as follows:

default where the -values are scalar weights.

The composability and interoperability extend to the construction of the regularization term. We can build each part of the regularization before combining them:


Having built our framework within the SimPEG combo-objective function also means that our PGI approach can be readily incorporated into any of the geophysical survey types that are supported in the source code, including electrical and electromagnetic methods (Astic and Oldenburg, 2019) as well as magnetic-vector inversions and gravity gradiometry (Astic et al., 2020). This makes the construction of the objective function for any specific joint problem relatively easy for the user.

To the second point, directives are functions that are evaluated at the beginning or end of each iteration in the optimization. They are the mechanism we use to make updates to components of the inversion, including the data misfit weights (equation 31) or smallness weights and reference model (equations 17 and 18), as well as for evaluating the target misfits and stopping criteria for the inversion.

Being open-source, working within the SimPEG framework allows us to be transparent about our implementation and to distribute the base code. This also offers us a unique opportunity to provide reproducible examples. Scripts and Jupyter notebooks to reproduce the examples presented in this study are available through Github at (Astic, 2020).

4 Example: The DO- kimberlite pipe

In this section, we present an illustration of how a joint PGI approach of several geophysical surveys, along with petrophysical data, can improve upon the Tikhonov inversions and the individual PGI results of each geophysical dataset. For that purpose, we compare 3D inversions results of synthetic gravity and magnetic data, based on the DO- kimberlite pipe, composed of two different facies. We show that for this case study, the geophysical models obtained from Tikhonov inversions, or even from a PGI approach for each survey, only lead to a binary resolution: host rock or kimberlite pipe. We then demonstrate that using this full joint PGI framework allows distinguishing two clear, distinct kimberlite facies from the background host rock.

4.1 Setup

Figure 3: Setup: TKC synthetic example: (a) Three cross-sections through the synthetic geological model. The dots represent the data locations for the gravity and magnetic survey; (b) Cross-plot and histograms of the physical properties of the synthetic model; (c) Synthetic ground gravity data; (d) Synthetic total amplitude magnetic data.

To illustrate the gains obtained with the joint PGI approach, we use simulated surface gravity and airborne magnetic data modelled from a synthetic model of the DO- Kimberlite pipe (Northwest Territories, Canada), part of what was formerly known as the Tli Kwi Cho (TKC) kimberlite deposit (Jansen and Witherly, 2004) (Fig. 3). The pipe has two distinctive kimberlite units of interest that are embedded in a background consisting of a granitic basement covered by a thin layer of till (Fig. 3a). The first unit is a pyroclastic and volcanoclastic kimberlite (called PK/VK), which has a weak magnetic susceptibility () and a very high density contrast (). The second unit is an hypabyssal kimberlite (called HK) which has a strong magnetic susceptibility () and a weak density contrast () (Fig. 3b).

The inversion mesh is a tensor mesh with 375 442 active cells; each has a pair of density-magnetic susceptibility values. The smallest cells are

m in all directions. We forward modelled the data over a grid of 961 receivers at the surface (Fig. 3

c and d). We added Gaussian noise to the gravity and magnetic data with standard deviations of

mGal and nT everywhere. These values are input into the data weighting matrices .

For each inversion, we added bound constraints so that the sought density contrast values are null or negative only, and the magnetic susceptibility contrast values are null or positive. We used the sensitivity of each survey to define the weights. Sensitivity-based weighting is a common practice for potential fields inversions (Li and Oldenburg, 1996, 1998; Portniaguine and Zhdanov, 2002; Mehanee et al., 2005). Note that the sensitivity of the gravity survey is used to weigh the density contrast model, while the sensitivity of the magnetic survey is used to weigh the magnetic-susceptibility model. The initial model is the background half-space for all inversions.

4.2 Tikhonov inversions

Figure 4: TKC gravity and magnetic Tikhonov inversion results. (a) Three cross-sections through the recovered density model; (b) Three cross-sections through the recovered magnetic susceptibility model; (c) Cross-plot of the two Tikhonov inversions of the gravity and magnetic data respectively. The points are coloured using the density and the magnetic-susceptibility contrast values (white, BCKGRD: background), density contrast only (blue), magnetic-susceptibility contrast only (orange), both density and magnetic susceptibility contrast (purple)). The same colourmap is used in (d); (d) Three cross-sections coloured based on the combination of density and magnetic-susceptibility contrast recovered by Tikhonov inversions.

We first run the individual inversions of the gravity and magnetic data using the well-established Tikhonov approach. The results, shown in Fig. 4, are relatively smooth. The gravity inversion (Fig. 4a) provides an approximate outline of the pipe. The magnetic inversion, on the other hand, shows a magnetic body centred on HK, but it is too diffuse to delineate a shape (Fig. 4b). Looking at the cross-plot (Fig. 4c), we observe the expected continuous Gaussian-like distribution of the model parameters (in the region allowed by the bound constraints). Petrophysical signatures are not reproduced, with notably the highest density contrast values being associated with the highest magnetic susceptibility values. In both cases, the two anomalous units are indistinguishable from each other. To highlight this, we show an overlap of the two inversions in Fig. 4(d). We coloured each point relative to its density - magnetic susceptibility values (also shown in Fig. 4c); white for both null values, blue-scale for only a density contrast, red-scale for a magnetic susceptibility only and a purple-scale for both). This highlights that combining both models does not show structures that seem closer to the ground truth. Thresholding would give highly variable results, depending on the values chosen to delineate units.

We then move to a PGI approach, by including petrophysical information, to invert each geophysical dataset individually. We assert what gains are made before moving to a full joint PGI.

4.3 Individual PGI

Figure 5: Results of the individual PGI. (a) Three cross-sections through the density model recovered using the petrophysical signature of PK/VK; (b) Cross-sections through the magnetic susceptibility model obtained using the HK petrophysical signature; (c) Cross-plot of the inverted models. The 2D distribution in the background has been determined by combining the two 1D distributions used for density and magnetic susceptibility PGI, respectively. With only one anomalous unit in each case, this still gives four possible combinations; (d) Cross-sections through the quasi-geological model built from the density and magnetic susceptibility models, see cross-plot.

We apply the framework developed in Astic and Oldenburg (2019) to invert each geophysical dataset individually with the PGI approach. Results are shown in Fig. 5. For the prior petrophysical distribution, we use the true value for the means and proportions of each unit. For the petrophysical noise level, we assign standard deviations of of the highest known mean value for each physical property except for the background for which we assign (see the 1D contour plots on the left and the bottom for each physical property respectively in 5d). The GMM parameters are held fixed.

With that approach, we found that both magnetic and gravity data can be explained, individually, by assuming a single unit, either PK/VK or HK. For concision, we choose to show here the gravity result recovered using only the petrophysical signature of the PK/VK unit (Fig. 5a), which is the most responsible for the gravity anomaly; for magnetic inversion, we show the model obtained by only using the petrophysical signature of HK, which is the unit that is the most responsible for the magnetic response (Fig. 5b). The additional models (gravity inversion with HK density signature and magnetic inversion with PK/VK magnetic signature) are shown in Appendix C where we also demonstrate the discrepancy in the recovered volumes of each unit between the magnetic and gravity inversions. The recovered volumes of each unit across inversions are not in agreement. For example, explaining the gravity anomaly with only a body with the same density as HK leads to a very large body (as PK/VK density contrast is four times more important), bigger than the volume of that same unit recovered through the magnetic inversion. The same reasoning applies to the PK/VK body.

The gravity PGI using the PK/VK unit petrophysical signature (Fig. 5a) gives useful information about the depth and delineation of the pipe, more than the Tikhonov inversion. The magnetic PGI using the HK petrophysical signature (Fig. 5b) places a body around the HK unit location, but miss its elongated shape. From the petrophysical perspective, both the gravity signature of PK/VK and the magnetic signature of HK are individually well reproduced, but not their combination. The density-magnetic susceptibility cross-plot for the individual PGI (Fig. 5c) is very far from the desired petrophysical distributions. Assuming even just one anomalous unit for each inversion as we did, there are still different combinations of physical properties. In this specific case, there is a clear anomalous body that has both the density signature of PK/VK and the magnetic signature of HK (Figs 5c and d), that we identify as ’undefined’ in the figure. It does not correspond to any unit signature and occupies a large volume. This hinders our ability to resolve two clear kimberlites from the inversions. Therefore we move to a joint inversion approach to take advantage of the density-magnetic susceptibility relationship in the inversion and finally delineate two distinct kimberlite facies.

4.4 Joint PGI

Figure 6: Results of the joint PGI approach. (a) Three cross-sections through the recovered density contrast model; (b) Cross-sections through the magnetic susceptibility contrast model; (c) Cross-plot of the inverted models. The colour has been determined by the clustering obtained from this framework joint inversion process. In the background the prior joint petrophysical distribution; (d) Cross-sections through the resulting geological model from the joint inversion process.

We now apply our extended PGI framework to the joint inversion of the gravity, magnetic and petrophysical data all-at-once. The parameters are again used to include the proper sensitivity weighting for each method. We use the same uncertainties that we used for the individual PGI. The non-diagonal elements of the covariance matrices are set to null, which just means we assume no correlations between the density and magnetic-susceptibility variations within a single rock unit. The GMM parameters are held fixed. The results are shown in Fig. 6.

The improvement is significant. The joint inversion succeeds in recovering two distinct kimberlite facies that reproduce the provided petrophysical signatures. The surface outline is well recovered but now with two distinct facies, contrary to the individual PGI. The vertical extension is comparable to the truth. We also get a sense of the elongated shape and tilt of the HK unit (Fig. 6d). From the petrophysical perspective, the signature of each unit is well reproduced. The HK magnetic-susceptibility of HK is slightly overestimated but still within the acceptable margins defined by the petrophysical noise levels we assigned (Fig. 6c). The discrepancy in the volume of each unit observed for the individual PGI has been resolved by providing the algorithm with the full 2-dimensional petrophysical distributions.

4.5 Comment on joint inversions of potential fields data with limited petrophysical information

Figure 7: Results of the joint PGI without providing the means of the physical properties for the kimberlite facies. (a) Three cross-sections through the recovered density contrast model; (b) Cross-sections through the magnetic susceptibility contrast model; (c) Cross-plot of the inverted models. The colour has been determined by the clustering obtained from this framework joint inversion process. In the background the prior joint petrophysical distribution; (d) Cross-sections through the resulting geological model from the joint inversion process.

In Astic and Oldenburg (2019), we emphasized the benefits of learning a GMM during the inversion process to compensate for uncertain or unknown petrophysical information. This allows running our PGI framework without providing physical property mean values. This is performed by the petrophysical characterization process (Fig. 1), which updates the GMM parameters at each iteration based on the current geophysical model. While we generalize the learning process of the GMM parameters to a multidimensional case in Appendix A, we have not made use of it in the present study so far. In our current understanding, it is not as efficient for potential fields as it appeared for electromagnetics inversions. Our hypothesis to explain this is related to the multiple illuminations of the anomalous bodies, in addition to the nonlinearity of the problems, that electromagnetics methods offer compared to potential fields methods. The volume over contrast ratio is the quantity of importance in potential fields. The four individual PGI provide an excellent example of it. Each dataset, gravity or magnetic, could be fitted by either reproducing the signature of PK/VK or HK. The difference in physical property contrast was compensated by a difference in the volume of the recovered body.

We illustrate this point by running a joint PGI without providing physical property mean values nor proportions for the kimberlite units. The means and proportions of the kimberlite units are learned iteratively through the inversion. This is done by setting the and parameters to null for the kimberlite facies (see Appendix A). We fix the mean of the physical properties for the background to its true value (null). We keep the covariances similar to what we used previously. They define our petrophysical noise and how spread each petrophysical signature can be.

We find a volume made of a single kimberlite facies, with physical properties within the prescribed ranges, that fit both the gravity and magnetic data. The resulting geophysical model, learned GMM and geological model are shown in Fig. 7. It thus appears to us that, contrary to the electromagnetics examples provided in Astic and Oldenburg (2019), a good enough guess of the physical properties signature is necessary to perform reasonable PGI of potential fields data. We direct the reader to our DO- field case study (Astic et al., 2020) for an example. This aspect will require more thorough research.

4.6 TKC example summary

In all of the inversions, we have shown here, the near-surface features such as the east tip of the PK/VK unit are relatively well recovered. From the petrophysical perspective, the density-magnetic susceptibility cross-plots for the standard Tikhonov inversions are very removed from the expected distributions represented in contours. Individual PGI give us more information about the depth and delineation of the main PK/VK body. Individual datasets can both be reproduced using a single kimberlite facies. The geological models are, however, incompatible once put together and the petrophysical signature of each unit is only partially recovered. Finally, the joint inversion distinguishes the two units in the recovered models and even starts to give us information about the elongated shape and tilt of the HK unit. A correct estimation of the physical properties contrast appears to be essential for the success of the inversions of potential fields data with PGI.

5 Discussion

We have expanded our PGI framework to carry out joint inversions. We have focused on the methodology and demonstrating the capabilities of the PGI joint framework on an example whose ground-truth is known. A challenging case study that jointly inverts three geophysical surveys (gravity, gravity gradiometry and magnetic-vector inversion) with five model parameters can be found in Astic et al. (2020)). We have so far produced examples of joint linear problems. In our previous work (Astic and Oldenburg, 2019)

, we applied the PGI approach to multiple electromagnetics nonlinear problems (magnetotelluric, direct-current resistivity and a field RESOLVE dataset (frequency-domain electromagnetics)). We plan to implement this approach for performing multi-physics inversions with electromagnetic methods.

Oldenburg et al. (2019) demonstrated a joint Tikhonov-style inversion of several electromagnetics surveys also using the SimPEG package, but without physical property information and with no specific strategy for individually fitting the various data misfits.

We proposed a strategy to balance any number of geophysical data misfits along with a coupling term. In our experiments, this strategy appeared to be critical to fitting multiple surveys as well as petrophysical data. On the contrary, the initial misfit weights played a minor role in comparison. The robustness of our re-weighting strategies for joint inversions remains to be tested with nonlinear geophysical problems.

The extension to several physical properties also offers the opportunity to include relationships that link physical properties together. Linear relationships are straight-forward to implement with Gaussian distributions through the covariance matrix, which can define tilted, elongated probability distributions. An example of such a linear trend can be found in our case study (Astic et al., 2020). In Appendix B, we discuss how to account for nonlinear relationships. Such nonlinear relationships are found, for example, between density and seismic velocities (Onizawa et al., 2002). Our framework is flexible enough that different relationships can be included for each rock unit. While very modest, our example shown in Fig. 9

is, to our best knowledge, the first one in the literature with such diverse relationships in a single inversion. For the moment, our framework assumes that those nonlinear relationships are given. An interesting avenue of research would be to develop the mathematics for the learning of those nonlinear relationships, along with the other GMM parameters (such as defined in Appendix

A). While we think it might be possible for nonlinear geophysical problems with multiple sources, our first assumption is that it might encounter similar issues like the ones mentioned above for the learning of the GMM parameters in potential fields inversion.

Our PGI framework is composed of three regularized optimization problems (Fig. 1). In Astic and Oldenburg (2019), we have laid the mathematical foundation of the framework, with an emphasis on the inclusion and the learning of the petrophysical signature (process in Fig. 1). In the current study, we focused on the coupling of several geophysical surveys with various physical properties, thus extending the process in Fig. 1. The third and last process, the geological identification, is still one where there is much room for advancement. The next step would be to add geological information, for example from a borehole. This could be used, for example, in our DO- study to adjust the depth of the HK unit. The inclusion and dynamic extrapolation of borehole lithology information through the inversion is part of our current research.

6 Conclusion

We have expanded our framework to include petrophysical and geological information as a coupling term to perform multi-physics joint inversions. We have described our strategies for handling multiple geophysical target misfits as well as a petrophysical and geological target misfit. We have presented our efforts to make the implementation modulable, extensible and shareable. Finally, we have demonstrated through the DO- synthetic example the gains that can be made by including various types of information into a single inversion. Only a joint approach for inverting the potential fields datasets allows us to delineate two kimberlite facies and to reproduce their petrophysical signature.

7 Acknowledgments

We sincerely thank the open-source software SimPEG community whose work has considerably facilitated the research presented here. Special thanks go to Dominique Fournier, for the implementation of the potential fields operators, Seogi Kang, for contributions to the Mapping module (Kang et al., 2015) and Rowan Cockett for his work in the creation and maintenance of SimPEG. We also thank Condor Consulting Inc., Peregrine Diamonds Ltd. and Kennecott for making the DO- example available for our research.


  • P. Afnimar, K. Koketsu, and K. Nakagawa (2002) Joint inversion of refraction and gravity data for the three-dimensional topography of a sediment–basement interface. Geophysical Journal International 151 (1), pp. 243–254. External Links: Document Cited by: §1.
  • U. Ascher and C. Greif (2011) A first course in numerical methods. edition, Society for Industrial and Applied Mathematics, Philadelphia, PA. External Links: Document, Link, Cited by: §3.4.
  • T. Astic, D. Fournier, and D. W. Oldenburg (2020) Joint inversion of potential fields data over the do-27 kimberlite using a gaussian mixture model prior. Submitted to Interpretation. Cited by: §3.4, §3.5, §4.5, §5, §5.
  • T. Astic and D. W. Oldenburg (2019) A framework for petrophysically and geologically guided geophysical inversion using a dynamic Gaussian mixture model prior. Geophysical Journal International 219 (3), pp. 1989–2012. External Links: ISSN 0956-540X, Document, Link, Cited by: Appendix A, Appendix A, §1, §2.4, §2.5, §2.5, §2.6, §2.6, §3.2, §3.2, §3.3, §3.5, §4.3, §4.5, §4.5, §5, §5, Algorithm 1.
  • T. Astic (2020) Collection of scripts for forward modeling and joint inversion of potential fields data. Note: 2020-01-31 External Links: Document, Link Cited by: §3.5.
  • N. Bournas, A. Prikhodko, K. Kwan, J. Legault, V. Polianichko, and S. Treshchev (2018) A new approach for kimberlite exploration using helicopter-borne tdem data. In SEG Technical Program Expanded Abstracts 2018, pp. 1853–1857. External Links: Document, Link, Cited by: §1.
  • R. Cockett, S. Kang, L. J. Heagy, A. Pidlisecky, and D. W. Oldenburg (2015) SimPEG: An open source framework for simulation and gradient based parameter estimation in geophysical applications. Computers and Geosciences 85, pp. 142–154. External Links: Document, ISSN 00983004, Link Cited by: §1, §2.3, §3.5.
  • A. P. Dempster, N. M. Laird, and D. B. Rubin (1977) Maximum likelihood from incomplete data via the {EM} algorithm. Journal of the Royal Statistical Society, series B 39 (1), pp. 1–38. Note: With discussion External Links: ISSN 0035-9246, Link Cited by: Appendix A, Appendix A.
  • S. G. R. Devriese, K. Davis, and D. W. Oldenburg (2017) Inversion of airborne geophysics over the do-27/do-18 kimberlites - part 1: potential fields. Interpretation 5 (3), pp. T299–T311. External Links: Document, Link, Cited by: §1.
  • J. Doetsch, N. Linde, I. Coscia, S. A. Greenhalgh, and A. G. Green (2010) Zonation for 3d aquifer characterization based on joint inversions of multimethod crosshole geophysical data. Geophysics 75 (6), pp. G53–G64. External Links: Document, Link, Cited by: §1.
  • L. A. Gallardo and M. A. Meju (2003) Characterization of heterogeneous near-surface materials by joint 2d inversion of dc resistivity and seismic data. Geophysical Research Letters 30 (13). External Links: Document Cited by: §1.
  • L. A. Gallardo and M. A. Meju (2004) Joint two-dimensional dc resistivity and seismic travel time inversion with cross-gradients constraints. Journal of Geophysical Research: Solid Earth 109 (B3), pp. . External Links: Document, Link, Cited by: §3.2.
  • L. A. Gallardo and M. A. Meju (2011) Structure-coupled multiphysics imaging in geophysical sciences. Reviews of Geophysics 49 (1). Note: doi: 10.1029/2010RG000330 External Links: Document, ISSN 8755-1209, Link Cited by: §1.
  • J. Giraud, E. Pakyuz-Charrier, M. Jessell, M. Lindsay, R. Martin, and V. Ogarko (2017) Uncertainty reduction in joint inversion using geologically conditioned petrophysical constraints. Geophysics 82 (6), pp. 1–61. External Links: Document, ISSN 0016-8033, Link Cited by: §1, §1.
  • D. Grana and E. Della Rossa (2010) Probabilistic petrophysical-properties estimation integrating statistical rock physics with seismic inversion. Geophysics 75 (3), pp. O21. External Links: Document, ISBN 00168033, ISSN 00168033 Cited by: §1.
  • D. Grana, T. Fjeldstad, and H. Omre (2017) Bayesian gaussian mixture linear inversion for geophysical inverse problems. Mathematical Geosciences 49 (4), pp. 493–515. External Links: ISSN 1874-8953, Document, Link Cited by: §1, §1.
  • E. Haber and D. Oldenburg (1997) Joint inversion: a structural approach. Inverse Problems 13 (1), pp. 63. External Links: Link, Document Cited by: §1.
  • L. J. Heagy, R. Cockett, S. Kang, G. K. Rosenkjaer, and D. W. Oldenburg (2017) A framework for simulation and inversion in electromagnetics. Computers & Geosciences 107, pp. 1–19. External Links: ISSN 0098-3004, Document, Link Cited by: §1, §3.5.
  • J. Jansen and K. Witherly (2004) The tli kwi cho kimberlite complex, northwest territories, canada: a geophysical case study. In 2004 SEG Annual Meeting, Cited by: §1, §4.1.
  • M. D. Jegen, R. W. Hobbs, P. Tarits, and A. Chave (2009) Joint inversion of marine magnetotelluric and gravity data incorporating seismic constraints: preliminary results of sub-basalt imaging off the faroe shelf. Earth and Planetary Science Letters 282 (1-4), pp. 47–55. External Links: Document Cited by: §1.
  • J. Kamm, I. A. Lundin, M. Bastani, M. Sadeghi, and L. B. Pedersen (2015) Joint inversion of gravity, magnetic, and petrophysical data — a case study from a gabbro intrusion in boden, sweden. Geophysics 80 (5), pp. B131–B152. External Links: Document, Link, Cited by: §1.
  • S. Kang, R. Cockett, L. J. Heagy, and D. W. Oldenburg (2015) Moving between dimensions in electromagnetic inversions. In SEG Technical Program Expanded Abstracts 2015, pp. 5000–5004. External Links: Document, Link, Cited by: §2.3, §7.
  • S. Kang, D. Fournier, and D. W. Oldenburg (2017) Inversion of airborne geophysics over the do-27/do-18 kimberlites — part 3: induced polarization. Interpretation 5 (3), pp. T327–T340. External Links: Document, Link, Cited by: §1.
  • P. Keating and P. Sailhac (2004) Use of the analytic signal to identify magnetic anomalies due to kimberlite pipes. Geophysics 69 (1), pp. 180–190. External Links: Document, Link, Cited by: §1.
  • P. G. Lelièvre, C. G. Farquharson, and C. A. Hurich (2012) Joint inversion of seismic traveltimes and gravity data on unstructured grids with application to mineral exploration. Geophysics 77 (1), pp. K1–K15. External Links: Document, Link, Cited by: §1, §3.2.
  • Y. G. Li and D. W. Oldenburg (1996) 3-D inversion of magnetic data. Geophysics 61 (2), pp. 394–408 (English). External Links: Document, ISBN 0016-8033, ISSN 1070485X Cited by: Appendix B, §4.1.
  • Y. Li, A. Melo, C. Martinez, and J. Sun (2019) Geology differentiation: a new frontier in quantitative geophysical interpretation in mineral exploration. The Leading Edge 38 (1), pp. 60–66. External Links: Document, Link, Cited by: §2.3.
  • Y. Li and D. W. Oldenburg (1998) 3-d inversion of gravity data. Geophysics 63 (1), pp. 109–119. External Links: Document, Link, Cited by: §4.1.
  • J. Macnae (1995) Applications of geophysics for the detection and exploration of kimberlites and lamproites. Journal of Geochemical Exploration 53 (1), pp. 213 – 243. Note: Diamond Exploration: Into the 21st Century External Links: ISSN 0375-6742, Document, Link Cited by: §1.
  • P. R. McGillivray and D. W. Oldenburg (1990) Methods for calculating fréchet derivatives and sensitivities for the non-linear inverse problem: a comparative study. Geophysical Prospecting 38 (5), pp. 499–524. External Links: Document, Link, Cited by: §3.4.
  • S. Mehanee, N. Golubev, and M. S. Zhdanov (2005) Weighted regularized inversion of magnetotelluric data. In SEG Technical Program Expanded Abstracts 1998, pp. 481–484. External Links: Document, Link, Cited by: §4.1.
  • A. T. Melo, J. Sun, and Y. Li (2017) Geophysical inversions applied to 3d geology characterization of an iron oxide copper-gold deposit in brazil. Geophysics 82 (5), pp. K1–K13. External Links: Document, Link, Cited by: §1.
  • D. W. Oldenburg, L. J. Heagy, S. Kang, and R. Cockett (2019) 3D electromagnetic modelling and inversion: a case for open source. Exploration Geophysics 0 (0), pp. 1–13. External Links: Document, Link, Cited by: §5.
  • D. W. Oldenburg, Y. Li, and R. G. Ellis (1997) Inversion of geophysical data over a copper gold porphyry deposit: a case history for mt. milligan. Geophysics 62 (5), pp. 1419–1431. External Links: Document, Link, Cited by: §1.
  • D. W. Oldenburg and Y. Li (2005) Inversion for applied geophysics: a tutorial. In D. K. Butler, ed., Near-Surface Geophysics, D. K. Butler (Ed.), pp. 89–150. External Links: Document, Link, Cited by: §2.2.
  • S. Onizawa, H. Mikada, H. Watanabe, and S. Sakashita (2002) A method for simultaneous velocity and density inversion and its application to exploration of subsurface structure beneath izu-oshima volcano, japan. Earth, Planets and Space 54 (8), pp. 803–817. External Links: ISSN 1880-5981, Document, Link Cited by: §1, §2.3, §5.
  • H. Paasche and J. Tronicke (2007)

    Cooperative inversion of 2d geophysical data sets: a zonal approach based on fuzzy c-means cluster analysis

    Geophysics 72 (3), pp. A35–A39. External Links: Document, Link, Cited by: §1.
  • R. L. Parker (1977) Understanding inverse theory. Annual Review of Earth and Planetary Sciences 5 (1), pp. 35–64. External Links: Document, Link, Cited by: §2.2, §2.6.
  • O. Portniaguine and M. S. Zhdanov (2002) 3‐d magnetic inversion with data compression and image focusing. Geophysics 67 (5), pp. 1532–1541. External Links: Document, Link, Cited by: §4.1.
  • P. Shamsipour, D. Marcotte, and M. Chouteau (2012) 3D stochastic joint inversion of gravity and magnetic data. Journal of Applied Geophysics 79, pp. 27 – 37. External Links: ISSN 0926-9851, Document, Link Cited by: §1.
  • A. Sosa, A. A. Velasco, L. Velazquez, M. Argaez, and R. Romero (2013) Constrained optimization framework for joint inversion of geophysical data sets. Geophysical Journal International 195 (3), pp. 1745–1762. External Links: ISSN 0956-540X, Document, Link, Cited by: §3.2.
  • M. D. Stefano, F. G. Andreasi, S. Re, M. Virgilio, and F. F. Snyder (2011) Multiple-domain, simultaneous joint inversion of geophysical data with application to subsalt imaging. Geophysics 76 (3), pp. R69–R80. External Links: Document, Link, Cited by: §1.
  • J. Sun and Y. Li (2015) Multidomain petrophysically constrained inversion and geology differentiation using guided fuzzy c-means clustering. Geophysics 80 (4), pp. ID1. External Links: Document, Link, /gsw/content_public/journal/geophysics/80/4/10.1190_geo2014-0049.1/2/geo2014-0049.pdf Cited by: §1.
  • J. Sun and Y. Li (2016) Joint inversion of multiple geophysical data using guided fuzzy c-means clustering. Geophysics 81 (3), pp. ID37–ID57. External Links: Document, Link, Cited by: §1, §3.2.
  • J. Sun and Y. Li (2017) Joint inversion of multiple geophysical and petrophysical data using generalized fuzzy clustering algorithms. Geophysical Journal International 208 (2), pp. 1201–1216. External Links: Document, ISSN 1365246X Cited by: §1, §1.
  • J. Sun and Y. Li (2018) Magnetization clustering inversion — part 1: building an automated numerical optimization algorithm. Geophysics 83 (5), pp. J61–J73. External Links: Document, Link, Cited by: §3.2.
  • A. N. Tikhonov and V. Y. Arsenin (1977) Solutions of ill-posed problems. V. H. Winston & Sons, Washington, D.C.: John Wiley & Sons, New York. Note: Translated from the Russian, Preface by translation editor Fritz John, Scripta Series in Mathematics External Links: MathReview (M. Z. Nashed) Cited by: §2.2.

Appendix A Updating the Gaussian mixture model in multi-dimensions

In this section, we generalize the learning of the GMM parameters presented in Astic and Oldenburg (2019)

to multidimensional Gaussian distributions, representing multiple physical properties at once. The main difference comes from the fact that the means are now vectors (scalar in 1D) and that the scalar variances in 1D become covariance matrices. We define the problem as a Maximum A Posteriori estimate of the GMM means, proportions and covariance matrices:


We choose to follow a conjugate prior approach for the choice of the prior distributions

. We use the Maximum A Posteriori Expectation-Maximization (MAP-EM) algorithm to estimate the parameters

(Dempster et al., 1977).

The computation of the responsibilities for the E-step of the MAP-EM algorithm (Dempster et al., 1977) stay the same as in Astic and Oldenburg (2019), except for the dimension of the parameters:


The conjugate prior for the proportions follow a Dirichlet distribution:

Thus giving the update:

where is the volume of the cell and is the volume of the active mesh. This allows the estimates to be mesh-independent by using volumetric proportions instead of cell counts.

The semi-conjugate for the means and covariance matrices follow respectively a Gaussian and Inverse-Wishart distributions:


The update to the means become:


The update to the covariance matrices takes the form of:


Appendix B Working with nonlinear petrophysical relationships

Figure 8: Gaussian mixture with various polynomial relationships: one linear (no addition required), one quadratic and one cubic. The 1D probability distribution for each physical property, and the respective histogram of each unit, are projected on the left and bottom panels.

While linear trends can be accounted through the covariance matrix to obtain elongated clusters, it is also possible to account for nonlinear relationships between physical properties in our framework (Fig. 8). This is achieved by composing the Gaussian function with the nonlinear relationship between the physical properties of the particular rock unit . This corresponds to using for each rock unit in the GMM in equation 14.

In this section, we present a joint inversion of two linear problems whose respective physical properties have nonlinear petrophysical relationships between each other. For simplicity, the physics of the two problems are the same and is based on the examples previously used in Li and Oldenburg (1996). The models are discretized on a 1D mesh defined on the interval and divided into 100 cells. Both datasets consist of data points at various frequencies equally distributed from to evaluated according to:


Each model presents three distinct units. The background for both models is set to while the two other units are linked through a quadratic and a cubic relationship respectively.

We inverted the two datasets using the presented joint PGI framework, including a priori knowledge of the petrophysical distributions and relationships. The result can be seen in the first row of Figs 9(a) to (c). Note that the petrophysical relationships are well reproduced. This allows recovering in great detail the two models.

For comparison, we also jointly inverted the datasets but without the polynomial relationships (by merely fitting a Gaussian distribution to each unit). The result can be seen in the second row of Figs 9(d) to (f). While the overall structures are well recovered, we miss some details of the models. The background is not as flat as with the full information, the lower tip of the model of problem is completely missed.

We also inverted both datasets individually using the classic Tikhonov inversion. The result is shown in the last row of Figs 9(g) to (i). Results are smoother, as expected. The background presents even more variations, but the overall structures are recovered. Same as for the joint inversion without the polynomial relationships, the details are missing.

Figure 9: Linear joint inversions with various types of physical properties relationships (no trend, quadratic, cubic). Panels a) to c) show the inversion result with PGI using the known nonlinear relationships. The first panel a) shows the result for the first problem, b) for the second problem. The panel c) shows the cross-plot of the models over the contour of the GMM with nonlinear relationships. Panels d) to f) show the result with PGI without nonlinear relationships. In panel f), we show the used GMM without nonlinear relationships. Panels g) to i) show the Tikhonov inversions result.

Appendix C Additional models obtained by individual PGI for the TKC synthetic case study

Figure 10: (a) magnetic PGI using PK/VK magnetic signature. Note that the recovered volume is much larger than the volume of PK/VK recovered from the gravity inversion with this unit density contrast (Fig. 5c); (b) gravity PGI using the HK density contrast signature. Note again the larger volume compared to Fig. 5(d).

Appendix D Algorithm

Initialization :

  • Input:

    • Initial geophysical model , GMM parameters and geological model .

  • Parameters:

    • Objective function: data’s noise , and .parameters.

    • Localized prior: specific for available locations , weights .

    • GMM prior weights: for the means, variances and proportions.

    • Optimization: -cooling factor (), sufficient decrease rate ().

  • Output:

    • , , .

while  and  do

6       543 Objective Function Descent Step:
  • Compute a model perturbation using an inexact Gauss-Newton.

  • Line search with Wolfe condition to find an that satisfy a sufficient decrease of

  • Return .

Update Petrophysical Distribution
  • Fit a new GMM on such as in equations 35, 40 and 42 until no sufficient increase of the posterior is observed.

  • Compute the membership of the current model as in equation 16 using .

  • Update and according to 17 and 18 respectively using .

Update regularization weights:
if  then
7               Decrease :
8        else if  and  then
9               Increase :
10        if (optional) and and  then
11               Include in Smoothness
13       12 Update geophysical data misfit weights:
if  then
14               update according to equation 31
16 end while
Algorithm 1 PGI extended from Astic and Oldenburg (2019) for joint inversion