HELP with ERT Crosshole Inversión. #820

aldefern · 2025-02-17T16:29:44Z

aldefern
Feb 17, 2025

Hello. I am working with some resistivity measurement data taken in two wells 15 meters apart from each other. In each well I have 36 electrodes separated by 0.7 m and the measurements have been taken in a "roll-along" manner. I have created the grid correctly but I cannot do the inversion because I cannot square the electrode positions. I am leaving you the notebook and the data to see if anyone has time and can take a look at it. Thanks in advance.
PS: If you think I could improve something, please let me know without any problem. This is for a specific time, in the future I will do a timelapse comparing it with other measurements.
ERT_Borehole .zip

Answered by halbmy

Feb 18, 2025

Your mesh is having the outer boundary with marker 2 and the inner with marker 1 (by default it is the other way),

so you need to tell which is the inversion and which is background by

mgr = ert.ERTManager(data, verbose=True)
mgr.setMesh(mesh)
mgr.inv.setRegularization(1, background=False)
mgr.inv.setRegularization(2, background=True)
model = mgr.invert(lam=20, verbose=True)

Then everything works well technically

However, something is wrong with your data file data.dat as for a lot of data identical indices are occurring:

# a b m n ...
1	2	1	3	
2	3	1	3	
3	3	2	4	
3	4	3	5	
4	5	3	5	
5	5	4	6	
5	6	5	7	
6	7	5	8	
7	8	6	8

There is something wrong.

View full answer

aldefern · 2025-02-17T19:10:53Z

aldefern
Feb 17, 2025
Author

Ok. I've made some changes and now it recognizes the electrodes. The problem now is that I cannot do the inversion.
HERE IS THE CODE:

import numpy as np
import pandas as pd
import pygimli as pg
from pygimli.physics import ert
from pygimli.viewer.mpl import showDataContainerAsMatrix
import matplotlib.pyplot as plt
import pygimli.meshtools as mt

# =============================
# 1️⃣ Cargar Datos de los Pozos
# =============================

# Nombre del archivo con los datos
file_name = "N1N3-005.txt"

# Cargar datos del Pozo N1 (primeras 33 filas)
df_N1 = pd.read_csv(file_name, delim_whitespace=True, skiprows=1, header=None, usecols=range(4, 31), nrows=33)

# Cargar datos del Pozo N3 (siguientes 33 filas)
df_N3 = pd.read_csv(file_name, delim_whitespace=True, skiprows=34, header=None, usecols=range(4, 31), nrows=33)

# Definir los nombres de las columnas
column_names = [
    "Spa.1", "Spa.2", "Spa.3", "Spa.4", "Rho", "Dev.", "M", "Sp", "Vp", "In", "Time",
    "Spa.5", "Spa.6", "Spa.7", "Spa.8", "Stack", "Rs-Check", "Vab", "Pab", "Rab", "Name",
    "Channel", "Overload", "Tx-Bat", "Rx-Bat", "Temp.", "Date"
]

# Asignar los nombres de columnas
df_N1.columns = column_names
df_N3.columns = column_names

# =============================
# 2️⃣ Cálculo de Resistividad Aparente
# =============================

factor_geom = 8.79642  # Factor geométrico

# Resistividad aparente para cada pozo
rhoa_N1 = df_N1["Vp"].values / df_N1["In"] 
rhoa_N3 = df_N3["Vp"].values / df_N3["In"] 

# Filtro de resistividades negativas
rhoa_N1[rhoa_N1 < 1] = 1
rhoa_N3[rhoa_N3 < 1] = 1



# Graficar resistividad aparente
plt.plot(rhoa_N1*factor_geom, 'r-', label="Pozo N1")
plt.title("Resistividad aparente - Pozo 1")
plt.show()

plt.plot(rhoa_N3*factor_geom, 'b-', label="Pozo N3")
plt.title("Resistividad aparente - Pozo 3")
plt.show()

# =============================
# 3️⃣ Preparación de Datos para Inversión
# =============================

# Extraer electrodos de cada pozo
M_N1, N_N1, A_N1, B_N1 = [np.abs(df_N1[col].values) for col in ["Spa.1", "Spa.2", "Spa.3", "Spa.4"]]
M_N3, N_N3, A_N3, B_N3 = [np.abs(df_N3[col].values) for col in ["Spa.1", "Spa.2", "Spa.3", "Spa.4"]]

# Concatenar datos de ambos pozos
M, N, A, B = [np.concatenate((var_N1, var_N3)) for var_N1, var_N3 in zip([M_N1, N_N1, A_N1, B_N1], [M_N3, N_N3, A_N3, B_N3])]
R = np.concatenate((rhoa_N1, rhoa_N3))





# Crear DataContainerERT para la inversión
data_container = pg.DataContainerERT()

data_container["a"], data_container["b"], data_container["m"], data_container["n"], data_container["r"]= A, B, M, N, R
data_container["k"] = np.ones_like(M)*factor_geom
data_container["rhoa"] = data_container["r"] * data_container["k"]
# Definir un error relativo del 5% con un mínimo de 1e-6 para evitar ceros
# Definir el error como un porcentaje fijo dentro del rango deseado
err_min = 0.03   # 3%
err_max = 0.0315 # 3.15%

