Joint_Inversion.ipynb

Kernel: Python 3 (system-wide)

Synthetic model - Mansfield area

We will be inverting gravity and magnetic data generated for the a synthetic model inspired from geological mreasurements collected in the Mansfield area (Victoria, Australia) as shown below. Details about the model used here can be found in Giraud et al. (2017) and Giraud et al. (2019) (geophysics and petrophysics) and Pakyuz-Charrier (2017, 2018) (geology).

Density contrast and magnetic susceptility

The density contrast and magnetic susceptibilitity models used here are shown in Figure 1 (top). The corresponding geophysical measurements are shown in Figure 1 (bottom).

Structure term

The entropy used for the structure term constraint in presented in Figure 2.

Geological model

The geological reference model is shown in Figure 1. The corresponding simulated probabilistic geological model is shown in Figure 4. The geological model was derived by Pakyuz-Charrier (2018).

Modelling petrophysical measurements

The petrophysical measurements are modelled using a Gaussian mixture model, as shown in Figure 3. It presents 2 ambiguities. Lithologies 1 and 2, and 5 and 6, respectively, have the same density contrasts. The same applies to lithologies 3 and 5, and 2 and 4, respectively, for magnetic susceptibility.

Constrained potential field inversions

A step-by-step guided process from unconstrained inversion to geologically constrained joint inversion

September 2019, AEGC Perth

Clément Barrière and Jérémie Giraud

Disclaimer

all material is provided without any warranty whatsoever. While we make every effort to ensure that this Jupyter notbeook is accurate and up to date, there remain the possibility of human error, and utilisation of the code provided here may by no means provide exact, certain answers to any modelling problem. We do not guarantee, and accepts no legal liability whatsoever arising from or connected to, the accuracy or completeness of any material contained in this Jupyter notebook or obtained through usage of it. The information here is subject to be changed without any notice/legal procedure. You are reading the text and code provided here and using this notebook at your own risks.

In this notebook we propose to run a series of inversions with increasing constraints. The aim of this exercise is to understand the influence of constraints onto inversion. These constraints range from the simple application to depth weight to counteract the decay of the sensitivity with depth to the application of spatially varying petrophysical constraints in joint inversion.

Introductory note

The Python code proposed in this practical tutorial provides a simplified version of the 3D parallel inversion platform Tomofast-x (Martin et al., 2013, Martin et al., 2018, Giraud et al., 2018, Giraud et al., 2019), which we call 'Tomoslow-2D'. While the implementation the clustering algorithm and the gradient constraints differs, the philosophy remains the same and results obtained with the code shown here present much similar features. It has been developed for educational and prototyping purposes.

The aim of this jupyter notebook is to provide a high-level understanding of the mechanics of constrained and joint inversion codes similar in philosophy to the implementation we follow. In the current implementation, the cell-size is chosen to be unity (1m x 1m), but the problem is scalable to any dimension.

Acknowledgements

The authors acknowledge Mark Jessell's support throughout the project.

       figure 1: True model for the Density contrast and the Magnetic suceptibility with the corresponding data set 

                                      figure 2: Shannon entropy of the synthetic model

        figure 3: Crossplot of Gaussian Mixture distribution of the Density contrast and the magnetic suceptibility

As one of the pre-requisites to run this code, we first need to load libraries for visualisation and numerical computations.

import matplotlib.pylab as plt
from scipy.sparse.linalg import lsqr
import scipy.stats as ss
import numpy
import matplotlib.colors as colors

plt.rcParams.update({'font.size': 16})

Following this, we need to import modules that contain fonctions called by the code during inversion.
The module Calc_Jac contains code that is necessary to calculates the Jacobian matrix of the problem considered.
The module Joint_inversion contains code to calculate depth weighting and functions to apply constraints to inversion.

The file names we define will be used and explained later in the notebook.

%load_ext autoreload
%autoreload 1

%aimport Calc_Jac
%aimport Joint_inversion
%aimport PC

filename_grav = "mod_grav.txt"
filename_mag = "mag_susc_mod.txt"
filename_geol_grav = "geol_weights.txt"
filename_geol_mag = "geol_weights_mag.txt"
filename_geol_gravmag = "geol_only_weights.txt"

FORMULATION OF THE COST FUNCTION

We formulate the inverse problem in a least-squre sense as follows: $\Large\theta (m)= \mid\mid (d-g (m)) \mid \mid^{2} _{2}+ \mid \mid W_{m} (m-m_{p}) \mid \mid^{2} _{2} + \alpha \mid \mid W_{H} \nabla_{m} \mid \mid^{2} _{2} + \mid \mid W_{p} P(m) \mid \mid^{2} _{2} \quad (1)$

Where the different terms mean:

• $m$ is the invered model.

• $d$ refers to the geophysical data.

• $m_{p}$ stands for the prior model.

• $g(m)$ refers to the respective model.

• $\mid \mid(d-g (m)) \mid \mid^{2} _{2}$ represent the data term, (Lines and Treitel, 1984).

• $\mid \mid W_{m} (m-m_{p}) \mid \mid^{2} _{2}$ represent the model term, (Hoerl and Kennard, 1970).

• $\alpha \mid \mid W_{H} \nabla_{m} \mid \mid^{2} _{2}$ is the structure term which stands for local gradient regularization, where $\nabla_{m}$ represent the model gradient and $\alpha$ represent the damping gradient constraints, (Li and Oldenburg, 1996, 1998).

• $\mid \mid W_{p} P(m) \mid \mid^{2} _{2}$ represent the petrophysics term or the the clustering constraint. Here $P(m)$ is the petrophysical distribution which is given by the following equation: $P(m) = log(\Sigma^{n_{f}}_{k=1} \omega_{k}N(m\mid \mu_{k},\sigma_{k}))$ where $\omega_{k} = \dfrac{1}{n_{f}}$ evrywhere if no geology is available or $\omega_{k} = \psi_{k,i}$ in the i-th cell otherwise (Giraud et al, 2017).
  Here, matrices $W_{m}$, $W_{H}$ and $W_{p}$ are spherical matrices.

Solver

We solve the inverse problem using the LSQR algorithm (Paige and Saunders, 1982).
In the following part of the code, A et B refers respectively to a matrix and a vector wich will be used to resolve the $Am=B$ in the least-square sense using the LSQR algorithm (Paige and Saunders, 1982).

Then the code use the lsqr function from the scipy.sparse.linalg library to solve the system of equations $Am=B$ . The LSQR function retruns the solution model called $m_{inv}$

For density contrast inversion, the LSQR function parameters are left to their default values except for the number of iterations wich is set to 100 to have a good compromise between speed and accuracy. The magnetic inversion requires a smaller tolerance to get a better approximated solution. More information about the way the LSQR function is built is available on the website: https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.sparse.linalg.lsqr.html

$Am=B$ can be rewritten as $m=(A^{T}A)^{-1}A^{T}B$. The calculation of an exact solution is limited to small models as calculating $(A^{T}A)^{-1}$ increases rapidly with the size of the model. For large models, the LSQR algorithm provides an approximated solution. In this notebook, the system is quite small so an exact solution could have been calculated but using the LSQR algorithm allow us to work on bigger models if necessary.

1. UNCONSTRAINED GRAVITY INVERSION

                            Figure 4: Probability of the different lithologies of the synthetic model

The simplest inversion we perform is constrained only by the model smallness term. The function to optimize is given as:

$\quad\theta (m)= \mid \mid W_{d} (d-g (m)) \mid \mid^{2} _{2}+ \mid \mid W_{m} (m-m_{p}) \mid \mid^{2} _{2} \quad (2)$

