DC-EM joint inversion#
Joint inversion is an important method to improve resolution properties by combining different methods. In the easiest case, the methods have the same subsurface parameter, for example direct current (DC) and electromagnetic (EM) measurements. In this example, we illustrate how two modelling operators can be combined by the JointInversion framework, using a vertical electric sounding (VES) and electromagnetic frequency sounding (FDEM).
A similar case has been documented by [Günther, 2013].
# We import the numpy and matplotlib import numpy as np import matplotlib.pyplot as plt # Next we import pyGIMLi and the modelling operators for block models import pygimli as pg from pygimli.physics.em import HEMmodelling from pygimli.physics.ves import VESModelling from pygimli.frameworks import JointModelling # For block models we need the Marquardt-Levenberg inversion scheme from pygimli.frameworks import MarquardtInversion from pygimli.viewer.mpl import drawModel1D
First we create a synthetic model and error models to be used later.
Part 1: Electromagnetic sounding#
We set forward operator and generate synthetic data with noise.
Result is inphase and outphase secondary fields divided by the primary field as a function of frequency. See inversion result & data fit figure below.
We first set up the independent EM inversion and run the model.
Part 2: Vertical Electric Sounding#
We set up the (DC) forward operator and generate synthetic data plus noise.
We set up the independent DC inversion and let it run.
Part 3: Joint inversion#
We create a the joint forward operator using the Joint inversion framework.
fDCEM = JointModelling([fDC, fEM]) fDCEM.setData([dataDC, dataEM]) # just for sizes!
Inversion is just as for the single inversions. The data vector is created by concatenating both data vectors. This is also done for the relative error.
The final output of the inversion is plotted for every method. Most-important measure is the chi-squared misfit that should be close to 1.
for inv in [invEM, invDC, invDCEM]: inv.echoStatus() print([invEM.chi2(), invDC.chi2(), invDCEM.chi2()]) # chi-square values
[1.0012335898454139, 1.1101824803379556, 1.0891389114801069]
We finally plot the inverted models along with data and model responses.
fig, (ax1, ax2, ax3) = plt.subplots(figsize=(10, 5), ncols=3) drawModel1D(ax1, synThk, synRes, plot='semilogx', color='C0', label="synth") drawModel1D(ax1, model=modelEM, color='C1', label="EM") drawModel1D(ax1, model=modelDC, color='C2', label="DC") drawModel1D(ax1, model=modelDCEM, color='C3', label="DC-EM") ax1.legend() ax1.set_xlim((10., 1000.)) ax1.set_ylim((40., 0.)) ax1.grid(which='both') ax2.semilogy(dataEM[0:nf], freq, 'x', color="C0", label='syn IP') ax2.semilogy(dataEM[nf:nf*2], freq, 'o', color="C0", label='syn OP') ax2.semilogy(invEM.response[0:nf], freq, '--', color="C1", label='EM') ax2.semilogy(invEM.response[nf:nf*2], freq, '--', color="C1") ax2.semilogy(invDCEM.response[na:na+nf], freq, ':', color="C3", label='DCEM') ax2.semilogy(invDCEM.response[na+nf:na+nf*2], freq, '2:', color="C3") 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.loglog(dataDC, ab2, 'x-', label='syn', color="C0") ax3.loglog(invDC.response, ab2, '-', label='DC', color="C2") ax3.loglog(invDCEM.response[0:na], ab2, '-', label='DCEM', color="C3") 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()
<matplotlib.legend.Legend object at 0x2e30c59f0>
All three inversions are able to reveal the subsurface structures. EM fails to describe the first layer-resistivity and also its thickness, for which DC does a better job. Both are similarly away from the synthetic model regarding the resistivity and upper depth of the third layer. EM can better resolve the good conductor at depth as expected. The joint inversion result combines the resolution properties of both methods and yields a result that is very close to the synthetic.
# 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.