Heagy et al., 2017 1D FDEM and TDEM inversions#

Here, we perform a 1D inversion using both the frequency and time domain codes. The forward simulations are conducted on a cylindrically symmetric mesh. The source is a point magnetic dipole source.

This example is used in the paper

Lindsey J. Heagy, Rowan Cockett, Seogi Kang, Gudni K. Rosenkjaer, Douglas W. Oldenburg, A framework for simulation and inversion in electromagnetics, Computers & Geosciences, Volume 107, 2017, Pages 1-19, ISSN 0098-3004, http://dx.doi.org/10.1016/j.cageo.2017.06.018.

This example is on figshare: https://doi.org/10.6084/m9.figshare.5035175

This example was updated for SimPEG 0.14.0 on January 31st, 2020 by Joseph Capriotti

(a) Recovered Models, (b) FDEM observed vs. predicted, (c) TDEM observed vs. predicted
min skin depth =  158.11388300841895 max skin depth =  500.0
max x  1267.687908603637 min z  -1242.6879086036365 max z  1242.687908603637
SimPEG.InvProblem will set Regularization.reference_model to m0.
SimPEG.InvProblem will set Regularization.reference_model to m0.
SimPEG.InvProblem will set Regularization.reference_model to m0.

                        SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.
                        ***Done using same Solver, and solver_opts as the Simulation3DMagneticFluxDensity 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.01e+00  7.90e+02  0.00e+00  7.90e+02    5.69e+02      0
   1  4.01e+00  1.98e+02  3.13e+01  3.23e+02    5.12e+02      0
   2  4.01e+00  5.58e+01  3.46e+01  1.94e+02    9.49e+01      0
   3  1.00e+00  5.43e+01  3.20e+01  8.64e+01    6.75e+01      0   Skip BFGS
   4  1.00e+00  2.10e+01  4.77e+01  6.88e+01    1.10e+02      0
   5  1.00e+00  1.37e+01  4.59e+01  5.97e+01    1.87e+01      0
   6  2.51e-01  1.05e+01  4.69e+01  2.22e+01    3.38e+01      0   Skip BFGS
   7  2.51e-01  5.76e+00  5.70e+01  2.00e+01    5.87e+01      0
------------------------- STOP! -------------------------
1 : |fc-fOld| = 0.0000e+00 <= tolF*(1+|f0|) = 7.9079e+01
1 : |xc-x_last| = 3.7089e-01 <= tolX*(1+|x0|) = 2.4026e+00
0 : |proj(x-g)-x|    = 5.8690e+01 <= tolG          = 1.0000e-01
0 : |proj(x-g)-x|    = 5.8690e+01 <= 1e3*eps       = 1.0000e-02
0 : maxIter   =      20    <= iter          =      8
------------------------- DONE! -------------------------
min diffusion distance  114.18394340269786 max diffusion distance  510.64611877484225
SimPEG.InvProblem will set Regularization.reference_model to m0.
SimPEG.InvProblem will set Regularization.reference_model to m0.
SimPEG.InvProblem will set Regularization.reference_model to m0.

                        SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.
                        ***Done using same Solver, and solver_opts as the Simulation3DMagneticFluxDensity 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  3.29e+00  1.64e+03  0.00e+00  1.64e+03    7.01e+02      0
   1  3.29e+00  1.15e+03  8.05e+01  1.41e+03    2.20e+03      0
   2  3.29e+00  8.36e+01  6.15e+01  2.86e+02    3.00e+02      0
   3  8.21e-01  6.64e+01  5.34e+01  1.10e+02    2.44e+02      0   Skip BFGS
   4  8.21e-01  3.94e+01  5.68e+01  8.60e+01    2.04e+02      1
   5  8.21e-01  1.69e+01  5.73e+01  6.40e+01    1.99e+02      0   Skip BFGS
   6  2.05e-01  5.36e+00  5.41e+01  1.65e+01    1.43e+01      0
------------------------- STOP! -------------------------
1 : |fc-fOld| = 0.0000e+00 <= tolF*(1+|f0|) = 1.6433e+02
1 : |xc-x_last| = 5.6398e-01 <= tolX*(1+|x0|) = 2.4026e+00
0 : |proj(x-g)-x|    = 1.4283e+01 <= tolG          = 1.0000e-01
0 : |proj(x-g)-x|    = 1.4283e+01 <= 1e3*eps       = 1.0000e-02
0 : maxIter   =      20    <= iter          =      7
------------------------- DONE! -------------------------

import discretize
from SimPEG import (
    maps,
    utils,
    data_misfit,
    regularization,
    optimization,
    inversion,
    inverse_problem,
    directives,
)
import numpy as np
from SimPEG.electromagnetics import frequency_domain as FDEM, time_domain as TDEM, mu_0
import matplotlib.pyplot as plt
import matplotlib

try:
    from pymatsolver import Pardiso as Solver
except ImportError:
    from SimPEG import SolverLU as Solver


