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

1Introduction¶

Mineral deposits, or other geologic features, are characterized by different physical properties, and hence multiple geophysical surveys can be used to delineate them. For example, kimberlites have signatures that depend upon density, magnetic susceptibility, and electrical conductivity. They are often discovered through data collected in airborne surveys and appear as circular low gravity anomalies, with high magnetic responses, and sometimes negative electromagnetic transient responses Macnae, 1995Keating & Sailhac, 2004Bournas et al., 2018. Although a joint interpretation of several individually inverted datasets can significantly improve our understanding of the subsurface Oldenburg et al., 1997Devriese et al., 2017Kang et al., 2017Giuseppe et al., 2014Paasche et al., 2006Paasche, 2016Martinez & Li, 2015Melo et al., 2017, multiple cases studies have shown that multi-physics inversions can reveal information that was not accessible through individual geophysical dataset inversions Doetsch et al., 2010Jegen et al., 2009Kamm et al., 2015Lelièvre & Farquharson, 2016. An extensive compilation of integrated imaging methods and their applications can be found in Moorkamp et al. (2016). Multi-physics inversions require a coupling term that mathematically describes a relationship between the different physical property models responsible for the geophysical data. Coupling methods generally use one or a combination of structural or physical property relationships.

The first frameworks for joint inversion focused on linking geophysical models through their structural similarities. Haber & 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 & Meju (2003) with the introduction of the concept of cross-gradient between geophysical models. This approach has become commonly used, and both Gallardo & Meju (2011) and Meju & Gallardo (2016) provide in-depth reviews of the method and its application. However, this strategy has several limitations: 1) Meju & Gallardo (2016) points out that “not all physical property distributions in the subsurface will be structurally coincident”; 2) it is unable to reproduce documented or expected petrophysical information Sun & Li, 2017. These drawbacks can be overcome by using other coupling methods.

The second coupling approach uses physical property relationships to link geophysical models. Some of the earliest works used empirical constitutive formulae as their physical properties constraint Afnimar et al., 2002Hoversten et al., 2006Chen et al., 2007. De Stefano et al. (2011) combined this approach with the above mentioned cross-gradient method for sub-salt imaging. Moorkamp et al. (2011) compared the constitutive relationship and the cross-gradient approaches on a 3-D synthetic example combining magnetotelluric, gravity, and seismic data. They concluded that, overall, a cross-gradient approach was preferable compared to using constitutive equations because deviations from the constitutive relations resulted in artefacts in the inverted models; in those situations, the cross-gradient method gave consistent satisfactory results. They also pointed out that the cross-gradient method relies on fewer assumptions about the models than the constitutive equations. Some stochastic frameworks have also been proposed that leverage geostatistical tools to define relationships between physical properties. Chen & Hoversten (2012) used a similar coupling approach as in Bosch (2004) Moorkamp et al., 2016 by building a rock-physics model from borehole data to jointly invert seismic and controlled-source electromagnetic data. Shamsipour et al. (2012) used the geostatistical techniques of cokriging and conditional simulation to jointly invert gravity and magnetic data assuming that the auto- and cross-covariances of the density and magnetic susceptibility follow a linear model of coregionalization. On the deterministic side, of which the framework we present belongs, recent frameworks use clustering techniques such as the fuzzy C-means (FCM) algorithm, which was first used in Paasche & Tronicke (2007) and further expanded in Lelièvre et al. (2012). This approach adds a clustering term to the objective function, which allows more flexible relationships between physical properties. Beyond the addition of the FCM term to the objective function, Sun & Li (2015) introduced an iterative update to the cluster centres throughout an inversion for a single physical property, a technique they called guided FCM; this introduced the concept of uncertainties for physical properties into the inversion. In Sun & Li (2016), they generalized this work to consider multiple physical properties, and further in Sun & Li (2017), they added tools to their approach to consider various types of correlations between physical properties (linear, quadratic, etc.). Giraud et al. (2017) represented the petrophysical information as a fixed GMM and focused on reducing uncertainties in stochastic geological modelling by linking potential field inversions and geological models through petrophysical information. To this end, they added, to the Tikhonov objective function, a sum of least-squares differences between the GMM function, evaluated at the current model, and reference values representing the likelihood of their prior knowledge. In Giraud et al. (2019), they modified their formulation of this coupling term to work with a least-squares difference between the log-likelihood of the GMM and their reference values. In both formulations, they required extensive and fixed quantitative petrophysical and geological information. In practice, that information is not often available and might be only qualitative.

In Astic & Oldenburg (2019), we presented a petrophysically and geologically guided inversion (PGI) framework that generalized concepts presented in previously published researches. We showed how geological and petrophysical knowledge represented as a univariate GMM could be incorporated in a voxel-based geophysical inversion through a single smallness term. Each contrasting geological unit was represented by a univariate Gaussian distribution, which summarized its physical property signature. Geological information was included in the GMM through its proportions in a manner similar to that of Giraud et al. (2017). The log-likelihood of the GMM was then used to regularize the geophysical inversion; this is analogous to the approach taken by Grana & Della Rossa (2010) and Grana et al. (2017). Incorporating both petrophysical and geological information into a single smallness term in the regularization had several advantages. First, it did not require adding a term in the classic Tikhonov formulation. Second, this approach brought the petrophysical data to the same level as the geophysical data, which allowed us to define a misfit, with a target value, between the geophysical model and the petrophysical and geological data. The iteration steps were decomposed into a suite of cyclic optimization problems across the geophysical, petrophysical, and geological data. The petrophysical step formalized the idea of learning the physical property mean values described in Sun & Li (2015) and generalized it to enable the variances and proportions of the GMM to be learned as well during the inversion. These updates to the GMM parameters allowed us to work with partial petrophysical information. The geological step built at each iteration a "quasi-geology model" Li et al., 2019, based on the current geophysical and petrophysical models. We applied the PGI approach on synthetic and field data, but the analysis was restricted to single datasets and a single physical property.

