input sensors and shots locations and pick times for first arrival tomography

[rachelPython]

Problem description

I am not so experienced with seismic analysis. I use to use REFLEX-W and some Matlab codes to calculate tomography by first arrival Now I have a .PCK file which I managed to read to a matrix. the matrix is set of 7 columns as follows : 0: pick first arrival times 1,2,3 : x,y,z sensor location 4,5,6 : x,y,z source (shot) location I know how to insert a velocity model how do i adjust this to the input that is requested in gimly package? the output I want is a tomograpgy plot

Your environment

Operating system: Windows Python version: 3.6 pyGIMLi version: 1.0.9+14.g82028afc Way of installation: Conda package

Actual behavior

the guide as i understand tells me to get the shot and sensor both on the same input inside the mesh also i noticed there is a problem to insert repeating values for shot locations for example

91 #150 #question

[florian-wagner]

Hi @rachelPython,

the previous issue #165 should help. We use the unified data format, which is in analogy to the format used for ERT (https://gitlab.com/resistivity-net/bert#the-unified-data-format). First line contains the number of positions (can be shots, receivers, or both), the function is defined in the second section with column s for shot, g for geophone and t for your measured traveltime. The crosshole example explains how to construct such a file on-the-fly in Python.

[halbmy]

Actually, when using REFLEXW, you should be able to export your first-arrivals into a TOM (Reflex Tomography) file that looks like:

     69.15         0      7.20      7.50      0.00     12.75
     69.57         0      7.20      7.50      0.00     13.00
     70.19         0      7.20      7.50      0.00     13.25

with the first column being the traveltime, the 3rd/4th and 5th/6th are x/y or x/z coordinates of shot and geophone, at least I did it that way. This can be read by

from pygimli.physics.traveltime.importData import readTOMfile
data = readTOMfile(filename)

Actually, once there was a manager class Tomography (inherited from Refraction) that could be initialized by .TOM files directly (tomo = Tomography(filename)) but it seems outdated. We have to clean up those managers in pg v1.1.

[rachelPython]

Thank you ! In the crosshole example it is made on synthetic data and manually calculating and inserting all locations and traveltimes. this looks very long and comlex for me if i have this information from a file (.PCK or .TOM) as you mentioned: _data_ = readTOMfile(filename) i see there are two meshe's : one for the model and second one for the source-sensor geometry and time arrival data how do i insert the data to the model ?

[halbmy]

(sorry I mixed up two issues which is why I had to edit and delete comments)

After reading in the data you can create a Refraction manager and do inversion

from pygimli.physics import Refraction
data = readTOMfile(filename)
ra = Refraction(data)

But before you can do the inversion you should create a mesh (as the automatic mesh generation is made for refraction data) e.g. a regular mesh using pg.createMesh2D. Just have a look at the crosshole example.

[rachelPython]

Thank you. I have managed to move my data to the unified data format. The Refraction is working! :-) Now I am trying to do the inversion and i see that i am on the boundries of the CPU and might need to use GPU instead . is there a way to define this inversion function to work with GPU ?

[florian-wagner]

Hi @rachelPython,

GPU applications with pygimli do, to our knowledge, not exist. I would first try to reduce the parameters of your inverse problem (i.e., play with the mesh size). How many traveltimes are you inverting? How many cells does your model have? How much RAM do you have available?

[rachelPython]

I didnt specify the mesh and let it be the deafult value Data: Sensors: 65 data: 739 (traveltimes) Mesh: Nodes: 26945 Cells: 53658 Boundaries: 80602 I have 16GB RAM The run time was very long (~about 10 minutes), and ended with an error :

Traceback

(most recent call last): File "", line 97, in File "C:\Users\Owner\Anaconda3\envs\obspy\lib\site-packages\pygimli\physics\traveltime\refraction.py", line 442, in invert slowness = self.inv.run() RuntimeError: C:/msys64/home/Guenther.T/src/gimli/gimli36/src/inversion.h: 638 double GIMLI::Inversion::getPhiD(const Vec&) const [with ModelValType = double; GIMLI::Inversion::Vec = GIMLI::Vector] getPhiD == inf

[halbmy]

From a resolution point of view, 50000 mesh cells is far too much compared to the relatively low number of data. RAM requirement should be low, the computation times can be higher but should be managable.

