Note
Click here to download the full example code
DC-EM Joint inversion¶
This is an old script from early pyGIMLi jointly inverting direct current (DC) and electromagnetic (EM) soundings on the modelling abstraction level. Note that this is not recommended as a basis for programming, because there is a dedicated framework for classical joint inversion. However, it explains what happens under the hood in the much simpler script that follows.
The case has been documented by [Gunther13].
import numpy as np
import matplotlib.pyplot as plt
import pygimli as pg
from pygimli.viewer.mpl import drawModel1D
First, we define a modelling class that calls two other classes and pastes their results to one vector.
class DCEM1dModelling(pg.core.ModellingBase):
"""Modelling jointing DC and EM 1Dforward operators."""
def __init__(self, nlay, ab2, mn2, freq, coilspacing, verbose=False):
"""Init number of layers, AB/2, MN/2, frequencies & coil spacing."""
pg.core.ModellingBase.__init__(self, verbose)
self.nlay_ = nlay
self.fDC_ = pg.core.DC1dModelling(nlay, ab2, mn2, verbose)
self.fEM_ = pg.core.FDEM1dModelling(nlay, freq, coilspacing, verbose)
self.mesh_ = pg.meshtools.createMesh1DBlock(nlay)
self.setMesh(self.mesh_)
def response(self, model):
"""Return concatenated response of DC and EM FOPs."""
return pg.cat(self.fDC_(model), self.fEM_(model))
The actual script starts here. There are some options to play with
First we create a synthetic model.
We first set up EM forward operator and generate synthetic data with noise
coilspacing = 50.
nf = 10
freq = pg.Vector(nf, 110.)
for i in range(nf-1):
freq[i+1] = freq[i] * 2.
fEM = pg.core.FDEM1dModelling(nlay, freq, coilspacing)
dataEM = fEM(model)
dataEM += pg.randn(len(dataEM), seed=1234) * noiseEM
We define model transformations: logarithms and log with upper+lower bounds
transRhoa = pg.trans.TransLog()
transThk = pg.trans.TransLog()
transRes = pg.trans.TransLogLU(1., 1000.)
transEM = pg.trans.Trans()
fEM.region(0).setTransModel(transThk)
fEM.region(1).setTransModel(transRes)
We set up the independent EM inversion and run the model.
Next we set up the DC forward operator and generate synthetic data with noise
ab2 = pg.Vector(20, 3.)
na = len(ab2)
mn2 = pg.Vector(na, 1.0)
for i in range(na-1):
ab2[i+1] = ab2[i] * 1.3
fDC = pg.core.DC1dModelling(nlay, ab2, mn2)
dataDC = fDC(model)
dataDC *= 1. + pg.randn(len(dataDC), seed=1234) * noiseDC / 100.
fDC.region(0).setTransModel(transThk)
fDC.region(1).setTransModel(transRes)
# We set up the independent DC inversion and let it run.
invDC = pg.core.Inversion(dataDC, fDC, transRhoa, verbose)
modelDC = pg.Vector(nlay*2-1, 20.)
invDC.setModel(modelDC)
invDC.setRelativeError(noiseDC/100.)
invDC.setLambda(lamDC)
invDC.setMarquardtScheme(0.9)
modelDC = invDC.run()
respDC = invDC.response()
Next we create a the joint forward operator (see class above).
fDCEM = DCEM1dModelling(nlay, ab2, mn2, freq, coilspacing)
fDCEM.region(0).setTransModel(transThk)
fDCEM.region(1).setTransModel(transRes)
We setup the joint inversion combining, transformations, data and errors.
transData = pg.trans.TransCumulative()
transData.add(transRhoa, na)
transData.add(transEM, nf*2)
invDCEM = pg.core.Inversion(pg.cat(dataDC, dataEM), fDCEM, transData, verbose)
modelDCEM = pg.Vector(nlay * 2 - 1, 20.)
invDCEM.setModel(modelDCEM)
err = pg.cat(dataDC * noiseDC / 100., pg.Vector(len(dataEM), noiseEM))
invDCEM.setAbsoluteError(err)
invDCEM.setLambda(lamDCEM)
invDCEM.setMarquardtScheme(0.9)
modelDCEM = invDCEM.run()
respDCEM = invDCEM.response()
The results of the inversion are plotted for comparison.
for inv in [invEM, invDC, invDCEM]:
inv.echoStatus()
print([invEM.chi2(), invDC.chi2(), invDCEM.chi2()]) # chi-square values
[1.0574072316294534, 1.1414069871330406, 1.1706466870746692]
%% We finally plot the results
fig = plt.figure(1, figsize=(10, 5))
ax1 = fig.add_subplot(131)
drawModel1D(ax1, thk, res, plot='semilogx', color='blue')
drawModel1D(ax1, modelEM[0:nlay-1], modelEM[nlay-1:nlay*2-1], color='green')
drawModel1D(ax1, modelDC[0:nlay-1], modelDC[nlay-1:nlay*2-1], color='cyan')
drawModel1D(ax1, modelDCEM[0:nlay-1], modelDCEM[nlay-1:nlay*2-1],
color='red')
ax1.legend(('syn', 'EM', 'DC', 'JI'))
ax1.set_xlim((10., 1000.))
ax1.set_ylim((40., 0.))
ax1.grid(which='both')
ax2 = fig.add_subplot(132)
ax2.semilogy(dataEM[0:nf], freq, 'bx', label='syn IP')
ax2.semilogy(dataEM[nf:nf*2], freq, 'bo', label='syn OP')
ax2.semilogy(respEM[0:nf], freq, 'g--', label='EM')
ax2.semilogy(respEM[nf:nf*2], freq, 'g--')
ax2.semilogy(respDCEM[na:na+nf], freq, 'r:', label='DCEM')
ax2.semilogy(respDCEM[na+nf:na+nf*2], freq, 'r:')
ax2.set_ylim((min(freq), max(freq)))
ax2.set_xlabel("IP/OP in %")
ax2.set_ylabel("$f$ in Hz")
ax2.yaxis.set_label_position("right")
ax2.grid(which='both')
ax2.legend(loc="best")
ax3 = fig.add_subplot(133)
ax3.loglog(dataDC, ab2, 'bx-', label='syn')
ax3.loglog(respDC, ab2, 'c-', label='DC')
ax3.loglog(respDCEM[0:na], ab2, 'r:', label='DCEM')
# ax3.axis('tight')
ax3.set_ylim((max(ab2), min(ab2)))
ax3.grid(which='both')
ax3.set_xlabel(r"$\rho_a$ in $\Omega$m")
ax3.set_ylabel("AB/2 in m")
ax3.yaxis.set_ticks_position("right")
ax3.yaxis.set_label_position("right")
ax3.legend(loc="best")
pg.wait()

# Günther, T. (2013): On Inversion of Frequency Domain Electromagnetic Data in
# Salt Water Problems - Sensitivity and Resolution. Ext. Abstr., 19th European
# Meeting of Environmental and Engineering Geophysics, Bochum, Germany.