2D crosshole ERT inversion

Inversion of 2D crosshole field data.

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 Kuras et al. (2009) in the frame of the ALERT project.

data = pg.getExampleData("ert/crosshole2d.dat")
print(data)

Data: Sensors: 144 data: 1256, nonzero entries: ['a', 'b', 'err', 'm', 'n', 'r', 'valid']

There are 144 electrodes, each 16 in nine 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"]
print(np.sort(data["rhoa"])[-5:])

[181.85035034 192.9622767  210.49617346 258.10205765 537.70069201]

We see that there is a single value above 500 Ohmm and a few other high ones. Therefore we delete the highest value. We plot the data in form of a crossplot between the A and M electrodes (with forward and backward sign).

data.remove(data["rhoa"] > 200)
m = np.sign(data["m"] - data["n"]) * data["m"]
showDataContainerAsMatrix(data, m, data["a"], "rhoa", cMap="Spectral_r")

(<matplotlib.axes._subplots.AxesSubplot object at 0x7fa06fa58280>, <matplotlib.colorbar.Colorbar object at 0x7fa074a930d0>)

We determine the x and z positions and create a regular grid with a spacing of 5xm that contains the electrodes as nodes. This is not necessary but improves quality of the forward response. Around the boreholes there is 0.5m space and all mesh cells have the marker 2.

ex = np.unique(pg.x(data))
ez = np.unique(pg.z(data))
dx = 0.05
nb = 8
xmin, xmax = min(ex) - nb*dx, max(ex) + nb*dx
zmin, zmax = min(ez) - nb*dx, 0
x = np.arange(xmin, xmax+.001, dx)
z = np.arange(zmin, zmax+.001, dx)
grid = mt.createGrid(x, z, marker=2)
ax, cb = pg.show(grid)
ax.plot(pg.x(data), pg.z(data), "mx")
print(grid)

Mesh: Nodes: 3977 Cells: 3840 Boundaries: 7816

In order to ensure correct boundary conditions, we append a triangular mesh outside that is automatically treated as background because they have a region marker of 1 and so-called world boundary conditions, i.e. Neumann BC at the Earth’s surface and mixed boundary conditions at the other boundaries.

mesh = mt.appendTriangleBoundary(grid, marker=1, boundary=5, worldMarkers=1)
pg.show(mesh, markers=True)

(<matplotlib.axes._subplots.AxesSubplot object at 0x7fa014c4fd00>, <matplotlib.colorbar.Colorbar object at 0x7fa06fe104c0>)

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 half as much weight into the vertical direction.

data["err"] = ert.estimateError(data)
mgr = ert.Manager(data)
mgr.inv.setRegularization(correlationLengths=[1, 0.5])
mgr.invert(mesh=mesh, verbose=True)
mgr.showResult(cMin=15, cMax=200)

fop: <pygimli.physics.ert.ertModelling.ERTModelling object at 0x7fa05564c680>
Data transformation: <pygimli.core._pygimli_.RTransLogLU object at 0x7fa029fc71d0>
Model transformation: <pygimli.core._pygimli_.RTransLog object at 0x7fa029d8e400>
min/max (data): 23.39/193
min/max (error): 3%/3%
min/max (start model): 68.62/68.62
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
inv.iter 1 ... chi² = 6.07 (dPhi = 92.58%) lam: 20
--------------------------------------------------------------------------------
inv.iter 2 ... chi² = 2.8 (dPhi = 39.28%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 3 ... chi² = 2.64 (dPhi = 1.53%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 4 ... chi² = 2.45 (dPhi = 0.79%) lam: 20.0
################################################################################
#                 Abort criteria reached: dPhi = 0.79 (< 2.0%)                 #
################################################################################

(<matplotlib.axes._subplots.AxesSubplot object at 0x7fa0224ea160>, <matplotlib.colorbar.Colorbar object at 0x7fa0554ac9a0>)

References

Kuras, O., Pritchard, J., Meldrum, P. I., Chambers, J. E., Wilkinson, P. B., Ogilvy, R. D.,and Wealthall, G. P. (2009). Monitoring hydraulic processes with Automated time-Lapse Electrical Resistivity Tomography (ALERT). Compte Rendus Geosciences - Special issue on Hydrogeophysics, 341(10-11):868–885.