At any rate, there is probably another problem, I guess due to very small distances and/or identical points. We can only help if you attach your data file and script.

[rachelPython]

Thank you very much for all your advice, this was very helpful . the problem was converting the traveltimes from [msec] to [sec] See my datafile attached in SGT format ‏‏sgtfile.txt this is a surface to borehole data, I believe i am getting something wrong with the geometry? this is my script - tomo.txt

[florian-wagner]

Hi @rachelPython,

a few things:

• ra.showVA and ra.showData currently assume all sources and receivers along the surface. Similarly, ert.showData shows classical electrical pseudosections for surface electrodes only. We cannot account for all acquisition geometries. So I am afraid you will need to find your own way to display the data, but that is just a visualization problem and the inversion will still work.
• Our convention is that the z-coordinate below the surface is negative.
• Although your approach with creating a regular grid would work (as in the crosshole example), you can think about using an unstructured mesh and cutting off the lower left and right parts of the model. This will save you inversion parameters and hence memory and computing time and will not affect the results, as no ray path is going to travel there.

sgtfile.sgt.txt tomo.py.txt

[rachelPython]

Thank you ! :) I would like to compare my inversion result to REFLEX-W result to see if it makes sense In order to do that I need to know the exact velocity model that comes as an input to the inversion function -

vest = ra.invert(lam=1100, vtop=300, vbottom=2000, zWeight=1.0)

I am trying to understand is it a gradient between vtop and vbottom ? also, is zWeight a smoothing parameter ?

[halbmy]

Sorry for the late response due to a conference that we all attended.

The starting model is, as usual in refraction tomography, a gradient model with increasing velocity. It is a simple interpolation of the depth below the surface to the interval [vTop, vBottom] so that these two velocities mark the extreme values for the shallowest and deepest cell. See also the source code in https://github.com/gimli-org/gimli/blob/7e1ba6f0ca757cf3ebf210afcb4f032bc7672de6/python/pygimli/physics/traveltime/ratools.py#L52

Note that the return is the slowness instead of velocity. I guess ReflexW will do something similar, but as the discretization is different, the initial velocity model will be slightly different. You can have a look at it by

ra = Refraction(data)
ra.invert(maxIter=0)  # do no inversion
ra.showResult()

for the BG_test.sgt data file of #165 as shown.

The resulting inversion model is of course very different from the starting model and can only be compared to the ReflexW result after some intepolation (see other issues and examples on that).

The classical smoothing parameter is called lam, but zWeight is probably more important: it defines the purely vertical roughness penalty in relation to a purely horizontal gradient and controls how "layered" you expect the model to be.

[florian-wagner]

@rachelPython Can we close this issue?

[GGDRriedel]

I would like to reopen this. I have a case where shots and geophones do not lign up. As far as I can see, the unified data format does not cover that, does it?

[carsten-forty2]

You can add any sensor position or sensor order in the unified data format as long as the indieces in the sensor token list 's' and 'g' match those sensors.

An ascending order is only needed if you want automatic mesh generation from the sensors in this data file. However, you can allways sort the DataContainer regarding the sensor positions.

[GGDRriedel]

Where is this documented?

I too am digging through the koeningssee example right now.

From that and the documentation I still have these questions left:

Is the flag 'y' actually the 'z' aka height from the ERT examples? Is this line even relevant to the input or just a comment? Clearly the lower

#s g t line is read by the program, is it? Do they have to be unique?

What other data formats are supported?

[carsten-forty2]

https://www.pygimli.org/gimliapi/classGIMLI_1_1DataContainer.html#a3d07be3931623a5ebef2bf14d06d4f50 Link to the api documentation of the DataContainer, they are usually not very verbosely documented and you need to play with it to see how they act. Feel free to support us by writing down your findings or improve these documentation.

'x', 'y', 'z' are carthesian coordinates where the last ('y' for 2D and 'z' for 3D) is the negative depth.

'# s g t' are data entries and don't need to be unique

This is just the description of the unified file format and it depends on the application if these informations are interpreted or not.

Regarding the file formats: Any new development is in the upcoming 1.1 version aka fw_cleaning branch. There is also generic load function for traveltime problems.

Btw. its better to open new issues for new questions.