PF: Gravity: Inversion LinearΒΆ

Create a synthetic block model and invert with a compact norm

  • ../../../_images/sphx_glr_plot_inversion_linear_0011.png
  • ../../../_images/sphx_glr_plot_inversion_linear_002.png
  • ../../../_images/sphx_glr_plot_inversion_linear_003.png

Out:

Begin calculation of forward operator: z
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 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
=============================== Projected GNCG ===============================
  #     beta     phi_d     phi_m       f      |proj(x-g)-x|  LS    Comment
-----------------------------------------------------------------------------
x0 has any nan: 0
   0  6.04e+07  1.73e+05  0.00e+00  1.73e+05    1.00e+02      0
   1  3.02e+07  1.71e+05  1.36e-05  1.72e+05    1.00e+02      0
   2  1.51e+07  1.70e+05  5.37e-05  1.71e+05    1.00e+02      0   Skip BFGS
   3  7.55e+06  1.67e+05  2.09e-04  1.68e+05    1.00e+02      0   Skip BFGS
   4  3.77e+06  1.61e+05  7.96e-04  1.64e+05    1.00e+02      0   Skip BFGS
   5  1.89e+06  1.50e+05  2.90e-03  1.56e+05    1.00e+02      0   Skip BFGS
   6  9.43e+05  1.33e+05  9.80e-03  1.42e+05    1.00e+02      0   Skip BFGS
   7  4.72e+05  1.07e+05  2.93e-02  1.21e+05    9.98e+01      0   Skip BFGS
   8  2.36e+05  7.81e+04  7.41e-02  9.56e+04    9.94e+01      0   Skip BFGS
   9  1.18e+05  5.15e+04  1.54e-01  6.97e+04    9.91e+01      0   Skip BFGS
  10  5.90e+04  3.15e+04  2.74e-01  4.77e+04    9.83e+01      0   Skip BFGS
  11  2.95e+04  1.76e+04  4.41e-01  3.06e+04    9.71e+01      0   Skip BFGS
  12  1.47e+04  8.60e+03  6.55e-01  1.83e+04    9.57e+01      0   Skip BFGS
  13  7.37e+03  3.64e+03  8.88e-01  1.02e+04    9.37e+01      0   Skip BFGS
  14  3.69e+03  1.40e+03  1.10e+00  5.44e+03    9.04e+01      0   Skip BFGS
  15  1.84e+03  5.52e+02  1.25e+00  2.86e+03    8.54e+01      0   Skip BFGS
  16  9.21e+02  2.67e+02  1.36e+00  1.52e+03    7.75e+01      0   Skip BFGS
Reached starting chifact with l2-norm regularization: Start IRLS steps...
eps_p: 0.17248982883996306 eps_q: 0.17248982883996306
Eps_p: 0.11499321922664203
Eps_q: 0.11499321922664203
delta phim:    inf
  17  4.61e+02  1.75e+02  1.85e+00  1.03e+03    4.10e+01      0   Skip BFGS
Beta search step
  18  7.78e+02  1.75e+02  1.85e+00  1.62e+03    6.69e+01      0   Skip BFGS
Beta search step
  19  1.22e+03  1.75e+02  1.85e+00  2.44e+03    8.53e+01      0
Beta search step
  20  1.00e+03  1.75e+02  1.85e+00  2.03e+03    8.01e+01      0
Eps_p: 0.07666214615109469
Eps_q: 0.07666214615109469
delta phim: 1.901e+01
  21  1.00e+03  1.99e+02  2.27e+00  2.47e+03    7.90e+01      0   Skip BFGS
Beta search step
  22  8.24e+02  1.99e+02  2.27e+00  2.07e+03    6.21e+01      0
Eps_p: 0.05110809743406313
Eps_q: 0.05110809743406313
delta phim: 7.569e-01
  23  8.24e+02  1.98e+02  2.73e+00  2.45e+03    8.00e+01      0   Skip BFGS
Beta search step
  24  6.83e+02  1.98e+02  2.73e+00  2.07e+03    6.54e+01      0
Eps_p: 0.034072064956042085
Eps_q: 0.034072064956042085
delta phim: 7.149e-01
  25  6.83e+02  1.99e+02  3.14e+00  2.34e+03    8.16e+01      0   Skip BFGS
Eps_p: 0.022714709970694722
Eps_q: 0.022714709970694722
delta phim: 5.869e-01
  26  6.83e+02  2.11e+02  3.39e+00  2.52e+03    8.45e+01      0
