EM: FDEM: 1D: InversionΒΆ

Here we will create and run a FDEM 1D inversion.

  • ../../../_images/sphx_glr_plot_inversion_1D_001.png
  • ../../../_images/sphx_glr_plot_inversion_1D_002.png
  • ../../../_images/sphx_glr_plot_inversion_1D_003.png

Out:

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.20e+03  1.64e+03  0.00e+00  1.64e+03    2.97e+03      0
   1  4.20e+03  9.28e+01  6.73e-03  1.21e+02    3.47e+02      0
   2  8.40e+02  6.48e+00  1.20e-02  1.66e+01    6.42e+01      0   Skip BFGS
   3  8.40e+02  1.64e+00  1.37e-02  1.32e+01    1.02e+00      0   Skip BFGS
   4  1.68e+02  1.63e+00  1.37e-02  3.93e+00    1.14e+01      0   Skip BFGS
   5  1.68e+02  1.45e+00  1.41e-02  3.82e+00    9.15e-01      0
   6  3.36e+01  1.33e+00  1.45e-02  1.82e+00    2.66e+00      0   Skip BFGS
------------------------- STOP! -------------------------
1 : |fc-fOld| = 1.9976e+00 <= tolF*(1+|f0|) = 1.6459e+02
1 : |xc-x_last| = 1.2197e-01 <= tolX*(1+|x0|) = 1.9349e+00
0 : |proj(x-g)-x|    = 2.6643e+00 <= tolG          = 1.0000e-01
0 : |proj(x-g)-x|    = 2.6643e+00 <= 1e3*eps       = 1.0000e-02
1 : maxIter   =       6    <= iter          =      6
------------------------- DONE! -------------------------

import numpy as np
from SimPEG import (Mesh, Maps, Utils, DataMisfit, Regularization,
                    Optimization, Inversion, InvProblem, Directives)
import SimPEG.EM as EM
import matplotlib.pyplot as plt
try:
    from pymatsolver import Pardiso as Solver
except ImportError:
    from SimPEG import SolverLU as Solver


def run(plotIt=True):

    cs, ncx, ncz, npad = 5., 25, 15, 15
    hx = [(cs, ncx), (cs, npad, 1.3)]
    hz = [(cs, npad, -1.3), (cs, ncz), (cs, npad, 1.3)]
    mesh = Mesh.CylMesh([hx, 1, hz], '00C')

    layerz = -100.

    active = mesh.vectorCCz < 0.
    layer = (mesh.vectorCCz < 0.) & (mesh.vectorCCz >= layerz)
    actMap = Maps.InjectActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
    mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * actMap
    sig_half = 2e-2
    sig_air = 1e-8
    sig_layer = 1e-2
    sigma = np.ones(mesh.nCz)*sig_air
    sigma[active] = sig_half
    sigma[layer] = sig_layer
    mtrue = np.log(sigma[active])

    if plotIt:
        fig, ax = plt.subplots(1, 1, figsize=(3, 6))
        plt.semilogx(sigma[active], mesh.vectorCCz[active])
        ax.set_ylim(-500, 0)
        ax.set_xlim(1e-3, 1e-1)
        ax.set_xlabel('Conductivity (S/m)', fontsize=14)
        ax.set_ylabel('Depth (m)', fontsize=14)
        ax.grid(color='k', alpha=0.5, linestyle='dashed', linewidth=0.5)

    rxOffset = 10.
    bzi = EM.FDEM.Rx.Point_b(
        np.array([[rxOffset, 0., 1e-3]]),
        orientation='z',
        component='imag'
    )

    freqs = np.logspace(1, 3, 10)
    srcLoc = np.array([0., 0., 10.])

    srcList = [EM.FDEM.Src.MagDipole([bzi], freq, srcLoc, orientation='Z')
               for freq in freqs]

    survey = EM.FDEM.Survey(srcList)
    prb = EM.FDEM.Problem3D_b(mesh, sigmaMap=mapping, Solver=Solver)
    prb.pair(survey)

    std = 0.05
    survey.makeSyntheticData(mtrue, std)

    survey.std = std
    survey.eps = np.linalg.norm(survey.dtrue)*1e-5

    if plotIt:
        fig, ax = plt.subplots(1, 1, figsize=(6, 6))
        ax.semilogx(freqs, survey.dtrue[:freqs.size], 'b.-')
        ax.semilogx(freqs, survey.dobs[:freqs.size], 'r.-')
        ax.legend(('Noisefree', '$d^{obs}$'), fontsize=16)
        ax.set_xlabel('Time (s)', fontsize=14)
        ax.set_ylabel('$B_z$ (T)', fontsize=16)
        ax.set_xlabel('Time (s)', fontsize=14)
        ax.grid(color='k', alpha=0.5, linestyle='dashed', linewidth=0.5)

    dmisfit = DataMisfit.l2_DataMisfit(survey)
    regMesh = Mesh.TensorMesh([mesh.hz[mapping.maps[-1].indActive]])
    reg = Regularization.Tikhonov(regMesh)
    opt = Optimization.InexactGaussNewton(maxIter=6)
    invProb = InvProblem.BaseInvProblem(dmisfit, reg, opt)

    # Create an inversion object
    beta = Directives.BetaSchedule(coolingFactor=5, coolingRate=2)
    betaest = Directives.BetaEstimate_ByEig(beta0_ratio=1e0)
    inv = Inversion.BaseInversion(invProb, directiveList=[beta, betaest])
    m0 = np.log(np.ones(mtrue.size)*sig_half)
    reg.alpha_s = 1e-3
    reg.alpha_x = 1.
    prb.counter = opt.counter = Utils.Counter()
    opt.LSshorten = 0.5
    opt.remember('xc')

    mopt = inv.run(m0)

    if plotIt:
        fig, ax = plt.subplots(1, 1, figsize=(3, 6))
        plt.semilogx(sigma[active], mesh.vectorCCz[active])
        plt.semilogx(np.exp(mopt), mesh.vectorCCz[active])
        ax.set_ylim(-500, 0)
        ax.set_xlim(1e-3, 1e-1)
        ax.set_xlabel('Conductivity (S/m)', fontsize=14)
        ax.set_ylabel('Depth (m)', fontsize=14)
        ax.grid(color='k', alpha=0.5, linestyle='dashed', linewidth=0.5)
        plt.legend(['$\sigma_{true}$', '$\sigma_{pred}$'], loc='best')



if __name__ == '__main__':
    run()
    plt.show()

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

Generated by Sphinx-Gallery