Magnetic Amplitude inversion on a TreeMesh

In this example, we demonstrate the use of magnetic amplitude inversion on 3D TreeMesh for the inversion of Total Magnetic Intensity (TMI) data affected by remanence. The original idea must be credited to Shearer and Li (2005) @ CSM

First we invert the TMI for an equivalent source layer, from which we recover 3-component magnetic data. This data is then transformed to amplitude

Secondly, we invert the non-linear inverse problem with SimPEG.Directives.UpdateSensitivityWeights. We also uses the SimPEG.Regularization.Sparse to apply sparsity assumption in order to improve the recovery of a cube prism.

import scipy as sp
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import NearestNDInterpolator
from SimPEG import (Mesh, Directives, Maps,
                    InvProblem, Optimization, DataMisfit,
                    Inversion, Utils, Regularization)

import SimPEG.PF as PF
from SimPEG.Utils import mkvc

# sphinx_gallery_thumbnail_number = 4

Setup

Define the survey and model parameters

First we need to define the direction of the inducing field As a simple case, we pick a vertical inducing field of magnitude 50,000 nT.

# We will assume a vertical inducing field
H0 = (50000., 90., 0.)

# The magnetization is set along a different direction (induced + remanence)
M = np.array([45., 90.])

# Block with an effective susceptibility
chi_e = 0.05

# Create grid of points for topography
# Lets create a simple Gaussian topo and set the active cells
[xx, yy] = np.meshgrid(np.linspace(-200, 200, 50), np.linspace(-200, 200, 50))
b = 100
A = 50
zz = A*np.exp(-0.5*((xx/b)**2. + (yy/b)**2.))
topo = np.c_[Utils.mkvc(xx), Utils.mkvc(yy), Utils.mkvc(zz)]

# Create and array of observation points
xr = np.linspace(-100., 100., 20)
yr = np.linspace(-100., 100., 20)
X, Y = np.meshgrid(xr, yr)
Z = A*np.exp(-0.5*((X/b)**2. + (Y/b)**2.)) + 5

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

# Here how the topography looks with a quick interpolation, just a Gaussian...
tri = sp.spatial.Delaunay(topo)
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1, projection='3d')
ax.plot_trisurf(
    topo[:, 0], topo[:, 1], topo[:, 2],
    triangles=tri.simplices, cmap=plt.cm.Spectral
)
ax.scatter3D(rxLoc[:, 0], rxLoc[:, 1], rxLoc[:, 2], c='k')
plt.show()
../../../_images/sphx_glr_plot_nonLinear_Amplitude_001.png

Inversion Mesh

Here, we create a TreeMesh with base cell size of 5 m. We reated a small utility function to center the mesh around points and to figure out the outer most dimension for adequate padding distance. The second stage allows to refine the mesh around points or surfaces (point assumed to follow an horiontal interface such as topo)

# Create a mesh
h = [5, 5, 5]
padDist = np.ones((3, 2)) * 100
nCpad = [4, 4, 2]

# Get extent of points
limx = np.r_[topo[:, 0].max(), topo[:, 0].min()]
limy = np.r_[topo[:, 1].max(), topo[:, 1].min()]
limz = np.r_[topo[:, 2].max(), topo[:, 2].min()]

# Get center of the mesh
midX = np.mean(limx)
midY = np.mean(limy)
midZ = np.mean(limz)

nCx = int(limx[0]-limx[1]) / h[0]
nCy = int(limy[0]-limy[1]) / h[1]
nCz = int(limz[0]-limz[1]+int(np.min(np.r_[nCx, nCy])/3)) / h[2]

# Figure out full extent required from input
extent = np.max(np.r_[nCx * h[0] + padDist[0, :].sum(),
                      nCy * h[1] + padDist[1, :].sum(),
                      nCz * h[2] + padDist[2, :].sum()])

maxLevel = int(np.log2(extent/h[0]))+1

# Number of cells at the small octree level
# For now equal in 3D
nCx, nCy, nCz = 2**(maxLevel), 2**(maxLevel), 2**(maxLevel)

# Define the mesh and origin
mesh = Mesh.TreeMesh([np.ones(nCx)*h[0],
                      np.ones(nCx)*h[1],
                      np.ones(nCx)*h[2]])

# Set origin
mesh.x0 = np.r_[-nCx*h[0]/2.+midX, -nCy*h[1]/2.+midY, -nCz*h[2]/2.+midZ]

