Layered initial model for refraction seismic inversion

[GiuliadePasqualeCEAZA]

Problem description

Hello,
I was surfing a bit of the API reference on the pygimli website cause I am interested into set a traveltime iversion with an initial layered model. Do you know if it is possible to set an other modelrather than the gradient one as starting point of the inversion? How could I set up this problem?

thank you in advance
Cheers
Giulia

Your environment

Please provide the output of print(pygimli.Report()) here. If that does not
work, please give provide some additional information on your:

Operating system: e.g. Windows, Linux or Mac?
Python version: e.g. 3.9, 3.10, etc.?
pyGIMLi version: Output of print(pygimli.__version__)
Way of installation: e.g. Conda package, manual compilation from source, etc.

Steps to reproduce

Tell us how to reproduce this issue. Ideally, you could paste the code that produces the error:

import pygimli as pg
...

Expected behavior

Tell us what should happen or what you want to achieve.

Actual behavior

Tell us what happens instead and/or provide the output of your script.

Paste your script output here.

If possible, please add one or more labels to your issue, e.g. if you expect that your issue is rather a question than a problem with the code, please add the label "question".

[halbmy]

Of course you can, you just have to provide a starting model vector and for this you will need the mesh. For a strictly layered mesh you will have to include these as boundaries. If you have good knowledge on the layer boundaries and velocities you should consider creating a mesh with different regions and then inv.setRegularization(nr, start=slowness_nr).

[halbmy]

Would be interesting to create an example for it. There might also be an issue with automatic background regions, so better use inv.setRegularization(background=False) before.

[GiuliadePasqualeCEAZA]

Thank you for the rapid answer: if you do create an example let me know to check if I'm doing it right
cheers
Giulia

[halbmy]

Maybe you can start doing it, we can help you and finally make an example together? That would be great!

[GiuliadePasqualeCEAZA]

Agrred: I was just trying something fairly easy adapting the 2D Refraction modelling and inversion example you have on the website:

###############################################################################################

import numpy as np
import pygimli as pg
import pygimli.meshtools as mt
import pygimli.physics.traveltime as tt


#1) Create synthetic dataset based on a two-layers model
geom = mt.createWorld(start=[-10., 80], end=[115, 135], layers=[120])
pg.show(geom)

mesh = mt.createMesh(geom, quality=34.3, area=3, smooth=[1, 10])
ax, _ = pg.show(mesh)

numberGeophones = 48
sensors = np.linspace(0., 117.5, numberGeophones)
scheme = tt.createRAData(sensors)

vp = np.array(mesh.cellMarkers())
vp[vp == 1] = 250
vp[vp == 2] = 1300

ax, _ = pg.show(mesh, vp, colorBar=True, logScale=False, label='v in m/s')
pg.viewer.mpl.drawSensors(ax, scheme.sensors(), diam=1.0, facecolor='white', edgecolor='black')

data = tt.simulate(slowness=1.0 / vp, scheme=scheme, mesh=mesh,
                   noiseLevel=0.001, noiseAbs=0.001, seed=1337, verbose=True)
tt.show(data)


#2) Setting inversion mesh and inversion parameters data
mgr      = tt.TravelTimeManager(data)
geom_Inv = mt.createWorld(start=[0.0, 110], end=[115, 135], layers=[120])
pg.show(geom)
mesh_Inv = mt.createMesh(geom_Inv, quality=30.2, area=5, smooth[1, 10])
pg.show(mesh_Inv)
mgr.applyMesh(mesh_Inv)
mgr.inv.setRegularization(background=False)
#here finally I set the initial model through this function cause the inv.setRegularization(nr, start=slowness_nr) gave me an error
mgr.fop.setRegionProperties(regionNr=1, startModel=100)      
mgr.fop.setRegionProperties(regionNr=2, startModel=1000)
vest = mgr.invert(data=data, mesh=mesh_Inv, secNodes=2,maxIter=10, verbose=True)

###############################################################################

the error mesage I get at this point is:

 CRITICAL - <class 'pygimli.physics.traveltime.TravelTimeManager.TravelTimeManager'>._ensureError(C:\Users\Giulia\anaconda3\envs\pg\lib\site-packages\pygimli\frameworks\methodManager.py:382)