First the code applies a depth weighting (DW) to the Jacobian matrix using a power law (Li and Oldenburg, 1998). This balances the relative importance of the different model-cells with depht via the following equation:
$\large\quad w(z) = \dfrac{1}{(z-z_{0})^{\beta/2}} \quad (3) \quad$,
where $z$ represent the depth of the cell, $\beta$ is usually set accordingly to the physics of the problem, and $z_{0}$ depends on the cell size (Li and Oldenburg, 1996, 1998).

LSQR formulation for least-square inversion

The LSQR algorithms approximates the solution of the following problem, where $m$ is the quantity inverted for:

We can therefore solve Equation (2) by setting:

• The matrix "A" as $\ A = \begin{bmatrix} G \\ W_{grav}I \end{bmatrix}$ with $G$ the Jacobian matrix and $W_{grav}I$ an indentity matrix mulitplied by a specific weight to match the size of the vector $B$.

• The vector "B" as $\ B = \begin{bmatrix} D_{g} \\ m_{p} \end{bmatrix}$ with $D_{g}$ being the data vector and $m_{p}$ the prior model wich initialised to 0.

The equation (2) can now be written as: $\begin{bmatrix} G \\ W_{grav}I \end{bmatrix}m = \begin{bmatrix} D_{g} \\ m_{p} \end{bmatrix}$

In the following lines, the code build these elements.

Note:

Several function used in the code are not detailed here to make the notebook more readable. Theses functions are:

• set_mod creates a python matrix from data stored in a text file.

• grav_jac calculates Jacobian matrix for the gravity inverse problem.

• mag_jac calculate Jacobian matrix for the magnetic inverse problem..

• get_dim gets the dimensions of a matrix stored in a text file.

• posx generates a vector of positions along the horizontal axis used in the calculation of the Jacobian matrix.

• Calc_DW calculates depth weighting. It creates matrices containing the depth weighting values for the gravity and the magnetic models.

• calc_grad calculates the gradients of the model.

• lsqr solves a system of equations in a least-squre framework using the LSQR algorithm (Paige and Saunders, 1982)

pos_grav = Calc_Jac.posx(filename_grav, start=0.5, interval=1)
line , column = Calc_Jac.get_dim(filename_grav)

Gij, Dg = Calc_Jac.grav_jac(filename_grav,pos_grav)
Id_grav = numpy.eye(line*column,line*column)

$Cm^{mag}$ and $Cm^{grav}$ contain the value of the DW that will be used to precondition the Jacobian matrix. From the gravity and magnetic problems, depth weighting is calcualted using equation (3) as follows:

• $\beta = 0.8\quad$ and $\quad z_{0} =0.5\quad$ for gravity (Boulanger and Chouteau, 2001);

• $\beta = 2\quad$ and $\quad z_{0} =0.5\quad$ for magnetic (Li and Chouteau, 1998).

Cm_grav , Cm_mag = Calc_Jac.calc_DW(filename_grav)

k = 0
for j in range(1,line+1):
    for i in range(1,column+1):
        Cm_mag[k][k] = 1/(j-0.5) ** (2/2)
        Cm_grav[k][k] = 1/(j - 0.5) ** (0.8/2)
        k +=1

### Zero model as prior model

#mp = numpy.zeros((line*column,))

### Vertical gradient as prior model

mp = numpy.reshape(Calc_Jac.calc_vertical_gradient(filename_grav),(line*column,))

for i in range(line * column):
    mp[i] *= Cm_grav[i][i]

Depth weighting is then applied to the Jacobian matrix $G$ (Gij):

for i in range(len(Gij)):
    for j in range(len(Gij[0])):
        Gij[i][j] /= Cm_grav[j][j]

As explained above, we will run a serie of inversions with increasing constraints. Under certain circumstances, these constraints might have competing influence on inversion due to the non-uniqueness of the solution and the topography of the cost function. Therefore, we need to implement different weights to strengthen or reduce the effect of a given contraint.
The first weight we add here is grav_w for the unconstrained density contrast inversion.
The next step is to buid the matrix $A$ and the vector $B$: here the matrix $W_{grav}I$ is multiplied by a weight: $grav_w$

grav_w = 1e-7 # Suggested value: 1e-7

run_already = False
A = numpy.concatenate((Gij,numpy.dot(Id_grav,grav_w)))
B = numpy.concatenate((Dg, numpy.dot(mp,grav_w)))

m_inv = lsqr(A, B, damp=0.0, atol=1e-08, btol=1e-08, conlim=1e08, iter_lim=100)

for h in range(line * column):
    m_inv[0][h] /= Cm_grav[h][h]

After inversion we use a number of plotting function for quality control and visualisation from the pyplot library:

We also define the data misfit and the model misfit as follow:

$\large data \quad misfit = \sqrt{\dfrac{\Sigma (d_{i}^{obs}-d_{i}^{calc})^{2}}{\Sigma (d_{i}^{obs})^{2}}}\quad$ $\large \quad \quad model \quad misfit = \sqrt{\dfrac{\Sigma (m_{i}^{true}-m_{i}^{inv})^{2}}{\Sigma (m_{i}^{true})^{2}}}$

These two values allow us to keep track of the difference between the inverted models and the true one.

rm_inv = numpy.reshape(m_inv[0], (line,column))
misfitgrav = Calc_Jac.model_misfitgrav(m_inv[0], filename_grav)
props = dict(boxstyle='round', facecolor='white', alpha=0.9)

plt.figure(figsize=(15,10))
plt.imshow(rm_inv, cmap='jet')
cbar = plt.colorbar(shrink=0.5)
cbar.set_label('$Kg.m^{-3}$', labelpad=-35, y=1.15, rotation=0)
plt.title('Unconstrained Density contrast model')
plt.xlabel('Horizontal position')
plt.ylabel('Depth')
plt.text(0.5, 2.5, r'Model misfit = ' + str(misfitgrav), bbox=props)
plt.show()

if not run_already:
    for i in range(len(Gij)):
        for j in range(len(Gij[0])):
            Gij[i][j] *= Cm_grav[j][j]
        
mul = numpy.matmul(Gij, m_inv[0])

rms = Calc_Jac.calc_RMS(Gij,Dg, m_inv[0])
        
plt.figure(figsize=(15,10))
plt.title('Gravity anomaly')
line_up, = plt.plot(mul, label='Line 1')
line_down, = plt.plot(Dg, label='Line 2')
plt.legend([line_up, line_down], ['Calculated data', 'Measured data'])
plt.xlabel('Horizontal position')
plt.ylabel('Gravity anomaly (gal)')
plt.grid(b=True, which='major', color='0.6', linestyle='-')
plt.minorticks_on()
plt.grid(b=True, which='minor', color='0.95', linestyle='--')
plt.text(20, 0.00225, r'rms = ' + str(rms), bbox = props)
plt.show()

run_already = True

Exercise: What happens if you change grav_w ?

2. ADDING THE STRUCTURE TERM

In this section, we utilise the information entropy (Shannon, 1948, Wellmann and Regenauer-Lieb, 2012) to adjust the strength of the smoothness term locally. Schweizer et al. (2017) define information entropy as follow: information entropy "quantifies the amount of missing information and hence, the uncertainty at a descrite location". For the $i^{th}$ cell, it is given by (Shannon, 1948):

$$\large H^{i} = H(\psi^{i}) = -\Sigma_{k=1}^{L} \psi^{i}_{k}log(\psi^{i}_{k}) \quad (4)$$

Here, we normalise $H$ and calculate $W_{H}$ such that it reaches its maximum values where uncertainty is high and conversely where entropy is low. We calculate $W_{H}$ (see eq. 1) as follows (Giraud et al, 2019):

$$\large W^{i}_{H}=\dfrac{max H-H^{i}}{maxH-minH} \quad (5)$$

