# Linear Magnetic InversionΒΆ

Objective:

In this tutorial we will create a simple magnetic problem from scratch using the SimPEG framework.

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.

The remaining of this notebook goes through all the important components of a 3D magnetic experiment. From mesh creation, topography, data and inverse problem.

Enjoy.

In [1]:

from SimPEG import Mesh
from SimPEG.Utils import mkvc, surface2ind_topo
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
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

In [2]:

# First we need to define the direction of the inducing field
# As a simple case, we pick a vertical inducing field of magnitude 50,000nT.
# From old convention, field orientation is given as an azimuth from North
# (positive clockwise) and dip from the horizontal (positive downward).
H0 = (60000.,90.,0.)

# Create a mesh
dx    = 5.

hxind = [(dx,5,-1.3), (dx, 10), (dx,5,1.3)]
hyind = [(dx,5,-1.3), (dx, 10), (dx,5,1.3)]
hzind = [(dx,5,-1.3),(dx, 10)]

mesh = Mesh.TensorMesh([hxind, hyind, hzind], 'CCC')

# Get index of the center
midx = int(mesh.nCx/2)
midy = int(mesh.nCy/2)

# Lets create a simple Gaussian topo and set the active cells
[xx,yy] = np.meshgrid(mesh.vectorNx,mesh.vectorNy)
zz = -np.exp( ( xx**2 + yy**2 )/ 75**2 ) + mesh.vectorNz[-1]

topo = np.c_[mkvc(xx),mkvc(yy),mkvc(zz)] # We would usually load a topofile

actv = surface2ind_topo(mesh,topo,'N') # Go from topo to actv cells
actv = np.asarray([inds for inds, elem in enumerate(actv, 1) if elem], dtype = int) - 1
#nC   = mesh.nC
#actv = np.asarray(range(mesh.nC))

# Create active map to go from reduce space to full
actvMap = Maps.InjectActiveCells(mesh, actv, -100)
nC = len(actv)

# Create and array of observation points
xr = np.linspace(-20., 20., 20)
yr = np.linspace(-20., 20., 20)
X, Y = np.meshgrid(xr, yr)

# Let just put the observation above the topo
Z = -np.exp( ( X**2 + Y**2 )/ 75**2 ) + mesh.vectorNz[-1] + 5.
#Z = np.ones(shape(X)) * mesh.vectorCCz[-1]

# Create a MAGsurvey
rxLoc = np.c_[mkvc(X.T), mkvc(Y.T), mkvc(Z.T)]
rxLoc = PF.BaseMag.RxObs(rxLoc)
srcField = PF.BaseMag.SrcField([rxLoc],param = H0)
survey = PF.BaseMag.LinearSurvey(srcField)



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 look 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 [3]:

# We can now create a susceptibility model and generate data
model = np.zeros((mesh.nCx,mesh.nCy,mesh.nCz))
model[(midx-2):(midx+2),(midy-2):(midy+2),-6:-2] = 0.02
model = mkvc(model)
model = model[actv]

# Create active map to go from reduce set to full
actvMap = Maps.InjectActiveCells(mesh, actv, -100)

# Creat reduced identity map
idenMap = Maps.IdentityMap(nP = nC)

# Create the forward model operator
prob = PF.Magnetics.MagneticIntegral(mesh, chiMap=idenMap, actInd=actv)

# Pair the survey and problem
survey.pair(prob)

# Compute linear forward operator and compute some data
d = prob.fields(model)

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!!


In [4]:

# Plot the model
m_true = actvMap * model
m_true[m_true==-100] = np.nan
plt.figure()
ax = plt.subplot(212)
mesh.plotSlice(m_true, ax = ax, normal = 'Y', ind=midy, grid=True, clim = (0., model.max()), pcolorOpts={'cmap':'viridis'})
plt.title('A simple block model.')
plt.xlabel('x'); plt.ylabel('z')

# We can now generate data
data = d + np.random.randn(len(d)) # We add some random Gaussian noise (1nT)
wd = np.ones(len(data))*1. # Assign flat uncertainties

plt.subplot(221)
plt.imshow(d.reshape(X.shape), extent=[xr.min(), xr.max(), yr.min(), yr.max()])
plt.title('True data.')
plt.colorbar()

