TKC Magnetics

Where are the diamonds: Using Earth’s Potential Fields (Magnetics)

In this tutorial we will run a magnetic inverse problem from an input file. The synthetic example used is based on the TKC kimberlite deposit, Northwest Territories. We are attempting to use the magnetic field to image the various rock units making up the diamond deposit. We will skip the model creation for now, at it is quite involved and out of scope for this notebook.

In [1]:
## First we need to load all the libraries and set up the path
## for the input files. Same files as used by the online tutorial
import matplotlib.pyplot as plt
from SimPEG import Mesh
from SimPEG import Utils
from SimPEG import Maps
from SimPEG import Regularization
from SimPEG import DataMisfit
from SimPEG import Optimization
from SimPEG import InvProblem
from SimPEG import Directives
from SimPEG import Inversion
from SimPEG import PF
from SimPEG.Utils.io_utils import download
import matplotlib
import matplotlib.colors as colors
import scipy as sp
import numpy as np
import time as tm
import os
import shutil
from ipywidgets.widgets import interact, IntSlider
%matplotlib inline
In [2]:
psep = os.path.sep
url = 'https://storage.googleapis.com/simpeg/tkc_synthetic/potential_fields/'
cloudfiles = ['MagData.obs', 'Mesh.msh',
              'Initm.sus', 'SimPEG_PF_Input.inp']

downloads = download([url+f for f in cloudfiles], folder='./MagTKC/', overwrite=True)
downloads = dict(zip(cloudfiles, downloads))
input_file = downloads['SimPEG_PF_Input.inp']
overwriting MagData.obs
overwriting Mesh.msh
overwriting Initm.sus
overwriting SimPEG_PF_Input.inp
Downloading https://storage.googleapis.com/simpeg/tkc_synthetic/potential_fields/MagData.obs
Downloading https://storage.googleapis.com/simpeg/tkc_synthetic/potential_fields/Mesh.msh
Downloading https://storage.googleapis.com/simpeg/tkc_synthetic/potential_fields/Initm.sus
Downloading https://storage.googleapis.com/simpeg/tkc_synthetic/potential_fields/SimPEG_PF_Input.inp
Download completed!
In [3]:
# Read in the input file which included all parameters at once (mesh, topo, model, survey, inv param, etc.)
driver = PF.MagneticsDriver.MagneticsDriver_Inv(input_file)
In [4]:
# We already have created a mesh, model and survey for this example.
# All the elements are stored in the driver, and can be accessed like this:
mesh = driver.mesh
survey = driver.survey

Setup

We are using the integral form of the magnetostatic problem. In the absence of free-currents or changing magnetic field, magnetic material can give rise to a secondary magnetic field according to:

\[\vec b = \frac{\mu_0}{4\pi} \int_{V} \vec M \cdot \nabla \nabla \left(\frac{1}{r}\right) \; dV\]

Where \(\mu_0\) is the magnetic permealitity of free-space, \(\vec M\) is the magnetization per unit volume and \(r\) defines the distance between the observed field \(\vec b\) and the magnetized object. Assuming a purely induced response, the strenght of magnetization can be written as:

\[\vec M = \mu_0 \kappa \vec H_0\]

where \(\vec H\) is an external inducing magnetic field, and \(\kappa\) the magnetic susceptibility of matter. As derived by Sharma 1966, the integral can be evaluated for rectangular prisms such that:

\[\vec b(P) = \mathbf{T} \cdot \vec H_0 \; \kappa\]

Where the tensor matrix \(\bf{T}\) relates the three components of magnetization \(\vec M\) to the components of the field \(\vec b\):

\[\begin{split}\mathbf{T} = \begin{pmatrix} T_{xx} & T_{xy} & T_{xz} \\ T_{yx} & T_{yy} & T_{yz} \\ T_{zx} & T_{zy} & T_{zz} \end{pmatrix}\end{split}\]

In general, we discretize the earth into a collection of cells, each contributing to the magnetic data such that:

\[\vec b(P) = \sum_{j=1}^{nc} \mathbf{T}_j \cdot \vec H_0 \; \kappa_j\]

giving rise to a linear problem.

In [5]:
# We did not include a topography file in this example as the information
# about inactive cells is already captured in our starting model.
# Line 6 of the input file specifies a VALUE to be used as inactive flag.

# Get the active cells
actv = driver.activeCells
nC = len(actv) # Number of active cells
ndv = -100

