2.5D DC inversion of with Iterative Reweighted Least SquaresΒΆ

This is an example for 2.5D DC Inversion with Iterative Reweighted Least Squares (IRLS). Earth includes a topography, and below the topography conductive and resistive cylinders are embedded. User is promoted to try different p, qx, qz. For instance a set of paraemters (default):

  • p=0 (sparse model, m)
  • qx=2 (smooth model, m in x-direction)
  • qz=2 (smooth model, m in z-direction)

But if you want share edges of the model, you can try:

  • p=0 (sparse model, m)
  • qx=0 (smooth model, m in x-direction)
  • qz=2 (smooth model, m in z-direction)
  • ../../../_images/sphx_glr_plot_dipoledipole_2_5Dinversion_irls_001.png
  • ../../../_images/sphx_glr_plot_dipoledipole_2_5Dinversion_irls_002.png
  • ../../../_images/sphx_glr_plot_dipoledipole_2_5Dinversion_irls_003.png
  • ../../../_images/sphx_glr_plot_dipoledipole_2_5Dinversion_irls_004.png

Out:

Using a seed of:  37
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
   0  4.26e+00  4.59e+03  0.00e+00  4.59e+03    9.35e+02      0
   1  1.32e+00  5.52e+02  8.04e+01  6.59e+02    1.22e+02      0
   2  7.21e-01  1.15e+02  1.87e+02  2.50e+02    2.35e+01      0   Skip BFGS
Reached starting chifact with l2-norm regularization: Start IRLS steps...
eps_p: 2.29859284353 eps_q: 2.29859284353
Eps_p: 1.91549403627
Eps_q: 1.91549403627
delta phim:    inf
   3  7.21e-01  6.23e+01  3.90e+02  3.44e+02    1.77e+01      0   Skip BFGS
Beta search step
   4  4.19e-01  6.23e+01  3.90e+02  2.26e+02    1.11e+01      0
Beta search step
   5  7.00e-01  6.23e+01  3.90e+02  3.36e+02    1.70e+01      0   Skip BFGS
Beta search step
   6  1.09e+00  6.23e+01  3.90e+02  4.88e+02    3.08e+01      0
Beta search step
   7  7.35e-01  6.23e+01  3.90e+02  3.49e+02    1.81e+01      0
Beta search step
   8  1.18e+00  6.23e+01  3.90e+02  5.22e+02    3.41e+01      0   Skip BFGS
Beta search step
   9  6.59e-01  6.23e+01  3.90e+02  3.20e+02    1.58e+01      0
Beta search step
  10  1.07e+00  6.23e+01  3.90e+02  4.81e+02    3.01e+01      0   Skip BFGS
Beta search step
  11  6.60e-01  6.23e+01  3.90e+02  3.20e+02    1.58e+01      0
Beta search step
  12  1.07e+00  6.23e+01  3.90e+02  4.81e+02    3.02e+01      0   Skip BFGS
Beta search step
  13  6.73e-01  6.23e+01  3.90e+02  3.25e+02    1.62e+01      0
Beta search step
  14  1.09e+00  6.23e+01  3.90e+02  4.89e+02    3.09e+01      0   Skip BFGS
Beta search step
  15  6.77e-01  6.23e+01  3.90e+02  3.27e+02    1.63e+01      0
Beta search step
  16  1.10e+00  6.23e+01  3.90e+02  4.91e+02    3.11e+01      0   Skip BFGS
Beta search step
  17  6.76e-01  6.23e+01  3.90e+02  3.26e+02    1.63e+01      0
Beta search step
  18  1.10e+00  6.23e+01  3.90e+02  4.90e+02    3.11e+01      0   Skip BFGS
Beta search step
  19  6.77e-01  6.23e+01  3.90e+02  3.26e+02    1.63e+01      0
Beta search step
  20  1.10e+00  6.23e+01  3.90e+02  4.91e+02    3.11e+01      0   Skip BFGS
Beta search step
  21  6.77e-01  6.23e+01  3.90e+02  3.26e+02    1.63e+01      0
Beta search step
  22  1.10e+00  6.23e+01  3.90e+02  4.91e+02    3.11e+01      0   Skip BFGS
Beta search step
  23  6.76e-01  6.23e+01  3.90e+02  3.26e+02    1.63e+01      0
Beta search step
  24  1.10e+00  6.23e+01  3.90e+02  4.93e+02    3.13e+01      0   Skip BFGS
Beta search step
  25  6.59e-01  6.23e+01  3.90e+02  3.19e+02    1.58e+01      0
Beta search step
  26  1.08e+00  6.23e+01  3.90e+02  4.85e+02    3.05e+01      0   Skip BFGS
Beta search step
  27  7.11e-01  6.23e+01  3.90e+02  3.40e+02    1.74e+01      0
Beta search step
  28  1.16e+00  6.23e+01  3.90e+02  5.17e+02    3.37e+01      0   Skip BFGS
Beta search step
  29  6.56e-01  6.23e+01  3.90e+02  3.18e+02    1.57e+01      0
Beta search step
  30  1.07e+00  6.23e+01  3.90e+02  4.81e+02    3.02e+01      0   Skip BFGS
Beta search step
  31  7.10e-01  6.23e+01  3.90e+02  3.40e+02    1.74e+01      0
