Parametric 1D Inversion of Sounding Data#

Here we use the module simpeg.electromangetics.static.resistivity to invert DC resistivity sounding data and recover the resistivities and layer thicknesses for a 1D layered Earth. In this tutorial, we focus on the following:

How to define sources and receivers from a survey file

How to define the survey

Defining a model that consists of resistivities and layer thicknesses

For this tutorial, we will invert sounding data collected over a layered Earth using a Wenner array. The end product is layered Earth model which explains the data.

Import modules#

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

from discretize import TensorMesh

from simpeg import (
    maps,
    data,
    data_misfit,
    regularization,
    optimization,
    inverse_problem,
    inversion,
    directives,
    utils,
)
from simpeg.electromagnetics.static import resistivity as dc
from simpeg.utils import plot_1d_layer_model

mpl.rcParams.update({"font.size": 16})

# sphinx_gallery_thumbnail_number = 2

Define File Names#

Here we provide the file paths to assets we need to run the inversion. The Path to the true model is also provided for comparison with the inversion results. These files are stored as a tar-file on our google cloud bucket: “https://storage.googleapis.com/simpeg/doc-assets/dcr1d.tar.gz”

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

overwriting /home/vsts/work/1/s/tutorials/05-dcr/dcr1d.tar.gz
Downloading https://storage.googleapis.com/simpeg/doc-assets/dcr1d.tar.gz
   saved to: /home/vsts/work/1/s/tutorials/05-dcr/dcr1d.tar.gz
Download completed!

Load Data, Define Survey and Plot#

Here we load the observed data, define the DC survey geometry and plot the data values.

# Load data
dobs = np.loadtxt(str(data_filename))

A_electrodes = dobs[:, 0:3]
B_electrodes = dobs[:, 3:6]
M_electrodes = dobs[:, 6:9]
N_electrodes = dobs[:, 9:12]
dobs = dobs[:, -1]

# Define survey
unique_tx, k = np.unique(np.c_[A_electrodes, B_electrodes], axis=0, return_index=True)
n_sources = len(k)
k = np.sort(k)
k = np.r_[k, len(k) + 1]

source_list = []
for ii in range(0, n_sources):
    # MN electrode locations for receivers. Each is an (N, 3) numpy array
    M_locations = M_electrodes[k[ii] : k[ii + 1], :]
    N_locations = N_electrodes[k[ii] : k[ii + 1], :]
    receiver_list = [
        dc.receivers.Dipole(
            M_locations,
            N_locations,
            data_type="apparent_resistivity",
        )
    ]

    # AB electrode locations for source. Each is a (1, 3) numpy array
    A_location = A_electrodes[k[ii], :]
    B_location = B_electrodes[k[ii], :]
    source_list.append(dc.sources.Dipole(receiver_list, A_location, B_location))

# Define survey
survey = dc.Survey(source_list)

# Plot apparent resistivities on sounding curve as a function of Wenner separation
# parameter.
electrode_separations = 0.5 * np.sqrt(
    np.sum((survey.locations_a - survey.locations_b) ** 2, axis=1)
)

fig = plt.figure(figsize=(11, 5))
mpl.rcParams.update({"font.size": 14})
ax1 = fig.add_axes([0.15, 0.1, 0.7, 0.85])
ax1.semilogy(electrode_separations, dobs, "b")
ax1.set_xlabel("AB/2 (m)")
ax1.set_ylabel(r"Apparent Resistivity ($\Omega m$)")
plt.show()

Assign Uncertainties#

Inversion with SimPEG requires that we define standard deviation on our data. This represents our estimate of the noise in our data. For DC sounding data, a relative error is applied to each datum. For this tutorial, the relative error on each datum will be 2.5%.

std = 0.025 * dobs

Define Data#

Here is where we define the data that are inverted. The data are defined by the survey, the observation values and the standard deviation.

data_object = data.Data(survey, dobs=dobs, standard_deviation=std)

Defining the Starting Model and Mapping#

In this case, the model consists of parameters which define the respective resistivities and thickness for a set of horizontal layer. Here, we choose to define a model consisting of 3 layers.

# Define the resistivities and thicknesses for the starting model. The thickness
# of the bottom layer is assumed to extend downward to infinity so we don't
# need to define it.
resistivities = np.r_[1e3, 1e3, 1e3]
layer_thicknesses = np.r_[50.0, 50.0]

