Can add surface wave inversion into pygimli?

[xichaoqiang]

Problem description

Can surface wave dispersion curve inversion be added to Pygimli? Pygimli has a useful inversion framework. Surface wave inversion is also a common method in geophysical exploration. Pysurf96 is a program for the forward modeling of the surface wave dispersion (Rayleigh wave Love wave,group velocity or phase velocity).

[halbmy]

Great idea! I always had surface waves in mind to add to the methods for near-surface characterization. I guess you are talking about 1D for the moment?

As long as you can describe your forward algorithm in a way that a model vector yields a forward response vector (fop.response()) as described in the tutorials, inversion should be straightforward and can be followed by any other laterally constraint, joint or coupled inversion.

If you need help just contact us here, maybe first describing what you have by an example.

[xichaoqiang]

Yes,It is dispersion 1D inversion.

The attachment is pysurf96 package.The package contains two files, a Fortran77 code from Computer programs in seismology by R. Hermann and a python wrapper file.

The code shown in the below is some examples of dispersion curve forward.

I think it is easy to modify in Pygimli.

import numpy as np
from pysurf96 import surf96

# Define the velocity model in km and km/s
thickness = np.array([5., 23., 8., 0]) 
vs = np.array([2, 3.6, 3.8, 3.3])# shear wave velocity 
vp = vs * 1.73 #compressional wave velocity 
rho = vp * .32 + .77 #density (g/cm³)

# Periods we are interested in
periods = np.linspace(1., 20., 20) #periods=1./frequency  (s)

#fundamental mode rayleigh wave phase velocity dipsersion curve
velocities1 = surf96(thickness, vp, vs, rho, periods, wave='rayleigh', mode=1, velocity='phase', flat_earth=True)

#1st high-mode rayleigh wave phase velocity dipsersion curve
velocities2 = surf96(thickness, vp, vs, rho, periods, wave='rayleigh', mode=2, velocity='phase', flat_earth=True)

#fundamental mode love wave group velocity dipsersion curve
velocities3= surf96(thickness, vp, vs, rho, periods, wave='love', mode=1, velocity='group', flat_earth=True)

pysurf96-master.zip

[florian-wagner]

@xichaoqiang To avoid misunderstandings: We are happy to help you implementing surface wave inversion in pyGIMLi, but we are not gonna do it ourselves anytime soon. So I am closing this issue for now.

P.S.: You can use forward operators from a different package and invert with pygimli. You "just" have to write a wrapper and populate the response and optionally createJacobian methods. This is detailed in the pyGIMLi paper.

[xichaoqiang]

@xichaoqiang To avoid misunderstandings: We are happy to help you implementing surface wave inversion in pyGIMLi, but we are not gonna do it ourselves anytime soon. So I am closing this issue for now.

P.S.: You can use forward operators from a different package and invert with pygimli. You "just" have to write a wrapper and populate the response and optionally createJacobian methods. This is detailed in the pyGIMLi paper.

I think we need reopen this issue.I think we can achieve dispersion curve inversion with pyGIMLi.

I have wirte a fop.response(model) and fop.createJacboain(model) .

The response(model) function is OK.

But met a problem in inversion. A error message shows that.

""" coeff = inv.run() RuntimeError: C:/msys64/home/Guenther.T/src/gimli/gimli36/src/vector.h: 659 GIMLI::Vector& GIMLI::Vector::operator-=(const GIMLI::Vector&) [with ValueType = double] 71 != 99 """

The attachment is my code and file use in example. pygimli_dcinv.zip

EVODCINV package need install by "pip install evodcinv"

[florian-wagner]

Hi @xichaoqiang,

you are on the right track, but there are a number of small mistakes:

• line 45: createJacboain --> createJacobian
• line 46: freq --> self.freq
• line 47: rang --> range
• lines 52 + 53: forward --> self.forward

Please note that the finite-difference based calculation of the sensitivities that you do in this method is the pyGIMLi default, i.e. you can leave it out for the beginning. The error message tells you that there is a mismatch between the measurements and the model response.

