Note
Click here to download the full example code
EM: TDEM: 1D: InversionΒΆ
Here we will create and run a TDEM 1D inversion.

Out:
SimPEG.Survey assigned new std of 5.00%
SimPEG.InvProblem will set Regularization.mref to m0.
SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.
***Done using same Solver and solverOpts as the 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 2.81e+03 2.60e+03 0.00e+00 2.60e+03 3.13e+03 0
1 2.81e+03 2.61e+02 7.46e-02 4.71e+02 4.03e+02 0
2 5.63e+02 1.07e+02 1.08e-01 1.68e+02 2.85e+02 0 Skip BFGS
3 5.63e+02 1.96e+01 1.68e-01 1.14e+02 1.28e+01 0 Skip BFGS
4 1.13e+02 1.86e+01 1.69e-01 3.77e+01 7.46e+01 0 Skip BFGS
5 1.13e+02 9.75e+00 1.97e-01 3.19e+01 2.95e+00 0
------------------------- STOP! -------------------------
1 : |fc-fOld| = 5.7587e+00 <= tolF*(1+|f0|) = 2.6055e+02
1 : |xc-x_last| = 2.0393e-01 <= tolX*(1+|x0|) = 3.0149e+00
0 : |proj(x-g)-x| = 2.9460e+00 <= tolG = 1.0000e-01
0 : |proj(x-g)-x| = 2.9460e+00 <= 1e3*eps = 1.0000e-02
1 : maxIter = 5 <= iter = 5
------------------------- DONE! -------------------------
/Users/lindseyjh/git/simpeg/simpeg/examples/09-tdem/plot_inversion_1D.py:72: UserWarning: Data has no positive values, and therefore cannot be log-scaled.
ax[0].loglog(rx.times, survey.dtrue, 'b.-')
/Users/lindseyjh/git/simpeg/simpeg/examples/09-tdem/plot_inversion_1D.py:73: UserWarning: Data has no positive values, and therefore cannot be log-scaled.
ax[0].loglog(rx.times, survey.dobs, 'r.-')
/Users/lindseyjh/git/simpeg/simpeg/examples/09-tdem/plot_inversion_1D.py:92: UserWarning: Matplotlib is currently using agg, which is a non-GUI backend, so cannot show the figure.
plt.show()
import numpy as np
from SimPEG import (Mesh, Maps, SolverLU, DataMisfit, Regularization,
Optimization, InvProblem, Inversion, Directives, Utils)
import SimPEG.EM as EM
import matplotlib.pyplot as plt
def run(plotIt=True):
cs, ncx, ncz, npad = 5., 25, 15, 15
hx = [(cs, ncx), (cs, npad, 1.3)]
hz = [(cs, npad, -1.3), (cs, ncz), (cs, npad, 1.3)]
mesh = Mesh.CylMesh([hx, 1, hz], '00C')
active = mesh.vectorCCz < 0.
layer = (mesh.vectorCCz < 0.) & (mesh.vectorCCz >= -100.)
actMap = Maps.InjectActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * actMap
sig_half = 2e-3
sig_air = 1e-8
sig_layer = 1e-3
sigma = np.ones(mesh.nCz)*sig_air
sigma[active] = sig_half
sigma[layer] = sig_layer
mtrue = np.log(sigma[active])
rxOffset = 1e-3
rx = EM.TDEM.Rx.Point_dbdt(
np.array([[rxOffset, 0., 30]]),
np.logspace(-5, -3, 31),
'z'
)
src = EM.TDEM.Src.MagDipole([rx], loc=np.array([0., 0., 80]))
survey = EM.TDEM.Survey([src])
prb = EM.TDEM.Problem3D_e(mesh, sigmaMap=mapping)
prb.Solver = SolverLU
prb.timeSteps = [(1e-06, 20), (1e-05, 20), (0.0001, 20)]
prb.pair(survey)
# create observed data
std = 0.05
survey.dobs = survey.makeSyntheticData(mtrue, std)
survey.std = std
survey.eps = 1e-5*np.linalg.norm(survey.dobs)
dmisfit = DataMisfit.l2_DataMisfit(survey)
regMesh = Mesh.TensorMesh([mesh.hz[mapping.maps[-1].indActive]])
reg = Regularization.Tikhonov(regMesh, alpha_s=1e-2, alpha_x=1.)
opt = Optimization.InexactGaussNewton(maxIter=5, LSshorten=0.5)
invProb = InvProblem.BaseInvProblem(dmisfit, reg, opt)
# Create an inversion object
beta = Directives.BetaSchedule(coolingFactor=5, coolingRate=2)
betaest = Directives.BetaEstimate_ByEig(beta0_ratio=1e0)
inv = Inversion.BaseInversion(invProb, directiveList=[beta, betaest])
m0 = np.log(np.ones(mtrue.size)*sig_half)
prb.counter = opt.counter = Utils.Counter()
opt.remember('xc')
mopt = inv.run(m0)
if plotIt:
fig, ax = plt.subplots(1, 2, figsize=(10, 6))
ax[0].loglog(rx.times, survey.dtrue, 'b.-')
ax[0].loglog(rx.times, survey.dobs, 'r.-')
ax[0].legend(('Noisefree', '$d^{obs}$'), fontsize=16)
ax[0].set_xlabel('Time (s)', fontsize=14)
ax[0].set_ylabel('$B_z$ (T)', fontsize=16)
ax[0].set_xlabel('Time (s)', fontsize=14)
ax[0].grid(color='k', alpha=0.5, linestyle='dashed', linewidth=0.5)
plt.semilogx(sigma[active], mesh.vectorCCz[active])
plt.semilogx(np.exp(mopt), mesh.vectorCCz[active])
ax[1].set_ylim(-600, 0)
ax[1].set_xlim(1e-4, 1e-2)
ax[1].set_xlabel('Conductivity (S/m)', fontsize=14)
ax[1].set_ylabel('Depth (m)', fontsize=14)
ax[1].grid(color='k', alpha=0.5, linestyle='dashed', linewidth=0.5)
plt.legend(['$\sigma_{true}$', '$\sigma_{pred}$'])
if __name__ == '__main__':
run()
plt.show()
Total running time of the script: ( 0 minutes 11.605 seconds)