# Define a mesh for plotting and regularization.
mesh = TensorMesh([(np.r_[layer_thicknesses, layer_thicknesses[-1]])], "0")
print(mesh)

# Define model. We are inverting for the layer resistivities and layer thicknesses.
# Since the bottom layer extends to infinity, it is not a model parameter for
# which we need to invert. For a 3 layer model, there is a total of 5 parameters.
# For stability, our model is the log-resistivity and log-thickness.
starting_model = np.r_[np.log(resistivities), np.log(layer_thicknesses)]

# Since the model contains two different properties for each layer, we use
# wire maps to distinguish the properties.
wire_map = maps.Wires(("rho", mesh.nC), ("t", mesh.nC - 1))
resistivity_map = maps.ExpMap(nP=mesh.nC) * wire_map.rho
layer_map = maps.ExpMap(nP=mesh.nC - 1) * wire_map.t

TensorMesh: 3 cells

                    MESH EXTENT             CELL WIDTH      FACTOR
dir    nC        min           max         min       max      max
---   ---  ---------------------------  ------------------  ------
 x      3          0.00        150.00     50.00     50.00    1.00

Define the Physics#

Here we define the physics of the problem. The data consists of apparent resistivity values. This is defined here.

simulation = dc.simulation_1d.Simulation1DLayers(
    survey=survey,
    rhoMap=resistivity_map,
    thicknessesMap=layer_map,
)

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.
# Within the data misfit, the residual between predicted and observed data are
# normalized by the data's standard deviation.
dmis = data_misfit.L2DataMisfit(simulation=simulation, data=data_object)

# Define the regularization on the parameters related to resistivity
mesh_rho = TensorMesh([mesh.h[0].size])
reg_rho = regularization.WeightedLeastSquares(
    mesh_rho, alpha_s=0.01, alpha_x=1, mapping=wire_map.rho
)

# Define the regularization on the parameters related to layer thickness
mesh_t = TensorMesh([mesh.h[0].size - 1])
reg_t = regularization.WeightedLeastSquares(
    mesh_t, alpha_s=0.01, alpha_x=1, mapping=wire_map.t
)

# Combine to make regularization for the inversion problem
reg = reg_rho + reg_t

# Define how the optimization problem is solved. Here we will use an inexact
# Gauss-Newton approach that employs the conjugate gradient solver.
opt = optimization.InexactGaussNewton(maxIter=50, maxIterCG=30)

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

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.

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

# Set the rate of reduction in trade-off parameter (beta) each time the
# the inverse problem is solved. And set the number of Gauss-Newton iterations
# for each trade-off paramter value.
beta_schedule = directives.BetaSchedule(coolingFactor=5.0, coolingRate=3.0)

# Options for outputting recovered models and predicted data for each beta.
save_iteration = directives.SaveOutputEveryIteration(save_txt=False)

# Setting a stopping criteria for the inversion.
target_misfit = directives.TargetMisfit(chifact=0.1)

# The directives are defined in a list
directives_list = [
    starting_beta,
    beta_schedule,
    target_misfit,
]

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, directiveList=directives_list)

# Run the inversion
recovered_model = inv.run(starting_model)

Running inversion with SimPEG v0.23.1.dev10+gf697d2455
simpeg.InvProblem will set Regularization.reference_model to m0.
simpeg.InvProblem will set Regularization.reference_model to m0.

                    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
============================ Inexact Gauss Newton ============================
  #     beta     phi_d     phi_m       f      |proj(x-g)-x|  LS    Comment
