3D DC Resistivity Inversion

Author: Devin C. Cowan

Keywords: DC resistivity, 3D inversion, weighted least-squares, tree mesh.

Summary: Here, DC resistivity data are inverted to recover a log-conductivity model using a weighted least-squares inversion approach. Most of the functionality required for this tutorial was presented in the 3D Forward Simulation and 2.5D Inversion tutorials. We strongly urge the reader to work through the aforementioned tutorials prior to working through this one, as we will avoid significant repitition of content.

Learning Objectives:

• Demonstrate 3D DC resistivity inversion with SimPEG.
• Discuss practical aspects of 3D DC resistivity inversion.
• 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¶

Surface topography is defined as an (N, 3) numpy.ndarray.

Load DC Resistivity Data¶

Option A: DCIP3D formatted data

The read_dcip3d_ubc utility function can be used to load DC/IP data stored in a UBC-GIF formatted file. This function outputs a SimPEG Data object.

Option B: Load XYZ formatted data

For this tutorial, the observed data are organized within an XYZ formatted file. We can use the read_dcip_xyz utility function to load data in this format. This function outputs a SimPEG Data object. The data are normalized voltages in units V/A.

The dict_headers input argument can be used to import other data columns that are not required to define the SimPEG data object. In this case, we import a column that defines the survey line associated with each datum. We will use this to plot pseudosections for individual survey lines.

Plot Observed Data in 3D Pseudosection¶

Here we use the plot_3d_pseudosection utility function to represent the observed data in pseudosection. Although our original data are normalized voltages, we use the apparent_resistivity_from_voltage utility function to convert the data to apparent resistivities.

Plot Observed Data in 2D Pseudosection¶

As demonstrated in the 3D Forward Simulation tutorial, we can parse the 3D survey into a set of 2D survey lines if we know which survey line is associated with each datum. Then we can plot individual pseudosections for each survey line. We use the convert_survey_3d_to_2d_lines to generate the list of 2D surveys. And we use the plot_2d_pseudosection utility function to plot the pseudosection for the desired survey line.

Assign Uncertainties¶

Approaches for applying reasonable uncertainties to normalized voltage and apparent resistivity data were presented in the 2.5D Inversion tutorial. Here, we apply uncertainties of 1e-7 V/A + 10 % to the normalized voltage data being inverted.

Design a (Tree) Mesh¶

Please visit the 2.5D Forward Simulation tutorial to see the best practices for mesh design. The same rules for defining appropriate meshes for 2.5D simulations and inversion apply to 3D simulations and inversion.

Tutorial mesh: Here, a minimum cell width of 25 m (or 1/4 the minimum electrode spacing) is used within our survey region. The largest electrode spacing was 1000 m, so a the padding was extended at least 3000 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 surface 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.

Mapping from the Model to the Mesh¶

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 resisitivities 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 3D DC resistivity was discussed in the 3D Forward Simulation tutorial. Here, we use the Simulation3DNodal 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. For reasonable cases, we can compute and store the Jacobian explicitly for 3D problems by setting storeJ to True. By doing this, we drastically reduce the run-time of the inversion. For larger problems being run on computers with insufficient RAM, you may need to set storeJ to False and invert the data without storing the sensitivities explicitly.

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, y and z are used to balance the smallness and smoothness terms. And 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.

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-4 and 1e-1. Warning: should not be used if storeJ set to False in the simulation class.

• 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. Warning: should not be used if storeJ set to False in the simulation class.

• 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

Analyzing Inversion Outputs¶

Plot Normalized Data Misfit in 3D Pseudosection¶

To see how well the recovered model reproduces the observed data, it is a good idea to compare the predicted and observed data. Here, we plot the normalized data misfit to determine whether the observed data are fit evenly. Correlated signatures in the normalized data misfits indicate we are overfitting/underfitting certain data. And that we must modify our uncertainties and re-run the inversion.

Plot 2D Pseudosections for Individual Survey Lines¶

For large datasets comprised of many survey lines, interpretation can be challenging if we plot every datum. Here, we plot the observed data, predicted data and normalized data misfit for a specified survey line.

Plot True and Recovered Conductivity Model¶