ERT inversion of data measured on a standing tree#

This example is about inversion of ERT data around a hollow lime tree. It applies the BERT methodology of Günther et al. (2006) to trees, first documented by Göcke et al. (2008) and Martin & Günther (2013). It has already been used in the BERT tutorial (Günther & Rücker, 2009-2023) in chapter 3.3 (closed geometries).

generate circular geometry
geometric factor generation
circular data representation
inversion

import pygimli as pg
import pygimli.meshtools as mt
from pygimli.physics import ert

Get some example data with topography, typically by a call like data = ert.load(“filename.dat”) that supports various file formats

data = pg.getExampleFile('ert/hollow_limetree.ohm', load=True, verbose=True)
print(data)

Data: Sensors: 24 data: 264, nonzero entries: ['a', 'b', 'i', 'm', 'n', 'r', 'u', 'valid']

It uses 24 electrodes around the circumference and applies a dipole-dipole (AB-MN) array. For each of the 24 dipoles being used as current injection (AB), 21 potential dipoles (from 3-4 to 23-24) can be measured of which 10 are measured, so that we end up in a total of 24*11=264 measurements. Apart from current (‘a’ and ‘b’) and potential (‘m’, ‘n’) electrodes, the file contains current (‘i’), voltage (‘u’) and resistance (‘r’) vectors for each of the 264 data.

We first generate the geometry by creating a close polygon. Between each two electrodes, we place three additional nodes whose positions are interpolated using a spline.

plc = mt.createPolygon(data.sensors(), isClosed=True,
                       addNodes=3, interpolate='spline')
ax, cb = pg.show(plc)

From this geometry, we create a triangular mesh with a quality factor, a maximum triangle area and post-smoothing.

mesh = mt.createMesh(plc, quality=34.3, area=2e-4, smooth=[10, 1])
print(mesh)
ax, _ = pg.show(mesh)
for i, s in enumerate(data.sensors()):
    ax.text(s.x(), s.y(), str(i+1), zorder=100)

Mesh: Nodes: 881 Cells: 1660 Boundaries: 2540

We first create the geometric factors to multiply the resistances to obtain “mean” resistivities for every quadrupol. We do this numerically using a refined mesh with quadratic shape functions using a constant resistivity of 1. The inverse of the modelled resistances is the geometric factor so that all apparent resistivities become 1. We then generate apparent resistivity, store it in the data and show it.

data["k"] = ert.createGeometricFactors(data, mesh=mesh, numerical=True)
data["rhoa"] = data["r"] * data["k"]
ax, cb = ert.show(data, "rhoa", circular=True)

/Users/fwagner/git/gimli_dev/source/pygimli/physics/ert/visualization.py:347: RuntimeWarning: invalid value encountered in cast
  ne = np.array(xe/de, dtype=int)

Note that we use the option circular. The values have almost rotation symmetry with lower values for shallow penetrations and higher ones for larger current-voltage distances.

We estimate a measuring error using default values (3%) and feed it into the ERT Manager.

data.estimateError()
mgr = ert.Manager(data)
mgr.invert(mesh=mesh, verbose=True)

fop: <pygimli.physics.ert.ertModelling.ERTModelling object at 0x29b949c10>
Data transformation: <pygimli.core._pygimli_.RTransLogLU object at 0x29199ea20>
Model transformation: <pygimli.core._pygimli_.RTransLog object at 0x29edcda30>
min/max (data): 41.73/1121
min/max (error): 3%/3%
min/max (start model): 205/205
--------------------------------------------------------------------------------
inv.iter 0 ... chi² =  436.39
--------------------------------------------------------------------------------
inv.iter 1 ... chi² =   88.06 (dPhi = 78.77%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 2 ... chi² =    3.54 (dPhi = 90.80%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 3 ... chi² =    2.17 (dPhi = 12.95%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 4 ... chi² =    2.11 (dPhi = 0.51%) lam: 20.0
################################################################################
#                Abort criterion reached: dPhi = 0.51 (< 2.0%)                 #
################################################################################

1660 [199.61760447433892,...,285.1723834882798]

We achieve a fairly good data fit with a chi-square value of 2 and compare data and model response, both are looking very similar.

ax = mgr.showFit(circular=True)

Finally, we have a look at the resulting resistivity distribution. It clearly shows cavity inside as high resistivity.

ax, cb = mgr.showResult(cMin=100, cMax=500, coverage=1)

References#

Günther, T., Rücker, C. & Spitzer, K. (2006): Three-dimensional modeling & inversion of dc resistivity data incorporating topography – II: Inversion. Geophys. J. Int. 166, 506-517, doi:10.1111/j.1365-246X.2006.03011.x
Göcke, L., Rust, S., Weihs, U., Günther, T. & Rücker, C. (2008): Combining sonic and electrical impedance tomography for the nondestructive testing of trees. - Western Arborist, Spring/2008: 1-11.
Martin, T. & Günther, T. (2013): Complex Resistivity Tomography (CRT) for fungus detection on standing oak trees. European Journal of Forest Research 132(5), 765-776, doi:10.1007/s10342-013-0711-4

Gallery generated by Sphinx-Gallery