# Refine the mesh around topography
# Get extent of points
F = NearestNDInterpolator(topo[:, :2], topo[:, 2])
zOffset = 0
# Cycle through the first 3 octree levels
for ii in range(3):

    dx = mesh.hx.min()*2**ii

    nCx = int((limx[0]-limx[1]) / dx)
    nCy = int((limy[0]-limy[1]) / dx)

    # Create a grid at the octree level in xy
    CCx, CCy = np.meshgrid(
        np.linspace(limx[1], limx[0], nCx),
        np.linspace(limy[1], limy[0], nCy)
    )

    z = F(mkvc(CCx), mkvc(CCy))

    # level means number of layers in current OcTree level
    for level in range(int(nCpad[ii])):

        mesh.insert_cells(
            np.c_[mkvc(CCx), mkvc(CCy), z-zOffset],
            np.ones_like(z)*maxLevel-ii,
            finalize=False
        )

        zOffset += dx

mesh.finalize()

# Define an active cells from topo
actv = Utils.surface2ind_topo(mesh, topo)
nC = int(actv.sum())

Forward modeling data

We can now generate TMI data

# Convert the inclination declination to vector in Cartesian
M_xyz = Utils.matutils.dip_azimuth2cartesian(np.ones(nC)*M[0], np.ones(nC)*M[1])

# Get the indicies of the magnetized block
ind = Utils.ModelBuilder.getIndicesBlock(
    np.r_[-20, -20, -10], np.r_[20, 20, 25],
    mesh.gridCC,
)[0]

# Assign magnetization value, inducing field strength will
# be applied in by the :class:`SimPEG.PF.Magnetics` problem
model = np.zeros(mesh.nC)
model[ind] = chi_e

# Remove air cells
model = model[actv]

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

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

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

# Compute some data and add some random noise
data = prob.fields(model)

# Split the data in components
nD = rxLoc.shape[0]

std = 5  # nT
data += np.random.randn(nD)*std
wd = np.ones(nD)*std

# Assigne data and uncertainties to the survey
survey.dobs = data
survey.std = wd


# Plot the model and data
plt.figure(figsize=(8, 8))
ax = plt.subplot(2, 1, 1)
im = Utils.PlotUtils.plot2Ddata(
        rxLoc, data, ax=ax, contourOpts={"cmap": "RdBu_r"}
)
plt.colorbar(im[0])
ax.set_title('Predicted data.')
plt.gca().set_aspect('equal', adjustable='box')

# Plot the vector model
ax = plt.subplot(2, 1, 2)

# Create active map to go from reduce set to full
actvPlot = Maps.InjectActiveCells(mesh, actv, np.nan)
mesh.plotSlice(
    actvPlot*model, ax=ax, normal='Y', ind=66,
    pcolorOpts={"vmin": 0., "vmax": 0.01}, grid=True,
)
ax.set_xlim([-200, 200])
ax.set_ylim([-100, 75])
ax.set_xlabel('x')
ax.set_ylabel('y')
plt.gca().set_aspect('equal', adjustable='box')

plt.show()
../../../_images/sphx_glr_plot_nonLinear_Amplitude_002.png

Out:

Begin forward: M=H0, Rx type= tmi
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 %

Equivalent Source

We first need to convert the TMI data into amplitude. We do this by for an effective susceptibility layer, from which we can forward component data

# Get the active cells for equivalent source is the top only
surf = Utils.modelutils.surface_layer_index(mesh, topo)
nC = np.count_nonzero(surf)  # Number of active cells

# Create active map to go from reduce set to full
surfMap = Maps.InjectActiveCells(mesh, surf, np.nan)

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

# Create static map
prob = PF.Magnetics.MagneticIntegral(
        mesh, chiMap=idenMap, actInd=surf,
        parallelized=False, equiSourceLayer=True)

prob.solverOpts['accuracyTol'] = 1e-4

# Pair the survey and problem
if survey.ispaired:
    survey.unpair()
survey.pair(prob)


# Create a regularization function, in this case l2l2
reg = Regularization.Sparse(
    mesh, indActive=surf, mapping=Maps.IdentityMap(nP=nC), scaledIRLS=False
)
reg.mref = np.zeros(nC)

# Specify how the optimization will proceed, set susceptibility bounds to inf
opt = Optimization.ProjectedGNCG(
    maxIter=20, lower=-np.inf, upper=np.inf, maxIterLS=20,
    maxIterCG=20, tolCG=1e-3
)

# Define misfit function (obs-calc)
dmis = DataMisfit.l2_DataMisfit(survey)
dmis.W = 1./survey.std

