PF: Magnetic: Inversion LinearΒΆ

Create a synthetic block model and invert with a compact norm

  • ../../../_images/sphx_glr_plot_inversion_linear_0012.png
  • ../../../_images/sphx_glr_plot_inversion_linear_0021.png
  • ../../../_images/sphx_glr_plot_inversion_linear_0031.png

Out:

Begin calculation of forward operator: ind
Done 0.0 %
Done 10.0 %
Done 20.0 %
Done 30.0 %
Done 40.0 %
Done 50.0 %
Done 60.0 %
Done 70.0 %
Done 80.0 %
Done 90.0 %
Done 100% ...forward operator completed!!

SimPEG.DataMisfit.l2_DataMisfit assigning default eps of 1e-5 * ||dobs||

    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
=============================== Projected GNCG ===============================
  #     beta     phi_d     phi_m       f      |proj(x-g)-x|  LS    Comment
-----------------------------------------------------------------------------
x0 has any nan: 0
   0  1.05e+07  2.87e+05  1.71e-06  2.87e+05    6.27e+01      0
   1  2.73e+06  1.05e+04  8.39e-04  1.28e+04    5.59e+01      0
   2  9.53e+05  1.01e+03  1.58e-03  2.51e+03    3.77e+01      0   Skip BFGS
   3  5.36e+05  3.20e+02  1.70e-03  1.23e+03    3.76e+01      0   Skip BFGS
   4  3.48e+05  2.50e+02  1.69e-03  8.38e+02    6.01e+01      1   Skip BFGS
Reached starting chifact with l2-norm regularization: Start IRLS steps...
eps_p: 0.010564457919211121 eps_q: 0.010564457919211121
Eps_p: 0.008803714932675935
Eps_q: 0.008803714932675935
delta phim: 1.722e+02
   5  3.48e+05  1.81e+02  2.31e-03  9.85e+02    3.63e+01      0
Eps_p: 0.00733642911056328
Eps_q: 0.00733642911056328
delta phim: 8.227e+03
   6  3.48e+05  1.81e+02  2.43e-03  1.03e+03    3.63e+01      7   Skip BFGS
Eps_p: 0.0061136909254694005
Eps_q: 0.0061136909254694005
delta phim: 3.137e-01
   7  3.48e+05  2.11e+02  2.42e-03  1.05e+03    5.99e+01      1
Eps_p: 0.005094742437891167
Eps_q: 0.005094742437891167
delta phim: 1.229e-01
   8  3.48e+05  1.85e+02  2.53e-03  1.07e+03    3.26e+01      0
Eps_p: 0.0042456186982426395
Eps_q: 0.0042456186982426395
delta phim: 7.787e-02
   9  3.48e+05  1.85e+02  2.51e-03  1.06e+03    3.32e+01      6   Skip BFGS
Eps_p: 0.0035380155818688663
Eps_q: 0.0035380155818688663
delta phim: 2.482e-01
  10  3.48e+05  1.96e+02  2.44e-03  1.04e+03    5.99e+01      1
Eps_p: 0.0029483463182240553
Eps_q: 0.0029483463182240553
delta phim: 1.326e-01
  11  3.48e+05  1.88e+02  2.34e-03  1.00e+03    5.58e+01      0
Eps_p: 0.0024569552651867127
Eps_q: 0.0024569552651867127
delta phim: 6.595e-02
  12  3.48e+05  1.89e+02  2.23e-03  9.64e+02    2.85e+01      3   Skip BFGS
Eps_p: 0.0020474627209889273
Eps_q: 0.0020474627209889273
delta phim: 2.098e-01
  13  3.48e+05  1.87e+02  2.12e-03  9.24e+02    5.70e+01      0
Eps_p: 0.0017062189341574396
Eps_q: 0.0017062189341574396
delta phim: 7.250e-02
  14  3.48e+05  1.90e+02  1.98e-03  8.79e+02    6.01e+01      2   Skip BFGS
Eps_p: 0.0014218491117978663
Eps_q: 0.0014218491117978663
delta phim: 1.746e-01
  15  3.48e+05  1.88e+02  1.83e-03  8.26e+02    5.18e+01      0
Eps_p: 0.0011848742598315554
Eps_q: 0.0011848742598315554
delta phim: 1.117e-01
  16  3.48e+05  1.89e+02  1.71e-03  7.84e+02    2.95e+01      0   Skip BFGS
