Inversion with structural constraints#

Any inversion suffers from ambiguity and missing resolution at depth. Structural information can provide valuable information about geological boundaries in the subsurface. Such information can come from wave methods like seismics (Tanner et al. 2020), ground-penetrating radar (Doetsch et al. 2012, Jiang et al. 2020) or even piece-wise from boreholes (Wunderlich et al. 2018). This can be done in both 2D (Jiang et al. 2020) using lines or in 3D using facets (Doetsch et al. 2012).

We demonstrate on a very simple example from a bedrock detection how a structural interface from a seismic refraction can lead to models that are far easier to interpret.

# We import the numpy and matplotlib
import numpy as np
# Next we import pyGIMLi and the ERT and mesh tools modules
import pygimli as pg
import pygimli.meshtools as mt
from pygimli.physics import ert

First we load the data and the interface

sphinx_gallery_thumbnail_number = 4

data = pg.getExampleData("ert/struct.dat")
xz = pg.getExampleData("ert/struct.txt")
ert.show(data)
plot 4 structural constraints
(<Axes: >, <matplotlib.colorbar.Colorbar object at 0x7f71513c0340>)

We estimate an error and do a normal inversion using the ERT manager. As we expect more or less layered structures, we choose a vertical weighting factor for the regularization of 0.5.

data.estimateError()
mgr = ert.Manager(data)
mgr.invert(paraDepth=150, zWeight=0.5, verbose=True)
fop: <pygimli.physics.ert.ertModelling.ERTModelling object at 0x7f71206c72e0>
Data transformation: <pygimli.core._pygimli_.RTransLogLU object at 0x7f7120e732e0>
Model transformation (cumulative):
         0 <pygimli.core._pygimli_.RTransLogLU object at 0x7f7122fce6e0>
min/max (data): 55.69/508
min/max (error): 3%/3%
min/max (start model): 254/254
--------------------------------------------------------------------------------
inv.iter 0 ... chi² =  343.06
--------------------------------------------------------------------------------
inv.iter 1 ... chi² =    9.41 (dPhi = 96.93%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 2 ... chi² =    7.32 (dPhi = 21.07%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 3 ... chi² =    6.16 (dPhi = 14.83%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 4 ... chi² =    5.10 (dPhi = 15.48%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 5 ... chi² =    3.27 (dPhi = 30.19%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 6 ... chi² =    0.42 (dPhi = 66.64%) lam: 20.0


################################################################################
#                  Abort criterion reached: chi² <= 1 (0.42)                   #
################################################################################

803 [61.6412762770448,...,1465.6683757223302]

We plot the inversion result using a defined color scale. Additionally, we plot the measured line on top of it and see that generally there is a change from conductive overburden to resistive bedrock.

ax, cb = mgr.showResult(cMin=50, cMax=1500)
ax.plot(xz[:, 0], xz[:, 1], "k-")
plot 4 structural constraints
[<matplotlib.lines.Line2D object at 0x7f7151630d30>]

We now want to include this line into the inversion mesh. Therefore we generate the geometry of the inversion, a polygone from the line and merge both objects before we send it to the mesh generation.

plc = mt.createParaMeshPLC(data, paraDepth=150, boundary=1)
line = mt.createPolygon(xz, marker=1)
plc += line
mesh = mt.createMesh(plc, quality=34.3)
pg.show(mesh, markers=True, showMesh=True)
plot 4 structural constraints
(<Axes: >, <matplotlib.colorbar.Colorbar object at 0x7f7122ca2080>)

We see the outer boundaries have markers of -1 (homogeneous Neumann, no-flow) and -2 (mixed) important to define the boundary conditions of the forward problem. The interface between the inversion domain (yellow) and the outer region (blue) and the added interface are >0 and do not affect the forward. We now use the generated mesh in a new ERT manager by setting it.

mgr = ert.Manager(data)
mgr.setMesh(mesh)
mgr.invert(verbose=True)
ax, cb = mgr.showResult(cMin=50, cMax=1500)
plot 4 structural constraints
fop: <pygimli.physics.ert.ertModelling.ERTModelling object at 0x7f715093d670>
Data transformation: <pygimli.core._pygimli_.RTransLogLU object at 0x7f713f19d3f0>
Model transformation: <pygimli.core._pygimli_.RTransLog object at 0x7f711d320cc0>
min/max (data): 55.69/508
min/max (error): 3%/3%
min/max (start model): 254/254
--------------------------------------------------------------------------------
inv.iter 0 ... chi² =  343.06
--------------------------------------------------------------------------------
inv.iter 1 ... chi² =    9.96 (dPhi = 96.63%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 2 ... chi² =    7.87 (dPhi = 19.95%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 3 ... chi² =    5.80 (dPhi = 23.98%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 4 ... chi² =    2.16 (dPhi = 50.37%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 5 ... chi² =    0.70 (dPhi = 42.68%) lam: 20.0


################################################################################
#                  Abort criterion reached: chi² <= 1 (0.70)                   #
################################################################################

As a result, we see sharp resistivity jumps across the line in the center of the profile, i.e. the resistivity structure follows the given hints. At other positions like at the beginning of the profile, it does not show such strong contrasts indicating that the bedrock interface is probably more shallow and there is likely a weathered zone.

References Doetsch, J., Linde, N., Pessognelli, M., Green, A.G. & Günther, T. (2012): Constraining 3-D electrical resistance tomography with GPR data for improved aquifer characterization. Journal of Applied Geophysics 78, 68-76, doi:10.1016/j.jappgeo.2011.04.008. Wunderlich, T., Fischer, P., Wilken, D., Hadler, H., Erkul, E., Mecking, R., Günther, T., Heinzelmann, M., Vött, A. & Rabbel, W. (2018): Constraining Electric Resistivity Tomography by Direct Push Electric Conductivity logs and vibracores: An exemplary study of the Fiume Morto silted riverbed (Ostia Antica, Western Italy). Geophysics 83(3), B87-B103, doi:10.1190/geo2016-0660.1. Tanner, D.C., Buness, H., Igel, J., Günther, T., Gabriel, G., Skiba, P., Plenefisch, T., Gestermann, N. & Walter, T. (2020): Fault Detection. in: Tanner, C.D. & Brandes, C. (Eds.): Understanding Faults, 380p., Elsevier, 81-146, doi:10.1016/B978-0-12-815985-9.00003-5. Jiang, C., Igel, J., Dlugosch, R., Müller-Petke, M., Günther, T., Helms, J., Lang, J. & Winsemann (2020): Magnetic resonance tomography constrained by ground-penetrating radar for improved hydrogeophysical characterisation, Geophysics 85(6), JM13-JM26, doi:10.1190/geo2020-0052.1. Norooz, R., Olsson, P.-I., Dahlin, T., Günther, T. & Bernstone, C. (2021): A geoelectrical pre-study of Älvkarleby test embankment dam: 3D forward modelling and effects of structural constraints on the 3D inversion model of zoned embankment dams. J. Appl. Geophys. 191, 104355, doi:10.1016/j.jappgeo.2021.104355.

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

Gallery generated by Sphinx-Gallery