This study extends the PGI framework to perform multi-physics inversions, involving several physical properties. We show how geological and petrophysical information represented as a multivariate GMM can be used to couple multiple voxel-based geophysical inversions through a single smallness term in the regularization. The updates to the means, covariances and proportions of the GMM are extended to inversions with multiple physical properties, which expands the work of Sun & Li (2016). Tools for handling various types of relationships between physical properties are designed; this further develops ideas presented in Sun & Li (2017). These capabilities, and the gains that they generate, are demonstrated on inversions of synthetic gravity and magnetic data.

In this paper, we start by reviewing the key concepts of the PGI approach and generalize them to the case of multiple physical properties and governing equations for performing multi-physics inversions. We then delineate our strategy to address the numerical challenges of finding a solution to the inverse problem that fits each dataset and, at the same time, adequately fits the petrophysical data. Finally, we demonstrate the advantages of performing joint inversions, with various levels of prior knowledge, by using a synthetic model of the DO-27 Tli Kwi Cho kimberlite pipe, Northwest Territories, Canada Jansen & Witherly, 2004.

2Background and motivation for a multi-physics PGI framework¶

In this section, we present the conventions we use for notation, provide background on the univariate PGI framework, and motivate its extension to multiple physical properties.

2.1Notation 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 $\mathbf{m}$.

• Bold uppercase symbols designate matrices, such as the weight matrix $\mathbf{W}$.

For a multi-physics inversion, the geophysical model is likely to contain multiple physical properties. Several surveys might be associated with the same physical property (e.g. gravity and gravity gradiometry) or one survey might depend on several physical properties (e.g. electromagnetic surveys depend on both electrical conductivity and magnetic permeability). We thus adopt the following notations for the geophysical model $\mathbf{m}$ indices:

A row of the $\mathbf{M}$ matrix represents all the $q$ physical properties that live in the same location. A column represents a single physical property at all the $n$ cells of the mesh. For clarity, we are consistent throughout this study with the index notation. The index $i$ always refers to the cell number, from 1 to $n$. The vector $\mathbf{m}_i$ then denotes all of the physical properties at the $i$^th cell:

Likewise, we denote the vector model for a single physical property on the whole mesh with the index $p \in \{1..q\}$ with a superscript:

Lastly, the geophysical model at iteration $t$ of an inversion is denoted with parentheses $\mathbf{m}^{(t)}$.

2.2The Tikhonov inverse problem and its PGI augmentation¶

The geophysical inverse problem can be posed as an optimization process where the goal is to find a geophysical model $\mathbf{m}$ that minimizes an objective function Φ. Using the same formulation of the inverse problem as in Oldenburg & Li (2005), the geophysical optimization problem takes the form:

In equation (4), the vector $\mathbf{m}$ is the geophysical model, which represents physical properties on a mesh. The term $\Phi_d$ is the geophysical data misfit. The regularization is composed of a smallness term $\Phi_s$, that penalizes model-values different from a reference model, and smoothness terms $\Phi_v$, which penalize variations between adjacent cells; those terms are weighted by positive scaling parameters $\left\{\alpha\right\}$. The trade-off parameter β is a positive scalar that adjusts the relative weighting between the regularization and the data misfit. A value of β is sought so that the data misfit $\Phi_d$ is below an acceptable target misfit $\Phi_d^*$ Parker, 1977. The scaling parameters $\alpha_s$ and $\{\alpha_{v}\}$ weight the relative importance of the smallness and smoothness terms. In the Tikhonov inversion, first introduced in Tikhonov & Arsenin (1977), each term of the objective function takes a least-squares form. In particular, the smallness term, which is essential in the PGI framework, can be written as:

where $\mathbf{m}_{\text{ref}}$ is a reference model and the matrix $\mathbf{W}_s$ represent local weights.

The PGI approach can be considered as an augmentation of the well-established Tikhonov inversion. This is detailed in Astic & Oldenburg (2019), where we start from a probabilistic formulation of the least-squares inverse problem Tarantola, 2005, to then include petrophysical and geological information in the form of a GMM. The resulting term is analogous to the smallness, but the reference model $\mathbf{m}_{\text{ref}}$ and the weights $\mathbf{W}_{s}$ are updated at each iteration. The minimization of the geophysical objective function is labelled Process 1 in the PGI framework (Fig. 1).

2.3Multivariate GMM: modelling multiple physical properties¶

To extend the approach presented in Astic & Oldenburg (2019), we represent the petrophysical signature of each geological unit $j$ ($j=1..c$) as a multivariate Gaussian probability distribution, denoted by $\mathcal{N}$. The Gaussian function representing the probability distribution of the $q$ physical properties of interest for each unit is defined by its mean $\mathbf{\mu}_j$ (vector of size $q$), and its covariance matrix $\mathbf{\Sigma}_j$ (matrix of size $q\times q$), plus its proportion $\pi_j$.

The multivariate Gaussian Mixture Model (GMM) simply sums the Gaussian probability distribution representing each known unit, weighted by their proportion:

where the variable Θ holds the GMM global variables $\Theta= \left\{\pi_j, \mathbf{\mu}_j, \mathbf{\Sigma}_j\right\}_{j=1..c}$. With some modifications, the GMM can also represent nonlinear relationships, such as polynomial as presented in Onizawa et al. (2002). We present those modifications in Appendix 1, along with an example of an inversion with various nonlinear relationships between two physical properties.

The geological classification (or membership) is denoted $\mathbf{z}$. It is defined as the most probable geological unit, given a set of values $\mathbf{x}$ for the $q$ physical properties:

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

2.4Motivation for simultaneously inverting multiple physical properties¶

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 (7). The bottom and left panels represent the marginal GMM probability distribution for each physical property individually. We notice that, while all three units are distinct in the two-dimensional space, they can overlap significantly when only considering one physical property at the time. For physical property 1, rock units 1 and 2 are distinct while rock unit 3 is indistinguishable from rock unit 1. For physical property 2, rock units 1 and 3 are now distinct while unit 2 is indistinguishable from either rock units 1 or 3. This highlights that it is only by jointly inverting for several physical properties that we might be able to uniquely identify three rock units. Units that might not be distinguishable in one survey may be in another one, and by simultaneously working with both physical properties in an inversion, we are able to explore the 2 (or multi)-dimensional physical property space in the centre panel of Fig. 2. 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.

