2.5D DC-IP inversion of Dipole Dipole array with TopographyΒΆ

This is an example for 2.5D DC-IP Inversion. For DC inversion, a resisistivity model (Ohm-m) is generated having conductive and resistive cylinders; they are respectively located right and left sides of the subsurface. For IP inversion, a chargeability model (V/V) is generated having a chargeable cyinder located at the center. Default survey_type is dipole-dipole, but this can be changed. to ‘pole-dipole’, ‘dipole-pole’, and ‘pole-pole’. By running DC and IP simulations synthetic DC and IP data are generated, respectively. Following two-stage approach (Oldenburg et al, 1999), first DC data is inverted to recover a resistivity model. Then by using the obtained resistivity model, sensitivity function is formed and used for subsequent IP inversion to recover a chargeability model.

  • ../../../_images/sphx_glr_plot_dcip_2_5Dinversion_001.png
  • ../../../_images/sphx_glr_plot_dcip_2_5Dinversion_002.png
  • ../../../_images/sphx_glr_plot_dcip_2_5Dinversion_003.png
  • ../../../_images/sphx_glr_plot_dcip_2_5Dinversion_004.png
  • ../../../_images/sphx_glr_plot_dcip_2_5Dinversion_005.png
  • ../../../_images/sphx_glr_plot_dcip_2_5Dinversion_006.png
  • ../../../_images/sphx_glr_plot_dcip_2_5Dinversion_007.png

Out:

/Users/lindseyjh/git/python_symlinks/SimPEG/EM/Static/DC/IODC.py:232: UserWarning: code under construction - API might change in the future
  "code under construction - API might change in the future"
/Users/lindseyjh/git/python_symlinks/SimPEG/EM/Static/DC/IODC.py:504: UserWarning: dx is set to 2.5 m (samllest electrode spacing (10.0) / 4)
  "dx is set to {} m (samllest electrode spacing ({}) / {})".format(dx, a, ncell_per_dipole)
/Users/lindseyjh/git/python_symlinks/SimPEG/EM/Static/DC/IODC.py:509: UserWarning: dz (1.25 m) is set to dx (2.5 m) / 2
  "dz ({} m) is set to dx ({} m) / {}".format(dz, dx, 2)
Using a seed of:  37
/Users/lindseyjh/git/simpeg/simpeg/examples/12-ip/plot_dcip_2_5Dinversion.py:107: UserWarning: Matplotlib is currently using agg, which is a non-GUI backend, so cannot show the figure.
  plt.show()
SimPEG.Survey assigned new std of 5.00%
SimPEG.Survey assigned new std of 5.00%
/Users/lindseyjh/git/simpeg/simpeg/examples/12-ip/plot_dcip_2_5Dinversion.py:156: UserWarning: Matplotlib is currently using agg, which is a non-GUI backend, so cannot show the figure.
  plt.show()
/Users/lindseyjh/git/simpeg/simpeg/examples/12-ip/plot_dcip_2_5Dinversion.py:164: UserWarning: Matplotlib is currently using agg, which is a non-GUI backend, so cannot show the figure.
  plt.show()
/Users/lindseyjh/git/simpeg/simpeg/examples/12-ip/plot_dcip_2_5Dinversion.py:176: UserWarning: Matplotlib is currently using agg, which is a non-GUI backend, so cannot show the figure.
  plt.show()
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
=============================== Projected GNCG ===============================
  #     beta     phi_d     phi_m       f      |proj(x-g)-x|  LS    Comment
-----------------------------------------------------------------------------
x0 has any nan: 0
   0  2.31e+00  4.59e+03  0.00e+00  4.59e+03    9.35e+02      0
   1  2.31e+00  1.76e+02  6.08e-01  1.77e+02    1.05e+02      0
------------------------- STOP! -------------------------
1 : |fc-fOld| = 0.0000e+00 <= tolF*(1+|f0|) = 4.5906e+02
1 : |xc-x_last| = 1.5091e+01 <= tolX*(1+|x0|) = 2.6208e+01
0 : |proj(x-g)-x|    = 1.0550e+02 <= tolG          = 1.0000e-01
0 : |proj(x-g)-x|    = 1.0550e+02 <= 1e3*eps       = 1.0000e-02
0 : maxIter   =      15    <= iter          =      2
------------------------- DONE! -------------------------
/Users/lindseyjh/git/simpeg/simpeg/examples/12-ip/plot_dcip_2_5Dinversion.py:228: UserWarning: Matplotlib is currently using agg, which is a non-GUI backend, so cannot show the figure.
  plt.show()
/Users/lindseyjh/git/simpeg/simpeg/examples/12-ip/plot_dcip_2_5Dinversion.py:240: UserWarning: Matplotlib is currently using agg, which is a non-GUI backend, so cannot show the figure.
  plt.show()
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***
Approximated diag(JtJ) with linear operator
model has any nan: 0
=============================== Projected GNCG ===============================
  #     beta     phi_d     phi_m       f      |proj(x-g)-x|  LS    Comment