In the following lines, we compute the gradients of the model in the horizontal and vertical directions ($grad_{x}$ and $grad_{y}$, respectively) using a forward difference scheme. We then perform element-wise multiplication of the gradients with $W_{H}$ to adjust the strength of the smoothness constrained applied during in inversion. This process relaxes the smoothness constraints in uncertain parts of the model while maintaining it is where uncertainty is lower.
We also define $w_{entropy}$ (w_entropy) to weight this structure (i.e. smoothess) term constraint in the cost function.

After this, we build the matrix A and vector B. We build A by concatenating the Jacobian matrix with w_entropy multiplying $grad_{x}$ and $grad_{y}$ and B by concatenating $d$ with all-zeros vectors indicating that the smoothness term should be minimized.

w_entropy = 5e-7 # Suggested value: 5e-7
gradx, grady = Joint_inversion.calc_grad(filename_grav, filename_geol_grav, filename_geol_mag, filename_geol_gravmag, mag = False,verbose=False)
empty_vec = numpy.zeros((line*column,))

Gij, Dg = Calc_Jac.grav_jac(filename_grav,pos_grav)

A = numpy.concatenate((A, numpy.dot(gradx,w_entropy), numpy.dot(grady,w_entropy)))
B = numpy.concatenate((B, numpy.dot(empty_vec,w_entropy), numpy.dot(empty_vec,w_entropy)))

m_inv = lsqr(A, B, damp=0.0, atol=1e-08, btol=1e-08, conlim=1e08, iter_lim=100)
for h in range(line * column):
    m_inv[0][h] /= Cm_grav[h][h]
    
rm_inv = numpy.reshape(m_inv[0], (line, column))
misfitgrav = Calc_Jac.model_misfitgrav(m_inv[0], filename_grav)
props = dict(boxstyle='round', facecolor='white', alpha=0.9)

plt.figure(figsize=(15,10))
plt.imshow(rm_inv, cmap='jet')
cbar2 = plt.colorbar(shrink=0.5)
cbar2.set_label('$kg.m^{-3}$', labelpad=-35, y=1.15, rotation=0)
plt.title('Density contrast model with structure term')
plt.xlabel('Horizontal position')
plt.ylabel('Depth')
plt.text(0.5, 2.5, r'Model misfit = ' + str(misfitgrav), bbox=props)
plt.show()

mul = numpy.matmul(Gij, m_inv[0])

rms = Calc_Jac.calc_RMS(Gij,Dg, m_inv[0])

plt.figure(figsize=(15,10))
plt.title('Gravity anomaly')
line_up, = plt.plot(mul, label='Line 1')
line_down, = plt.plot(Dg, label='Line 2')
plt.legend([line_up, line_down], ['Calculated data', 'Measured data'])
plt.xlabel('Horizontal position')
plt.ylabel('Gravity anomaly (gal)')
plt.grid(b=True, which='major', color='0.6', linestyle='-')
plt.minorticks_on()
plt.grid(b=True, which='minor', color='0.95', linestyle='--')
props = dict(boxstyle='round', facecolor='white', alpha=0.5) 
plt.text(20, 0.00225, r'rms = ' + str(rms), bbox = props) 
plt.show() 

Exercise: What happens if you change w_entropy ?

3. ADDING PETROPHYSICAL CONSTRAINTS

In this part of the code we add the petrophysical constraints term to the cost function (eq. 1):

$\theta (m)= ... + \mid \mid W_{p} P(m) \mid \mid^{2} _{2}$

Single domain inversions

Gravity modelling

As shown in Figure 1, lithologies [1 and 2] and [5 and 6] have the same density contrast. We can therefore model density contrasts using a mixture of 4 Gaussians. The individual Gaussians are parameterised as:

• mean values of: 0, 100, 200 and 300 $kg/m^3$

• standard deviation of 40 $kg/m^3$.

Magnetic modelling

Using the same logic as for gravity, lithologies [2 and 4] and [3 and 5] have the same magnetic susceptibility. We can therefore model magnetic susceptibility with a mixture of 4 Gaussians. The Gaussians are parameterised as:

• mean values of: 0, 0.025, 0.05 and 0.0075 SI

• standard deviation of 0.0001 SI (suggested value as starting point for tests in this notebook)

An exemple of such Gaussian mixture model is presented on the figure 5.

    Figure 5: Exemple of a Gaussian mixture model for the density contrast and the magnetic suceptibility in 1D and 2D

Use of probabilistic geological modelling

The results of probabilistic geological modelling we use consists of the observation probability of lithologies. This is useful to adapt the weight of the different lithologies in the Gaussian mixture model accordingly with their observation probability.

Implementation

We optimize of the petrophysical constraints term using a k-nearest neighbour algorithm, in which the distance between the different cluster is calculated using the Gaussian mixture model (Mahalanobis distance). The implementation is explained below.

In what follows we perform inversion petrophysical constraints. We start with gravity inversion, followed by the inversion of magnetic data.

means_grav = [0,100,200,300]
standard_dev_grav = [40,40,40,40]

The following lines compute a first resulting model with $A$ and $B$ from part 1. We use this model as a prior model to compute the value of the Gaussian mixture. It is represented by the term pdfs, calculated by the function PDF.

line, column = Calc_Jac.get_dim(filename_grav)

### Zero model as prior model

#mp = numpy.zeros((line*column,))

### Vertical gradient as prior model

mp = numpy.reshape(Calc_Jac.calc_vertical_gradient(filename_grav),(line*column,))

for i in range(line * column):
    mp[i] *= Cm_grav[i][i]

Id_grav = numpy.eye(line*column,line*column)

for i in range(len(Gij)):
    for j in range(len(Gij[0])):
        Gij[i][j] /= Cm_grav[j][j]

A = numpy.concatenate((Gij,numpy.dot(Id_grav,grav_w)))
B = numpy.concatenate((Dg, numpy.dot(mp,grav_w)))

m_inv = lsqr(A, B, damp=0.0, atol=1e-08, btol=1e-08, conlim=1e08, iter_lim=100)

for h in range(line * column):
    m_inv[0][h] /= Cm_grav[h][h]


xs = m_inv[0]
nb_gauss, pdfs = PC.PDF(xs, means_grav, standard_dev_grav)

Next, the code enters a loop described below.

For each cell of the model we determine the maximum of the 4 Gaussians. For each model-cell, the value of the center of the best matching Gaussian in the mixture is stored in the vector $P_{vector}$ (P_vector). This vector is then multiplied by $weight$ (weight). When geological information is availabe, this value depending on the the probability of the different lithologies. An identity matrix ($pfdp$) (pfdp) is also created to match the size of $P_{vector}$. In this part of the code, we define $pgrav_w$ (pgrav_w) to adjust the strength of the petrophysical constraint during inversion.

Then we concatenate $A$ with $pfdp$ and $B$ with the $P_{vector}$ and solve the system using LSQR.

At each iteration, $A$ and $B$ the resulting model is used to create a new vector $P_{vector}$ for an other instance of the loop.

At the end of the loop, plotting function are to visualise the solution. As we can observed after executing these lines, the values cluster around the centers of the Gaussians of the mixture model.

Calc_Jac.hide_toggle(for_next=True)

Toggle show/hide next cell

entropy = False
verbose = True
use_geology = False

pgrav_w = 7e-7 # Suggested value: 7e-7

weight = Calc_Jac.set_mod(filename_geol_grav)

mod_pre = m_inv[0]
diff_mod = 0
for i in range(line * column):
    diff_mod += mp[i] - mod_pre[i]
diff_mod = diff_mod ** 2

stop = 0