# Create the default L2 inverse problem from the above objects
invProb = InvProblem.BaseInvProblem(dmis, reg, opt)

# Specify how the initial beta is found
betaest = Directives.BetaEstimate_ByEig()

# Target misfit to stop the inversion,
# try to fit as much as possible of the signal, we don't want to lose anything
IRLS = Directives.Update_IRLS(f_min_change=1e-3, minGNiter=1,
                              beta_tol=1e-1)
update_Jacobi = Directives.UpdatePreconditioner()
# Put all the parts together
inv = Inversion.BaseInversion(invProb,
                              directiveList=[betaest, IRLS, update_Jacobi])

# Run the equivalent source inversion
mstart = np.ones(nC)*1e-4
mrec = inv.run(mstart)

Out:

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

    SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.
    ***Done using same Solver and solverOpts as the problem***
Begin forward: M=H0, Rx type= tmi
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 %
Approximated diag(JtJ) with linear operator
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  5.96e+04  9.64e+03  5.05e-02  1.27e+04    7.18e+06      0
   1  2.98e+04  4.26e+03  3.68e-02  5.35e+03    9.90e+04      0
   2  1.49e+04  2.47e+03  7.96e-02  3.65e+03    7.25e+04      0   Skip BFGS
   3  7.45e+03  1.19e+03  1.40e-01  2.23e+03    4.80e+04      0   Skip BFGS
   4  3.73e+03  5.03e+02  2.04e-01  1.26e+03    2.91e+04      0   Skip BFGS
Reached starting chifact with l2-norm regularization: Start IRLS steps...
eps_p: 0.00745756694134104 eps_q: 0.00745756694134104
delta phim: 4.430e+01
   5  1.86e+03  1.98e+02  2.60e-01  6.82e+02    1.66e+04      0   Skip BFGS
Beta search step
   5  4.34e+03  1.98e+02  2.60e-01  1.33e+03    5.51e+03      0   Skip BFGS
Beta search step
   5  3.41e+03  1.98e+02  2.60e-01  1.08e+03    2.85e+03      0
Beta search step
   5  5.35e+03  1.98e+02  2.60e-01  1.59e+03    1.45e+04      0
Beta search step
   5  3.66e+03  1.98e+02  2.60e-01  1.15e+03    6.01e+02      0
delta phim: 2.762e-01
   6  3.66e+03  1.93e+02  2.61e-01  1.15e+03    6.09e-02      0
------------------------- STOP! -------------------------
1 : |fc-fOld| = 4.4403e-02 <= tolF*(1+|f0|) = 1.2653e+03
1 : |xc-x_last| = 1.6684e-04 <= tolX*(1+|x0|) = 1.0085e-01
1 : |proj(x-g)-x|    = 6.0870e-02 <= tolG          = 1.0000e-01
0 : |proj(x-g)-x|    = 6.0870e-02 <= 1e3*eps       = 1.0000e-02
0 : maxIter   =      20    <= iter          =      6
------------------------- DONE! -------------------------

Forward Amplitude Data

Now that we have an equialent source layer, we can forward model alh three components of the field and add them up: \(|B| = \sqrt{( Bx^2 + Bx^2 + Bx^2 )}\)

prob.forwardOnly = True
prob.rx_type = 'xyz'
prob._G = None
prob.modelType = 'amplitude'
prob.model = mrec
pred = prob.fields(mrec)

bx = pred[:nD]
by = pred[nD:2*nD]
bz = pred[2*nD:]

bAmp = (bx**2. + by**2. + bz**2.)**0.5


# Plot the layer model and data
plt.figure(figsize=(8, 8))
ax = plt.subplot(2, 2, 1)
im = Utils.PlotUtils.plot2Ddata(
        rxLoc, invProb.dpred, ax=ax, contourOpts={"cmap": "RdBu_r"}
)
plt.colorbar(im[0])
ax.set_title('Predicted data.')
plt.gca().set_aspect('equal', adjustable='box')

ax = plt.subplot(2, 2, 2)
im = Utils.PlotUtils.plot2Ddata(
        rxLoc, bAmp, ax=ax, contourOpts={"cmap": "RdBu_r"}
)
plt.colorbar(im[0])
ax.set_title('Calculated amplitude')
plt.gca().set_aspect('equal', adjustable='box')

# Plot the equivalent layer model
ax = plt.subplot(2, 1, 2)
mesh.plotSlice(
    surfMap*mrec, ax=ax, normal='Y', ind=66,
    pcolorOpts={"vmin": 0., "vmax": 0.01}, grid=True
)
ax.set_xlim([-200, 200])
ax.set_ylim([-100, 75])
ax.set_xlabel('x')
ax.set_ylabel('y')
plt.gca().set_aspect('equal', adjustable='box')

