Maps: ComboMapsΒΆ

Invert synthetic magnetic data with variable background values and a single block anomaly buried at depth. We will use the Sum Map to invert for both the background values and an heterogeneous susceptibiilty model.

1

  • ../../../_images/sphx_glr_plot_sumMap_001.png
  • ../../../_images/sphx_glr_plot_sumMap_002.png

Out:

Begin forward: M=H0, Rx type= tmi
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 %
Begin forward: M=H0, Rx type= tmi
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 %
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***
Approximated diag(JtJ) with linear operator
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.73e+13  4.55e+06  1.72e-06  3.43e+07    7.47e-03      0
   1  8.64e+12  4.71e+06  7.81e-10  4.72e+06    6.26e+01      0
   2  4.32e+12  4.69e+06  3.08e-09  4.70e+06    6.37e+01      0
   3  2.16e+12  4.63e+06  1.22e-08  4.66e+06    6.35e+01      0   Skip BFGS
   4  1.08e+12  4.53e+06  4.75e-08  4.58e+06    6.32e+01      0   Skip BFGS
   5  5.40e+11  4.34e+06  1.80e-07  4.44e+06    6.29e+01      0   Skip BFGS
   6  2.70e+11  4.00e+06  6.55e-07  4.17e+06    6.28e+01      0   Skip BFGS
   7  1.35e+11  3.44e+06  2.19e-06  3.73e+06    6.27e+01      0   Skip BFGS
   8  6.75e+10  2.66e+06  6.40e-06  3.09e+06    6.25e+01      0   Skip BFGS
   9  3.37e+10  1.80e+06  1.55e-05  2.33e+06    6.21e+01      0   Skip BFGS
  10  1.69e+10  1.08e+06  3.06e-05  1.60e+06    6.16e+01      0   Skip BFGS
  11  8.43e+09  6.09e+05  5.03e-05  1.03e+06    6.04e+01      0   Skip BFGS
  12  4.22e+09  3.41e+05  7.24e-05  6.47e+05    5.87e+01      0   Skip BFGS
  13  2.11e+09  2.11e+05  9.38e-05  4.08e+05    5.73e+01      0   Skip BFGS
  14  1.05e+09  1.55e+05  1.12e-04  2.73e+05    5.59e+01      0   Skip BFGS
  15  5.27e+08  1.31e+05  1.28e-04  1.98e+05    5.41e+01      0   Skip BFGS
  16  2.64e+08  1.15e+05  1.50e-04  1.54e+05    5.17e+01      0   Skip BFGS
  17  1.32e+08  9.51e+04  2.05e-04  1.22e+05    5.03e+01      0   Skip BFGS
  18  6.59e+07  6.95e+04  3.46e-04  9.23e+04    4.95e+01      0   Skip BFGS
  19  3.29e+07  4.19e+04  6.46e-04  6.32e+04    4.93e+01      0   Skip BFGS
  20  1.65e+07  2.02e+04  1.11e-03  3.85e+04    4.83e+01      0   Skip BFGS
  21  8.24e+06  7.90e+03  1.62e-03  2.13e+04    4.83e+01      0   Skip BFGS
  22  4.12e+06  2.73e+03  2.05e-03  1.12e+04    4.81e+01      0   Skip BFGS
  23  2.06e+06  9.44e+02  2.34e-03  5.77e+03    4.79e+01      0   Skip BFGS
  24  1.03e+06  4.02e+02  2.52e-03  3.00e+03    4.83e+01      0   Skip BFGS
  25  5.15e+05  2.42e+02  2.62e-03  1.59e+03    4.77e+01      0   Skip BFGS
Reached starting chifact with l2-norm regularization: Start IRLS steps...
eps_p: 0.010155126601143325 eps_q: 0.010155126601143325
eps_p: 0.011973340358565595 eps_q: 0.011973340358565595
delta phim: 3.476e+02
  26  2.57e+05  1.93e+02  3.70e-03  1.15e+03    3.09e+01      0   Skip BFGS
delta phim: 6.290e+03
  27  2.57e+05  1.87e+02  3.94e-03  1.20e+03    2.60e+01      0
delta phim: 2.550e-01
  28  2.57e+05  1.87e+02  4.07e-03  1.24e+03    2.19e+01      0
delta phim: 1.632e-01
  29  2.57e+05  1.95e+02  4.10e-03  1.25e+03    5.09e+01      0
delta phim: 1.555e-01
  30  2.57e+05  1.99e+02  4.10e-03  1.25e+03    2.44e+01      0
delta phim: 1.065e-01
  31  2.57e+05  2.01e+02  4.00e-03  1.23e+03    3.21e+01      3
delta phim: 2.078e-01
  32  2.57e+05  2.12e+02  3.89e-03  1.21e+03    2.10e+01      0
delta phim: 2.829e-03
  33  2.57e+05  2.14e+02  3.73e-03  1.17e+03    5.81e+01      3   Skip BFGS
delta phim: 1.683e-01
  34  2.57e+05  2.19e+02  3.54e-03  1.13e+03    2.63e+01      0
delta phim: 4.538e-02
  35  2.57e+05  2.20e+02  3.37e-03  1.09e+03    3.18e+01      4   Skip BFGS
Beta search step
  35  2.11e+05  2.20e+02  3.37e-03  9.31e+02    3.20e+01      0
delta phim: 1.447e-01
  36  2.11e+05  2.07e+02  3.20e-03  8.82e+02    4.69e+01      0   Skip BFGS
