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.11/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,
)
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)
/home/vsts/work/1/s/examples/03-magnetics/plot_inv_mag_MVI_VectorAmplitude.py:105: BreakingChangeWarning:
Since SimPEG v0.25.0, the 'get_indices_block' function returns a single array with the cell indices, instead of a tuple with a single element. This means that we don't need to unpack the tuple anymore to access to the cell indices.
If you were using this function as in:
ind = get_indices_block(p0, p1, mesh.cell_centers)[0]
Make sure you update it to:
ind = get_indices_block(p0, p1, mesh.cell_centers)
To hide this warning, add this to your script or notebook:
import warnings
from simpeg.utils import BreakingChangeWarning
warnings.filterwarnings(action='ignore', category=BreakingChangeWarning)
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, cg_maxiter=20, cg_rtol=1e-3
)
# 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)
Running inversion with SimPEG v0.25.0
================================================= Projected GNCG =================================================
# beta phi_d phi_m f |proj(x-g)-x| LS iter_CG CG |Ax-b|/|b| CG |Ax-b| Comment
-----------------------------------------------------------------------------------------------------------------
0 1.34e+05 1.19e+04 0.00e+00 1.19e+04 0 inf inf
1 1.34e+05 7.00e+03 1.50e-02 9.01e+03 2.43e+03 0 20 2.13e-01 4.89e+04
2 6.69e+04 4.77e+03 3.83e-02 7.33e+03 2.02e+03 0 13 5.87e-04 5.19e+01
3 3.35e+04 2.79e+03 8.04e-02 5.48e+03 1.86e+03 0 14 6.19e-04 3.79e+01
4 1.67e+04 1.39e+03 1.39e-01 3.72e+03 1.74e+03 0 15 9.53e-04 4.31e+01
5 8.36e+03 6.18e+02 2.03e-01 2.32e+03 1.62e+03 0 18 8.18e-04 2.49e+01
6 4.18e+03 2.67e+02 2.61e-01 1.36e+03 1.49e+03 0 20 1.58e-03 3.01e+01
Reached starting chifact with l2-norm regularization: Start IRLS steps...
irls_threshold 0.008897945940906758
7 4.18e+03 4.93e+02 3.71e-01 2.04e+03 1.56e+03 0 18 5.24e-04 1.02e+01
8 4.18e+03 7.49e+02 4.52e-01 2.64e+03 1.56e+03 0 13 9.41e-04 1.89e+01
9 2.69e+03 6.57e+02 6.57e-01 2.42e+03 1.20e+03 0 14 5.10e-04 3.57e+00
10 1.83e+03 6.17e+02 8.92e-01 2.24e+03 1.53e+03 0 12 9.01e-04 1.47e+01
11 1.28e+03 5.71e+02 1.08e+00 1.95e+03 1.57e+03 0 11 9.32e-04 1.90e+01
12 1.28e+03 6.49e+02 1.06e+00 2.00e+03 1.70e+03 0 9 9.21e-04 2.79e+01
13 8.72e+02 5.29e+02 1.15e+00 1.53e+03 9.50e+02 0 11 9.02e-04 5.47e+00
14 8.72e+02 5.37e+02 1.07e+00 1.47e+03 1.47e+03 0 10 4.60e-04 5.17e+00
15 8.72e+02 5.35e+02 1.00e+00 1.41e+03 1.37e+03 0 10 3.35e-04 2.90e+00
16 8.72e+02 5.32e+02 9.51e-01 1.36e+03 1.33e+03 0 9 9.80e-04 7.17e+00
Reach maximum number of IRLS cycles: 10
------------------------- STOP! -------------------------
1 : |fc-fOld| = 1.6611e+01 <= tolF*(1+|f0|) = 1.1929e+03
1 : |xc-x_last| = 2.9300e-02 <= tolX*(1+|x0|) = 1.0290e-01
0 : |proj(x-g)-x| = 1.3253e+03 <= tolG = 1.0000e-01
0 : |proj(x-g)-x| = 1.3253e+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.643 seconds)
Estimated memory usage: 448 MB