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

Thibaut Astic; Lindsey J. Heagy; Douglas W Oldenburg

doi:10.1093/gji/ggaa378

Article

Contents

1Working with nonlinear petrophysical relationships¶

Figure 1:GMM with various polynomial relationships: one linear (no addition required), one quadratic and one cubic. In the main panel, thicker contour lines indicate higher probabilities. 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. 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:

d_j = \int_0^1 e^{-jx}\cos(2\pi j x) m(x) dx, ~j=1, 3, .., 59.

(1)

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.

Figure 2: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.

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.

\mathcal{P}(\Theta|\mathbf{m}) \propto \mathcal{M}(\mathbf{m}|\Theta)\mathcal{P}(\Theta).

(2)

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:

n_{ij}^{(k)} = \frac{\mathcal{P}(z_i=j)^{(k-1)}\mathcal{N}(\mathbf{m}_i|{\mathbf{\mu}_j}^{(k-1)}, {\mathbf{\Sigma}_j}^{(k-1)})}{ \sum_{t=1}^c \mathcal{P}(z_i=t)^{(k-1)} \mathcal{N}({\mathbf{m}}_i|{\mathbf{\mu}_t}^{(k-1)}, {\mathbf{\Sigma}_j}^{(k-1)})} .

(3)

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

{\pi}^{(k)}_j = \frac{V_{j}^{(k)}+\zeta_j {{\pi}_j}_{\text{prior}}V}{V(1+\sum_{t=1}^c \zeta_t {{\pi}_t}_{\text{prior}})} ,

(4)

with:

V_{j}^{(k)} = \sum^n_{i=1} v_i n_{ij}^{(k)},

(5)

and

V=\sum^n_{i=1} v_i,

(6)

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$ :

{\mu_j^p}^{(k)}=\frac{V_{j}^{(k)}{\bar{{m}^p_j}}^{(k)} + \kappa^p_j {\pi_j}_{\text{prior}} V {\mu^p_j}_{\text{prior}}}{V_{j}^{(k)}+\kappa_j^p {\pi_j}_{\text{prior}} V} ,

(7)

with:

{\bar{\mathbf{m}}}_j^{(k)} = \frac{\sum^n_{i=1} v_i n_{ij}^{(k)} \mathbf{m}_i}{V_{j}^{(k)}},

(8)

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$ :

{\mathbf{\Sigma}_j}^{(k)} = \frac{{{V_{j}^{(k)}} {\mathbf{\Sigma}_{\bar{\mathbf{m}}}}_j}^{(k)} + \nu_j {\pi_j}_{prior} V {\mathbf{\Sigma}_j}_{prior}} {{V_{j}^{(k)}} + \nu_j {\pi_j}_{prior} V} ,

(9)

with:

{\mathbf{\Sigma}_{\bar{\mathbf{m}}}}_j^{(k)} =\frac{1}{{V_{j}^{(k)}}} \sum_{i=1}^{n} v_i n_{ij}^{(k)}(\mathbf{m}_i-\bar{\mathbf{m}_j}^{(k)})(\mathbf{m}_i-\bar{\mathbf{m}_j}^{(k)})^\top,

(10)

3Algorithm¶

Algorithm 1

InitializationInitialization

PGI extended from Astic & Oldenburg (2019) for joint inversion

Input:
- Initial geophysical model $\mathbf{m}^{(0)}$ , GMM parameters $\Theta^{(0)}$ and geological model $\mathbf{z}^{(0)}$ .
Parameters:
- Objective function: data’s noise $\left\{{\mathbf{W}_d}_p\right\}_{p=1..q}$ , $\beta^{(0)}$ and $\left\{\alpha\right\}$ .parameters.
- Localized prior: specific $\mathcal{P}(z_i)$ for available locations $i \in \{1..n\}$ , weights $\left\{\mathbf{w}\right\}$ .
- GMM prior weights: $\left\{\mathbf{\kappa}_j, \nu_j, \zeta_j\right\}_{j=1..c}$ for the means, variances and proportions.
- Optimization: β-cooling factor γ ( $> 1$ ), sufficient decrease rate τ ( $\leq 1$ ), tolerance on target misfit ε.
Output:
- $\mathbf{m}$ , Θ, $\mathbf{z}$ .

