2.5D DC Resistivity Inversion

Keywords: DC resistivity, 2.5D inversion, sparse norm, tree mesh.

Summary: Here we invert DC resistivity data on a tree mesh to recover the subsurface distribution of electric properties. To demonstrate a range of functionality within SimPEG, we apply two approaches:

1. Weighted least-squares inversion, where we invert normalized voltage data to recover a log-conductivity model
2. Iteratively re-weighted least-squares (IRLS) inversion, where we invert apparent resistivity data to recover a log-resistivity model

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. The iteratively re-weighted least-squares approach is able to recover sparse and/or blocky structures. Because this tutorial focusses primarily on inversion-related functionality, we urge the reader to become familiar with functionality explained in the 2.5D Forward Simulation tutorial before working through this one.

Learning Objectives:

• Introduce geophysical inversion with SimPEG.
• Assigning appropriate uncertainties to normalized voltage and apparent resistivity data.
• Designing a suitable mesh for 2.5D DC resistivity inversion.
• Choosing suitable parameters for the inversion.
• Specifying directives that are applied throughout the inversion.
• Apply the sensitivity weighting commonly used when inverting magnetic data.
• Inversion with weighted least-squares and sparse-norm regularizations.
• Analyzing inversion outputs.

Import Modules¶

Here, we import all of the functionality required to run the notebook for the tutorial exercise. All of the functionality specific to DC resistivity is imported from simpeg.electromagnetics.static.resistivity. We also import some useful utility functions from simpeg.utils. 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.

Load Tutorial Data and Plot¶

For most geophysical inversion projects, a reasonable inversion result can be obtained so long as the practitioner has observed data, the survey geometry and topography. For this tutorial, the observed data and topography files are provided. Here, we download and import the observed data and topography into the SimPEG framework.

Load the Topography¶

True surface topography is defined as an (N, 3) numpy.ndarray. For the 2.5D problem geometry however, topography is an (N, 2) numpy.ndarray, where the first coordinate represent along-line position and the second coordinate represents the vertical position. In this tutorial, we assume the topography and electrode locations are defined according to the 2.5D geometry.

Load DC Resistivity Data¶

Option A: DCIP2D formatted data

For this tutorial, the observed data are organized with a UBC-GIF DCIP2D formatted data file. We can use the read_dcip2d_ubc utility function to load data in this format. This function outputs a SimPEG Data object. The data are normalized voltages in units V/A.

Option B: Survey from ABMN electrode locations

If you have CSV-formatted data containing the ABMN electrode locations, you will need to:

• load the file into a numpy.ndarray
• extract the data column and the A, B, M and N electrode locations
• transform the electrode locations to the 2.5D geometry if necessary; i.e. (N, 2) numpy.ndarray
• use the generate_survey_from_abmn_locations to generate a survey object automatically.
• Define a Data object attached to the survey and observed data.

E.g. for a file containing electrode locations already formatted for a 2.5D geometry:

data_array = np.loadtxt(data_filename, skiprows=1)

dobs = data_array[:, -1]
A = data_array[:, 0:2]
B = data_array[:, 2:4]
M = data_array[:, 4:6]
N = data_array[:, 6:8]

survey = generate_survey_from_abmn_locations(
    locations_a=A,
    locations_b=B,
    locations_m=M,
    locations_n=M,
    data_type='volt'
)

dc_data = data.Data(survey, dobs=dobs)

Plot DC Resistivity Data in Pseudo-Section¶

Here we use the plot_pseudosection utility function to plot the normalized voltage data in pseudosection. We also use the apparent_resistivity_from_voltage utility function to convert the data to apparent resistivities, which are also plotted in pseudosection.

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 normalized voltage data, we generally apply a percent uncertainty and a very small floor uncertainty to all data. Differences in electrode spacing and subsurface conductivity result in measured voltages that span many orders of magnitude. A percent uncertainty ensures all data are fit equally. Depending on the quality of the data, we may choose a percent uncertainty between 2% - 10%. The floor uncertainty ensures stability when there are zero-crossings or erroneously small voltages. Here, we apply uncertainties of 1e-7 V/A + 5 %.