plt.subplot(222)
plt.imshow(data.reshape(X.shape), extent=[xr.min(), xr.max(), yr.min(), yr.max()])
plt.title('Data + Noise')
plt.colorbar()

plt.tight_layout()

In [5]:

# Create distance weights from our linera forward operator
wr = np.sum(prob.G**2.,axis=0)**0.5
wr = ( wr/np.max(wr) )

wr_FULL = actvMap * wr
wr_FULL[wr_FULL==-100] = np.nan
plt.figure()
ax = plt.subplot()
mesh.plotSlice(wr_FULL, ax = ax, normal = 'Y', ind=midx, grid=True, clim = (0, wr.max()),pcolorOpts={'cmap':'viridis'})
plt.title('Distance weighting')
plt.xlabel('x');plt.ylabel('z')


Once we have our problem, we can use the inversion tools in SimPEG to run our inversion:

In [6]:

#survey.makeSyntheticData(data, std=0.01)
survey.dobs= data
survey.std = wd
survey.mtrue = model

# Create a regularization
reg = Regularization.Sparse(mesh, indActive=actv, mapping=idenMap)
reg.cell_weights = wr
reg.norms = [0, 1, 1, 1]
reg.eps_p = 1e-3
reg.eps_1 = 1e-3
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-3, minGNiter=3)
update_Jacobi = Directives.Update_lin_PreCond()
inv = Inversion.BaseInversion(invProb, directiveList=[betaest, IRLS, update_Jacobi])

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

SimPEG.DataMisfit.l2_DataMisfit assigning default eps of 1e-5 * ||dobs||