Make sure that self.forward always gives back the same length of frequencies. I assume that this is not the case, but I don't know ECODCINV.

[xichaoqiang]

Hi @xichaoqiang,

you are on the right track, but there are a number of small mistakes:

• line 45: createJacboain --> createJacobian
• line 46: freq --> self.freq
• line 47: rang --> range
• lines 52 + 53: forward --> self.forward

Please note that the finite-difference based calculation of the sensitivities that you do in this method is the pyGIMLi default, i.e. you can leave it out for the beginning. The error message tells you that there is a mismatch between the measurements and the model response.

Make sure that self.forward always gives back the same length of frequencies. I assume that this is not the case, but I don't know ECODCINV.

ECODCINV is a dispersion curve package ,which inversion method is PSO algorithm.

The dispersion curve forward method of my code is base on it.

I correct my mistake and,make sure length of frequencies equals to velocity.

[florian-wagner]

ECODCINV is a dispersion curve package ,which inversion method is PSO algorithm.

The dispersion curve forward method of my code is base on it.

Yes, I understand that, but I mean I don't know in detail what the function does, which you are using within self.forward. But this has to be clear to you, not to me. Just make sure that it always produces the same length response vector.

(P.S.: The PSO inversion algorithm in ECODCINV is irrelevant as you are only using its forward capabilities here.)

[xichaoqiang]

ECODCINV是一种色散曲线包，其反演方法是PSO算法。 我的代码的色散曲线前向方法基于它。

是的，我理解这一点，但我的意思是我不详细地知道你正在使用的功能是什么self.forward。但这一点必须清楚，而不是我。只需确保它始终生成相同长度的响应向量。

（PS：ECODCINV中的PSO反演算法无关紧要，因为此处仅使用其前向功能。）

Yes,the result of my start model is not right.Thanks. I will correct it.

[xichaoqiang]

It works. It is briliant.

How can it get the inversion information of each iteration?

red point is observe data.

Blue line in Fig.1 is dispersion curve of initial model

Blue line in Fig.2 is dispersion curve of inverted model

[florian-wagner]

Hey @xichaoqiang,

glad to hear that it works!

There is some output in the console (which you may not see in Spyder or Jupyter Notebooks). With pygimli 1.1, we will have better control to access the information from each iteration.

[xichaoqiang]

嘿@xichaoqiang，

很高兴听到它有效！

控制台中有一些输出（您可能在Spyder或Jupyter笔记本中看不到）。使用pygimli 1.1，我们将有更好的控制来访问每次迭代的信息。

OK,thanks.I'll try the new version.

[carsten-forty2]

If you miss console output you can init the inversion with verbose=True or debug=True

1.1. is not yet ready, however there is an unstable branch for called fw_cleaning. You can try but there are some slight api and concept changes, so its likely that you need to adapt your scripts.

[xichaoqiang]

If you miss console output you can init the inversion with verbose=True or debug=True

1.1. is not yet ready, however there is an unstable branch for called fw_cleaning. You can try but there are some slight api and concept changes, so its likely that you need to adapt your scripts.

Thank you.I only get the information of final iteration.

Can't I get the inforamtion of each iteration now like tomography?

""" 4: Model: min = 195.269; max = 738.592 4: Response: min = 201.333; max = 734.856 4: rms/rrms(data, Response) = 1.97832/0.94894% 4: chi^2(data, Response, error, log) = 391.375 4: Phi = 7827.5+217.002*51.2=7827.5 1.9783195244976863 0.9489398145678535 391.3748141008751 6 [195.2693454832814, 268.9549553771356, 342.72085151520116, 517.7852642621755, 610.0854987811659, 738.5919511674316] """

[carsten-forty2]

One reason for the development of 1.1. is to move the main inversion loop from the c++ core into pure python so we have more control and possibilities. Until we finished 1.1., its a little extensive to collect informations per iteration.

You can manually loop with inv.oneStep() or if you set the inv = RInversion(..., dosave=True), you can get a vector of models for each iteration with inv.modelHistory(). With this models you can calculate responses and fitting values manually.

[xichaoqiang]

Thank you.

[xichaoqiang]

For the stability of the inversion, can I set the step value of each iteration not to exceed a certain percentage of the model?

I know we can set upper and bottom value of inversion model.

[halbmy]

I guess you are referring to the Line Search, i.e. determining optimum step length of the model update. This is done internally using a quadratic function of step lengths 0 (old model), 1 (full update) and 0.3 with a fallback to a small (0.1) step length. So far you cannot change or deactivate that, but I am confident it is going to work. In pg v1.1 we will enable choosing the linesearch method more specifically.

[xichaoqiang]

I guess you are referring to the Line Search, i.e. determining optimum step length of the model update. This is done internally using a quadratic function of step lengths 0 (old model), 1 (full update) and 0.3 with a fallback to a small (0.1) step length. So far you cannot change or deactivate that, but I am confident it is going to work. In pg v1.1 we will enable choosing the linesearch method more specifically.

Hello, I have three question about surface wave inversion with version 1.1.

1.How can I get the data information of each iteration with the newer version,so I can determine whether it is overfitting?

2.Can I set the lower and the upper bounds in the inversion.

3.Can I set the reference model of the object function?

[halbmy]

As there are several ways to Rome: Would you like to share your code (or at least the pygimli-related parts)? Then we could help you.

Are you still using pysurf96? Or anything different or an own open-source code? Would be interesting to have this as an example for pyGIMLi...

[xichaoqiang]

Thanks for your reply.I would like to add this as an example for pyGIMLi.

Recently I am using disba as a tool for forward modeling of my dispersion curve. This is an open source package of python. https://github.com/keurfonluu/disba.