plt.show()
../../../_images/sphx_glr_plot_nonLinear_Amplitude_003.png

Out:

Begin forward: M=H0, Rx type= xyz
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 %

Amplitude Inversion

Now that we have amplitude data, we can invert for an effective susceptibility. This is a non-linear inversion.

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

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

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

# Create the forward model operator
prob = PF.Magnetics.MagneticIntegral(
    mesh, chiMap=idenMap, actInd=actv,
    modelType='amplitude', rx_type='xyz'
)
prob.model = mstart
# Change the survey to xyz components
surveyAmp = PF.BaseMag.LinearSurvey(survey.srcField)

# Pair the survey and problem
surveyAmp.pair(prob)
# Create a regularization function, in this case l2l2
wr = np.sum(prob.G**2., axis=0)**0.5
wr = (wr/np.max(wr))
# Re-set the observations to |B|
surveyAmp.dobs = bAmp
surveyAmp.std = wd

# Create a sparse regularization
reg = Regularization.Sparse(mesh, indActive=actv, mapping=idenMap)
reg.norms = np.c_[1, 0, 0, 0]
reg.mref = np.zeros(nC)
reg.cell_weights = wr
# Data misfit function
dmis = DataMisfit.l2_DataMisfit(surveyAmp)
dmis.W = 1./surveyAmp.std

# Add directives to the inversion
opt = Optimization.ProjectedGNCG(maxIter=30, lower=0., upper=1.,
                                 maxIterLS=20, maxIterCG=20,
                                 tolCG=1e-3)

invProb = InvProblem.BaseInvProblem(dmis, reg, opt)

# Here is the list of directives
betaest = Directives.BetaEstimate_ByEig()

# Specify the sparse norms
IRLS = Directives.Update_IRLS(f_min_change=1e-3,
                              minGNiter=1, coolingRate=1,
                              betaSearch=False)

# Special directive specific to the mag amplitude problem. The sensitivity
# weights are update between each iteration.
update_SensWeight = Directives.UpdateSensitivityWeights()
update_Jacobi = Directives.UpdatePreconditioner()

# Put all together
inv = Inversion.BaseInversion(
    invProb, directiveList=[
        betaest, IRLS, update_SensWeight, update_Jacobi
        ]
)

# Invert
mrec_Amp = inv.run(mstart)

Out:

Begin forward: M=H0, Rx type= xyz
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 %
SimPEG.DataMisfit.l2_DataMisfit assigning default eps of 1e-5 * ||dobs||

    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  2.23e+08  1.47e+04  3.64e-06  1.55e+04    5.77e+01      0
   1  1.11e+08  1.11e+04  8.91e-06  1.21e+04    8.34e+01      0
   2  5.57e+07  8.34e+03  2.64e-05  9.81e+03    7.41e+01      0   Skip BFGS
   3  2.78e+07  5.72e+03  5.98e-05  7.38e+03    6.80e+01      0   Skip BFGS
   4  1.39e+07  3.50e+03  1.16e-04  5.11e+03    7.62e+01      0   Skip BFGS
   5  6.96e+06  1.88e+03  1.99e-04  3.27e+03    1.07e+02      0   Skip BFGS
   6  3.48e+06  8.78e+02  3.00e-04  1.92e+03    1.32e+02      0   Skip BFGS
   7  1.74e+06  3.68e+02  3.95e-04  1.06e+03    1.31e+02      0   Skip BFGS
Reached starting chifact with l2-norm regularization: Start IRLS steps...
eps_p: 0.00785326494414669 eps_q: 0.00785326494414669
delta phim: 4.815e+01
   8  8.70e+05  1.54e+02  5.60e-04  6.41e+02    1.03e+02      0   Skip BFGS
delta phim: 3.968e+03
   9  1.94e+06  8.12e+01  7.53e-04  1.54e+03    6.08e+01      0   Skip BFGS
delta phim: 7.986e-01
  10  1.50e+06  2.53e+02  6.77e-04  1.27e+03    8.39e+01      0
delta phim: 1.737e-01
  11  1.50e+06  2.08e+02  8.01e-04  1.41e+03    4.61e+01      0
delta phim: 4.740e-01
  12  1.18e+06  2.41e+02  8.70e-04  1.27e+03    7.54e+01      0
delta phim: 3.663e-01
  13  1.18e+06  2.07e+02  1.01e-03  1.40e+03    4.59e+01      0