# Create active map to return to full space after the inversion
actvMap = Maps.InjectActiveCells(mesh, actv, ndv)

# Create a reduced identity map for the inverse problem
idenMap = Maps.IdentityMap(nP=nC)

Forward system:

Now that we have all our spatial components, we can create our linear system. For a single location and single component of the data, the system would looks like this:

\[\begin{split}b_x = \begin{bmatrix} T_{xx}^1 &... &T_{xx}^{nc} & T_{xy}^1 & ... & T_{xy}^{nc} & T_{xz}^1 & ... & T_{xz}^{nc}\\ \end{bmatrix} \begin{bmatrix} \mathbf{M}_x \\ \mathbf{M}_y \\ \mathbf{M}_z \end{bmatrix} \\\end{split}\]

where each of \(T_{xx},\;T_{xy},\;T_{xz}\) are [nc x 1] long. For the \(y\) and \(z\) component, we need the two other rows of the tensor \(\mathbf{T}\). In our simple induced case, the magnetization direction \(\mathbf{M_x,\;M_y\;,Mz}\) are known and assumed to be constant everywhere, so we can reduce the size of the system such that:

\[\vec{\mathbf{d}}_{\text{pred}} = (\mathbf{T\cdot M})\; \kappa\]

In most geophysical surveys, we are not collecting all three components, but rather the magnitude of the field, or \(Total\;Magnetic\;Intensity\) (TMI) data. Because the inducing field is really large, we will assume that the anomalous fields are parallel to \(H_0\):

\[d^{TMI} = \hat H_0 \cdot \vec d\]

We then end up with a much smaller system:

\[d^{TMI} = \mathbf{F\; \kappa}\]

where \(\mathbf{F} \in \mathbb{R}^{nd \times nc}\) is our \(forward\) operator.

In [6]:
# Now that we have a model and a survey we can build the linear system ...
# Create the forward model operator (the argument forwardOnly=False store the forward matrix to memory)
prob = PF.Magnetics.MagneticIntegral(mesh, chiMap=idenMap, actInd=actv, forwardOnly=False, rtype = 'tmi')

# Pair the survey and problem
survey.pair(prob)
In [7]:
# Fist time that we ask for predicted data,
# the dense matrix T is calculated.
# This is generally the bottleneck of the integral formulation in terms of cost
d = prob.fields(driver.m0)

# Add noise to the data and assign uncertainties
data = d + np.random.randn(len(d)) # We add some random Gaussian noise (1nT)
wd = np.ones(len(data))*1. # Assign flat uncertainties

survey.dobs = data
survey.std = wd

# [OPTIONAL] You can write the observations to UBC format here
#PF.Magnetics.writeUBCobs('MAG_Synthetic_data.obs',survey,data)
fig = PF.Magnetics.plot_obs_2D(survey.srcField.rxList[0].locs,d=data ,title='Mag Obs')
Begin calculation of forward operator: ind
Done 0.0 %
Done 10.0 %
Done 20.0 %
Done 30.0 %
Done 40.0 %
Done 50.0 %
Done 60.0 %
Done 70.0 %
Done 80.0 %
Done 90.0 %
Done 100% ...forward operator completed!!

../../_images/case-studies_PF_TKC_Mag_9_1.png

Inverse problem

We have generated synthetic data, we now what to see if we can solve the inverse. Using the usual formulation, we seek a model that can reproduce the data, let’s say a least-squares measure of the form:

The inverse problem is hard because we don’t have great data coverage, and the Earth is big, and there is usually noise in the data. So we need to add something to regularize it. The simplest way to do it is to penalize solutions that won’t make sense geologically, for example to assume that the model is small. The usual smooth inversion function use an l2-norm measure:

The full objective function to be minimized can be written as:

which will yield our usual small and smooth models.

We propose a fancier regularization function that can allow to recover sparse and blocky solutions. Starting with the well known Ekblom norm:

where \(x_i\) denotes some function of the model parameter, and \(\epsilon\) is a small value to avoid singularity as \(m\rightarrow0\). For p=2, we get the usual least-squares measure and we recover the regularization presented above. For \(p \leq 1\), the function becomes non-linear which requires some tweaking.

We can linearize the function by updating the penality function iteratively, commonly known as an Iterative Re-weighted Least-Squares (IRLS) method:

where we added the superscript \(\square^{(k)}\) to denote the IRLS iterations. The weights \(r(x)\) are computed from model values obtained at a previous iteration such that:

where \({r}(x) \in \mathbb{R}^{nc}\).

In matrix form, our objective function simply becomes:

where the IRLS weights \(\mathbf{R}_s\) and \(\mathbf{R}_i\) are diagonal matrices defined as:

we added two scaling parameters \(\eta_p\) and \(\eta_q\) for reasons that we won’t dicuss here, but turn out to be important to get stable solves.

In order to initialize the IRLS and get an estimate for the stabilizing parameters \(\epsilon_p\) and \(\epsilon_q\), we first invert with the smooth \(l_2\)-norm. The whole IRLS process is implemented with a directive added to the inversion workflow (see below).

In [8]:
# It is potential fields, so we will need to push the inverison down
# Create distance weights from our linera forward operator
wr = np.sum(prob.G**2.,axis=0)**0.5
wr = ( wr/np.max(wr) )

reg = Regularization.Sparse(mesh, indActive=actv, mapping=idenMap)
reg.norms = driver.lpnorms
reg.eps_p = 5e-3
reg.eps_q = 5e-3
reg.cell_weights = wr

dmis = DataMisfit.l2_DataMisfit(survey)
dmis.W = 1/wd

# Add directives to the inversion
opt = Optimization.ProjectedGNCG(maxIter=100 ,lower=0.,upper=1., maxIterLS = 20, maxIterCG= 10, tolCG = 1e-3)
invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
betaest = Directives.BetaEstimate_ByEig()

# Here is where the norms are applied
# Use pick a treshold parameter empirically based on the distribution of model
# parameters (run last cell to see the histogram before and after IRLS)
IRLS = Directives.Update_IRLS(f_min_change = 1e-4, minGNiter=2)
update_Jacobi = Directives.Update_lin_PreCond()
inv = Inversion.BaseInversion(invProb, directiveList=[betaest,IRLS,update_Jacobi])

m0 = np.ones(idenMap.nP)*1e-4

SimPEG.DataMisfit.l2_DataMisfit assigning default eps of 1e-5 * ||dobs||
In [9]:
# Run inversion...
mrec = inv.run(m0)
SimPEG.InvProblem will set Regularization.mref to m0.

    SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.
    ***Done using same Solver and solverOpts as the problem***
model has any nan: 0
=============================== Projected GNCG ===============================
  #     beta     phi_d     phi_m       f      |proj(x-g)-x|  LS    Comment
-----------------------------------------------------------------------------
x0 has any nan: 0
   0  1.43e+09  2.87e+05  0.00e+00  2.87e+05    1.45e+02      0
   1  7.17e+08  2.44e+05  1.40e-05  2.54e+05    1.44e+02      0
   2  3.59e+08  2.13e+05  4.63e-05  2.29e+05    1.43e+02      0   Skip BFGS
   3  1.79e+08  1.68e+05  1.37e-04  1.93e+05    1.42e+02      0   Skip BFGS
   4  8.97e+07  1.17e+05  3.41e-04  1.48e+05    1.40e+02      0   Skip BFGS
   5  4.48e+07  7.23e+04  7.00e-04  1.04e+05    1.37e+02      0   Skip BFGS
   6  2.24e+07  4.05e+04  1.20e-03  6.74e+04    1.33e+02      0   Skip BFGS
   7  1.12e+07  2.13e+04  1.80e-03  4.14e+04    1.29e+02      0   Skip BFGS
   8  5.60e+06  1.03e+04  2.48e-03  2.42e+04    1.24e+02      0   Skip BFGS
   9  2.80e+06  4.48e+03  3.17e-03  1.34e+04    1.22e+02      0   Skip BFGS
  10  1.40e+06  1.81e+03  3.82e-03  7.15e+03    1.20e+02      0   Skip BFGS
  11  7.01e+05  6.84e+02  4.35e-03  3.73e+03    1.19e+02      0   Skip BFGS
  12  3.50e+05  2.57e+02  4.75e-03  1.92e+03    1.17e+02      0   Skip BFGS
  13  1.75e+05  9.70e+01  5.04e-03  9.80e+02    1.14e+02      0   Skip BFGS
Convergence with smooth l2-norm regularization: Start IRLS steps...
L[p qx qy qz]-norm : [ 0.  1.  1.  1.  1.]
eps_p: 0.005 eps_q: 0.005
Regularization decrease: 8.382e+01
  14  1.75e+05  3.88e+01  5.25e-03  9.59e+02    4.07e+01      0   Skip BFGS
  15  1.75e+05  3.49e+01  4.84e-03  8.82e+02    1.31e+02      0
Regularization decrease: 1.453e-01
  16  4.83e+05  3.06e+01  4.84e-03  2.37e+03    5.17e+01      0
  17  4.83e+05  2.94e+02  4.15e-03  2.30e+03    6.04e+01      0
Regularization decrease: 1.722e-01
  18  3.31e+05  1.23e+02  4.15e-03  1.50e+03    1.05e+02      0
  19  3.31e+05  1.56e+02  4.01e-03  1.48e+03    1.43e+02      1   Skip BFGS
Regularization decrease: 8.029e-02
  20  3.31e+05  8.46e+01  4.01e-03  1.41e+03    6.73e+01      0
  21  3.31e+05  9.35e+01  3.98e-03  1.41e+03    6.64e+01      4   Skip BFGS
Regularization decrease: 7.463e-02
  22  1.98e+05  1.42e+02  3.98e-03  9.29e+02    1.36e+02      0
  23  1.98e+05  4.99e+01  4.14e-03  8.68e+02    1.31e+02      0
Regularization decrease: 8.676e-03
  24  3.04e+05  5.49e+01  4.09e-03  1.30e+03    5.99e+01      3   Skip BFGS
  25  3.04e+05  1.00e+02  3.93e-03  1.30e+03    6.56e+01      1
Regularization decrease: 7.008e-02
  26  2.35e+05  1.10e+02  3.93e-03  1.03e+03    1.42e+02      0
  27  2.35e+05  9.85e+01  3.93e-03  1.02e+03    1.44e+02      0
Regularization decrease: 2.251e-02
  28  1.92e+05  1.03e+02  3.84e-03  8.41e+02    1.08e+02      2
  29  1.92e+05  1.02e+02  3.83e-03  8.38e+02    7.06e+01      2   Skip BFGS
Regularization decrease: 5.241e-03
  30  2.33e+05  6.96e+01  3.83e-03  9.63e+02    1.12e+02      0   Skip BFGS
  31  2.33e+05  7.86e+01  3.69e-03  9.38e+02    1.45e+02      0
Regularization decrease: 4.601e-02
  32  2.33e+05  8.26e+01  3.65e-03  9.33e+02    1.35e+02      3   Skip BFGS
  33  2.33e+05  8.61e+01  3.62e-03  9.30e+02    7.51e+01      2
Reach maximum number of IRLS cycles: 10
------------------------- STOP! -------------------------
1 : |fc-fOld| = 0.0000e+00 <= tolF*(1+|f0|) = 2.8740e+04
1 : |xc-x_last| = 2.1463e-02 <= tolX*(1+|x0|) = 1.0168e-01
0 : |proj(x-g)-x|    = 7.5108e+01 <= tolG          = 1.0000e-01
0 : |proj(x-g)-x|    = 7.5108e+01 <= 1e3*eps       = 1.0000e-02
0 : maxIter   =     100    <= iter          =     34
------------------------- DONE! -------------------------
In [10]:
# Get the final model back to full space and plot!!
m_lp = actvMap*mrec
m_lp[m_lp==ndv] = np.nan

# Get the smooth model aslo
m_l2 = actvMap*IRLS.l2model
m_l2[m_l2==ndv] = np.nan

m_true = actvMap*driver.m0
m_true[m_true==ndv] = np.nan
#[OPTIONAL] Save both models to file
#Mesh.TensorMesh.writeModelUBC(mesh,'SimPEG_MAG_l2l2.sus',m_l2)
#Mesh.TensorMesh.writeModelUBC(mesh,'SimPEG_MAG_lplq.sus',m_lp)

# Plot the recoverd models
vmin, vmax = 0., 0.015

mesh = Mesh.TensorMesh([mesh.hx, mesh.hy, mesh.hz],x0="CCN")

def slide(s,normal):

    if normal == "Z":
        fig = plt.figure(figsize=(10*1.2, 8))
    else:
        fig = plt.figure(figsize=(10*1.2, 4))

    ax1 = plt.subplot(2,2,3)
    dat = mesh.plotSlice(m_lp, ax = ax1, normal=normal, ind=s, clim=np.r_[vmin, vmax],pcolorOpts={'cmap':'viridis'})
#     plt.colorbar(dat[0])
    plt.gca().set_aspect('equal')
    plt.title('Compact model')

    if normal == "Z":
        plt.xlim(-600, 600)
        plt.ylim(-600, 600.)
    else:
        plt.xlim(-600, 600)
        plt.ylim(-500, 0.)

    ax2 = plt.subplot(2,2,1)
    dat = mesh.plotSlice(m_l2, ax = ax2, normal=normal, ind=s, clim=np.r_[vmin, vmax],pcolorOpts={'cmap':'viridis'})
#     plt.colorbar(dat[0])
    plt.gca().set_aspect('equal')
    plt.title('Smooth model')

    if normal == "Z":
        plt.xlim(-600, 600)
        plt.ylim(-600, 600.)
    else:
        plt.xlim(-600, 600)
        plt.ylim(-500, 0.)

    ax2.set_xticklabels([])

    ax2 = plt.subplot(1,2,2)
    dat = mesh.plotSlice(m_true, ax = ax2, normal=normal, ind=s, clim=np.r_[vmin, vmax],pcolorOpts={'cmap':'viridis'})
#     plt.colorbar(dat[0])
    plt.gca().set_aspect('equal')
    plt.title('True model')

    pos =  ax2.get_position()

    ax2.yaxis.set_visible(False)
    if normal == "Z":
        plt.xlim(-600, 600)
        plt.ylim(-600, 600.)
        ax2.set_position([pos.x0 -0.04 , pos.y0,  pos.width, pos.height])
    else:
        plt.xlim(-600, 600)
        plt.ylim(-500, 0.)

    pos =  ax2.get_position()
    cbarax = fig.add_axes([pos.x0 + 0.375 , pos.y0 + 0.05,  pos.width*0.1, pos.height*0.75])  ## the parameters are the specified position you set
    cb = fig.colorbar(dat[0],cax=cbarax, orientation="vertical", ax = ax2, ticks=np.linspace(vmin,vmax, 4))
    cb.set_label("Susceptibility (SI)",size=12)

    #{OPTIONAL} Save the figure to png
    #fig.savefig('PF_Compact.png',dpi = 150)

interact(slide, s=(0,36), normal=['Y','Z','X'])

# interact(lambda ind: viz(m_l2, ind, normal="Z"), ind=IntSlider(min=0, max=32,step=1, value=28))

Out[10]:
<function __main__.slide>
../../_images/case-studies_PF_TKC_Mag_13_1.png
In [11]:
# Lets compare the distribution of model parameters and model gradients
plt.figure(figsize=[15,5])
ax = plt.subplot(121)
plt.hist(IRLS.l2model,100)
plt.plot((reg.eps_p,reg.eps_p),(0,mesh.nC),'r--')
plt.yscale('log', nonposy='clip')
plt.title('Hist model - l2-norm')
plt.legend(['$\epsilon_p$'])

ax = plt.subplot(122)
plt.hist(reg.regmesh.cellDiffxStencil*IRLS.l2model,100)
plt.plot((reg.eps_q,reg.eps_q),(0,mesh.nC),'r--')
plt.yscale('log', nonposy='clip')
plt.title('Hist model gradient - l2-norm')
plt.legend(['$\epsilon_q$'])

# Lets look at the distribution of model parameters and model gradients
plt.figure(figsize=[15,5])
ax = plt.subplot(121)
plt.hist(mrec,100)
plt.plot((reg.eps_p,reg.eps_p),(0,mesh.nC),'r--')
plt.yscale('log', nonposy='clip')
plt.title('Hist model - lp-norm')
plt.legend(['$\epsilon_p$'])

ax = plt.subplot(122)
plt.hist(reg.regmesh.cellDiffxStencil*mrec,100)
plt.plot((reg.eps_q,reg.eps_q),(0,mesh.nC),'r--')
plt.yscale('log', nonposy='clip')
plt.title('Hist model gradient - lp-norm')
plt.legend(['$\epsilon_q$'])
Out[11]:
<matplotlib.legend.Legend at 0x12236c9d0>
../../_images/case-studies_PF_TKC_Mag_14_1.png
../../_images/case-studies_PF_TKC_Mag_14_2.png
In [ ]: