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.

• 
• 
• 
• 
• 

/home/vsts/work/1/s/examples/04-dcip/plot_inv_dcip_dipoledipole_parametric_inversion.py:50: 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
/home/vsts/work/1/s/simpeg/base/pde_simulation.py:490: 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.


Running inversion with SimPEG v0.23.1.dev10+gf697d2455
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.
/home/vsts/work/1/s/simpeg/base/pde_simulation.py:490: 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 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
/usr/share/miniconda/envs/simpeg-test/lib/python3.10/site-packages/pymatsolver/solvers.py:415: FutureWarning:

In Future pymatsolver v0.4.0, passing a vector of shape (n, 1) to the solve method will return an array with shape (n, 1), instead of always returning a flattened array. This is to be consistent with numpy.linalg.solve broadcasting.

   0  0.00e+00  7.19e+03  0.00e+00  7.19e+03    3.08e+04      0
   1  0.00e+00  3.85e+03  7.75e-03  3.85e+03    2.76e+03      0
   2  0.00e+00  1.96e+03  7.54e-01  1.96e+03    9.18e+03      0   Skip BFGS
   3  0.00e+00  1.36e+03  7.43e-01  1.36e+03    5.51e+02      0
   4  0.00e+00  1.05e+03  1.62e+00  1.05e+03    5.14e+03      8
   5  0.00e+00  1.01e+03  2.07e+00  1.01e+03    4.87e+03      4
   6  0.00e+00  9.27e+02  1.06e+00  9.27e+02    9.00e+02      1
   7  0.00e+00  4.77e+02  2.47e+02  4.77e+02    6.13e+03      5
   8  0.00e+00  3.95e+02  2.82e+02  3.95e+02    1.73e+03      1
   9  0.00e+00  2.77e+02  2.41e+02  2.77e+02    3.16e+03      2
  10  0.00e+00  1.88e+02  2.64e+02  1.88e+02    8.14e+02      0
------------------------- STOP! -------------------------
1 : |fc-fOld| = 8.9462e+01 <= tolF*(1+|f0|) = 7.1862e+02
1 : |xc-x_last| = 2.1051e+00 <= tolX*(1+|x0|) = 1.0246e+01
0 : |proj(x-g)-x|    = 8.1440e+02 <= tolG          = 1.0000e-01
0 : |proj(x-g)-x|    = 8.1440e+02 <= 1e3*eps       = 1.0000e-02
1 : maxIter   =      10    <= iter          =     10
------------------------- DONE! -------------------------

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


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,
    )

    # 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()

Estimated memory usage: 294 MB