Note
Go to the end to download the full example code.
2D inversion of Loop-Loop EM Data#
In this example, we consider a single line of loop-loop EM data at 30kHz with 3 different coil separations [0.32m, 0.71m, 1.18m]. We will use only Horizontal co-planar orientations (vertical magnetic dipole), and look at the real and imaginary parts of the secondary magnetic field.
We use the simpeg.maps.Surject2Dto3D mapping to invert for a 2D model
and perform the forward modelling in 3D.
import numpy as np
import matplotlib.pyplot as plt
import time
try:
from pymatsolver import Pardiso as Solver
except ImportError:
from simpeg import SolverLU as Solver
import discretize
from simpeg import (
maps,
optimization,
data_misfit,
regularization,
inverse_problem,
inversion,
directives,
Report,
)
from simpeg.electromagnetics import frequency_domain as FDEM
Setup#
Define the survey and model parameters
sigma_surface = 10e-3
sigma_deep = 40e-3
sigma_air = 1e-8
coil_separations = [0.32, 0.71, 1.18]
freq = 30e3
print("skin_depth: {:1.2f}m".format(500 / np.sqrt(sigma_deep * freq)))
skin_depth: 14.43m
Define a dipping interface between the surface layer and the deeper layer
z_interface_shallow = -0.25
z_interface_deep = -1.5
x_dip = np.r_[0.0, 8.0]
def interface(x):
interface = np.zeros_like(x)
interface[x < x_dip[0]] = z_interface_shallow
dipping_unit = (x >= x_dip[0]) & (x <= x_dip[1])
x_dipping = (-(z_interface_shallow - z_interface_deep) / x_dip[1]) * (
x[dipping_unit]
) + z_interface_shallow
interface[dipping_unit] = x_dipping
interface[x > x_dip[1]] = z_interface_deep
return interface
Forward Modelling Mesh#
Here, we set up a 3D tensor mesh which we will perform the forward simulations on.
Note
In practice, a smaller horizontal discretization should be used to improve accuracy, particularly for the shortest offset (eg. you can try 0.25m).
csx = 0.5 # cell size for the horizontal direction
csz = 0.125 # cell size for the vertical direction
pf = 1.3 # expansion factor for the padding cells
npadx = 7 # number of padding cells in the x-direction
npady = 7 # number of padding cells in the y-direction
npadz = 11 # number of padding cells in the z-direction
core_domain_x = np.r_[-11.5, 11.5] # extent of uniform cells in the x-direction
core_domain_z = np.r_[-2.0, 0.0] # extent of uniform cells in the z-direction
# number of cells in the core region
ncx = int(np.diff(core_domain_x) / csx)
ncz = int(np.diff(core_domain_z) / csz)
# create a 3D tensor mesh
mesh = discretize.TensorMesh(
[
[(csx, npadx, -pf), (csx, ncx), (csx, npadx, pf)],
[(csx, npady, -pf), (csx, 1), (csx, npady, pf)],
[(csz, npadz, -pf), (csz, ncz), (csz, npadz, pf)],
]
)
# set the origin
mesh.x0 = np.r_[
-mesh.h[0].sum() / 2.0, -mesh.h[1].sum() / 2.0, -mesh.h[2][: npadz + ncz].sum()
]
print("the mesh has {} cells".format(mesh.nC))
mesh.plot_grid()
the mesh has 34200 cells
<Axes3D: xlabel='x1', ylabel='x2', zlabel='x3'>
Inversion Mesh#
Here, we set up a 2D tensor mesh which we will represent the inversion model on
inversion_mesh = discretize.TensorMesh([mesh.h[0], mesh.h[2][mesh.cell_centers_z <= 0]])
inversion_mesh.x0 = [-inversion_mesh.h[0].sum() / 2.0, -inversion_mesh.h[1].sum()]
inversion_mesh.plot_grid()
<Axes: xlabel='x1', ylabel='x2'>
Mappings#
Mappings are used to take the inversion model and represent it as electrical conductivity on the inversion mesh. We will invert for log-conductivity below the surface, fixing the conductivity of the air cells to 1e-8 S/m
# create a 2D mesh that includes air cells
mesh2D = discretize.TensorMesh([mesh.h[0], mesh.h[2]], x0=mesh.x0[[0, 2]])
active_inds = mesh2D.gridCC[:, 1] < 0 # active indices are below the surface
mapping = (
maps.Surject2Dto3D(mesh)
* maps.InjectActiveCells( # populates 3D space from a 2D model
mesh2D, active_inds, sigma_air
)
* maps.ExpMap( # adds air cells
nP=inversion_mesh.nC
) # takes the exponential (log(sigma) --> sigma)
)
True Model#
Create our true model which we will use to generate synthetic data for
m_true = np.log(sigma_deep) * np.ones(inversion_mesh.nC)
interface_depth = interface(inversion_mesh.gridCC[:, 0])
m_true[inversion_mesh.gridCC[:, 1] > interface_depth] = np.log(sigma_surface)
fig, ax = plt.subplots(1, 1)
cb = plt.colorbar(inversion_mesh.plot_image(m_true, ax=ax, grid=True)[0], ax=ax)
cb.set_label(r"$\log(\sigma)$")
ax.set_title("true model")
ax.set_xlim([-10, 10])
ax.set_ylim([-2, 0])
(-2.0, 0.0)
Survey#
Create our true model which we will use to generate synthetic data for
src_locations = np.arange(-11, 11, 0.5)
src_z = 0.25 # src is 0.25m above the surface
orientation = "z" # z-oriented dipole for horizontal co-planar loops
# reciever offset in 3D space
rx_offsets = np.vstack([np.r_[sep, 0.0, 0.0] for sep in coil_separations])
# create our source list - one source per location
source_list = []
for x in src_locations:
src_loc = np.r_[x, 0.0, src_z]
rx_locs = src_loc - rx_offsets
rx_real = FDEM.Rx.PointMagneticFluxDensitySecondary(
locations=rx_locs, orientation=orientation, component="real"
)
rx_imag = FDEM.Rx.PointMagneticFluxDensitySecondary(
locations=rx_locs, orientation=orientation, component="imag"
)
src = FDEM.Src.MagDipole(
receiver_list=[rx_real, rx_imag],
location=src_loc,
orientation=orientation,
frequency=freq,
)
source_list.append(src)
# create the survey and problem objects for running the forward simulation
survey = FDEM.Survey(source_list)
prob = FDEM.Simulation3DMagneticFluxDensity(
mesh, survey=survey, sigmaMap=mapping, solver=Solver
)
Set up data for inversion#
Generate clean, synthetic data. Later we will invert the clean data, and assign a standard deviation of 0.05, and a floor of 1e-11.
t = time.time()
data = prob.make_synthetic_data(
m_true, relative_error=0.05, noise_floor=1e-11, add_noise=False
)
dclean = data.dclean
print("Done forward simulation. Elapsed time = {:1.2f} s".format(time.time() - t))
def plot_data(data, ax=None, color="C0", label=""):
if ax is None:
fig, ax = plt.subplots(1, 3, figsize=(15, 5))
# data is [re, im, re, im, ...]
data_real = data[0::2]
data_imag = data[1::2]
for i, offset in enumerate(coil_separations):
ax[i].plot(
src_locations,
data_real[i :: len(coil_separations)],
color=color,
label="{} real".format(label),
)
ax[i].plot(
src_locations,
data_imag[i :: len(coil_separations)],
"--",
color=color,
label="{} imag".format(label),
)
ax[i].set_title("offset = {:1.2f}m".format(offset))
ax[i].legend()
ax[i].grid(which="both")
ax[i].set_ylim(np.r_[data.min(), data.max()] + 1e-11 * np.r_[-1, 1])
ax[i].set_xlabel("source location x (m)")
ax[i].set_ylabel("Secondary B-Field (T)")
plt.tight_layout()
return ax
ax = plot_data(dclean)
Done forward simulation. Elapsed time = 12.61 s
Set up the inversion#
We create the data misfit, simple regularization
(a least-squares-style regularization, simpeg.regularization.LeastSquareRegularization)
The smoothness and smallness contributions can be set by including
alpha_s, alpha_x, alpha_y as input arguments when the regularization is
created. The default reference model in the regularization is the starting
model. To set something different, you can input an mref into the
regularization.
We estimate the trade-off parameter, beta, between the data
misfit and regularization by the largest eigenvalue of the data misfit and
the regularization. Here, we use a fixed beta, but could alternatively
employ a beta-cooling schedule using simpeg.directives.BetaSchedule
dmisfit = data_misfit.L2DataMisfit(simulation=prob, data=data)
reg = regularization.WeightedLeastSquares(inversion_mesh)
opt = optimization.InexactGaussNewton(maxIterCG=10, remember="xc")
invProb = inverse_problem.BaseInvProblem(dmisfit, reg, opt)
betaest = directives.BetaEstimate_ByEig(beta0_ratio=0.05, n_pw_iter=1, seed=1)
target = directives.TargetMisfit()
directiveList = [betaest, target]
inv = inversion.BaseInversion(invProb, directiveList=directiveList)
print("The target misfit is {:1.2f}".format(target.target))
The target misfit is 264.00
Run the inversion#
We start from a half-space equal to the deep conductivity.
Running inversion with SimPEG v0.22.1
simpeg.InvProblem will set Regularization.reference_model to m0.
simpeg.InvProblem will set Regularization.reference_model to m0.
simpeg.InvProblem will set Regularization.reference_model to m0.
simpeg.InvProblem will set Regularization.reference_model to m0.
simpeg.InvProblem will set Regularization.reference_model to m0.
simpeg.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.
***Done using same Solver, and solver_opts as the Simulation3DMagneticFluxDensity problem***
model has any nan: 0
============================ Inexact Gauss Newton ============================
# beta phi_d phi_m f |proj(x-g)-x| LS Comment
-----------------------------------------------------------------------------
x0 has any nan: 0
0 3.20e+00 1.42e+04 0.00e+00 1.42e+04 3.21e+03 0
1 3.20e+00 1.48e+03 5.42e+00 1.50e+03 4.61e+02 0
------------------------- STOP! -------------------------
1 : |fc-fOld| = 0.0000e+00 <= tolF*(1+|f0|) = 1.4191e+03
1 : |xc-x_last| = 6.2345e+00 <= tolX*(1+|x0|) = 1.3056e+01
0 : |proj(x-g)-x| = 4.6099e+02 <= tolG = 1.0000e-01
0 : |proj(x-g)-x| = 4.6099e+02 <= 1e3*eps = 1.0000e-02
0 : maxIter = 20 <= iter = 2
------------------------- DONE! -------------------------
Inversion Complete. Elapsed Time = 219.17 s
Plot the predicted and observed data#
fig, ax = plt.subplots(1, 3, figsize=(15, 5))
plot_data(dclean, ax=ax, color="C0", label="true")
plot_data(invProb.dpred, ax=ax, color="C1", label="predicted")
array([<Axes: title={'center': 'offset = 0.32m'}, xlabel='source location x (m)', ylabel='Secondary B-Field (T)'>,
<Axes: title={'center': 'offset = 0.71m'}, xlabel='source location x (m)', ylabel='Secondary B-Field (T)'>,
<Axes: title={'center': 'offset = 1.18m'}, xlabel='source location x (m)', ylabel='Secondary B-Field (T)'>],
dtype=object)
Plot the recovered model#
fig, ax = plt.subplots(1, 2, figsize=(12, 5))
# put both plots on the same colorbar
clim = np.r_[np.log(sigma_surface), np.log(sigma_deep)]
# recovered model
cb = plt.colorbar(
inversion_mesh.plot_image(mrec, ax=ax[0], clim=clim)[0],
ax=ax[0],
)
ax[0].set_title("recovered model")
cb.set_label(r"$\log(\sigma)$")
# true model
cb = plt.colorbar(
inversion_mesh.plot_image(m_true, ax=ax[1], clim=clim)[0],
ax=ax[1],
)
ax[1].set_title("true model")
cb.set_label(r"$\log(\sigma)$")
# # uncomment to plot the true interface
# x = np.linspace(-10, 10, 50)
# [a.plot(x, interface(x), 'k') for a in ax]
[a.set_xlim([-10, 10]) for a in ax]
[a.set_ylim([-2, 0]) for a in ax]
plt.tight_layout()
plt.show()
Print the version of SimPEG and dependencies#
Report()
Moving Forward#
If you have suggestions for improving this example, please create a pull request on the example in SimPEG
- You might try:
improving the discretization
changing beta
changing the noise model
playing with the regulariztion parameters
…
Total running time of the script: (3 minutes 57.973 seconds)
Estimated memory usage: 2935 MB