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.

  • ../../../_images/sphx_glr_plot_dipoledipole_parametric_inversion_001.png
  • ../../../_images/sphx_glr_plot_dipoledipole_parametric_inversion_002.png
  • ../../../_images/sphx_glr_plot_dipoledipole_parametric_inversion_003.png
  • ../../../_images/sphx_glr_plot_dipoledipole_parametric_inversion_004.png
  • ../../../_images/sphx_glr_plot_dipoledipole_parametric_inversion_005.png

Out:

Compute fields
SimPEG.DataMisfit.l2_DataMisfit assigning default eps of 1e-5 * ||dobs||
SimPEG.InvProblem will set Regularization.mref to m0.

    SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.
    ***Done using same Solver and solverOpts as the 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
Compute fields
   0  0.00e+00  2.74e+03  0.00e+00  2.74e+03    9.20e+03      0
Compute fields
   1  0.00e+00  1.86e+03  1.76e-02  1.86e+03    8.75e+02      0
Compute fields
   2  0.00e+00  9.05e+02  3.30e+00  9.05e+02    2.82e+03      0   Skip BFGS
Compute fields
   3  0.00e+00  6.06e+02  2.11e+00  6.06e+02    2.23e+02      0
Compute fields
Compute fields
Compute fields
Compute fields
Compute fields
Compute fields
   4  0.00e+00  4.42e+02  8.59e+00  4.42e+02    2.05e+03      5   Skip BFGS
Compute fields
Compute fields
Compute fields
   5  0.00e+00  4.24e+02  6.10e+00  4.24e+02    3.47e+02      2
Compute fields
Compute fields
Compute fields
Compute fields
Compute fields
Compute fields
Compute fields
   6  0.00e+00  3.38e+02  1.02e+01  3.38e+02    1.54e+03      6
Compute fields
Compute fields
Compute fields
   7  0.00e+00  3.27e+02  7.95e+00  3.27e+02    1.37e+03      2
Compute fields
Compute fields
Compute fields
Compute fields
Compute fields
Compute fields
Compute fields
   8  0.00e+00  3.26e+02  9.82e+00  3.26e+02    1.37e+03      6
Compute fields
Compute fields
Compute fields
Compute fields
   9  0.00e+00  2.32e+02  3.43e+01  2.32e+02    7.60e+02      3   Skip BFGS
Compute fields
Compute fields
Compute fields
Compute fields
Compute fields
Compute fields
Compute fields
  10  0.00e+00  1.67e+02  1.54e+01  1.67e+02    2.11e+02      6
------------------------- STOP! -------------------------
1 : |fc-fOld| = 6.5617e+01 <= tolF*(1+|f0|) = 2.7367e+02
1 : |xc-x_last| = 3.7212e+00 <= tolX*(1+|x0|) = 1.0246e+01
0 : |proj(x-g)-x|    = 2.1067e+02 <= tolG          = 1.0000e-01
0 : |proj(x-g)-x|    = 2.1067e+02 <= 1e3*eps       = 1.0000e-02
1 : maxIter   =      10    <= iter          =     10
------------------------- DONE! -------------------------

from SimPEG import DC, Mesh
from SimPEG import (Maps, Utils, DataMisfit, Regularization,
                    Optimization, Inversion, InvProblem, 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., 200.
    ymin, ymax = 0., 0.
    zmin, zmax = 0, 0
    endl = np.array([[xmin, ymin, zmin], [xmax, ymax, zmax]])
    # Generate DC survey object
    survey = DC.Utils.gen_DCIPsurvey(endl, survey_type=survey_type, dim=2,
                                     a=10, b=10, n=10)
    survey.getABMN_locations()
    survey = IO.from_ambn_locations_to_survey(
        survey.a_locations, survey.b_locations,
        survey.m_locations, survey.n_locations,
        survey_type, data_dc_type='volt'
    )

    # Obtain 2D TensorMesh
    mesh, actind = IO.set_mesh()
    # Flat topography
    actind = Utils.surface2ind_topo(
        mesh, np.c_[mesh.vectorCCx, mesh.vectorCCx*0.]
    )
    survey.drapeTopo(mesh, actind, option="top")
    # Use Exponential Map: m = log(rho)
    actmap = Maps.InjectActiveCells(
        mesh, indActive=actind, valInactive=np.log(1e8)
    )
    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.plotImage(
        temp, grid=False, ax=ax, gridOpts={'alpha': 0.2},
        clim=(10, 1000),
        pcolorOpts={"cmap": "viridis", "norm": colors.LogNorm()}
    )
    ax.plot(
        survey.electrode_locations[:, 0],
        survey.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.plotImage(
        temp, grid=False, ax=ax, gridOpts={'alpha': 0.2},
        clim=(10, 1000),
        pcolorOpts={"cmap": "viridis", "norm": colors.LogNorm()}
    )
    ax.plot(
        survey.electrode_locations[:, 0],
        survey.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.Problem2D_N(
        mesh, rhoMap=mapping, storeJ=True,
        Solver=Solver
    )
    # Pair problem with survey
    try:
        prb.pair(survey)
    except:
        survey.unpair()
        prb.pair(survey)

    # Make synthetic DC data with 5% Gaussian noise
    dtrue = survey.makeSyntheticData(mtrue, std=0.05, force=True)

    # Show apparent resisitivty pseudo-section
    IO.plotPseudoSection(
        data=survey.dobs/IO.G, data_type='apparent_resistivity'
    )

    # Show apparent resisitivty histogram
    fig = plt.figure()
    out = hist(survey.dobs/IO.G, bins=20)
    plt.show()
    # Set uncertainty
    # floor
    eps = 10**(-3.2)
    # percentage
    std = 0.05
    dmisfit = DataMisfit.l2_DataMisfit(survey)
    uncert = abs(survey.dobs) * std + eps
    dmisfit.W = 1./uncert

    # Map for a regularization
    mesh_1d = Mesh.TensorMesh([parametric_block.nP])
    # Related to inversion
    reg = Regularization.Simple(mesh_1d, alpha_x=0.)
    opt = Optimization.InexactGaussNewton(maxIter=10)
    invProb = InvProblem.BaseInvProblem(dmisfit, reg, opt)
    beta = Directives.BetaSchedule(coolingFactor=5, coolingRate=2)
    betaest = Directives.BetaEstimate_ByEig(beta0_ratio=1e0)
    target = Directives.TargetMisfit()
    updateSensW = Directives.UpdateSensitivityWeights()
    update_Jacobi = Directives.UpdatePreconditioner()
    invProb.beta = 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
    vmin, vmax = rho.min(), rho.max()
    fig, ax = plt.subplots(2, 1, figsize=(20, 6))
    out1 = mesh.plotImage(
            rho_true, clim=(10, 1000),
            pcolorOpts={"cmap": "viridis", "norm": colors.LogNorm()},
            ax=ax[0]
    )
    out2 = mesh.plotImage(
        rho_est, clim=(10, 1000),
        pcolorOpts={"cmap": "viridis", "norm": colors.LogNorm()},
        ax=ax[1]
    )
    out = [out1, out2]
    for i in range(2):
        ax[i].plot(
            survey.electrode_locations[:, 0],
            survey.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("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()

Total running time of the script: ( 0 minutes 48.839 seconds)

Gallery generated by Sphinx-Gallery