Eps_p: 0.0009873952165262963
Eps_q: 0.0009873952165262963
delta phim: 1.373e-01
  17  3.48e+05  1.90e+02  1.61e-03  7.51e+02    4.17e+01      0
Eps_p: 0.000822829347105247
Eps_q: 0.000822829347105247
delta phim: 8.836e-02
  18  3.48e+05  1.95e+02  1.49e-03  7.13e+02    6.15e+01      2   Skip BFGS
Eps_p: 0.0006856911225877058
Eps_q: 0.0006856911225877058
delta phim: 1.491e-01
  19  3.48e+05  1.90e+02  1.39e-03  6.76e+02    1.03e+01      0
Eps_p: 0.0005714092688230882
Eps_q: 0.0005714092688230882
delta phim: 4.737e-02
  20  3.48e+05  1.90e+02  1.35e-03  6.58e+02    9.42e+00      0
Eps_p: 0.00047617439068590686
Eps_q: 0.00047617439068590686
delta phim: 1.355e-02
  21  3.48e+05  1.88e+02  1.29e-03  6.38e+02    1.67e+01      0
Eps_p: 0.0003968119922382557
Eps_q: 0.0003968119922382557
delta phim: 1.511e-01
  22  3.48e+05  1.89e+02  1.22e-03  6.13e+02    2.47e+01      0
Eps_p: 0.00033067666019854646
Eps_q: 0.00033067666019854646
delta phim: 2.169e-01
  23  3.48e+05  1.90e+02  1.14e-03  5.87e+02    2.22e+01      0
Eps_p: 0.0002755638834987887
Eps_q: 0.0002755638834987887
delta phim: 2.334e-01
  24  3.48e+05  1.93e+02  1.05e-03  5.59e+02    2.99e+01      0
Eps_p: 0.00022963656958232393
Eps_q: 0.00022963656958232393
delta phim: 1.935e-01
  25  3.48e+05  1.96e+02  9.50e-04  5.27e+02    3.20e+01      0
Eps_p: 0.00019136380798526995
Eps_q: 0.00019136380798526995
delta phim: 1.540e-01
  26  3.48e+05  2.01e+02  8.44e-04  4.95e+02    3.06e+01      0
Eps_p: 0.00015946983998772497
Eps_q: 0.00015946983998772497
delta phim: 1.170e-01
  27  3.48e+05  2.05e+02  7.49e-04  4.65e+02    9.93e+00      0   Skip BFGS
Eps_p: 0.00013289153332310416
Eps_q: 0.00013289153332310416
delta phim: 3.135e-02
  28  3.48e+05  2.07e+02  6.67e-04  4.39e+02    9.59e+00      0   Skip BFGS
Eps_p: 0.00011074294443592014
Eps_q: 0.00011074294443592014
delta phim: 2.117e-03
  29  3.48e+05  2.09e+02  5.97e-04  4.16e+02    1.04e+01      0   Skip BFGS
Eps_p: 9.228578702993345e-05
Eps_q: 9.228578702993345e-05
delta phim: 2.288e-02
  30  3.48e+05  2.10e+02  5.37e-04  3.97e+02    1.03e+01      0   Skip BFGS
Eps_p: 7.690482252494455e-05
Eps_q: 7.690482252494455e-05
delta phim: 3.368e-02
  31  3.48e+05  2.10e+02  4.85e-04  3.79e+02    1.01e+01      0   Skip BFGS
Eps_p: 6.408735210412046e-05
Eps_q: 6.408735210412046e-05
delta phim: 2.727e-02
  32  3.48e+05  2.10e+02  4.39e-04  3.63e+02    9.82e+00      0   Skip BFGS
Eps_p: 5.340612675343372e-05
Eps_q: 5.340612675343372e-05
delta phim: 2.375e-02
  33  3.48e+05  2.09e+02  3.97e-04  3.48e+02    9.05e+00      0
Eps_p: 4.450510562786143e-05
Eps_q: 4.450510562786143e-05
delta phim: 1.425e-02
  34  3.48e+05  2.09e+02  3.60e-04  3.34e+02    8.78e+00      0
Eps_p: 3.7087588023217865e-05
Eps_q: 3.7087588023217865e-05
delta phim: 4.267e-03
  35  3.48e+05  2.08e+02  3.27e-04  3.22e+02    8.62e+00      0