data_container["err"] = np.clip(0.03 * data_container["rhoa"], err_min, err_max)





# =============================
# 4️⃣ Definir Posiciones de Electrodos
# =============================

num_electrodos = 36  # Número de electrodos por pozo
sep_electrodos = 0.7  # Separación entre electrodos (m)
dist_pozos = 15  # Distancia entre los pozos (m)

# Generar posiciones de electrodos
electrode_positions = [pg.RVector3(-1*dist_pozos/2,-1*i * sep_electrodos,0 ) for i in range(num_electrodos)] + \
                      [pg.RVector3(dist_pozos/2, -1*i * sep_electrodos,0) for i in range(num_electrodos)]

# Asignar posiciones al contenedor de datos
data_container.setSensorPositions(electrode_positions)
print(data_container)

# =============================
# 5️⃣ Visualización de la Matriz de Datos
# =============================

m = data_container["m"]
ax, cb = showDataContainerAsMatrix(data_container, m, data_container["a"], "rhoa", cMap="Spectral_r")

# =============================
# 6️⃣ Creación de Malla para la Inversión
# =============================

# Parámetros para la malla
n_electrodos = 36
sep_pozos = 15

# Coordenadas de los electrodos
z_electrodos = np.linspace(0, -(n_electrodos - 1) * sep_electrodos, n_electrodos)
x_pozos = [-sep_pozos / 2, sep_pozos / 2]

data_x = np.concatenate([np.full(n_electrodos, x) for x in x_pozos])
data_z = np.tile(z_electrodos, 2)

datax = pg.RVector(data_x)
dataz = pg.RVector(data_z)
data = np.column_stack((data_x, data_z))

# Crear la malla
ex, ez = np.unique(datax), np.unique(dataz)
dx, nb = 0.7, 4
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=1)
grid = mt.appendTriangleBoundary(grid, marker=2, boundary=4, worldMarkers=True)

# Marcar la base de la malla
for b in grid.boundaries():
    if np.isclose(b.center().y(), -1.7):
        b.setMarker(1)

# Visualizar la malla
ax, cb = pg.show(grid, showMesh=True, markers=True)
ax.plot(pg.x(data_container), data_z, "mx")

print(grid)

# =============================
# 7️⃣ Resistivity Inversion
# =============================

# Create inversion object
inv = ert.ERTManager(data_container, verbose=True)

# Run inversion with the created mesh
model = inv.invert(mesh=grid, lam=20)

# Visualize the inversion
fig, ax = plt.subplots(figsize=(8, 6))
pg.show(grid, model, ax=ax, label="Resistivity (Ωm)", cMap="Spectral_r")
plt.title("Resistivity Model Obtained by Inversion")
plt.show()

Here is the error I get:

RuntimeError                              Traceback (most recent call last)
Cell In[104], [line 9](vscode-notebook-cell:?execution_count=104&line=9)
      [6](vscode-notebook-cell:?execution_count=104&line=6) inv = ert.ERTManager(data_container, verbose=True)
      [8](vscode-notebook-cell:?execution_count=104&line=8) # Run inversion with the created mesh
----> [9](vscode-notebook-cell:?execution_count=104&line=9) model = inv.invert(mesh=grid, lam=20)
     [11](vscode-notebook-cell:?execution_count=104&line=11) # Visualize the inversion
     [12](vscode-notebook-cell:?execution_count=104&line=12) fig, ax = plt.subplots(figsize=(8, 6))