while diff_mod >= 1e-3 and stop != 20:


    P_vector = []

    for k in range(line * column):
        L = []
        for l in range(nb_gauss):
            
            if use_geology:
                L.append(pdfs[l][k] * weight[k][l])
            else:
                L.append(pdfs[l][k])
                
        max_pdf = L.index(numpy.amax(L))
        P_vector.append(Cm_grav[k][k] * pgrav_w * means_grav[max_pdf])

    B = numpy.concatenate((B, P_vector))

    pdfp = numpy.eye(line * column, line * column)

    pdfp = numpy.dot(pdfp,pgrav_w)

    A = numpy.concatenate((A, pdfp))



    m_inv = lsqr(A, B, damp=0.0, atol=1e-08, btol=1e-08, conlim=1e08, iter_lim=100, x0=xs)

    if not entropy:

        A = numpy.split(A, [column, line*column+column])
        A = numpy.concatenate((A[0], A[1]))

        B= numpy.split(B, [column, line*column+column])
        B= numpy.concatenate((B[0], B[1]))

    for h in range(line * column):
        m_inv[0][h] /= Cm_grav[h][h]

    rm_inv = numpy.reshape(m_inv[0], (line, column))

    xs = m_inv[0]

    nb_gauss, pdfs = PC.PDF(xs, means_grav, standard_dev_grav)

    if verbose:
        
        misfitgrav = Calc_Jac.model_misfitgrav(m_inv[0], filename_grav)
        props = dict(boxstyle='round', facecolor='white', alpha=0.9)
        
        plt.figure(figsize=(15,10))
        plt.imshow(rm_inv, cmap='jet')
        cbar3 = plt.colorbar(shrink=0.5)
        cbar3.set_label('$kg.m^{-3}$', labelpad=-25, y=1.15, rotation=0)

        if entropy:
            plt.title('Density contrast model with petrophysical contraints and structure term')
        else:
            plt.title('Density contrast model with petrophysical contraints')
        plt.xlabel('Horizontal position')
        plt.ylabel('Depth')
        plt.text(0.5, 2.5, r'Model misfit = ' + str(misfitgrav), bbox=props)
        plt.draw()
        plt.pause(0.0001)
        plt.clf()

    diff_mod = 0
    for i in range(line * column):
        diff_mod += m_inv[0][i] - mod_pre[i]
    diff_mod = diff_mod ** 2

    mod_pre = m_inv[0]
    stop += 1

<Figure size 432x288 with 0 Axes>

<Figure size 432x288 with 0 Axes>

Petrophysical constraints using probabilistic geological modelling

When probabilistic geological modelling results are available, we can condition the petrophysical using the observation probability of the different rock types.
These probabilities are obtained using the Monte Carlo Uncertainty Estimation method (MCUE). In MCUE, a series of geological models honouring geological measurements is generated. These models are then used to calculate a probabilistic model comprised of the observation probability of the different lithologies in every cell of the model (Pakyuz-Charrier et al, 2018).

When available, such probabilities can be used to weight the different Gaussians of the mixture model representing the petrpohysical measurements in the area accordingly with the lithology they correspond to (Giraud et al., 2017, 2019). In such case, the petrophysical constraints applied vary across the model accordingly with probabilistic geological modelling results.

The probabilities used here are shown below.

        Figure 6: Observation probabilities of lithologies labelled 1 (basement) through 6 (most recent sediments)

Calc_Jac.hide_toggle(for_next=True)

Toggle show/hide next cell

pc=True
entropy = False
use_geology = False

grav_w = 1e-7 # Suggested value: 1e-7
pgrav_w = 7e-7 # Suggested value: 7e-7
w_entropy = 5e-7 # Suggested value: 5e-7

line, column = Calc_Jac.get_dim(filename_grav)
Cm_grav, Cm_mag = Calc_Jac.calc_DW(filename_grav)

pos_grav = Calc_Jac.posx(filename_grav,start=0.5, interval=1)

Gij, Dg = Calc_Jac.grav_jac(filename_grav,pos_grav)
Id_grav = numpy.eye(line*column,line*column)

### Zero model as prior model

#mp = numpy.zeros((line*column,))

### Vertical gradient as prior model

mp = numpy.reshape(Calc_Jac.calc_vertical_gradient(filename_grav),(line*column,))

for i in range(line * column):
    mp[i] *= Cm_grav[i][i]

for i in range(len(Gij)):
    for j in range(len(Gij[0])):
        Gij[i][j] /= Cm_grav[j][j]

A = numpy.concatenate((Gij,numpy.dot(Id_grav,grav_w)))
B = numpy.concatenate((Dg, numpy.dot(mp,grav_w)))

if entropy:
    gradx, grady = Joint_inversion.calc_grad(filename_grav, filename_geol_grav, filename_geol_mag, filename_geol_gravmag, mag = False,verbose=False)
    empty_vec = numpy.zeros((line*column,))

    A = numpy.concatenate((A, numpy.dot(gradx,w_entropy), numpy.dot(grady,w_entropy)))
    B = numpy.concatenate((B, numpy.dot(empty_vec,w_entropy), numpy.dot(empty_vec,w_entropy)))

if pc:

    means_grav = [0,100,200,300]
    standard_dev_grav = [40,40,40,40]

    m_inv = lsqr(A, B, damp=0.0, atol=1e-08, btol=1e-08, conlim=1e08, iter_lim=100)

    for h in range(line * column):
        m_inv[0][h] /= Cm_grav[h][h]


    xs = m_inv[0]
    nb_gauss, pdfs = PC.PDF(xs, means_grav, standard_dev_grav)

    mod_pre = m_inv[0]
    diff_mod = 0
    for i in range(line * column):
        diff_mod += mp[i] - mod_pre[i]
    diff_mod = diff_mod ** 2

    stop = 0

    while diff_mod >= 1e-3 and stop != 20:

        P_vector = []

        for k in range(line * column):
            L = []
            for l in range(nb_gauss):

                if use_geology:
                    L.append(pdfs[l][k] * weight[k][l])
                else:
                    L.append(pdfs[l][k])

            max_pdf = L.index(numpy.amax(L))
            P_vector.append(Cm_grav[k][k] * pgrav_w * means_grav[max_pdf])

        B = numpy.concatenate((B, P_vector))

        pdfp = numpy.eye(line * column, line * column)

        pdfp = numpy.dot(pdfp,pgrav_w)

        A = numpy.concatenate((A, pdfp))



        m_inv = lsqr(A, B, damp=0.0, atol=1e-08, btol=1e-08, conlim=1e08, iter_lim=100, x0=xs)

        if not entropy:

            A = numpy.split(A, [column, line*column+column])
            A = numpy.concatenate((A[0], A[1]))

            B= numpy.split(B, [column, line*column+column])
            B= numpy.concatenate((B[0], B[1]))

        else:
            A = numpy.split(A, [column, line*column*3+column])
            A = numpy.concatenate((A[0], A[1]))

            B = numpy.split(B, [column, line*column*3+column])
            B = numpy.concatenate((B[0], B[1]))

        for h in range(line * column):
            m_inv[0][h] /= Cm_grav[h][h]

        rm_inv = numpy.reshape(m_inv[0], (line, column))
        m_inv_grav_p = rm_inv # to keep it for potential use at next stage. 

        xs = m_inv[0]

        nb_gauss, pdfs = PC.PDF(xs, means_grav, standard_dev_grav)

        if verbose:
            
            misfitgrav = Calc_Jac.model_misfitgrav(m_inv[0], filename_grav)
            props = dict(boxstyle='round', facecolor='white', alpha=0.9)
            
            plt.figure(figsize=(15,10))
            plt.imshow(rm_inv, cmap='jet')
            cbar4 = plt.colorbar(shrink=0.5)
            cbar4.set_label('$kg.m^{-3}$', labelpad=-30, y=1.15, rotation=0)

            if entropy:
                plt.title('Density contrast model with petrophysical contraints and structure term')
            else:
                plt.title('Density contrast model with petrophysical contraints')

            plt.xlabel('Horizontal position')
            plt.ylabel('Depth')
            plt.text(0.5, 2.5, r'Model misfit = ' + str(misfitgrav), bbox=props)
            plt.draw()
            plt.pause(0.0001)
            plt.clf()

        diff_mod = 0
        for i in range(line * column):
            diff_mod += m_inv[0][i] - mod_pre[i]
        diff_mod = diff_mod ** 2

        mod_pre = m_inv[0]
        stop += 1

