Note
Go to the end to download the full example code.
Heagy et al., 2017 1D FDEM and TDEM inversions#
Here, we perform a 1D inversion using both the frequency and time domain codes. The forward simulations are conducted on a cylindrically symmetric mesh. The source is a point magnetic dipole source.
This example is used in the paper
Lindsey J. Heagy, Rowan Cockett, Seogi Kang, Gudni K. Rosenkjaer, Douglas W. Oldenburg, A framework for simulation and inversion in electromagnetics, Computers & Geosciences, Volume 107, 2017, Pages 1-19, ISSN 0098-3004, http://dx.doi.org/10.1016/j.cageo.2017.06.018.
This example is on figshare: https://doi.org/10.6084/m9.figshare.5035175
This example was updated for SimPEG 0.14.0 on January 31st, 2020 by Joseph Capriotti
min skin depth = 158.11388300841895 max skin depth = 500.0
max x 1267.687908603637 min z -1242.6879086036365 max z 1242.687908603637
/usr/share/miniconda/envs/simpeg-test/lib/python3.11/site-packages/scipy/sparse/_construct.py:543: FutureWarning:
Input has data type int64, but the output has been cast to float64. In the future, the output data type will match the input. To avoid this warning, set the `dtype` parameter to `None` to have the output dtype match the input, or set it to the desired output data type.
Note: In Python 3.11, this warning can be generated by a call of scipy.sparse.diags(), but the code indicated in the warning message will refer to an internal call of scipy.sparse.diags_array(). If that happens, check your code for the use of diags().
Running inversion with SimPEG v0.25.1.dev9+g471344c9a
============================ Inexact Gauss Newton ============================
# beta phi_d phi_m f |proj(x-g)-x| LS Comment
-----------------------------------------------------------------------------
0 3.89e+00 1.67e+03 0.00e+00 1.67e+03
1 3.89e+00 4.06e+02 6.60e+01 6.63e+02 1.19e+03 0
2 3.89e+00 1.10e+02 7.16e+01 3.89e+02 1.07e+03 0
3 3.89e+00 1.09e+02 6.57e+01 3.64e+02 1.96e+02 0 Skip BFGS
4 9.72e-01 4.82e+01 9.58e+01 1.41e+02 1.31e+02 0
5 9.72e-01 3.19e+01 9.42e+01 1.23e+02 2.31e+02 0
6 9.72e-01 2.81e+01 9.41e+01 1.20e+02 5.80e+01 0 Skip BFGS
7 2.43e-01 2.10e+01 1.09e+02 4.75e+01 9.63e+01 0
8 2.43e-01 1.63e+01 1.07e+02 4.23e+01 1.26e+02 0
9 2.43e-01 1.50e+01 1.08e+02 4.12e+01 5.23e+01 1 Skip BFGS
10 6.07e-02 1.26e+01 1.19e+02 1.98e+01 4.15e+01 0
11 6.07e-02 1.23e+01 1.17e+02 1.94e+01 3.46e+01 0
12 6.07e-02 1.23e+01 1.17e+02 1.94e+01 4.22e+00 1 Skip BFGS
13 1.52e-02 1.22e+01 1.19e+02 1.40e+01 7.37e+00 2
14 1.52e-02 1.21e+01 1.20e+02 1.40e+01 1.42e+01 2 Skip BFGS
15 1.52e-02 1.20e+01 1.20e+02 1.38e+01 1.73e+01 0
16 3.80e-03 1.20e+01 1.21e+02 1.24e+01 1.95e+00 3
17 3.80e-03 1.19e+01 1.23e+02 1.24e+01 4.85e+00 2 Skip BFGS
18 3.80e-03 1.19e+01 1.24e+02 1.23e+01 9.89e+00 0
19 9.49e-04 1.19e+01 1.25e+02 1.20e+01 8.26e-01 3 Skip BFGS
20 9.49e-04 1.18e+01 1.27e+02 1.20e+01 5.34e+00 2 Skip BFGS
------------------------- STOP! -------------------------
1 : |fc-fOld| = 1.5570e-02 <= tolF*(1+|f0|) = 1.6677e+02
1 : |xc-x_last| = 1.6095e-01 <= tolX*(1+|x0|) = 2.4026e+00
0 : |proj(x-g)-x| = 8.6662e+00 <= tolG = 1.0000e-01
0 : |proj(x-g)-x| = 8.6662e+00 <= 1e3*eps = 1.0000e-02
1 : maxIter = 20 <= iter = 20
------------------------- DONE! -------------------------
min diffusion distance 114.18394344131536 max diffusion distance 510.6461189475449
Running inversion with SimPEG v0.25.1.dev9+g471344c9a
============================ Inexact Gauss Newton ============================
# beta phi_d phi_m f |proj(x-g)-x| LS Comment
-----------------------------------------------------------------------------
0 3.24e+00 3.24e+03 0.00e+00 3.24e+03
1 3.24e+00 2.12e+03 1.55e+02 2.62e+03 1.42e+03 0
2 3.24e+00 1.70e+02 1.19e+02 5.55e+02 4.15e+03 0
3 3.24e+00 1.41e+02 1.05e+02 4.81e+02 5.55e+02 0 Skip BFGS
4 8.09e-01 6.80e+01 1.22e+02 1.67e+02 4.86e+02 0 Skip BFGS
5 8.09e-01 2.45e+01 1.13e+02 1.16e+02 5.08e+02 0
6 8.09e-01 1.88e+01 1.10e+02 1.08e+02 6.64e+01 1
7 2.02e-01 1.68e+01 1.09e+02 3.89e+01 7.03e+01 2 Skip BFGS
8 2.02e-01 1.64e+01 1.09e+02 3.85e+01 9.19e+01 1 Skip BFGS
9 2.02e-01 1.56e+01 1.13e+02 3.85e+01 1.39e+02 0 Skip BFGS
10 5.06e-02 1.21e+01 1.11e+02 1.77e+01 1.54e+02 0
11 5.06e-02 1.19e+01 1.09e+02 1.74e+01 1.75e+01 2
12 5.06e-02 1.08e+01 1.11e+02 1.64e+01 7.18e+01 0
13 1.26e-02 1.06e+01 1.12e+02 1.20e+01 2.74e+01 3 Skip BFGS
14 1.26e-02 1.05e+01 1.14e+02 1.20e+01 4.42e+01 2
15 1.26e-02 1.02e+01 1.13e+02 1.16e+01 4.80e+01 0
16 3.16e-03 9.98e+00 1.17e+02 1.03e+01 4.70e+00 4
------------------------- STOP! -------------------------
1 : |fc-fOld| = 2.2767e-01 <= tolF*(1+|f0|) = 3.2376e+02
1 : |xc-x_last| = 2.4002e-01 <= tolX*(1+|x0|) = 2.4026e+00
0 : |proj(x-g)-x| = 4.6974e+00 <= tolG = 1.0000e-01
0 : |proj(x-g)-x| = 4.6974e+00 <= 1e3*eps = 1.0000e-02
0 : maxIter = 20 <= iter = 16
------------------------- DONE! -------------------------
import discretize
from simpeg import (
maps,
utils,
data_misfit,
regularization,
optimization,
inversion,
inverse_problem,
directives,
)
import numpy as np
from simpeg.electromagnetics import frequency_domain as FDEM, time_domain as TDEM, mu_0
import matplotlib.pyplot as plt
import matplotlib
def run(plotIt=True, saveFig=False):
# Set up cylindrically symmeric mesh
cs, ncx, ncz, npad = 10.0, 15, 25, 13 # padded cyl mesh
hx = [(cs, ncx), (cs, npad, 1.3)]
hz = [(cs, npad, -1.3), (cs, ncz), (cs, npad, 1.3)]
mesh = discretize.CylindricalMesh([hx, 1, hz], "00C")
# Conductivity model
layerz = np.r_[-200.0, -100.0]
layer = (mesh.cell_centers_z >= layerz[0]) & (mesh.cell_centers_z <= layerz[1])
active = mesh.cell_centers_z < 0.0
sig_half = 1e-2 # Half-space conductivity
sig_air = 1e-8 # Air conductivity
sig_layer = 5e-2 # Layer conductivity
sigma = np.ones(mesh.shape_cells[2]) * sig_air
sigma[active] = sig_half
sigma[layer] = sig_layer
# Mapping
actMap = maps.InjectActiveCells(mesh, active, np.log(1e-8), nC=mesh.shape_cells[2])
mapping = maps.ExpMap(mesh) * maps.SurjectVertical1D(mesh) * actMap
mtrue = np.log(sigma[active])
# ----- FDEM problem & survey ----- #
rxlocs = utils.ndgrid([np.r_[50.0], np.r_[0], np.r_[0.0]])
bzr = FDEM.Rx.PointMagneticFluxDensitySecondary(rxlocs, "z", "real")
bzi = FDEM.Rx.PointMagneticFluxDensitySecondary(rxlocs, "z", "imag")
freqs = np.logspace(2, 3, 5)
srcLoc = np.array([0.0, 0.0, 0.0])
print(
"min skin depth = ",
500.0 / np.sqrt(freqs.max() * sig_half),
"max skin depth = ",
500.0 / np.sqrt(freqs.min() * sig_half),
)
print(
"max x ",
mesh.cell_centers_x.max(),
"min z ",
mesh.cell_centers_z.min(),
"max z ",
mesh.cell_centers_z.max(),
)
source_list = [
FDEM.Src.MagDipole([bzr, bzi], freq, srcLoc, orientation="Z") for freq in freqs
]
surveyFD = FDEM.Survey(source_list)
prbFD = FDEM.Simulation3DMagneticFluxDensity(
mesh, survey=surveyFD, sigmaMap=mapping
)
rel_err = 0.03
dataFD = prbFD.make_synthetic_data(mtrue, relative_error=rel_err, add_noise=True)
dataFD.noise_floor = np.linalg.norm(dataFD.dclean) * 1e-5
# FDEM inversion
np.random.seed(1)
dmisfit = data_misfit.L2DataMisfit(simulation=prbFD, data=dataFD)
regMesh = discretize.TensorMesh([mesh.h[2][mapping.maps[-1].active_cells]])
reg = regularization.WeightedLeastSquares(regMesh)
opt = optimization.InexactGaussNewton(cg_maxiter=10)
invProb = inverse_problem.BaseInvProblem(dmisfit, reg, opt)
# Inversion Directives
beta = directives.BetaSchedule(coolingFactor=4, coolingRate=3)
betaest = directives.BetaEstimate_ByEig(beta0_ratio=1.0, random_seed=518936)
target = directives.TargetMisfit()
directiveList = [beta, betaest, target]
inv = inversion.BaseInversion(invProb, directiveList=directiveList)
m0 = np.log(np.ones(mtrue.size) * sig_half)
reg.alpha_s = 5e-1
reg.alpha_x = 1.0
prbFD.counter = opt.counter = utils.Counter()
opt.remember("xc")
moptFD = inv.run(m0)
# TDEM problem
times = np.logspace(-4, np.log10(2e-3), 10)
print(
"min diffusion distance ",
1.28 * np.sqrt(times.min() / (sig_half * mu_0)),
"max diffusion distance ",
1.28 * np.sqrt(times.max() / (sig_half * mu_0)),
)
rx = TDEM.Rx.PointMagneticFluxDensity(rxlocs, times, "z")
src = TDEM.Src.MagDipole(
[rx],
waveform=TDEM.Src.StepOffWaveform(),
location=srcLoc, # same src location as FDEM problem
)
surveyTD = TDEM.Survey([src])
prbTD = TDEM.Simulation3DMagneticFluxDensity(
mesh, survey=surveyTD, sigmaMap=mapping
)
prbTD.time_steps = [(5e-5, 10), (1e-4, 10), (5e-4, 10)]
rel_err = 0.03
dataTD = prbTD.make_synthetic_data(mtrue, relative_error=rel_err, add_noise=True)
dataTD.noise_floor = np.linalg.norm(dataTD.dclean) * 1e-5
# TDEM inversion
dmisfit = data_misfit.L2DataMisfit(simulation=prbTD, data=dataTD)
regMesh = discretize.TensorMesh([mesh.h[2][mapping.maps[-1].active_cells]])
reg = regularization.WeightedLeastSquares(regMesh)
opt = optimization.InexactGaussNewton(cg_maxiter=10)
invProb = inverse_problem.BaseInvProblem(dmisfit, reg, opt)
# directives
beta = directives.BetaSchedule(coolingFactor=4, coolingRate=3)
betaest = directives.BetaEstimate_ByEig(beta0_ratio=1.0, random_seed=518936)
target = directives.TargetMisfit()
directiveList = [beta, betaest, target]
inv = inversion.BaseInversion(invProb, directiveList=directiveList)
m0 = np.log(np.ones(mtrue.size) * sig_half)
reg.alpha_s = 5e-1
reg.alpha_x = 1.0
prbTD.counter = opt.counter = utils.Counter()
opt.remember("xc")
moptTD = inv.run(m0)
# Plot the results
if plotIt:
plt.figure(figsize=(10, 8))
ax0 = plt.subplot2grid((2, 2), (0, 0), rowspan=2)
ax1 = plt.subplot2grid((2, 2), (0, 1))
ax2 = plt.subplot2grid((2, 2), (1, 1))
fs = 13 # fontsize
matplotlib.rcParams["font.size"] = fs
# Plot the model
# z_true = np.repeat(mesh.cell_centers_z[active][1:], 2, axis=0)
# z_true = np.r_[mesh.cell_centers_z[active][0], z_true, mesh.cell_centers_z[active][-1]]
activeN = mesh.nodes_z <= 0.0 + cs / 2.0
z_true = np.repeat(mesh.nodes_z[activeN][1:-1], 2, axis=0)
z_true = np.r_[mesh.nodes_z[activeN][0], z_true, mesh.nodes_z[activeN][-1]]
sigma_true = np.repeat(sigma[active], 2, axis=0)
ax0.semilogx(sigma_true, z_true, "k-", lw=2, label="True")
ax0.semilogx(
np.exp(moptFD),
mesh.cell_centers_z[active],
"bo",
ms=6,
markeredgecolor="k",
markeredgewidth=0.5,
label="FDEM",
)
ax0.semilogx(
np.exp(moptTD),
mesh.cell_centers_z[active],
"r*",
ms=10,
markeredgecolor="k",
markeredgewidth=0.5,
label="TDEM",
)
ax0.set_ylim(-700, 0)
ax0.set_xlim(5e-3, 1e-1)
ax0.set_xlabel("Conductivity (S/m)", fontsize=fs)
ax0.set_ylabel("Depth (m)", fontsize=fs)
ax0.grid(which="both", color="k", alpha=0.5, linestyle="-", linewidth=0.2)
ax0.legend(fontsize=fs, loc=4)
# plot the data misfits - negative b/c we choose positive to be in the
# direction of primary
ax1.plot(freqs, -dataFD.dobs[::2], "k-", lw=2, label="Obs (real)")
ax1.plot(freqs, -dataFD.dobs[1::2], "k--", lw=2, label="Obs (imag)")
dpredFD = prbFD.dpred(moptTD)
ax1.loglog(
freqs,
-dpredFD[::2],
"bo",
ms=6,
markeredgecolor="k",
markeredgewidth=0.5,
label="Pred (real)",
)
ax1.loglog(
freqs, -dpredFD[1::2], "b+", ms=10, markeredgewidth=2.0, label="Pred (imag)"
)
ax2.loglog(times, dataTD.dobs, "k-", lw=2, label="Obs")
ax2.loglog(
times,
prbTD.dpred(moptTD),
"r*",
ms=10,
markeredgecolor="k",
markeredgewidth=0.5,
label="Pred",
)
ax2.set_xlim(times.min() - 1e-5, times.max() + 1e-4)
# Labels, gridlines, etc
ax2.grid(which="both", alpha=0.5, linestyle="-", linewidth=0.2)
ax1.grid(which="both", alpha=0.5, linestyle="-", linewidth=0.2)
ax1.set_xlabel("Frequency (Hz)", fontsize=fs)
ax1.set_ylabel("Vertical magnetic field (-T)", fontsize=fs)
ax2.set_xlabel("Time (s)", fontsize=fs)
ax2.set_ylabel("Vertical magnetic field (T)", fontsize=fs)
ax2.legend(fontsize=fs, loc=3)
ax1.legend(fontsize=fs, loc=3)
ax1.set_xlim(freqs.max() + 1e2, freqs.min() - 1e1)
ax0.set_title("(a) Recovered Models", fontsize=fs)
ax1.set_title("(b) FDEM observed vs. predicted", fontsize=fs)
ax2.set_title("(c) TDEM observed vs. predicted", fontsize=fs)
plt.tight_layout(pad=1.5)
if saveFig is True:
plt.savefig("example1.png", dpi=600)
if __name__ == "__main__":
run(plotIt=True, saveFig=True)
plt.show()
Total running time of the script: (0 minutes 54.887 seconds)
Estimated memory usage: 321 MB