3Extension of the PGI framework to multi-physics inversions¶

3.1Definition of the GMM prior¶

Astic & Oldenburg (2019) developed a GMM smallness prior to include petrophysical and geological information in inversions involving only one type of geophysical survey with only one physical property to recover. In equation (8), we propose a generalized version of the GMM smallness prior that is designed to couple multiple physical properties and incorporate geological information. Note that now the means are vectors that are 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 $i$) and lithologically (index $j$) dependent.

where:

• $c$ is the number of distinct rock units.

• $n$ is the number of active cells in the mesh.

• $\mathbf{m}_i$ represents the physical property values at location $i$.

• $\mathcal{P}(z_{i}=j)$ is the a priori probability of observing rock unit $j$ at location $i$. It can be either constant over the whole area, then denoted by $\pi_j$, or locally determined by a priori geological knowledge.

• $\mathbf{\mu}_j$ contains the means of the physical properties of rock unit $j$.

• $\mathbf{\Sigma}_j$ is the covariance matrix of the physical properties of rock unit $j$.

• $\mathbf{W}_i^{-1}$ is a weighting term at location $i$, used for example to include depth or sensitivity weighting. We define it from weights $\left\{{w}_{i, p}, ~i=1..n, ~p=1..q\right\}$; this is a scalar value for each cell $i$ and physical property $p$. $\mathbf{W}_i$ is defined as a diagonal matrix made of the combination of all the weights at the specific location $i$: $\mathbf{W}_i = \text{diag}(\mathbf{w}_{i})$, with the same notation convention as for the model $\mathbf{m}$. This allows the weighting to be different for each physical property. We can thus weight each physical property according to the survey on which it depends.

• Θ holds the GMM global variables $\Theta= \left\{\pi_j, \mathbf{\mu}_j, \mathbf{\Sigma}_j\right\}_{j=1..c}$.

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

3.2The multi-physics PGI geophysical objective function¶

In Astic & Oldenburg (2019), we demonstrated how to use the negative log-likelihood of a univariate GMM as the smallness term in the Tikhonov inverse problem. The resulting smallness term could be approximated by a least-squares misfit between the current model and a reference model $\mathbf{m}_{\text{ref}}$, which was updated at each iteration, as was the smallness matrix $\mathbf{W}_s$. Those dynamic reference model and smallness matrix updates were determined based on the current geophysical model and the geological and petrophysical prior information. The goal of the least-squares approximation was to enable the use of the PGI framework with compiled codes working with the Tikhonov formulation.

Generalizing the result obtained in Astic & Oldenburg (2019) to multiple physical properties, we use the negative log-likelihood of the GMM defined in equation (8) to obtain a single smallness term that couples all of the model parameters. That smallness term can be approximated by the following least-squares misfit:

with:

where $\mathbf{\Sigma}^{-1/2}$ is the upper triangular matrix from the Cholesky decomposition of the precision matrix $\mathbf{\Sigma}^{-1}$.

Note that our implementation can handle either the log-likelihood of the GMM or its least-squares approximation. When the petrophysical signature of the rock units are not known, it is possible to learn the parameters of the GMM Astic & Oldenburg, 2019 (Fig. 1, Process 2). The extension of that learning process to multiple physical properties can be found in Appendix 2.

3.3Petrophysical target misfit¶

The PGI smallness expresses a misfit between the petrophysical and geological information and the geophysical model. To measure the goodness of fit and define a stopping criterion for the petrophysical misfit, Astic & Oldenburg (2019) defined a measure $\Phi_{\text{petro}}$ and its target value $\Phi_{\text{petro}}^*$. This measure is similar to the PGI smallness term but without the weights $\mathbf{W}_i$ in equation (12). The same approach can be taken here to define the value $\Phi_{\text{petro}}$ for the multivariate case:

with:

where $z_i$ and $\mathbf{m}_{\text{ref}}$ are the same as in equations (10) and (11).

Looking at the term $\Phi_{\text{petro}}$ (equation (13)) from a probabilistic point of view Tarantola, 2005Astic & Oldenburg, 2019, each variable $\mathbf{m}_i$ follows a multivariate Gaussian of dimension $q$, with mean $\mathbf{m}_{\text{ref}}(\Theta, z_i)$, and covariance matrix $\left(\mathbf{\tilde{W}}_{s}(\Theta, z_i)^\top \mathbf{\tilde{W}}_{s}(\Theta, z_i)\right)^{-1}$. Thus, the variable $\mathbf{\tilde{W}}_{s_i}(\Theta, z_i)(\mathbf{m}_i-\mathbf{m}_{\text{ref}_i})$ follows a multivariate Gaussian variable with mean $\mathbf{0}$ and an identity covariance matrix. Thus, the sum in $\Phi_{\text{petro}}$ follows a chi-squared distribution and we can apply Pearson’s chi-squared test Pearson, 1900. The target misfit value $\Phi_{\text{petro}}^*$ is defined as the expectation of $\Phi_{\text{petro}}$:

with $n$ being the number of active cells in the mesh and $q$ being the number of physical properties. This generalizes the result obtained in Astic & Oldenburg (2019) ($q=1$). This is a similar approach to the definition of a target misfit for the geophysical data as given in Parker (1977).

Our algorithm stops when all of the target misfits, geophysical and petrophysical, are achieved. In the next section, we present our strategy for handling multiple geophysical data misfits as well as an additional petrophysical misfit, each with its target value we seek to reach.

4Numerical considerations for reaching multiple target misfits¶

We have multiple geophysical data misfits that we wish to fit. The inclusion of petrophysical and geological data with PGI adds another data misfit term that also needs to reach its target misfit (section 3.3). In this section, we provide our strategies for choosing and dynamically adjusting 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 3.

4.1Objective function with multiple geophysical data misfits¶

The objective function we seek to minimize for the multi-physics inversion process (Fig. 1, Process 1) takes the form:

with:

The data misfit term $\Phi_d(\mathbf{m})$ now contains multiple geophysical data misfits, each defined as a weighted least-squares norm. $\Phi_d^k$ is the data misfit of the $k$^th survey, where the forward operator $\mathbb{F}^k$ generates the predicted data for that survey from $\mathbf{m}^{\{k\}}$. The notation $\mathbf{m}^{\{k\}}$ denotes the subset of model parameters associated with the $k$^th 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 $k$^th survey is symbolized by $\mathbf{d}_{\text{obs}}^k$ and the uncertainty on those measurements by the matrix $\mathbf{W}_{d}^k$. Each data misfit $\Phi_d^k$ is weighted by a scaling parameter $\chi_k$. Those $\{\chi\}$ scaling parameters are important for balancing the various geophysical data misfits and finding a solution that fits all of them. Our approach for updating these parameters is developed later in this section.

The regularization is still composed of the smallness and smoothness terms. The smallness term $\Phi_s$ is our coupling term, which is defined in equation (9). The smoothness terms, one for each direction and physical property, are represented by $\Phi_{v, p}$. In the smoothness terms (equation (18)), the smoothness operators (usually first or second-order difference) are represented by the matrix $\mathbf{L}_v$, and weights (sensitivity or depth) are represented by the matrix $\mathbf{W}_{v,p}$.

Equation 16 is an intricate objective function that is the sum of many quadratic regularization terms, each of which is multiplied by an adjustable constant. Finding values for these constants and carrying out a nonlinear inversion to produce a model that acceptably fits the data, and is a good candidate for representing information from the true geology model, is numerically challenging.

In the following sections, we present our approach for estimating and updating these parameters throughout the inversion. The smoothness scaling parameters $\{\alpha_{v, p}\}$ are the only values that we keep constant. The scaling parameter $\alpha_s$ weights the importance of the petrophysical misfit term, and the trade-off parameter β influences the importance of the regularization (which contains the petrophysical misfit) relative to the geophysical data misfits. The geophysical data misfit scaling parameters $\{\chi\}$ are used to adjust the relative importance of each geophysical dataset. Values of β, $\alpha_s$, and $\{\chi\}$ are sought so that each geophysical data misfit $\Phi_d^k$ is below or equal to its target misfit ${\Phi_d^k}^*$, along with a value of the petrophysical data misfit $\Phi_{\text{petro}}$ that is less than or equal to its target misfit $\Phi_{\text{petro}}^*$.

4.2The regularization scaling parameters β and $\{\alpha\}$¶

Two types of scaling parameters act on the regularization terms; they are the trade-off parameter β and the $\{\alpha\}$ parameters. In our implementation, we keep the $\{\alpha_{v,p}\}$ parameters acting on the smoothness terms constant while we update β and $\alpha_s$ to reach a suitable solution to the PGI problem. Next, we outline our strategies for each of these parameters.

4.2.1Fixed parameters: the smoothness scaling parameters $\{\alpha_{v, p}\}$¶

In the smoothness terms, the scaling parameters $\{\alpha_{v, p}\}$ control the relative importance of spatial derivative terms in the regularization. Each set of assigned values will yield different outcomes. This is often a way in which model space can be explored (e.g. preferential smoothness in some directions, see Williams (2008)Lelièvre et al. (2009)). They are generally specified a priori, and we keep them fixed in our objective function throughout the inversion process. In addition to the common practical considerations for choosing the smoothness parameters Oldenburg & Li, 2005Williams, 2006, we use the $\{\alpha_{v, p}\}$ to weight each physical property, by dividing it by the square of its expected maximum amplitude (available through the GMM means if provided). This helps equalize the contribution of the smoothness terms to the objective function value when parameters have widely different scales (like density, log- electrical conductivity and magnetic susceptibility contrasts). The scaling of the physical properties in our extended smallness term (equation (9)) is taken care of by the covariance matrices of the GMM.

4.2.2Adjusted parameters: β and $\alpha_s$¶

In the case of a single geophysical data misfit with a petrophysical misfit, Astic & Oldenburg (2019) developed a strategy for cooling β and warming $\alpha_s$ to find a solution to the inverse problem that reaches the target values of both misfits. This approach is still appropriate in the multi-physics inversion framework and is what we use in this study (step 6 in algorithm 1). For the multi-physics case, we alter it in the following way: when all geophysical data misfits are equal or below their target value, the strategy for warming the scaling parameter $\alpha_s$ on the coupling term is:

4.3Balancing the geophysical data misfits with the scaling parameters $\{\chi\}$¶

Our goal is to develop a strategy for scaling multiple geophysical data misfits so that each geophysical dataset is adequately fit. We propose a strategy where the scaling parameters $\left\{\chi_k\right\}_{(k=1..r)}$ in equation (17) are successively updated. Our approach has a heuristic foundation and does not incur the significant computational cost often associated with optimization-based approaches, and generalize to any number of geophysical data misfits. Before presenting its details, we first outline some strategies that others have taken in addressing this problem.

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

Several approaches for weighting multiple geophysical data misfits have been proposed in the recent literature. Some frameworks do not follow any prescribed strategy for updating the scaling parameters of the geophysical data misfits. This is the case for the approaches proposed by Sun & Li (2016)Sun & Li (2017) and Sosa et al. (2013); both keep those scaling parameters constant. Sun & Li (2017) use values of unity, while Sosa et al. (2013) normalize each geophysical data misfit by its number of data. In our experience, keeping the weights constant has led us to overfit some surveys while underfitting others. Other frameworks have adopted the approach of running their joint inversions for multiple combinations of parameters. For three geophysical data misfits, Moorkamp et al. (2011) adopted a manual check-and-guess approach to adjust the parameters. For two geophysical data misfits, Giraud et al. (2019) ran a subset of their inversion hundreds of times for various combinations of scaling parameters before choosing values based on the L-curve principle Hansen & O'Leary, 1993Hansen, 2000Santos & Bassrei, 2007. They then manually "fine-tuned" those values using the full joint inversion problem. To avoid the issue of having to choose multiple appropriate scaling parameters, Bijani et al. (2017) developed a compromise between deterministic and stochastic optimizations for joint inversions. They adopted a "Pareto Multi-Objective Global Optimization" strategy with genetic algorithms that generate populations of candidate models that "simultaneously minimize multiple objectives in a Pareto-optimal sense," rather than working with a fully aggregated objective function. This approach was still computationally expensive and limited to small 2-D studies in the paper.