File c:\Users\Cuenta\anaconda3\envs\pg\Lib\site-packages\pygimli\frameworks\methodManager.py:843, in MeshMethodManager.invert(self, data, mesh, startModel, **kwargs)
    [837](file:///C:/Users/Cuenta/anaconda3/envs/pg/Lib/site-packages/pygimli/frameworks/methodManager.py:837) # if "blockyModel" in kwargs:  # redundant, now in fw
    [838](file:///C:/Users/Cuenta/anaconda3/envs/pg/Lib/site-packages/pygimli/frameworks/methodManager.py:838) #     self.fw.blockyModel = kwargs["blockyModel"]
    [839](file:///C:/Users/Cuenta/anaconda3/envs/pg/Lib/site-packages/pygimli/frameworks/methodManager.py:839) # if "stopAtChi1" in kwargs:
    [840](file:///C:/Users/Cuenta/anaconda3/envs/pg/Lib/site-packages/pygimli/frameworks/methodManager.py:840) #     self.fw._stopAtChi1 = kwargs["stopAtChi1"]
    [842](file:///C:/Users/Cuenta/anaconda3/envs/pg/Lib/site-packages/pygimli/frameworks/methodManager.py:842) self.preRun(**kwargs)
--> [843](file:///C:/Users/Cuenta/anaconda3/envs/pg/Lib/site-packages/pygimli/frameworks/methodManager.py:843) self.fw.run(dataVals, errorVals, **kwargs)
    [844](file:///C:/Users/Cuenta/anaconda3/envs/pg/Lib/site-packages/pygimli/frameworks/methodManager.py:844) self.postRun(**kwargs)
    [845](file:///C:/Users/Cuenta/anaconda3/envs/pg/Lib/site-packages/pygimli/frameworks/methodManager.py:845) return self.paraModel(self.fw.model)

File c:\Users\Cuenta\anaconda3\envs\pg\Lib\site-packages\pygimli\frameworks\inversion.py:1422, in ClassicInversion.run(self, dataVals, errorVals, **kwargs)
   [1420](file:///C:/Users/Cuenta/anaconda3/envs/pg/Lib/site-packages/pygimli/frameworks/inversion.py:1420) # self.inv.start()  # start is reset() and run() so better run?
   [1421](file:///C:/Users/Cuenta/anaconda3/envs/pg/Lib/site-packages/pygimli/frameworks/inversion.py:1421) self.inv.setMaxIter(0)
-> [1422](file:///C:/Users/Cuenta/anaconda3/envs/pg/Lib/site-packages/pygimli/frameworks/inversion.py:1422) self.inv.start()
   [1423](file:///C:/Users/Cuenta/anaconda3/envs/pg/Lib/site-packages/pygimli/frameworks/inversion.py:1423) self.maxIter = maxIterTmp
   [1424](file:///C:/Users/Cuenta/anaconda3/envs/pg/Lib/site-packages/pygimli/frameworks/inversion.py:1424) if self.verbose:

RuntimeError: ./core/src/inversion.cpp:95		double GIMLI::RInversion::getPhiD(const Vec&) const  getPhiD == inf

0 replies

halbmy · 2025-02-18T08:11:41Z

halbmy
Feb 18, 2025
Maintainer

An infinite data objective (phiD) is mostly due to zero error values (or already infinite input). You can get additional output in form of files by adding debug=True to the inversion and have a look at the data.

I didn't run the whole code but it all comes down to the inversion line. Can you please save and attach the mesh and the data container (both with mesh.save('mesh.bms') and datacontainer.save('data.dat')?

0 replies

aldefern · 2025-02-18T08:54:01Z

aldefern
Feb 18, 2025
Author

Datos.zip
Here it is. I also think it is reading the [a b m n] wrong. They should not be integer numbers

0 replies

aldefern · 2025-02-18T09:45:31Z

aldefern
Feb 18, 2025
Author

I have changed a couple of thing regarding the electrodes' positions. Here is a zip with everything (data, code, response, data_container)
Allfiles.zip

0 replies

halbmy · 2025-02-18T12:27:07Z

halbmy
Feb 18, 2025
Maintainer

Your mesh is having the outer boundary with marker 2 and the inner with marker 1 (by default it is the other way),

so you need to tell which is the inversion and which is background by

mgr = ert.ERTManager(data, verbose=True)
mgr.setMesh(mesh)
mgr.inv.setRegularization(1, background=False)
mgr.inv.setRegularization(2, background=True)
model = mgr.invert(lam=20, verbose=True)

Then everything works well technically

However, something is wrong with your data file data.dat as for a lot of data identical indices are occurring:

# a b m n ...
1	2	1	3	
2	3	1	3	
3	3	2	4	
3	4	3	5	
4	5	3	5	
5	5	4	6	
5	6	5	7	
6	7	5	8	
7	8	6	8

There is something wrong.

0 replies

halbmy · 2025-05-23T12:15:55Z

halbmy
May 23, 2025
Maintainer

Does it work?

0 replies

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

HELP with ERT Crosshole Inversión. #820

Uh oh!

{{title}}

Uh oh!

Replies: 6 comments

Uh oh!

{{title}}

Uh oh!

Uh oh!

{{editor}}'s edit

{{editor}}'s edit

Uh oh!

Uh oh!

{{title}}

Uh oh!

Uh oh!

{{title}}

Uh oh!

Uh oh!

{{title}}

Uh oh!

Uh oh!

{{title}}

Uh oh!

Uh oh!

{{title}}

Uh oh!

Select a reply

Uh oh!

HELP with ERT Crosshole Inversión. #820

Uh oh!

aldefern Feb 17, 2025

Replies: 6 comments

Uh oh!

Uh oh!

aldefern Feb 17, 2025 Author

Here is the error I get:

Uh oh!

halbmy Feb 18, 2025 Maintainer

Uh oh!

aldefern Feb 18, 2025 Author

Uh oh!

aldefern Feb 18, 2025 Author

Uh oh!

halbmy Feb 18, 2025 Maintainer

Uh oh!

halbmy May 23, 2025 Maintainer

aldefern
Feb 17, 2025

aldefern
Feb 17, 2025
Author

halbmy
Feb 18, 2025
Maintainer

aldefern
Feb 18, 2025
Author

aldefern
Feb 18, 2025
Author

halbmy
Feb 18, 2025
Maintainer

halbmy
May 23, 2025
Maintainer