.. 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_parametric.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_parametric.py: Parametric 1D Inversion of Sounding Data ======================================== Here we use the module *simpeg.electromangetics.static.resistivity* to invert DC resistivity sounding data and recover the resistivities and layer thicknesses for a 1D layered Earth. In this tutorial, we focus on the following: - How to define sources and receivers from a survey file - How to define the survey - Defining a model that consists of resistivities and layer thicknesses 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 22-25 Import modules -------------- .. GENERATED FROM PYTHON SOURCE LINES 25-53 .. 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 54-62 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 62-81 .. 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 82-88 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 88-139 .. code-block:: Python # Load data dobs = np.loadtxt(str(data_filename)) 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_parametric_001.png :alt: plot inv 1 dcr sounding parametric :srcset: /content/tutorials/05-dcr/images/sphx_glr_plot_inv_1_dcr_sounding_parametric_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 140-148 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.5%. .. GENERATED FROM PYTHON SOURCE LINES 148-152 .. code-block:: Python std = 0.025 * dobs .. GENERATED FROM PYTHON SOURCE LINES 153-159 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 159-162 .. code-block:: Python data_object = data.Data(survey, dobs=dobs, standard_deviation=std) .. GENERATED FROM PYTHON SOURCE LINES 163-170 Defining the Starting Model and Mapping --------------------------------------- In this case, the model consists of parameters which define the respective resistivities and thickness for a set of horizontal layer. Here, we choose to define a model consisting of 3 layers. .. GENERATED FROM PYTHON SOURCE LINES 170-193 .. code-block:: Python # Define the resistivities and thicknesses for the starting model. The thickness # of the bottom layer is assumed to extend downward to infinity so we don't # need to define it. resistivities = np.r_[1e3, 1e3, 1e3] layer_thicknesses = np.r_[50.0, 50.0] # Define a mesh for plotting and regularization. mesh = TensorMesh([(np.r_[layer_thicknesses, layer_thicknesses[-1]])], "0") print(mesh) # Define model. We are inverting for the layer resistivities and layer thicknesses. # Since the bottom layer extends to infinity, it is not a model parameter for # which we need to invert. For a 3 layer model, there is a total of 5 parameters. # For stability, our model is the log-resistivity and log-thickness. starting_model = np.r_[np.log(resistivities), np.log(layer_thicknesses)] # Since the model contains two different properties for each layer, we use # wire maps to distinguish the properties. wire_map = maps.Wires(("rho", mesh.nC), ("t", mesh.nC - 1)) resistivity_map = maps.ExpMap(nP=mesh.nC) * wire_map.rho layer_map = maps.ExpMap(nP=mesh.nC - 1) * wire_map.t .. rst-class:: sphx-glr-script-out .. code-block:: none TensorMesh: 3 cells MESH EXTENT CELL WIDTH FACTOR dir nC min max min max max --- --- --------------------------- ------------------ ------ x 3 0.00 150.00 50.00 50.00 1.00 .. GENERATED FROM PYTHON SOURCE LINES 194-200 Define the Physics ------------------ Here we define the physics of the problem. The data consists of apparent resistivity values. This is defined here. .. GENERATED FROM PYTHON SOURCE LINES 200-207 .. code-block:: Python simulation = dc.simulation_1d.Simulation1DLayers( survey=survey, rhoMap=resistivity_map, thicknessesMap=layer_map, ) .. GENERATED FROM PYTHON SOURCE LINES 208-218 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 218-247 .. 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 on the parameters related to resistivity mesh_rho = TensorMesh([mesh.h[0].size]) reg_rho = regularization.WeightedLeastSquares( mesh_rho, alpha_s=0.01, alpha_x=1, mapping=wire_map.rho ) # Define the regularization on the parameters related to layer thickness mesh_t = TensorMesh([mesh.h[0].size - 1]) reg_t = regularization.WeightedLeastSquares( mesh_t, alpha_s=0.01, alpha_x=1, mapping=wire_map.t ) # Combine to make regularization for the inversion problem reg = reg_rho + reg_t # Define how the optimization problem is solved. Here we will use an inexact # Gauss-Newton approach that employs the conjugate gradient solver. opt = optimization.InexactGaussNewton(maxIter=50, maxIterCG=30) # Define the inverse problem inv_prob = inverse_problem.BaseInvProblem(dmis, reg, opt) .. GENERATED FROM PYTHON SOURCE LINES 248-255 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 255-278 .. code-block:: Python # Defining a starting value for the trade-off parameter (beta) between the data # misfit and the regularization. starting_beta = directives.BetaEstimate_ByEig(beta0_ratio=1e1) # Set the rate of reduction in trade-off parameter (beta) each time the # the inverse problem is solved. And set the number of Gauss-Newton iterations # for each trade-off paramter value. beta_schedule = directives.BetaSchedule(coolingFactor=5.0, coolingRate=3.0) # Options for outputting recovered models and predicted data for each beta. save_iteration = directives.SaveOutputEveryIteration(save_txt=False) # Setting a stopping criteria for the inversion. target_misfit = directives.TargetMisfit(chifact=0.1) # The directives are defined in a list directives_list = [ starting_beta, beta_schedule, target_misfit, ] .. GENERATED FROM PYTHON SOURCE LINES 279-285 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 285-292 .. code-block:: Python # Here we combine the inverse problem and the set of directives inv = inversion.BaseInversion(inv_prob, directiveList=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 simpeg.InvProblem will set Regularization.reference_model to m0. simpeg.InvProblem will set Regularization.reference_model to m0. /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 ============================ Inexact Gauss Newton ============================ # beta phi_d phi_m f |proj(x-g)-x| LS Comment ----------------------------------------------------------------------------- x0 has any nan: 0 0 1.44e+04 1.48e+03 0.00e+00 1.48e+03 5.94e+03 0 1 1.44e+04 7.15e+02 3.64e-04 7.20e+02 6.58e+02 0 2 1.44e+04 7.09e+02 3.66e-04 7.14e+02 1.43e+01 0 3 2.88e+03 7.09e+02 3.62e-04 7.10e+02 7.70e+02 0 Skip BFGS 4 2.88e+03 6.95e+02 3.25e-03 7.05e+02 9.55e+01 0 5 2.88e+03 6.95e+02 3.21e-03 7.04e+02 1.08e+02 0 6 5.77e+02 6.95e+02 3.28e-03 6.97e+02 7.53e+02 1 7 5.77e+02 6.32e+02 6.39e-02 6.69e+02 3.15e+02 0 8 5.77e+02 6.26e+02 7.06e-02 6.66e+02 4.70e+02 1 9 1.15e+02 6.24e+02 7.17e-02 6.32e+02 7.96e+02 2 10 1.15e+02 5.98e+02 1.39e-01 6.14e+02 9.27e+02 3 11 1.15e+02 5.68e+02 2.61e-01 5.98e+02 5.71e+02 2 12 2.31e+01 5.26e+02 3.16e-01 5.34e+02 1.02e+03 3 13 2.31e+01 4.44e+02 1.14e+00 4.70e+02 1.36e+03 2 14 2.31e+01 2.69e+02 2.37e+00 3.23e+02 6.66e+02 1 15 4.61e+00 1.98e+02 2.91e+00 2.12e+02 6.99e+02 2 16 4.61e+00 1.42e+02 4.29e+00 1.62e+02 8.64e+02 2 17 4.61e+00 1.22e+02 7.35e+00 1.56e+02 1.91e+03 1 Skip BFGS 18 9.23e-01 6.83e+01 6.49e+00 7.43e+01 3.83e+02 0 19 9.23e-01 4.95e+01 7.84e+00 5.67e+01 6.38e+02 2 20 9.23e-01 3.86e+01 8.51e+00 4.65e+01 3.05e+02 2 21 1.85e-01 3.02e+01 9.77e+00 3.20e+01 5.00e+02 2 22 1.85e-01 1.99e+01 1.45e+01 2.26e+01 5.11e+02 1 Skip BFGS 23 1.85e-01 6.40e+00 2.58e+01 1.12e+01 3.51e+02 0 Skip BFGS 24 3.69e-02 3.65e+00 2.22e+01 4.46e+00 4.25e+01 0 25 3.69e-02 3.48e+00 2.00e+01 4.22e+00 1.18e+02 2 ------------------------- STOP! ------------------------- 1 : |fc-fOld| = 0.0000e+00 <= tolF*(1+|f0|) = 1.4803e+02 1 : |xc-x_last| = 1.6464e-01 <= tolX*(1+|x0|) = 1.4182e+00 0 : |proj(x-g)-x| = 1.1810e+02 <= tolG = 1.0000e-01 0 : |proj(x-g)-x| = 1.1810e+02 <= 1e3*eps = 1.0000e-02 0 : maxIter = 50 <= iter = 26 ------------------------- DONE! ------------------------- .. GENERATED FROM PYTHON SOURCE LINES 293-296 Examining the Results --------------------- .. GENERATED FROM PYTHON SOURCE LINES 296-329 .. code-block:: Python # Define true model and layer thicknesses true_model = np.r_[1e3, 4e3, 2e2] true_layers = np.r_[100.0, 100.0] # Plot true model and recovered model fig = plt.figure(figsize=(5, 5)) x_min = np.min([np.min(resistivity_map * recovered_model), np.min(true_model)]) x_max = np.max([np.max(resistivity_map * recovered_model), np.max(true_model)]) ax1 = fig.add_axes([0.2, 0.15, 0.7, 0.7]) plot_1d_layer_model(true_layers, true_model, ax=ax1, plot_elevation=True, color="b") plot_1d_layer_model( layer_map * recovered_model, resistivity_map * recovered_model, ax=ax1, plot_elevation=True, color="r", ) ax1.set_xlabel(r"Resistivity ($\Omega m$)") ax1.set_xlim(0.9 * x_min, 1.1 * x_max) ax1.legend(["True Model", "Recovered Model"]) # Plot the true and apparent resistivities on a sounding curve fig = plt.figure(figsize=(11, 5)) ax1 = fig.add_axes([0.2, 0.05, 0.6, 0.8]) ax1.semilogy(electrode_separations, dobs, "b") ax1.semilogy(electrode_separations, inv_prob.dpred, "r") ax1.set_xlabel("AB/2 (m)") ax1.set_ylabel(r"Apparent Resistivity ($\Omega m$)") ax1.legend(["True Sounding Curve", "Predicted Sounding Curve"]) plt.show() .. rst-class:: sphx-glr-horizontal * .. image-sg:: /content/tutorials/05-dcr/images/sphx_glr_plot_inv_1_dcr_sounding_parametric_002.png :alt: plot inv 1 dcr sounding parametric :srcset: /content/tutorials/05-dcr/images/sphx_glr_plot_inv_1_dcr_sounding_parametric_002.png :class: sphx-glr-multi-img * .. image-sg:: /content/tutorials/05-dcr/images/sphx_glr_plot_inv_1_dcr_sounding_parametric_003.png :alt: plot inv 1 dcr sounding parametric :srcset: /content/tutorials/05-dcr/images/sphx_glr_plot_inv_1_dcr_sounding_parametric_003.png :class: sphx-glr-multi-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 32.766 seconds) **Estimated memory usage:** 288 MB .. _sphx_glr_download_content_tutorials_05-dcr_plot_inv_1_dcr_sounding_parametric.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_parametric.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_inv_1_dcr_sounding_parametric.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_inv_1_dcr_sounding_parametric.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_