Vertical Electrical Sounding

[ferryatm]

datatest.xlsx

Please help me, I'm new to this and really want to learn pygimli for VES. I have Schlumberger VES field data, can you make me a 1D VES inversion code? like starting from synthetic data and then calling up the field data that I have in Excel. I plan to compare the results with the curves in IPI2Win and Res1Dinv. Please help me just this once :(

[halbmy]

Please have a look at the tutorial https://www.pygimli.org/_tutorials_auto/3_inversion/plot_2-dc1dblock.html

I created a script

import pandas as pd
import matplotlib.pyplot as plt
import pygimli as pg
from pygimli.physics.ves import VESModelling
from pygimli.frameworks import MarquardtInversion
from pygimli.viewer.mpl import drawModel1D

data = pd.read_excel("datatest.xlsx")
fop = VESModelling(ab2=data["AB/2 (m)"],
                   mn2=data["MN/2 (m)"],
                   nLayers=2)
inv = MarquardtInversion(fop=fop, verbose=True)

rhoa = data["Rho (Ohm.m)"].values  # pandas->numpy
model = inv.run(rhoa, relativeError=0.05, verbose=True)

fig, ax = plt.subplots(ncols=2)
ax[0].loglog(rhoa, ab2, "x", label="data")
ax[0].loglog(inv.response, ab2, "-", label="response")
ax[0].invert_yaxis()
ax[0].legend()
ax[0].grid()
drawModel1D(model=model, ax=ax[1])

[ferryatm]

Why do the resistivity and depth graphs in the tutorial have two orange and blue curves? sorry, I still have a lot to learn

[ferryatm]

and can the inversion use Gauss-Newton method?

[halbmy]

Yes, the inversion uses a Gauss-Newton method by default.

[halbmy]

The blue crosses are the measured (or synthetic modelled) data and the orange line is the model response. In case of a synthetic experiment, one can plot both the synthetic (blue) along with the inverted (orange) model. A label is missing here. Please read the tutorial in detail and give feedback which part is not sufficiently explained.

[halbmy]

Well, the focus of pyGIMLi is not so much on Vertical Electrical Soundings and there seem to be little people working on it. Maybe you can start a github project and create some notebooks with models that you have in mind, and we can have a look and help you.

[ferryatm]

maybe we'll put that aside for now. returning to the code you provided, why does the output that appears, for the resulting depth, never exceed 30 meters? even I have tried 10 different trajectories and the depth is still no more than 30 meters? doesn't AB/2 = 100 indicate a depth of 100 meters? I also tried the IPI2win software and it can detect depths of up to 80 meters for 8 layers.

[halbmy]

AB/2=100m indicates a maximum depth of 30m according to the text books. But you can use a different starting model. After all, you are able to fit your data, but I guess with too many layers than necessary.

[ferryatm]

Sorry, I think I really have to learn more about this. Can you help me again to adjust the code? Maybe we should adjust the model initials or other parameters? Please help me once again. My lecturer has tried the IPI2win and Res1Dinv software, the results are as shown in the picture. However, the PyGIMLi results are always inaccurate or dissimilar. I've been trying for a long time but still. Please help me for just this 1 data. I was asked by my lecturer to compare the results of PyGIMLi with software so that in the future PyGIMLi can be recommended for geophysical processing.

Data on Sheet Dawuhan 2:

MN/2 (m) | AB/2 (m) | K | R (Ohm) | Rho (Ohm.m)

0.5 | 1.5 | 6.29 | 2.663 | 16.75027
0.5 | 2 | 11.79 | 1.425 | 16.80075
0.5 | 2.5 | 18.86 | 0.7492 | 14.129912
0.5 | 3 | 27.5 | 0.5317 | 14.62175
0.5 | 4 | 49.5 | 0.3454 | 17.0973
0.5 | 5 | 77.79 | 0.2341 | 18.210639
0.5 | 6 | 112.36 | 0.1685 | 18.93266
1 | 6 | 55 | 0.3277 | 18.93266
1 | 7 | 75.43 | 0.2435 | 19.29370252
1 | 8 | 99 | 0.1935 | 20.12281196
1 | 9 | 125.71 | 0.1529 | 20.19062751
1 | 10 | 155.57 | 0.1238 | 20.23107692
1 | 12 | 224.71 | 0.08969 | 21.17088088
1 | 15 | 352 | 0.06077 | 22.47006893
5 | 15 | 62.86 | 0.2601 | 22.47006893
5 | 20 | 117.86 | 0.1383 | 22.40156154
5 | 25 | 188.57 | 0.0898 | 23.27226286
5 | 30 | 275 | 0.06066 | 22.92581458
10 | 30 | 125.71 | 0.1391 | 22.92581458
10 | 40 | 235.71 | 0.06951 | 21.48093177
10 | 50 | 377.14 | 0.05463 | 27.01230324
10 | 60 | 550 | 0.03517 | 25.36079579
20 | 60 | 251.43 | 0.07155 | 25.36079579
20 | 70 | 353.57 | 0.05723 | 28.5256335
20 | 80 | 471.43 | 0.04651 | 30.91004939
20 | 90 | 605 | 0.03933 | 33.54404752
20 | 100 | 754.29 | 0.03496 | 37.17456691
30 | 100 | 476.67 | 0.05171 | 37.17456691
30 | 120 | 707.14 | 0.04037 | 43.05441709
30 | 150 | 1131.43 | 0.03059 | 52.19882499

