MT: 1D: InversionΒΆ

Forward model 1D MT data. Setup and run a MT 1D inversion.

../../../_images/sphx_glr_plot_inversion_MT1D_001.png

Out:

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  2.41e+02  4.93e+04  0.00e+00  4.93e+04    8.88e+03      0
   1  2.41e+02  6.71e+03  2.97e+01  1.39e+04    2.00e+03      0
   2  2.41e+02  4.63e+03  2.84e+01  1.15e+04    4.44e+02      0   Skip BFGS
   3  2.41e+02  4.43e+03  2.67e+01  1.09e+04    2.92e+02      0   Skip BFGS
   4  6.03e+01  4.15e+03  2.48e+01  5.64e+03    1.25e+03      0
   5  6.03e+01  2.13e+03  4.30e+01  4.73e+03    2.17e+02      0
   6  6.03e+01  2.05e+03  4.29e+01  4.64e+03    1.54e+02      0
   7  6.03e+01  1.83e+03  4.56e+01  4.58e+03    2.06e+02      0   Skip BFGS
   8  1.51e+01  1.65e+03  4.76e+01  2.37e+03    4.60e+02      0   Skip BFGS
   9  1.51e+01  8.11e+02  8.46e+01  2.09e+03    9.76e+02      0
  10  1.51e+01  2.95e+02  8.14e+01  1.52e+03    1.74e+02      0
  11  1.51e+01  1.62e+02  8.38e+01  1.43e+03    3.03e+02      0   Skip BFGS
  12  3.77e+00  1.55e+02  8.32e+01  4.68e+02    1.97e+02      1
  13  3.77e+00  9.13e+01  9.12e+01  4.35e+02    2.18e+02      0
  14  3.77e+00  8.52e+01  9.05e+01  4.26e+02    2.01e+02      1
  15  3.77e+00  5.45e+01  9.15e+01  3.99e+02    1.55e+02      0   Skip BFGS
------------------------- STOP! -------------------------
1 : |fc-fOld| = 0.0000e+00 <= tolF*(1+|f0|) = 4.9305e+03
1 : |xc-x_last| = 2.3205e-01 <= tolX*(1+|x0|) = 5.6454e+00
0 : |proj(x-g)-x|    = 1.5454e+02 <= tolG          = 1.0000e-01
0 : |proj(x-g)-x|    = 1.5454e+02 <= 1e3*eps       = 1.0000e-02
0 : maxIter   =      25    <= iter          =     16
------------------------- DONE! -------------------------

import matplotlib.pyplot as plt

import numpy as np
import SimPEG as simpeg
from SimPEG.EM import NSEM

np.random.seed(1983)


def run(plotIt=True):

    # Setup the forward modeling
    # Setting up 1D mesh and conductivity models to forward model data.
    # Frequency
    nFreq = 26
    freqs = np.logspace(2, -3, nFreq)
    # Set mesh parameters
    ct = 10
    air = simpeg.Utils.meshTensor([(ct, 25, 1.4)])
    core = np.concatenate((
        np.kron(simpeg.Utils.meshTensor([(ct, 15, -1.2)]), np.ones((8, ))),
        simpeg.Utils.meshTensor([(ct, 5)])))
    bot = simpeg.Utils.meshTensor([(core[0], 20, -1.4)])
    x0 = -np.array([np.sum(np.concatenate((core, bot)))])
    # Make the model
    m1d = simpeg.Mesh.TensorMesh([np.concatenate((bot, core, air))], x0=x0)

    # Setup model variables
    active = m1d.vectorCCx < 0.
    layer1 = (m1d.vectorCCx < -500.) & (m1d.vectorCCx >= -800.)
    layer2 = (m1d.vectorCCx < -3500.) & (m1d.vectorCCx >= -5000.)
    # Set the conductivity values
    sig_half = 1e-2
    sig_air = 1e-8
    sig_layer1 = .2
    sig_layer2 = .2
    # Make the true model
    sigma_true = np.ones(m1d.nCx) * sig_air
    sigma_true[active] = sig_half
    sigma_true[layer1] = sig_layer1
    sigma_true[layer2] = sig_layer2
    # Extract the model
    m_true = np.log(sigma_true[active])
    # Make the background model
    sigma_0 = np.ones(m1d.nCx) * sig_air
    sigma_0[active] = sig_half
    m_0 = np.log(sigma_0[active])

    # Set the mapping
    actMap = simpeg.Maps.InjectActiveCells(
        m1d, active, np.log(1e-8), nC=m1d.nCx)
    mappingExpAct = simpeg.Maps.ExpMap(m1d) * actMap

    # Setup the layout of the survey, set the sources and the connected receivers
    # Receivers
    rxList = []
    rxList.append(
        NSEM.Rx.Point_impedance1D(simpeg.mkvc(np.array([-0.5]), 2).T, 'real'))
    rxList.append(
        NSEM.Rx.Point_impedance1D(simpeg.mkvc(np.array([-0.5]), 2).T, 'imag'))
    # Source list
    srcList = []
    for freq in freqs:
            srcList.append(NSEM.Src.Planewave_xy_1Dprimary(rxList, freq))
    # Make the survey
    survey = NSEM.Survey(srcList)
    survey.mtrue = m_true

    # Set the problem
    problem = NSEM.Problem1D_ePrimSec(
        m1d, sigmaPrimary=sigma_0, sigmaMap=mappingExpAct)
    problem.pair(survey)

    # Forward model data
    # Project the data
    survey.dtrue = survey.dpred(m_true)
    survey.dobs = (
        survey.dtrue + 0.01 *
        abs(survey.dtrue) *
        np.random.randn(*survey.dtrue.shape))

    # Assign uncertainties
    std = 0.025  # 5% std
    survey.std = np.abs(survey.dobs * std)
    # Assign the data weight
    Wd = 1. / survey.std

    # Setup the inversion proceedure
    # Define a counter
    C = simpeg.Utils.Counter()
    # Optimization
    opt = simpeg.Optimization.InexactGaussNewton(maxIter=25)
    opt.counter = C
    opt.LSshorten = 0.1
    opt.remember('xc')
    # Data misfit
    dmis = simpeg.DataMisfit.l2_DataMisfit(survey)
    dmis.W = Wd
    # Regularization - with a regularization mesh
    regMesh = simpeg.Mesh.TensorMesh([m1d.hx[active]], m1d.x0)
    reg = simpeg.Regularization.Tikhonov(regMesh)
    reg.mrefInSmooth = True
    reg.alpha_s = 1e-2
    reg.alpha_x = 1.
    reg.mrefInSmooth = True
    # Inversion problem
    invProb = simpeg.InvProblem.BaseInvProblem(dmis, reg, opt)
    invProb.counter = C
    # Beta schedule
    beta = simpeg.Directives.BetaSchedule()
    beta.coolingRate = 4.
    beta.coolingFactor = 4.
    # Initial estimate of beta
    betaest = simpeg.Directives.BetaEstimate_ByEig(beta0_ratio=10.)
    # Target misfit stop
    targmis = simpeg.Directives.TargetMisfit()
    targmis.target = survey.nD
    # Create an inversion object
    directives = [beta, betaest, targmis]
    inv = simpeg.Inversion.BaseInversion(invProb, directiveList=directives)

    # Run the inversion
    mopt = inv.run(m_0)

    if plotIt:
        fig = NSEM.Utils.dataUtils.plotMT1DModelData(problem, [mopt])
        fig.suptitle('Target - smooth true')
        fig.axes[0].set_ylim([-10000, 500])

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

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

Gallery generated by Sphinx-Gallery