.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "content/user-guide/tutorials/04-magnetics/plot_inv_2a_magnetics_induced.py"
.. LINE NUMBERS ARE GIVEN BELOW.
.. only:: html
.. note::
:class: sphx-glr-download-link-note
:ref:`Go to the end <sphx_glr_download_content_user-guide_tutorials_04-magnetics_plot_inv_2a_magnetics_induced.py>`
to download the full example code.
.. rst-class:: sphx-glr-example-title
.. _sphx_glr_content_user-guide_tutorials_04-magnetics_plot_inv_2a_magnetics_induced.py:
Sparse Norm Inversion for Total Magnetic Intensity Data on a Tensor Mesh
========================================================================
Here we invert total magnetic intensity (TMI) data to recover a magnetic
susceptibility model. We formulate the inverse problem as an iteratively
re-weighted least-squares (IRLS) optimization problem. For this tutorial, we
focus on the following:
- Defining the survey from xyz formatted data
- Generating a mesh based on survey geometry
- Including surface topography
- Defining the inverse problem (data misfit, regularization, optimization)
- Specifying directives for the inversion
- Setting sparse and blocky norms
- Plotting the recovered model and data misfit
Although we consider TMI data in this tutorial, the same approach
can be used to invert other types of geophysical data.
.. GENERATED FROM PYTHON SOURCE LINES 25-28
Import modules
--------------
.. GENERATED FROM PYTHON SOURCE LINES 28-53
.. code-block:: Python
import os
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import tarfile
from discretize import TensorMesh
from discretize.utils import active_from_xyz
from simpeg.potential_fields import magnetics
from simpeg.utils import plot2Ddata, model_builder
from simpeg import (
maps,
data,
inverse_problem,
data_misfit,
regularization,
optimization,
directives,
inversion,
utils,
)
# sphinx_gallery_thumbnail_number = 3
.. GENERATED FROM PYTHON SOURCE LINES 54-63
Load Data and Plot
------------------
File paths for assets we are loading. To set up the inversion, we require
topography and field observations. The true model defined on the whole mesh
is loaded to compare with the inversion result. These files are stored as a
tar-file on our google cloud bucket:
"https://storage.googleapis.com/simpeg/doc-assets/magnetics.tar.gz"
.. GENERATED FROM PYTHON SOURCE LINES 63-83
.. code-block:: Python
# storage bucket where we have the data
data_source = "https://storage.googleapis.com/simpeg/doc-assets/magnetics.tar.gz"
# download the data
downloaded_data = utils.download(data_source, overwrite=True)
# unzip the tarfile
tar = tarfile.open(downloaded_data, "r")
tar.extractall()
tar.close()
# path to the directory containing our data
dir_path = downloaded_data.split(".")[0] + os.path.sep
# files to work with
topo_filename = dir_path + "magnetics_topo.txt"
data_filename = dir_path + "magnetics_data.obs"
.. rst-class:: sphx-glr-script-out
.. code-block:: none
Downloading https://storage.googleapis.com/simpeg/doc-assets/magnetics.tar.gz
saved to: /home/vsts/work/1/s/tutorials/04-magnetics/magnetics.tar.gz
Download completed!
.. GENERATED FROM PYTHON SOURCE LINES 84-91
Load Data and Plot
------------------
Here we load and plot synthetic TMI data. Topography is generally
defined as an (N, 3) array. TMI data is generally defined with 4 columns:
x, y, z and data.
.. GENERATED FROM PYTHON SOURCE LINES 91-124
.. code-block:: Python
topo_xyz = np.loadtxt(str(topo_filename))
dobs = np.loadtxt(str(data_filename))
receiver_locations = dobs[:, 0:3]
dobs = dobs[:, -1]
# Plot
fig = plt.figure(figsize=(6, 5))
v_max = np.max(np.abs(dobs))
ax1 = fig.add_axes([0.1, 0.1, 0.75, 0.85])
plot2Ddata(
receiver_locations,
dobs,
ax=ax1,
ncontour=30,
clim=(-v_max, v_max),
contourOpts={"cmap": "bwr"},
)
ax1.set_title("TMI Anomaly")
ax1.set_xlabel("x (m)")
ax1.set_ylabel("y (m)")
ax2 = fig.add_axes([0.85, 0.05, 0.05, 0.9])
norm = mpl.colors.Normalize(vmin=-np.max(np.abs(dobs)), vmax=np.max(np.abs(dobs)))
cbar = mpl.colorbar.ColorbarBase(
ax2, norm=norm, orientation="vertical", cmap=mpl.cm.bwr
)
cbar.set_label("$nT$", rotation=270, labelpad=15, size=12)
plt.show()
.. image-sg:: /content/user-guide/tutorials/04-magnetics/images/sphx_glr_plot_inv_2a_magnetics_induced_001.png
:alt: TMI Anomaly
:srcset: /content/user-guide/tutorials/04-magnetics/images/sphx_glr_plot_inv_2a_magnetics_induced_001.png
:class: sphx-glr-single-img
.. GENERATED FROM PYTHON SOURCE LINES 125-134
Assign Uncertainty
------------------
Inversion with SimPEG requires that we define standard deviation on our data.
This represents our estimate of the noise in our data. For magnetic inversions,
a constant floor value is generall applied to all data. For this tutorial, the
standard deviation on each datum will be 2% of the maximum observed magnetics
anomaly value.
.. GENERATED FROM PYTHON SOURCE LINES 134-139
.. code-block:: Python
maximum_anomaly = np.max(np.abs(dobs))
std = 0.02 * maximum_anomaly * np.ones(len(dobs))
.. GENERATED FROM PYTHON SOURCE LINES 140-148
Defining the Survey
-------------------
Here, we define survey that will be used for the simulation. Magnetic
surveys are simple to create. The user only needs an (N, 3) array to define
the xyz locations of the observation locations, the list of field components
which are to be modeled and the properties of the Earth's field.
.. GENERATED FROM PYTHON SOURCE LINES 148-174
.. code-block:: Python
# Define the component(s) of the field we are inverting as a list. Here we will
# invert total magnetic intensity data.
components = ["tmi"]
# Use the observation locations and components to define the receivers. To
# simulate data, the receivers must be defined as a list.
receiver_list = magnetics.receivers.Point(receiver_locations, components=components)
receiver_list = [receiver_list]
# Define the inducing field H0 = (intensity [nT], inclination [deg], declination [deg])
inclination = 90
declination = 0
strength = 50000
source_field = magnetics.sources.UniformBackgroundField(
receiver_list=receiver_list,
amplitude=strength,
inclination=inclination,
declination=declination,
)
# Define the survey
survey = magnetics.survey.Survey(source_field)
.. GENERATED FROM PYTHON SOURCE LINES 175-181
Defining the Data
-----------------
Here is where we define the data that is inverted. The data is defined by
the survey, the observation values and the standard deviations.
.. GENERATED FROM PYTHON SOURCE LINES 181-185
.. code-block:: Python
data_object = data.Data(survey, dobs=dobs, standard_deviation=std)
.. GENERATED FROM PYTHON SOURCE LINES 186-192
Defining a Tensor Mesh
----------------------
Here, we create the tensor mesh that will be used to invert TMI data.
If desired, we could define an OcTree mesh.
.. GENERATED FROM PYTHON SOURCE LINES 192-199
.. code-block:: Python
dh = 5.0
hx = [(dh, 5, -1.3), (dh, 40), (dh, 5, 1.3)]
hy = [(dh, 5, -1.3), (dh, 40), (dh, 5, 1.3)]
hz = [(dh, 5, -1.3), (dh, 15)]
mesh = TensorMesh([hx, hy, hz], "CCN")
.. GENERATED FROM PYTHON SOURCE LINES 200-208
Starting/Reference Model and Mapping on Tensor Mesh
---------------------------------------------------
Here, we would create starting and/or reference models for the inversion as
well as the mapping from the model space to the active cells. Starting and
reference models can be a constant background value or contain a-priori
structures. Here, the background is 1e-4 SI.
.. GENERATED FROM PYTHON SOURCE LINES 208-224
.. code-block:: Python
# Define background susceptibility model in SI. Don't make this 0!
# Otherwise the gradient for the 1st iteration is zero and the inversion will
# not converge.
background_susceptibility = 1e-4
# Find the indecies of the active cells in forward model (ones below surface)
active_cells = active_from_xyz(mesh, topo_xyz)
# Define mapping from model to active cells
nC = int(active_cells.sum())
model_map = maps.IdentityMap(nP=nC) # model consists of a value for each cell
# Define starting model
starting_model = background_susceptibility * np.ones(nC)
.. GENERATED FROM PYTHON SOURCE LINES 225-231
Define the Physics
------------------
Here, we define the physics of the magnetics problem by using the simulation
class.
.. GENERATED FROM PYTHON SOURCE LINES 233-234
Define the problem. Define the cells below topography and the mapping
.. GENERATED FROM PYTHON SOURCE LINES 234-244
.. code-block:: Python
simulation = magnetics.simulation.Simulation3DIntegral(
survey=survey,
mesh=mesh,
model_type="scalar",
chiMap=model_map,
active_cells=active_cells,
engine="choclo",
)
.. GENERATED FROM PYTHON SOURCE LINES 245-252
.. tip::
Since SimPEG v0.22.0 we can use `Choclo
<https://www.fatiando.org/choclo>`_ as the engine for running the magnetic
simulations, which results in faster and more memory efficient runs. Just
pass ``engine="choclo"`` when constructing the simulation.
.. GENERATED FROM PYTHON SOURCE LINES 255-264
Define Inverse Problem
----------------------
The inverse problem is defined by 3 things:
1) Data Misfit: a measure of how well our recovered model explains the field data
2) Regularization: constraints placed on the recovered model and a priori information
3) Optimization: the numerical approach used to solve the inverse problem
.. GENERATED FROM PYTHON SOURCE LINES 264-292
.. code-block:: Python
# Define the data misfit. Here the data misfit is the L2 norm of the weighted
# residual between the observed data and the data predicted for a given model.
# Within the data misfit, the residual between predicted and observed data are
# normalized by the data's standard deviation.
dmis = data_misfit.L2DataMisfit(data=data_object, simulation=simulation)
# Define the regularization (model objective function)
reg = regularization.Sparse(
mesh,
active_cells=active_cells,
mapping=model_map,
reference_model=starting_model,
gradient_type="total",
)
# Define sparse and blocky norms p, qx, qy, qz
reg.norms = [0, 0, 0, 0]
# Define how the optimization problem is solved. Here we will use a projected
# Gauss-Newton approach that employs the conjugate gradient solver.
opt = optimization.ProjectedGNCG(
maxIter=20, lower=0.0, upper=1.0, maxIterLS=20, maxIterCG=10, tolCG=1e-3
)
# Here we define the inverse problem that is to be solved
inv_prob = inverse_problem.BaseInvProblem(dmis, reg, opt)
.. GENERATED FROM PYTHON SOURCE LINES 293-300
Define Inversion Directives
---------------------------
Here we define any directives that are carried out during the inversion. This
includes the cooling schedule for the trade-off parameter (beta), stopping
criteria for the inversion and saving inversion results at each iteration.
.. GENERATED FROM PYTHON SOURCE LINES 300-335
.. code-block:: Python
# Defining a starting value for the trade-off parameter (beta) between the data
# misfit and the regularization.
starting_beta = directives.BetaEstimate_ByEig(beta0_ratio=5)
# Options for outputting recovered models and predicted data for each beta.
save_iteration = directives.SaveOutputEveryIteration(save_txt=False)
# Defines the directives for the IRLS regularization. This includes setting
# the cooling schedule for the trade-off parameter.
update_IRLS = directives.UpdateIRLS(
f_min_change=1e-4,
max_irls_iterations=30,
cooling_factor=1.5,
misfit_tolerance=1e-2,
)
# Updating the preconditioner if it is model dependent.
update_jacobi = directives.UpdatePreconditioner()
# Setting a stopping criteria for the inversion.
target_misfit = directives.TargetMisfit(chifact=1)
# Add sensitivity weights
sensitivity_weights = directives.UpdateSensitivityWeights(every_iteration=False)
# The directives are defined as a list.
directives_list = [
sensitivity_weights,
starting_beta,
save_iteration,
update_IRLS,
update_jacobi,
]
.. GENERATED FROM PYTHON SOURCE LINES 336-342
Running the Inversion
---------------------
To define the inversion object, we need to define the inversion problem and
the set of directives. We can then run the inversion.
.. GENERATED FROM PYTHON SOURCE LINES 342-352
.. code-block:: Python
# Here we combine the inverse problem and the set of directives
inv = inversion.BaseInversion(inv_prob, directives_list)
# Print target misfit to compare with convergence
# print("Target misfit is " + str(target_misfit.target))
# Run the inversion
recovered_model = inv.run(starting_model)
.. rst-class:: sphx-glr-script-out
.. code-block:: none
Running inversion with SimPEG v0.23.1.dev28+g33c4221d8
simpeg.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.
***Done using the default solver Pardiso and no solver_opts.***
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.97e+05 2.10e+04 0.00e+00 2.10e+04 1.39e+02 0
1 1.31e+05 4.80e+02 1.01e-02 1.80e+03 1.30e+02 0
Reached starting chifact with l2-norm regularization: Start IRLS steps...
irls_threshold 0.0013525622846263725
2 1.31e+05 2.59e+02 1.77e-02 2.59e+03 1.90e+02 0 Skip BFGS
3 7.89e+04 3.83e+02 2.11e-02 2.05e+03 1.17e+02 0
4 9.63e+04 2.57e+02 2.79e-02 2.95e+03 1.89e+02 0
5 5.32e+04 4.69e+02 2.98e-02 2.05e+03 1.14e+02 0
6 6.89e+04 2.29e+02 3.87e-02 2.89e+03 1.86e+02 0
7 4.13e+04 3.87e+02 3.98e-02 2.03e+03 1.17e+02 0
8 6.18e+04 1.69e+02 4.24e-02 2.79e+03 1.80e+02 0
9 7.99e+04 2.31e+02 3.06e-02 2.68e+03 1.69e+02 0
10 1.01e+05 2.42e+02 2.25e-02 2.50e+03 1.59e+02 0
11 1.23e+05 2.56e+02 1.80e-02 2.47e+03 1.39e+02 0
12 1.46e+05 2.71e+02 1.44e-02 2.38e+03 1.32e+02 0
13 1.75e+05 2.65e+02 1.18e-02 2.33e+03 1.40e+02 0
14 2.10e+05 2.64e+02 9.91e-03 2.35e+03 1.52e+02 0
15 2.54e+05 2.61e+02 8.19e-03 2.34e+03 1.50e+02 0
16 3.08e+05 2.60e+02 7.04e-03 2.43e+03 1.52e+02 0
17 3.68e+05 2.66e+02 5.71e-03 2.37e+03 1.35e+02 0
18 4.32e+05 2.76e+02 5.02e-03 2.45e+03 1.33e+02 0
19 5.00e+05 2.83e+02 4.39e-03 2.48e+03 1.32e+02 0
20 5.85e+05 2.78e+02 3.53e-03 2.34e+03 1.31e+02 0
------------------------- STOP! -------------------------
1 : |fc-fOld| = 1.3683e+02 <= tolF*(1+|f0|) = 2.1023e+03
1 : |xc-x_last| = 1.7318e-02 <= tolX*(1+|x0|) = 1.0219e-01
0 : |proj(x-g)-x| = 1.3132e+02 <= tolG = 1.0000e-01
0 : |proj(x-g)-x| = 1.3132e+02 <= 1e3*eps = 1.0000e-02
1 : maxIter = 20 <= iter = 20
------------------------- DONE! -------------------------
.. GENERATED FROM PYTHON SOURCE LINES 353-356
Recreate True Model
-------------------
.. GENERATED FROM PYTHON SOURCE LINES 356-369
.. code-block:: Python
background_susceptibility = 0.0001
sphere_susceptibility = 0.01
true_model = background_susceptibility * np.ones(nC)
ind_sphere = model_builder.get_indices_sphere(
np.r_[0.0, 0.0, -45.0], 15.0, mesh.cell_centers
)
ind_sphere = ind_sphere[active_cells]
true_model[ind_sphere] = sphere_susceptibility
.. GENERATED FROM PYTHON SOURCE LINES 370-373
Plotting True Model and Recovered Model
---------------------------------------
.. GENERATED FROM PYTHON SOURCE LINES 373-424
.. code-block:: Python
# Plot True Model
fig = plt.figure(figsize=(9, 4))
plotting_map = maps.InjectActiveCells(mesh, active_cells, np.nan)
ax1 = fig.add_axes([0.08, 0.1, 0.75, 0.8])
mesh.plot_slice(
plotting_map * true_model,
normal="Y",
ax=ax1,
ind=int(mesh.shape_cells[1] / 2),
grid=True,
clim=(np.min(true_model), np.max(true_model)),
pcolor_opts={"cmap": "viridis"},
)
ax1.set_title("Model slice at y = 0 m")
ax2 = fig.add_axes([0.85, 0.1, 0.05, 0.8])
norm = mpl.colors.Normalize(vmin=np.min(true_model), vmax=np.max(true_model))
cbar = mpl.colorbar.ColorbarBase(
ax2, norm=norm, orientation="vertical", cmap=mpl.cm.viridis, format="%.1e"
)
cbar.set_label("SI", rotation=270, labelpad=15, size=12)
plt.show()
# Plot Recovered Model
fig = plt.figure(figsize=(9, 4))
plotting_map = maps.InjectActiveCells(mesh, active_cells, np.nan)
ax1 = fig.add_axes([0.08, 0.1, 0.75, 0.8])
mesh.plot_slice(
plotting_map * recovered_model,
normal="Y",
ax=ax1,
ind=int(mesh.shape_cells[1] / 2),
grid=True,
clim=(np.min(recovered_model), np.max(recovered_model)),
pcolor_opts={"cmap": "viridis"},
)
ax1.set_title("Model slice at y = 0 m")
ax2 = fig.add_axes([0.85, 0.1, 0.05, 0.8])
norm = mpl.colors.Normalize(vmin=np.min(recovered_model), vmax=np.max(recovered_model))
cbar = mpl.colorbar.ColorbarBase(
ax2, norm=norm, orientation="vertical", cmap=mpl.cm.viridis, format="%.1e"
)
cbar.set_label("SI", rotation=270, labelpad=15, size=12)
plt.show()
.. rst-class:: sphx-glr-horizontal
*
.. image-sg:: /content/user-guide/tutorials/04-magnetics/images/sphx_glr_plot_inv_2a_magnetics_induced_002.png
:alt: Model slice at y = 0 m
:srcset: /content/user-guide/tutorials/04-magnetics/images/sphx_glr_plot_inv_2a_magnetics_induced_002.png
:class: sphx-glr-multi-img
*
.. image-sg:: /content/user-guide/tutorials/04-magnetics/images/sphx_glr_plot_inv_2a_magnetics_induced_003.png
:alt: Model slice at y = 0 m
:srcset: /content/user-guide/tutorials/04-magnetics/images/sphx_glr_plot_inv_2a_magnetics_induced_003.png
:class: sphx-glr-multi-img
.. GENERATED FROM PYTHON SOURCE LINES 425-428
Plotting Predicted Data and Misfit
----------------------------------
.. GENERATED FROM PYTHON SOURCE LINES 428-468
.. code-block:: Python
# Predicted data with final recovered model
dpred = inv_prob.dpred
# Observed data | Predicted data | Normalized data misfit
data_array = np.c_[dobs, dpred, (dobs - dpred) / std]
fig = plt.figure(figsize=(17, 4))
plot_title = ["Observed", "Predicted", "Normalized Misfit"]
plot_units = ["nT", "nT", ""]
ax1 = 3 * [None]
ax2 = 3 * [None]
norm = 3 * [None]
cbar = 3 * [None]
cplot = 3 * [None]
v_lim = [np.max(np.abs(dobs)), np.max(np.abs(dobs)), np.max(np.abs(data_array[:, 2]))]
for ii in range(0, 3):
ax1[ii] = fig.add_axes([0.33 * ii + 0.03, 0.11, 0.25, 0.84])
cplot[ii] = plot2Ddata(
receiver_list[0].locations,
data_array[:, ii],
ax=ax1[ii],
ncontour=30,
clim=(-v_lim[ii], v_lim[ii]),
contourOpts={"cmap": "bwr"},
)
ax1[ii].set_title(plot_title[ii])
ax1[ii].set_xlabel("x (m)")
ax1[ii].set_ylabel("y (m)")
ax2[ii] = fig.add_axes([0.33 * ii + 0.27, 0.11, 0.01, 0.84])
norm[ii] = mpl.colors.Normalize(vmin=-v_lim[ii], vmax=v_lim[ii])
cbar[ii] = mpl.colorbar.ColorbarBase(
ax2[ii], norm=norm[ii], orientation="vertical", cmap=mpl.cm.bwr
)
cbar[ii].set_label(plot_units[ii], rotation=270, labelpad=15, size=12)
plt.show()
.. image-sg:: /content/user-guide/tutorials/04-magnetics/images/sphx_glr_plot_inv_2a_magnetics_induced_004.png
:alt: Observed, Predicted, Normalized Misfit
:srcset: /content/user-guide/tutorials/04-magnetics/images/sphx_glr_plot_inv_2a_magnetics_induced_004.png
:class: sphx-glr-single-img
.. rst-class:: sphx-glr-timing
**Total running time of the script:** (0 minutes 32.526 seconds)
**Estimated memory usage:** 288 MB
.. _sphx_glr_download_content_user-guide_tutorials_04-magnetics_plot_inv_2a_magnetics_induced.py:
.. only:: html
.. container:: sphx-glr-footer sphx-glr-footer-example
.. container:: sphx-glr-download sphx-glr-download-jupyter
:download:`Download Jupyter notebook: plot_inv_2a_magnetics_induced.ipynb <plot_inv_2a_magnetics_induced.ipynb>`
.. container:: sphx-glr-download sphx-glr-download-python
:download:`Download Python source code: plot_inv_2a_magnetics_induced.py <plot_inv_2a_magnetics_induced.py>`
.. container:: sphx-glr-download sphx-glr-download-zip
:download:`Download zipped: plot_inv_2a_magnetics_induced.zip <plot_inv_2a_magnetics_induced.zip>`
.. only:: html
.. rst-class:: sphx-glr-signature
`Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_