.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "content/tutorials/05-dcr/plot_inv_1_dcr_sounding_irls.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_content_tutorials_05-dcr_plot_inv_1_dcr_sounding_irls.py: Sparse 1D Inversion of Sounding Data ==================================== Here we use the module *simpeg.electromangetics.static.resistivity* to invert DC resistivity sounding data and recover a 1D electrical resistivity model. In this tutorial, we focus on the following: - How to define sources and receivers from a survey file - How to define the survey - 1D inversion of DC resistivity data with iteratively re-weighted least-squares 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. .. GENERATED FROM PYTHON SOURCE LINES 21-24 Import modules -------------- .. GENERATED FROM PYTHON SOURCE LINES 24-51 .. code-block:: Python 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 .. GENERATED FROM PYTHON SOURCE LINES 52-60 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" .. GENERATED FROM PYTHON SOURCE LINES 60-79 .. code-block:: Python # 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" .. rst-class:: sphx-glr-script-out .. code-block:: none 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! .. GENERATED FROM PYTHON SOURCE LINES 80-86 Load Data, Define Survey and Plot --------------------------------- Here we load the observed data, define the DC survey geometry and plot the data values. .. GENERATED FROM PYTHON SOURCE LINES 86-138 .. code-block:: Python # Load data dobs = np.loadtxt(str(data_filename)) # Extract source and receiver electrode locations and the observed data 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() .. image-sg:: /content/tutorials/05-dcr/images/sphx_glr_plot_inv_1_dcr_sounding_irls_001.png :alt: plot inv 1 dcr sounding irls :srcset: /content/tutorials/05-dcr/images/sphx_glr_plot_inv_1_dcr_sounding_irls_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 139-146 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%. .. GENERATED FROM PYTHON SOURCE LINES 146-150 .. code-block:: Python std = 0.02 * np.abs(dobs) .. GENERATED FROM PYTHON SOURCE LINES 151-157 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. .. GENERATED FROM PYTHON SOURCE LINES 157-161 .. code-block:: Python data_object = data.Data(survey, dobs=dobs, standard_deviation=std) .. GENERATED FROM PYTHON SOURCE LINES 162-168 Defining a 1D Layered Earth (1D Tensor Mesh) -------------------------------------------- Here, we define the layer thicknesses for our 1D simulation. To do this, we use the TensorMesh class. .. GENERATED FROM PYTHON SOURCE LINES 168-177 .. code-block:: Python # Define layer thicknesses layer_thicknesses = 5 * np.logspace(0, 1, 25) # Define a mesh for plotting and regularization. mesh = TensorMesh([(np.r_[layer_thicknesses, layer_thicknesses[-1]])], "0") print(mesh) .. rst-class:: sphx-glr-script-out .. code-block:: none TensorMesh: 26 cells MESH EXTENT CELL WIDTH FACTOR dir nC min max min max max --- --- --------------------------- ------------------ ------ x 26 0.00 546.90 5.00 50.00 1.10 .. GENERATED FROM PYTHON SOURCE LINES 178-190 Define a Starting and Reference Model ------------------------------------- Here, we create starting and/or reference models for the 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 log(1000) Ohm meters. Define log-resistivity values for each layer since our model is the log-resistivity. Don't make the values 0! Otherwise the gradient for the 1st iteration is zero and the inversion will not converge. .. GENERATED FROM PYTHON SOURCE LINES 190-197 .. code-block:: Python # Define model. A resistivity (Ohm meters) or conductivity (S/m) for each layer. starting_model = np.log(2e2 * np.ones((len(layer_thicknesses) + 1))) # Define mapping from model to active cells. model_map = maps.IdentityMap(nP=len(starting_model)) * maps.ExpMap() .. GENERATED FROM PYTHON SOURCE LINES 198-203 Define the Physics ------------------ Here we define the physics of the problem using the Simulation1DLayers class. .. GENERATED FROM PYTHON SOURCE LINES 203-211 .. code-block:: Python simulation = dc.simulation_1d.Simulation1DLayers( survey=survey, rhoMap=model_map, thicknesses=layer_thicknesses, ) .. GENERATED FROM PYTHON SOURCE LINES 212-222 Define Inverse Problem ---------------------- The inverse problem is defined by 3 things: 1) Data Misfit: a measure of how well our recovered model explains the field data 2) Regularization: constraints placed on the recovered model and a priori information 3) Optimization: the numerical approach used to solve the inverse problem .. GENERATED FROM PYTHON SOURCE LINES 222-245 .. code-block:: Python # 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 (model objective function). Here, 'p' defines the # the norm of the smallness term and 'q' defines the norm of the smoothness # term. reg = regularization.Sparse(mesh, mapping=model_map) reg.reference_model = starting_model p = 0 q = 0 reg.norms = [p, q] # Define how the optimization problem is solved. Here we will use an inexact # Gauss-Newton approach that employs the conjugate gradient solver. opt = optimization.ProjectedGNCG(maxIter=100, maxIterLS=20, maxIterCG=20, tolCG=1e-3) # Define the inverse problem inv_prob = inverse_problem.BaseInvProblem(dmis, reg, opt) .. GENERATED FROM PYTHON SOURCE LINES 246-253 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. .. GENERATED FROM PYTHON SOURCE LINES 253-279 .. code-block:: Python # Apply and update sensitivity weighting as the model updates update_sensitivity_weights = directives.UpdateSensitivityWeights() # Reach target misfit for L2 solution, then use IRLS until model stops changing. IRLS = directives.UpdateIRLS(max_irls_iterations=40, f_min_change=1e-5) # Defining a starting value for the trade-off parameter (beta) between the data # misfit and the regularization. starting_beta = directives.BetaEstimate_ByEig(beta0_ratio=20) # Update the preconditionner update_Jacobi = directives.UpdatePreconditioner() # Options for outputting recovered models and predicted data for each beta. save_iteration = directives.SaveOutputEveryIteration(save_txt=False) # The directives are defined as a list. directives_list = [ update_sensitivity_weights, IRLS, starting_beta, update_Jacobi, save_iteration, ] .. GENERATED FROM PYTHON SOURCE LINES 280-286 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. .. GENERATED FROM PYTHON SOURCE LINES 286-293 .. code-block:: Python # Here we combine the inverse problem and the set of directives inv = inversion.BaseInversion(inv_prob, directives_list) # Run the inversion recovered_model = inv.run(starting_model) .. rst-class:: sphx-glr-script-out .. code-block:: none Running inversion with SimPEG v0.23.0 /home/vsts/work/1/s/simpeg/simulation.py:197: DefaultSolverWarning: Using the default solver: Pardiso. If you would like to suppress this notification, add warnings.filterwarnings('ignore', simpeg.utils.solver_utils.DefaultSolverWarning) to your script. simpeg.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv. ***Done using same Solver, and solver_opts as the Simulation1DLayers problem*** model has any nan: 0 =============================== Projected GNCG =============================== # beta phi_d phi_m f |proj(x-g)-x| LS Comment ----------------------------------------------------------------------------- x0 has any nan: 0 0 1.57e-03 4.32e+04 6.45e+02 4.32e+04 3.57e+03 0 1 7.84e-04 1.11e+04 1.17e+05 1.12e+04 6.18e+03 2 2 3.92e-04 2.60e+03 3.90e+05 2.75e+03 6.05e+03 0 3 1.96e-04 8.49e+01 3.65e+05 1.56e+02 4.20e+02 0 4 9.80e-05 3.48e+01 4.76e+05 8.14e+01 1.23e+02 0 Skip BFGS Reached starting chifact with l2-norm regularization: Start IRLS steps... irls_threshold 9.846781913625243 5 9.80e-05 2.21e+01 7.13e+04 2.91e+01 1.82e+02 0 Skip BFGS 6 1.56e-04 1.26e+01 1.41e+05 3.47e+01 2.66e+02 0 Skip BFGS 7 2.51e-04 1.23e+01 1.53e+05 5.08e+01 2.70e+02 2 Skip BFGS 8 3.67e-04 1.55e+01 1.30e+05 6.32e+01 2.58e+02 0 Skip BFGS 9 6.12e-04 1.12e+01 1.35e+05 9.41e+01 1.48e+02 0 Skip BFGS 10 9.00e-04 1.52e+01 1.27e+05 1.30e+02 4.61e+01 0 Skip BFGS 11 1.15e-03 2.05e+01 1.26e+05 1.65e+02 4.05e+01 0 Skip BFGS 12 1.46e-03 2.49e+01 1.32e+05 2.17e+02 5.34e+01 0 Skip BFGS 13 9.10e-04 3.02e+01 1.38e+05 1.56e+02 5.19e+01 0 14 5.68e-04 2.47e+01 1.77e+05 1.25e+02 2.46e+01 0 15 7.00e-04 2.18e+01 2.21e+05 1.77e+02 5.73e+01 0 Skip BFGS 16 4.48e-04 2.80e+01 2.28e+05 1.30e+02 4.21e+01 0 17 2.88e-04 2.29e+01 2.98e+05 1.09e+02 1.96e+01 0 18 3.69e-04 2.02e+01 3.76e+05 1.59e+02 5.64e+01 0 Skip BFGS 19 4.73e-04 2.64e+01 3.87e+05 2.09e+02 5.31e+01 0 20 2.91e-04 3.14e+01 4.12e+05 1.51e+02 5.48e+01 0 Skip BFGS 21 1.79e-04 2.54e+01 5.29e+05 1.20e+02 3.16e+01 0 22 2.22e-04 2.14e+01 6.66e+05 1.69e+02 5.57e+01 0 Skip BFGS 23 2.77e-04 2.60e+01 7.03e+05 2.21e+02 5.44e+01 0 24 1.73e-04 3.01e+01 7.40e+05 1.58e+02 5.97e+01 0 25 1.08e-04 2.46e+01 8.81e+05 1.19e+02 4.73e+01 0 26 1.40e-04 1.97e+01 1.05e+06 1.67e+02 8.40e+01 0 27 1.73e-04 2.17e+01 1.09e+06 2.10e+02 1.02e+02 0 28 2.14e-04 2.37e+01 1.08e+06 2.56e+02 1.34e+02 0 29 2.65e-04 2.59e+01 1.03e+06 2.98e+02 1.62e+02 0 30 1.69e-04 2.86e+01 9.17e+05 1.83e+02 1.09e+02 0 31 1.07e-04 2.48e+01 8.36e+05 1.15e+02 6.48e+01 0 Skip BFGS 32 1.31e-04 2.23e+01 7.46e+05 1.20e+02 7.38e+01 0 Skip BFGS 33 1.59e-04 2.32e+01 6.04e+05 1.19e+02 8.08e+01 0 34 1.93e-04 2.39e+01 4.83e+05 1.17e+02 8.34e+01 0 Skip BFGS 35 2.35e-04 2.42e+01 3.89e+05 1.16e+02 8.16e+01 0 Skip BFGS 36 2.86e-04 2.42e+01 3.17e+05 1.15e+02 7.86e+01 0 Skip BFGS 37 3.48e-04 2.41e+01 2.61e+05 1.15e+02 7.53e+01 0 Skip BFGS 38 4.24e-04 2.39e+01 2.16e+05 1.15e+02 7.18e+01 0 Skip BFGS 39 5.15e-04 2.37e+01 1.79e+05 1.16e+02 6.80e+01 0 Skip BFGS 40 6.27e-04 2.35e+01 1.49e+05 1.17e+02 6.55e+01 0 Skip BFGS 41 7.63e-04 2.33e+01 1.24e+05 1.18e+02 6.34e+01 0 Skip BFGS 42 9.28e-04 2.32e+01 1.03e+05 1.19e+02 6.27e+01 0 Skip BFGS 43 1.13e-03 2.31e+01 8.60e+04 1.20e+02 6.47e+01 0 Skip BFGS 44 1.37e-03 2.30e+01 7.16e+04 1.21e+02 6.53e+01 0 Skip BFGS Reach maximum number of IRLS cycles: 40 ------------------------- STOP! ------------------------- 1 : |fc-fOld| = 0.0000e+00 <= tolF*(1+|f0|) = 4.3194e+03 1 : |xc-x_last| = 4.1521e-02 <= tolX*(1+|x0|) = 2.8016e+00 0 : |proj(x-g)-x| = 6.5261e+01 <= tolG = 1.0000e-01 0 : |proj(x-g)-x| = 6.5261e+01 <= 1e3*eps = 1.0000e-02 0 : maxIter = 100 <= iter = 45 ------------------------- DONE! ------------------------- .. GENERATED FROM PYTHON SOURCE LINES 294-297 Examining the Results --------------------- .. GENERATED FROM PYTHON SOURCE LINES 297-328 .. code-block:: Python # Define true model and layer thicknesses true_model = np.r_[1e3, 4e3, 2e2] true_layers = np.r_[100.0, 100.0] # Extract Least-Squares model l2_model = inv_prob.l2model # Plot true model and recovered model fig = plt.figure(figsize=(6, 4)) x_min = np.min(np.r_[model_map * recovered_model, model_map * l2_model, true_model]) x_max = np.max(np.r_[model_map * recovered_model, model_map * l2_model, true_model]) ax1 = fig.add_axes([0.2, 0.15, 0.7, 0.7]) plot_1d_layer_model(true_layers, true_model, ax=ax1, color="k") plot_1d_layer_model(layer_thicknesses, model_map * l2_model, ax=ax1, color="b") plot_1d_layer_model(layer_thicknesses, model_map * recovered_model, ax=ax1, color="r") ax1.set_xlabel(r"Resistivity ($\Omega m$)") ax1.set_xlim(0.9 * x_min, 1.1 * x_max) ax1.legend(["True Model", "L2-Model", "Sparse Model"]) # Plot the true and apparent resistivities on a sounding curve fig = plt.figure(figsize=(11, 5)) ax1 = fig.add_axes([0.2, 0.1, 0.6, 0.8]) ax1.semilogy(electrode_separations, dobs, "k") ax1.semilogy(electrode_separations, simulation.dpred(l2_model), "b") ax1.semilogy(electrode_separations, simulation.dpred(recovered_model), "r") ax1.set_xlabel("AB/2 (m)") ax1.set_ylabel(r"Apparent Resistivity ($\Omega m$)") ax1.legend(["True Sounding Curve", "Predicted (L2-Model)", "Predicted (Sparse)"]) plt.show() .. rst-class:: sphx-glr-horizontal * .. image-sg:: /content/tutorials/05-dcr/images/sphx_glr_plot_inv_1_dcr_sounding_irls_002.png :alt: plot inv 1 dcr sounding irls :srcset: /content/tutorials/05-dcr/images/sphx_glr_plot_inv_1_dcr_sounding_irls_002.png :class: sphx-glr-multi-img * .. image-sg:: /content/tutorials/05-dcr/images/sphx_glr_plot_inv_1_dcr_sounding_irls_003.png :alt: plot inv 1 dcr sounding irls :srcset: /content/tutorials/05-dcr/images/sphx_glr_plot_inv_1_dcr_sounding_irls_003.png :class: sphx-glr-multi-img .. rst-class:: sphx-glr-script-out .. code-block:: none /home/vsts/work/1/s/simpeg/utils/plot_utils.py:354: UserWarning: Attempt to set non-positive xlim on a log-scaled axis will be ignored. /home/vsts/work/1/s/tutorials/05-dcr/plot_inv_1_dcr_sounding_irls.py:315: UserWarning: Attempt to set non-positive xlim on a log-scaled axis will be ignored. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 38.391 seconds) **Estimated memory usage:** 288 MB .. _sphx_glr_download_content_tutorials_05-dcr_plot_inv_1_dcr_sounding_irls.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_inv_1_dcr_sounding_irls.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_inv_1_dcr_sounding_irls.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_inv_1_dcr_sounding_irls.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_