.. 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_2_dcr2d_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_2_dcr2d_irls.py: 2.5D DC Resistivity Inversion with Sparse Norms =============================================== Here we invert a line of DC resistivity data to recover an electrical conductivity model. We formulate the inverse problem as a least-squares optimization problem. For this tutorial, we focus on the following: - Defining the survey - Generating a mesh based on survey geometry - Including surface topography - Defining the inverse problem (data misfit, regularization, directives) - Applying sensitivity weighting - Plotting the recovered model and data misfit .. GENERATED FROM PYTHON SOURCE LINES 20-23 Import modules -------------- .. GENERATED FROM PYTHON SOURCE LINES 23-61 .. code-block:: Python import os import numpy as np import matplotlib as mpl import matplotlib.pyplot as plt from matplotlib.colors import LogNorm import tarfile from discretize import TreeMesh from discretize.utils import mkvc, refine_tree_xyz, active_from_xyz from SimPEG.utils import model_builder from SimPEG import ( maps, data_misfit, regularization, optimization, inverse_problem, inversion, directives, utils, ) from SimPEG.electromagnetics.static import resistivity as dc from SimPEG.electromagnetics.static.utils.static_utils import ( plot_pseudosection, apparent_resistivity_from_voltage, ) from SimPEG.utils.io_utils.io_utils_electromagnetics import read_dcip2d_ubc try: from pymatsolver import Pardiso as Solver except ImportError: from SimPEG import SolverLU as Solver mpl.rcParams.update({"font.size": 16}) # sphinx_gallery_thumbnail_number = 3 .. GENERATED FROM PYTHON SOURCE LINES 62-71 Define File Names ----------------- Here we provide the file paths to assets we need to run the inversion. The path to the true model conductivity and chargeability models are 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/dcr2d.tar.gz" .. GENERATED FROM PYTHON SOURCE LINES 71-91 .. code-block:: Python # storage bucket where we have the data data_source = "https://storage.googleapis.com/simpeg/doc-assets/dcr2d.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 topo_filename = dir_path + "topo_xyz.txt" data_filename = dir_path + "dc_data.obs" .. rst-class:: sphx-glr-script-out .. code-block:: none overwriting /home/ssoler/simpeg/tutorials/05-dcr/dcr2d.tar.gz Downloading https://storage.googleapis.com/simpeg/doc-assets/dcr2d.tar.gz saved to: /home/ssoler/simpeg/tutorials/05-dcr/dcr2d.tar.gz Download completed! .. GENERATED FROM PYTHON SOURCE LINES 92-100 Load Data, Define Survey and Plot --------------------------------- Here we load the observed data, define the DC and IP survey geometry and plot the data values using pseudo-sections. **Warning**: In the following example, the observations file is assumed to be sorted by sources .. GENERATED FROM PYTHON SOURCE LINES 100-105 .. code-block:: Python # Load data topo_xyz = np.loadtxt(str(topo_filename)) dc_data = read_dcip2d_ubc(data_filename, "volt", "general") .. GENERATED FROM PYTHON SOURCE LINES 106-115 Plot Observed Data in Pseudo-Section ------------------------------------ Here, we demonstrate how to plot 2D data in pseudo-section. First, we plot the actual data (voltages) in pseudo-section as a scatter plot. This allows us to visualize the pseudo-sensitivity locations for our survey. Next, we plot the data as apparent conductivities in pseudo-section with a filled contour plot. .. GENERATED FROM PYTHON SOURCE LINES 115-151 .. code-block:: Python # Plot voltages pseudo-section fig = plt.figure(figsize=(12, 5)) ax1 = fig.add_axes([0.1, 0.15, 0.75, 0.78]) plot_pseudosection( dc_data, plot_type="scatter", ax=ax1, scale="log", cbar_label="V/A", scatter_opts={"cmap": mpl.cm.viridis}, ) ax1.set_title("Normalized Voltages") plt.show() # Get apparent conductivities from volts and survey geometry apparent_conductivities = 1 / apparent_resistivity_from_voltage( dc_data.survey, dc_data.dobs ) # Plot apparent conductivity pseudo-section fig = plt.figure(figsize=(12, 5)) ax1 = fig.add_axes([0.1, 0.15, 0.75, 0.78]) plot_pseudosection( dc_data.survey, apparent_conductivities, plot_type="contourf", ax=ax1, scale="log", cbar_label="S/m", mask_topography=True, contourf_opts={"levels": 20, "cmap": mpl.cm.viridis}, ) ax1.set_title("Apparent Conductivity") plt.show() .. rst-class:: sphx-glr-horizontal * .. image-sg:: /content/tutorials/05-dcr/images/sphx_glr_plot_inv_2_dcr2d_irls_001.png :alt: Normalized Voltages :srcset: /content/tutorials/05-dcr/images/sphx_glr_plot_inv_2_dcr2d_irls_001.png :class: sphx-glr-multi-img * .. image-sg:: /content/tutorials/05-dcr/images/sphx_glr_plot_inv_2_dcr2d_irls_002.png :alt: Apparent Conductivity :srcset: /content/tutorials/05-dcr/images/sphx_glr_plot_inv_2_dcr2d_irls_002.png :class: sphx-glr-multi-img .. GENERATED FROM PYTHON SOURCE LINES 152-160 Assign Uncertainties -------------------- Inversion with SimPEG requires that we define the uncertainties on our data. This represents our estimate of the standard deviation of the noise in our data. For DC data, the uncertainties are 10% of the absolute value. .. GENERATED FROM PYTHON SOURCE LINES 160-163 .. code-block:: Python dc_data.standard_deviation = 0.05 * np.abs(dc_data.dobs) .. GENERATED FROM PYTHON SOURCE LINES 164-169 Create Tree Mesh ------------------ Here, we create the Tree mesh that will be used invert the DC data .. GENERATED FROM PYTHON SOURCE LINES 169-217 .. code-block:: Python dh = 4 # base cell width dom_width_x = 3200.0 # domain width x dom_width_z = 2400.0 # domain width z nbcx = 2 ** int(np.round(np.log(dom_width_x / dh) / np.log(2.0))) # num. base cells x nbcz = 2 ** int(np.round(np.log(dom_width_z / dh) / np.log(2.0))) # num. base cells z # Define the base mesh hx = [(dh, nbcx)] hz = [(dh, nbcz)] mesh = TreeMesh([hx, hz], x0="CN") # Mesh refinement based on topography mesh = refine_tree_xyz( mesh, topo_xyz[:, [0, 2]], octree_levels=[0, 0, 4, 4], method="surface", finalize=False, ) # Mesh refinement near transmitters and receivers. First we need to obtain the # set of unique electrode locations. electrode_locations = np.c_[ dc_data.survey.locations_a, dc_data.survey.locations_b, dc_data.survey.locations_m, dc_data.survey.locations_n, ] unique_locations = np.unique( np.reshape(electrode_locations, (4 * dc_data.survey.nD, 2)), axis=0 ) mesh = refine_tree_xyz( mesh, unique_locations, octree_levels=[4, 4], method="radial", finalize=False ) # Refine core mesh region xp, zp = np.meshgrid([-600.0, 600.0], [-400.0, 0.0]) xyz = np.c_[mkvc(xp), mkvc(zp)] mesh = refine_tree_xyz( mesh, xyz, octree_levels=[0, 0, 2, 8], method="box", finalize=False ) mesh.finalize() .. rst-class:: sphx-glr-script-out .. code-block:: none /home/ssoler/simpeg/tutorials/05-dcr/plot_inv_2_dcr2d_irls.py:182: DeprecationWarning: The surface option is deprecated as of `0.9.0` please update your code to use the `TreeMesh.refine_surface` functionality. It will be removed in a future version of discretize. /home/ssoler/mambaforge/envs/simpeg-test/lib/python3.8/site-packages/scipy/interpolate/_interpolate.py:698: RuntimeWarning: invalid value encountered in divide /home/ssoler/simpeg/tutorials/05-dcr/plot_inv_2_dcr2d_irls.py:203: DeprecationWarning: The radial option is deprecated as of `0.9.0` please update your code to use the `TreeMesh.refine_points` functionality. It will be removed in a future version of discretize. /home/ssoler/simpeg/tutorials/05-dcr/plot_inv_2_dcr2d_irls.py:210: DeprecationWarning: The box option is deprecated as of `0.9.0` please update your code to use the `TreeMesh.refine_bounding_box` functionality. It will be removed in a future version of discretize. .. GENERATED FROM PYTHON SOURCE LINES 218-226 Project Surveys to Discretized Topography ----------------------------------------- It is important that electrodes are not model as being in the air. Even if the electrodes are properly located along surface topography, they may lie above the discretized topography. This step is carried out to ensure all electrodes like on the discretized surface. .. GENERATED FROM PYTHON SOURCE LINES 226-246 .. code-block:: Python # Create 2D topography. Since our 3D topography only changes in the x direction, # it is easy to define the 2D topography projected along the survey line. For # arbitrary topography and for an arbitrary survey orientation, the user must # define the 2D topography along the survey line. topo_2d = np.unique(topo_xyz[:, [0, 2]], axis=0) # Find cells that lie below surface topography ind_active = active_from_xyz(mesh, topo_2d) # Extract survey from data object survey = dc_data.survey # Shift electrodes to the surface of discretized topography survey.drape_electrodes_on_topography(mesh, ind_active, option="top") # Reset survey in data object dc_data.survey = survey .. GENERATED FROM PYTHON SOURCE LINES 247-255 Starting/Reference Model and Mapping on Tree Mesh --------------------------------------------------- Here, we would create starting and/or reference models for the DC 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 natural log of 0.01 S/m. .. GENERATED FROM PYTHON SOURCE LINES 255-268 .. code-block:: Python # Define conductivity model in S/m (or resistivity model in Ohm m) air_conductivity = np.log(1e-8) background_conductivity = np.log(1e-2) active_map = maps.InjectActiveCells(mesh, ind_active, np.exp(air_conductivity)) nC = int(ind_active.sum()) conductivity_map = active_map * maps.ExpMap() # Define model starting_conductivity_model = background_conductivity * np.ones(nC) .. GENERATED FROM PYTHON SOURCE LINES 269-274 Define the Physics of the DC Simulation --------------------------------------- Here, we define the physics of the DC resistivity problem. .. GENERATED FROM PYTHON SOURCE LINES 274-280 .. code-block:: Python # Define the problem. Define the cells below topography and the mapping simulation = dc.simulation_2d.Simulation2DNodal( mesh, survey=survey, sigmaMap=conductivity_map, solver=Solver, storeJ=True ) .. GENERATED FROM PYTHON SOURCE LINES 281-291 Define DC 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 291-328 .. 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(data=dc_data, simulation=simulation) # Define the regularization (model objective function). Here, 'p' defines the # the norm of the smallness term, 'qx' defines the norm of the smoothness # in x and 'qz' defines the norm of the smoothness in z. regmap = maps.IdentityMap(nP=int(ind_active.sum())) reg = regularization.Sparse( mesh, active_cells=ind_active, reference_model=starting_conductivity_model, mapping=regmap, gradient_type="total", alpha_s=0.01, alpha_x=1, alpha_y=1, ) reg.reference_model_in_smooth = True # Include reference model in smoothness p = 0 qx = 1 qz = 1 reg.norms = [p, qx, qz] # Define how the optimization problem is solved. Here we will use an inexact # Gauss-Newton approach. opt = optimization.InexactGaussNewton(maxIter=40) # Here we define the inverse problem that is to be solved inv_prob = inverse_problem.BaseInvProblem(dmis, reg, opt) .. GENERATED FROM PYTHON SOURCE LINES 329-336 Define DC 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 336-363 .. code-block:: Python # Apply and update sensitivity weighting as the model updates update_sensitivity_weighting = directives.UpdateSensitivityWeights() # Reach target misfit for L2 solution, then use IRLS until model stops changing. update_IRLS = directives.Update_IRLS( max_irls_iterations=25, minGNiter=1, chifact_start=1.0 ) # 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) # Options for outputting recovered models and predicted data for each beta. save_iteration = directives.SaveOutputEveryIteration(save_txt=False) # Update preconditioner update_jacobi = directives.UpdatePreconditioner() directives_list = [ update_sensitivity_weighting, update_IRLS, starting_beta, save_iteration, update_jacobi, ] .. GENERATED FROM PYTHON SOURCE LINES 364-370 Running the DC 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 370-377 .. code-block:: Python # Here we combine the inverse problem and the set of directives dc_inversion = inversion.BaseInversion(inv_prob, directiveList=directives_list) # Run inversion recovered_conductivity_model = dc_inversion.run(starting_conductivity_model) .. rst-class:: sphx-glr-script-out .. code-block:: none SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv. ***Done using same Solver, and solver_opts as the Simulation2DNodal 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.90e+03 3.43e+04 0.00e+00 3.43e+04 6.52e+03 0 1 9.50e+02 5.06e+03 2.89e+00 7.81e+03 8.67e+02 0 2 4.75e+02 2.28e+03 4.74e+00 4.53e+03 3.95e+02 0 Skip BFGS 3 2.38e+02 1.07e+03 6.54e+00 2.62e+03 2.40e+02 0 Skip BFGS 4 1.19e+02 4.74e+02 8.28e+00 1.46e+03 1.40e+02 0 Skip BFGS 5 5.94e+01 2.13e+02 9.80e+00 7.95e+02 7.89e+01 0 Skip BFGS Reached starting chifact with l2-norm regularization: Start IRLS steps... irls_threshold 2.4994906088963162 6 2.97e+01 1.01e+02 1.68e+01 6.01e+02 3.08e+01 0 Skip BFGS 7 5.92e+01 6.79e+01 2.13e+01 1.33e+03 1.15e+02 0 8 4.42e+01 1.82e+02 2.34e+01 1.22e+03 4.99e+01 0 9 3.15e+01 1.99e+02 2.83e+01 1.09e+03 2.10e+02 0 Skip BFGS 10 2.55e+01 1.56e+02 2.83e+01 8.77e+02 3.94e+01 0 11 2.55e+01 1.22e+02 2.85e+01 8.49e+02 2.35e+01 0 Skip BFGS 12 3.98e+01 1.21e+02 2.76e+01 1.22e+03 9.52e+01 0 Skip BFGS 13 2.79e+01 2.06e+02 2.51e+01 9.08e+02 4.63e+01 0 14 2.79e+01 1.39e+02 2.48e+01 8.32e+02 4.11e+01 0 15 2.79e+01 1.29e+02 2.34e+01 7.82e+02 5.45e+01 0 Skip BFGS 16 2.79e+01 1.26e+02 2.18e+01 7.36e+02 3.43e+01 0 Skip BFGS 17 2.79e+01 1.23e+02 2.05e+01 6.95e+02 3.15e+01 0 Skip BFGS 18 2.79e+01 1.22e+02 1.92e+01 6.60e+02 3.23e+01 0 Skip BFGS 19 2.79e+01 1.22e+02 1.82e+01 6.29e+02 3.28e+01 0 Skip BFGS 20 2.79e+01 1.23e+02 1.75e+01 6.13e+02 3.35e+01 0 Skip BFGS 21 2.79e+01 1.27e+02 1.70e+01 6.02e+02 3.37e+01 0 Skip BFGS 22 2.79e+01 1.31e+02 1.64e+01 5.91e+02 3.41e+01 0 Skip BFGS 23 2.79e+01 1.36e+02 1.59e+01 5.80e+02 3.45e+01 0 Skip BFGS 24 2.79e+01 1.42e+02 1.53e+01 5.70e+02 3.49e+01 0 Skip BFGS 25 2.79e+01 1.47e+02 1.48e+01 5.61e+02 3.54e+01 0 Skip BFGS 26 2.29e+01 1.52e+02 1.43e+01 4.80e+02 2.19e+01 0 Skip BFGS 27 2.29e+01 1.39e+02 1.45e+01 4.70e+02 2.89e+01 0 28 2.29e+01 1.40e+02 1.42e+01 4.64e+02 3.24e+01 0 29 2.29e+01 1.41e+02 1.39e+01 4.58e+02 3.44e+01 0 30 2.29e+01 1.41e+02 1.36e+01 4.53e+02 3.62e+01 0 Skip BFGS Reach maximum number of IRLS cycles: 25 ------------------------- STOP! ------------------------- 1 : |fc-fOld| = 0.0000e+00 <= tolF*(1+|f0|) = 3.4278e+03 1 : |xc-x_last| = 3.8056e-01 <= tolX*(1+|x0|) = 3.1119e+01 0 : |proj(x-g)-x| = 3.6201e+01 <= tolG = 1.0000e-01 0 : |proj(x-g)-x| = 3.6201e+01 <= 1e3*eps = 1.0000e-02 0 : maxIter = 40 <= iter = 31 ------------------------- DONE! ------------------------- .. GENERATED FROM PYTHON SOURCE LINES 378-381 Recreate True Conductivity Model -------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 381-398 .. code-block:: Python true_background_conductivity = 1e-2 true_conductor_conductivity = 1e-1 true_resistor_conductivity = 1e-3 true_conductivity_model = true_background_conductivity * np.ones(len(mesh)) ind_conductor = model_builder.get_indices_sphere( np.r_[-120.0, -180.0], 60.0, mesh.gridCC ) true_conductivity_model[ind_conductor] = true_conductor_conductivity ind_resistor = model_builder.get_indices_sphere(np.r_[120.0, -180.0], 60.0, mesh.gridCC) true_conductivity_model[ind_resistor] = true_resistor_conductivity true_conductivity_model[~ind_active] = np.NaN .. GENERATED FROM PYTHON SOURCE LINES 399-402 Plotting True and Recovered Conductivity Model ---------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 402-449 .. code-block:: Python # Get L2 and sparse recovered model in base 10 l2_conductivity = conductivity_map * inv_prob.l2model l2_conductivity[~ind_active] = np.NaN recovered_conductivity = conductivity_map * recovered_conductivity_model recovered_conductivity[~ind_active] = np.NaN # Plot True Model norm = LogNorm(vmin=1e-3, vmax=1e-1) fig = plt.figure(figsize=(9, 15)) ax1 = 3 * [None] ax2 = 3 * [None] title_str = [ "True Conductivity Model", "Smooth Recovered Model", "Sparse Recovered Model", ] plotting_model = [ true_conductivity_model, l2_conductivity, recovered_conductivity, ] for ii in range(0, 3): ax1[ii] = fig.add_axes([0.14, 0.75 - 0.3 * ii, 0.68, 0.2]) mesh.plot_image( plotting_model[ii], ax=ax1[ii], grid=False, range_x=[-700, 700], range_y=[-600, 0], pcolor_opts={"norm": norm}, ) ax1[ii].set_xlim(-600, 600) ax1[ii].set_ylim(-600, 0) ax1[ii].set_title(title_str[ii]) ax1[ii].set_xlabel("x (m)") ax1[ii].set_ylabel("z (m)") ax2[ii] = fig.add_axes([0.84, 0.75 - 0.3 * ii, 0.03, 0.2]) cbar = mpl.colorbar.ColorbarBase(ax2[ii], norm=norm, orientation="vertical") cbar.set_label(r"$\sigma$ (S/m)", rotation=270, labelpad=15, size=12) plt.show() .. image-sg:: /content/tutorials/05-dcr/images/sphx_glr_plot_inv_2_dcr2d_irls_003.png :alt: True Conductivity Model, Smooth Recovered Model, Sparse Recovered Model :srcset: /content/tutorials/05-dcr/images/sphx_glr_plot_inv_2_dcr2d_irls_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 450-453 Plotting Predicted DC Data and Misfit ------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 453-488 .. code-block:: Python # Predicted data from recovered model dpred = inv_prob.dpred dobs = dc_data.dobs std = dc_data.standard_deviation # Plot fig = plt.figure(figsize=(9, 13)) data_array = [np.abs(dobs), np.abs(dpred), (dobs - dpred) / std] plot_title = ["Observed Voltage", "Predicted Voltage", "Normalized Misfit"] plot_units = ["V/A", "V/A", ""] scale = ["log", "log", "linear"] ax1 = 3 * [None] cax1 = 3 * [None] cbar = 3 * [None] cplot = 3 * [None] for ii in range(0, 3): ax1[ii] = fig.add_axes([0.15, 0.72 - 0.33 * ii, 0.65, 0.21]) cax1[ii] = fig.add_axes([0.81, 0.72 - 0.33 * ii, 0.03, 0.21]) cplot[ii] = plot_pseudosection( survey, data_array[ii], "contourf", ax=ax1[ii], cax=cax1[ii], scale=scale[ii], cbar_label=plot_units[ii], mask_topography=True, contourf_opts={"levels": 25, "cmap": mpl.cm.viridis}, ) ax1[ii].set_title(plot_title[ii]) plt.show() .. image-sg:: /content/tutorials/05-dcr/images/sphx_glr_plot_inv_2_dcr2d_irls_004.png :alt: Observed Voltage, Predicted Voltage, Normalized Misfit :srcset: /content/tutorials/05-dcr/images/sphx_glr_plot_inv_2_dcr2d_irls_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (2 minutes 46.683 seconds) **Estimated memory usage:** 9 MB .. _sphx_glr_download_content_tutorials_05-dcr_plot_inv_2_dcr2d_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_2_dcr2d_irls.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_inv_2_dcr2d_irls.py ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_