-----------------------------------------------------------------------------
x0 has any nan: 0
   0  2.78e+04  6.80e+02  0.00e+00  6.80e+02    1.45e+04      0
   1  2.78e+04  4.28e+02  5.61e-04  4.43e+02    2.55e+03      0
   2  5.57e+03  3.26e+02  1.80e-03  3.36e+02    3.03e+03      0   Skip BFGS
   3  5.57e+03  2.96e+02  5.82e-03  3.28e+02    7.79e+03      0   Skip BFGS
   4  1.11e+03  1.97e+02  6.12e-03  2.04e+02    4.05e+03      0
   5  1.11e+03  1.51e+02  1.31e-02  1.65e+02    3.45e+03      1
   6  2.23e+02  1.04e+02  1.85e-02  1.08e+02    3.92e+03      0
   7  2.23e+02  6.96e+01  5.74e-02  8.24e+01    6.68e+03      0
------------------------- STOP! -------------------------
1 : |fc-fOld| = 0.0000e+00 <= tolF*(1+|f0|) = 6.8057e+01
1 : |xc-x_last| = 8.0313e-02 <= tolX*(1+|x0|) = 1.0000e-01
0 : |proj(x-g)-x|    = 6.6764e+03 <= tolG          = 1.0000e-01
0 : |proj(x-g)-x|    = 6.6764e+03 <= 1e3*eps       = 1.0000e-02
0 : maxIter   =      15    <= iter          =      8
------------------------- DONE! -------------------------
/Users/lindseyjh/git/simpeg/simpeg/examples/12-ip/plot_dcip_2_5Dinversion.py:296: UserWarning: Matplotlib is currently using agg, which is a non-GUI backend, so cannot show the figure.
  plt.show()

from SimPEG import DC, IP
from SimPEG import Maps, Utils
import matplotlib.pyplot as plt
from matplotlib import colors
import numpy as np
from pylab import hist
try:
    from pymatsolver import Pardiso as Solver
except ImportError:
    from SimPEG import SolverLU as Solver