while $\text{any}({{\Phi_d^k}}>{{\Phi_d^k}^*}, k=1..r)$ and $\Phi_{\text{petro}}>\Phi_{\text{petro}}^*$ do

Objective Function Descent Step:

Compute a model perturbation $\delta \mathbf{m}$ using an inexact Gauss-Newton.
Line search with Wolfe condition to find an η that satisfy a sufficient decrease of Φ.
Return $\mathbf{m}^{(t)} = \mathbf{m}^{(t-1)}+\eta \delta \mathbf{m}$ .

Update Petrophysical Distribution

Fit a new GMM $\Theta^{(t)}$ on $\mathbf{m}^{(t)}$ such as in equations (4), (7) and (9) until no sufficient increase of the posterior is observed.

Classification:

Compute the membership $\mathbf{z}^{(t)}$ of the current model $\mathbf{m}^{(t)}$ as in equation (10) using $\Theta^{(t)}$ .
Update $\mathbf{m}_{\text{ref}}$ and $\mathbf{W}_s$ according to equations (11) and (12) respectively using $\mathbf{z}^{(t)}$ .

Update regularization weights:
$\text{all}({\Phi_d^k}^{(t)}\geq\text{max}((1+\epsilon){\Phi_d^k}^*, \tau {\Phi_d^k}^{(t-1)}), k=1..r)$ Decrease β: $\beta^{(t)}=\frac{\beta^{(t-1)}}{\gamma}$ . $\text{all}({\Phi_d^k}^{(t)}\leq{\Phi_d^k}^*, k=1..r)$ and $\Phi_{\text{petro}}>\Phi_{\text{petro}}^*$ Increase $\alpha_s$ : $\alpha_s^{(t)}=\alpha_s^{(t-1)} \times \text{median}(\frac{{{\Phi_d^k}^*}}{{{{{\Phi_d^k}}}^{(t)}}}, k=1..r)$ (equation (19)). (optional) $\text{all}({\Phi_d^k}^{(t)}\leq{\Phi_d^k}^*)$ and $\Phi_{\text{petro}}>\Phi_{\text{petro}}^*$ and $\mathbf{z}^{(t)}==\mathbf{z}^{(t-1)}$ Include $\mathbf{m}_{\text{ref}}$ in Smoothness.

Update geophysical data misfit weights:
$\text{any}({\Phi_d^k}^{(t)}\leq {\Phi_d^k}^*, k=1..r)$ update $\left\{\chi\right\}$ according to equations (20) to (22).

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¶

(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. %sc); (b) gravity PGI using the HK density contrast signature. Note again the larger volume compared to Fig. %s(d). — Figure 3:(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. 6c); (b) gravity PGI using the HK density contrast signature. Note again the larger volume compared to Fig. 6(d).

Figure 4:Results of the multi-physics PGI without providing the means of the physical properties for the kimberlite facies; a single kimberlite body is enough to meet the petrophysical requirements (the spread set by the covariances matrices) and reproduce the geophysical datasets; (a) Plan map, East-West and North-South cross-sections through the recovered density contrast model; (b) Plan map, East-West and North-South 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 and side panels, we show the learned petrophysical GMM distribution; (d) Plan map, East-West and North-South cross-sections through the resulting quasi-geology model.

References¶

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

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

[cite-Li1996] Li, Y., & Oldenburg, D. W. (1996). 3-D inversion of magnetic data. Geophysics, 61(2), 394–408. 10.1190/1.1443968

[cite-ggz389] 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

[cite-ExpectationMaximization] 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

[cite-Murphy2012] Murphy, K. P. (2012). Machine Learning: A Probabilistic Perspective. The MIT Press.

[cite-SeogiMapping] 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