Welcome to gimli Discussions!

[halbmy]

👋 Welcome!

We’re using Discussions as a place to connect with other members of our community without opening issues. We hope that you:

• Ask questions you’re wondering about.
• Share ideas.
• Engage with other community members.
• Welcome others and are open-minded. Remember that this is a community we
  build together 💪.

To get started, comment below with an introduction of yourself and tell us about what you do with this community.

[Mouctarmmc]

Hi
What if we try to input real data and test the model. From a csv file with x positions, depth, rho, IP, and sp

[halbmy]

I don't understand. You have real data. And you have a model? But the model consists of points? Note that a model is a discretized (e.g. by triangles, quadrangles, tetrahedra, hexahedra etc.) subsurface representation with piece-wise constant values. Please have a look at the modelling examples.

[Mouctarmmc]

How to proceed

[halbmy]

Please explain in a bit more detail so that we can understand your problem.

[Mouctarmmc]

I have a .dat file of resistivity and ip values of schlumberger array. How can I insert into the model and appreciate the results. NB :it's a tomography

[halbmy]

If your model is already in a mesh with piece-wise constant resistivity you have to create it, otherwise create a new one and make interpolation from your point data. Please have a look at the examples and tutorials on mesh generation, interpolation and ERT simulation. Don't hesitate to come back with a specific question (e.g. with code parts, images, files etc.).

[LuFrag]

Hello,

I am trying to use PyGimli to invert a 1-D GPR waveform for a research project, I am a real beginner.

I already have a FDTD forward operator, but I am struggling to set up the inversion and the start model using regions.

Thanks in advance for all your help !

First of all, I want to try a blocky inversion considering a two-layer model (with absorbing boundary conditions at the end). Each layer is only charaterized by its dielectric permittivity (eps) for now. So, I have two thicknesses and two parameters to invert.

class myFwd_Block(pg.Modelling):
    def __init__(self, n_layers):
        """Initialize the model."""
        self.n_lay = n_layers
        self.mesh1d = pg.meshtools.createMesh1DBlock(self.n_lay)
        super().__init__()
        self.setMesh(self.mesh1d)
    
    def response(self, model):
        """Forward response."""
        A = fwd_block(model,self.n_lay)
        return A
        
    def createStartModel(self,data):
        
        th = pg.Vector(self.n_lay,1)
        ep = pg.Vector(self.n_lay,12)
        stmo = pg.cat(th,ep)
        return stmo

I am using a createStartModel function because I did not know how to set up regions like thicknesses and eps. The I set up the forward model and the inversion framework using a model vector composed of my thicknesses and eps.

# Setting the forward model into the inversion scheme
thk =np.ones(n_layers)
fop = myFwd_Block(n_layers)
epsilon = np.ones(n_layers) * 4
sig = np.ones(n_layers) * 2e-4
model_base = pg.cat(thk,epsilon)
response = fop.response(model_base)
plt.plot(response)

inv = pg.Inversion()
inv.setForwardOperator(fop)
inv._inv.setMarquardtScheme(0.5)

# inv._inv.setLambda(100)
model = inv.run(data, relativeError,\
                verbose=True,maxIter=20,lam=10,debug=True)

After setting my forward model, I launch the inversion. The start model is correctly launched I guess. But then the number of parameters (4 at the begininng) goes to 3 during the first iteration ! An then it crashes.

1. Is it because the last layer is considered semi-infinite inside PyGimli framework ? Maybe I don't need this assumption because of the absorbing boundary conditions of my model ?
2. Is it possible to better set up my starting model using the RegionManager in myFwd_Block function ? Could you please give me some hint on how to do that ?
3. If I would like to invert the permittivity and conductivity, is there a lot of modifications to do in the workflow and in the forward operator ?

30/05/22 - 15:42:34 - pyGIMLi - INFO - Starting inversion.
30/05/22 - 15:42:34 - pyGIMLi - INFO - Created startmodel from forward operator: 4 [1.0, 1.0, 12.0, 12.0]
fop: <__main__.myFwd_Block object at 0x000001C6027F0540>
Data transformation: <pgcore._pygimli_.RTransLin object at 0x000001C603558220>
Model transformation: <pgcore._pygimli_.RTransLog object at 0x000001C6031AB680>
min/max (data): -0.23/0.50
min/max (error): -5.3e+04%/1.2e+06%
min/max (start model): 1.00/12.00
min/max(dweight) = 1000/1000
min/max(dweight) = 1000/1000