Eps_p: 3.090632335268156e-05
Eps_q: 3.090632335268156e-05
delta phim: 6.279e-03
  36  3.48e+05  2.06e+02  2.98e-04  3.10e+02    9.00e+00      0
Eps_p: 2.5755269460567964e-05
Eps_q: 2.5755269460567964e-05
delta phim: 1.480e-02
  37  3.48e+05  2.04e+02  2.71e-04  2.99e+02    9.88e+00      0
Eps_p: 2.1462724550473305e-05
Eps_q: 2.1462724550473305e-05
delta phim: 2.612e-02
  38  3.48e+05  2.02e+02  2.48e-04  2.88e+02    8.89e+00      0
Eps_p: 1.7885603792061087e-05
Eps_q: 1.7885603792061087e-05
delta phim: 3.307e-02
  39  3.48e+05  1.98e+02  2.26e-04  2.77e+02    9.84e+00      0
Eps_p: 1.4904669826717574e-05
Eps_q: 1.4904669826717574e-05
delta phim: 3.476e-02
  40  3.48e+05  1.94e+02  2.07e-04  2.66e+02    1.03e+01      0
Eps_p: 1.2420558188931312e-05
Eps_q: 1.2420558188931312e-05
delta phim: 4.981e-02
  41  3.48e+05  1.89e+02  1.90e-04  2.55e+02    9.30e+00      0   Skip BFGS
Eps_p: 1.035046515744276e-05
Eps_q: 1.035046515744276e-05
delta phim: 3.831e-02
  42  3.48e+05  1.85e+02  1.83e-04  2.49e+02    1.06e+01      0   Skip BFGS
Eps_p: 8.625387631202301e-06
Eps_q: 8.625387631202301e-06
delta phim: 4.096e-02
  43  3.48e+05  1.83e+02  1.71e-04  2.43e+02    1.03e+01      0   Skip BFGS
Eps_p: 7.187823026001918e-06
Eps_q: 7.187823026001918e-06
delta phim: 1.287e-02
  44  3.48e+05  1.82e+02  1.58e-04  2.37e+02    1.37e+01      0
Reach maximum number of IRLS cycles: 40
------------------------- STOP! -------------------------
1 : |fc-fOld| = 0.0000e+00 <= tolF*(1+|f0|) = 2.8716e+04
1 : |xc-x_last| = 3.2683e-03 <= tolX*(1+|x0|) = 1.0075e-01
0 : |proj(x-g)-x|    = 1.3720e+01 <= tolG          = 1.0000e-01
0 : |proj(x-g)-x|    = 1.3720e+01 <= 1e3*eps       = 1.0000e-02
0 : maxIter   =     100    <= iter          =     45
------------------------- DONE! -------------------------

import matplotlib.pyplot as plt
import numpy as np

from SimPEG import Mesh
from SimPEG import Utils
from SimPEG import Maps
from SimPEG import Regularization
from SimPEG import DataMisfit
from SimPEG import Optimization
from SimPEG import InvProblem
from SimPEG import Directives
from SimPEG import Inversion
from SimPEG import PF