delta phim: 9.432e-03
  37  2.11e+05  2.10e+02  3.02e-03  8.48e+02    5.93e+01      2   Skip BFGS
delta phim: 9.755e-02
  38  2.11e+05  2.15e+02  2.80e-03  8.07e+02    3.52e+01      0
delta phim: 5.375e-02
  39  2.11e+05  2.16e+02  2.67e-03  7.79e+02    3.53e+01      9   Skip BFGS
delta phim: 1.107e-01
  40  2.11e+05  2.19e+02  2.51e-03  7.48e+02    6.06e+01      1
delta phim: 2.167e-02
  41  2.11e+05  2.13e+02  2.36e-03  7.13e+02    2.86e+01      0
delta phim: 6.997e-02
  42  2.11e+05  2.15e+02  2.26e-03  6.93e+02    3.48e+01      2   Skip BFGS
delta phim: 7.203e-02
  43  2.11e+05  2.11e+02  2.15e-03  6.66e+02    2.77e+01      0
delta phim: 4.761e-02
  44  2.11e+05  2.12e+02  2.06e-03  6.48e+02    6.02e+01      1   Skip BFGS
delta phim: 4.355e-02
  45  2.11e+05  2.06e+02  1.96e-03  6.20e+02    2.85e+01      0
Reach maximum number of IRLS cycles: 20
------------------------- STOP! -------------------------
1 : |fc-fOld| = 0.0000e+00 <= tolF*(1+|f0|) = 3.4255e+06
1 : |xc-x_last| = 1.1552e-03 <= tolX*(1+|x0|) = 1.0075e-01
0 : |proj(x-g)-x|    = 2.8474e+01 <= tolG          = 1.0000e-03
0 : |proj(x-g)-x|    = 2.8474e+01 <= 1e3*eps       = 1.0000e-03
0 : maxIter   =     100    <= iter          =     46
------------------------- DONE! -------------------------

from SimPEG import (
    Mesh, Utils, Maps, Regularization,
    DataMisfit, Optimization, InvProblem,
    Directives, Inversion, PF
)
import numpy as np
import matplotlib.pyplot as plt


def run(plotIt=True):

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

    # 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
    model = np.zeros(mesh.nC)

    # Change values in half the domain
    model[mesh.gridCC[:,0] < 0] = 0.01

    # Add a block in half-space
    model = Utils.ModelBuilder.addBlock(mesh.gridCC, model, np.r_[-10,-10,20], np.r_[10,10,40], 0.05)

    model = Utils.mkvc(model)
    model = model[actv]


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

    # 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

    # Plot the data
    rxLoc = survey.srcField.rxList[0].locs

    # Creat a homogenous maps for the two domains
    domains = [mesh.gridCC[actv,0] < 0, mesh.gridCC[actv,0] >= 0]
    homogMap = Maps.SurjectUnits(domains)

    # Create a wire map for a second model space, voxel based
    wires = Maps.Wires(('homo', len(domains)), ('hetero', len(actv)))

    # Create Sum map
    sumMap = Maps.SumMap([homogMap*wires.homo, wires.hetero])

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

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

    # Make depth weighting
    wr = np.zeros(sumMap.shape[1])

    # Take the cell number out of the scaling.
    # Want to keep high sens for large volumnes
    scale = Utils.sdiag(np.r_[Utils.mkvc(1./homogMap.P.sum(axis=0)),np.ones_like(actv)])

    for ii in range(survey.nD):
        wr += ((prob.G[ii, :]*prob.chiMap.deriv(np.ones(sumMap.shape[1])*1e-4)*scale)/survey.std[ii])**2.

    # Scale the model spaces independently
    wr[wires.homo.index] /= (np.max((wires.homo*wr)))
    wr[wires.hetero.index] /= (np.max(wires.hetero*wr))
    wr = wr**0.5

    ## Create a regularization
    # For the homogeneous model
    regMesh = Mesh.TensorMesh([len(domains)])

    reg_m1 = Regularization.Sparse(regMesh, mapping=wires.homo)
    reg_m1.cell_weights = wires.homo*wr
    reg_m1.norms = np.c_[0, 2, 2, 2]
    reg_m1.mref = np.zeros(sumMap.shape[1])

    # Regularization for the voxel model
    reg_m2 = Regularization.Sparse(mesh, indActive=actv, mapping=wires.hetero)
    reg_m2.cell_weights = wires.hetero*wr
    reg_m2.norms = np.c_[0, 1, 1, 1]
    reg_m2.mref =  np.zeros(sumMap.shape[1])

    reg = reg_m1 + reg_m2

    # 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, tolG=1e-3, eps=1e-6)
    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-3, minGNiter=1)
    update_Jacobi = Directives.UpdatePreconditioner()
    inv = Inversion.BaseInversion(invProb,
                                  directiveList=[IRLS, betaest, update_Jacobi])

    # Run the inversion
    m0 = np.ones(sumMap.shape[1])*1e-4  # Starting model
    prob.model = m0
    mrecSum = inv.run(m0)
    if plotIt:

        mesh.plot_3d_slicer(actvMap * model, aspect="equal", zslice=30, pcolorOpts={"cmap":'inferno_r'}, transparent='slider')

        mesh.plot_3d_slicer(actvMap * sumMap * mrecSum, aspect="equal", zslice=30, pcolorOpts={"cmap":'inferno_r'}, transparent='slider')



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

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

Gallery generated by Sphinx-Gallery