def run(plotIt=True, survey_type="dipole-dipole"):
    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 = DC.Utils.gen_DCIPsurvey(endl, survey_type=survey_type, dim=2,
                                     a=10, b=10, n=10)
    survey_dc.getABMN_locations()
    survey_dc = IO.from_ambn_locations_to_survey(
        survey_dc.a_locations, survey_dc.b_locations,
        survey_dc.m_locations, survey_dc.n_locations,
        survey_type, data_dc_type='volt', data_ip_type='volt'
    )

    # Obtain 2D TensorMesh
    mesh, actind = IO.set_mesh()
    topo, mesh1D = DC.Utils.genTopography(mesh, -10, 0, its=100)
    actind = Utils.surface2ind_topo(mesh, np.c_[mesh1D.vectorCCx, topo])
    survey_dc.drapeTopo(mesh, actind, option="top")

    # Build conductivity and chargeability model
    blk_inds_c = Utils.ModelBuilder.getIndicesSphere(
        np.r_[60., -25.], 12.5, mesh.gridCC
    )
    blk_inds_r = Utils.ModelBuilder.getIndicesSphere(
        np.r_[140., -25.], 12.5, mesh.gridCC
    )
    blk_inds_charg = Utils.ModelBuilder.getIndicesSphere(
        np.r_[100., -25], 12.5, mesh.gridCC
    )
    sigma = np.ones(mesh.nC)*1./100.
    sigma[blk_inds_c] = 1./10.
    sigma[blk_inds_r] = 1./1000.
    sigma[~actind] = 1./1e8
    rho = 1./sigma
    charg = np.zeros(mesh.nC)
    charg[blk_inds_charg] = 0.1

    # Show the true conductivity model
    if plotIt:
        fig, axs = plt.subplots(2,1, figsize=(12, 6))
        temp_rho = rho.copy()
        temp_rho[~actind] = np.nan
        temp_charg = charg.copy()
        temp_charg[~actind] = np.nan

        out1 = mesh.plotImage(
            temp_rho, grid=True, ax=axs[0], gridOpts={'alpha': 0.2},
            clim=(10, 1000),
            pcolorOpts={"cmap": "viridis", "norm": colors.LogNorm()}
        )
        out2 = mesh.plotImage(
            temp_charg, grid=True, ax=axs[1], gridOpts={'alpha': 0.2},
            clim=(0, 0.1),
            pcolorOpts={"cmap": "magma"}
        )
        for i in range(2):
            axs[i].plot(
                survey_dc.electrode_locations[:, 0],
                survey_dc.electrode_locations[:, 1], 'kv'
            )
            axs[i].set_xlim(IO.grids[:, 0].min(), IO.grids[:, 0].max())
            axs[i].set_ylim(-IO.grids[:, 1].max(), IO.grids[:, 1].min())
            axs[i].set_aspect('equal')
        cb = plt.colorbar(out1[0], ax=axs[0])
        cb.set_label("Resistivity (ohm-m)")
        cb = plt.colorbar(out2[0], ax=axs[1])
        cb.set_label("Chargeability")

        plt.show()

    # Use Exponential Map: m = log(rho)
    actmap = Maps.InjectActiveCells(
        mesh, indActive=actind, valInactive=np.log(1e8)
    )
    mapping = Maps.ExpMap(mesh) * actmap

    # Generate mtrue_dc for resistivity
    mtrue_dc = np.log(rho[actind])

    # 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_dc)
    except:
        survey_dc.unpair()
        prb.pair(survey_dc)

    # Make synthetic DC data with 5% Gaussian noise
    dtrue_dc = survey_dc.makeSyntheticData(mtrue_dc, std=0.05, force=True)
    IO.data_dc = dtrue_dc

    # Generate mtrue_ip for chargability
    mtrue_ip = charg[actind]
    # Generate 2.5D DC problem
    # "N" means potential is defined at nodes
    prb_ip = IP.Problem2D_N(
        mesh, etaMap=actmap, storeJ=True, rho=rho,
        Solver=Solver
    )
    survey_ip = IP.from_dc_to_ip_survey(survey_dc, dim="2.5D")
    prb_ip.pair(survey_ip)

    dtrue_ip = survey_ip.makeSyntheticData(mtrue_ip, std=0.05)

    IO.data_ip = dtrue_ip

    # Show apparent resisitivty pseudo-section
    if plotIt:
        IO.plotPseudoSection(
            data_type='apparent_resistivity', scale='log',
            cmap='viridis'
        )
        plt.show()

    # Show apparent chargeability pseudo-section
    if plotIt:
        IO.plotPseudoSection(
            data_type='apparent_chargeability', scale='linear',
            cmap='magma'
        )
        plt.show()

    # Show apparent resisitivty histogram
    if plotIt:
        fig = plt.figure(figsize=(10, 4))
        ax1 = plt.subplot(121)
        out = hist(np.log10(abs(IO.voltages)), bins=20)
        ax1.set_xlabel("log10 DC voltage (V)")
        ax2 = plt.subplot(122)
        out = hist(IO.apparent_resistivity, bins=20)
        ax2.set_xlabel("Apparent Resistivity ($\Omega$m)")
        plt.tight_layout()
        plt.show()

    # Set initial model based upon histogram
    m0_dc = np.ones(actmap.nP)*np.log(100.)
    # Set uncertainty
    # floor
    eps_dc = 10**(-3.2)
    # percentage
    std_dc = 0.05

    mopt_dc, pred_dc = DC.run_inversion(
        m0_dc, survey_dc, actind, mesh, std_dc, eps_dc,
        beta0_ratio=1e0,
        use_sensitivity_weight=True
        )

    # Convert obtained inversion model to resistivity
    # rho = M(m), where M(.) is a mapping

    rho_est = mapping*mopt_dc
    rho_est[~actind] = np.nan
    rho_true = rho.copy()
    rho_true[~actind] = np.nan

    # show recovered conductivity
    if plotIt:
        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_dc.electrode_locations[:, 0],
                survey_dc.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')
        plt.tight_layout()
        plt.show()

    # Show apparent resisitivty histogram
    if plotIt:
        fig = plt.figure(figsize=(10, 4))
        ax1 = plt.subplot(121)
        out = hist(np.log10(abs(IO.voltages_ip)), bins=20)
        ax1.set_xlabel("log10 IP voltage (V)")
        ax2 = plt.subplot(122)
        out = hist(IO.apparent_chargeability, bins=20)
        ax2.set_xlabel("Apparent Chargeability (V/V)")
        plt.tight_layout()
        plt.show()


    # Set initial model based upon histogram
    m0_ip = np.ones(actmap.nP)*1e-10
    # Set uncertainty
    # floor
    eps_ip = 10**(-4)
    # percentage
    std_ip = 0.05
    # Clean sensitivity function formed with true resistivity
    prb_ip._Jmatrix = None
    # Input obtained resistivity to form sensitivity
    prb_ip.rho = mapping*mopt_dc
    mopt_ip, _ = IP.run_inversion(
        m0_ip, survey_ip, actind, mesh, std_ip, eps_ip,
        upper=np.Inf, lower=0.,
        beta0_ratio=1e0,
        use_sensitivity_weight=True
    )

    # Convert obtained inversion model to chargeability
    # charg = M(m), where M(.) is a mapping for cells below topography

    charg_est = actmap*mopt_ip
    charg_est[~actind] = np.nan
    charg_true = charg.copy()
    charg_true[~actind] = np.nan

    # show recovered chargeability
    if plotIt:
        fig, ax = plt.subplots(2, 1, figsize=(20, 6))
        out1 = mesh.plotImage(
                charg_true, clim=(0, 0.1),
                pcolorOpts={"cmap": "magma"},
                ax=ax[0]
        )
        out2 = mesh.plotImage(
            charg_est, clim=(0, 0.1),
            pcolorOpts={"cmap": "magma"},
            ax=ax[1]
        )
        out = [out1, out2]
        for i in range(2):
            ax[i].plot(
                survey_dc.electrode_locations[:, 0],
                survey_dc.electrode_locations[:, 1], 'rv'
            )
            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')
        plt.tight_layout()
        plt.show()

if __name__ == '__main__':
    survey_type = 'dipole-dipole'
    run(survey_type=survey_type, plotIt=True)

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

Gallery generated by Sphinx-Gallery