[ferryatm]

Datasheet Plotting.xlsx

Or maybe you can use Excel. Please help me once again in improving the code so that it is at least similar to the software. Try for sheet 'Dawuhan 2'

[halbmy]

Maybe this issue will help you: #673

[ferryatm]

just the same

[ferryatm]

Can we set the initial or start model?

[halbmy]

Can you please explain your problem in detail:

1. Are you technically able to invert your sounding curve or are there problems reading the data etc.?
2. Can you fit the data appropriately and if so how many layers do you need for it?
   For your sounding curve, do do not need more than five layers and I am pretty sure that the last layer at 100m depth is beyond what you can resolve with a maximum AB/2 of 150m.

[ferryatm]

Maybe you can see my Google Colab. The resulting output cannot detect up to the maximum depth and only reaches 30 m. I want to adapt the results obtained from my lecturer in this Madiun city government project report. It can reach 7 to 9 layers with a depth of up to 98 m. You can see the report in this pdf.

intip.in/GoogleColabTA

Laporan madiun 24.pdf

[ferryatm]

This is for my thesis btw. Sorry if I seem to ask a lot of questions, because I was previously interested in PyGIMLi and wanted to try it for VES processing

[halbmy]

Again, every geophysics text book will with a maximum AB/2 of 150m you can NOT resolve a layer at a depth of 100m accordingly. Like in the example in Figure 4 of your thesis

you see that the sounding curve is saturating but there is no indication of a decrease. So your model is one of the equivalent models, but there is no substantation of the last two layers. According to Occams razor, the simplest model that explains the data within error bounds should be chosen and this is a three-layer case, if you remove the first data point, even a two-layer case.

[ferryatm]

OK, I'll try again. What about my Google Colab? Is there something wrong? And also can the start model be defined or initiated? or can it only be set by default?

[halbmy]

I cannot access the colab notebook

[ferryatm]

https://colab.research.google.com/drive/1x4s849gFNAzhIaUmtaPcM12kl0MF0syV?usp=sharing

…

________________________________ Dari: Thomas Günther ***@***.***> Dikirim: Kamis, 13 Juni 2024 16.26 Kepada: gimli-org/gimli ***@***.***> Cc: FERRY ATH THAARIQ MUDHOFIR ***@***.***>; Author ***@***.***> Subjek: Re: [gimli-org/gimli] Vertical Electrical Sounding (Discussion #708) I cannot access the colab notebook — Reply to this email directly, view it on GitHub<#708 (reply in thread)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AYMQCYRNP4PDAJUP6SO3OS3ZHFQSXAVCNFSM6AAAAABHDGYB4WVHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM4TONRRGIZDI>. You are receiving this because you authored the thread.Message ID: ***@***.***>

[halbmy]

The collection of different types of synthetic models looks really good. Would be good to have it openly visible somewhere, e.g. a github repository that renders jupyter notebooks.

Nevertheless the number of layers of the last model is too large to be resolved.

[ferryatm]

thank you, but my question has not been answered, regarding can the initial or start model be defined or initiated? or can it only be set by default? If so, can you provide an example?

[ferryatm]

I use this but there is an error ''MarquardtInversion' object has no attribute 'setModel''

Extract data from the dataframe

ab2 = data["AB/2 (m)"].values
mn2 = data["MN/2 (m)"].values

Set up the forward operator for VES

nLayers = 7
fop = VESModelling(ab2, mn2, nLayers=nLayers)
inv = MarquardtInversion(fop=fop, verbose=True)

Extract resistivity data

rhoa = data["Rho (Ohm.m)"].values

Set initial model (example values, modify as needed)

initial_model = np.array([100, 80, 60, 50, 40, 30, 20])

Set the initial model for the inversion

inv.setModel(initial_model)

Run inversion

model = inv.run(rhoa, relativeError=0.05, lam=10, maxIter=100, verbose=True)

Plot results

fig, ax = plt.subplots(ncols=2, figsize=(12, 6))

Data and response

ax[0].loglog(ab2, rhoa, "x", label="Data", markersize=8)
ax[0].loglog(ab2, inv.response, "-", label="Response")
ax[0].legend()
ax[0].grid()
ax[0].set_xlabel("AB/2 (m)")
ax[0].set_ylabel("Apparent Resistivity (Ohm.m)")
ax[0].set_title("Field Data and Model Response")

Draw the 1D model

drawModel1D(model=model, ax=ax[1])
ax[1].set_xlabel("Resistivity (Ohm.m)")
ax[1].set_ylabel("Depth (m)")
ax[1].set_title("1D Resistivity Model")
ax[1].grid(True)

plt.tight_layout()
plt.show()

[ferryatm]

how to calculate the uncertainty value of the parameters of the Levenberg-Marquardt algorithm? can you provide a code example?

[halbmy]

see #746