.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "content/examples/04-dcip/plot_inv_dcip_dipoledipole_parametric_inversion.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_examples_04-dcip_plot_inv_dcip_dipoledipole_parametric_inversion.py: Parametric DC inversion with Dipole Dipole array ================================================ This is an example for a parametric inversion with a DC survey. Resistivity structure of the subsurface is parameterized as following parameters: - sigma_background: background conductivity - sigma_block: block conductivity - block_x0: horizotontal location of the block (center) - block_dx: width of the block - block_y0: depth of the block (center) - block_dy: thickness of the block User is promoted to try different initial values of the parameterized model. .. GENERATED FROM PYTHON SOURCE LINES 18-230 .. rst-class:: sphx-glr-horizontal * .. image-sg:: /content/examples/04-dcip/images/sphx_glr_plot_inv_dcip_dipoledipole_parametric_inversion_001.png :alt: True resistivity model :srcset: /content/examples/04-dcip/images/sphx_glr_plot_inv_dcip_dipoledipole_parametric_inversion_001.png :class: sphx-glr-multi-img * .. image-sg:: /content/examples/04-dcip/images/sphx_glr_plot_inv_dcip_dipoledipole_parametric_inversion_002.png :alt: Initial resistivity model :srcset: /content/examples/04-dcip/images/sphx_glr_plot_inv_dcip_dipoledipole_parametric_inversion_002.png :class: sphx-glr-multi-img * .. image-sg:: /content/examples/04-dcip/images/sphx_glr_plot_inv_dcip_dipoledipole_parametric_inversion_003.png :alt: plot inv dcip dipoledipole parametric inversion :srcset: /content/examples/04-dcip/images/sphx_glr_plot_inv_dcip_dipoledipole_parametric_inversion_003.png :class: sphx-glr-multi-img * .. image-sg:: /content/examples/04-dcip/images/sphx_glr_plot_inv_dcip_dipoledipole_parametric_inversion_004.png :alt: plot inv dcip dipoledipole parametric inversion :srcset: /content/examples/04-dcip/images/sphx_glr_plot_inv_dcip_dipoledipole_parametric_inversion_004.png :class: sphx-glr-multi-img * .. image-sg:: /content/examples/04-dcip/images/sphx_glr_plot_inv_dcip_dipoledipole_parametric_inversion_005.png :alt: True resistivity model, Recovered resistivity model :srcset: /content/examples/04-dcip/images/sphx_glr_plot_inv_dcip_dipoledipole_parametric_inversion_005.png :class: sphx-glr-multi-img .. rst-class:: sphx-glr-script-out .. code-block:: none /home/vsts/work/1/s/examples/04-dcip/plot_inv_dcip_dipoledipole_parametric_inversion.py:55: UserWarning: code under construction - API might change in the future dx is set to 2.5 m (samllest electrode spacing (10.0) / 4) dz (1.25 m) is set to dx (2.5 m) / 2 SimPEG.InvProblem will set Regularization.reference_model to m0. SimPEG.InvProblem will set Regularization.reference_model to m0. SimPEG.InvProblem will set Regularization.reference_model to m0. 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 0.00e+00 7.27e+03 0.00e+00 7.27e+03 3.06e+04 0 1 0.00e+00 3.95e+03 7.73e-03 3.95e+03 2.76e+03 0 2 0.00e+00 2.02e+03 7.95e-01 2.02e+03 9.49e+03 0 Skip BFGS 3 0.00e+00 1.36e+03 7.83e-01 1.36e+03 6.05e+02 0 4 0.00e+00 1.01e+03 1.64e+00 1.01e+03 5.19e+03 8 5 0.00e+00 9.07e+02 2.78e+00 9.07e+02 6.02e+03 3 6 0.00e+00 8.79e+02 3.18e+00 8.79e+02 5.78e+03 3 ------------------------------------------------------------------ 0 : ft = 9.1111e+02 <= alp*descent = 8.7897e+02 1 : maxIterLS = 10 <= iterLS = 10 ------------------------- End Linesearch ------------------------- The linesearch got broken. Boo. | .. code-block:: Python from SimPEG.electromagnetics.static import resistivity as DC, utils as DCutils from discretize import TensorMesh from discretize.utils import active_from_xyz from SimPEG import ( maps, utils, data_misfit, regularization, optimization, inversion, inverse_problem, directives, ) import matplotlib.pyplot as plt from matplotlib import colors import numpy as np from pylab import hist try: from pymatsolver import PardisoSolver as Solver except ImportError: from SimPEG import SolverLU as Solver def run( plotIt=True, survey_type="dipole-dipole", rho_background=1e3, rho_block=1e2, block_x0=100, block_dx=10, block_y0=-10, block_dy=5, ): np.random.seed(1) # Initiate I/O class for DC IO = DC.IO() # Obtain ABMN locations xmin, xmax = 0.0, 200.0 ymin, ymax = 0.0, 0.0 zmin, zmax = 0, 0 endl = np.array([[xmin, ymin, zmin], [xmax, ymax, zmax]]) # Generate DC survey object survey = DCutils.generate_dcip_survey( endl, survey_type=survey_type, dim=2, a=10, b=10, n=10 ) survey = IO.from_abmn_locations_to_survey( survey.locations_a, survey.locations_b, survey.locations_m, survey.locations_n, survey_type, data_dc_type="volt", ) # Obtain 2D TensorMesh mesh, actind = IO.set_mesh() # Flat topography actind = active_from_xyz( mesh, np.c_[mesh.cell_centers_x, mesh.cell_centers_x * 0.0] ) survey.drape_electrodes_on_topography(mesh, actind, option="top") # Use Exponential Map: m = log(rho) parametric_block = maps.ParametricBlock(mesh, slopeFact=1e2) mapping = maps.ExpMap(mesh) * parametric_block # Set true model # val_background,val_block, block_x0, block_dx, block_y0, block_dy mtrue = np.r_[np.log(1e3), np.log(10), 100, 10, -20, 10] # Set initial model m0 = np.r_[ np.log(rho_background), np.log(rho_block), block_x0, block_dx, block_y0, block_dy, ] rho = mapping * mtrue rho0 = mapping * m0 # Show the true conductivity model fig = plt.figure(figsize=(12, 3)) ax = plt.subplot(111) temp = rho.copy() temp[~actind] = np.nan out = mesh.plot_image( temp, grid=False, ax=ax, grid_opts={"alpha": 0.2}, pcolor_opts={"cmap": "viridis", "norm": colors.LogNorm(10, 1000)}, ) ax.plot( survey.unique_electrode_locations[:, 0], survey.unique_electrode_locations[:, 1], "k.", ) ax.set_xlim(IO.grids[:, 0].min(), IO.grids[:, 0].max()) ax.set_ylim(-IO.grids[:, 1].max(), IO.grids[:, 1].min()) cb = plt.colorbar(out[0]) cb.set_label("Resistivity (ohm-m)") ax.set_aspect("equal") ax.set_title("True resistivity model") plt.show() # Show the true conductivity model fig = plt.figure(figsize=(12, 3)) ax = plt.subplot(111) temp = rho0.copy() temp[~actind] = np.nan out = mesh.plot_image( temp, grid=False, ax=ax, grid_opts={"alpha": 0.2}, pcolor_opts={"cmap": "viridis", "norm": colors.LogNorm(10, 1000)}, ) ax.plot( survey.unique_electrode_locations[:, 0], survey.unique_electrode_locations[:, 1], "k.", ) ax.set_xlim(IO.grids[:, 0].min(), IO.grids[:, 0].max()) ax.set_ylim(-IO.grids[:, 1].max(), IO.grids[:, 1].min()) cb = plt.colorbar(out[0]) cb.set_label("Resistivity (ohm-m)") ax.set_aspect("equal") ax.set_title("Initial resistivity model") plt.show() # Generate 2.5D DC problem # "N" means potential is defined at nodes prb = DC.Simulation2DNodal( mesh, survey=survey, rhoMap=mapping, storeJ=True, solver=Solver ) # Make synthetic DC data with 5% Gaussian noise data = prb.make_synthetic_data(mtrue, relative_error=0.05, add_noise=True) # Show apparent resisitivty pseudo-section IO.plotPseudoSection(data=data.dobs / IO.G, data_type="apparent_resistivity") # Show apparent resisitivty histogram fig = plt.figure() out = hist(data.dobs / IO.G, bins=20) plt.show() # Set standard_deviation # floor eps = 10 ** (-3.2) # percentage relative = 0.05 dmisfit = data_misfit.L2DataMisfit(simulation=prb, data=data) uncert = abs(data.dobs) * relative + eps dmisfit.standard_deviation = uncert # Map for a regularization mesh_1d = TensorMesh([parametric_block.nP]) # Related to inversion reg = regularization.WeightedLeastSquares(mesh_1d, alpha_x=0.0) opt = optimization.InexactGaussNewton(maxIter=10) invProb = inverse_problem.BaseInvProblem(dmisfit, reg, opt) target = directives.TargetMisfit() invProb.beta = 0.0 inv = inversion.BaseInversion(invProb, directiveList=[target]) prb.counter = opt.counter = utils.Counter() opt.LSshorten = 0.5 opt.remember("xc") # Run inversion mopt = inv.run(m0) # Convert obtained inversion model to resistivity # rho = M(m), where M(.) is a mapping rho_est = mapping * mopt rho_true = rho.copy() # show recovered conductivity fig, ax = plt.subplots(2, 1, figsize=(20, 6)) out1 = mesh.plot_image( rho_true, pcolor_opts={"cmap": "viridis", "norm": colors.LogNorm(10, 1000)}, ax=ax[0], ) out2 = mesh.plot_image( rho_est, pcolor_opts={"cmap": "viridis", "norm": colors.LogNorm(10, 1000)}, ax=ax[1], ) out = [out1, out2] for i in range(2): ax[i].plot( survey.unique_electrode_locations[:, 0], survey.unique_electrode_locations[:, 1], "kv", ) ax[i].set_xlim(IO.grids[:, 0].min(), IO.grids[:, 0].max()) ax[i].set_ylim(-IO.grids[:, 1].max(), IO.grids[:, 1].min()) cb = plt.colorbar(out[i][0], ax=ax[i]) cb.set_label(r"Resistivity ($\Omega$m)") ax[i].set_xlabel("Northing (m)") ax[i].set_ylabel("Elevation (m)") ax[i].set_aspect("equal") ax[0].set_title("True resistivity model") ax[1].set_title("Recovered resistivity model") plt.tight_layout() plt.show() if __name__ == "__main__": run() plt.show() .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 23.439 seconds) **Estimated memory usage:** 9 MB .. _sphx_glr_download_content_examples_04-dcip_plot_inv_dcip_dipoledipole_parametric_inversion.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_dcip_dipoledipole_parametric_inversion.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_inv_dcip_dipoledipole_parametric_inversion.py ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_