2D/3D cross-borehole ERT forward and inversion using pygimli and pyBERT

[XueyuanK]

Hello, 1) I wonder whether pyGIMLi and pyBERT (not BERT with cmd line) can solve the large-scale 3D cross-borehole ERT and IP forward and inversion problem. 2) Because both pyGIMLi and pyBERT only provide ERT scheme (such as dipole-dipole) for surface electrodes (e.g., pb.createData), I created the bipole-bipole scheme for borehole electrodes from other ERT software (i.e., AGI Earthimager) and then import this scheme into pygimli with the function “. resize” and “.set”, like this:

read the scheme for borehole from agi software

read_scheme=np.loadtxt('schemedata.dat') #schemedata.dat is created by AGI software a_scheme=numpy2gvec(read_scheme[:,0]-1) #transform the numpy array into pgimli vector b_scheme=numpy2gvec(read_scheme[:,1]-1) m_scheme=numpy2gvec(read_scheme[:,2]-1) n_scheme=numpy2gvec(read_scheme[:,3]-1)

reset the scheme

scheme.resize(len(a_scheme)) scheme.set('a',a_scheme) scheme.set('b',b_scheme) scheme.set('m',m_scheme) scheme.set('n',n_scheme) scheme.set('valid',pg.RVector(len(a_scheme),1)) Besides, I changed the position of electrodes with “setSensorPositions”. pyGIMLi can solve the forward problem here, but it will go wrong when running the inversion. I wonder whether the mesh I built is not suitable for inversion. Both the forward and inversion use the meshERT: meshERT=mt.createGrid(x=np.linspace(-10,20,61),y=np.linspace(-50,0,101)) 3) Besides, when in 2D mode, should the borehole electrodes coordinates be on z or y?

Thanks in advance! Bests, Xueyuan

code&scheme.zip

[halbmy]

I had a look into your scripts.

1. The line map = np.array(scheme.dataMap()) can never work as the dataMap is sort of a dictionary (you can find its keys/datafields with DM.keys() and access them it with DM['rhoa']) but it the line is not needed anyway. (Besides better not use the Python keyword map for a variable name).

2. Both for forward calculation and inversion you will need a larger mesh to avoid problems from boundary conditions. The function pg.meshtools.appendTriangleBoundary() helps doing that. (I suggest adding an option to ert.simulate that does it automatically).

3. You should not use the same mesh for forward calculation and inversion as this is considered an inverse crime. At any rate, only invert for an inner part and use the outer box as (prolongated) background.

4. Instead of using a random model I suggest creating a geostatistical model by res = pg.utils.geostatistics.generateGeostatisticalModel(meshERT) * 10 + 100

5. A noiseLevel=0 is currently indicating "no noise added" which is finally spoiling the inversion by infinite weighting. Better use a few percent of relative noise.

[XueyuanK]

Thanks for your suggestions! I have tried the pg.meshtools.appendTriangleBoundary() to help avoid the boundary condition, prepared coarser mesh for inversion and set noiseLevel=1. However, the inversion still can not work properly.

I wonder whether it is due to the "Setting topography=1" in the screenshot below, but I did not set any topography here.

Besides, it seems that the report shows the presence of negative resistivity, and I try to use data.removeInvalid() to remove it. However, the number of data points has not decreased. "found neg. resp, save and abort. 236 -119.453 -119.036 10 16 13 19"

[halbmy]

The algorithm says "found Neumann domain" which means that all boundaries are of homogeneous Neumann (no-flow) type (and therefore non-trivial). However, in the function appendTriangleBoundary https://github.com/gimli-org/gimli/blob/bc43dbcc50d707ee616f72ea13adce559b8c765c/python/pygimli/meshtools/grid.py#L63 we use mixed boundaries:

        poly.createEdge(n1, n2, pg.MARKER_BOUND_HOMOGEN_NEUMANN)
        poly.createEdge(n2, n3, pg.MARKER_BOUND_MIXED)
        poly.createEdge(n3, n4, pg.MARKER_BOUND_MIXED)
        poly.createEdge(n4, n1, pg.MARKER_BOUND_MIXED)

So something goes wrong here. We can only help if you provide the script.

[XueyuanK]

