2.5D IP Inversion

Author: Devin C. Cowan

Keywords: induced polarization, 2.5D inversion, sparse norm, tree mesh.

Summary: Here we invert IP data to recover the subsurface chargeability distribution on a tree mesh. We demonstrate two inversion approaches for recovering a chargeability model:

1. Weighted least-squares inversion
2. Iteratively re-weighted least-squares (IRLS) inversion

Almost all of the details of defining and carrying out weighted least-squares and iteratively weighted least squares inversion were covered in the 2.5D DC Resistivity Inversion tutorial. And functionality specific to defining the IP forward simulation was covered in the 2.5D IP Forward Simulation tutorial. As we will focus primarily on content specific to IP inversion, we strongly urge the reader to work through the aforementioned tutorials prior to working through this one.

Learning Objectives:

• Assigning appropriate uncertainties to IP data
• Designing a mesh for IP inversion
• Obtaining a background conductivity/resistivity model for the IP inversion
• Analyzing inversion results

Import Modules¶

Here, we import all of the functionality required to run the notebook for the tutorial exercise. All of the functionality specific to IP is imported from simpeg.electromagnetics.static.induced_polarization. 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 and Plot Tutorial Data¶

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 the IP 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 apparent chargeabilities in V/V.

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='apparent_chargeability'
)

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

Plot IP Data in Pseudo-Section¶

Here we use the plot_pseudosection utility function to represent the observed apparent chargeability data in pseudosection. In this case, the apparent chargeabilities are represented in V/V.

In the case of secondary voltage data, you would need to convert the values to apparent chargeabilities. This accomplished by normalizing the secondary IP voltages by the DC voltages.

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 secondary 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 secondary voltage IP data 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.

For apparent chargeability data, we generally apply a floor uncertainty to all data. This is because apparent chargeability data is obtained by normalizing the secondary voltage defining the IP signature by the DC voltage. This effectively removes the influence of electrode spacing on the amplitude of the datum. Here, we apply a floor uncertainty of 5e-3 V/V to all data.

Design a (Tree) Mesh¶

IP inversion is carried out on the same mesh that is used for DC resistivity inversion. Here, we reproduce the survey-based meshing approach that was used in the 2.5D Inversion of DC Resistivity Data tutorial. The best-practices for meshes design for DC/IP problems is provided in the 2.5D Forward Simulation of DC Resistivity Data tutorial.

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¶

Use the active_from_xyz utility function to obtain the indices of the active mesh cells from topography (e.g. cells below surface).

Project Electrodes to Discretized Topography¶

Surface IP 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 surface IP 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.

Define the Background Conductivity/Resistivity¶

In order to invert IP data, we require the background conductivity/resistivity model. In practice, a conductivity/resistivity model is recovered by inverting DC resistivity data; see the 2.5D Inversion of DC Resistivity Data tutorial. Since the conductivity/resistivity model greatly impacts the sensitivities for the IP problem, IP inversion is most successful when the recovered conductivity/resistivity model accurately characterizes the true subsurface distribution of electrical properties.

For the tutorial, we will demonstrate the importance of using an appropriate background conductivity model for IP inversion. For our weighted least-squares inversion approach, we use a halfspace as the background conductivity model. For the iteratively re-weighted least-squares inversion approach, we use the true subsurface conductivity model. In both cases, the model is defined on the active mesh cells, and we must use the simpeg.maps.InjectActiveCells mapping to map from the model space to the conductivities on the entire mesh.

Weighted Least-Squares Inversion¶

Here, a weighted least-squares inversion approach is used to apparent chargeability data to recover an intrinsic chargeability model.

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.

SimPEG uses a linearized formulation for simulating IP data; see Simulation2DCellCentered or Simulation2DNodal. In this formulation, any standard definition of the chargeability can be inverted. And the units of the recovered chargeability values is the same as the data units; e.g. intrinsic chargeability (V/V or mV/V) or integrated chargeability (ms). If you are inverting secondary voltages however, the chargeability model must represent intrinsic chargeabilities ($0 \leq \eta \leq 1$) in V/V.

In this tutorial, the model parameters are the subsurface chargeabilities ($0 \leq \eta \leq 1$) in V/V. We use the simpeg.maps.InjectActiveCells map to project the active cell chargeabilities to the entire mesh.

Starting and Reference Models for IP Inversion¶

Here, we would create starting and/or reference models for the IP 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 the 1e-6 V/V. 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.

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

Since we are using a background conductivity model, we use sigma to define the background conductivity on the entire mesh. Because the IP problem has been linearized, we store the sensitivities for the inversion using storeJ.

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. Length scales along x and y are used to balance the smallness and smoothness terms. Here, length scales of 5 are used to emphasize more smoothness in the recovered model. 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 chargeabilities 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 ProjectedGNCG class to solve the optimization problem using projected Gauss-Newton with conjugate gradietn 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. Here, the lower property is set to 0 to ensure recovered chargeability values are positive.

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 applied while the inversion is running. Here, we apply common directives for weighted least-squares inversion of IP data and describe their roles. These are:

• UpdateSensitivityWeights: Apply sensitivity weighting to counteract the natural tendency of IP inversion to place materials near the electrodes. Since the IP problem has been linearized, we do not need to update the sensitivities and we set every_iteration=False. For IP 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-4 and 1e-1.

• UpdatePreconditioner: Apply Jacobi preconditioner when solving optimization problem to reduce the number of conjugate gradient iterations. We set update_every_iteration=False because the ideal preconditioner because the forward problem and regularization are not 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 data misfit. As we can see, the recovered model reproduces the observed data quite well. And the data misfit map indicates we are not overfitting certain electrodes relative to others. We applied a uniform floor uncertainty to all data. However, if the uncertainty on each datum is different, we may choose to plot the normalized data misfit; i.e. the misfit normalized by the uncertainties.

Plot the Recovered Model¶

As we can see, the inversion recovers a smooth model. We recover a confined chargeable structure centered at (-160, 50). We also recover an extreme chargeable feature at depth at x < 0 m. Where sensitivity is low, we expect the inversion to default to the reference model values (i.e. zero) given the regularization we chose. As we will see, the second feature is an artifact of choosing an inappropriate background conductivity model

Iteratively Re-weighted Least-Squares Inversion¶

Here, we provide a step-by-step best-practices approach for IRLS inversion of apparent chargeability data on a tree mesh to recover a chargeability model. Many of the steps are the same as our previous approach. As a result, we will avoid repeating information whenever possible. It should be noted that once the IRLS inversion is complete, we can extract the least-squares model the is recovered prior to performing the IRLS step. This will be done to compare least-squares inversion results for poor and appropriate background conductivity models.

Because we are changing the background conductivity, the sensitivities used for the L2 inversion cannot be reused for the IRLS inversion. Here, we define the simulation that can be used for the IRLS inversion, where the true background conductivity is true.

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. We use the same length scales that were used in the weighted least-squares inversion. Here, we apply a 0-norm on the smallness and a 2-norm on the smoothness.

Define the Optimization¶

Define the Inverse Problem¶

Provide Inversion Directives¶

Here, we create common directives for IRLS inversion of IP 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. Note that every_iteration has been set to True for the Jacobi preconditioner. This is because the IRLS regularization is model-dependent.

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¶

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. As we will see, the L2-model extracted from the IRLS inversion is a great improvement over the L2-model we recovered via weighted least-squares inversion. And the only difference is the background conductivity model.