Sometimes we have additional information on the model that can be expressed in form of an equation. Let us consider a simple vertical electrical sounding (VES) that is inverted for a four-layer case with unknown layer thicknesses.

Assume we know the depth of the layer boundary, e.g. from a borehole, in this case between the third and fourth layer be $$z_2$$. So we can formulate this by an equation $$d_1+d_2+d_3=z_2$$. This equation can be added to the inverse problem in the style of Lagrangian, i.e. by an additional regularization term.

We use the LSQR inversion framework as used by Wagner et al. (2019) to constrain the sum of water, ice, air and rock fraction to be 1.

import numpy as np
import pygimli as pg
from pygimli.frameworks.lsqrinversion import LSQRInversion
from pygimli.physics.ves import VESModelling
from pygimli.viewer.mpl import drawModel1D


We set up a synthetic model of four layers and show the sounding curve

nlay = 4  # number of layers
lam = 200.  # (initial) regularization parameter
errPerc = 3.  # relative error of 3 percent
ab2 = np.logspace(-1, 2, 50)  # AB/2 distance (current electrodes)
mn2 = ab2 / 3.  # MN/2 distance (potential electrodes)
f = VESModelling(ab2=ab2, mn2=mn2, nLayers=nlay)
synres = [100., 500., 20., 800.]  # synthetic resistivity
synthk = [0.5, 3.5, 6.]  # synthetic thickness (nlay-th layer is infinite)
rhoa = f(synthk+synres)
rhoa = rhoa * (pg.randn(len(rhoa)) * errPerc / 100. + 1.)
pg.plt.loglog(rhoa, ab2, "x-")
pg.plt.grid(True)


Next, we set up an inversion instance with log transformation on data and model side and run the inversion.

inv = LSQRInversion(fop=f, verbose=True)
tLog = pg.trans.TransLog()  #
inv.dataTrans = tLog
inv.modelTrans = tLog
startModel = pg.cat(pg.Vector(nlay-1, 8), pg.Vector(nlay, pg.median(rhoa)))
inv.inv.setMarquardtScheme()
error = pg.Vector(len(rhoa), errPerc/100)
model1 = inv.run(rhoa, error, lam=1000, startModel=startModel)
print(model1)
print(inv.chi2(), inv.relrms(), pg.sum(inv.model[:nlay-1]))

fop: <pygimli.physics.ves.vesModelling.VESModelling object at 0x7f8031337100>
Data transformation: <pygimli.core._pygimli_.RTransLog object at 0x7f801dbae750>
Model transformation (cumulative):
0 <pygimli.core._pygimli_.RTransLogLU object at 0x7f80313a73a0>
1 <pygimli.core._pygimli_.RTransLogLU object at 0x7f803130c0a0>
min/max (data): 88.33/312
min/max (error): 3%/3%
min/max (start model): 8/140
--------------------------------------------------------------------------------
Jacobian MT:( 0 -- 7 ) / 7 ...
inv.iter 0 ... chi² =  190.10
--------------------------------------------------------------------------------
inv.iter 1 ... Jacobian MT:( 0 -- 7 ) / 7 ...
chi² =  174.06 (dPhi = 8.44%) lam: 1000.0
--------------------------------------------------------------------------------
inv.iter 2 ... Jacobian MT:( 0 -- 7 ) / 7 ...
chi² =  170.01 (dPhi = 2.33%) lam: 1000.0
--------------------------------------------------------------------------------
inv.iter 3 ... Jacobian MT:( 0 -- 7 ) / 7 ...
chi² =  166.37 (dPhi = 2.14%) lam: 1000.0
--------------------------------------------------------------------------------
inv.iter 4 ... Jacobian MT:( 0 -- 7 ) / 7 ...
chi² =  154.59 (dPhi = 7.08%) lam: 1000.0
--------------------------------------------------------------------------------
inv.iter 5 ... Jacobian MT:( 0 -- 7 ) / 7 ...
chi² =  133.15 (dPhi = 13.87%) lam: 1000.0
--------------------------------------------------------------------------------
inv.iter 6 ... Jacobian MT:( 0 -- 7 ) / 7 ...
chi² =   58.54 (dPhi = 56.03%) lam: 1000.0
--------------------------------------------------------------------------------
inv.iter 7 ... Jacobian MT:( 0 -- 7 ) / 7 ...
chi² =   13.69 (dPhi = 76.61%) lam: 1000.0
--------------------------------------------------------------------------------
inv.iter 8 ... Jacobian MT:( 0 -- 7 ) / 7 ...
chi² =    1.54 (dPhi = 88.74%) lam: 1000.0
--------------------------------------------------------------------------------
inv.iter 9 ... Jacobian MT:( 0 -- 7 ) / 7 ...
chi² =    0.85 (dPhi = 44.65%) lam: 1000.0