For apparent resistivity data, we also apply a percent uncertainty and a very small floor uncertainty to all data. A percent uncertainty will fit conductive and resistive structures equally; unlike a floor uncertainty which will prioritize fitting more resistive structures. Depending on the quality of the data, we may choose a percent uncertainty between 2% - 10%. The floor uncertainty ensures stability when there are zero-crossings or erroneously small voltages.

Design a (Tree) Mesh¶

The same rules for defining appropriate meshes for forward simulation of DC resistivity apply to inversion. Please visit the 2.5D Forward Simulation tutorial to see the best practices for mesh design.

Tutorial mesh: Here, a minimum cell width of 4 m (or 1/10 the minimum electrode spacing) is used within our survey region. The largest electrode spacing was 400 m, so a the padding was extended at least 1200 m from the survey region. Using the refine_surface method, we refine the tree mesh where there is significant topography. And using the refine_points methods, we refine based on electrodes locations. Visit the tree mesh API to see additional refinement methods.

Define the Active Cells¶

Simulated geophysical data are dependent on the subsurface distribution of physical property values. As a result, the cells lying below the surface topography are commonly referred to as ‘active cells’. And air cells, whose physical property values are fixed, are commonly referred to as ‘inactive cells’. Here, the discretize active_from_xyz utility function is used to find the indices of the active cells using the mesh and surface topography. The output quantity is a bool array.

Project Electrodes to Discretized Topography¶

Surface DC resistivity data will not be modeled accurately if the electrodes are modeled as living above or below the surface. It is especially problematic when electrodes are modeled as living in the air. Prior to inverting DC resistivity data, we must project the electrodes from their true elevation to the surface of the discretized topography. This is done using the drape_electrodes_on_topography method.

Weighted Least-Squares Inversion¶

Here, a weighted least-squares inversion approach is used to invert normalized voltage data to recover a log-conductivity model. We invert for log-conductivity because

1. It naturally enforces positive conductivity values in the recovered model
2. Electrical conductivity is a physical property that spans many orders of magnitude

Mapping from the Model to the Mesh¶

In SimPEG, the term ‘model’ is not synonymous with the physical property values defined on the mesh. For whatever model we choose, we must define a mapping from the set of model parameters (a 1D numpy.ndarray) to the physical property values of all cells in the mesh. Mappings are created using the simpeg.maps module.

For DC resistivity, the model parameters define the subsurface electrical conductivity, and the electrical conductivities of the cells in the air are given a fixed value. Here, the model defines the log-conductivity values for all active cells. We use the simpeg.maps.ExpMap to map from the model parameters to the conductivity values for all active cells. Then we use the simpeg.maps.InjectActiveCells map to project the active cell conductivities to the entire mesh. As explained in the 2.5D Forward Simulation tutorial, air cells are given a fixed value of 1e-8 S/m to ensure the forward problem is well-conditionned. And the $*$ operator is used to combine the separate mapping objects into a single mapping.

Starting/Reference Models¶

The starting model defines a reasonable starting point for the inversion and does not necessarily represent an initial estimate of the true model. 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. For DC resistivity inversion, the starting model is frequently a halfspace estimated from the set of apparent resistivities/conductivities.

The reference model is used to include a-priori information. By default, the starting model is set as the reference model. The impact of the reference model on the inversion will be discussed in another tutorial.

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 2.5D DC resistivity was discussed in the 2.5D Forward Simulation tutorial. Here, we use the Simulation2DNodal simulation which solves for the electric potential on mesh nodes.

Since our model consists of log-conductivities, we use sigmaMap to set the mapping from the model to the mesh cells. Except for in extreme cases, we can compute and store the Jacobian explicitly for 2.5D problems. By doing so, we drastically reduce the run-time of the inversion.

Define the 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.

Define the Regularization¶

