1D Inversion of for a Single Sounding#

Here we use the module simpeg.electromangetics.frequency_domain_1d to invert frequency 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 frequency domain data for a single sounding. The end product is layered Earth model which explains the data. The survey consisted of a vertical magnetic dipole source located 30 m above the surface. The receiver measured the vertical component of the secondary field at a 10 m offset from the source in ppm.

Import modules#

import os
import tarfile
import numpy as np
import matplotlib.pyplot as plt

from discretize import TensorMesh

import simpeg.electromagnetics.frequency_domain as fdem
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

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/em1dfm.tar.gz”

# storage bucket where we have the data
data_source = "https://storage.googleapis.com/simpeg/doc-assets/em1dfm.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 + "em1dfm_data.txt"

Downloading https://storage.googleapis.com/simpeg/doc-assets/em1dfm.tar.gz
   saved to: /home/vsts/work/1/s/tutorials/07-fdem/em1dfm.tar.gz
Download completed!

Load Data and Plot#

Here we load and plot the 1D sounding data. In this case, we have the secondary field response in ppm for a set of frequencies.

# Load field data
# dobs = np.loadtxt(str(data_filename))
dobs = np.loadtxt(str(data_filename), skiprows=1)

# Define receiver locations and observed data
frequencies = dobs[:, 0]
dobs = mkvc(dobs[:, 1:].T)

fig, ax = plt.subplots(1, 1, figsize=(7, 7))
ax.loglog(frequencies, np.abs(dobs[0::2]), "k-o", lw=3)
ax.loglog(frequencies, np.abs(dobs[1::2]), "k:o", lw=3)
ax.set_xlabel("Frequency (Hz)")
ax.set_ylabel("|Hs/Hp| (ppm)")
ax.set_title("Magnetic Field as a Function of Frequency")
ax.legend(["Real", "Imaginary"])

Magnetic Field as a Function of Frequency

<matplotlib.legend.Legend object at 0x7ff8ee7481d0>

Defining the Survey#

Here we demonstrate a general way to define the receivers, sources and survey. The survey consisted of a vertical magnetic dipole source located 30 m above the surface. The receiver measured the vertical component of the secondary field at a 10 m offset from the source in ppm.

source_location = np.array([0.0, 0.0, 30.0])
moment = 1.0

receiver_location = np.array([10.0, 0.0, 30.0])
receiver_orientation = "z"
data_type = "ppm"

# Receiver list
receiver_list = []
receiver_list.append(
    fdem.receivers.PointMagneticFieldSecondary(
        receiver_location,
        orientation=receiver_orientation,
        data_type=data_type,
        component="real",
    )
)
receiver_list.append(
    fdem.receivers.PointMagneticFieldSecondary(
        receiver_location,
        orientation=receiver_orientation,
        data_type=data_type,
        component="imag",
    )
)

# Define source list
source_list = []
for freq in frequencies:
    source_list.append(
        fdem.sources.MagDipole(
            receiver_list=receiver_list,
            frequency=freq,
            location=source_location,
            orientation="z",
            moment=moment,
        )
    )

# Survey
survey = fdem.survey.Survey(source_list)

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.

# 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, noise_floor=uncertainties)

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.

# 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")

Define a Starting and/or Reference Model and the Mapping#

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.

# 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()

Define the Physics using a Simulation Object#

simulation = fdem.Simulation1DLayered(
    survey=survey, thicknesses=inv_thicknesses, sigmaMap=model_mapping
)

Define Inverse Problem#

The inverse problem is defined by 3 things:

Data Misfit: a measure of how well our recovered model explains the field data

Regularization: constraints placed on the recovered model and a priori information

Optimization: the numerical approach used to solve the inverse problem

# 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)

# Define the regularization (model objective function)
reg_map = maps.IdentityMap(nP=mesh.nC)
reg = regularization.Sparse(mesh, mapping=reg_map, alpha_s=0.025, alpha_x=1.0)