To limit the number of multiple runs of the same inversion, Lelièvre et al. (2012) devised a rigorously defined, but computationally expensive, strategy for dynamically balancing two geophysical data misfits. Their approach relied on adjusting first the trade-off parameter until the two data misfits are Pareto-optimal. Next, it adjusted the relative weights of the two surveys to fit both geophysical surveys. It then reinforced the importance of the coupling term before going into another round of adjustments of the trade-off and surveys weights parameters. The approach developed in Astic & Oldenburg (2019) for a single geophysical data misfit, but with a petrophysical data misfit, is related to the work of Lelièvre et al. (2012). Both focused first on fitting the geophysical data misfit terms and then adjusting the coupling term. On the contrary, the strategy presented in Gallardo & Meju (2004) favoured the cross-gradients coupling over the geophysical misfits.

Here we define a practical, computationally inexpensive, heuristic 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 thereby generalize the work of Lelièvre et al. (2012). For the full algorithm, the reader can refer to Appendix 3.

4.3.2Strategy for updating $\{\chi\}$¶

We now define our strategy for weighting the multiple geophysical data misfits to reach all the target values. We use the scaling parameters $\left\{\chi_k\right\}_{(k=1..r)}$ defined in equation (17). We dynamically update each geophysical misfit scaling parameter based on its current misfit and target value, compared to the other surveys.

We start with a set of initial scaling parameters $\{\chi\}$ that sums to unity. To ensure the progress of all data misfits, while limiting the possibilities of overfitting any given term, we update the scaling parameters $\{\chi\}$ during the inversion. Our approach is philosophically similar to what is proposed in Astic & Oldenburg (2019) for β and $\alpha_s$ in order to balance the geophysical data misfit and the petrophysical misfit at each iteration in the inversion, and generalizes ideas proposed in Lelièvre et al. (2012) to more than two geophysical data misfits. 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 scaling parameters of the other geophysical data misfit terms. We then normalize the sum of the scaling parameters to be equal to unity again; this is to keep the importance of the total $\Phi_d$ term relatively similar before and after adjusting the scaling parameters $\{\chi\}$. If several surveys are below their respective targets, we simply use the median of the ratios to warm the scaling parameters of the still unfit surveys. Thus, at any iteration $(t)$ of the geophysical inverse problem, if an ensemble of $\left\{k_f\right\}$ surveys has reached their respective targets, we warm the scaling parameters of the remaining $\left\{k_u\right\}$ surveys that are not yet fit as (step 7 in algorithm 1):

then we normalize the sum:

An example of convergence curves for the data misfits and evolution of the dynamic scaling parameters is proposed in Fig. 8 for a multi-physics PGI with full petrophysical information.

Our strategy has proven to be insensitive to the initialization of the scaling parameters $\{\chi\}$ for linear problems. To demonstrate this point, we show in Fig. 3 the evolution of the scaling parameters $\{\chi\}$ for three multi-physics PGI with full petrophysical information. The synthetic example presented in section 6.4 is run with various initializations $\{\chi_0\}$. The outcomes of all these three PGI were similar to the result we show in Fig. 7. The scaling parameters χ associated with the magnetic and gravity misfits, respectively, all finish at approximately the same value even though the initializations are very different. The final χ scaling values are about 0.8 for the gravity data misfit and around 0.2 for the magnetic data misfit. This is an appealing property as it reduces the need to fine-tune the initialization of the scaling parameters $\{\chi\}$, which can be costly for large-scale inversions.

5Numerical implementation¶

We implemented our framework as part of the open-source software SimPEG Cockett et al., 2015Heagy et al., 2017. As such, we are able to share both the software environment and the scripts to reproduce the examples shown in this paper Astic, 2020. In this section, we highlight some key points of our implementation to encourage the use of this work and future collaborations. A more detailed tutorial is provided in Appendix 4, which lays out a pseudocode sketch of the implementation.

SimPEG is designed to be a modular, extensible framework for simulations and inversions of geophysical data. In particular, two features enabled us to focus our implementation efforts on the PGI framework, while using tools provided by the open-source community (such as the forward operators):

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

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

The first point enables 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 $\mathbf{W}_d$, and methods, including functions to evaluate the misfit given a model as well as derivatives for use in the optimization routines. To construct a gravity and magnetic joint inversion, we first define each misfit term independently and then 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.

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 update to components of the inversion, including the data misfit scaling parameters (equation (20)), smallness weights and reference model (equations (11) and (12)), and for evaluating the target misfits and stopping criteria for the inversion.

6Example: The DO-27 kimberlite pipe¶

In this section, we illustrate the joint PGI approach on synthetic gravity and magnetic data based on the DO-27 kimberlite pipe Jansen & Witherly, 2004, which is composed of two different kimberlite facies. We compare standard Tikhonov inversions of the individual geophysical datasets and independent PGIs of the gravity and magnetic data with the multi-physics PGI approach. Both the Tikhonov inversions and single-physics PGIs produce models that enable only a binary distinction: kimberlite or host rock. Only the multi-physics PGI allows us to identify the two kimberlite facies as distinct from the background host rock.

Scripts and Jupyter notebooks to reproduce the examples presented in this study are available through GitHub at https://github.com/simpeg-research/Astic-2020-JointInversion Astic, 2020.

6.1Setup¶

The DO-27 Kimberlite pipe (Northwest Territories, Canada) is part of a complex known as the Tli Kwi Cho (TKC) kimberlite cluster Jansen & Witherly, 2004 (Fig. 4). The pipe has two distinctive kimberlite units that are embedded in a background consisting of a granitic basement covered by a thin layer of till (Fig. 4a). The first pipe unit is a pyroclastic and volcanoclastic kimberlite (called PK/VK), which has a weak magnetic susceptibility and a very high negative density contrast. The second unit is a hypabyssal kimberlite (called HK), which has a strong magnetic susceptibility and a weak negative density contrast. The Tikhonov inversions of the field gravity and magnetic datasets have been documented in Devriese et al. (2017).