Hello, Sorry for the late reply. Here is the test script I used. 2d_borehole_inversion.zip

[halbmy]

Sorry for the late feedback, I just have been too busy. Actually this is a nice example. I took your code and simplified it a bit.

First, you should use to use mt.createWorld() instead of mt.createRectangle() for the synthetic model, because this will make sure that the proper boundary conditions (upper surface Neumann, otherwise mixed) will be used and that the computation is much more accurate due to secondary potential (which cannot be used in a Neumann body that would be produced by mt.createRectangle().

Second, you should give the markers in the grid (inversion region) a marker>1. pg.show(mesh, markers=True) shows that the outside region has the marker 1 but the inner region has marker 0. By default, the smallest region is used as background and all other for inversion, which will declare the inner region to background. If you use: inner = pg.createGrid(np.arange(-10, 20, 0.5), -np.arange(0.0, 30, 0.5), marker=2) all will be working well.

# try to use the pygimli to forward the ERT simulation
# 2D cross-borehole electrodes forward and inversion
import numpy as np
import pygimli as pg
import pygimli.meshtools as mt
import pybert as pb

scheme = pb.createData(elecs=pg.utils.grange(start=0, end=25, dx=1.0),
                       schemeName='dd')
# change the electrodes position
pos = np.array(scheme.sensorPositions())  # get xyz coor for electrodes
# borehole 1
pos[0:10, 0] = 0
for i in range(1, 11, 1):
    pos[i-1, 1] = -2*i

# surface
for i in range(0, 6, 1):
    pos[i+10, 0] = 2*i

# borehole 2
pos[16:26, 0] = 10
for i in range(1, 11, 1):
    pos[i+15, 1] = -2*i

scheme.setSensorPositions(pos)  # reset the position for electrodes
pos_xy = pos[:, 0:2]

# read the scheme for borehole from agi software
read_scheme = np.loadtxt('schemedata.dat')  # created by AGI software
scheme.resize(read_scheme.shape[0])
for i, tok in enumerate(['a', 'b', 'm', 'n']):
    scheme.set(tok, read_scheme[:, i]-1)

scheme.set('valid', pg.Vector(scheme.size(), 1))

# create the mesh for forward
world = mt.createWorld(start=[-40, -100], end=[50, 0], marker=0)
c0 = mt.createCircle(pos=(5.0, -10.0), radius=4, segments=25, marker=1)
geom = mt.mergePLC([world, c0])
for sen in pos_xy:
    geom.createNode(sen)  # create new node for electrodes

mesh_fwd = mt.createMesh(geom, quality=34, area=.5)
model = np.array([10, 2300])[mesh_fwd.cellMarkers()]
ert = pb.ERTManager()
data = ert.simulate(mesh_fwd, res=model, scheme=scheme,
                    noiseLevel=0.01, noiseAbs=1e-6)
data.save('data.dat', 'a b m n rhoa k err')
# plot apparent resistivity over a and m index
if 0:
    pg.mplviewer.plotDataContainerAsMatrix(data, 'a', 'm', 'rhoa')
    raise SystemExit

# inversion
inner = pg.createGrid(np.arange(-10, 20, 0.5), -np.arange(0.0, 30, 0.5),
                      marker=2)

mesh = mt.appendTriangleBoundary(inner, xbound=30, ybound=82)
pg.show(mesh, markers=True)
mesh.save("mesh.bms")
mesh.exportBoundaryVTU('boundary.vtu')
# %% inversion
model = ert.invert(data, mesh=mesh, lam=20)
ert.paraDomain.save("pd.bms")
model.save("model.vec")
ert.showResult()
pg.wait()

[XueyuanK]

Hello, Thanks for your feedback! Both the forward and inversion work well!! 1）However, I still have some questions about the mesh for forward model. Actually, I try to use pygimli as a geophysical forward model for coupled hydrogeophysical inversion. I need to couple the hydrogeological model (e.g. MODFLOW and MT3DMS) and the geophysical model (i.e., pygimli). Therefore, I need to input the resistivity values calculated from the MT3DMS (or TOUGH) and petrophysical relationship (e.g., Archie’s law) into the pygimli in the form of rectangular grids. So, I try to use mesh_fwd=mt.createGrid(x=np.linspace(-40,50,91),y=np.linspace(-100,0,101)) to create mesh for forward model. However, the forward model cannot work properly. As you mentioned above, maybe it is because the boundary conditions are set incorrectly in mt.createGrid. So, how to set proper BCs for forward modeling? 2）Besides, I would like to use pyglimi to solve 3D ERT cross-borehole forward and inversion problem (also in coupled hydrophysical inversion framework). I try to use mt.createGrid and mt.appendTetrahedronBoundary to create a grids for forward modeling, but something went wrong. So, I wonder how to create proper hexahedral grids for 3D modeling. inner_fwd=mt.createGrid(x=np.linspace(-10,30,41),y=np.linspace(0,15,16),z=np.linspace(-20,0,21)) mesh_fwd=mt.appendTetrahedronBoundary(inner_fwd, xbound=75, ybound=16,zbound=80) #error here 3）Questions about parameter setting: how to properly set the and noiseAbs in forward modeling? And lamda in the inversion? 2D3D_ERT_modeling.zip

Attached pleased the scripts. Thanks for your help! Bests, Xueyuan

[halbmy]

1. the hydro models are usually not big enough for a good forward calculation. Therefore mt.appendTriangleBoundary is always a good idea to extend the grid with a triangle mesh. As a side effects, the boundary conditions (top=Neumann, others=Mixed) will be set automatically. Otherwise you you would have to search for all boundaries at xmin/xmax/zmax to set them

   for b in grid.boundaries():
   if b.center().x() == xmin or ....
       b.setMarker(pg.MARKER_BC_MIXED)...

2. Unlike appendTriangleBoundary, appendTetrahedronBoundary is not full stable and is sometimes resulting in additional nodes on the boundary. I cannot tell you without seeing the error.
3. The absolute error represents higher errors for lower voltage gains. You can take the default values, determine some from your field data (by reciprocal analysis), have a look into literature or just set them to zero. Particularly for synthetic data, the regularization parameter lam needs to be set so that the data are fitted within noise. There is also chi2=1 optimization scheme for adjusting that.

[halbmy]

Please note that the basic functionality of BERT has been moved into the ert module, so instead of pb.createData or pb.ERTManager() you should write

from pygimli.physics import ert
data = ert.createData(...
mgr = ert.ERTManager()

[simondehout]

Hi Thomas! Thank you for the information and solution code. I've tried to run it in spyder without importing BERT but I receive the following error: AttributeError: 'DataContainerERT' object has no attribute 'fop' The error comes from invert() function on line 64. I've imported it from VESManager as it is the only module which has the invert() function (I believe). My code is as follows:

`# try to use the pygimli to forward the ERT simulation
# 2D cross-borehole electrodes forward and inversion
import numpy as np
import pygimli as pg
import pygimli.meshtools as mt
from pygimli.physics import ert

scheme = ert.createData(elecs=pg.utils.grange(start=0, end=25, dx=1.0),
                       schemeName='dd')
# change the electrodes position
pos = np.array(scheme.sensorPositions())  # get xyz coor for electrodes
# borehole 1
pos[0:10, 0] = 0
for i in range(1, 11, 1):
    pos[i-1, 1] = -2*i

# surface
for i in range(0, 6, 1):
    pos[i+10, 0] = 2*i

# borehole 2
pos[16:26, 0] = 10
for i in range(1, 11, 1):
    pos[i+15, 1] = -2*i

scheme.setSensorPositions(pos)  # reset the position for electrodes
pos_xy = pos[:, 0:2]

# read the scheme for borehole from agi software
read_scheme = np.loadtxt('schemedata.dat')  # created by AGI software
scheme.resize(read_scheme.shape[0])
for i, tok in enumerate(['a', 'b', 'm', 'n']):
    scheme.set(tok, read_scheme[:, i]-1)

scheme.set('valid', pg.Vector(scheme.size(), 1))

# create the mesh for forward
world = mt.createWorld(start=[-40, -100], end=[50, 0], marker=0)
c0 = mt.createCircle(pos=(5.0, -10.0), radius=4, segments=25, marker=1)
geom = mt.mergePLC([world, c0])
for sen in pos_xy:
    geom.createNode(sen)  # create new node for electrodes

mesh_fwd = mt.createMesh(geom, quality=34, area=.5)
model = np.array([10, 2300])[mesh_fwd.cellMarkers()]
mgr = ert.ERTManager()
data = mgr.simulate(mesh_fwd, res=model, scheme=scheme,
                    noiseLevel=0.01, noiseAbs=1e-6)
data.save('data.dat', 'a b m n rhoa k err')
# plot apparent resistivity over a and m index
if 0:
    pg.mplviewer.plotDataContainerAsMatrix(data, 'a', 'm', 'rhoa')
    raise SystemExit

# inversion
inner = mt.createGrid(np.arange(-10, 20, 0.5), -np.arange(0.0, 30, 0.5),
                      marker=2)

mesh = mt.appendTriangleBoundary(inner, xbound=30, ybound=82)
pg.show(mesh, markers=True)
mesh.save("mesh.bms")
mesh.exportBoundaryVTU('boundary.vtu')
# %% inversion
model = ert.VESManager.invert(data, mesh=mesh, lam=20)
ert.paraDomain.save("pd.bms")
model.save("model.vec")
ert.showResult()
pg.wait()`

Let me know whether I might a problem with my inputs or the version of pygimli (1.3.1) Thank you!

[halbmy]

I would need the file schemedata.dat to check the script.

[halbmy]

At least it will not work to use VESManager.invert(data) with a DataContainerERT as the invert call now accepts only a data vector. However, it should work like this, so I will have to do some modifications.

[simondehout]

Hi Thank you for the reply. I used the your code (April 16, 2019) and the schemedata file given by the initial post. I'm trying to see whether I can get the same result/figure as you posted previously.

Thanks in advance,

[simondehout]

I actually get this image as output. But it seems quite different compared to yours. My corrected code looks like this.

# try to use the pygimli to forward the ERT simulation
# 2D cross-borehole electrodes forward and inversion
import numpy as np

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

scheme = ert.createData(elecs=pg.utils.grange(start=0, end=25, dx=1.0),
                       schemeName='dd')
# change the electrodes position
pos = np.array(scheme.sensorPositions())  # get xyz coor for electrodes
# borehole 1
pos[0:10, 0] = 0
for i in range(1, 11, 1):
    pos[i-1, 1] = -2*i

# surface
for i in range(0, 6, 1):
    pos[i+10, 0] = 2*i

# borehole 2
pos[16:26, 0] = 10
for i in range(1, 11, 1):
    pos[i+15, 1] = -2*i

scheme.setSensorPositions(pos)  # reset the position for electrodes
pos_xy = pos[:, 0:2]

# read the scheme for borehole from agi software
read_scheme = np.loadtxt('schemedata.dat')  # created by AGI software
scheme.resize(read_scheme.shape[0])
for i, tok in enumerate(['a', 'b', 'm', 'n']):
    scheme.set(tok, read_scheme[:, i]-1)

scheme.set('valid', pg.Vector(scheme.size(), 1))

# create the mesh for forward
world = mt.createWorld(start=[-40, -100], end=[50, 0], marker=0)
c0 = mt.createCircle(pos=(5.0, -10.0), radius=4, segments=25, marker=1)
geom = mt.mergePLC([world, c0])
for sen in pos_xy:
    geom.createNode(sen)  # create new node for electrodes

mesh_fwd = mt.createMesh(geom, quality=34, area=.5)
model = np.array([10, 25])[mesh_fwd.cellMarkers()]
mgr = ert.ERTManager()
data = mgr.simulate(mesh_fwd, res=model, scheme=scheme,
                    noiseLevel=0.01, noiseAbs=1e-6)
data.save('data.dat', 'a b m n rhoa k err')
# plot apparent resistivity over a and m index
if 0:
    pg.mplviewer.plotDataContainerAsMatrix(data, 'a', 'm', 'rhoa')
    raise SystemExit

# inversion
inner = mt.createGrid(np.arange(-10, 20, 0.5), -np.arange(0.0, 30, 0.5),
                      marker=2)

mesh = mt.appendTriangleBoundary(inner, xbound=30, ybound=82)
pg.show(mesh, markers=True)
mesh.save("mesh.bms")
mesh.exportBoundaryVTU('boundary.vtu')
# %% inversion
model = mgr.invert(data, mesh=mesh, lam=20)
mgr.paraDomain.save("pd.bms")
model.save("model.vec")
mgr.showResult()
pg.wait()

[halbmy]

Your inversion mesh does not make sense (or appendTriangleMarker needs to be modified): If you change the y vector to np.arange(-30.0, .1, 0.5) it will make sense:

[halbmy]

My result is the following:

Note that I changed the area argument in the forward modelling from .5 to 3 because it caused an extremely fine mesh.

[halbmy]

With blockyModel=True it becomes even better:

[simondehout]

Great! Thank you. It looks better indeed. Could you maybe share with me the full version of your code to generate the last figure? Thanks a lot :)

[halbmy]

"""2D cross-borehole electrodes forward and inversion."""
import numpy as np

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

scheme = ert.createData(elecs=pg.utils.grange(start=0, end=25, dx=1.0),
                        schemeName='dd')
# change the electrodes position
pos = np.array(scheme.sensorPositions())  # get xyz coor for electrodes
# borehole 1
pos[0:10, 0] = 0
for i in range(1, 11, 1):
    pos[i-1, 1] = -2*i
# surface
for i in range(0, 6, 1):
    pos[i+10, 0] = 2*i
# borehole 2
pos[16:26, 0] = 10
for i in range(1, 11, 1):
    pos[i+15, 1] = -2*i

scheme.setSensorPositions(pos)  # reset the position for electrodes
pos_xy = pos[:, 0:2]
# read the scheme for borehole from agi software
read_scheme = np.loadtxt('schemedata.dat')  # created by AGI software
scheme.resize(read_scheme.shape[0])
for i, tok in enumerate(['a', 'b', 'm', 'n']):
    scheme.set(tok, read_scheme[:, i]-1)

scheme.set('valid', pg.Vector(scheme.size(), 1))
# create the mesh for forward
world = mt.createWorld(start=[-40, -100], end=[50, 0], marker=0)
c0 = mt.createCircle(pos=(5.0, -10.0), radius=4, segments=25, marker=1)
geom = mt.mergePLC([world, c0])
for sen in pos_xy:
    geom.createNode(sen)  # create new node for electrodes

pg.show(geom)
# %%
mesh_fwd = mt.createMesh(geom, quality=34, area=3)
pg.show(mesh_fwd)
# %%
model = np.array([10, 25])[mesh_fwd.cellMarkers()]
mgr = ert.ERTManager()
data = mgr.simulate(mesh_fwd, res=model, scheme=scheme,
                    noiseLevel=0.01, noiseAbs=1e-6)
data.save('data.dat', 'a b m n rhoa k err')
# plot apparent resistivity over a and m index
if 0:
    pg.mplviewer.plotDataContainerAsMatrix(data, 'a', 'm', 'rhoa')
    raise SystemExit

# %% inversion
inner = mt.createGrid(np.arange(-10, 20, 0.5), np.arange(-30.0, .1, 0.5),
                      marker=2)
mesh = mt.appendTriangleBoundary(inner, xbound=30, ybound=82)
ax, _ = pg.show(mesh, markers=True)
ax.plot(pg.x(scheme), pg.y(scheme), "r.")
# %%
mesh.save("mesh.bms")
mesh.exportBoundaryVTU('boundary.vtu')
# %% inversion
model = mgr.invert(data, mesh=mesh, lam=20, blockyModel=True)
mgr.showResult()
# %%
mgr.paraDomain.save("pd.bms")
model.save("model.vec")
pg.wait()

[kafahabibullah1453]

Please note that the basic functionality of BERT has been moved into the ert module, so instead of pb.createData or pb.ERTManager() you should write

from pygimli.physics import ert
data = ert.createData(...
mgr = ert.ERTManager()

Excuse me sir, is there any reason, why the functionality from BERT is moved to pygimli?

[halbmy]

Pure ERT (and simple IP) modelling and inversion is so widespread that it should be part of pyGIMLi without the need for other packages. The pyBERT project is mainly for spectral IP (time or frequency domain) and for extensive examples.