# reference model
reg.reference_model = starting_model

# Define sparse and blocky norms p, q
reg.norms = [0, 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=50, maxIterLS=20, cg_maxiter=30, cg_rtol=1e-3)

# Define the inverse problem
inv_prob = inverse_problem.BaseInvProblem(dmis, reg, opt)

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.

# 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)

# Directive for the IRLS
update_IRLS = directives.UpdateIRLS(max_irls_iterations=30, irls_cooling_factor=1.5)

# 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,
]

/home/vsts/work/1/s/simpeg/directives/_directives.py:1865: FutureWarning:

SaveEveryIteration.save_txt has been deprecated, please use SaveEveryIteration.on_disk. It will be removed in version 0.26.0 of SimPEG.

/home/vsts/work/1/s/simpeg/directives/_directives.py:1866: FutureWarning:

SaveEveryIteration.save_txt has been deprecated, please use SaveEveryIteration.on_disk. It will be removed in version 0.26.0 of SimPEG.

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.

# 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)

Running inversion with SimPEG v0.25.1.dev6+gb840df108
================================================= Projected GNCG =================================================
  #     beta     phi_d     phi_m       f      |proj(x-g)-x|  LS   iter_CG   CG |Ax-b|/|b|  CG |Ax-b|   Comment
-----------------------------------------------------------------------------------------------------------------
2.65e+02  2.81e+02  0.00e+00  2.81e+02                         0           inf          inf
2.65e+02  4.45e+01  2.07e-01  9.94e+01    7.91e+01      0      8        6.10e-04     4.82e-02
1.33e+02  8.35e+00  2.05e-01  3.55e+01    3.68e+01      0      8        6.35e-04     2.34e-02
Reached starting chifact with l2-norm regularization: Start IRLS steps...
irls_threshold 2.7179393704559085
1.33e+02  7.90e+00  1.98e-01  3.42e+01    8.19e+00      0      8        8.86e-04     7.26e-03
2.17e+02  1.32e+01  1.54e-01  4.66e+01    1.87e+01      0      7        4.89e-04     9.14e-03
1.63e+02  1.43e+01  1.60e-01  4.03e+01    6.51e+00      0      7        5.83e-04     3.79e-03
1.18e+02  1.36e+01  1.80e-01  3.48e+01    3.48e+00      0      7        6.88e-04     2.40e-03
8.75e+01  1.24e+01  1.93e-01  2.92e+01    3.57e+00      0      8        5.28e-04     1.89e-03
6.82e+01  1.08e+01  1.86e-01  2.35e+01    3.70e+00      0      8        2.86e-04     1.06e-03
6.82e+01  9.61e+00  1.54e-01  2.01e+01    5.74e+00      0      8        9.22e-04     5.29e-03
6.82e+01  8.64e+00  1.29e-01  1.75e+01    6.47e+00      0      9        4.45e-04     2.88e-03
1.08e+02  9.08e+00  8.82e-02  1.86e+01    1.87e+01      0      9        9.44e-04     1.77e-02
1.08e+02  8.80e+00  6.67e-02  1.60e+01    1.26e+01      0      10       4.30e-04     5.42e-03
1.69e+02  9.02e+00  4.31e-02  1.63e+01    3.22e+01      0      9        7.02e-04     2.26e-02
1.69e+02  8.47e+00  3.23e-02  1.39e+01    2.03e+01      0      10       3.03e-04     6.13e-03
2.69e+02  8.42e+00  2.24e-02  1.44e+01    3.78e+01      0      9        4.90e-04     1.85e-02
4.28e+02  8.48e+00  1.52e-02  1.50e+01    3.82e+01      0      9        5.31e-04     2.03e-02
6.81e+02  8.57e+00  1.01e-02  1.55e+01    3.77e+01      0      9        2.43e-04     9.18e-03
1.08e+03  8.66e+00  6.70e-03  1.59e+01    7.12e+01      0      9        3.13e-04     2.23e-02
1.70e+03  8.75e+00  4.41e-03  1.63e+01    4.87e+01      0      9        4.39e-04     2.14e-02
2.67e+03  8.85e+00  2.90e-03  1.66e+01    3.24e+01      0      9        5.57e-04     1.81e-02
4.18e+03  8.94e+00  1.91e-03  1.69e+01    3.22e+01      0      9        5.10e-04     1.64e-02
6.52e+03  9.04e+00  1.26e-03  1.72e+01    3.21e+01      0      9        5.34e-04     1.71e-02
6.52e+03  8.60e+00  8.90e-04  1.44e+01    1.25e+01      0      10       1.28e-04     1.60e-03
1.03e+04  8.51e+00  6.01e-04  1.47e+01    3.47e+01      0      9        8.93e-04     3.10e-02
1.64e+04  8.52e+00  4.00e-04  1.51e+01    3.33e+01      0      9        7.81e-04     2.60e-02
2.60e+04  8.58e+00  2.65e-04  1.55e+01    3.29e+01      0      9        7.54e-04     2.48e-02
4.11e+04  8.65e+00  1.75e-04  1.58e+01    3.27e+01      0      9        7.23e-04     2.36e-02
6.49e+04  8.74e+00  1.15e-04  1.62e+01    3.25e+01      0      9        6.69e-04     2.17e-02
1.02e+05  8.83e+00  7.57e-05  1.66e+01    3.23e+01      0      9        5.85e-04     1.89e-02
1.60e+05  8.93e+00  4.98e-05  1.69e+01    3.22e+01      0      9        5.20e-04     1.67e-02
2.49e+05  9.02e+00  3.27e-05  1.72e+01    3.21e+01      0      9        5.28e-04     1.69e-02
2.49e+05  8.59e+00  2.32e-05  1.44e+01    1.25e+01      0      10       1.30e-04     1.62e-03
Reach maximum number of IRLS cycles: 30
------------------------- STOP! -------------------------
: |fc-fOld| = 8.1682e-02 <= tolF*(1+|f0|) = 2.8196e+01
: |xc-x_last| = 1.2955e-01 <= tolX*(1+|x0|) = 1.2741e+00
: |proj(x-g)-x|    = 1.2459e+01 <= tolG          = 1.0000e-01
: |proj(x-g)-x|    = 1.2459e+01 <= 1e3*eps       = 1.0000e-02
: maxIter   =      50    <= iter          =     32
------------------------- DONE! -------------------------