For this example, we use simulated surface gravity and airborne magnetic data modelled from a synthetic model of the DO-27 pipe. The forward and inversion mesh is a tensor mesh with 375 442 active cells; each has a pair of density-magnetic susceptibility values. The smallest cells are cubes with a 10 m edge length. All chosen values for the surveys and geological units are consistent with observations documented in Devriese et al. (2017). For the PK/VK unit, we assume a magnetic susceptibility of $5\cdot10^{-3}~\text{SI}$ and a density contrast with the background of $-0.8~\text{g/cm}^3$. For the HK unit, the magnetic susceptibility is set to $2\cdot10^{-2}~\text{SI}$ and the density contrast to $-0.2~\text{g/cm}^3$ (Fig. 4b). We forward modelled the data over a grid of 961 receivers, at the surface for the gravity survey and at the height of 20 m for the airborne magnetic survey (Fig. 4c and d). We added unbiased Gaussian noise to the gravity and magnetic data with standard deviations of 0.01 mGal and 1 nT, respectively. These standard deviations are input into the data weighting matrices $\left\{{\mathbf{W}_d^k, ~k=1, 2}\right\}$.

For each inversion, we added bound constraints so that the sought density contrast values are null or negative, and the magnetic susceptibility contrast values are null or positive. We used the sensitivity of each survey to define the $\left\{\mathbf{w}_{ip}, ~i=1..n, p=1..q\right\}$ weights. Each physical property is weighted by the sensitivity of its associated survey. Sensitivity-based weighting is a common practice for potential fields inversions Li & Oldenburg, 1996Li & Oldenburg, 1998Portniaguine & Zhdanov, 2002Mehanee et al., 2005. The initial model is the background half-space for all inversions.

6.2Tikhonov inversions¶

We first run the individual inversions of the gravity and magnetic data using the well-established Tikhonov approach described in Section 2.2. The results, shown in Fig. 5, are relatively smooth. The gravity inversion (Fig. 5a) provides an approximate outline of the pipe. The magnetic inversion (Fig. 5b) shows a body centred on HK, but it is too diffuse to delineate a shape. Fig. 5(c) shows the cross-plot of the recovered density and magnetic susceptibility models. Each point is coloured based on both its density and magnetic susceptibility values: white when both contrasts are low, with a blue-scale for a significant density contrast only, with a red-scale for only a significant magnetic susceptibility, and with a purple-scale when both contrasts are significant. 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 strongest density contrasts being co-located with the highest magnetic susceptibility values (mesh cells coloured in purple in Fig. 5c and d); this is in contradiction with the setup where PK/VK, the unit with a low density, is distinct from HK, the unit with high magnetic susceptibility. In both gravity and magnetic inversions, the two anomalous units are indistinguishable from each other. To highlight this, we show an overlap of the two inversions in Fig. 5(d). We coloured each point relative to its density and magnetic susceptibility values, as in Fig. 5(c). This juxtaposition highlights that combining both models does not show structures that seem closer to the ground truth. Post-inversions classification would give highly variable results, depending on the thresholds chosen to delineate units.

We next move to a PGI approach and include petrophysical information. We start by inverting each geophysical dataset individually, and we assess what gains are made before moving to a multi-physics PGI.

6.3Single-physics PGIs¶

We apply the PGI framework developed in Astic & Oldenburg (2019) to invert each geophysical dataset individually. The results are shown in Fig. 6. For the prior petrophysical distribution, we use the true value for the means of each unit. For the petrophysical noise levels, we assign standard deviations of $3.5\%$ of the highest known mean value for each physical property except for the background for which we assign $1.75\%$ (see the 1D left and bottom panels for the distribution of each physical property, respectively, in Fig. 6c). For the proportions, we also used the true values. We acknowledge that proportion values could be difficult to estimate in practice. However, in our experiments, the values of the global proportions have not had a significant impact on the inversion result. The use of locally varying proportions can, however, guide the reproduction of particular features Giraud et al., 2017Astic & Oldenburg, 2019. All the GMM parameters are held fixed in those single-physics inversions.

In carrying out the inversions, we found that both magnetic and gravity data can be explained individually by assuming a single unit, either PK/VK or HK. Each dataset, gravity or magnetic, can be fit by either reproducing the signature of PK/VK or HK, or any value in between. The difference in physical property contrast is compensated for by a difference in the volume of the recovered body. The two kimberlite facies are thus indistinguishable when we consider one geophysical survey at the time. Adding a third cluster in either inversion to represent the second kimberlite facies does not help, as it only gives the algorithm more "choices" that are not supported by the data. For conciseness, we choose to show here the gravity result recovered using only the petrophysical signature of the PK/VK unit (Fig. 6a), which is the most responsible for the gravity anomaly; for the 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. 6b). The additional models (gravity inversion with HK’s density signature and magnetic inversion with PK/VK’s magnetic signature) are shown in Appendix 5; these demonstrate the discrepancy in the recovered volumes of each unit between the magnetic and gravity inversions. For example, explaining the gravity anomaly with only a body with the same density as HK leads to a very large body, 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. 6a) gives useful information about the depth and delineation of the pipe that was not available from the Tikhonov inversion. The magnetic PGI using the HK petrophysical signature (Fig. 6b) places a body around the HK unit location but misses its elongated shape. From the petrophysical perspective, both the gravity signature of PK/VK and the magnetic signature of HK are individually well reproduced. However, the combination of the density and magnetic susceptibility contrasts recovered by the individual PGIs (cross-plot in Fig. 6c) is very far from the desired multidimensional petrophysical distributions. Even by assuming just two units for each inversion (background and kimberlite) as we did, there are still four different combinations of density and magnetic susceptibility values. In this specific case, looking at Figs 6c) and d), there is: 1) a cluster representing the background with both weak density and magnetic susceptibility contrasts (coloured in white); 2) a cluster with a large density contrast and a low susceptibility (coloured in blue); this is close to the petrophysical signature of the PK/VK unit; 3) a cluster with high magnetic susceptibility and very small density contrasts (coloured in orange); This would be the HK unit; 4) a cluster that has both high magnetic susceptibility and large density contrasts, that we identify as ’undefined’ in the figures. This last cluster does not correspond to any unit signature and occupies a large volume. This hinders our ability to resolve two clear kimberlite facies from the inversions. Therefore, this motivates us to move to a multi-physics inversion approach to take advantage of the density-magnetic susceptibility relationships in the inversion and finally delineate two distinct kimberlite facies.

