.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "content/tutorials/08-tdem/plot_inv_1_em1dtm.py"
.. LINE NUMBERS ARE GIVEN BELOW.
.. only:: html
.. note::
:class: sphx-glr-download-link-note
Click :ref:`here <sphx_glr_download_content_tutorials_08-tdem_plot_inv_1_em1dtm.py>`
to download the full example code
.. rst-class:: sphx-glr-example-title
.. _sphx_glr_content_tutorials_08-tdem_plot_inv_1_em1dtm.py:
1D Inversion of Time-Domain Data for a Single Sounding
======================================================
Here we use the module *SimPEG.electromangetics.time_domain_1d* to invert
time domain data and recover a 1D electrical conductivity model.
In this tutorial, we focus on the following:
- How to define sources and receivers from a survey file
- How to define the survey
- Sparse 1D inversion of with iteratively re-weighted least-squares
For this tutorial, we will invert 1D time domain data for a single sounding.
The end product is layered Earth model which explains the data. The survey
consisted of a horizontal loop with a radius of 6 m, located 20 m above the
surface. The receiver measured the vertical component of the magnetic flux
at the loop's centre.
.. GENERATED FROM PYTHON SOURCE LINES 23-26
Import modules
--------------
.. GENERATED FROM PYTHON SOURCE LINES 26-53
.. code-block:: default
import os
import tarfile
import numpy as np
import matplotlib.pyplot as plt
from discretize import TensorMesh
import SimPEG.electromagnetics.time_domain as tdem
from SimPEG.utils import mkvc, plot_1d_layer_model
from SimPEG import (
maps,
data,
data_misfit,
inverse_problem,
regularization,
optimization,
directives,
inversion,
utils,
)
plt.rcParams.update({"font.size": 16, "lines.linewidth": 2, "lines.markersize": 8})
# sphinx_gallery_thumbnail_number = 2
.. GENERATED FROM PYTHON SOURCE LINES 54-62
Download Test Data File
-----------------------
Here we provide the file path to the data we plan on inverting.
The path to the data file is stored as a
tar-file on our google cloud bucket:
"https://storage.googleapis.com/simpeg/doc-assets/em1dtm.tar.gz"
.. GENERATED FROM PYTHON SOURCE LINES 62-80
.. code-block:: default
# storage bucket where we have the data
data_source = "https://storage.googleapis.com/simpeg/doc-assets/em1dtm.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
data_filename = dir_path + "em1dtm_data.txt"
.. rst-class:: sphx-glr-script-out
.. code-block:: none
Downloading https://storage.googleapis.com/simpeg/doc-assets/em1dtm.tar.gz
saved to: /home/vsts/work/1/s/tutorials/08-tdem/em1dtm.tar.gz
Download completed!
.. GENERATED FROM PYTHON SOURCE LINES 81-87
Load Data and Plot
------------------
Here we load and plot the 1D sounding data. In this case, we have the B-field
response to a step-off waveform.
.. GENERATED FROM PYTHON SOURCE LINES 87-102
.. code-block:: default
# Load field data
dobs = np.loadtxt(str(data_filename), skiprows=1)
times = dobs[:, 0]
dobs = mkvc(dobs[:, -1])
fig = plt.figure(figsize=(7, 7))
ax = fig.add_axes([0.15, 0.15, 0.8, 0.75])
ax.loglog(times, np.abs(dobs), "k-o", lw=3)
ax.set_xlabel("Times (s)")
ax.set_ylabel("|B| (T)")
ax.set_title("Observed Data")
.. image-sg:: /content/tutorials/08-tdem/images/sphx_glr_plot_inv_1_em1dtm_001.png
:alt: Observed Data
:srcset: /content/tutorials/08-tdem/images/sphx_glr_plot_inv_1_em1dtm_001.png
:class: sphx-glr-single-img
.. rst-class:: sphx-glr-script-out
.. code-block:: none
Text(0.5, 1.0, 'Observed Data')
.. GENERATED FROM PYTHON SOURCE LINES 103-110
Defining the Survey
-------------------
Here we demonstrate a general way to define the receivers, sources, waveforms and survey.
For this tutorial, we define a single horizontal loop source as well
a receiver which measures the vertical component of the magnetic flux.
.. GENERATED FROM PYTHON SOURCE LINES 110-147
.. code-block:: default
# Source loop geometry
source_location = np.array([0.0, 0.0, 20.0])
source_orientation = "z" # "x", "y" or "z"
source_current = 1.0 # peak current amplitude
source_radius = 6.0 # loop radius
# Receiver geometry
receiver_location = np.array([0.0, 0.0, 20.0])
receiver_orientation = "z" # "x", "y" or "z"
# Receiver list
receiver_list = []
receiver_list.append(
tdem.receivers.PointMagneticFluxDensity(
receiver_location, times, orientation=receiver_orientation
)
)
# Define the source waveform.
waveform = tdem.sources.StepOffWaveform()
# Sources
source_list = [
tdem.sources.CircularLoop(
receiver_list=receiver_list,
location=source_location,
waveform=waveform,
current=source_current,
radius=source_radius,
)
]
# Survey
survey = tdem.Survey(source_list)
.. GENERATED FROM PYTHON SOURCE LINES 148-154
Assign Uncertainties and Define the Data Object
-----------------------------------------------
Here is where we define the data that are inverted. The data are defined by
the survey, the observation values and the uncertainties.
.. GENERATED FROM PYTHON SOURCE LINES 154-162
.. code-block:: default
# 5% of the absolute value
uncertainties = 0.05 * np.abs(dobs) * np.ones(np.shape(dobs))
# Define the data object
data_object = data.Data(survey, dobs=dobs, standard_deviation=uncertainties)
.. GENERATED FROM PYTHON SOURCE LINES 163-169
Defining a 1D Layered Earth (1D Tensor Mesh)
--------------------------------------------
Here, we define the layer thicknesses for our 1D simulation. To do this, we use
the TensorMesh class.
.. GENERATED FROM PYTHON SOURCE LINES 169-177
.. code-block:: default
# Layer thicknesses
inv_thicknesses = np.logspace(0, 1.5, 25)
# Define a mesh for plotting and regularization.
mesh = TensorMesh([(np.r_[inv_thicknesses, inv_thicknesses[-1]])], "0")
.. GENERATED FROM PYTHON SOURCE LINES 178-190
Define a Starting and Reference Model
-------------------------------------
Here, we 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 starting model is log(0.1) S/m.
Define log-conductivity values for each layer since our model is the
log-conductivity. Don't make the values 0!
Otherwise the gradient for the 1st iteration is zero and the inversion will
not converge.
.. GENERATED FROM PYTHON SOURCE LINES 190-198
.. code-block:: default
# Define model. A resistivity (Ohm meters) or conductivity (S/m) for each layer.
starting_model = np.log(0.1 * np.ones(mesh.nC))
# Define mapping from model to active cells.
model_mapping = maps.ExpMap()
.. GENERATED FROM PYTHON SOURCE LINES 199-202
Define the Physics using a Simulation Object
--------------------------------------------
.. GENERATED FROM PYTHON SOURCE LINES 202-207
.. code-block:: default
simulation = tdem.Simulation1DLayered(
survey=survey, thicknesses=inv_thicknesses, sigmaMap=model_mapping
)
.. GENERATED FROM PYTHON SOURCE LINES 208-218
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 218-243
.. code-block:: default
# 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.
# The weighting is defined by the reciprocal of the uncertainties.
dmis = data_misfit.L2DataMisfit(simulation=simulation, data=data_object)
dmis.W = 1.0 / uncertainties
# Define the regularization (model objective function)
reg_map = maps.IdentityMap(nP=mesh.nC)
reg = regularization.Sparse(mesh, mapping=reg_map, alpha_s=0.01, alpha_x=1.0)
# set reference model
reg.mref = starting_model
# Define sparse and blocky norms p, q
reg.norms = [1, 0]
# Define how the optimization problem is solved. Here we will use an inexact
# Gauss-Newton approach that employs the conjugate gradient solver.
opt = optimization.ProjectedGNCG(maxIter=100, maxIterLS=20, maxIterCG=30, tolCG=1e-3)
# Define the inverse problem
inv_prob = inverse_problem.BaseInvProblem(dmis, reg, opt)
.. GENERATED FROM PYTHON SOURCE LINES 244-251
Define Inversion Directives
---------------------------
Here we define any directiveas 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 251-282
.. code-block:: default
# Defining a starting value for the trade-off parameter (beta) between the data
# misfit and the regularization.
starting_beta = directives.BetaEstimate_ByEig(beta0_ratio=1e1)
# Update the preconditionner
update_Jacobi = directives.UpdatePreconditioner()
# Options for outputting recovered models and predicted data for each beta.
save_iteration = directives.SaveOutputEveryIteration(save_txt=False)
# Directives for the IRLS
update_IRLS = directives.Update_IRLS(
max_irls_iterations=30, minGNiter=1, coolEpsFact=1.5, update_beta=True
)
# Updating the preconditionner if it is model dependent.
update_jacobi = directives.UpdatePreconditioner()
# Add sensitivity weights
sensitivity_weights = directives.UpdateSensitivityWeights()
# The directives are defined as a list.
directives_list = [
sensitivity_weights,
starting_beta,
save_iteration,
update_IRLS,
update_jacobi,
]
.. GENERATED FROM PYTHON SOURCE LINES 283-289
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 289-297
.. code-block:: default
# Here we combine the inverse problem and the set of directives
inv = inversion.BaseInversion(inv_prob, directives_list)
# Run the inversion
recovered_model = inv.run(starting_model)
.. rst-class:: sphx-glr-script-out
.. code-block:: none
SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.
***Done using same Solver, and solver_opts as the Simulation1DLayered 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.26e+03 2.00e+03 0.00e+00 2.00e+03 4.02e+02 0
1 6.30e+02 4.83e+02 3.48e-01 7.02e+02 4.66e+02 1
2 3.15e+02 4.27e+02 6.11e-01 6.19e+02 1.18e+03 0
3 1.57e+02 1.66e+01 3.48e-01 7.14e+01 1.13e+02 0
Reached starting chifact with l2-norm regularization: Start IRLS steps...
irls_threshold 2.3469232142123606
4 7.87e+01 3.53e+00 3.62e-01 3.20e+01 1.29e+01 0 Skip BFGS
5 4.47e+02 1.65e+00 4.44e-01 2.00e+02 1.05e+02 0 Skip BFGS
6 1.05e+03 5.73e+00 3.93e-01 4.20e+02 1.17e+02 0
7 8.03e+02 2.00e+01 3.66e-01 3.13e+02 2.80e+01 0
8 8.03e+02 1.42e+01 3.89e-01 3.27e+02 7.19e+00 0
9 8.03e+02 1.44e+01 3.95e-01 3.32e+02 8.39e+00 0
10 8.03e+02 1.47e+01 3.96e-01 3.33e+02 1.00e+01 0
Minimum decrease in regularization.End of IRLS
------------------------- STOP! -------------------------
1 : |fc-fOld| = 0.0000e+00 <= tolF*(1+|f0|) = 2.0047e+02
1 : |xc-x_last| = 2.0349e-01 <= tolX*(1+|x0|) = 1.2741e+00
0 : |proj(x-g)-x| = 1.0015e+01 <= tolG = 1.0000e-01
0 : |proj(x-g)-x| = 1.0015e+01 <= 1e3*eps = 1.0000e-02
0 : maxIter = 100 <= iter = 11
------------------------- DONE! -------------------------
.. GENERATED FROM PYTHON SOURCE LINES 298-300
Plotting Results
---------------------
.. GENERATED FROM PYTHON SOURCE LINES 300-347
.. code-block:: default
# Load the true model and layer thicknesses
true_model = np.array([0.1, 1.0, 0.1])
true_layers = np.r_[40.0, 40.0, 160.0]
# Extract Least-Squares model
l2_model = inv_prob.l2model
print(np.shape(l2_model))
# Plot true model and recovered model
fig = plt.figure(figsize=(8, 9))
x_min = np.min(
np.r_[model_mapping * recovered_model, model_mapping * l2_model, true_model]
)
x_max = np.max(
np.r_[model_mapping * recovered_model, model_mapping * l2_model, true_model]
)
ax1 = fig.add_axes([0.2, 0.15, 0.7, 0.7])
plot_1d_layer_model(true_layers, true_model, ax=ax1, show_layers=False, color="k")
plot_1d_layer_model(
mesh.h[0], model_mapping * l2_model, ax=ax1, show_layers=False, color="b"
)
plot_1d_layer_model(
mesh.h[0], model_mapping * recovered_model, ax=ax1, show_layers=False, color="r"
)
ax1.set_xlim(0.01, 10)
ax1.grid()
ax1.set_title("True and Recovered Models")
ax1.legend(["True Model", "L2-Model", "Sparse Model"])
plt.gca().invert_yaxis()
# Plot predicted and observed data
dpred_l2 = simulation.dpred(l2_model)
dpred_final = simulation.dpred(recovered_model)
fig = plt.figure(figsize=(7, 7))
ax1 = fig.add_axes([0.15, 0.15, 0.8, 0.75])
ax1.loglog(times, np.abs(dobs), "k-o")
ax1.loglog(times, np.abs(dpred_l2), "b-o")
ax1.loglog(times, np.abs(dpred_final), "r-o")
ax1.grid()
ax1.set_xlabel("times (Hz)")
ax1.set_ylabel("|Hs/Hp| (ppm)")
ax1.set_title("Predicted and Observed Data")
ax1.legend(["Observed", "L2-Model", "Sparse"], loc="upper right")
plt.show()
.. rst-class:: sphx-glr-horizontal
*
.. image-sg:: /content/tutorials/08-tdem/images/sphx_glr_plot_inv_1_em1dtm_002.png
:alt: True and Recovered Models
:srcset: /content/tutorials/08-tdem/images/sphx_glr_plot_inv_1_em1dtm_002.png
:class: sphx-glr-multi-img
*
.. image-sg:: /content/tutorials/08-tdem/images/sphx_glr_plot_inv_1_em1dtm_003.png
:alt: Predicted and Observed Data
:srcset: /content/tutorials/08-tdem/images/sphx_glr_plot_inv_1_em1dtm_003.png
:class: sphx-glr-multi-img
.. rst-class:: sphx-glr-script-out
.. code-block:: none
(26,)
.. rst-class:: sphx-glr-timing
**Total running time of the script:** ( 0 minutes 23.479 seconds)
**Estimated memory usage:** 18 MB
.. _sphx_glr_download_content_tutorials_08-tdem_plot_inv_1_em1dtm.py:
.. only:: html
.. container:: sphx-glr-footer sphx-glr-footer-example
.. container:: sphx-glr-download sphx-glr-download-python
:download:`Download Python source code: plot_inv_1_em1dtm.py <plot_inv_1_em1dtm.py>`
.. container:: sphx-glr-download sphx-glr-download-jupyter
:download:`Download Jupyter notebook: plot_inv_1_em1dtm.ipynb <plot_inv_1_em1dtm.ipynb>`
.. only:: html
.. rst-class:: sphx-glr-signature
`Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_