All error values need to be larger then 0. Either give and err argument or fill dataContainer  with a valid 'err'  -82.54067435292752 18.79883003318321
Traceback (most recent call last):
  File "C:\Users\Giulia\Desktop\OneDrive\CEAZA\CODES\Exemplos_Pygimli\2D_Refraction_seismic_layered-Initial-Model.py", line 40, in <module>
    vest = mgr.invert(data=data, mesh=mesh_Inv, secNodes=2,maxIter=10, verbose=True)
  File "C:\Users\Giulia\anaconda3\envs\pg\lib\site-packages\pygimli\physics\traveltime\TravelTimeManager.py", line 254, in invert
    slowness = super().invert(data, mesh, **kwargs)
  File "C:\Users\Giulia\anaconda3\envs\pg\lib\site-packages\pygimli\frameworks\methodManager.py", line 818, in invert
    errorVals = self._ensureError(self.fop.data, dataVals)
  File "C:\Users\Giulia\anaconda3\envs\pg\lib\site-packages\pygimli\frameworks\methodManager.py", line 382, in _ensureError
    pg.critical("All error values need to be larger then 0. Either"
  File "C:\Users\Giulia\anaconda3\envs\pg\lib\site-packages\pygimli\core\logger.py", line 301, in critical
    raise Exception(_msg(*args))
Exception: All error values need to be larger then 0. Either give and err argument or fill dataContainer  with a valid 'err'  -82.54067435292752 18.79883003318321

##############

which somehow is related to the data error, but when I look into the error of the data are actually fine (between 0.001 and 0.0010886418269230769). So not sure what I am doing wrong

thank you in advance
Cheers
Giulia

[halbmy]

Your velocity distribution looks like that

but your sensors are still at the surface (y=z=0), so they are outside of the modelling domain. I guess they are supposed to be at the surface? Then set their y position to 130, or consider moving your model so that z=0 is the surface.

[GiuliadePasqualeCEAZA]

Great: thank you for the hint!

I modify my "world" and now it start the inversion but I think the problem is that is still assigning the region with smaller marker to background (even thought I run the command " mgr.inv.setRegularization(background=False)"). So the inversion mesh results doesn't appears as the one I set in the codes and the results did not reflect the inicial model.
I write here the code I used: any idea/help/suggestions?

thank you in advance for your time
cheers
Giulia

import numpy as np
import pygimli as pg
import pygimli.meshtools as mt
import pygimli.physics.traveltime as tt


#1) Create synthetic dataset based on a two-layers model
geom = mt.createWorld(start=[-10., -60.], end=[120., 0.], layers=[-20.])
pg.show(geom)

mesh = mt.createMesh(geom, quality=34.3, area=3, smooth=[1, 10])
ax, _ = pg.show(mesh)

numberGeophones = 48
sensors = np.linspace(0., 117.5, numberGeophones)
scheme = tt.createRAData(sensors)

vp = np.array(mesh.cellMarkers())
vp[vp == 1] = 250
vp[vp == 2] = 1300

ax, _ = pg.show(mesh, vp, colorBar=True, logScale=False, label='v in m/s')
pg.viewer.mpl.drawSensors(ax, scheme.sensors(), diam=1.0, facecolor='white', edgecolor='black')

data = tt.simulate(slowness=1.0 / vp, scheme=scheme, mesh=mesh,
                   noiseLevel=0.001, noiseAbs=0.001, seed=1337, verbose=True)
tt.show(data)


#2) Setting inversion mesh and inversion parametersdata
mgr      = tt.TravelTimeManager(data)
geom_Inv = mt.createWorld(start=[0.0, -40.], end=[115, 0.], layers=[-20.])
pg.show(geom_Inv)
mesh_Inv = mt.createMesh(geom_Inv, quality=30.2, area=5, smooth=[1, 10])
pg.show(mesh_Inv)
mgr.inv.setRegularization(background=False)
mgr.applyMesh(mesh_Inv)
mgr.fop.setRegionProperties(regionNr=1, startModel=1/100.)
mgr.fop.setRegionProperties(regionNr=2, startModel=1/1000.)
vest = mgr.invert(data=data, useGradient=False, secNodes=2, maxIter=10, verbose=True)
np.testing.assert_array_less(mgr.inv.inv.chi2(), 1.1)
ax, _ = mgr.showResult(cMin=min(vp), cMax=max(vp), logScale=False)
ax, _ = pg.show(geom_Inv, ax=ax, fillRegion=False, regionMarker=False)

ax, _ = pg.show(mgr.mesh, vest, label="Velocity [m/s]")

[halbmy]

I did a few things:

1. I first created the data container and added the sensors to the world thus ensuring a better forward quality (and by the way a much nicer mesh)
2. I used setMesh instead of applyMesh and set the regularization options AFTER setting the mesh.
   This lead to a quite nice result with a chi2 below 3 which can definitely be further improved
   

numberGeophones = 48
sensors = np.linspace(0., 117.5, numberGeophones)
scheme = tt.createRAData(sensors)

#1) Create synthetic dataset based on a two-layers model
geom = mt.createWorld(start=[-10., -60.], end=[120., 0.], layers=[-20.])
# pg.show(geom)
for sen in scheme.sensors():
    geom.createNode(sen)

mesh = mt.createMesh(geom, quality=34.3, area=4, smooth=[1, 10])
ax, _ = pg.show(mesh, markers=True, showMesh=True)

vp = np.array(mesh.cellMarkers())
vp[vp == 1] = 250
vp[vp == 2] = 1300

ax, _ = pg.show(mesh, vp, showMesh=True, colorBar=True, logScale=False, label='v in m/s')
pg.viewer.mpl.drawSensors(ax, scheme.sensors(), diam=1.0, facecolor='white', edgecolor='black')

data = tt.simulate(slowness=1.0 / vp, scheme=scheme, mesh=mesh,
                   noiseLevel=0.001, noiseAbs=0.001, seed=1337, verbose=True)
tt.show(data)

#2) Setting inversion mesh and inversion parametersdata
mgr      = tt.TravelTimeManager(data)
geom_Inv = mt.createWorld(start=[0.0, -40.], end=[115, 0.], layers=[-20.])
# pg.show(geom_Inv)
mesh_Inv = mt.createMesh(geom_Inv, quality=30.2, area=5, smooth=[1, 10])
# pg.show(mesh_Inv)
mgr.setMesh(mesh_Inv)
mgr.inv.setRegularization(background=False)
mgr.fop.setRegionProperties(regionNr=1, startModel=1/100.)
mgr.fop.setRegionProperties(regionNr=2, startModel=1/1000.)
vest = mgr.invert(data=data, useGradient=False, secNodes=2, maxIter=10, verbose=True)
ax, _ = mgr.showResult(cMin=min(vp), cMax=max(vp), logScale=False)

[halbmy]

If you make the model big enough (the last sensor was outside of the inversion region) and add the sensors you can fit the data with chi^2=1 with a practically homogeneous model

geom_Inv = mt.createWorld(start=[0.0, -40.], end=[120, 0.], layers=[-20.])
for sen in data.sensors():
    geom_Inv.createNodeWithCheck(sen)
mesh_Inv = mt.createMesh(geom_Inv, quality=34, area=5, smooth=[1, 10])
mgr.setMesh(mesh_Inv, secNodes=3)

[GiuliadePasqualeCEAZA]

Thank you so much!
Just run the code on my laptop but the results look quite different...

Could it be something with my version of pygimli (1.4.1)?
cheers
Giulia

[GiuliadePasqualeCEAZA]

Here the code I run:

import numpy as np
import pygimli as pg
import pygimli.meshtools as mt
import pygimli.physics.traveltime as tt

numberGeophones = 48
sensors = np.linspace(0., 117.5, numberGeophones)
scheme = tt.createRAData(sensors)

#1) Create synthetic dataset based on a two-layers model
geom = mt.createWorld(start=[-10., -60.], end=[120., 0.], layers=[-20.])
# pg.show(geom)
for sen in scheme.sensors():
    geom.createNode(sen)

mesh = mt.createMesh(geom, quality=34.3, area=4, smooth=[1, 10])
ax, _ = pg.show(mesh, markers=True, showMesh=True)

vp = np.array(mesh.cellMarkers())
vp[vp == 1] = 250
vp[vp == 2] = 1300

ax, _ = pg.show(mesh, vp, showMesh=True, colorBar=True, logScale=False, label='v in m/s')
pg.viewer.mpl.drawSensors(ax, scheme.sensors(), diam=1.0, facecolor='white', edgecolor='black')

data = tt.simulate(slowness=1.0 / vp, scheme=scheme, mesh=mesh,
                   noiseLevel=0.001, noiseAbs=0.001, seed=1337, verbose=True)
tt.show(data)

#2) Setting inversion mesh and inversion parametersdata
mgr      = tt.TravelTimeManager(data)
geom_Inv = mt.createWorld(start=[0.0, -40.], end=[120, 0.], layers=[-20.])
for sen in data.sensors():
    geom_Inv.createNodeWithCheck(sen)
mesh_Inv = mt.createMesh(geom_Inv, quality=34, area=5, smooth=[1, 10])
mgr.setMesh(mesh_Inv, secNodes=3)
mgr.inv.setRegularization(background=False)
mgr.fop.setRegionProperties(regionNr=1, startModel=1/100.)
mgr.fop.setRegionProperties(regionNr=2, startModel=1/1000.)
vest = mgr.invert(data=data, useGradient=False, secNodes=2, maxIter=10, verbose=True)
ax, _ = mgr.showResult(cMin=min(vp), cMax=max(vp), logScale=False)

[makeabhishek]

I tried the code provided by @halbmy. It gave me similar results but not exactly same.

Sorry for jumping here, but I also have a query related to this question. How can we fix the slowness of a layer so it does not change during inversion. the function to do that doesn't give correct result.
mgr.inv.setRegularization(1,fix=1/2300)

[andieie]

Hello ! Sorry to also cut in. I am running a similar example but instead using
mgr.inv.setRegularization(regionNr=1, single=True) to get back a single value for every layer. My problem is that when I pass a startModel it gives a bad and singular result for all of the regions. This does not happen for ERT but I am running into this problem in TT. Maybe it is similar to what @makeabhishek is experiencing. If I don't set a startModel, the results look good.

[halbmy]

Do you set a starting value in the inv.setRegularization() or in the inv.run() call?

[andieie]

I set it in the inv.setRegularization() as Carsten suggested on Mattermost. It looks more like this:

mgr.fop.setRegionProperties(regionNr=1, single=True, startModel=1/1000)

It works well with ERT but I am having trouble getting sound results with physics.traveltime.TravelTimeManager.

[makeabhishek]

Hi @andieie
If you can share your inversion code snippet or better a full code, that would be easier to follow.

[andieie]

I set up the same example with the synthetic data from @GiuliadePasqualeCEAZA above. and then the type of inversion I am doing is this :

#1) Create synthetic dataset based on a two-layers model
geom = mt.createWorld(start=[-10., -60.], end=[120., 0.], layers=[-20.])
pg.show(geom)

mesh = mt.createMesh(geom, quality=34.3, area=3, smooth=[1, 10])
ax, _ = pg.show(mesh)

numberGeophones = 48
sensors = np.linspace(0., 117.5, numberGeophones)
scheme = tt.createRAData(sensors)

vp = np.array(mesh.cellMarkers())
vp[vp == 1] = 250
vp[vp == 2] = 1300
ax, _ = pg.show(mesh, vp, colorBar=True, logScale=False, label='v in m/s')
pg.viewer.mpl.drawSensors(ax, scheme.sensors(), diam=1.0, facecolor='white', edgecolor='black')

data = tt.simulate(slowness=1.0 / vp, scheme=scheme, mesh=mesh,
                  noiseLevel=0.001, noiseAbs=0.001, seed=1337, verbose=True)
tt.show(data)

#2) Setting inversion mesh and inversion parametersdata
mgr  = tt.TravelTimeManager(data)
plc= mt.createWorld(start=[0.0, -40.], end=[115, 0.], layers=[-20.])
mesh_Inv = mt.createMesh(plc, quality=30.2, area=5, smooth=[1, 10])
geom_Inv = mt.appendTriangleBoundary(mesh_Inv, marker=0, xbound=10, ybound=10, isSubSurface=False)
mgr.setMesh(geom_Inv)
mgr.inv.setRegularization(background=False)
mgr.inv.setRegularization(1, single=True, startModel=1/1000.)
mgr.inv.setRegularization(2, single=True, startModel=1/2000.)
vest = mgr.invert(data=data, useGradient=False, secNodes=2
ax, _ = mgr.showResult(cMin=min(vp), cMax=max(vp), logScale=False)

I guess I don't need to create a larger mesh/append boundary, but I was just curious as to why this doesn't give the correct result.

[makeabhishek]

May be if you can just use

mgr.inv.setRegularization(background=False)
mgr.inv.setRegularization(1, startModel=1/1000.)
mgr.inv.setRegularization(2, startModel=1/2000.)

Or
mgr.inv.setRegularization(1, limits=[1/1000, 1/1100])