Building constraints matrix
constraint matrix of size(nBounds x nModel) 3 x 3
C:/msys64/home/halbm/gimli/gimli/core/src/inversion.cpp:43		void GIMLI::RInversion::checkConstraints()  resize model 4 to fit constrain size: 3
check Jacobian: wrong dimensions: (0x0) should be (1000x3)  force: 1
jacobian size invalid, forced recalc: 1

[halbmy]

1. Yes, in a block model inversion (of both thickness and values) the last layer is always infinite, so you need nlay-1 thickness values.

2. You can also pass a starting model directly to the inv.run(startModel= call

3. No, you specify two parameters by pg.meshtools.createMesh1DBlock(self.n_lay, 2) and modifiy the starting model accordingly

[LuFrag]

Hello again,

Thanks for your quick answer 👍

I modified my forward operator to add a semi-infinite layer at the end of my 1-D model. Moreover the start model is also OK.

The inversion using a Marquardt-Levemberg scheme runs : the results below are from an inversion from the two-layer model.

The inversion is highly sensible to the starting model (not a suprise !). I would like to define boundaries for each parameter of the inversion, but it seems that the pg.trans.Log, pg.transLogLU commands are not working as in the example showed on Transfrom2021.

inv = pg.Inversion()
inv.setForwardOperator(fop)
transModel = pg.trans.TransLog(1) # > 1 epsilon (and what about the thickness ?)
inv.transModel = transModel
inv._inv.setMarquardtScheme(0.9)

lamb = [250,300]

for l in lamb:
    model = inv.run(data, relativeError,\
                    verbose=True,maxIter=10,lam=l)

1. Is there a way to verify the implementation of these transformations ? Or am I doing something wrong ?

2. Is there an user-defined way to impose differents bounds and variations on thicknesses and other parameters ? Do I have to use Block1DInversion class ?

3. For the classical Inversion framework (Occam inversion), is there also a semi-infinite layer at the bottom of the model ?

Here is some results of the Occam inversion using a model with 20 layers with thickness of 20 cm.

4. Is the runchi1() method still available to find the better regularization parameter ? I am having problems because the inversion is not reaching chi2=1 for most of the runs, it only quits because dPhi < 2% or even goes to a negative dPhi.

5. Is there a way to verify that the data and model transformation are well aplied in this case ? I am only trying :

 inv = pg.Inversion()
inv.setForwardOperator(fop)
#inv.inv.setLambda(100)
transModel = pg.trans.TransLog(1) # > 1 epsilon
inv.transModel = transModel

damps = np.logspace(1,3,10)

for l in damps:
    model = inv.run(data, relativeError, startModel=12,verbose=True,maxIter=50,lam=l)

Thanks again !

[halbmy]

Thickness and epsilon values should of course be treated differently, with different limits etc. So the best way is a region-specific behaviour. The block model generation automatically creates nProperties+1 regions.

>>> mesh=pg.meshtools.createMesh1DBlock(2)  # 2 layers, 1 property
>>> print(mesh.cellMarkers())
3 [0, 1, 1]
>>> mesh=pg.meshtools.createMesh1DBlock(3, 2)  # 3 layers, 2 properties
>>> print(mesh.cellMarkers())
8 [0, 0, 1, 1, 1, 2, 2, 2]

For region-specific bounds, you should have a look at the 2022 Transform tutorials (https://github.com/gimli-org/transform2022).

inv.setRegularization(0, limits=[0, 10])  # thickness
inv.setRegularization(1, limits=[3, 100])  # epsilon

3. In the Occam inversion there is usually also a semi-infinite layer at the end but it is up to your forward operator how this is done internally.
4. For the Marquardt inversion (where lambda is successively decreased), the inv.runChi1() method does not make much sense. With a large enough value it should converge. You might setup the lambda factor or the deltaPhi stopping criterion.

[LuFrag]

Hello,

Thanks again for your answers and help.

The region specific bounds worked nicely for the Marquadt-Levemberg scheme! As you suggested, I just fine-tuned the lambda factor (using the trace of the J.T*J matrix) and let deltaPhi as default.

inv = pg.Inversion(fop=fop)
inv.setRegularization(0, limits = [1,2]) # Limits for region (0): thicknesses
inv.setRegularization(1, limits = [1,80]) # Limits for region (1) : epsr

inv._inv.setMarquardtScheme(0.9)

model = inv.run(data, relativeError,\
                verbose=True,maxIter=50,lam=500)

thk_inv,eps_inv = np.array(model[0:n_layers-1]),np.array(model[n_layers-1:])

# Retrieving inversion information verbose
model_h = np.array(inv.modelHistory)
chi2_h = np.array(inv.chi2History)

# Choosing starting damping factor
J = pg.utils.gmat2numpy(fop.jacobian())

trace_J = np.trace(np.matmul(J.T,J))

Departing from a starting model with [thk1 = 0.5 m, epsr1 = epsr2 = 20], the inversion converges to the simulated model [thk1 = 1 m, epsr1 = 4, epsr2 = 16].

Concerning the classical inversion framework (smooth, Occam-Type), I saw the example for a smooth VES inversion (https://www.pygimli.org/_tutorials_auto/3_inversion/plot_4-dc1dsmooth.html). The runChi1() method seems to work for the class pg.core.Inversion but it is not recognized by the class pg.frameworks.inversion.

    model =inv.runChi1()

AttributeError: 'Inversion' object has no attribute 'runChi1'

1. Is it possible to use runChi1() ? How can I correctly set up the pg.Inversion class to be able to run this method ? I am surely doing something wrong because I don't even call my data, error vectors...

I tried to make a L-curve using phiData and phiModel to figure out the regularization factor from the elbow of the curve, but it seems that my data misift is not improving with the range of regularization factors chosen (1e-9 to 1e+9).

2. Maybe I have to try another approach or even redefine my cost function ?

Thanks !

[halbmy]

The current pg.Inversion does not support runChi1(), you would have to use the underlying core inversion object. However:

The runChi1 method step by step decreases or increases the regularization parameter until the data are fitted within noise, using a global regularization scheme with a fixed objective function. For a Marquardt type inversion (where the regularization is only local and the regularization parameter is varied), this does NOT make any sense as mentioned above. The same holds for the L-curve.

[LuFrag]

Thanks for your answer. Now the inversion works.

I have a new question.

I am working now on some ERT data acquired using reciprocal werner schlumberger mode (WSR). When I plot the apparent resistivity section, the quadripoles are identified as Wenner-Alpha type (WA). Is it normal because I am using a unconventional scheme ?

How can I be sure that the a,b,m,n electrodes are correctly taken in into accound ? IN WSR, the m, n are the outer electrodes and a,b are the inner electrodes, the contrary of wenner-schlumberger.

You will find the data attached.

ERT48_gapout.zip

Thnaks again :)

[halbmy]

Yes, that's normal. Actually we could name it WS (Wenner-Schlumberger) if it not a classical Wenner but this is just a label. We could also exchange the inner and outer spacings in the label, but I don't see a strong need as it is just a label. The pseudosection is displayed correctly.

To check whether all is correct, you can save the data using data.save("out.dat", "a b m n rhoa") and check whether the indices point to the correct positions. Or use

ex = pg.x(data)
np.column_stack([ex[data[tok]] for tok in 'abmn'])

[LuFrag]

Thanks for your answer, I checked the electrodes and everything's fine.

I was wondering how the regularization constraints type of PyGimli relate to the smoothness-contraints (blocky and smooth) applied on Res2DInv ? Maybe it is a silly question, but I will try my chance anyway.

I found the topic "Change regularization for crosshole traveltime tomograpy #355", but it is not answered yet. For a blocky inversion, you propose to use the option blockyModel=True.

I tried something like below, but the inversion results does not change from the default options of the ERTManager.invert():

data = pg.physics.ert.load(filename)

ert = pg.physics.ert.ERTManager()

mod = ert.invert(data,mesh=mesh,lam=5,verbose=True,blockyModel=True)

Even if I try the code below, the results still the same as the default inversion parameters :

ert.inv.setRegularization(cType=2,zWeight=0.2)

mod = ert.invert(data,mesh=mesh,lam=5,verbose=True)

I am missing something, but what ? Maybe the regularization effects are not visible because I have a lot of apparent resistivity data points to constraint the inversion ?

Let me know if you need my data file.

Thanks again for your help,
Luis

[ZSSSSZJU]

Hi！
I'm trying to use pyGIMLI to study the concentration distribution of underground pollutants through complex resistivity inversion imaging
Since I have just started my research, I hope you can forgive me for asking questions that may seem naive.
question desription：
In order to solve the resistivity distribution at a certain depth, I try to use vertical electrode distribution，
When using function：

velec=[[15,i] for i in range (-20,1)]
scheme = pygimli.physics.createData(elecs=velec,schemeName='dd')

It seems correct, but in the function:showERTdata,An error has occurred：
Something gone wrong while drawing data. Try fallback with equidistant electrodes
However, from the final inversion results, it seems to basically conform to the preset model（）
I wonder why the program reports an error and how to debug it
i truly appreciate your asistance and look forward to you reply!

[halbmy]

The message is not an error but just a warning. showERTdata is designed for surface ERT layouts and it detects that this is obviously none and uses electrode numbers instead of horizontal distances. So all is right.

[ZSSSSZJU]

Sincerely thank you for your quick reply and pay tribute to your masterpiece pyGIMLI

[ZSSSSZJU]

Hi！Thank you again for your patient answer
now i Now I have a new problem，although it may seem too easy for you.
problem discription:
Generally speaking, most problems can be solved by arranging electrodes on the surface，But as stated in my last question, I want to use two rows of electrodes arranged vertically to realize an inversion，just like the image below:

It seems that joint inversion using the data of two rows of electrodes is a simple and direct method, which involves the assembly of Jacobian matrix and others. Is there any simple and direct instruction in pyGIMLI? Or I need to complete this process manually.

[halbmy]

Joint inversion is indeed a simple task for pyGIMLi but here you don't need it. You just need to create your data file holding all the electrodes and appropriate quadrupoles. The latter is of course a bit complicating, but there are several papers on designing crosshole ERT arrays.

[ZSSSSZJU]

Well, in fact, the problem may be more complex. The two rows of electrodes are in the form of dipole dipole, and both are complex resistivity. So at first, I thought it would be great if there were a simpler way to perform joint inversion by assembling the Jacobian matrix. What functions need to be called after arranging the data files? Can complex resistivity also be realized?

[halbmy]

The two rows of electrodes are in the form of dipole dipole, and both are complex resistivity.

I can't get any meaning of this sentence. The electrode positions have nothing to do with the array chosen and of course both are independent on the resistivity distribution. Crosshole ERT has nothing to do with joint inversion.

As the examples demonstrate, of course you can use complex resistivity, but you should start with real-valued first. In case you don't mean IP but complex structures, this is of course all up to you.

[ZSSSSZJU]

Please forgive my mistake. I want to express IP data. Now that I have solved the problem, I have a small question about geometric factors:

I noticed that in the calculation of apparent resistivity, the geometric factor is usually assumed to be isotropic and homogeneous

My question is:

When measuring the apparent resistivity, is any electrode distribution applicable to this geometric factor calculation? For example, do you need to consider the difference of boundary conditions when calculating the geometric factor between electrodes distributed on the surface and those distributed underground?

Thank you for your patience every time!

[halbmy]

The geometric factor is based on homogeneous isotropic resistivity. In case of trivial topography (unbounded domain, no topography) it is computed analytically, otherwise numerically. Generally it can be computed for any type of electrode array, but there might be arrays (e.g. null arrays associated with infinite geometric factors) where the concept does not make sense.

[LuFrag]

Hello everyone,

I am trying to invert traveltime refraction seismics data and ERT data using the same mesh. This is not (yet) a joint inversion, I just want to use the same unstructured mesh to invert data separately on a conventional way.

1. Do the electrodes and geophones must be collocated for the inversion to better perform in each case ? In my case, the shots were done at the electrode locations.

Here is my protocol to create the common mesh:

From a open topography polyline (created with createPolygon, IsClosed=False), I used the function appendTriangleBoundary from meshtools to create a first mesh with only one region. This mesh was used for the traveltime inversion. Something like the example in the API documentation :

surf = mt.createPolygon([[0, 0],[5, 3], [10, 1]], boundaryMarker=-1, addNodes=5, interpolate='spline')
m = mt.appendTriangleBoundary(surf, xbound=4, ybound=4, marker=2, addNodes=5, isSubSurface=True)
ax, _ = pg.show(m, markers=True, showBoundaries=True, showMesh=True, ax=axs[1][2])

Then, I used the function appendTriangleBoundary again to create a « far away » zone on the first mesh. This new composed mesh was used for the ERT inversion, since I can retrieve the inverted ERT data on the same region (background excluded) as the traveltime inverted velocities.

2. Is it the best way to assure that both inversion were done on the same mesh ? How does the mesh have to be related to geophones or electrode positions ? The sensors must be at the nodes ?

I tried to invert the travel time data using the ERT mesh, but it does not seem to recognize the boundaries or the regions. Here is the error I get :

RuntimeError                              Traceback (most recent call last)
~\AppData\Local\Temp\ipykernel_7188\15074670.py in <cell line: 4>()
      2 #vel = donnees_sismique_refrac.invert(maxIter=8,verbose=True)
      3 #vel = donnees_sismique_refrac.invert(mesh=mesh,maxIter=8,verbose=True,lam=840)
----> 4 vel = donnees_sismique_refrac.invert(mesh=mesh,maxIter=8,verbose=True)

~\Anaconda3\envs\pg\lib\site-packages\pygimli\physics\traveltime\TravelTimeManager.py in invert(self, data, useGradient, vTop, vBottom, secNodes, **kwargs)
    239             self.fop._useGradient = None
    240 
--> 241         slowness = super().invert(data, mesh, **kwargs)
    242         velocity = 1.0 / slowness
    243         self.fw.model = velocity

~\Anaconda3\envs\pg\lib\site-packages\pygimli\frameworks\methodManager.py in invert(self, data, mesh, startModel, **kwargs)
    757 
    758         self.preRun(**kwargs)
--> 759         self.fw.run(dataVals, errorVals, **kwargs)
    760         self.postRun(**kwargs)
    761         return self.paraModel(self.fw.model)

~\Anaconda3\envs\pg\lib\site-packages\pygimli\frameworks\inversion.py in run(self, dataVals, errorVals, **kwargs)
    551             self._preStep(0, self)
    552 
--> 553         self.inv.start()
    554         self.maxIter = maxIterTmp
    555 

~\Anaconda3\envs\pg\lib\site-packages\pygimli\physics\traveltime\modelling.py in createJacobian(self, par)
     82             pg.critical("no mesh")
     83 
---> 84         return self._core.createJacobian(par)
     85 
     86     def response(self, par):

RuntimeError: C:/msys64/home/halbm/gimli/gimli/core/src/sparsematrix.h:342		GIMLI::SparseMapMatrix<ValueType, IndexType>::MatElement GIMLI::SparseMapMatrix<ValueType, IndexType>::Aux::operator[](IndexType) [with ValueType = double; IndexType = long long unsigned int; GIMLI::SparseMapMatrix<ValueType, IndexType>::MatElement = GIMLI::MatrixElement<double, long long unsigned int, std::map<std::pair<long long unsigned int, long long unsigned int>, double, std::less<std::pair<long long unsigned int, long long unsigned int> >, std::allocator<std::pair<const std::pair<long long unsigned int, long long unsigned int>, double> > > >]  idx = 18446744073709551615, 11 maxcol = 2961 stype: 0

Thanks again for your help and support, let me know if you need more information,

Cheers,
Luis

[ZSSSSZJU]

Hello again！
I want to ask a simple question, as you can see, im trying to use datacontainer to read information from the file of sensor position and quadrupoles-measurement

and all steps are fine (Input and output of sensors and quadrupoles) until the file is saved, The problem lies in the output of quadrupoles, The following is a brief code

scheme=pg.DataContainerERT()
for i in range(sensor_num+2,sensor_num+scheme_num+2,1):
    a=np.append(a,int(data[i].split()[0]))
    b=np.append(b,int(data[i].split()[1]))
    m=np.append(m,int(data[i].split()[2]))
    n=np.append(n,int(data[i].split()[3]))
    
scheme["a"]=a
scheme["b"]=b
scheme["m"]=m
scheme["n"]=n

now i print the scheme and the quadrupoles, the output is alright

print(scheme)
print(scheme["a"])
print(scheme["b"])
print(scheme["m"])
print(scheme["n"])

output:

Data: Sensors: 10 data: 7, nonzero entries: ['a', 'b', 'm', 'n']
[1 2 3 4 5 6 7]
[2 3 4 5 6 7 8]
[3 4 5 6 7 8 9]
[ 4  5  6  7  8  9 10]

but after i saved the file,it turns out that the data in 'a', 'b','m','n' are all different with the number in the original file(1 larger)
scheme.save('schemedemo.txt')

I'm sorry to bother you with such a simple question，but it did make me puzzled. Thank you again for your patience and pay tribute to your great work“pygmili”

[halbmy]

You're not the only one to be puzzled.

The unified data format was created (for ERT) long time ago when we have been working with Matlab instead of Python. Therefore we counted from one in the file (and 0 meaning infinite electrode). However, in Python it makes more sense to count from zero. But when creating data containers in Python you will have to count from zero but it will be saved to file as one.

[ZSSSSZJU]

Thank you for your answer. It seems that there is no problem with the data in the datacontainer. Just need to modify it when save it. Thank you very much

[ZSSSSZJU]

Hello again!
I have been reading the 'naive IP data inversion' in the example carefully these days，There are two little issues that I don't quite understand.
Q1：
in the example(line268-274):

dmag = np.abs(data_rcomplex)
dpha = np.arctan2(data_rcomplex.imag, data_rcomplex.real) * 1000
fig, axes = plt.subplots(1, 2, figsize=(20 / 2.54, 10 / 2.54))
k = np.array(ert.createGeometricFactors(scheme))
ert.showERTData(
    scheme, vals=dmag * k, ax=axes[0], label=r'$|\rho_a|~[\Omega$m]')

it seems to mean the 'dmag' is the impedance we measured （Unit: Ohms）however, the 'dmag' is extracted from "data['rho']" simulated by ert.simulate and in ERT , data['rho'] usually means resistivity （Unit: Ohms*m）it makes me confused
Q2
after i finished complex resistivity modeling and saved the simulated data, I found that compared to ERT modeling, CR modeling seems to generate a new array named ['iperr'] , And its value is equivalent to one percent of the phase value (without error), The following are the codes I use when simulating：
data = ert.simulate(mesh, scheme=scheme, res=rhomap, noiseLevel=1, noiseAbs=1e-5)
Finally, I would like to ask what is the role of argument"seed=XXX" in ert.simulate, i didn't find it in the GIIMLI documentation
I believe these are very simple questions for you. Thank you for your patient answers every time！. Good day！

[halbmy]

Q1: data["rhoa"] means apparent resistivity and this is where dmag comes from. There is no data["rho"]. However, we set the geometric factors deliberately to 1m so that it is actually the impedance. This is just a scaling.

Q2: to describe the IP error we use the token iperr which is default 1% but you can pass a phase error by ert.simulate(..., phiErr=...

Q3: seed sets the random generator to a specific state so that the noise is always the same. See numpy.rand API doc.

[ZSSSSZJU]

Dear Professor！First of all, I would like to express the apology for my unclear expression of data['rhoa'] last time
Secondly, I have a small problem that I'm not very clear about：
Problem Description ：When I tried to perform CRmodeling (the synthetic model is as follows),

rhomap = [
[1, pg.utils.complex.toComplex(100, 0 / 1000)]
[2, pg.utils.complex.toComplex(400, 0 / 1000)],
[3, pg.utils.complex.toComplex(200, 0 / 1000)],
[4, pg.utils.complex.toComplex(100, -300 / 1000)]
]

It's very easy to get apparent complex resistivity by ert.simulate().


However, some of the phia are positive, Even if I remove the error, it's still the case. which means that the waveform of the potential is more advanced than the waveform of the current and behaves more like inductive reactance. This is not common in my previous laboratory SIP research. and as you can see, I did not set an region with a phi larger than 0. So I would like to know why a positive value occurs.
Sorry for my frequent questions, but your patient answers have been a great help to me. Thank you very much

[halbmy]

The intrinsic phase (phi) is always negative (resistivity phase, or positive for conductivity phase) for physical reasons, and this is what you measure in the lab as long as the probe is rather homogeneous.

In contrast, the apparent phase can be positive as a result of negative sensitivities, which are the weighting factors. Typically, it is reported that such an effect occurs for chargeable material in good conductors (typically shallow layers) but I can well image that for dipole-dipole data it can also occur for anomalies deeper than the penetration depth as the sensitivity is negative there, too:
https://github.com/gimli-org/notebooks/blob/main/ERT/arrays/ERTArrays.ipynb

Dahlin, T. & Loke, M.H. (2015): Negative apparent chargeability in time-domain induced polarisation data. J. Appl. Geoph. 123, http://dx.doi.org/10.1016/j.jappgeo.2015.08.012

[ZSSSSZJU]

Dear Professor！
First of all, I want to express my gratitude to you for your meticulous answer last time about the positive phia. After receiving your reply, I carefully reviewed the relevant documents, I understand that sensitivity distribution can cause phia to be positive (when the sensitivity of a negative-phi region is negative).
Recently, I have spent a lot of time studying example of naive CRtomo. I find it really convenient and intuitive to implement CRTomo (kemna 2000) with pygimli. However, I still have some confusion and would like to ask you questions. I divided the question into two parts, one is about questions in the demonstration, and the other is about problems in the code I created myself.
Q1
in this demonstration, geometric factors are set to 1, and potential electrodes are switched:
m = scheme['m'], n = scheme['n'], scheme['m'] = n, scheme['n'] = m
scheme.set('k', [1 for x in range(scheme.size())])
As I asked before, this is equivalent to creating the impedance which is used in inversion. Is this because the magnitude difference between rho and phi is too large due to a lagre geometric factor, resulting in distorted inversion results. In my own code, I set geometric factor to - 1 to avoid switching potential electrodes. scheme.set('k', [-1 for x in range(scheme.size())])
Q2
I noticed that when tranforming jacobian to log-log sensitivities by chain rule
J = J0 / np.exp(response[:, np.newaxis]) * np.exp(m_old)[np.newaxis, :]
It seems that J is composed of the real and imaginary parts of resistivity(ρ’ and ρ‘’) because both response and m_old are in the form of ln(magnitude)+j*phase, So will there be a mathematical error when calculating the delta_m? however the result seems to be fine.
term1 = J.conj().T.dot(WdTwd).dot(J) + lam * Wm.T.dot(Wm)
term2 = J.conj().T.dot(WdTwd).dot(d - response) - lam * Wm.T.dot(Wm).dot(m_old)
Q3
in this demonstration, No noise added to the forward data. data = ert.simulate(mesh, res=rhomap, scheme=scheme,) , Then we added noise to the forward data: d_rlog = np.log(np.exp(d_rlog.real) + noise_magnitude) + 1j * (d_rlog.imag + noise_phase), But during inversion, it still uses d = np.log(data_rcomplex) instead of d_rlog, and the error matrix Wd only contains magnitude error.it makes me puzzled.
The above are some of my problems with the code in the example. Next, I will expound some of the problems encountered in the code I created myself
I created two python files, one responsible for creating a synthetic model simulating underground contaminant plumes and generating a txt file(field_data.txt) containing sensor positions, measurement scheme, impedance and phase that i simulated (k is set to -1). Another python file is responsible for reading the txt file and inverting it. The synthesis model and field_data.txt are as follows：

[1, pg.utils.complex.toComplex(100, 0 / 1000)],
[2, pg.utils.complex.toComplex(200, 0 / 1000)],
[3, pg.utils.complex.toComplex(300, 0 / 1000)],
[4, pg.utils.complex.toComplex(50, -300 / 1000)]
field_data.txt
The mesh for inversion is as follows：

Before inversion, I output the forward data read from the txt file, and the result indicates that the reading is normal and there are no errors：

Inversion result (just 1 iteration, lam=100) and the misfit are as follows：


So I tried a second iteration, The code is as follows

**My questions are: **
Q4
Is my code for the second iteration correct? The reason I ask you this question is not because I am too lazy to check my code. In fact, during the first iteration, chi2 (calculated by chi_2=((d - response).conj().T.dot(WdTwd).dot(d - response)).real/len(d)) had dropped from 73.914 to 0.687, which surprised me with such a fast convergence rate. As a result, the model update was very small during the second iteration, so I was unable to determine whether the code I wrote was executing the next iteration, or just repeating the first step due to some parameter errors.
My questions are all above
Such a lengthy question may cause some inconvenience to you, but I did spend a lot of time studying and ’trial and error‘ to solve the problem myself. Thank you very much for your patience and care all the time, and I would like to extend my most sincere gratitude to you!

[halbmy]

Q1: there should be no difference whether you use k=-1 or k=1 and switch M and N

Q2: No, I think that's correct. The (real-valued) model vector consist of real and imaginary parts but that's implemented.

Q3: Yes, you're right. We should use the noisified data.

Q4: Not sure and no time for a deep look.

[ZSSSSZJU]

Dear professor！
Thank you very much for taking the time to respond last time！
Regarding the last question I asked last time, I have tried many times since then. Finally, I found that the reason for low-chi2 was that I overestimated the error during inversion. It seems that error estimation plays a very important role in inversion.
Q1: so when carrying the real world field data inversion(just like the example"field dara ERT"), How to correctly estimate errors without knowing the artificial manual noise
After using noisified data, I rewrote the error estimation section, Finally, after three iterations, the model update tends to converge, with chi-2 decreasing from 87225 to 2.96, and the final result's RRMS being 0.0389：

Where chi_2 and RRMS are calculated by the following formula
chi_2_=((d - response).conj().T.dot(WdTwd).dot(d - response)).real/len(d)
RRMS=np.sqrt((((d-response)/d).conj().T.dot((d-response)/d))/len(d))
Q2: From the perspective of chi-square error and RRMS, whether the inversion result is correct enough
Q3: What should I do to eliminate (or optimize?) some abnormal regions (black rectangles) in the inversion results, I know that increasing dipole lengths can help to improve resolution for deeper model, However, these abnormal regions are relatively shallow

Thank you for your patient guidance and answers! Best wishes！

[halbmy]

For synthetic data, you should know your error and fit your model accordingly.

For field data, you should have a look at the misfit distribution, e.g. as done in
https://www.pygimli.org/_tutorials_auto/3_inversion/plot_8-regionWise.html
If this shown uncorrelated noise with no systematics, your data are perfectly fitted. If there is systematics, there could be still structures that are not yet resolved. But it could also be 3D effects that are not explained by a 2D model.
For additional improvements you might test all options for the inversion in general or specific regions.

[ZSSSSZJU]

Dear professor, I suddenly came up with a question about ERT in the 2D domain:

In the two examples of '2D crosshole ERT inversion' and the example 'ERT field data with topography', there is no setting on whether there is a topography effect. At the same time, in the input files (in class of Datacontainer) of these two examples, the coordinates of the sensors are both in the form of (x, 0, z). So, how does pygimli determine whether the input file is a data form crosshole or data with topography

[ZSSSSZJU]

I tried again on the above issue and found that if you input a txt document (Fig below) with sensor coordinates in the form of (x, 0, y) like the example '2D crosshole ERT inversion', when performing ert. simulate (), an error message will appear: Warning! Swap YZ coordinates for sensor positions to meet mesh dimensions.
Does this mean that during simulation, if it is a Crosshole measurement method, the sensor coordinates must be filled in as (x, y, 0). In addition, I noticed that the example '2D crosshole ERT inversion' involves meshing before inversion, including the setting of geometric dimensions and boundary conditions. In the example 'ERT field data with topography', this step is not present. Does this mean that when we read an input file with sensor coordinates in the form of (x, 0, z) and do mod=mgr.invert (data, lam=10, verbose=True) without providing a correct grid, it will be considered by the program to be an inversion with topography, resulting in an error

[halbmy]

The z coordinate could indeed be a topography or a (negative) depth below the surface (including topography). The decision is made on the mesh, whether the upper surface is constant zero (flat earth, electrodes possibly in the subsurface) or not (topography).

A 2D inversion is done on a 2D mesh containing two dimensions, i.e. x and y positions. If the data file has y=0 and z not, the dimensions are swapped so that it matches the mesh. But this is done automatically.

[ZSSSSZJU]

Dear Professor, I have two questions regarding 2D-ERT that I would like to discuss with you.
Q1
Based on my understanding, the parameter zweight is used to control the smoothness in the depth direction during inversion. When zweight is reduced, the inversion results tend to show more distinct layering structures. This is somewhat similar to αx/αz in CRTomo (Kemna, 2000). I wonder if there is any specific mathematical expression expression for zweight in pyGIMLI?
Q2
During the imaging process, the parameter coverage provides an approximate confidence region for the inversion results by evaluating the sensitivity. However, probably due to my lack of skills, I did not find the pg.core.coverageDCtrans ( ) function in the function coverage ( ), which looks like a logarithmic sensitivity calculation based on the chain rule. I wonder what is the specific calculation process of pg.core.coverageDCtrans ( ) ? In addition, for the CR ( complex resistivity ) problem, the elements in the Jacobi matrix should also be complex numbers, which seems to be less intuitive than ERT. At this time, how should Coverage be calculated ? ( Maybe only use the amplitude of CR ? )
I am sorry that these two questions may not be specific, thank you for your patience !

[halbmy]

Q1: yes, zWeight is like az/ax (which is however only for regular grid), for details see Coscia et al. (2011)

Q2: The function is in the C++ core and sums up the absolute values of the jacobian and scales by the inner derivative. For the complex problem, the Jacobian consists of two real-valued matrices. The coverage should, however, be computed differently.

[ZSSSSZJU]

Thanks I think my problem has been solved. Thank you very much.

[halbmy]

For new topics we should create a new discussion with a proper title so that is better found.

[DustinRayMIZ]

Hey all, I was wondering if anyone knows of a way to save the simulated traveltimes for refraction tomography synthetic/forward models to a file. I can't seem to find anything about it in the documentation. I would like to try and compare the synthetic model I've gotten with a tomographic inversion using that same simulated data in Refrainv.py and compare them!