6.4Multi-physics PGI with petrophysical information¶

We now apply our multi-physics PGI framework to jointly invert the gravity, magnetic data with petrophysical information (means, covariances, proportions for each unit, background, PK/VK, HK, of the GMM). The parameters $\left\{\mathbf{w}_{ij}\right\}_{i=1..n, j=1..c}$ are again used to include the appropriate sensitivity weighting for each method and physical property. We use the same uncertainties that we used for the individual PGIs. The off-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 still held fixed. The results are shown in Fig. 7.

The improvement is significant. The joint inversion succeeds in recovering two distinct kimberlite facies that reproduce the provided petrophysical signatures. The quasi-geology model (Fig. 7d) is geologically consistent and does not introduce erroneous structures. The surface outline of the pipe is well recovered. The vertical extension is similar to that of the true model. We also now have indications of the elongated shape and tilt of the HK unit. The magnetic susceptibility of HK is slightly overestimated but still within the acceptable margins defined by the petrophysical noise levels we assigned (Fig. 7c).

To illustrate the behaviour of our heuristic strategy for the update of the objective function scaling parameters, we provide in Fig. 8 the convergence curves of the three data misfits (gravity, magnetic, and petrophysical), and the evolution of the various dynamic scaling parameters, for the multi-physics PGI with full petrophysical information. All target values are reached after 25 iterations, and the PGI stops.

This result was obtained by providing the petrophysical means of the rock units in the GMM. In the next inversion, we devise our approach for using the multi-physics PGI framework when quantitative information is not available.

6.5Multi-physics PGI with limited information¶

We have illustrated the gains made by the multi-physics PGI framework when extensive and quantitative a priori information is provided. We now investigate how to perform multi-physics inversions when a priori information to design the coupling term is not available. Astic & Oldenburg (2019) emphasized the benefits of learning a GMM during the inversion process to compensate for uncertain, or unknown, petrophysical information. At each iteration, the GMM parameters are determined by running a Maximum A Posteriori Expectation-Maximization (MAP-EM) clustering algorithm Dempster et al., 1977. The MAP-EM algorithm estimates compromise values for the GMM parameters based on the prior GMM parameters, weighted by confidence parameters in this prior knowledge, and the current geophysical model. We generalize the learning process of the GMM parameters to a multidimensional case in Appendix 2.

In the next example, we demonstrate how learning the means of the kimberlite units iteratively through the inversion allows us to still perform multi-physics inversions without providing physical property mean values. This is done by acting on the $\{\mathbf{\kappa}\}$ confidence parameters in the means. A confidence κ value of zero indicates that the mean is fully learned from the inversion, while an infinite confidence fixes the mean to its prior value. In all the inversions with limited information, we fix the means of the background for both physical properties to their true values (zero); this is a usual assumption in Tikhonov inversions of potential fields that the background has a zero contrast. We keep the covariances of the GMM fixed and similar to what we used previously. The covariance matrices define our petrophysical noise levels and how spread each petrophysical signature can be.

6.5.1Employing qualitative information with PGI¶

Once the Tikhonov inversions have been run (see Fig. 5), it already appears likely that the gravity and magnetic anomalies are mostly generated by two distinct bodies, as the centres of the two recovered anomalous bodies (density and magnetic susceptibility) are at different locations. Lacking quantitative petrophysical information, the multi-physics PGI framework allows us to formulate the following "interpreter’s assumption": one kimberlite unit is responsible for the gravity anomaly, while a second one is responsible for the magnetic response. Employing this assumption is made possible in our framework by defining the confidences in the means of each unit $\{\mathbf{\kappa}\}$ as vectors. This allows us to act on each physical property mean value of each unit. For example, for the kimberlite unit that is assumed to be responsible for the magnetic response, we set the confidence κ in its magnetic susceptibility to zero; the MAP-EM algorithm decides its value at each iteration based solely on the current geophysical model. On the contrary, its mean density contrast is kept fixed at zero by setting the confidence κ in this mean value to infinity. The same procedure is applied for the kimberlite unit that is assumed to be responsible for the gravity response. The initialization of the density contrast and magnetic susceptibility mean values, for the kimberlite units responsible of the gravity and magnetic response respectively, has little impact on the inversion result, so long as the initial guess is reasonable. It is also common practice, in general, to run clustering algorithms multiple times from various initializations before choosing a specific outcome Dempster et al., 1977Murphy, 2012. In the specific result shown in Fig. 9, the density contrast and magnetic susceptibility mean values for the respective kimberlite rock units were initialized at $-1 \text{ g/cm}^3$ and 0.1 SI. Similar results were obtained with other initializations ($-0.4 \text{ g/cm}^3$ and 0.01 SI etc., but not with $0 \text{ g/cm}^3$ and 0 SI for all units). Because of the weak dependency of the result with regard to the initialization, we choose not to show the initial value in Fig. 10, which presents the evolution of the means throughout the multi-physics PGI with qualitative information.

The result of the multi-physics PGI, with no petrophysical information but with the assumption of distinct low-density and magnetized units, is shown in Fig. 9. Three distinct clusters (background, low-density unit, magnetized unit) are well recovered. The final learned means are respectively $-0.33 \text{ g/cm}^3$ for the low-density unit and $1.1 \cdot10^{-2}$ SI for the magnetized unit.

This multi-physics inversion has several advantages over any of the single-physics inversions. First, by bringing in a qualitative, geologic assumption, we are able to delineate two units by avoiding the overlap of low density and high susceptibility anomalies. Second, we get a sense of the dip of the HK unit. None of those two achievements was reached by the Tikhonov inversions or the single-physics PGIs, even with petrophysical information.

The means of the GMM are learned iteratively following the constraints defined by our assumptions: the background has a fixed contrast of zero in both physical properties, one rock unit is responsible for the gravity response with a null magnetic contrast, and one rock unit is responsible for the magnetic response with a null density contrast. The evolution of the estimations of the means throughout the inversion is shown in Fig. 10. All target values are reached after 36 iterations, and the PGI stops.

