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
  • plot sumMap
  • plot sumMap
/home/vsts/work/1/s/simpeg/optimization.py:1422: FutureWarning:

InexactCG.tolCG has been deprecated, please use InexactCG.cg_atol. It will be removed in version 0.26.0 of SimPEG.

/home/vsts/work/1/s/simpeg/optimization.py:943: FutureWarning:

InexactCG.maxIterCG has been deprecated, please use InexactCG.cg_maxiter. It will be removed in version 0.26.0 of SimPEG.


Running inversion with SimPEG v0.24.1.dev29+g256406f92
model has any nan: 0
================================================= Projected GNCG =================================================
  #     beta     phi_d     phi_m       f      |proj(x-g)-x|  LS    iterCG   CG |Ax-b|/|b|  CG |Ax-b|   Comment
-----------------------------------------------------------------------------------------------------------------
x0 has any nan: 0
   0  5.71e+05  9.11e+06  4.27e-04  9.11e+06    6.29e+01      0       0          inf          inf
   1  2.85e+05  1.30e+05  7.63e-02  1.52e+05    3.49e+01      0      10       2.85e-05     3.92e+04
   2  1.43e+05  4.62e+04  2.32e-01  7.93e+04    4.95e+01      0      10       4.74e-03     2.14e+05
   3  7.14e+04  1.98e+04  3.61e-01  4.56e+04    5.52e+01      0      10       4.11e-03     6.30e+03
   4  3.57e+04  6.38e+03  4.74e-01  2.33e+04    4.82e+01      0      10       3.91e-04     4.09e+03
   5  1.78e+04  4.63e+03  5.09e-01  1.37e+04    3.56e+01      1      10       1.04e-02     5.30e+03
   6  8.92e+03  9.29e+02  5.99e-01  6.28e+03    4.86e+01      0      10       5.26e-04     5.54e+03
   7  4.46e+03  9.19e+02  6.00e-01  3.60e+03    5.29e+01      4      10       1.46e-01     2.22e+04
   8  2.23e+03  4.39e+02  6.45e-01  1.88e+03    3.36e+01      0      10       3.26e-02     3.73e+04
   9  1.11e+03  4.39e+02  6.45e-01  1.16e+03    3.54e+01      7      10       1.77e-01     4.71e+04
Reached starting chifact with l2-norm regularization: Start IRLS steps...
irls_threshold 0.010208932913106157
irls_threshold 0.012044322391514583
  10  1.11e+03  3.85e+02  9.01e-01  1.39e+03    3.42e+01      0      10       2.38e-01     7.66e+04
  11  1.11e+03  3.94e+02  9.58e-01  1.46e+03    6.17e+01      1      10       1.98e-01     6.49e+04
  12  1.11e+03  3.94e+02  1.00e+00  1.51e+03    6.17e+01     13      10       1.54e-01     1.82e+05
  13  1.11e+03  3.99e+02  1.02e+00  1.54e+03    3.19e+01      4      10       4.17e-02     4.92e+04
  14  1.11e+03  4.09e+02  1.03e+00  1.55e+03    3.69e+01      0      10       4.76e-01     1.19e+05
  15  1.11e+03  4.05e+02  1.02e+00  1.54e+03    6.08e+01      2      10       1.09e-02     1.79e+04
  16  1.11e+03  4.05e+02  9.93e-01  1.51e+03    6.08e+01     11      10       6.66e-02     7.00e+04
  17  1.11e+03  4.11e+02  9.49e-01  1.47e+03    3.21e+01      3      10       1.32e-01     1.39e+05
  18  1.11e+03  4.12e+02  9.00e-01  1.42e+03    3.13e+01      0      10       2.26e-01     6.90e+04
  19  1.11e+03  4.19e+02  8.41e-01  1.36e+03    6.11e+01      2      10       2.59e-01     5.80e+04
  20  1.11e+03  4.26e+02  7.90e-01  1.31e+03    2.88e+01      0      10       3.93e-02     4.11e+04
  21  1.11e+03  4.28e+02  7.29e-01  1.24e+03    3.39e+01      4      10       4.32e-01     3.39e+04
  22  1.11e+03  4.39e+02  6.76e-01  1.19e+03    5.35e+01      0      10       3.95e-02     1.89e+04
  23  9.05e+02  4.58e+02  6.06e-01  1.01e+03    6.18e+01      1      10       3.71e-01     6.20e+04
  24  7.37e+02  4.55e+02  5.51e-01  8.61e+02    3.73e+01      0      10       7.95e-03     1.28e+04
  25  6.00e+02  4.56e+02  4.93e-01  7.52e+02    3.74e+01      2      10       7.45e-03     1.46e+04
  26  6.00e+02  4.33e+02  4.35e-01  6.94e+02    6.04e+01      1      10       6.43e-03     1.36e+04
  27  6.00e+02  4.10e+02  3.69e-01  6.31e+02    3.21e+01      0      10       3.27e-02     2.69e+04
  28  6.00e+02  4.13e+02  3.15e-01  6.02e+02    3.70e+01      1      10       1.03e-01     7.59e+03
  29  6.00e+02  4.04e+02  2.60e-01  5.60e+02    5.99e+01      0      10       8.06e-03     7.65e+03