<Figure size 432x288 with 0 Axes>

<Figure size 432x288 with 0 Axes>

Exercise: What happens if you change entropy to True ?

4. MAGNETIC INVERSION

The inversion of magnetic data follows the same philosophy as the gravity data.
In terms of petrophysical properties, the centers of the gaussians are: 0, 0.025, 0.05, 0.075 SI.

Parameters of the magnetic inversion:

• filename_mag: the name of the file containing the magnetic data.

• mag_w: weight of the the magnetic data.

• pmag_w: The weights of the petrophysical constraints for the magnetic model.

• w_entropy: weight of the information entropy.

• pc: Boolean to add the petrophysical constraints to the system, set to False by default

• entropy: Boolean to add the entropy and the gradients of the model to the system, set to False by default

• Geol : Boolean the geological information of the model, set to False by default

Calc_Jac.hide_toggle(for_next=True)

Toggle show/hide next cell

pc = False
entropy = False
use_geology = False

mag_w = 1e-13 # Suggested value: 5e-13
pmag_w = 3e-13 # Suggested value: 5e-13
w_entropy = 0e-14 # Suggested value: 5e-14

means_mag = [0, 0.025, 0.05, 0.075]
standard_dev_mag = [0.02, 0.02, 0.02, 0.02]

cmap = plt.get_cmap('jet')
newcmap = Calc_Jac.truncate_colormap(cmap, 0.5, 0.95)

line, column = Calc_Jac.get_dim(filename_mag)

pos_mag = Calc_Jac.posx(filename_mag, start=0.5, interval = 1)
Db, Jij = Calc_Jac.mag_jac(filename_mag, pos_mag)

### Zero model as prior model

# mp = numpy.zeros((line*column,))

### Vertical gradient as prior model

mp = numpy.reshape(Calc_Jac.calc_vertical_gradient(filename_mag), (line * column,))

for i in range(line * column):
    mp[i] *= Cm_mag[i][i]


### Only with the vertical component
Kij = numpy.zeros((column, line * column))

Dm = numpy.zeros((column,))

for i in range(len(Jij)):
    for j in range(len(Jij[0])):
        Kij[i][j] = Jij[i][j][0]
        Kij[i][j] /= Cm_mag[j][j]

for e in range(len(Db)):
    Dm[e] = Db[e][0]

weight = Calc_Jac.set_mod(filename_geol_mag)
mean_mod = []


Id = numpy.eye(line * column, line * column)

A = numpy.concatenate((Kij, numpy.dot(Id, mag_w)))
B = numpy.concatenate((Dm, numpy.dot(mp,mag_w)))

if use_geology:
    
    probam = 0
    for i in range(line*column):
        for j in range(len(means_mag)):
            probam += weight[i][j] * means_mag[j]
        mean_mod.append(probam)
        probam = 0
    
    for i in range(line*column):
        mean_mod[i]*= Cm_mag[i][i]
    
    B = numpy.concatenate((Dm, numpy.dot(mean_mod,mag_w)))


if entropy:
    gradx, grady = Joint_inversion.calc_grad(filename_grav, filename_geol_grav, filename_geol_mag, filename_geol_gravmag, mag = True,verbose=False)
    empty_vec = numpy.zeros((line*column,))

    A = numpy.concatenate((A, numpy.dot(gradx,w_entropy), numpy.dot(grady,w_entropy)))
    B = numpy.concatenate((B, numpy.dot(empty_vec,w_entropy), numpy.dot(empty_vec,w_entropy)))

if pc:

    m_inv = lsqr(A, B, damp=0.0, atol=1e-20, btol=1e-20, conlim=1e20, iter_lim=1000)

    for h in range(line * column):
        m_inv[0][h] /= Cm_mag[h][h]

    xs = m_inv[0]
    nb_gauss, pdfs = PC.PDF(xs,means_mag, standard_dev_mag)



    mod_pre = m_inv[0]
    diff_mod = 0
    for i in range(line * column):
        diff_mod += mp[i] - mod_pre[i]
    diff_mod = diff_mod**2

    e = 0

    while diff_mod>=1e-5 and e==0:

        P_vector = []

        for k in range(line * column):
            L = []
            for l in range(nb_gauss):

                if use_geology:
                    L.append(pdfs[l][k] * weight[k][l])
                else:
                    L.append(pdfs[l][k])

            max_pdf = L.index(numpy.amax(L))
            P_vector.append(Cm_mag[k][k] * pmag_w * means_mag[max_pdf])

        B = numpy.concatenate((B, P_vector))

        pdfp = numpy.eye(line * column, line * column)
        pdfp = numpy.dot(pdfp, pmag_w)

        A = numpy.concatenate((A, pdfp))

        m_inv = lsqr(A, B, damp=0.0, atol=1e-20, btol=1e-20, conlim=1e20, iter_lim=1000)


        if not entropy:

            A = numpy.split(A, [column, line*column+column])
            A = numpy.concatenate((A[0], A[1]))

            B = numpy.split(B, [column, line*column+column])
            B = numpy.concatenate((B[0], B[1]))


        else:
            A = numpy.split(A, [column, line*column*3+column])
            A = numpy.concatenate((A[0], A[1]))

            B = numpy.split(B, [column, line*column*3+column])
            B = numpy.concatenate((B[0], B[1]))

        for h in range(line * column):
            m_inv[0][h] /= Cm_mag[h][h]

        rm_inv = numpy.reshape(m_inv[0], (line, column))

        xs = m_inv[0]
        nb_gauss, pdfs = PC.PDF(xs,means_mag, standard_dev_mag)

        if verbose:
            plt.figure(figsize=(15,10))
            plt.imshow(rm_inv, newcmap)
            cbar5 = plt.colorbar(shrink=0.5)
            cbar5.set_label('SI', labelpad=-50, y=1.15, rotation=0)
            if entropy:
                plt.title('Magnetic susc. model with petrophysical contraints and structure term')
            else:
                plt.title('Magnetic susc. model with petrophysical contraints')
            
            plt.xlabel('Horizontal position')
            plt.ylabel('Depth')
            plt.draw()
            plt.pause(0.0001)
            plt.clf()

        diff_mod = 0
        for i in range(line*column):
            diff_mod += m_inv[0][i] - mod_pre[i]
        diff_mod = diff_mod ** 2

        mod_pre = m_inv[0]
        e += 1

if not pc:

    m_inv = lsqr(A, B, damp=0.0, atol=1e-20, btol=1e-20, conlim=1e20, iter_lim=1000)
    for h in range(line*column):
        m_inv[0][h] /= Cm_mag[h][h]


rm_inv = numpy.reshape(m_inv[0], (line, column))
m_inv_mag_p=rm_inv; # to keep it for potential use at next stage. 