In [7]:

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.46e+10  4.13e+05  0.00e+00  4.13e+05    6.27e+01      0
1  7.32e+09  3.62e+05  1.62e-06  3.74e+05    6.26e+01      0
2  3.66e+09  3.27e+05  5.17e-06  3.46e+05    6.25e+01      0   Skip BFGS
3  1.83e+09  2.74e+05  1.57e-05  3.03e+05    6.23e+01      0   Skip BFGS
4  9.15e+08  2.09e+05  4.15e-05  2.47e+05    6.19e+01      0   Skip BFGS
5  4.58e+08  1.46e+05  9.06e-05  1.88e+05    6.10e+01      0   Skip BFGS
6  2.29e+08  9.78e+04  1.66e-04  1.36e+05    6.01e+01      0   Skip BFGS
7  1.14e+08  6.38e+04  2.71e-04  9.48e+04    5.82e+01      0   Skip BFGS
8  5.72e+07  3.89e+04  4.26e-04  6.33e+04    5.48e+01      0   Skip BFGS
9  2.86e+07  2.06e+04  6.53e-04  3.93e+04    5.08e+01      0   Skip BFGS
10  1.43e+07  9.12e+03  9.33e-04  2.25e+04    4.45e+01      0   Skip BFGS
11  7.15e+06  3.45e+03  1.21e-03  1.21e+04    3.50e+01      0   Skip BFGS
12  3.58e+06  1.19e+03  1.41e-03  6.24e+03    4.59e+01      0   Skip BFGS
13  1.79e+06  4.81e+02  1.55e-03  3.24e+03    4.36e+01      0   Skip BFGS
14  8.94e+05  2.73e+02  1.62e-03  1.72e+03    4.17e+01      0   Skip BFGS
15  4.47e+05  2.16e+02  1.67e-03  9.61e+02    3.41e+01      0   Skip BFGS
Convergence with smooth l2-norm regularization: Start IRLS steps...
L[p qx qy qz]-norm : [0 1 1 1]
eps_p: 0.001 eps_q: 0.000812873564105
Regularization decrease: 3.191e+01
16  4.47e+05  1.98e+02  1.69e-03  9.55e+02    4.43e+01      0   Skip BFGS
17  4.47e+05  2.13e+02  1.24e-03  7.65e+02    5.97e+01      0
18  4.47e+05  2.16e+02  1.21e-03  7.56e+02    5.94e+01      0   Skip BFGS
Regularization decrease: 2.527e-01
19  4.25e+05  2.10e+02  1.21e-03  7.24e+02    3.42e+01      0
20  4.25e+05  2.14e+02  1.07e-03  6.71e+02    5.91e+01      0
21  4.25e+05  2.12e+02  1.07e-03  6.68e+02    6.04e+01      0
Regularization decrease: 9.408e-02
22  4.04e+05  2.11e+02  1.07e-03  6.43e+02    2.81e+01      0
23  4.04e+05  2.13e+02  9.28e-04  5.88e+02    4.42e+01      0
24  4.04e+05  2.08e+02  9.27e-04  5.82e+02    3.26e+01      0
Regularization decrease: 1.260e-01
25  4.04e+05  2.08e+02  9.27e-04  5.82e+02    4.06e+01      0
26  4.04e+05  2.07e+02  7.88e-04  5.26e+02    3.68e+01      0
27  4.04e+05  2.02e+02  7.84e-04  5.18e+02    3.42e+01      0
Regularization decrease: 1.564e-01
28  4.04e+05  2.02e+02  7.84e-04  5.18e+02    3.52e+01      0
29  4.04e+05  1.99e+02  6.93e-04  4.79e+02    3.67e+01      0
30  4.04e+05  1.98e+02  6.88e-04  4.76e+02    4.25e+01      0
Regularization decrease: 8.832e-02
31  4.04e+05  1.98e+02  6.88e-04  4.76e+02    3.20e+01      0
32  4.04e+05  1.98e+02  6.15e-04  4.46e+02    4.01e+01      0
33  4.04e+05  1.97e+02  6.13e-04  4.45e+02    3.51e+01      0
Regularization decrease: 7.651e-02
34  4.04e+05  1.97e+02  6.13e-04  4.45e+02    3.78e+01      0
35  4.04e+05  1.99e+02  5.50e-04  4.21e+02    3.88e+01      0
36  4.04e+05  1.98e+02  5.49e-04  4.20e+02    3.54e+01      0
Regularization decrease: 1.255e-01
37  4.04e+05  1.98e+02  5.49e-04  4.20e+02    4.89e+01      0
38  4.04e+05  1.99e+02  5.27e-04  4.12e+02    3.69e+01      0
39  4.04e+05  1.99e+02  5.26e-04  4.11e+02    2.69e+01      0
Regularization decrease: 6.452e-02
40  4.04e+05  1.99e+02  5.26e-04  4.11e+02    3.58e+01      0
41  4.04e+05  1.99e+02  5.22e-04  4.10e+02    3.74e+01      0
42  4.04e+05  1.99e+02  5.21e-04  4.09e+02    3.06e+01      0
Regularization decrease: 1.349e-02
43  4.04e+05  1.99e+02  5.21e-04  4.09e+02    2.75e+01      0
44  4.04e+05  1.99e+02  5.18e-04  4.08e+02    3.45e+01      0
45  4.04e+05  1.99e+02  5.18e-04  4.08e+02    2.61e+01      0
Reach maximum number of IRLS cycles: 10
------------------------- STOP! -------------------------
1 : |fc-fOld| = 0.0000e+00 <= tolF*(1+|f0|) = 4.1262e+04
1 : |xc-x_last| = 7.8384e-04 <= tolX*(1+|x0|) = 1.0075e-01
0 : |proj(x-g)-x|    = 2.6113e+01 <= tolG          = 1.0000e-01
0 : |proj(x-g)-x|    = 2.6113e+01 <= 1e3*eps       = 1.0000e-02
0 : maxIter   =     100    <= iter          =     46
------------------------- DONE! -------------------------


Inversion has converged. We can plot sections through the model.

In [8]:

# Here is the recovered susceptibility model
ypanel = midx
zpanel = -4
m_l2 = actvMap * IRLS.l2model
m_l2[m_l2==-100] = np.nan

m_lp = actvMap * mrec
m_lp[m_lp==-100] = np.nan

m_true = actvMap * model
m_true[m_true==-100] = np.nan

plt.figure()

#Plot L2 model
ax = plt.subplot(231)
mesh.plotSlice(m_l2, ax = ax, normal = 'Z', ind=zpanel, grid=True, clim = (0., model.max()),pcolorOpts={'cmap':'viridis'})
plt.plot(([mesh.vectorCCx[0],mesh.vectorCCx[-1]]),([mesh.vectorCCy[ypanel],mesh.vectorCCy[ypanel]]),color='w')
plt.title('Plan l2-model.')
plt.gca().set_aspect('equal')
plt.ylabel('y')
ax.xaxis.set_visible(False)

