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 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)
[::::::::::::::::::::::::::::::::::::: 83% :::::::::::::::::::::::             ] 8193 of 9844 complete
[:::::::::::::::::::::::::::::::::::: 100% ::::::::::::::::::::::::::::::::::::] 9844 of 9844 complete
md5: 641890bb17cb2bdf052cbc348669dfd0
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 = pg.plt.subplots()
tt.drawFirstPicks(ax, data)
plot 04 koenigsee

We initialize the refraction manager.

mgr = tt.TravelTimeManager(data)

# Alternatively, one can plot a matrix plot of apparent velocities which is the
# more general function also making sense for crosshole data.

ax, cbar = mgr.showData()
plot 04 koenigsee

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 0x7f3369b32ea0>
Data transformation: <pygimli.core._pygimli_.RTrans object at 0x7f3379533130>
Model transformation (cumulative):
         0 <pygimli.core._pygimli_.RTransLogLU object at 0x7f33793b7e80>
min/max (data): 3.5e-04/0.03
min/max (error): 3%/3%
min/max (start model): 2.0e-04/0.002
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
inv.iter 2 ... chi² = 9.5 (dPhi = 29.31%) lam: 20
--------------------------------------------------------------------------------
inv.iter 3 ... chi² = 8.28 (dPhi = 12.52%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 4 ... chi² = 7.15 (dPhi = 13.18%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 5 ... chi² = 6.08 (dPhi = 14.59%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 6 ... chi² = 5.24 (dPhi = 13.37%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 7 ... chi² = 4.76 (dPhi = 8.56%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 8 ... chi² = 4.35 (dPhi = 7.67%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 9 ... chi² = 4.14 (dPhi = 4.32%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 10 ... chi² = 4.0 (dPhi = 3.18%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 11 ... chi² = 3.86 (dPhi = 3.0%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 12 ... chi² = 3.79 (dPhi = 1.6%) lam: 20.0
################################################################################
#                 Abort criteria reached: dPhi = 1.6 (< 2.0%)                  #
################################################################################

1204 [1127.028037028366,...,2611.3199864775293]

First have a look at the data fit. Plot the measured (crosses) and modelled (lines) traveltimes.

ax, cbar = mgr.showData(firstPicks=True, linewidth=0)
tt.drawFirstPicks(ax, data, mgr.inv.response, marker=None)
plot 04 koenigsee

Show resulting tomogram along with fit. You may want to save your results.

mgr.showResultAndFit()
mgr.saveResult()  # saves the results (mesh, velocity, vtk) in a folder
plot 04 koenigsee
'./20220829-12.29/TravelTimeManager'

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

ax, cbar = mgr.showResult(logScale=False, cMin=500, cMax=3000, cMap="plasma_r",
                          coverage=mgr.standardizedCoverage())
mgr.drawRayPaths(ax=ax, color="k", lw=0.3, alpha=0.5)

# mgr.coverage() yields the ray coverage in m and standardizedCoverage as 0/1
plot 04 koenigsee
<matplotlib.collections.LineCollection object at 0x7f33306c3760>

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 30.612 seconds)

Gallery generated by Sphinx-Gallery