Reach maximum number of IRLS cycles: 20
------------------------- STOP! -------------------------
1 : |fc-fOld| = 0.0000e+00 <= tolF*(1+|f0|) = 9.1102e+05
1 : |xc-x_last| = 3.5517e-03 <= tolX*(1+|x0|) = 1.0075e-01
0 : |proj(x-g)-x|    = 5.9937e+01 <= tolG          = 1.0000e-03
0 : |proj(x-g)-x|    = 5.9937e+01 <= 1e3*eps       = 1.0000e-03
0 : maxIter   =     100    <= iter          =     30
------------------------- DONE! -------------------------

from discretize import TensorMesh
from discretize.utils import active_from_xyz
from simpeg import (
    utils,
    maps,
    regularization,
    data_misfit,
    optimization,
    inverse_problem,
    directives,
    inversion,
)
from simpeg.potential_fields import magnetics
import numpy as np
import matplotlib.pyplot as plt


def run(plotIt=True):
    h0_amplitude, h0_inclination, h0_declination = (50000.0, 90.0, 0.0)

    # Create a mesh
    dx = 5.0

    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 = TensorMesh([hxind, hyind, hzind], "CCC")

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

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

    # Go from topo to array of indices of active cells
    actv = active_from_xyz(mesh, topo, "N")
    nC = int(actv.sum())
    # Create and array of observation points
    xr = np.linspace(-20.0, 20.0, 20)
    yr = np.linspace(-20.0, 20.0, 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.nodes_z[-1] + 5.0

    # Create a MAGsurvey
    rxLoc = np.c_[utils.mkvc(X.T), utils.mkvc(Y.T), utils.mkvc(Z.T)]
    rxLoc = magnetics.Point(rxLoc)
    srcField = magnetics.UniformBackgroundField(
        receiver_list=[rxLoc],
        amplitude=h0_amplitude,
        inclination=h0_inclination,
        declination=h0_declination,
    )
    survey = magnetics.Survey(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.model_builder.add_block(
        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)

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

    # Create the forward model operator
    prob = magnetics.Simulation3DIntegral(
        mesh,
        survey=survey,
        chiMap=idenMap,
        active_cells=actv,
        store_sensitivities="forward_only",
    )

    # Compute linear forward operator and compute some data
    data = prob.make_synthetic_data(
        model, relative_error=0.0, noise_floor=1, add_noise=True
    )

    # Create 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", nC))

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

    # Create the forward model operator
    prob = magnetics.Simulation3DIntegral(
        mesh, survey=survey, chiMap=sumMap, active_cells=actv, store_sensitivities="ram"
    )

    # Make sensitivity weighting
    # Take the cell number out of the scaling.
    # Want to keep high sens for large volumes
    wr = (
        prob.getJtJdiag(np.ones(sumMap.shape[1]))
        / np.r_[homogMap.P.T * mesh.cell_volumes[actv], mesh.cell_volumes[actv]] ** 2.0
    )
    # Scale the model spaces independently
    wr[wires.homo.index] /= np.max((wires.homo * wr)) * utils.mkvc(
        homogMap.P.sum(axis=0).flatten()
    )
    wr[wires.hetero.index] /= np.max(wires.hetero * wr)
    wr = wr**0.5

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

    reg_m1 = regularization.Sparse(regMesh, mapping=wires.homo)
    reg_m1.set_weights(weights=wires.homo * wr)

    reg_m1.norms = [0, 2]
    reg_m1.reference_model = np.zeros(sumMap.shape[1])

    # Regularization for the voxel model
    reg_m2 = regularization.Sparse(
        mesh, active_cells=actv, mapping=wires.hetero, gradient_type="components"
    )
    reg_m2.set_weights(weights=wires.hetero * wr)

    reg_m2.norms = [0, 0, 0, 0]
    reg_m2.reference_model = np.zeros(sumMap.shape[1])

    reg = reg_m1 + reg_m2

    # Data misfit function
    dmis = data_misfit.L2DataMisfit(simulation=prob, data=data)

    # Add directives to the inversion
    opt = optimization.ProjectedGNCG(
        maxIter=100,
        lower=0.0,
        upper=1.0,
        maxIterLS=20,
        maxIterCG=10,
        tolCG=1e-3,
        tolG=1e-3,
        eps=1e-6,
    )
    invProb = inverse_problem.BaseInvProblem(dmis, reg, opt)
    betaest = directives.BetaEstimate_ByEig(beta0_ratio=1e-2)

    # Here is where the norms are applied
    # Use pick a threshold parameter empirically based on the distribution of
    #  model parameters
    IRLS = directives.UpdateIRLS(f_min_change=1e-3)

    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,
            pcolor_opts={"cmap": "inferno_r"},
            transparent="slider",
        )

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


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

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

Estimated memory usage: 314 MB

Gallery generated by Sphinx-Gallery