# Vertica section
ax = plt.subplot(234)
mesh.plotSlice(m_l2, ax = ax, normal = 'Y', ind=midx, grid=True, clim = (0., model.max()),pcolorOpts={'cmap':'viridis'})
plt.plot(([mesh.vectorCCx[0],mesh.vectorCCx[-1]]),([mesh.vectorCCz[zpanel],mesh.vectorCCz[zpanel]]),color='w')
plt.plot(([mesh.vectorCCx[0],mesh.vectorCCx[-1]]),([Z.min(),Z.max()]),color='k')
plt.title('E-W l2-model.')
plt.gca().set_aspect('equal')
plt.xlabel('x')
plt.ylabel('z')

#Plot Lp model
ax = plt.subplot(232)
mesh.plotSlice(m_lp, ax = ax, normal = 'Z', ind=zpanel, grid=True, clim = (0., model.max()),pcolorOpts={'cmap':'viridis'})
plt.plot(([mesh.vectorCCx[0],mesh.vectorCCx[-1]]),([mesh.vectorCCy[ypanel],mesh.vectorCCy[ypanel]]),color='w')
plt.title('Plan lp-model.')
plt.gca().set_aspect('equal')
ax.xaxis.set_visible(False)
ax.yaxis.set_visible(False)

# Vertical section
ax = plt.subplot(235)
mesh.plotSlice(m_lp, ax = ax, normal = 'Y', ind=midx, grid=True, clim = (0., model.max()),pcolorOpts={'cmap':'viridis'})
plt.plot(([mesh.vectorCCx[0],mesh.vectorCCx[-1]]),([mesh.vectorCCz[zpanel],mesh.vectorCCz[zpanel]]),color='w')
plt.title('E-W lp-model.')
plt.gca().set_aspect('equal')
ax.yaxis.set_visible(False)
plt.xlabel('x')

#Plot True model
ax = plt.subplot(233)
mesh.plotSlice(m_true, ax = ax, normal = 'Z', ind=zpanel, grid=True, clim = (0., model.max()),pcolorOpts={'cmap':'viridis'})
plt.plot(([mesh.vectorCCx[0],mesh.vectorCCx[-1]]),([mesh.vectorCCy[ypanel],mesh.vectorCCy[ypanel]]),color='w')
plt.title('Plan true model.')
plt.gca().set_aspect('equal')
ax.xaxis.set_visible(False)
ax.yaxis.set_visible(False)

# Vertical section
ax = plt.subplot(236)
mesh.plotSlice(m_true, ax = ax, normal = 'Y', ind=midx, grid=True, clim = (0., model.max()),pcolorOpts={'cmap':'viridis'})
plt.plot(([mesh.vectorCCx[0],mesh.vectorCCx[-1]]),([mesh.vectorCCz[zpanel],mesh.vectorCCz[zpanel]]),color='w')
plt.title('E-W true model.')
plt.gca().set_aspect('equal')
plt.xlabel('x')
ax.yaxis.set_visible(False)



Great, we have a 3D model of susceptibility, but the job is not done yet. A VERY important step of the inversion workflow is to look at how well the model can predict the observed data. The figure below compares the observed, predicted and normalized residual.

In [9]:

# Plot predicted data and residual
plt.figure()
pred = prob.fields(mrec)  #: this is matrix multiplication!!

plt.subplot(221)
plt.imshow(data.reshape(X.shape))
plt.title('Observed data.')
plt.colorbar()

plt.subplot(222)
plt.imshow(pred.reshape(X.shape))
plt.title('Predicted data.')
plt.colorbar()

plt.subplot(223)
plt.imshow(data.reshape(X.shape) - pred.reshape(X.shape))
plt.title('Residual data.')
plt.colorbar()

plt.subplot(224)
plt.imshow( (data.reshape(X.shape) - pred.reshape(X.shape)) / wd.reshape(X.shape) )
plt.title('Normalized Residual')
plt.colorbar()

plt.tight_layout()


Good job! Hopefully we covered all the important points regarding the inversion of magnetic field data using the integral formulation.

Make sure you visit the notebook for the compact norms regularization.

Cheers!

In [10]: