Note
Go to the end to download the full example code.
Field data inversion (“Koenigsee”)#
This minimalistic example shows how to use the Refraction Manager to invert a field data set. Here, we consider the Koenigsee data set, which represents classical refraction seismics data set with slightly heterogeneous overburden and some high-velocity bedrock. The data file can be found in the pyGIMLi example data repository.
# We import pyGIMLi and the traveltime module.
import matplotlib.pyplot as plt
import pygimli as pg
import pygimli.physics.traveltime as tt
The helper function pg.getExampleData downloads the data set to a temporary location and loads it. Printing the data reveals that there are 714 data points using 63 sensors (shots and geophones) with the data columns s (shot), g (geophone), and t (traveltime). By default, there is also a validity flag.
data = pg.getExampleData("traveltime/koenigsee.sgt", verbose=True)
print(data)
Data: Sensors: 63 data: 714, nonzero entries: ['g', 's', 't', 'valid']
Let’s have a look at the data in the form of traveltime curves.
fig, ax = plt.subplots()
lines = tt.drawFirstPicks(ax, data)

We initialize the refraction manager.

Finally, we call the invert method and plot the result.The mesh is created based on the sensor positions on-the-fly.
mgr.invert(secNodes=3, paraMaxCellSize=5.0,
zWeight=0.2, vTop=500, vBottom=5000, verbose=1)
fop: <pygimli.physics.traveltime.modelling.TravelTimeDijkstraModelling object at 0x1717d9fd0>
Data transformation: Identity transform
Model transformation (cumulative):
0 Logarithmic LU transform, lower bound 0.0, upper bound 0.0
min/max (data): 3.5e-04/0.03
min/max (error): 3%/3%
min/max (start model): 2.0e-04/0.002
--------------------------------------------------------------------------------
inv.iter 0 ... chi² = 156.33
--------------------------------------------------------------------------------
inv.iter 1 ... chi² = 12.32 (dPhi = 91.49%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 2 ... chi² = 8.91 (dPhi = 27.39%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 3 ... chi² = 6.41 (dPhi = 27.46%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 4 ... chi² = 5.53 (dPhi = 13.07%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 5 ... chi² = 4.95 (dPhi = 10.00%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 6 ... chi² = 4.36 (dPhi = 10.88%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 7 ... chi² = 3.95 (dPhi = 8.51%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 8 ... chi² = 3.71 (dPhi = 5.44%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 9 ... chi² = 3.61 (dPhi = 2.35%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 10 ... chi² = 3.51 (dPhi = 2.41%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 11 ... chi² = 3.46 (dPhi = 1.19%) lam: 20.0
################################################################################
# Abort criterion reached: dPhi = 1.19 (< 2.0%) #
################################################################################
1090 [874.4940674147222,...,2692.581750977291]
Look at the fit between measured (crosses) and modelled (lines) traveltimes.
mgr.showFit(firstPicks=True)

You can plot only the model and customize with a bunch of keywords

You can play around with the gradient starting model (vTop and vBottom arguments) and the regularization strength lam and customize the mesh.
Total running time of the script: (0 minutes 10.716 seconds)