-----------------------------------------------------------------------------
x0 has any nan: 0
1.44e+04  1.48e+03  0.00e+00  1.48e+03    5.94e+03      0
1.44e+04  7.15e+02  3.64e-04  7.20e+02    6.58e+02      0
1.44e+04  7.09e+02  3.66e-04  7.14e+02    1.43e+01      0
2.88e+03  7.09e+02  3.63e-04  7.10e+02    7.70e+02      0   Skip BFGS
2.88e+03  6.95e+02  3.26e-03  7.04e+02    9.59e+01      0
2.88e+03  6.95e+02  3.23e-03  7.04e+02    1.10e+02      0
5.75e+02  6.95e+02  3.29e-03  6.97e+02    7.53e+02      1
5.75e+02  6.32e+02  6.42e-02  6.69e+02    3.20e+02      0
5.75e+02  6.25e+02  7.09e-02  6.66e+02    4.60e+02      1
1.15e+02  6.23e+02  7.22e-02  6.32e+02    7.99e+02      2
1.15e+02  5.97e+02  1.40e-01  6.13e+02    9.21e+02      3
1.15e+02  5.65e+02  2.64e-01  5.95e+02    5.66e+02      2
2.30e+01  5.25e+02  3.21e-01  5.32e+02    1.01e+03      3
2.30e+01  4.39e+02  1.15e+00  4.65e+02    1.33e+03      2
2.30e+01  2.62e+02  2.41e+00  3.17e+02    6.51e+02      1
4.60e+00  1.94e+02  2.97e+00  2.08e+02    6.81e+02      2
4.60e+00  1.42e+02  4.33e+00  1.62e+02    9.21e+02      2
4.60e+00  1.05e+02  7.38e+00  1.39e+02    1.67e+03      1   Skip BFGS
9.21e-01  5.57e+01  8.18e+00  6.32e+01    3.44e+02      0
9.21e-01  4.19e+01  8.40e+00  4.96e+01    4.12e+02      2
9.21e-01  3.26e+01  9.78e+00  4.16e+01    5.73e+02      2   Skip BFGS
1.84e-01  2.88e+01  1.23e+01  3.11e+01    5.77e+02      1
1.84e-01  2.24e+01  1.23e+01  2.47e+01    4.70e+02      2
1.84e-01  1.45e+01  1.61e+01  1.75e+01    6.30e+02      1   Skip BFGS
3.68e-02  6.40e+00  1.57e+01  6.98e+00    6.28e+01      0
3.68e-02  5.11e+00  1.88e+01  5.80e+00    3.74e+02      1
3.68e-02  3.08e+00  2.22e+01  3.89e+00    1.44e+02      1
7.36e-03  2.89e+00  2.24e+01  3.06e+00    1.63e+02      1   Skip BFGS
7.36e-03  2.65e+00  2.69e+01  2.85e+00    2.76e+02      0
------------------------- STOP! -------------------------
: |fc-fOld| = 0.0000e+00 <= tolF*(1+|f0|) = 1.4803e+02
: |xc-x_last| = 9.3033e-02 <= tolX*(1+|x0|) = 1.4182e+00
: |proj(x-g)-x|    = 2.7552e+02 <= tolG          = 1.0000e-01
: |proj(x-g)-x|    = 2.7552e+02 <= 1e3*eps       = 1.0000e-02
: maxIter   =      50    <= iter          =     29
------------------------- DONE! -------------------------

Examining the Results#

# Define true model and layer thicknesses
true_model = np.r_[1e3, 4e3, 2e2]
true_layers = np.r_[100.0, 100.0]

# Plot true model and recovered model
fig = plt.figure(figsize=(5, 5))

x_min = np.min([np.min(resistivity_map * recovered_model), np.min(true_model)])
x_max = np.max([np.max(resistivity_map * recovered_model), np.max(true_model)])

ax1 = fig.add_axes([0.2, 0.15, 0.7, 0.7])
plot_1d_layer_model(true_layers, true_model, ax=ax1, plot_elevation=True, color="b")
plot_1d_layer_model(
    layer_map * recovered_model,
    resistivity_map * recovered_model,
    ax=ax1,
    plot_elevation=True,
    color="r",
)
ax1.set_xlabel(r"Resistivity ($\Omega m$)")
ax1.set_xlim(0.9 * x_min, 1.1 * x_max)
ax1.legend(["True Model", "Recovered Model"])

# Plot the true and apparent resistivities on a sounding curve
fig = plt.figure(figsize=(11, 5))
ax1 = fig.add_axes([0.2, 0.05, 0.6, 0.8])
ax1.semilogy(electrode_separations, dobs, "b")
ax1.semilogy(electrode_separations, inv_prob.dpred, "r")
ax1.set_xlabel("AB/2 (m)")
ax1.set_ylabel(r"Apparent Resistivity ($\Omega m$)")
ax1.legend(["True Sounding Curve", "Predicted Sounding Curve"])
plt.show()