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

## 1Working with nonlinear petrophysical relationships¶

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. 1). This is achieved by composing the Gaussian function with the nonlinear relationship $P_j$ between the physical properties of the particular rock unit $j$. This corresponds to using $\mathcal{N}(P_j(\mathbf{m}_i)|\mathbf{\mu}_j, \mathbf{W}_i^{-1}\mathbf{\Sigma}_j\mathbf{W}_i^{-1})$ for each rock unit $j$ in the GMM in equation (8).

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 is the same and is based on the example previously used in Li & Oldenburg (1996). The models are discretized on a 1D mesh defined on the interval $[0, 1]$ and divided into 100 cells. Both datasets consist of 30 data points at various frequencies equally distributed from 1 to 59 evaluated according to:

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

We invert the two datasets using the multi-physics PGI framework, including *a priori* knowledge of the petrophysical distributions and relationships. The result can be seen in the first row of Figs 2(a) to (c). Note that the petrophysical relationships are well reproduced. This allows the recovery of details of the two models.

For comparison, we also jointly invert 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 2(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 2 is completely missed.

We also invert both datasets individually using the classic Tikhonov inversion. The result is shown in the last row of Figs 2(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.

## 2Updating the Gaussian mixture model in multi-dimensions¶

In this section, we generalize the learning of the GMM parameters presented in Astic & Oldenburg (2019) to multivariate Gaussian distributions, representing multiple physical properties. The main difference comes from the fact that the means are now vectors (they are scalars in 1-D) and that the scalar variances in 1D become covariance matrices. We define the problem as a Maximum A Posteriori (MAP) estimate of the GMM means, proportions and covariance matrices, using the MAP Expectation-Maximization (MAP-EM) algorithm Dempster *et al.*, 1977.

We choose to follow a semi-conjugate prior approach for the choice of the prior distributions $\mathcal{P}(\Theta)$ Astic & Oldenburg, 2019Murphy, 2012.

The computation of the responsibilities for the E-step of the MAP-EM algorithm stay the same as in Astic & Oldenburg (2019), except for the dimension of the parameters:

The update to the proportions stays similar to the univariate case Astic & Oldenburg, 2019:

with:

and

where $\zeta_j$ is the confidence in the prior proportion ${\pi_j}_{\text{prior}}$ of the rock unit $j$, $v_i$ is the volume of the $i$^{th} cell and $V$ is the volume of the active mesh. This allows the estimates to be mesh-independent by using volumetric proportions instead of cell counts.

The mean update proposed in Astic & Oldenburg (2019) generalizes to each physical property $p$:

with:

where $\kappa_j^p$ is the confidence in ${{\mu^p_j}_{\text{prior}}}$, which is the prior mean of the physical property $p$ of the rock unit $j$. The confidences $\{\mathbf{\kappa}\}$ are thus consider as vectors. This is an important tool that we use in section 6.5.1 to formulate a geologic assumption about the model.

The update to the covariance matrices is, with $\nu_j$ the confidence in the prior covariance matrix ${\mathbf{\Sigma}_j}_{\text{prior}}$ of the rock unit $j$:

with:

## 3Algorithm¶

## 4Pseudo-code for the implementation in SimPEG¶

Here, we provide an overview for the use of our implementation of the multi-physics PGI framework by other scientists. We outline the main components of the multi-physics PGI implementation using `SimPEG`

and other core tools in the Python ecosystem. As in the paper, we consider the multi-physics inversion of gravity and magnetic data.

We begin by creating a simulation `mesh`

and defining mappings that translate the model vector, which contains all of the parameters we will invert for, to physical properties on the mesh to be used in each forward simulation. The inversion model is a single vector; for an inversion with multiple physical parameters, we stack them and use the `Wires`

map to keep track of which indices in the vector correspond to each physical parameter. default model_mappings.py

Note that mappings can be composed; for example if we wanted to invert for log-susceptibility rather than susceptibility, then we would compose the `wires.susceptibility`

mapping with an `ExpMap`

instance. For examples, see Kang *et al.* (2015).

Next, we construct the forward simulations for both the gravity and magnetic data. The `survey`

objects contain the locations of the receivers as well as parameters defining the source field for the magnetics simulation (magnitude, inclination, and declination). default data_misfits.py

With both gravity and magnetic simulations defined, we now have the ability to compute predicted data given a model. The next step is to construct the data misfit term as described in equation (17). default combodatamisfit.py

The `mag_data`

and `grav_data`

objects contain the observed data as well as user-specified uncertainties, and the scalar `chi`

-values are initialized according to equation (22).

Next, we construct the GMM and regularization which consists of the petrophysical smallness term described in equation (9) and smoothness terms for both the density and susceptibility. To construct the GMM, we use the `WeightedGaussianMixtureModel`

object, which inherits from and extends the `sklearn.mixture.GaussianMixture`

object in Scikit-Learn Pedregosa *et al.*, 2011. The instantiated `gmm`

object has methods for fitting a GMM on a model and computing membership of a model as needed in steps 4 and 5 in Algorithm 1. default regularization.py

Having defined the components of the objective function, we now specify the optimization routine, in this case, inexact Gauss Newton with projections for bound constraints. default optimization.py

Finally, we assemble the inversion using `directives`

in SimPEG to orchestrate updates throughout the inversion. Each directive has an `initialize`

method which is called at the beginning of the inversion and can be used to set initial values of parameters, for example to initialize β and the α-values, as well as an `end_iteration`

method which is called after a model update (at the end of step 3 in Algorithm 1) and can be used to update parameters in the inversion. The updates in steps 4 through 7 in Algorithm 1 make use of the `directives`

functionality. default inverse_problem.py

## 5Additional models obtained by PGIs for the DO-27 synthetic case study¶

- Li, Y., & Oldenburg, D. W. (1996). 3-D inversion of magnetic data.
*Geophysics*,*61*(2), 394–408. 10.1190/1.1443968 - Astic, T., & Oldenburg, D. W. (2019). A framework for petrophysically and geologically guided geophysical inversion using a dynamic Gaussian mixture model prior.
*Geophysical Journal International*,*219*(3), 1989–2012. 10.1093/gji/ggz389 - Dempster, A. P., Laird, N. M., & Rubin, D. B. (1977). Maximum likelihood from incomplete data via the EM algorithm.
*Journal of the Royal Statistical Society, Series B*,*39*(1), 1–38. http://links.jstor.org/sici?sici=0035-9246(1977)39:1%7B%5C%25%7D3C1:MLFIDV%7B%5C%25%7D3E2.0.CO;2-Z%7B%5C&%7Dorigin=MSN - Murphy, K. P. (2012).
*Machine Learning: A Probabilistic Perspective*. The MIT Press. - Kang, S., Cockett, R., Heagy, L. J., & Oldenburg, D. W. (2015). Moving between dimensions in electromagnetic inversions.
*SEG Technical Program Expanded Abstracts 2015*, 5000–5004. 10.1190/segam2015-5930379.1