Plotting Results#

# Load the true model and layer thicknesses
true_model = np.array([0.1, 1.0, 0.1])
true_layers = np.r_[20.0, 40.0, 160.0]

# Extract Least-Squares model
l2_model = inv_prob.l2model

# 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.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=(11, 6))
ax1 = fig.add_axes([0.2, 0.1, 0.6, 0.8])
ax1.loglog(frequencies, np.abs(dobs[0::2]), "k-o")
ax1.loglog(frequencies, np.abs(dobs[1::2]), "k:o")
ax1.loglog(frequencies, np.abs(dpred_l2[0::2]), "b-o")
ax1.loglog(frequencies, np.abs(dpred_l2[1::2]), "b:o")
ax1.loglog(frequencies, np.abs(dpred_final[0::2]), "r-o")
ax1.loglog(frequencies, np.abs(dpred_final[1::2]), "r:o")
ax1.set_xlabel("Frequencies (Hz)")
ax1.set_ylabel("|Hs/Hp| (ppm)")
ax1.set_title("Predicted and Observed Data")
ax1.legend(
    [
        "Observed (real)",
        "Observed (imag)",
        "L2-Model (real)",
        "L2-Model (imag)",
        "Sparse (real)",
        "Sparse (imag)",
    ],
    loc="upper left",
)
plt.show()