To understand the role of the regularization in the inversion, please visit this online resource. Here, we use the WeightedLeastSquares regularization class to constrain the inversion result. Here, the scaling constants that balance the smallness and smoothness terms are set directly. We use the inverse square of the smallest cell size to balance the smallness and smoothness terms. Setting alpha_s to a very small value will recover the smoothest model. And the reference model is only applied to the smallness term.

By default, the regularization acts on the model parameters; which in the case are the log-conductivities of the active cells. So we need to specify which cells are active in the regularization. And 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.

Define the Optimization Algorithm¶

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.

Define the 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.

Provide Inversion Directives¶

Directives represent operations that are carried out during the inversion. Here, we apply common directives for weighted least-squares inversion of DC resistivity data and describe their roles. These are:

• UpdateSensitivityWeights: Apply sensitivity weighting to counteract the natural tendency of DC resistivity inversion to place materials near the electrodes. Since the problem is non-linear and the sensitivities are updated with every model, we set every_iteration=True. For DC resistivity inversion, we do not want to use the entire dynamic range of the sensitivities to generate our weighting. So we generally set threshold_value to a value betwewen 1e-5 and 1e-3.

• 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 2 or 3 works well for DC resistivity 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

Plot the Data Misfit¶

This step is necessary for determining whether the recovered model accurately reproduces observed anomalies. Here, we plot the observed data, predicted data for the recovered model, and the normalized data misfit. As we can see, the recovered model reproduces the observed data quite well. And the normalized misfit map indicates we are not overfitting certain electrodes relative to others.

If desired, one could also use the apparent_resistivity_from_voltage utility function compare observed and predicted data as apparent resistivities.

Plot the Recovered Model¶

As we can see, weighted least-squares regularization leads to the recovery of smooth models.

Iteratively Re-weighted Least-Squares Inversion¶

Here, we provide a step-by-step best-practices approach for IRLS inversion of apparent resistivity data on a tree mesh to recover a log-resistivity model. Many of the steps are the same as our previous approach. As a result, we will avoid repeating information whenever possible.

Define An Apparent Resistivity Survey¶

We cannot reuse a survey that was used to simulate a different data type. So we must define a new survey object for inverting apparent resistivity data. This can be done by extracting the ABMN electrode locations from the previous survey object and using the generate_survey_from_abmn_locations utility function.

Define the Data¶

Assign Uncertainties¶

Uncertainties for apparent resistivity data were discussed earlier. Here we apply uncertainties of 1e-3 $\Omega m$ + 5 % to all data.

Mapping from the Model to the Mesh¶

Here, the model defines the log-resistivity values for all active cells. We use the simpeg.maps.ExpMap to map from the model parameters to the resistivity values for all active cells. Then we use the simpeg.maps.InjectActiveCells map to project the active cell resisitivities to the entire mesh. As explained in the 2.5D Forward Simulation tutorial, air cells are given a fixed value of 1e8 $\Omega m$ to ensure the forward problem is well-conditionned. And the $*$ operator is used to combine the separate mapping objects into a single mapping.

Starting/Reference Models¶

The starting model is defined according to the median apparent resistivity value. The reference model is equal to the starting model.

Define the Forward Simulation¶

Define the Data Misfit¶

Define the Regularization¶

Here, we use the Sparse regularization class to constrain the inversion result using an IRLS approach. Length scales along x and y are used to balance the smallness and smoothness terms. Here, length scales of 100 are used to more strongly emphasize smoothness in the recovered model. 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 0-norm to the smallness and 2-norm to first-order smoothness along the x and y directions.

Define the Optimization Algorithm¶

Define the Inverse Problem¶

Provide Inversion 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 UpdateIRLS directive.

You will notice that we don’t use the BetaSchedule and TargetMisfit directives. Here, the beta cooling schedule is set in the UpdateIRLS 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¶

Plot True, L2 and IRLS Models¶

Here, we compare the models recovered from weighted least-squares and iteratively re-weighted least-squares inversion to the true model.