if verbose:
    misfitmag = Calc_Jac.model_misfitmag(m_inv[0], filename_mag)
    props = dict(boxstyle='round', facecolor='white', alpha=0.9)
    
    if pc and not entropy:
        plt.figure(figsize=(15,10))
        plt.imshow(rm_inv, newcmap)
        cbar6 = plt.colorbar(shrink=0.5)
        cbar6.set_label('SI', labelpad=-50, y=1.15, rotation=0)
        plt.title('Magnetic susc. model with petrophysical contraints')
        plt.xlabel('Horizontal position')
        plt.ylabel('Depth')
        plt.text(0.5, 2.5, r'Model misfit = ' + str(misfitmag), bbox=props)
        plt.show()

    if entropy and pc:
        plt.figure(figsize=(15,10))
        plt.imshow(rm_inv, newcmap)
        cbar7 = plt.colorbar(shrink=0.5)
        cbar7.set_label('SI', labelpad=-50, y=1.15, rotation=0)
        plt.title('Magnetic susc. model with petrophysical contraints and structure term')
        plt.xlabel('Horizontal position')
        plt.ylabel('Depth')
        plt.text(0.5, 2.5, r'Model misfit = ' + str(misfitmag), bbox=props)
        plt.show()

    if entropy and not pc:
        plt.figure(figsize=(15,10))
        plt.imshow(rm_inv, newcmap)
        cbar8 = plt.colorbar(shrink=0.5)
        cbar8.set_label('SI', labelpad=-50, y=1.15, rotation=0)
        plt.title('Magnetic susc. model with structure term')
        plt.xlabel('Horizontal position')
        plt.ylabel('Depth')
        plt.text(0.5, 2.5, r'Model misfit = ' + str(misfitmag), bbox=props)
        plt.show()

    if not pc and not entropy:
        plt.figure(figsize=(15,10))
        plt.imshow(rm_inv, newcmap)
        cbar9 = plt.colorbar(shrink=0.5)
        cbar9.set_label('SI', labelpad=-50, y=1.15, rotation=0)
        plt.title('Unconstrained magnetic susc. model')
        plt.xlabel('Horizontal position')
        plt.ylabel('Dephth')
        plt.text(0.5, 2.5, r'Model misfit = ' + str(misfitmag), bbox=props)
        plt.show()

for i in range(len(Jij)):
    for j in range(len(Jij[0])):
        Kij[i][j] *= Cm_mag[j][j]

mul = numpy.matmul(Kij, m_inv[0])

rms = Calc_Jac.calc_RMS(Kij,Dm, m_inv[0])

plt.figure(figsize=(15,10))
plt.title('Magnetic anomaly')
line_up, = plt.plot(mul, label='Line 1')
line_down, = plt.plot(Dm, label='Line 2')
plt.legend([line_up, line_down], ['Calculated data', 'Measured data'])
plt.xlabel('Horizontal position')
plt.ylabel('Magnetic anomaly (nT)')
plt.grid(b=True, which='major', color='0.6', linestyle='-')
plt.minorticks_on()
plt.grid(b=True, which='minor', color='0.95', linestyle='--')
props = dict(boxstyle='round', facecolor='white', alpha=0.5)
plt.text(20, 0.25e-13, r'rms = ' + str(rms), bbox=props)
plt.show()

Exercise: What happens if you change :

• 1: geol to True,

• 2: entropy to True

• 3: geol and entropy to True

5. JOINT INVERSION

Now that we have performed single-domain inversions using petrophysical constraints, we can combine them to perform joint inversion.
The method is very similar to single-domain inversion. The only difference is that we compute the Gaussian mixture model in cross-plot domain (see Figure 5), via the python function multivariate_normal.pdf form the scipy.stats library. This will retrun a $P_{vector}$ (P_vector) that will be used in the LSQR.

More information on the multivariate_nomal.pdf function is available on this link : http://lagrange.univ-lyon1.fr/docs/scipy/0.17.1/generated/scipy.stats.multivariate_normal.html

the code below offers the possibility to perform joint inversion using global petrophysical constraints (ie, no geological information), petrophysical constraints conditioned by probabilistic geological modelling, and structural constraints using geological uncertainty information.

Note: not all parameters or combinations of constraints may work well at the first attempt.

Joint_inversion parameters:

• filename_grav and filename_mag: the name of the files containing the gravimetic and magnetic data

• grav_w and mag_w: weights of the gravimetric and magnetic data, respectively.

• pgrav_w and pmag_w: The weights of the petrophysical constraints for the gravimetric and magnetic models

• w_entropy: weight for the structural term

• entropy: Boolean to add the entropy and the gradients of the model to the system, set to False by default

• verbose: Boolean to view or not the solution, set to True by default

• Geol : Boolean the geological information of the model, set to False by default

Calc_Jac.hide_toggle(for_next=True)

Toggle show/hide next cell

entropy = False
use_geology = True
verbose = True 

grav_w = 5e-9 # Suggested value: 1e-8 ; 2e-7 
mag_w = 5e-13 # Suggested value: 1e-14 ; 1e-13
pgrav_w = 2e-8 # Suggested value: 5e-8; 1e-8
pmag_w = 9e-13 # Suggested value: 9e-13
w_entropygrav = 0e-7 # Suggested value: 5e-7
w_entropymag = 0e-12 # Suggested value: 5e-12

means_onlygrav = [0, 100, 200, 300]
means_onlymag = [0, 0.025, 0.05,0.075]

means_grav = [0, 0, 100, 200, 300, 300]
means_mag = [0, 0.025, 0.05, 0.025, 0.05, 0.075]
standard_dev_grav = [40, 40, 40, 40, 40, 40]
standard_dev_mag = [0.0000135, 0.0000135, 0.0000135, 0.0000135, 0.0000135, 0.0000135]

weight_grav = Calc_Jac.set_mod(filename_geol_grav)
weight_mag = Calc_Jac.set_mod(filename_geol_mag)

cmap = plt.get_cmap('jet')
newcmap = Calc_Jac.truncate_colormap(cmap, 0.5, 0.95)

line, column = Calc_Jac.get_dim(filename_grav)
Cm_grav, Cm_mag = Calc_Jac.calc_DW(filename_grav)

pos_grav = Calc_Jac.posx(filename_grav, start=0.5, interval=1)

Gij, Dg = Calc_Jac.grav_jac(filename_grav, pos_grav)
Id_grav = numpy.eye(line * column, line * column)

### Zero model as prior model ###
# mpg = numpy.zeros((line*column,))
# mpm = numpy.zeros((line*column,))

### Vertical gradient as prior model ###
# mpg = numpy.reshape(Calc_Jac.calc_vertical_gradient(filename_grav), (line * column,))
# mpm = numpy.reshape(Calc_Jac.calc_vertical_gradient(filename_mag), (line * column,))

### Use separate petrophysical inversion results as prior model ### - careful, Jupyter might have a problem with this. 
mpg = numpy.reshape(m_inv_grav_p, (line * column,))
mpm = numpy.reshape(m_inv_mag_p, (line * column,))

for i in range(line * column):
    mpg[i] *= Cm_grav[i][i]
    mpm[i] *= Cm_mag[i][i]

for i in range(len(Gij)):
    for j in range(len(Gij[0])):
        Gij[i][j] /= Cm_grav[j][j]

A_grav = numpy.concatenate((Gij, numpy.dot(Id_grav, grav_w)))
B_grav = numpy.concatenate((Dg, numpy.dot(mpg, grav_w)))

pos_mag = Calc_Jac.posx(filename_mag, start=0.5, interval=1)
Db, Jij = Calc_Jac.mag_jac(filename_mag, pos_mag)

### Only with the vertical component
Kij = numpy.zeros((column, line * column))

Dm = numpy.zeros((column,))

for i in range(len(Jij)):
    for j in range(len(Jij[0])):
        Kij[i][j] = Jij[i][j][0]
        Kij[i][j] /= Cm_mag[j][j]

for e in range(len(Db)):
    Dm[e] = Db[e][0]

