1D Inversion for a Single Sounding

Author: Devin C. Cowan

Keywords: total magnetic intensity, integral formulation, inversion, sparse norm, tensor mesh, tree mesh.

Summary: Here we invert total magnetic intensity data to recover a susceptibility model. We demonstrate two approaches for recovering a susceptibility model:

1. Weighted least-squares inversion for a tensor mesh
2. Iteratively re-weighted least-squares (IRLS) inversion for a tree mesh

The weighted least-squares approach is a great introduction to geophysical inversion with SimPEG. One drawback however, is that it recovers smooth structures which may not be representative of the true model. To recover sparse and/or blocky structures, we demonstrate the iteratively re-weighted least-squares approach.

Learning Objectives: Because this tutorial focusses primarily on inversion-related functionality, we urge the reader to become familiar with functionality explained in the 1D Forward Simulation of Time-Domain EM Data for a Single Sounding tutorial before working through this one. For this tutorial, we focus on:

• How to carry out 1D geophysical inversion with SimPEG.
• How to assign appropriate uncertainties to TDEM data.
• Choosing suitable parameters for the inversion.
• Specifying directives that are applied throughout the inversion.
• Weighted least-squares, sparse-norm and parametric inversion.
• Analyzing inversion outputs.

Importing Modules¶

Here, we import all of the functionality required to run the notebook for the tutorial exercise. All of the functionality specific to the forward simulation of 1D time domain EM data are imported from the simpeg.electromagnetics.time_domain module. Classes required to define the data misfit, regularization, optimization, etc... are imported from elsewhere within SimPEG. We also import some useful utility functions from simpeg.utils. To generate the mesh used for the inversion, we use the discretize package.

Download and Extract the Tutorial Data¶

For this tutorial, the frequencies and observed data for 1D sounding are stored within a tar-file. Here, we download and extract the data file.

Load and Plot the Data¶

Here we load and plot the 1D sounding data. In this case, we have B-field data for a step-off excitation. The columns of the data file are: times (s) and Bz (T).

Defining the Survey¶

Here, we define the survey geometry. For a comprehensive description of constructing TDEM surveys in SimPEG, see the 1D Forward Simulation of Time-Domain EM Data for a Single Sounding tutorial.

Here, the survey consisted of a circular transmitter loop with a radius of 6 m located 20 m above the Earth’s surface. The receiver measured the vertical component of the secondary magnetic field at a the loop’s centre.

Assign Uncertainties¶

Inversion with SimPEG requires that we define the uncertainties on our data; that is, an estimate of the standard deviation of the noise on our data assuming it is uncorrelated Gaussian with zero mean. An online resource explaining uncertainties and their role in the inversion can be found here.

For off-time TDEM data, a percent uncertainty plus a small floor value is generally applied to all data. Depending on many factors, we may apply a percent uncertainty between 5% and 20%. The floor uncertainty ensure we do not try to overfit late time channels when the signal to noise become sufficiently small. For systems where multiple field directions are measured for the same source, we may not want to apply a uniform percent uncertainty to all data. Doing so may cause the inversion to overfit weaker components. In this case, the uncertainty may be a percent of the amplitude of the secondary field.

Defining the Data¶

The SimPEG Data class is required for inversion and connects the observed data, uncertainties and survey geometry.

Weighted Least-Squares Inversion¶

Here, we use the weighted least-squares inversion approach to recover the log-conductivities on a 1D layered Earth. We impose no a-priori information about the number of layers (geological units) or their thicknesses. Instead, we define a large number of layers with exponentially increasing thicknesses. And the depth, thickness and electrical properties of the Earth are inferred from the recovered model.

Defining a 1D Layered Earth¶

Let us assume we have a reasonable estimate of the regional conductivity within our area of interest. For the earliest time channel and the estimated conductivity, we compute the minimum diffusion distance:

The minimum layer thickness is some fraction of the minimum diffusion distance. Next, we use the latest time channel and the estimated conductivity to compute the maximum diffusion distance:

Starting from our minimum layer thickness, we continue to add layers with exponentially increasing thicknesses. We do so until the layers extend to some multiple of the maximum diffusion distance.

Model and Mapping to Layer Conductivities¶

