Note
Click here to download the full example code
VES inversion for a blocky model#
This tutorial shows how an built-in forward operator is used for inversion. A DC 1D (VES) modelling is used to generate data, noisify and invert them.
We import numpy, matplotlib and the 1D plotting function
import numpy as np
import matplotlib.pyplot as plt
import pygimli as pg
from pygimli.viewer.mpl import drawModel1D
some definitions before (model, data and error)
initialize the forward modelling operator
other ways are by specifying a Data Container or am/an/bm/bn distances
the forward operator can be called by f.response(model) or simply f(model)
create some transformations used for inversion
transThk = pg.trans.TransLog() # log-transform ensures thk>0
transRho = pg.trans.TransLogLU(1, 1000) # lower and upper bound
transRhoa = pg.trans.TransLog() # log transformation for data
set model transformation for thickness and resistivity
f.region(0).setTransModel(transThk) # 0=thickness
f.region(1).setTransModel(transRho) # 1=resistivity
generate start model values from median app. resistivity & spread
set up inversion
inv = pg.core.Inversion(rhoa, f, transRhoa, True) # data vector, fop, verbose
# could also be set by inv.setTransData(transRhoa)
set error model, regularization strength and Marquardt scheme
inv.setRelativeError(errPerc / 100.0) # alternative: setAbsoluteError in Ohmm
inv.setLambda(lam) # (initial) regularization parameter
inv.setMarquardtScheme(0.9) # decrease lambda by factor 0.9
model = f.createStartVector() # creates from region start value
model[nlay] *= 1.5 # change default model by changing 2nd layer resistivity
inv.setModel(model) #
run actual inversion and extract resistivity and thickness
rrms=3.21%, chi^2=1.136
show estimated&synthetic models and data with model response in 2 subplots
fig, ax = plt.subplots(ncols=2, figsize=(8, 6)) # two-column figure
drawModel1D(ax[0], synthk, synres, plot='semilogx', color='r')
drawModel1D(ax[0], thk, res, color='b')
ax[0].grid(True, which='both')
ax[0].set_ylabel('z (m)')
ax[0].set_xlabel(r'$\rho$ ($\Omega$m)')
ax[1].loglog(rhoa, ab2, 'rx-', label='data') # sounding curve
ax[1].loglog(inv.response(), ab2, 'b-', label='response')
ax[1].set_ylim((max(ab2), min(ab2))) # downwards according to penetration
ax[1].grid(True, which='both')
ax[1].set_xlabel(r'$\rho_a$ ($\Omega$m)')
ax[1].set_ylabel('AB/2 (m)')
ax[1].legend(loc='best')
plt.show()