Id = numpy.eye(line * column, line * column)

A_mag = numpy.concatenate((Kij, numpy.dot(Id, mag_w)))
B_mag = numpy.concatenate((Dm, numpy.dot(mpm, mag_w)))


if use_geology:
    mean_modmag = []
    mean_modgrav = []
    probam = 0
    probag = 0
    for i in range(line*column):
        for j in range(len(means_onlymag)):
            probam += weight_mag[i][j] * means_onlymag[j]
        mean_modmag.append(probam)
        probam = 0
        
        
        for k in range(len(means_onlygrav)):
            probag += weight_grav[i][k] * means_onlygrav[k]
        mean_modgrav.append(probag)
        probag = 0

    for i in range(line*column):
        mean_modmag[i] *= Cm_mag[i][i]
        mean_modgrav[i] *= Cm_grav[i][i]
    B_grav = numpy.concatenate((Dg, numpy.dot(mean_modgrav, grav_w)))
    B_mag = numpy.concatenate((Dm, numpy.dot(mean_modmag,mag_w)))




if entropy:
    gradx, grady = Joint_inversion.calc_grad(filename_grav, filename_geol)
    empty_vec = numpy.zeros((line*column,))
    A_grav = numpy.concatenate((A_grav, numpy.dot(gradx,w_entropygrav), numpy.dot(grady,w_entropygrav)))
    B_grav = numpy.concatenate((B_grav, numpy.dot(empty_vec,w_entropygrav), numpy.dot(empty_vec,w_entropygrav)))

    gradxm, gradym = Joint_inversion.calc_grad(filename_mag, filename_geol)
    A_mag = numpy.concatenate((A_mag, numpy.dot(gradxm,w_entropymag), numpy.dot(gradym,w_entropymag)))
    B_mag = numpy.concatenate((B_mag, numpy.dot(empty_vec,w_entropymag), numpy.dot(empty_vec,w_entropymag)))


m_inv_grav = lsqr(A_grav, B_grav, damp=0.0, atol=1e-08, btol=1e-08, conlim=1e08, iter_lim=100)
m_inv_mag = lsqr(A_mag, B_mag, damp=0.0, atol=1e-20, btol=1e-20, conlim=1e20, iter_lim=1000)


for h in range(line * column):
    m_inv_mag[0][h] /= Cm_mag[h][h]
    m_inv_grav[0][h] /= Cm_grav[h][h]

xs = m_inv_grav[0]
ys = m_inv_mag[0]

nb_gauss, pdfs = PC.PDF_2d(xs,ys, means_grav, means_mag, standard_dev_grav, standard_dev_mag)

weight_gravmag = Calc_Jac.set_mod(filename_geol_gravmag)

mod_pre = m_inv_grav[0]
diff_mod = 0
for i in range(line * column):
    diff_mod += mp[i] - mod_pre[i]
diff_mod = diff_mod ** 2

e = 0

while diff_mod >= 1e-3 and e != 20:

    G_vector = []
    M_vector = []
    Ind_vector = []

    for k in range(line * column):
        L = []
        for l in range(nb_gauss):
            if use_geology:
                L.append(pdfs[l][k] * weight_gravmag[k][l])
            else:
                L.append(pdfs[l][k])
        max_pdf = L.index(numpy.amax(L))
        Ind_vector.append(max_pdf)
        
        G_vector.append(Cm_grav[k][k] * pgrav_w * means_grav[max_pdf])
        M_vector.append(Cm_mag[k][k] * pmag_w * means_mag[max_pdf])
        
    B_grav = numpy.concatenate((B_grav, G_vector))
    B_mag = numpy.concatenate((B_mag, M_vector))

    pdfp_g = numpy.eye(line * column,line * column)
    pdfp_g = numpy.dot(pdfp_g, pgrav_w)
    A_grav = numpy.concatenate((A_grav, pdfp_g))

    pdfp_m = numpy.eye(line * column, line * column)
    pdfp_m = numpy.dot(pdfp_m, pmag_w)
    A_mag = numpy.concatenate((A_mag, pdfp_m))

    m_inv_grav = lsqr(A_grav, B_grav, damp=0.0, atol=1e-08, btol=1e-08, conlim=1e08, iter_lim=100)
    m_inv_mag = lsqr(A_mag, B_mag, damp=0.0, atol=1e-20, btol=1e-20, conlim=1e20, iter_lim=1000)


    if not entropy:

        A_grav = numpy.split(A_grav, [column,line*column+column])
        A_grav = numpy.concatenate((A_grav[0], A_grav[1]))
        A_mag = numpy.split(A_mag, [column, line*column+column])
        A_mag = numpy.concatenate((A_mag[0], A_mag[1]))

        B_grav = numpy.split(B_grav, [column, line*column+column])
        B_grav = numpy.concatenate((B_grav[0], B_grav[1]))
        B_mag = numpy.split(B_mag, [column, line*column+column])
        B_mag = numpy.concatenate((B_mag[0], B_mag[1]))

    else:
        A_grav = numpy.split(A_grav, [column, line*column*3+column])
        A_grav = numpy.concatenate((A_grav[0], A_grav[1]))
        A_mag = numpy.split(A_mag, [column, line*column*3+column])
        A_mag = numpy.concatenate((A_mag[0], A_mag[1]))


        B_grav = numpy.split(B_grav, [column, line*column*3+column])
        B_grav = numpy.concatenate((B_grav[0], B_grav[1]))
        B_mag = numpy.split(B_mag, [column, line*column*3+column])
        B_mag = numpy.concatenate((B_mag[0],B_mag[1]))

    for h in range(line * column):
        m_inv_grav[0][h] /= Cm_grav[h][h]
        m_inv_mag[0][h] /= Cm_mag[h][h]

    xs = m_inv_grav[0]
    ys = m_inv_mag[0]
    nb_gauss, pdfs = PC.PDF_2d(xs,ys, means_grav, means_mag, standard_dev_grav, standard_dev_mag)

    reTeta_grav = numpy.reshape(m_inv_grav[0], (line, column))
    rm_invmag = numpy.reshape(m_inv_mag[0], (line, column))

    if verbose:
        plt.figure(figsize=(15,10))
        plt.subplot(221)
        plt.imshow(reTeta_grav, cmap='jet')
        cbar10 = plt.colorbar(shrink=0.5)
        cbar10.set_label('$kg.m^{-3}$', labelpad=-30, y=1.15, rotation=0)
        plt.title('iteration number: '+ str(e+1))
        plt.xlabel('Horizontal position')
        plt.ylabel('Depth')

        plt.subplot(222)
        plt.imshow(rm_invmag, newcmap)
        cbar11 = plt.colorbar(shrink=0.5)
        cbar11.set_label('$SI$', labelpad=-40, y=1.15, rotation=0)
        plt.xlabel('Horizontal position')
        plt.ylabel('Depth')
        
        plt.subplot(212)
        plt.title('Cross plot of the value od the Magnetic and Density contrast inversion')
        plt.scatter(xs, ys, c=Ind_vector, alpha=0.5)
        plt.plot(means_grav, means_mag, 'ks')
        plt.xlabel('Density constrast value')
        plt.ylabel('Magnetic susc. value')
        plt.draw()
        plt.pause(0.0001)
        plt.clf()
        
    diff_mod = 0

    for i in range(line * column):
        diff_mod += m_inv_grav[0][i] - mod_pre[i]
    diff_mod = diff_mod ** 2

    mod_pre = m_inv_grav[0]
    e += 1

