Making 3D velocity models

[Ali1990dashti]

Dear community,

First and foremost, thanks for making open-source codes and for your support.
I recently started my work with pyGIMLi and want to make a simple 3D Refraction modeling and inversion with it. The 2D example is presented here. I am using pyGIMLi 1.4.0.
This is my code:

import numpy as np
import matplotlib.pyplot as plt
import pygimli as pg
from pygimli import meshtools as mt
from pygimli.physics import TravelTimeManager
from pygimli.viewer import showMesh

# dimensions of the geometry
left = -30
right = 30
depth = 50

numberGeophones = 30
zeros_column = np.zeros((numberGeophones, 1))
pos = np.linspace(-right, right, numberGeophones)
sens = np.hstack((pos.reshape((-1, 1)), zeros_column)) # it makes nodes at the sensor locations

# meshing block
geometry = mt.createWorld(start=[left, 0],
                       end=[right, -depth],
                       layers=[-5, -20]) # I considered three layers in the model
for sen in sens:
    geometry.createNode(sen)
mesh = mt.createMesh(geometry, 
                     area=1.0,
                     quality=33,
                     smooth=[2, 4]
                    )

extruded_mesh = mt.extrudeMesh(mesh, np.linspace(0, 10, 5)) # adding dimension Z in this case 
extruded_mesh.rotate([np.pi/2, 0, 0]) # rotating mesh to switch y/z direction and view the top of the extended mesh. 

sensors = np.linspace(-right, right, numberGeophones)
scheme = pg.physics.traveltime.createRAData(sensors)

velocities = [200, 400, 800] # material properties
mgr = TravelTimeManager()
vp = np.array(extruded_mesh.cellMarkers())
vp[vp == 1] = velocities[0]
vp[vp == 2] = velocities[1]
vp[vp == 3] = velocities[2]

data = mgr.simulate(slowness=1.0 / vp, scheme=scheme, mesh=extruded_mesh,
                    noiseLevel=0.001, noiseAbs=0.001, seed=1337,
                    verbose=True) # this block of codes so long to be run!
mgr.showData(data)

vest = mgr.invert(data, secNodes=2, paraMaxCellSize=15.0,
                  maxIter=100, verbose=True)
np.testing.assert_array_less(mgr.inv.inv.chi2(), 1.1)

And the following graph comes with this line of code:

ax, _ = mgr.showResult(cMin=min(vp), cMax=max(vp), logScale=False)
pg.show(geometry, ax=ax, fillRegion=False, regionMarker=False)

But as you can see the data only covers depths to -22 m:

To be honest I am not sure which mistakes I have done. Therefore, I am looking for hints from people who are using pyGIMLi.

Thanks

[halbmy]

I don't fully get it. You have a 2D profile of sensors, but you do a 3D synthetic model computation (why?). Finally you invert the data 2D as this is the default if you don't specify a mesh. By default, the mesh depth is derived from the shot-geophone distance and 22m is probably deep enough. I don't think that rays with such limited offsets and small contrasts will reach the third layer, but you can check when looking at the first-arrival traveltime curves slope. But anyway you can set the model depth by paraDepth.

[Ali1990dashti]

Thanks for dedicating time to my issue.
I thought maybe by the extrusion the sensor position will also changed to 3D. When I call print(scheme.sensorPositions().array()) it gives me an array of x, y and z values. This is why i thought sensor positions are 3D.
I very much appreciate if you let me know how I can implement sensors correctly in such model.

[halbmy]

Positions are always in the 3D space, this holds also for sensors. So of course, you can create a spread of 3D positions, but there is no ready-made createData function for it. Then you also would have to create a 3D inversion mesh and pass it. But I still don't understand where the problem is.

[Ali1990dashti]

Sorry for my naive questions and thanks for your help.
Is there any example available showing how to import a 3D mesh ad and use it for seismic inversion? My mesh has been made in GMSH and I defined the physical properties of the mesh, i.e. each volume has its specific name. then, I need to set the velocity of each volume and in the last step simulation and inversion. My mesh is here. My problem is that I do not know how to do it and very much appreciate it if you provide me with similar examples.