def run(plotIt=True, saveFig=False):
    # Set up cylindrically symmeric mesh
    cs, ncx, ncz, npad = 10.0, 15, 25, 13  # padded cyl mesh
    hx = [(cs, ncx), (cs, npad, 1.3)]
    hz = [(cs, npad, -1.3), (cs, ncz), (cs, npad, 1.3)]
    mesh = discretize.CylindricalMesh([hx, 1, hz], "00C")

    # Conductivity model
    layerz = np.r_[-200.0, -100.0]
    layer = (mesh.cell_centers_z >= layerz[0]) & (mesh.cell_centers_z <= layerz[1])
    active = mesh.cell_centers_z < 0.0
    sig_half = 1e-2  # Half-space conductivity
    sig_air = 1e-8  # Air conductivity
    sig_layer = 5e-2  # Layer conductivity
    sigma = np.ones(mesh.shape_cells[2]) * sig_air
    sigma[active] = sig_half
    sigma[layer] = sig_layer

    # Mapping
    actMap = maps.InjectActiveCells(mesh, active, np.log(1e-8), nC=mesh.shape_cells[2])
    mapping = maps.ExpMap(mesh) * maps.SurjectVertical1D(mesh) * actMap
    mtrue = np.log(sigma[active])

    # ----- FDEM problem & survey ----- #
    rxlocs = utils.ndgrid([np.r_[50.0], np.r_[0], np.r_[0.0]])
    bzr = FDEM.Rx.PointMagneticFluxDensitySecondary(rxlocs, "z", "real")
    bzi = FDEM.Rx.PointMagneticFluxDensitySecondary(rxlocs, "z", "imag")

    freqs = np.logspace(2, 3, 5)
    srcLoc = np.array([0.0, 0.0, 0.0])

    print(
        "min skin depth = ",
        500.0 / np.sqrt(freqs.max() * sig_half),
        "max skin depth = ",
        500.0 / np.sqrt(freqs.min() * sig_half),
    )
    print(
        "max x ",
        mesh.cell_centers_x.max(),
        "min z ",
        mesh.cell_centers_z.min(),
        "max z ",
        mesh.cell_centers_z.max(),
    )

    source_list = [
        FDEM.Src.MagDipole([bzr, bzi], freq, srcLoc, orientation="Z") for freq in freqs
    ]

    surveyFD = FDEM.Survey(source_list)
    prbFD = FDEM.Simulation3DMagneticFluxDensity(
        mesh, survey=surveyFD, sigmaMap=mapping, solver=Solver
    )
    rel_err = 0.03
    dataFD = prbFD.make_synthetic_data(mtrue, relative_error=rel_err, add_noise=True)
    dataFD.noise_floor = np.linalg.norm(dataFD.dclean) * 1e-5

    # FDEM inversion
    np.random.seed(1)
    dmisfit = data_misfit.L2DataMisfit(simulation=prbFD, data=dataFD)
    regMesh = discretize.TensorMesh([mesh.h[2][mapping.maps[-1].indActive]])
    reg = regularization.WeightedLeastSquares(regMesh)
    opt = optimization.InexactGaussNewton(maxIterCG=10)
    invProb = inverse_problem.BaseInvProblem(dmisfit, reg, opt)

    # Inversion Directives
    beta = directives.BetaSchedule(coolingFactor=4, coolingRate=3)
    betaest = directives.BetaEstimate_ByEig(beta0_ratio=1.0, seed=518936)
    target = directives.TargetMisfit()
    directiveList = [beta, betaest, target]

    inv = inversion.BaseInversion(invProb, directiveList=directiveList)
    m0 = np.log(np.ones(mtrue.size) * sig_half)
    reg.alpha_s = 5e-1
    reg.alpha_x = 1.0
    prbFD.counter = opt.counter = utils.Counter()
    opt.remember("xc")
    moptFD = inv.run(m0)

    # TDEM problem
    times = np.logspace(-4, np.log10(2e-3), 10)
    print(
        "min diffusion distance ",
        1.28 * np.sqrt(times.min() / (sig_half * mu_0)),
        "max diffusion distance ",
        1.28 * np.sqrt(times.max() / (sig_half * mu_0)),
    )
    rx = TDEM.Rx.PointMagneticFluxDensity(rxlocs, times, "z")
    src = TDEM.Src.MagDipole(
        [rx],
        waveform=TDEM.Src.StepOffWaveform(),
        location=srcLoc,  # same src location as FDEM problem
    )

    surveyTD = TDEM.Survey([src])
    prbTD = TDEM.Simulation3DMagneticFluxDensity(
        mesh, survey=surveyTD, sigmaMap=mapping, solver=Solver
    )
    prbTD.time_steps = [(5e-5, 10), (1e-4, 10), (5e-4, 10)]

    rel_err = 0.03
    dataTD = prbTD.make_synthetic_data(mtrue, relative_error=rel_err, add_noise=True)
    dataTD.noise_floor = np.linalg.norm(dataTD.dclean) * 1e-5

    # TDEM inversion
    dmisfit = data_misfit.L2DataMisfit(simulation=prbTD, data=dataTD)
    regMesh = discretize.TensorMesh([mesh.h[2][mapping.maps[-1].indActive]])
    reg = regularization.WeightedLeastSquares(regMesh)
    opt = optimization.InexactGaussNewton(maxIterCG=10)
    invProb = inverse_problem.BaseInvProblem(dmisfit, reg, opt)

    # directives
    beta = directives.BetaSchedule(coolingFactor=4, coolingRate=3)
    betaest = directives.BetaEstimate_ByEig(beta0_ratio=1.0, seed=518936)
    target = directives.TargetMisfit()
    directiveList = [beta, betaest, target]

    inv = inversion.BaseInversion(invProb, directiveList=directiveList)
    m0 = np.log(np.ones(mtrue.size) * sig_half)
    reg.alpha_s = 5e-1
    reg.alpha_x = 1.0
    prbTD.counter = opt.counter = utils.Counter()
    opt.remember("xc")
    moptTD = inv.run(m0)

    # Plot the results
    if plotIt:
        plt.figure(figsize=(10, 8))
        ax0 = plt.subplot2grid((2, 2), (0, 0), rowspan=2)
        ax1 = plt.subplot2grid((2, 2), (0, 1))
        ax2 = plt.subplot2grid((2, 2), (1, 1))

        fs = 13  # fontsize
        matplotlib.rcParams["font.size"] = fs

        # Plot the model
        # z_true = np.repeat(mesh.cell_centers_z[active][1:], 2, axis=0)
        # z_true = np.r_[mesh.cell_centers_z[active][0], z_true, mesh.cell_centers_z[active][-1]]
        activeN = mesh.nodes_z <= 0.0 + cs / 2.0
        z_true = np.repeat(mesh.nodes_z[activeN][1:-1], 2, axis=0)
        z_true = np.r_[mesh.nodes_z[activeN][0], z_true, mesh.nodes_z[activeN][-1]]
        sigma_true = np.repeat(sigma[active], 2, axis=0)

        ax0.semilogx(sigma_true, z_true, "k-", lw=2, label="True")

        ax0.semilogx(
            np.exp(moptFD),
            mesh.cell_centers_z[active],
            "bo",
            ms=6,
            markeredgecolor="k",
            markeredgewidth=0.5,
            label="FDEM",
        )
        ax0.semilogx(
            np.exp(moptTD),
            mesh.cell_centers_z[active],
            "r*",
            ms=10,
            markeredgecolor="k",
            markeredgewidth=0.5,
            label="TDEM",
        )
        ax0.set_ylim(-700, 0)
        ax0.set_xlim(5e-3, 1e-1)

        ax0.set_xlabel("Conductivity (S/m)", fontsize=fs)
        ax0.set_ylabel("Depth (m)", fontsize=fs)
        ax0.grid(which="both", color="k", alpha=0.5, linestyle="-", linewidth=0.2)
        ax0.legend(fontsize=fs, loc=4)

        # plot the data misfits - negative b/c we choose positive to be in the
        # direction of primary

        ax1.plot(freqs, -dataFD.dobs[::2], "k-", lw=2, label="Obs (real)")
        ax1.plot(freqs, -dataFD.dobs[1::2], "k--", lw=2, label="Obs (imag)")

        dpredFD = prbFD.dpred(moptTD)
        ax1.loglog(
            freqs,
            -dpredFD[::2],
            "bo",
            ms=6,
            markeredgecolor="k",
            markeredgewidth=0.5,
            label="Pred (real)",
        )
        ax1.loglog(
            freqs, -dpredFD[1::2], "b+", ms=10, markeredgewidth=2.0, label="Pred (imag)"
        )

        ax2.loglog(times, dataTD.dobs, "k-", lw=2, label="Obs")
        ax2.loglog(
            times,
            prbTD.dpred(moptTD),
            "r*",
            ms=10,
            markeredgecolor="k",
            markeredgewidth=0.5,
            label="Pred",
        )
        ax2.set_xlim(times.min() - 1e-5, times.max() + 1e-4)

        # Labels, gridlines, etc
        ax2.grid(which="both", alpha=0.5, linestyle="-", linewidth=0.2)
        ax1.grid(which="both", alpha=0.5, linestyle="-", linewidth=0.2)

        ax1.set_xlabel("Frequency (Hz)", fontsize=fs)
        ax1.set_ylabel("Vertical magnetic field (-T)", fontsize=fs)

        ax2.set_xlabel("Time (s)", fontsize=fs)
        ax2.set_ylabel("Vertical magnetic field (T)", fontsize=fs)

        ax2.legend(fontsize=fs, loc=3)
        ax1.legend(fontsize=fs, loc=3)
        ax1.set_xlim(freqs.max() + 1e2, freqs.min() - 1e1)

        ax0.set_title("(a) Recovered Models", fontsize=fs)
        ax1.set_title("(b) FDEM observed vs. predicted", fontsize=fs)
        ax2.set_title("(c) TDEM observed vs. predicted", fontsize=fs)

        plt.tight_layout(pad=1.5)

        if saveFig is True:
            plt.savefig("example1.png", dpi=600)


if __name__ == "__main__":
    run(plotIt=True, saveFig=True)
    plt.show()

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

Estimated memory usage: 119 MB

Gallery generated by Sphinx-Gallery