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
from pymatsolver import Pardiso as Solver

from SimPEG import (
    Mesh, Maps, Optimization,
    DataMisfit, Regularization, InvProblem,
    Inversion, Directives
from SimPEG.EM import FDEM
from SimPEG.Utils import mkvc


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., 8.]

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 = (
    )*(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.


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.]  # 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 = Mesh.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.hx.sum()/2., -mesh.hy.sum()/2., -mesh.hz[:npadz+ncz].sum()

print("the mesh has {} cells".format(mesh.nC))


the mesh has 34200 cells

Inversion Mesh

Here, we set up a 2D tensor mesh which we will represent the inversion model on

inversion_mesh = Mesh.TensorMesh([mesh.hx, mesh.hz[mesh.vectorCCz<=0]])
inversion_mesh.x0 = [-inversion_mesh.hx.sum()/2., -inversion_mesh.hy.sum()]


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 = Mesh.TensorMesh([mesh.hx, mesh.hz], x0=mesh.x0[[0, 2]])
active_inds = mesh2D.gridCC[:, 1] < 0  # active indices are below the surface

mapping = (
    Maps.Surject2Dto3D(mesh) *  # populates 3D space from a 2D model
    Maps.InjectActiveCells(mesh2D, active_inds, sigma_air) *  # adds air cells
    Maps.ExpMap(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.plotImage(m_true, ax=ax, grid=True)[0], ax=ax)
ax.set_title("true model")
ax.set_xlim([-10, 10])
ax.set_ylim([-2, 0])


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.] for sep in coil_separations])

# create our source list - one source per location
srcList = []
for x in src_locations:
    src_loc = np.r_[x, 0., src_z]
    rx_locs = src_loc - rx_offsets

    rx_real = FDEM.Rx.Point_bSecondary(
        locs=rx_locs, orientation=orientation, component='real'
    rx_imag = FDEM.Rx.Point_bSecondary(
        locs=rx_locs, orientation=orientation, component='imag'

    src = FDEM.Src.MagDipole(
        rxList=[rx_real, rx_imag], loc=src_loc, orientation=orientation,


# create the survey and problem objects for running the forward simulation
survey = FDEM.Survey(srcList)
prob = FDEM.Problem3D_b(mesh, sigmaMap=mapping, Solver=Solver)



Generate clean, synthetic data

t = time.time()
dclean = survey.dpred(m_true)
    "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):
            src_locations, data_real[i::len(coil_separations)],
            color=color, label="{} real".format(label)
            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].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)")

    return ax

ax = plot_data(dclean)


Done forward simulation. Elapsed time = 30.94 s

Set up data for inversion

We will invert the clean data, and assign a standard deviation of 0.03, and a floor of 1e-11.

survey.std = 0.03
survey.eps = 1e-11
survey.dobs = dclean

Set up the inversion

We create the data misfit, simple regularization (a Tikhonov-style regularization, SimPEG.Regularization.Simple) 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 = DataMisfit.l2_DataMisfit(survey)
reg = Regularization.Simple(inversion_mesh)
opt = Optimization.InexactGaussNewton(maxIterCG=10, remember="xc")
invProb = InvProblem.BaseInvProblem(dmisfit, reg, opt)

betaest = Directives.BetaEstimate_ByEig(beta0_ratio=0.25)
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 132.00

Run the inversion

We start from a half-space equal to the deep conductivity.

m0 = np.log(sigma_deep) * np.ones(inversion_mesh.nC)

t = time.time()
mrec = inv.run(m0)
print("\n Inversion Complete. Elapsed Time = {:1.2f} s".format(time.time() - t))


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  1.88e+00  3.66e+03  0.00e+00  3.66e+03    8.25e+02      0
   1  1.88e+00  3.95e+02  4.18e+01  4.74e+02    1.08e+02      0
------------------------- STOP! -------------------------
1 : |fc-fOld| = 0.0000e+00 <= tolF*(1+|f0|) = 3.6625e+02
1 : |xc-x_last| = 5.0176e+00 <= tolX*(1+|x0|) = 1.3056e+01
0 : |proj(x-g)-x|    = 1.0759e+02 <= tolG          = 1.0000e-01
0 : |proj(x-g)-x|    = 1.0759e+02 <= 1e3*eps       = 1.0000e-02
0 : maxIter   =      20    <= iter          =      2
------------------------- DONE! -------------------------

 Inversion Complete. Elapsed Time = 760.55 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")

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.plotImage(mrec, ax=ax[0], clim=clim)[0], ax=ax[0],
ax[0].set_title("recovered model")

# true model
cb = plt.colorbar(
    inversion_mesh.plotImage(m_true, ax=ax[1], clim=clim)[0], ax=ax[1],
ax[1].set_title("true model")

# # 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]


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: ( 13 minutes 13.444 seconds)

Gallery generated by Sphinx-Gallery