Recall from the 1D Forward Simulation of Time-Domain EM Data for a Single Sounding tutorial that the ‘model’ is not necessarily synonymous with physical property values. And that we need to define a mapping from the model to the set of input parameters required by the forward simulation. When inverting to recover electrical conductivities (or resistivities), it is best to use the log-value, as electrical conductivities of rocks span many order of magnitude.

Here, the model defines the log-conductivity values for a defined set of subsurface layers. And we use the simpeg.maps.ExpMap to map from the model parameters to the conductivity values required by the forward simulation.

Starting/Reference Models¶

The starting model defines a reasonable starting point for the inversion. Because electromagnetic problems are non-linear, your choice in starting model does have an impact on the recovered model. For DC resistivity inversion, we generally choose our starting model based on apparent resistivities. For the tutorial example, the apparent resistivities were near 1000 $\Omega m$. It should be noted that the starting model cannot be vector of zeros, otherwise the inversion will be unable to compute a gradient direction at the first iteration.

The reference model is used to include a-priori information. The impact of the reference model on the inversion will be discussed in another tutorial. The reference model for basic inversion approaches is either zero or equal to the starting model.

Notice that the length of the starting and reference models is equal to the number of model parameters!!!

Define the Forward Simulation¶

A simulation object defining the forward problem is required in order to predict data and calculate misfits for recovered models. A comprehensive description of the simulation object for 1D DC resistivity was discussed in the 1D Forward Simulation of Time-Domain EM Data for a Single Sounding tutorial. Here, we use the Simulation1DLayers which simulates the data according to a 1D Hankel transform solution.

The layer thicknesses are a static property of the simulation, and we set them using the thicknessess keyword argument. Since our model consists of log-conductivities, we use sigmaMap to set the mapping from the model to the layer conductivities.

Data Misfit¶

To understand the role of the data misfit in the inversion, please visit this online resource. Here, we use the L2DataMisfit class to define the data misfit. In this case, the data misfit is the L2 norm of the weighted residual between the observed data and the data predicted for a given model. When instantiating the data misfit object within SimPEG, we must assign an appropriate data object and simulation object as properties.

Regularization¶

To understand the role of the regularization in the inversion, please visit this online resource.

To define the regularization within SimPEG, we must define a 1D tensor mesh. Meshes are designed using the discretize package. Whereas layer thicknesses and our model are defined from our top-layer down, tensor meshes are defined from the bottom up. So to define a 1D tensor mesh for the regularization, we:

• add an extra layer to the end of our thicknesses so that the number of cells in the 1D mesh equals the number of model parameters
• reverse the order so that the model parameters in the regularization match up with the appropriate cell
• define the tensor mesh from the cell widths

By default, the regularization acts on the model parameters. In this case, the model parameters are the log-resistivities, not the electric resistivities!!! Here, we use the WeightedLeastSquares regularization class to constrain the inversion result. Here, length scale along x are used to balance the smallness and smoothness terms; yes, x is smoothness along the vertical direction. And the reference model is only applied to the smallness term. If we wanted to apply the regularization to a function of the model parameters, we would need to set an approprate mapping object using the mapping keyword argument.

Optimization¶

Here, we use the InexactGaussNewton class to solve the optimization problem using the inexact Gauss Newton with conjugate gradient solver. Reasonable default values have generally been set for the properties of each optimization class. However, the user may choose to set custom values; e.g. the accuracy tolerance for the conjugate gradient solver or the number of line searches.

Inverse Problem¶

We use the BaseInvProblem class to fully define the inverse problem that is solved at each beta (trade-off parameter) iteration. The inverse problem requires appropriate data misfit, regularization and optimization objects.

Inversion Directives¶

To understand the role of directives in the inversion, please visit this online resource. Here, we apply common directives for weighted least-squares inversion of gravity data and describe their roles. These are:

• UpdatePreconditioner: Apply Jacobi preconditioner when solving optimization problem to reduce the number of conjugate gradient iterations. We set update_every_iteration=True because the ideal preconditioner is model-dependent.

• BetaEstimate_ByEig: Compute and set starting trade-off parameter (beta) based on largest eigenvalues.