Beta search step
  27  5.65e+02  2.11e+02  3.39e+00  2.12e+03    7.03e+01      0
Eps_p: 0.015143139980463148
Eps_q: 0.015143139980463148
delta phim: 4.385e-01
  28  5.65e+02  2.05e+02  3.55e+00  2.21e+03    8.35e+01      0   Skip BFGS
Eps_p: 0.010095426653642098
Eps_q: 0.010095426653642098
delta phim: 3.146e-01
  29  5.65e+02  2.14e+02  3.67e+00  2.29e+03    8.60e+01      0
Beta search step
  30  4.57e+02  2.14e+02  3.67e+00  1.89e+03    6.65e+01      0
Eps_p: 0.006730284435761399
Eps_q: 0.006730284435761399
delta phim: 2.079e-01
  31  4.57e+02  2.13e+02  3.85e+00  1.97e+03    8.79e+01      0   Skip BFGS
Beta search step
  32  3.72e+02  2.13e+02  3.85e+00  1.64e+03    7.28e+01      0
Eps_p: 0.004486856290507599
Eps_q: 0.004486856290507599
delta phim: 2.013e-01
  33  3.72e+02  2.10e+02  3.87e+00  1.65e+03    8.72e+01      0   Skip BFGS
Eps_p: 0.0029912375270050658
Eps_q: 0.0029912375270050658
delta phim: 2.332e-01
  34  3.72e+02  2.14e+02  3.80e+00  1.63e+03    8.72e+01      0
Eps_p: 0.0019941583513367104
Eps_q: 0.0019941583513367104
delta phim: 2.106e-01
  35  3.72e+02  2.18e+02  3.71e+00  1.59e+03    8.84e+01      0
Beta search step
  36  3.06e+02  2.18e+02  3.71e+00  1.35e+03    6.69e+01      0
Eps_p: 0.0013294389008911402
Eps_q: 0.0013294389008911402
delta phim: 2.100e-01
  37  3.06e+02  2.13e+02  3.61e+00  1.32e+03    8.48e+01      0   Skip BFGS
Eps_p: 0.0008862926005940935
Eps_q: 0.0008862926005940935
delta phim: 2.012e-01
  38  3.06e+02  2.12e+02  3.50e+00  1.28e+03    8.54e+01      0
Eps_p: 0.0005908617337293956
Eps_q: 0.0005908617337293956
delta phim: 1.470e-01
  39  3.06e+02  2.06e+02  3.41e+00  1.25e+03    8.54e+01      0
Eps_p: 0.00039390782248626376
Eps_q: 0.00039390782248626376
delta phim: 1.221e-01
  40  3.06e+02  1.98e+02  3.34e+00  1.22e+03    8.76e+01      0
Eps_p: 0.0002626052149908425
Eps_q: 0.0002626052149908425
delta phim: 1.036e-01
  41  3.06e+02  1.94e+02  3.29e+00  1.20e+03    8.84e+01      0
Eps_p: 0.00017507014332722836
Eps_q: 0.00017507014332722836
delta phim: 8.939e-02
  42  3.06e+02  1.91e+02  3.26e+00  1.19e+03    9.02e+01      0
Eps_p: 0.0001167134288848189
Eps_q: 0.0001167134288848189
delta phim: 6.854e-02
  43  3.06e+02  1.90e+02  3.23e+00  1.18e+03    9.04e+01      0
Eps_p: 7.780895258987926e-05
Eps_q: 7.780895258987926e-05
delta phim: 3.062e-02
  44  3.06e+02  1.89e+02  3.22e+00  1.17e+03    9.09e+01      0
Eps_p: 5.187263505991951e-05
Eps_q: 5.187263505991951e-05
delta phim: 1.760e-02
  45  3.06e+02  1.88e+02  3.21e+00  1.17e+03    9.18e+01      0
Eps_p: 3.4581756706613005e-05
Eps_q: 3.4581756706613005e-05
delta phim: 5.778e-02
  46  3.06e+02  1.88e+02  3.20e+00  1.17e+03    9.21e+01      0
Eps_p: 2.3054504471075336e-05
Eps_q: 2.3054504471075336e-05
delta phim: 8.552e-02
  47  3.06e+02  1.88e+02  3.19e+00  1.17e+03    9.26e+01      0