delta phim: 4.073e-01
  14  9.41e+05  2.37e+02  1.10e-03  1.28e+03    6.93e+01      0
delta phim: 2.936e-01
  15  9.41e+05  2.06e+02  1.27e-03  1.40e+03    5.15e+01      0
delta phim: 2.757e-01
  16  7.63e+05  2.29e+02  1.39e-03  1.29e+03    6.37e+01      0
delta phim: 1.533e-01
  17  7.63e+05  2.00e+02  1.62e-03  1.44e+03    6.28e+01      0
delta phim: 8.857e-02
  18  7.63e+05  2.18e+02  1.84e-03  1.62e+03    6.76e+01      0
delta phim: 3.490e-02
  19  6.00e+05  2.44e+02  2.02e-03  1.45e+03    6.72e+01      0
delta phim: 9.743e-02
  20  6.00e+05  2.19e+02  2.23e-03  1.56e+03    6.05e+01      0
delta phim: 8.281e-02
  21  4.62e+05  2.53e+02  2.31e-03  1.32e+03    6.35e+01      0
delta phim: 7.610e-02
  22  3.73e+05  2.31e+02  2.48e-03  1.16e+03    5.82e+01      0   Skip BFGS
delta phim: 6.047e-02
  23  3.73e+05  2.20e+02  2.64e-03  1.20e+03    5.47e+01      0   Skip BFGS
delta phim: 1.812e-01
  24  2.94e+05  2.42e+02  2.73e-03  1.04e+03    6.85e+01      0
delta phim: 1.693e-01
  25  2.40e+05  2.26e+02  2.89e-03  9.21e+02    6.21e+01      0
delta phim: 2.618e-01
  26  2.40e+05  2.14e+02  3.02e-03  9.39e+02    5.84e+01      0   Skip BFGS
delta phim: 2.761e-01
  27  1.97e+05  2.24e+02  3.06e-03  8.27e+02    5.45e+01      0
Reach maximum number of IRLS cycles: 20
------------------------- STOP! -------------------------
1 : |fc-fOld| = 0.0000e+00 <= tolF*(1+|f0|) = 1.5549e+03
1 : |xc-x_last| = 1.7773e-03 <= tolX*(1+|x0|) = 1.0194e-01
0 : |proj(x-g)-x|    = 5.4464e+01 <= tolG          = 1.0000e-01
0 : |proj(x-g)-x|    = 5.4464e+01 <= 1e3*eps       = 1.0000e-02
0 : maxIter   =      30    <= iter          =     28
------------------------- DONE! -------------------------

Final Plot

Let’s compare the smooth and compact model Note that the recovered effective susceptibility block is slightly offseted to the left of the true model. This is due to the wrong assumption of a vertical magnetization. Important to remember that the amplitude inversion is weakly sensitive to the magnetization direction, but can still have an impact.

# Plot the layer model and data
plt.figure(figsize=(12, 8))
ax = plt.subplot(3, 1, 1)
im = Utils.PlotUtils.plot2Ddata(
        rxLoc, invProb.dpred, ax=ax, contourOpts={"cmap": "RdBu_r"}
 )
plt.colorbar(im[0])
ax.set_title('Predicted data.')
plt.gca().set_aspect('equal', adjustable='box')

# Plot the vector model
ax = plt.subplot(3, 1, 2)
im = mesh.plotSlice(
    actvPlot*invProb.l2model, ax=ax, normal='Y', ind=66,
    pcolorOpts={"vmin": 0., "vmax": 0.01}, grid=True,
)
plt.colorbar(im[0])
ax.set_xlim([-200, 200])
ax.set_ylim([-100, 75])
ax.set_xlabel('x')
ax.set_ylabel('y')
plt.gca().set_aspect('equal', adjustable='box')

# Plot the amplitude model
ax = plt.subplot(3, 1, 3)
im = mesh.plotSlice(
    actvPlot*mrec_Amp, ax=ax, normal='Y', ind=66,
    pcolorOpts={"vmin": 0., "vmax": 0.01}, grid=True
)
plt.colorbar(im[0])
ax.set_xlim([-200, 200])
ax.set_ylim([-100, 75])
ax.set_xlabel('x')
ax.set_ylabel('y')
plt.gca().set_aspect('equal', adjustable='box')
plt.show()
../../../_images/sphx_glr_plot_nonLinear_Amplitude_004.png

Total running time of the script: ( 1 minutes 55.484 seconds)

Gallery generated by Sphinx-Gallery