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

noiseEM = 1.  # absolute (per cent of primary signal)
noiseDC = 3.  # in per cent
lamEM, lamDC, lamDCEM = 300., 500., 500.  # regularization strength
verbose = False

First we create a synthetic model.

nlay = 3  # number of layers
thk = pg.Vector(nlay - 1, 15.0)  # 15m thickness each
res = pg.Vector(nlay, 200.0)  # 200 Ohmm
res[1] = 10.
res[2] = 50.
model = pg.cat(thk, res)  # paste together to one 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.

invEM = pg.core.Inversion(dataEM, fEM, transEM, True, True)
modelEM = pg.Vector(nlay * 2 - 1, 50.)
invEM.setModel(modelEM)
invEM.setAbsoluteError(noiseEM)
invEM.setLambda(lamEM)
invEM.setMarquardtScheme(0.9)
modelEM = invEM.run()
respEM = invEM.response()

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.