Beta search step
  32  1.16e+00  6.23e+01  3.90e+02  5.17e+02    3.36e+01      0   Skip BFGS
Beta search step
  33  6.58e-01  6.23e+01  3.90e+02  3.19e+02    1.58e+01      0
Beta search step
  34  1.08e+00  6.23e+01  3.90e+02  4.83e+02    3.04e+01      0   Skip BFGS
Beta search step
  35  7.09e-01  6.23e+01  3.90e+02  3.39e+02    1.73e+01      0
Beta search step
  36  1.16e+00  6.23e+01  3.90e+02  5.16e+02    3.36e+01      0   Skip BFGS
Beta search step
  37  6.58e-01  6.23e+01  3.90e+02  3.19e+02    1.58e+01      0
Beta search step
  38  1.08e+00  6.23e+01  3.90e+02  4.84e+02    3.05e+01      0   Skip BFGS
Beta search step
  39  7.05e-01  6.23e+01  3.90e+02  3.38e+02    1.72e+01      0
Beta search step
  40  1.16e+00  6.23e+01  3.90e+02  5.14e+02    3.34e+01      0   Skip BFGS
------------------------- STOP! -------------------------
1 : |fc-fOld| = 1.7663e+02 <= tolF*(1+|f0|) = 4.5906e+02
1 : |xc-x_last| = 0.0000e+00 <= tolX*(1+|x0|) = 2.6208e+01
0 : |proj(x-g)-x|    = 3.3382e+01 <= tolG          = 1.0000e-01
0 : |proj(x-g)-x|    = 3.3382e+01 <= 1e3*eps       = 1.0000e-02
1 : maxIter   =      40    <= iter          =     40
------------------------- DONE! -------------------------

from SimPEG import DC
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 Pardiso as Solver
except ImportError:
    from SimPEG import SolverLU as Solver


def run(plotIt=True, survey_type="dipole-dipole", p=0., qx=2., qz=2.):
    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()
    topo, mesh1D = DC.Utils.genTopography(mesh, -10, 0, its=100)
    actind = Utils.surface2ind_topo(mesh, np.c_[mesh1D.vectorCCx, topo])
    survey.drapeTopo(mesh, actind, option="top")

    # Build a conductivity 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
    )
    layer_inds = mesh.gridCC[:, 1] > -5.
    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

    # Show the true conductivity model
    if plotIt:
        fig = plt.figure(figsize=(12, 3))
        ax = plt.subplot(111)
        temp = rho.copy()
        temp[~actind] = np.nan
        out = mesh.plotImage(
            temp, grid=True, 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')
        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
    mtrue = 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)
    except:
        survey.unpair()
        prb.pair(survey)

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

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

    # Show apparent resisitivty histogram
    if plotIt:
        fig = plt.figure()
        out = hist(survey.dobs/IO.G, bins=20)
        plt.xlabel("Apparent Resisitivty ($\Omega$m)")
        plt.show()

    # Set initial model based upon histogram
    m0 = np.ones(actmap.nP)*np.log(100.)

    # 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
    regmap = Maps.IdentityMap(nP=int(actind.sum()))

    # Related to inversion
    reg = Regularization.Sparse(
        mesh, indActive=actind, mapping=regmap,
        gradientType = 'components'
    )
    #     gradientType = 'components'
    reg.norms = np.c_[p, qx, qz, 0.]
    IRLS = Directives.Update_IRLS(maxIRLSiter=20, minGNiter=1)

    opt = Optimization.InexactGaussNewton(maxIter=40)
    invProb = InvProblem.BaseInvProblem(dmisfit, reg, opt)
    beta = Directives.BetaSchedule(coolingFactor=5, coolingRate=2)
    betaest = Directives.BetaEstimate_ByEig(beta0_ratio=1e0)
    target = Directives.TargetMisfit()
    update_Jacobi = Directives.UpdatePreconditioner()
    inv = Inversion.BaseInversion(
        invProb, directiveList=[
            betaest, IRLS
        ]
        )
    prb.counter = opt.counter = Utils.Counter()
    opt.LSshorten = 0.5
    opt.remember('xc')

    # Run inversion
    mopt = inv.run(m0)

    rho_est = mapping*mopt
    rho_est_l2 = mapping*invProb.l2model
    rho_est[~actind] = np.nan
    rho_est_l2[~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(3, 1, figsize=(20, 9))
        out1 = mesh.plotImage(
                rho_true, clim=(10, 1000),
                pcolorOpts={"cmap": "viridis", "norm": colors.LogNorm()},
                ax=ax[0]
        )
        out2 = mesh.plotImage(
            rho_est_l2, clim=(10, 1000),
            pcolorOpts={"cmap": "viridis", "norm": colors.LogNorm()},
            ax=ax[1]
        )
        out3 = mesh.plotImage(
            rho_est, clim=(10, 1000),
            pcolorOpts={"cmap": "viridis", "norm": colors.LogNorm()},
            ax=ax[2]
        )

        out = [out1, out2, out3]
        titles = ["True", "L2", ("L%d, Lx%d, Lz%d")%(p, qx, qz)]
        for i in range(3):
            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[i].set_title(titles[i])
        plt.tight_layout()
        plt.show()


if __name__ == '__main__':
    run()

Total running time of the script: ( 2 minutes 33.688 seconds)

Gallery generated by Sphinx-Gallery