def run(plotIt=True):

    # Define the inducing field parameter
    H0 = (50000, 90, 0)

    # Create a mesh
    dx = 5.

    hxind = [(dx, 5, -1.3), (dx, 10), (dx, 5, 1.3)]
    hyind = [(dx, 5, -1.3), (dx, 10), (dx, 5, 1.3)]
    hzind = [(dx, 5, -1.3), (dx, 10)]

    mesh = Mesh.TensorMesh([hxind, hyind, hzind], 'CCC')

    # Get index of the center
    midx = int(mesh.nCx/2)
    midy = int(mesh.nCy/2)

    # Lets create a simple Gaussian topo and set the active cells
    [xx, yy] = np.meshgrid(mesh.vectorNx, mesh.vectorNy)
    zz = -np.exp((xx**2 + yy**2) / 75**2) + mesh.vectorNz[-1]

    # We would usually load a topofile
    topo = np.c_[Utils.mkvc(xx), Utils.mkvc(yy), Utils.mkvc(zz)]

    # Go from topo to actv cells
    actv = Utils.surface2ind_topo(mesh, topo, 'N')
    actv = np.asarray([inds for inds, elem in enumerate(actv, 1) if elem],
                      dtype=int) - 1

    # Create active map to go from reduce space to full
    actvMap = Maps.InjectActiveCells(mesh, actv, -100)
    nC = len(actv)

    # Create and array of observation points
    xr = np.linspace(-20., 20., 20)
    yr = np.linspace(-20., 20., 20)
    X, Y = np.meshgrid(xr, yr)

    # Move the observation points 5m above the topo
    Z = -np.exp((X**2 + Y**2) / 75**2) + mesh.vectorNz[-1] + 5.

    # Create a MAGsurvey
    rxLoc = np.c_[Utils.mkvc(X.T), Utils.mkvc(Y.T), Utils.mkvc(Z.T)]
    rxLoc = PF.BaseMag.RxObs(rxLoc)
    srcField = PF.BaseMag.SrcField([rxLoc], param=H0)
    survey = PF.BaseMag.LinearSurvey(srcField)

    # We can now create a susceptibility model and generate data
    # Here a simple block in half-space
    model = np.zeros((mesh.nCx, mesh.nCy, mesh.nCz))
    model[(midx-2):(midx+2), (midy-2):(midy+2), -6:-2] = 0.02
    model = Utils.mkvc(model)
    model = model[actv]

    # Create active map to go from reduce set to full
    actvMap = Maps.InjectActiveCells(mesh, actv, -100)

    # Creat reduced identity map
    idenMap = Maps.IdentityMap(nP=nC)

    # Create the forward model operator
    prob = PF.Magnetics.MagneticIntegral(mesh, chiMap=idenMap, actInd=actv)

    # Pair the survey and problem
    survey.pair(prob)

    # Compute linear forward operator and compute some data
    d = prob.fields(model)

    # Add noise and uncertainties
    # We add some random Gaussian noise (1nT)
    data = d + np.random.randn(len(d))
    wd = np.ones(len(data))*1.  # Assign flat uncertainties

    survey.dobs = data
    survey.std = wd
    survey.mtrue = model

    # Create sensitivity weights from our linear forward operator
    rxLoc = survey.srcField.rxList[0].locs
    wr = np.sum(prob.G**2., axis=0)**0.5
    wr = (wr/np.max(wr))

    # Create a regularization
    reg = Regularization.Sparse(mesh, indActive=actv, mapping=idenMap)
    reg.cell_weights = wr
    reg.mref = np.zeros(nC)
    reg.norms = np.c_[0, 1, 1, 1]
    # reg.eps_p, reg.eps_q = 1e-0, 1e-0

    # Data misfit function
    dmis = DataMisfit.l2_DataMisfit(survey)
    dmis.W = 1/wd

    # Add directives to the inversion
    opt = Optimization.ProjectedGNCG(maxIter=100, lower=0., upper=1.,
                                     maxIterLS=20, maxIterCG=10, tolCG=1e-3)
    invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
    betaest = Directives.BetaEstimate_ByEig(beta0_ratio=1e-1)

    # Here is where the norms are applied
    # Use pick a treshold parameter empirically based on the distribution of
    #  model parameters
    IRLS = Directives.Update_IRLS(
        f_min_change=1e-4, maxIRLSiter=40
    )
    saveDict = Directives.SaveOutputEveryIteration(save_txt=False)
    update_Jacobi = Directives.UpdatePreconditioner()
    inv = Inversion.BaseInversion(
        invProb, directiveList=[IRLS, betaest, update_Jacobi, saveDict]
    )

    # Run the inversion
    m0 = np.ones(nC)*1e-4  # Starting model
    mrec = inv.run(m0)

    if plotIt:
        # Here is the recovered susceptibility model
        ypanel = midx
        zpanel = -5
        m_l2 = actvMap * invProb.l2model
        m_l2[m_l2 == -100] = np.nan

        m_lp = actvMap * mrec
        m_lp[m_lp == -100] = np.nan

        m_true = actvMap * model
        m_true[m_true == -100] = np.nan

        # Plot the data
        PF.Magnetics.plot_obs_2D(rxLoc, d=d)

        plt.figure()

        # Plot L2 model
        ax = plt.subplot(321)
        mesh.plotSlice(m_l2, ax=ax, normal='Z', ind=zpanel,
                       grid=True, clim=(model.min(), model.max()))
        plt.plot(([mesh.vectorCCx[0], mesh.vectorCCx[-1]]),
                 ([mesh.vectorCCy[ypanel], mesh.vectorCCy[ypanel]]), color='w')
        plt.title('Plan l2-model.')
        plt.gca().set_aspect('equal')
        plt.ylabel('y')
        ax.xaxis.set_visible(False)
        plt.gca().set_aspect('equal', adjustable='box')

        # Vertica section
        ax = plt.subplot(322)
        mesh.plotSlice(m_l2, ax=ax, normal='Y', ind=midx,
                       grid=True, clim=(model.min(), model.max()))
        plt.plot(([mesh.vectorCCx[0], mesh.vectorCCx[-1]]),
                 ([mesh.vectorCCz[zpanel], mesh.vectorCCz[zpanel]]), color='w')
        plt.title('E-W l2-model.')
        plt.gca().set_aspect('equal')
        ax.xaxis.set_visible(False)
        plt.ylabel('z')
        plt.gca().set_aspect('equal', adjustable='box')

        # Plot Lp model
        ax = plt.subplot(323)
        mesh.plotSlice(m_lp, ax=ax, normal='Z', ind=zpanel,
                       grid=True, clim=(model.min(), model.max()))
        plt.plot(([mesh.vectorCCx[0], mesh.vectorCCx[-1]]),
                 ([mesh.vectorCCy[ypanel], mesh.vectorCCy[ypanel]]), color='w')
        plt.title('Plan lp-model.')
        plt.gca().set_aspect('equal')
        ax.xaxis.set_visible(False)
        plt.ylabel('y')
        plt.gca().set_aspect('equal', adjustable='box')

        # Vertical section
        ax = plt.subplot(324)
        mesh.plotSlice(m_lp, ax=ax, normal='Y', ind=midx,
                       grid=True, clim=(model.min(), model.max()))
        plt.plot(([mesh.vectorCCx[0], mesh.vectorCCx[-1]]),
                 ([mesh.vectorCCz[zpanel], mesh.vectorCCz[zpanel]]), color='w')
        plt.title('E-W lp-model.')
        plt.gca().set_aspect('equal')
        ax.xaxis.set_visible(False)
        plt.ylabel('z')
        plt.gca().set_aspect('equal', adjustable='box')

        # Plot True model
        ax = plt.subplot(325)
        mesh.plotSlice(m_true, ax=ax, normal='Z', ind=zpanel,
                       grid=True, clim=(model.min(), model.max()))
        plt.plot(([mesh.vectorCCx[0], mesh.vectorCCx[-1]]),
                 ([mesh.vectorCCy[ypanel], mesh.vectorCCy[ypanel]]), color='w')
        plt.title('Plan true model.')
        plt.gca().set_aspect('equal')
        plt.xlabel('x')
        plt.ylabel('y')
        plt.gca().set_aspect('equal', adjustable='box')

        # Vertical section
        ax = plt.subplot(326)
        mesh.plotSlice(m_true, ax=ax, normal='Y', ind=midx,
                       grid=True, clim=(model.min(), model.max()))
        plt.plot(([mesh.vectorCCx[0], mesh.vectorCCx[-1]]),
                 ([mesh.vectorCCz[zpanel], mesh.vectorCCz[zpanel]]), color='w')
        plt.title('E-W true model.')
        plt.gca().set_aspect('equal')
        plt.xlabel('x')
        plt.ylabel('z')
        plt.gca().set_aspect('equal', adjustable='box')

        # Plot convergence curves
        fig, axs = plt.figure(), plt.subplot()
        axs.plot(saveDict.phi_d, 'k', lw=2)
        axs.plot(
            np.r_[IRLS.iterStart, IRLS.iterStart],
            np.r_[0, np.max(saveDict.phi_d)], 'k:'
        )

        twin = axs.twinx()
        twin.plot(saveDict.phi_m, 'k--', lw=2)
        axs.text(
            IRLS.iterStart, 0,
            'IRLS Steps', va='bottom', ha='center',
            rotation='vertical', size=12,
            bbox={'facecolor': 'white'}
        )

        axs.set_ylabel('$\phi_d$', size=16, rotation=0)
        axs.set_xlabel('Iterations', size=14)
        twin.set_ylabel('$\phi_m$', size=16, rotation=0)

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

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

Gallery generated by Sphinx-Gallery