Note
Go to the end to download the full example code
3D crosshole ERT inversion#
In this example, we demonstrate the inversion of 3D crosshole field data. Instead of a regular grid or an irregular tetrahedral mesh, we use triangular prism mesh (triangles in x-y plane and regular along z). This is beneficial in cases of a predominant layering that can be accounted for by using the zWeight inversion parameter because the boundaries are perfectly vertical.
We import the used libraries pygimli, meshtools the ERT module and a function for displaying a data container.
import numpy as np
import pygimli as pg
import pygimli.meshtools as mt
from pygimli.physics import ert
from pygimli.viewer.mpl import showDataContainerAsMatrix
We load the data file from the example repository. It represents a crosshole data set published by Doetsch et al. (2010) in the frame of the RECORD project where boreholes were installed to monitor the exchange between river and groundwater (Coscia et al., 2008).
data = pg.getExampleData("ert/crosshole3d.dat")
print(data)
Data: Sensors: 36 data: 753, nonzero entries: ['a', 'b', 'm', 'n', 'r', 'valid']
There are 36 electrodes, each 9 in four boreholes, and 1256 data that are resistances only. Therefore we first compute the geometric factors and then the apparent resistivities, of which we plot the last few values.
data["k"] = ert.geometricFactors(data, dim=3)
data["rhoa"] = data["r"] * data["k"]
We plot the data in form of a crossplot between the A-B and M-N electrode combinations. Values are in the range between 100 and 500 Ohmmeters.
We first extract the borehole locations, i.e. the x and y positions of the electrodes. From these we create a rectangle with 40% boundary and marker 2 and add the borehole positions to it.
From this PLC, we create a mesh using a maximum cell size. We add an outer (modelling) boundary.
We create a vertical discretization vector with dense spacing in the range of the electrodes and a coarser discretization above and below.
[ -0. -0.8 -1.6 -2.4 -3.2 -4. -4.4 -4.8 -5.2 -5.6 -6. -6.4
-6.8 -7.2 -7.6 -8. -8.4 -8.8 -9.2 -9.6 -10. -10.8 -11.6 -12.4
-13.2 -14. -14.8]
From the 2d mesh and the z vector we create a 3D triangular prism mesh that obtains the marker of the 2D mesh. Additionally, we set all cells above or below to marker 1 which is by default the background region. In total we have 56k cells, of which most are background and less than 20k cells are inverted.
mesh = mt.createMesh3D(mesh2d, zVec, pg.core.MARKER_BOUND_HOMOGEN_NEUMANN,
pg.core.MARKER_BOUND_MIXED)
print(mesh)
for c in mesh.cells():
cd = -c.center().z() # center depth
if cd < dTop or cd > dBot:
c.setMarker(1)
mesh["region"] = pg.Vector(mesh.cellMarkers())
sli = mt.extract2dSlice(mesh, y=2.5)
ax, cb = pg.show(sli, "region", showMesh=True)
_ = ax.plot(pg.x(data), pg.z(data), "mo", markersize=1)
Mesh: Nodes: 11664 Cells: 21736 Boundaries: 3830
We estimate an error using default values, i.e. 3% relative error and an absolute error of 100uV at an assumed current of 100mA which is almost zero. Inversion is run with less weight into the vertical direction.
data.estimateError(relativeError=0.03, absoluteUError=1e-4)
# data["err"] = ert.estimateError(data) # the same
mgr = ert.Manager(data)
mgr.invert(mesh=mesh, zWeight=0.3, verbose=True)
fop: <pygimli.physics.ert.ertModelling.ERTModelling object at 0x7ff74d13f290>
Data transformation: <pygimli.core._pygimli_.RTransLogLU object at 0x7ff74d18a430>
Model transformation (cumulative):
0 <pygimli.core._pygimli_.RTransLogLU object at 0x7ff74ed0df60>
min/max (data): 82.22/548
min/max (error): 3%/3.1%
min/max (start model): 243/243
--------------------------------------------------------------------------------
inv.iter 0 ... chi² = 157.11
--------------------------------------------------------------------------------
inv.iter 1 ... chi² = 14.28 (dPhi = 89.85%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 2 ... chi² = 3.46 (dPhi = 60.78%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 3 ... chi² = 2.16 (dPhi = 21.95%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 4 ... chi² = 1.56 (dPhi = 7.91%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 5 ... chi² = 1.54 (dPhi = 0.47%) lam: 20.0
################################################################################
# Abort criterion reached: dPhi = 0.47 (< 2.0%) #
################################################################################
6015 [326.05337833337836,...,140.11615873813034]
Eventually, we are able to fit the data to a chi-square value close to 1, or about 3% RRMS. We visualize the result using the pyVista package with a clip. The electrodes are marked by magenta points. We mainly see a layering with highest resistivity on top and lowest at the bottom.
pd = mgr.paraDomain
pd["res"] = mgr.model
pl, _ = pg.show(pd, label="res", style="surface", cMap="Spectral_r", hold=True,
filter={"clip": dict(normal=[1, 1, 0], origin=[2, 2, -6])})
pl.add_points(data.sensors().array(), color="magenta")
pl.camera_position = "yz"
pl.camera.azimuth = 20
pl.camera.elevation = 20
pl.camera.zoom(1.2)
_ = pl.show()
We also have a closer look at a slice through the middle.
para = mgr.paraDomain
para["res"] = mgr.paraModel()
slice = mt.extract2dSlice(para, origin=[4, 4, 0], normal=[1, 1, 0])
pg.show(slice, "res", cMap="Spectral_r", cMin=100, cMax=500)
(<Axes: >, <matplotlib.colorbar.Colorbar object at 0x7ff757be7670>)
References#
Coscia, I., S. Greenhalgh, N. Linde, A. Green, T. Günter, J. Doetsch, and T. Vogt, 2010, A multi-borehole 3-D ERT monitoring system for aquifer characterization using river flood events as a natural tracer: Ext. Abstr. 16th Annual EAGE meeting of Environmental and Engineering Geophysics.
Doetsch, J., Linde, N., Coscia, I., Greenhalgh, S., & Green, A. (2010): Zonation for 3D aquifer characterization based on joint inversions of multimethod crosshole geophysical data. GEOPHYSICS (2010),75(6): G53.
Total running time of the script: (3 minutes 11.499 seconds)