Note
Go to the end to download the full example code.
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
/home/vsts/work/1/s/simpeg/utils/model_builder.py:37: BreakingChangeWarning:
Since SimPEG v0.25.0, the 'get_indices_block' function returns a single array with the cell indices, instead of a tuple with a single element. This means that we don't need to unpack the tuple anymore to access to the cell indices.
If you were using this function as in:
ind = get_indices_block(p0, p1, mesh.cell_centers)[0]
Make sure you update it to:
ind = get_indices_block(p0, p1, mesh.cell_centers)
To hide this warning, add this to your script or notebook:
import warnings
from simpeg.utils import BreakingChangeWarning
warnings.filterwarnings(action='ignore', category=BreakingChangeWarning)
Running inversion with SimPEG v0.25.1.dev6+gb840df108
================================================= Projected GNCG =================================================
# beta phi_d phi_m f |proj(x-g)-x| LS iter_CG CG |Ax-b|/|b| CG |Ax-b| Comment
-----------------------------------------------------------------------------------------------------------------
0 5.32e+05 9.10e+06 4.27e-04 9.10e+06 0 inf inf
1 5.32e+05 1.29e+05 7.89e-02 1.71e+05 6.29e+01 0 3 7.66e-04 1.05e+06
2 2.66e+05 4.21e+04 2.44e-01 1.07e+05 3.53e+01 0 9 5.79e-04 2.59e+04
3 1.33e+05 1.79e+04 3.71e-01 6.73e+04 4.92e+01 0 10 1.26e-02 1.77e+04
4 6.66e+04 5.67e+03 4.80e-01 3.76e+04 5.56e+01 0 9 4.52e-04 4.89e+03
5 3.33e+04 4.27e+03 5.12e-01 2.13e+04 4.73e+01 1 10 1.15e-02 5.59e+03
6 1.66e+04 8.45e+02 5.99e-01 1.08e+04 3.56e+01 0 10 4.85e-04 5.48e+03
7 8.32e+03 8.38e+02 6.00e-01 5.83e+03 4.84e+01 4 10 1.22e+00 1.74e+05
8 4.16e+03 4.12e+02 6.42e-01 3.08e+03 5.35e+01 0 10 3.20e-02 3.78e+04
9 2.08e+03 4.12e+02 6.42e-01 1.75e+03 3.33e+01 7 10 2.39e-01 6.67e+04
10 1.04e+03 3.61e+02 6.61e-01 1.05e+03 3.51e+01 0 10 7.19e-02 2.46e+04
Reached starting chifact with l2-norm regularization: Start IRLS steps...
irls_threshold 0.010219925373431072
irls_threshold 0.012367468447233831
11 1.04e+03 3.66e+02 9.08e-01 1.31e+03 5.26e+01 2 10 4.22e-01 5.69e+04
12 1.04e+03 3.70e+02 9.67e-01 1.37e+03 6.12e+01 3 10 1.71e-01 1.42e+05
13 1.04e+03 3.70e+02 1.01e+00 1.42e+03 3.64e+01 12 10 1.01e-01 9.47e+04
14 1.04e+03 3.71e+02 1.04e+00 1.45e+03 3.64e+01 4 10 1.04e-01 9.81e+04
15 1.04e+03 3.79e+02 1.01e+00 1.43e+03 5.43e+01 0 10 4.88e-01 1.72e+05
16 1.04e+03 3.79e+02 1.00e+00 1.42e+03 5.93e+01 4 10 2.87e-01 1.59e+05
17 1.04e+03 3.80e+02 9.77e-01 1.40e+03 2.98e+01 4 10 1.77e-01 2.94e+04
18 1.04e+03 4.01e+02 9.15e-01 1.35e+03 3.34e+01 0 10 3.02e-01 1.48e+05
19 1.04e+03 4.02e+02 8.76e-01 1.31e+03 6.23e+01 7 10 2.08e-02 4.35e+04
20 1.04e+03 4.07e+02 8.00e-01 1.24e+03 6.23e+01 2 10 4.24e-02 8.90e+04
21 1.04e+03 3.99e+02 7.26e-01 1.15e+03 3.48e+01 0 10 1.75e-01 1.44e+05
22 1.04e+03 4.09e+02 6.73e-01 1.11e+03 3.39e+01 2 10 1.01e-01 2.15e+04
23 1.04e+03 4.12e+02 5.94e-01 1.03e+03 6.15e+01 0 10 8.57e-03 9.98e+03
24 1.04e+03 4.16e+02 5.45e-01 9.83e+02 3.29e+01 3 10 3.87e-02 9.47e+03
25 1.04e+03 4.29e+02 4.66e-01 9.13e+02 3.52e+01 0 10 1.69e-02 1.16e+04
26 1.04e+03 4.31e+02 3.89e-01 8.35e+02 6.18e+01 0 10 7.33e-03 9.52e+03
27 1.04e+03 4.29e+02 3.32e-01 7.74e+02 6.16e+01 1 10 6.76e-03 8.29e+03
28 1.04e+03 4.22e+02 2.83e-01 7.16e+02 3.53e+01 0 10 1.50e-02 1.16e+04
29 1.04e+03 4.34e+02 2.36e-01 6.80e+02 4.81e+01 0 10 2.61e-01 9.25e+03
30 1.04e+03 4.28e+02 2.04e-01 6.40e+02 6.19e+01 0 10 4.87e-03 6.07e+03
Reach maximum number of IRLS cycles: 20
------------------------- STOP! -------------------------
1 : |fc-fOld| = 1.5172e+01 <= tolF*(1+|f0|) = 9.1019e+05
1 : |xc-x_last| = 4.1430e-03 <= tolX*(1+|x0|) = 1.0075e-01
0 : |proj(x-g)-x| = 6.1934e+01 <= tolG = 1.0000e-03
0 : |proj(x-g)-x| = 6.1934e+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,
cg_maxiter=10,
cg_rtol=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 25.829 seconds)
Estimated memory usage: 339 MB