The attachment is my inversion program but not test for my conda environment has been broken. :(

pygimili_dcinv.zip

[halbmy]

Thanks for the code. The dispa model takes the following layer properties:

• thickness
• p-wave velocity $v_p$
• s-wave velocity $v_s$
• density But you seem to invert only for $v_s$. I know it is hard to find all properties without additional knowledge, but in reality you usually don't know the other properties exactly.

[xichaoqiang]

Thanks for the code. The dispa model takes the following layer properties:

• thickness
• p-wave velocity $v_p$
• s-wave velocity $v_s$
• density But you seem to invert only for $v_s$. I know it is hard to find all properties without additional knowledge, but in reality you usually don't know the other properties exactly.

Yes. For near surface applications , surface wave velocity is the most sensitive to s-wave velocity than other parameters. Generally speaking, only inversion of the S-wave velocity is sufficient.

1. We can set the inversion model as a collection of many thin layers. (Thickenss)

2. P-wave velocity is not very sensitive to surface wave velocity. We can get the P-wave velocity by fixing the poisson’s ratio of layers in the inversion. (Vp)

3. Density is the least sensitive to surface wave velocity. So it is usually set to 2000 or other values based on prior information. (Density)

[xichaoqiang]

pygimli_dic.zip

The forward computation is OK.

Howerver, I met a error message in the inversion step.

Can you help me check my code? The code is in the attachment.

Here is the error message.

8/01/21 - 17:17:48 - pyGIMLi - INFO - Starting inversion.
min/max(dweight) = 2/2
18/01/21 - 17:17:48 - pyGIMLi - INFO - Set default startmodel to median(data values)=202.27420563453796
18/01/21 - 17:17:48 - pyGIMLi - INFO - Created startmodel from forward operator: 0
fop: <__main__.DCModelling object at 0x7fcf1c04b570>
Data transformation: <pygimli.core._pygimli_.RTransLin object at 0x7fcf04392e30>
Model transformation: <pygimli.core._pygimli_.RTransLog object at 0x7fcf04386d30>
min/max (data): 185/693
min/max (error): 0.072%/0.27%
Traceback (most recent call last):

  File "/home/xcq/下载/pygimli_dcinv_2021.py", line 58, in <module>
    coeff = inv.run(vr,relError,verbose=True)

  File "/home/xcq/.conda/envs/pyinv/lib/python3.7/site-packages/pygimli/frameworks/inversion.py", line 411, in run
    print("min/max (start model): {0}/{1}".format(pf(min(self.startModel)),

ValueError: min() arg is an empty sequence

[halbmy]

I've had a look at your code, it is really very good.

The main error came from the fact that the starting model did not match the parameter count as you did not set a mesh. I included the lines

self.mesh_ = pg.meshtools.createMesh1D(len(self.thk))
self.setMesh(self.mesh_)

to create a 1D mesh (of which only vp is inverted but not its thickness), e.g. an Occam type inversion. If you do so, the starting model is accepted and the inversion starts running. However, I get an error in disba

DispersionError: failed to find root for fundamental mode

which is beyond my expertise. I attached my script. pygimli_dcinv_2021TG.zip

Other comments:

• setting an empty Jacobian matrix does not make sense unless you create a function createJacobian where it is filled, if you comment this line it is computed by the perturbation mode
• Marquardt inversion makes more sense for a block inversion (including thickness), I commented this in favor of smoothness constraints (1st order, default for any 1D/2D/3D mesh)
• you should use the same error for noisifying the data and the error weighting
• I did a few more cosmetic operations like renaming forward to response

[xichaoqiang]

pygimli_dcinv_0120.zip

It can be run now.I can get a good inversion result.

Thank you very much.

1.How can I save the response result of each iteration,so I can determine whether it is overfitting?

2.Can I set the lower and the upper bounds in the inversion.

3.Can I set the reference model of the object function?

[halbmy]

An error of 0.1 is very unrealistic, I guess 10 is a more realistic magnitude but it works well with it:

1. There is a model history available through inv.modelHistory. So far there is no response history (but probably helpful), but you can foward model the individual responses. However you can also do other things in the inversion by setting inv.setPostStep(callMePerIteration) to a function def callMePerIteration:. However, taking an early model (by looking at inv.ch2History is in my opinion not a good choice, you should rather vary inv.run(lam= so that the final chi-square is one (in your case 50 is a good number yielding chi2=1.06).

[halbmy]

2. A lower bound can be set by f.modelTrans=pg.trans.TransLog(150), two-sided bounds by f.modelTrans=pg.trans.TransLogLU(150, 900)

[halbmy]

3. By default, the starting model is also a reference model (but we should allow disabling this by a switch). In case of smoothness constraints, this does not make a difference for a homogeneous model but in your case you are specifying the synthetic model plus some shift which does not make sense (therefore the results are so close to the synthetic). When trying to use a homogeneous starting model by inv.run(..., startModel=np.ones_like(vs)*300, I am running into the surf96 error. I suggest using a gradient model with a low gradient such as np.arange(len(vs))*10+200, then all runs well.

[halbmy]

Thank you again for this very interesting issue. We should include this as an example.

• We will have to use a different depth discretization (more layers) to avoid the inverse crime.
• How realistic is it to invert field measurements? Right now you assume knowing layer boundaries and Vp (e.g. from VSP).
• Will you have to adapt it so that the poisson ratio is fixed instead of the P-wave velocity?
• How important is the choice of the density?

[xichaoqiang]

Tanks for your reply.

1. Can I set the boundaries with numpy.array ,the length of boundary is equal to number of layer
2. I will adapt the code to fix the poisson ratio
3. We always neglected the density in surface wave inversion

[halbmy]

Every model region has its own model transformation. After calling mesh=mt.createMesh1D(10) you could set the markers (mesh.cell(i).setMarker(j)) or alternatively use createMesh2D(10, 1, markerType=1), i.e. a 2D mesh with just one row and increasing markers from 1 to 10.

[xichaoqiang]

pygimli_dcinv_20210120_possion.zip

In this code,we can fix the possion in the inversion step.

I think response history will be a good optition for me.

[halbmy]

I like the Poisson-based code and would be interested whether it can be applied to real data?

As written before the inversion should be either a Block-1D model with variable thickness or an Occam inversion with a higher number of predefined layers. Right now we assume knowing the geometry.

We will add the response history as this is a good idea.

[xichaoqiang]

I think it can be applied to real data.

I am also intrested in inverting dispersion curve like occam inversion or block-1D model.

I think it can be applied like other data in a similar way.