• BetaSchedule: Size reduction of the trade-off parameter at every beta iteration, and the number of Gauss-Newton iterations for each beta. In general, a coolingFactor between 1.5 and 2.5, and coolingRate of 3 works well for TDEM inversion. Cooling beta too quickly will result in portions of the model getting trapped in local minima. And we will not be finding the solution that minimizes the optimization problem if the cooling rate is too small.

• TargetMisfit: Terminates the inversion when the data misfit equals the target misfit. A chifact=1 terminates the inversion when the data misfit equals the number of data.

The directive objects are organized in a list. Upon starting the inversion or updating the recovered model at each iteration, the inversion will call each directive within the list in order. The order of the directives matters, and SimPEG will throw an error if directives are organized into an improper order. Some directives, like the BetaEstimate_ByEig are only used when starting the inversion. Other directives, like UpdatePreconditionner, are used whenever the model is updated.

Define and Run the Inversion¶

We define the inversion using the BaseInversion class. The inversion class must be instantiated with an appropriate inverse problem object and directives list. The run method, along with a starting model, is respondible for running the inversion. The output is a 1D numpy.ndarray containing the recovered model parameters

Inversion Outputs¶

Data Misfit¶

Recovered Model¶

Iteratively Re-weighted Least-Squares Inversion¶

Here, we use the iteratively reweighted least-squares (IRLS) inversion approach to recover sparse and/or blocky models on the set layers.

Define the Forward Simulation¶

Define the Data Misfit¶

Regularization¶

Here, we use the Sparse regularization class to constrain the inversion result using an IRLS approach. Here, the scaling constants that balance the smallness and smoothness terms are set directly. Equal emphasis on smallness and smoothness is generally applied by using the inverse square of the smallest cell dimension. The reference model is only applied to the smallness term; which is redundant for the tutorial example since we have set the reference model to an array of zeros. Here, we apply a 1-norm to the smallness term and a 1-norm to first-order smoothness along the x (vertical direction).

Optimization¶

Inverse Problem¶

Directives¶

Here, we create common directives for IRLS inversion of total magnetic intensity data and describe their roles. In additon to the UpdateSensitivityWeights, UpdatePreconditioner and BetaEstimate_ByEig (described before), inversion with sparse-norms requires the Update_IRLS directive.

You will notice that we don’t use the BetaSchedule and TargetMisfit directives. Here, the beta cooling schedule is set in the Update_IRLS directive using the coolingFactor and coolingRate properties. The target misfit for the L2 portion of the IRLS approach is set with the chifact_start property.

Define and Run the Inversion¶

Data Misfit and Recovered Model¶

Parametric Inversion¶

Here, we assume the subsurface is defined by a 3-layered Earth. However, the electrical properties and thicknesses of the layers are unknown. Here, we define our model to include log-conductivities and log-thicknesses. When including quantities that span different scales, it is frequently best to define the model in terms of log values so that each quantity influences the predicted data evenly.

Model and Mapping¶

For a 3-layered Earth model, the model consists of 2 log-thicknesses and 3 log-conductivities. Similar to the 1D Forward Simulation of Time Domain EM Data for a Single Sounding tutorial, need a mapping that extract log-thicknesses and log-resistivities from the model, and mappings that convert log-values to property values. For this, we require the simpeg.maps.Wires mapping and simpeg.maps.ExpMap mapping classes. Note that successive mappings can be chained together using the $*$ operator.

Starting and Reference Model¶

This problem is highly non-linear so it is important to have a reasonable estimate of the true model.

Forward Simulation¶

Because the layer thicknesses are part of the model, we define the thicknessesMap. Because we are working in terms of electrical resistivity, we must define the rhoMap.

Data Misfit¶

(Combo) Regularization¶

We need to define a regularization for each model parameter type. In this case, we have log-thicknesses and log-conductivities. For each model parameter type, we create a 1D tensor mesh with length equal to the number of parameters. In the mapping keyword argument, we used the wire map that extracts the specific model parameters from the model.

Using the $*$ operator, separate regularizations can be summed to form a regularization that is also a ComboObjectiveFunction. By setting the multipliers property, we can emphasize the relative contributions of the log-thicknesses and log-conductivities regularizations.

Optimization¶

Inverse Problem¶

Directives¶

Define and Run Inversion¶

Data Misfit and Recovered Model¶