EM: FDEM: 1D: InversionΒΆ

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

(Source code)

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
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):
    """
        EM: FDEM: 1D: Inversion
        =======================

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

    """

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