if verbose:
    misfitgrav = Calc_Jac.model_misfitgrav(m_inv_grav[0], filename_grav)
    misfitmag = Calc_Jac.model_misfitmag(m_inv_mag[0], filename_mag)
    props = dict(boxstyle='round', facecolor='white', alpha=0.9)
    
    plt.figure(figsize=(15,10))
    plt.subplot(221)
    plt.imshow(reTeta_grav, cmap='jet')
    cbar12 = plt.colorbar(shrink=0.5)
    cbar12.set_label('$kg.m^{-3}$', labelpad=-30, y=1.15, rotation=0)
    plt.title('Density constraint model - joint inversion')
    plt.xlabel('Horizontal position')
    plt.ylabel('Depth')
    plt.text(2, 5, r'Model misfit = ' + str(misfitgrav), bbox=props)

    plt.subplot(222)
    plt.imshow(rm_invmag, newcmap)
    cbar13 = plt.colorbar()
    cbar13.set_label('$SI$', labelpad=-40, y=1.15, rotation=0)
    plt.title('Magnetic susc. model - joint inversion')
    plt.xlabel('Horizontal position')
    plt.ylabel('Depth')
    plt.text(2, 5, r'Model misfit = ' + str(misfitmag), bbox=props)

    plt.subplot(212)
    plt.title('Cross plot of the value od the Magnetic and Density contrast inversion')
    plt.scatter(xs, ys, c=Ind_vector, alpha=0.5)
    plt.plot(means_grav, means_mag, 'ks')
    plt.xlabel('Density constrast value')
    plt.ylabel('Magnetic susc. value')
    plt.draw()
    plt.pause(0.0001)
    plt.clf()
    
    for i in range(len(Kij)):
        for j in range(len(Kij[0])):
            Kij[i][j] *= Cm_mag[j][j]
            Gij[i][j] *= Cm_grav[j][j]

    mulm = numpy.matmul(Kij, m_inv_mag[0])

    mulg = numpy.matmul(Gij, m_inv_grav[0])
    
    rmsmag = Calc_Jac.calc_RMS(Kij,Dm,  m_inv_mag[0])

    rmsgrav = Calc_Jac.calc_RMS(Gij,Dg, m_inv_grav[0])

    plt.figure(figsize=(15,10))
    plt.subplot(211)
    plt.title('Gravity anomaly')
    line_up, = plt.plot(mulg, label='Line 1')
    line_down, = plt.plot(Dg, label='Line 2')
    plt.legend([line_up, line_down], ['Calculated data', 'Measured data'])
    plt.xlabel('Horizontal position')
    plt.ylabel('Gravity anomaly (gal)')
    plt.grid(b=True, which='major', color='0.6', linestyle='-')
    plt.minorticks_on()
    plt.grid(b=True, which='minor', color='0.95', linestyle='--')
    props = dict(boxstyle='round', facecolor='white', alpha=0.5)
    plt.text(20, 0.00215, r'rms grav = ' + str(rmsgrav), bbox=props)

    plt.subplot(212)
    plt.title('Magnetic anomaly')
    line_up, = plt.plot(mulm, label='Line 1')
    line_down, = plt.plot(Dm, label='Line 2')
    plt.legend([line_up, line_down], ['Calculated data', 'Measured data'])
    plt.xlabel('Horizontal position')
    plt.ylabel('Magnetic anomaly (nT)')
    plt.grid(b=True, which='major', color='0.6', linestyle='-')
    plt.minorticks_on()
    plt.grid(b=True, which='minor', color='0.95', linestyle='--')
    props = dict(boxstyle='round', facecolor='white', alpha=0.5)
    plt.text(20, 0.1e-13, r'rms mag = ' + str(rmsmag), bbox=props)
    plt.show()
    
    del(mulg, mulm, rmsmag, rmsgrav)
    
del(Gij, Kij, Jij, A_mag, A_grav, B_grav, B_mag, m_inv_grav, m_inv_mag, Dg, Dm, Db, mpg, mpm, L)

<Figure size 432x288 with 0 Axes>

<Figure size 432x288 with 0 Axes>

<Figure size 432x288 with 0 Axes>

Exercise: What happens if you change :

• 1: geol to True,

• 2: entropy to True

• 3: geol and entropy to True

• 4: Prior model? How do you get a more geologically sensible model?

$REFERENCES$

• Boulanger, O. and Chouteau, M., 2001, Constraints in 3D gravity inversion, Geophys. Prospect., 49(2), 265–280, doi:10.1046/j.1365-2478.2001.00254.x, 2001.

• Giraud, J., E. Pakyuz-Charrier, M. Jessell, M. Lindsay, R. Martin, and V. Ogarko, 2017, Uncertainty reduction through geologically conditioned petrophysical constraints in joint inversion: GEOPHYSICS, 82, ID19-ID34.

• Giraud, J., Lindsay, M., Ogarko, V., Jessell, M., Martin, R. & Pakyuz- Charrier, E., 2019. Integration of geological uncertainty into geophysical inversion by means of local gradient regularization, Solid Earth, 10, 193–210.

• Giraud, J., Ogarko, V., Lindsay, M., Pakyuz-Charrier, E., Jessell, M., and Martin, R., 2019, Sensitivity of constrained joint inversions to geological and petrophysical input data uncertainties with posterior geological analysis, Geophysical Journal International, Volume 218, Issue 1, July 2019, Pages 666–688, https://doi.org/10.1093/gji/ggz152

• Li, Y. & Chouteau, M., 1998. Three-dimentional gravity modeling in all space: Survesy in GEOPHYSICS, 19, 10.1023/A:1006554408567.

• Li, Y. & Oldenburg D., 1996. 3-D inversion of magnetic data. GEOPHYSICS, 61, 10.1190/1.1443968.

• Li, Y. & Oldenburg D.,1998. 3-D inversion of gravity data. GEOPHYSICS, 63, 10.1190/1.1444302.

• Martin, R., Monteiller, V., Komatitsch, D., Perrouty, S., Jessell, M., Bonvalot, S. & Lindsay, M.D., 2013. Gravity inversion using wavelet-based compression on parallel hybrid CPU/GPU systems: application to southwest Ghana, Geophys. J. Int., 195, 1594–1619.

• Martin, R., Ogarko, V., Komatitsch, D. & Jessell, M., 2018. Parallel threedimensional electrical capacitance data imaging using a nonlinear inversion algorithm and Lp norm-based model, Measurement, 128, 428–445.

• Paige, C.C. & Saunders, M.A., 1982. LSQR: an algorithm for sparse linear equations and sparse least squares, ACM Trans. Math. Softw., 8, 43–71.

• Pakyuz-Charrier, E., & Intrepid Geophysics. (2017). Mansfield (Victoria, Australia) area original GeoModeller model and relevant MCUE outputs [Data set]. Zenodo. http://doi.org/10.5281/zenodo.848225

• Pakyuz-Charrier, E., Lindsay, M., Ogarko, V., Giraud, J., and Jessell, M.: Monte Carlo simulation for uncertainty estimation on structural data in implicit 3-D geological modeling, a guide for disturbance distribution selection and parameterization, Solid Earth, 9, 385-402, https://doi.org/10.5194/se-9-385-2018, 2018.

• Schweizer, D., Blum, P. and Butscher, C.: Uncertainty assessment in 3-D geological models of increasing complexity, Solid 10 Earth, 8(2), 515–530, doi:10.5194/se-8-515-2017, 2017.

• Shannon, C. E. E.: A Mathematical Theory of Communication, Bell Syst. Tech. J., 27(3), 379–423, doi:10.1002/j.1538- 7305.1948.tb01338.x, 1948.

• Wellmann, J. F. and Regenauer-Lieb, K.: Uncertainties have a meaning: Information entropy as a quality measure for 3-D geological models, Tectonophysics, 526–529, 207–216, doi:10.1016/j.tecto.2011.05.001, 2012.