Finally, we note that distinguishing two bodies, each mostly responsible for a particular geophysical response, is not automatically required from the geophysical datasets themselves. Indeed, we present in Appendix 5, Fig. 4 a multi-physics PGI with no petrophysical information and only two clusters, where the confidences $\{\mathbf{\kappa}\}$ are all zeros for the kimberlite unit; the background mean is kept fixed at zero contrast for both physical properties. A single anomalous body is able to fit both gravity and magnetic datasets, with a learned mean and an acceptable spread according to the set covariance matrices. This further highlights the gains made possible by the opportunities to simply incorporate a qualitative, geologic assumption within the inversion framework.

6.6DO-27 example summary¶

From the petrophysical perspective, the density-magnetic susceptibility cross-plot for the standard Tikhonov inversions is very different from the expected distribution (Fig. 7c). The smoothness of the recovered models and physical property distributions does not allow us to delineate and distinguish between the two kimberlite facies clearly. Using PGI, individual datasets can both be reproduced using a single kimberlite facies. The individual gravity PGI gives us more information about the depth and delineation of the main PK/VK body. The individual magnetic PGI yields a reasonable estimate for the depth of the HK unit but misses its elongated shape. The two individually recovered quasi-geology models are, however, incompatible when they are combined because of the significant overlap of the recovered PK/VK and HK units. The multi-physics PGI without petrophysical information produces a geological model that distinguishes between the two kimberlite facies. It also begins to give us information about the elongated shape and dip of the HK unit; this result was not achieved by any of the single-physics inversions, not even by the ones that included petrophysical information. However, the accuracy of the boundaries of the bodies is affected by the lack of petrophysical information. The result is improved by providing petrophysical information to the multi-physics PGI, which yields our best recovered model.

7Discussion¶

We have expanded the PGI framework developed in Astic & Oldenburg (2019) to carry out multi-physics joint inversions, and we have 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 data from multiple surveys as well as petrophysical data. Finally, we have used a synthetic example to demonstrate the capabilities of the multi-physics PGI framework.

With regards to the iterative learning of GMM means when limited information is available, considering the confidences $\{\mathbf{\kappa}\}$ as vectors is an important contribution of our framework and it advances the approach of Sun & Li (2016). In their approach, the updates to the means are controlled per unit only, without differentiating the physical properties that are well-known from the undocumented ones. They can either learn the means of a unit or keep it fixed, whereas the framework we present is capable of learning specific components of the GMM means for each unit, as demonstrated in the DO-27 example.

We demonstrated examples of multi-physics inversions with potential fields, which are linear problems. In previous works Astic & Oldenburg, 2019, the PGI approach was applied to nonlinear electromagnetic problems (magnetotelluric, direct-current resistivity, and a field frequency-domain electromagnetic dataset), but it considered only individual surveys depending on a single physical property. We plan to implement this approach for performing multi-physics inversions with electromagnetic methods. This will also be an opportunity to test the robustness of our reweighting strategies for multi-physics inversions with nonlinear geophysical problems. Areas such as the DO-27 kimberlite pipe Devriese et al., 2017Fournier et al., 2017Kang et al., 2017, with many different types of geophysical surveys available, are prime candidates for the application of the PGI approach to refine the image of the subsurface structures by integrating more datasets and physical properties in a single inversion.

As we apply the PGI framework to more complex problems, the handling of various types of relationships between physical properties is required. Linear relationships are straight-forward to implement with Gaussian distributions through the covariance matrix, which can define tilted, elongated probability distributions. In Appendix 1, 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 modest in size, our example in Fig. 2 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 2).

Our PGI framework is composed of three regularized optimization problems (Fig. 1). Astic & Oldenburg (2019) laid the mathematical foundation of the framework, with an emphasis on the inclusion and the learning of the petrophysical signatures (Process 2 in Fig. 1). In the current study, we focused on the coupling of several geophysical surveys with various physical properties, thus extending the Process 1 in Fig. 1. The third and last process, the geological identification, is still one where there is much room for advancement. Astic & Oldenburg (2019) showed, with a direct-current resistivity example, the efficacy of the proportions when they are locally set to zero or unity. Proportions values of zero or unity are "constraining," in the sense that they forbid local occurrences of certain units, rather than just favouring it. While intermediate values of the global proportions appear to have minimal impact on the inversions in our experiments, further studies are required to address the importance and the effects of intermediate (strictly between zero and unity) local proportions. The approach taken by Giraud et al. (2017) combines local proportions from stochastic geological modelling with various warm-started initial models. The integration, extrapolation, and learning of geological information, is part of our current active research. More types of geological information, such as dips, contacts or strikes, also need to be formalized within our framework approach. Combining the PGI smallness term with approaches including prior structural information in the smoothness terms Fournier & Oldenburg, 2019Lelièvre & Oldenburg, 2009Brown et al., 2012Yan et al., 2017Giraud et al., 2019 is also to be investigated.

This framework has been implemented as part of the open-source SimPEG (Simulation and Parameter Estimation in Geophysics) project Cockett et al., 2015. This has two major advantages. First, this enables the reproducibility of the approach by making freely available online, in GitHub repositories, the software environment and Python scripts to recreate the example (https://github.com/simpeg-research/Astic-2020-JointInversion, Astic (2020)). Second, by providing a common environment, it allows the implementation of the framework to be readily used with any type of geophysical surveys (such as EM surveys) or discretization (such as OcTree meshes) that are supported in the source code, and facilitate the collaboration with others researchers Oldenburg et al., 2019.

8Conclusion¶

We have expanded the PGI framework to use petrophysical and geological information, represented as a Gaussian Mixture Model, as a coupling term to perform multi-physics joint inversions. We described our strategies for handling multiple geophysical target misfits as well as a petrophysical target misfit. We presented our efforts to make the implementation modular, extensible, and shareable. Finally, we demonstrated, through the DO-27 kimberlite pipe 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 field datasets allowed us to delineate two kimberlite facies and to reproduce their petrophysical signature.