################################################################################
#                  Abort criterion reached: chi² <= 1 (0.85)                   #
################################################################################
7 [0.47756113568843633, 3.5984889865997918, 10.08337951145904, 100.02705937149561, 466.09368409548637, 33.40005059759289, 806.1065163478178]
0.8533334506972959 2.7797567795258313 14.159429633747267


To formulate the constraints, we need to set up a matrix for the left side and a vector for the right side of the equation. The layer thicknesses are the first values in the model vector, and setting 1 implements d1+d2+d3:

d1 d2 d3 r1 r2 r3 r4

G = [1 1 1 0 0 0 0] c = [z2] The constraints G * m = c are set by setParameterConstraints with a Lagrangian parameter that determines how well the equation needs to be fit.

G = pg.Matrix(rows=1, cols=len(startModel))
for i in range(3):
G.setVal(0, i, 1)

c = pg.Vector(1, pg.sum(synthk))
inv.setParameterConstraints(G, c, 100)
model2 = inv.run(rhoa, error, lam=1000, startModel=startModel)
print(model2)
print(inv.chi2(), inv.relrms(), pg.sum(inv.model[:nlay-1]))

fop: <pygimli.physics.ves.vesModelling.VESModelling object at 0x7f8031337100>
Data transformation: <pygimli.core._pygimli_.RTransLog object at 0x7f801dbae750>
Model transformation (cumulative):
0 <pygimli.core._pygimli_.RTransLogLU object at 0x7f803130fac0>
1 <pygimli.core._pygimli_.RTransLogLU object at 0x7f803130fac0>
min/max (data): 88.33/312
min/max (error): 3%/3%
min/max (start model): 8/140
--------------------------------------------------------------------------------
Jacobian MT:( 0 -- 7 ) / 7 ...
inv.iter 0 ... chi² =  190.10
--------------------------------------------------------------------------------
inv.iter 1 ... Jacobian MT:( 0 -- 7 ) / 7 ...
chi² =  167.58 (dPhi = 11.84%) lam: 1000.0
--------------------------------------------------------------------------------
inv.iter 2 ... Jacobian MT:( 0 -- 7 ) / 7 ...
chi² =  163.49 (dPhi = 2.44%) lam: 1000.0
--------------------------------------------------------------------------------
inv.iter 3 ... Jacobian MT:( 0 -- 7 ) / 7 ...
chi² =  145.41 (dPhi = 11.06%) lam: 1000.0
--------------------------------------------------------------------------------
inv.iter 4 ... Jacobian MT:( 0 -- 7 ) / 7 ...
chi² =   75.56 (dPhi = 48.03%) lam: 1000.0
--------------------------------------------------------------------------------
inv.iter 5 ... Jacobian MT:( 0 -- 7 ) / 7 ...
chi² =   42.94 (dPhi = 43.18%) lam: 1000.0
--------------------------------------------------------------------------------
inv.iter 6 ... Jacobian MT:( 0 -- 7 ) / 7 ...
chi² =   10.66 (dPhi = 75.17%) lam: 1000.0
--------------------------------------------------------------------------------
inv.iter 7 ... Jacobian MT:( 0 -- 7 ) / 7 ...
chi² =    1.14 (dPhi = 89.33%) lam: 1000.0
--------------------------------------------------------------------------------
inv.iter 8 ... Jacobian MT:( 0 -- 7 ) / 7 ...
chi² =    0.85 (dPhi = 25.19%) lam: 1000.0

################################################################################
#                  Abort criterion reached: chi² <= 1 (0.85)                   #
################################################################################
7 [0.4755776777228442, 3.7908826308108723, 5.737859354785321, 99.98585356290495, 462.2946464467531, 19.283180379893853, 775.9349467624263]
0.8513158854632327 2.7803159352775815 10.004319663319038


All three models are plotted together and are equivalently fitting the data. Due to the additional constraints, the model is much closer to the synthetic.

fig, ax = pg.plt.subplots()
drawModel1D(ax, synthk, synres, plot="semilogx", label="synth")
drawModel1D(ax, model=model1, label="unconstrained")
drawModel1D(ax, model=model2, label="constrained")
ax.set_ylim(15, 0)
ax.grid(True)
ax.legend()

<matplotlib.legend.Legend object at 0x7f801d065a50>


References#

Wagner, F.M., Mollaret, C., Günther, T., Kemna, A., Hauck, A. (2019):

Quantitative imaging of water, ice, and air in permafrost systems through petrophysical joint inversion of seismic refraction and electrical resistivity data. Geophys. J. Int. 219, 1866-1875. doi:10.1093/gji/ggz402.

Total running time of the script: (0 minutes 18.101 seconds)

Gallery generated by Sphinx-Gallery