Eps_p: 1.5369669647383556e-05
Eps_q: 1.5369669647383556e-05
delta phim: 8.073e-02
  48  3.06e+02  1.88e+02  3.19e+00  1.16e+03    9.36e+01      0
Eps_p: 1.0246446431589038e-05
Eps_q: 1.0246446431589038e-05
delta phim: 5.729e-02
  49  3.06e+02  1.87e+02  3.19e+00  1.16e+03    9.32e+01      0
Eps_p: 6.830964287726025e-06
Eps_q: 6.830964287726025e-06
delta phim: 4.321e-02
  50  3.06e+02  1.87e+02  3.19e+00  1.16e+03    9.29e+01      0
Eps_p: 4.55397619181735e-06
Eps_q: 4.55397619181735e-06
delta phim: 1.430e-02
  51  3.06e+02  1.87e+02  3.19e+00  1.16e+03    9.30e+01      0
Eps_p: 3.0359841278782333e-06
Eps_q: 3.0359841278782333e-06
delta phim: 2.275e-02
  52  3.06e+02  1.87e+02  3.19e+00  1.16e+03    9.27e+01      0
Eps_p: 2.023989418585489e-06
Eps_q: 2.023989418585489e-06
delta phim: 3.204e-03
  53  3.06e+02  1.87e+02  3.19e+00  1.16e+03    9.28e+01      0
Eps_p: 1.3493262790569927e-06
Eps_q: 1.3493262790569927e-06
delta phim: 1.091e-02
  54  3.06e+02  1.87e+02  3.19e+00  1.16e+03    9.25e+01      0
Eps_p: 8.995508527046618e-07
Eps_q: 8.995508527046618e-07
delta phim: 9.917e-03
  55  3.06e+02  1.87e+02  3.19e+00  1.16e+03    9.25e+01      0
Reach maximum number of IRLS cycles: 30
------------------------- STOP! -------------------------
1 : |fc-fOld| = 0.0000e+00 <= tolF*(1+|f0|) = 1.7309e+04
1 : |xc-x_last| = 7.2770e-04 <= tolX*(1+|x0|) = 1.0100e-01
0 : |proj(x-g)-x|    = 9.2543e+01 <= tolG          = 1.0000e-01
0 : |proj(x-g)-x|    = 9.2543e+01 <= 1e3*eps       = 1.0000e-02
0 : maxIter   =     100    <= iter          =     56
------------------------- DONE! -------------------------

import numpy as np
import matplotlib.pyplot as plt

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

    # Create a mesh
    dx = 5.

    hxind = [(dx, 5, -1.3), (dx, 15), (dx, 5, 1.3)]
    hyind = [(dx, 5, -1.3), (dx, 15), (dx, 5, 1.3)]
    hzind = [(dx, 5, -1.3), (dx, 7), (3.5, 1), (2, 5)]

    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(-30., 30., 20)
    yr = np.linspace(-30., 30., 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] + 0.1

    # Create a MAGsurvey
    rxLoc = np.c_[Utils.mkvc(X.T), Utils.mkvc(Y.T), Utils.mkvc(Z.T)]
    rxLoc = PF.BaseGrav.RxObs(rxLoc)
    srcField = PF.BaseGrav.SrcField([rxLoc])
    survey = PF.BaseGrav.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-5):(midx-1), (midy-2):(midy+2), -10:-6] = 0.75
    model[(midx+1):(midx+5), (midy-2):(midy+2), -10:-6] = -0.75
    model = Utils.mkvc(model)
    model = model[actv]

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

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

    # Create the forward model operator
    prob = PF.Gravity.GravityIntegral(mesh, rhoMap=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))*1e-3
    wd = np.ones(len(data))*1e-3  # 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.norms = np.c_[1, 0, 0, 0]

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

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

    # 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=30, coolEpsFact=1.5, beta_tol=1e-1,
    )
    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 = -7
        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

        vmin, vmax = mrec.min(), mrec.max()

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

        plt.figure()

        # Plot L2 model
        ax = plt.subplot(321)
        mesh.plotSlice(m_l2, ax=ax, normal='Z', ind=zpanel,
                       grid=True, clim=(vmin, vmax))
        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=(vmin, vmax))
        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=(vmin, vmax))
        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=(vmin, vmax))
        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=(vmin, vmax))
        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=(vmin, vmax))
        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, np.max(saveDict.phi_d)/2.,
            '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 29.290 seconds)

Gallery generated by Sphinx-Gallery