Note
Go to the end to download the full example code.
Magnetic inversion on a TreeMesh#
In this example, we demonstrate the use of a Magnetic Vector Inversion on 3D TreeMesh for the inversion of magnetic data.
The inverse problem uses the simpeg.regularization.VectorAmplitude
regularization.
from simpeg import (
data,
data_misfit,
directives,
maps,
inverse_problem,
optimization,
inversion,
regularization,
)
from simpeg import utils
from simpeg.utils import mkvc, sdiag
from discretize.utils import mesh_builder_xyz, refine_tree_xyz, active_from_xyz
from simpeg.potential_fields import magnetics
import numpy as np
import matplotlib.pyplot as plt
# sphinx_gallery_thumbnail_number = 3
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.
np.random.seed(1)
# We will assume a vertical inducing field
h0_amplitude, h0_inclination, h0_declination = (50000.0, 90.0, 0.0)
# 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.0 + (yy / b) ** 2.0))
topo = np.c_[utils.mkvc(xx), utils.mkvc(yy), utils.mkvc(zz)]
# Create an array of observation points
xr = np.linspace(-100.0, 100.0, 20)
yr = np.linspace(-100.0, 100.0, 20)
X, Y = np.meshgrid(xr, yr)
Z = A * np.exp(-0.5 * ((X / b) ** 2.0 + (Y / b) ** 2.0)) + 5
# Create a MAGsurvey
xyzLoc = np.c_[mkvc(X.T), mkvc(Y.T), mkvc(Z.T)]
rxLoc = magnetics.receivers.Point(xyzLoc)
srcField = magnetics.sources.UniformBackgroundField(
receiver_list=[rxLoc],
amplitude=h0_amplitude,
inclination=h0_inclination,
declination=h0_declination,
)
survey = magnetics.survey.Survey(srcField)
Inversion Mesh#
Here, we create a TreeMesh with base cell size of 5 m.
# Create a mesh
h = [5, 5, 5]
padDist = np.ones((3, 2)) * 100
mesh = mesh_builder_xyz(
xyzLoc, h, padding_distance=padDist, depth_core=100, mesh_type="tree"
)
mesh = refine_tree_xyz(
mesh, topo, method="surface", octree_levels=[2, 6], finalize=True
)
# Define an active cells from topo
actv = active_from_xyz(mesh, topo)
nC = int(actv.sum())
/usr/share/miniconda/envs/simpeg-test/lib/python3.10/site-packages/discretize/utils/mesh_utils.py:528: FutureWarning:
In discretize v1.0 the TreeMesh will change the default value of diagonal_balance to True, which will likely slightly change meshes you have previously created. If you need to keep the current behavior, explicitly set diagonal_balance=False.
/home/vsts/work/1/s/examples/03-magnetics/plot_inv_mag_MVI_VectorAmplitude.py:88: DeprecationWarning:
The surface option is deprecated as of `0.9.0` please update your code to use the `TreeMesh.refine_surface` functionality. It will be removed in a future version of discretize.
Forward modeling data#
We can now create a magnetization model and generate data.
model_azm_dip = np.zeros((mesh.nC, 2))
model_amp = np.ones(mesh.nC) * 1e-8
ind = utils.model_builder.get_indices_block(
np.r_[-30, -20, -10],
np.r_[30, 20, 25],
mesh.gridCC,
)[0]
model_amp[ind] = 0.05
model_azm_dip[ind, 0] = 45.0
model_azm_dip[ind, 1] = 90.0
# Remove air cells
model_azm_dip = model_azm_dip[actv, :]
model_amp = model_amp[actv]
model = sdiag(model_amp) * utils.mat_utils.dip_azimuth2cartesian(
model_azm_dip[:, 0], model_azm_dip[:, 1]
)
# Create reduced identity map
idenMap = maps.IdentityMap(nP=nC * 3)
# Create the simulation
simulation = magnetics.simulation.Simulation3DIntegral(
survey=survey, mesh=mesh, chiMap=idenMap, active_cells=actv, model_type="vector"
)
# Compute some data and add some random noise
d = simulation.dpred(mkvc(model))
std = 10 # nT
synthetic_data = d + np.random.randn(len(d)) * std
wd = np.ones(len(d)) * std
# Assign data and uncertainties to the survey
data_object = data.Data(survey, dobs=synthetic_data, standard_deviation=wd)
# Create a projection matrix for plotting later
actv_plot = maps.InjectActiveCells(mesh, actv, np.nan)
Inversion#
We can now attempt the inverse calculations.
# Create sensitivity weights from our linear forward operator
rxLoc = survey.source_field.receiver_list[0].locations
# This Mapping connects the regularizations for the three-component
# vector model
wires = maps.Wires(("p", nC), ("s", nC), ("t", nC))
m0 = np.ones(3 * nC) * 1e-4 # Starting model
# Create the regularization on the amplitude of magnetization
reg = regularization.VectorAmplitude(
mesh,
mapping=idenMap,
active_cells=actv,
reference_model_in_smooth=True,
norms=[0.0, 2.0, 2.0, 2.0],
gradient_type="total",
)
# Data misfit function
dmis = data_misfit.L2DataMisfit(simulation=simulation, data=data_object)
dmis.W = 1.0 / data_object.standard_deviation
# The optimization scheme
opt = optimization.ProjectedGNCG(
maxIter=20, lower=-10, upper=10.0, maxIterLS=20, maxIterCG=20, tolCG=1e-4
)
# The inverse problem
invProb = inverse_problem.BaseInvProblem(dmis, reg, opt)
# Estimate the initial beta factor
betaest = directives.BetaEstimate_ByEig(beta0_ratio=1e1)
# Add sensitivity weights
sensitivity_weights = directives.UpdateSensitivityWeights()
# Here is where the norms are applied
IRLS = directives.UpdateIRLS(
f_min_change=1e-3, max_irls_iterations=10, misfit_tolerance=5e-1
)
# Pre-conditioner
update_Jacobi = directives.UpdatePreconditioner()
inv = inversion.BaseInversion(
invProb, directiveList=[sensitivity_weights, IRLS, update_Jacobi, betaest]
)
# Run the inversion
mrec = inv.run(m0)
/home/vsts/work/1/s/simpeg/optimization.py:1422: FutureWarning:
InexactCG.tolCG has been deprecated, please use InexactCG.cg_atol. It will be removed in version 0.26.0 of SimPEG.
/home/vsts/work/1/s/simpeg/optimization.py:943: FutureWarning:
InexactCG.maxIterCG has been deprecated, please use InexactCG.cg_maxiter. It will be removed in version 0.26.0 of SimPEG.
Running inversion with SimPEG v0.24.1.dev29+g256406f92
simpeg.InvProblem will set Regularization.reference_model to m0.
simpeg.InvProblem will set Regularization.reference_model to m0.
simpeg.InvProblem will set Regularization.reference_model to m0.
simpeg.InvProblem will set Regularization.reference_model to m0.
model has any nan: 0
================================================= Projected GNCG =================================================
# beta phi_d phi_m f |proj(x-g)-x| LS iterCG CG |Ax-b|/|b| CG |Ax-b| Comment
-----------------------------------------------------------------------------------------------------------------
x0 has any nan: 0
0 1.28e+05 1.19e+04 0.00e+00 1.19e+04 2.43e+03 0 0 inf inf
1 6.38e+04 6.86e+03 1.61e-02 7.89e+03 2.02e+03 0 20 2.10e-01 4.82e+04
2 3.19e+04 4.63e+03 4.05e-02 5.92e+03 1.85e+03 0 20 7.33e-06 6.39e-01
3 1.60e+04 2.67e+03 8.41e-02 4.01e+03 1.74e+03 0 20 2.89e-05 1.74e+00
4 7.98e+03 1.32e+03 1.44e-01 2.46e+03 1.61e+03 0 20 8.09e-05 3.57e+00
5 3.99e+03 5.84e+02 2.07e-01 1.41e+03 1.48e+03 0 20 4.22e-04 1.25e+01
Reached starting chifact with l2-norm regularization: Start IRLS steps...
irls_threshold 0.008951116313856873
6 3.99e+03 2.52e+02 4.69e-01 2.12e+03 1.55e+03 0 20 1.73e-03 3.18e+01
7 3.99e+03 4.65e+02 5.57e-01 2.69e+03 1.55e+03 0 20 3.06e-04 5.73e+00
8 2.62e+03 7.08e+02 6.87e-01 2.51e+03 1.25e+03 0 20 1.25e-05 2.43e-01
9 1.80e+03 6.39e+02 1.03e+00 2.50e+03 1.53e+03 0 20 6.08e-06 4.56e-02
10 1.27e+03 6.10e+02 1.29e+00 2.25e+03 1.57e+03 0 20 2.28e-06 3.77e-02
11 1.27e+03 5.68e+02 1.29e+00 2.21e+03 1.70e+03 0 20 8.50e-07 1.73e-02
12 8.68e+02 6.46e+02 1.11e+00 1.61e+03 9.35e+02 0 20 7.40e-08 2.22e-03
13 8.68e+02 5.28e+02 1.14e+00 1.52e+03 1.47e+03 0 20 8.87e-07 5.34e-03
14 8.68e+02 5.36e+02 1.04e+00 1.44e+03 1.37e+03 0 20 1.11e-07 1.25e-03
15 8.68e+02 5.34e+02 9.70e-01 1.38e+03 1.32e+03 0 20 7.78e-08 6.73e-04
Reach maximum number of IRLS cycles: 10
------------------------- STOP! -------------------------
1 : |fc-fOld| = 0.0000e+00 <= tolF*(1+|f0|) = 1.1929e+03
1 : |xc-x_last| = 2.9759e-02 <= tolX*(1+|x0|) = 1.0290e-01
0 : |proj(x-g)-x| = 1.3245e+03 <= tolG = 1.0000e-01
0 : |proj(x-g)-x| = 1.3245e+03 <= 1e3*eps = 1.0000e-02
0 : maxIter = 20 <= iter = 16
------------------------- DONE! -------------------------
Final Plot#
Let’s compare the smooth and compact model
plt.figure(figsize=(12, 6))
ax = plt.subplot(2, 2, 1)
im = utils.plot_utils.plot2Ddata(xyzLoc, synthetic_data, ax=ax)
plt.colorbar(im[0])
ax.set_title("Predicted data.")
plt.gca().set_aspect("equal", adjustable="box")
for ii, (title, mvec) in enumerate(
[("True model", model), ("Smooth model", invProb.l2model), ("Sparse model", mrec)]
):
ax = plt.subplot(2, 2, ii + 2)
mesh.plot_slice(
actv_plot * mvec.reshape((-1, 3), order="F"),
v_type="CCv",
view="vec",
ax=ax,
normal="Y",
grid=True,
quiver_opts={
"pivot": "mid",
"scale": 8 * np.abs(mvec).max(),
"scale_units": "inches",
},
)
ax.set_xlim([-200, 200])
ax.set_ylim([-100, 75])
ax.set_title(title)
ax.set_xlabel("x")
ax.set_ylabel("z")
plt.gca().set_aspect("equal", adjustable="box")
plt.show()
print("END")
# Plot the final predicted data and the residual
# plt.figure()
# ax = plt.subplot(1, 2, 1)
# utils.plot_utils.plot2Ddata(xyzLoc, invProb.dpred, ax=ax)
# ax.set_title("Predicted data.")
# plt.gca().set_aspect("equal", adjustable="box")
#
# ax = plt.subplot(1, 2, 2)
# utils.plot_utils.plot2Ddata(xyzLoc, synthetic_data - invProb.dpred, ax=ax)
# ax.set_title("Data residual.")
# plt.gca().set_aspect("equal", adjustable="box")

END
Total running time of the script: (0 minutes 24.393 seconds)
Estimated memory usage: 418 MB