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)

[[<matplotlib.lines.Line2D object at 0x7ff3d28b0340>], [<matplotlib.lines.Line2D object at 0x7ff3d59fe7c0>], [<matplotlib.lines.Line2D object at 0x7ff3d5b210a0>], [<matplotlib.lines.Line2D object at 0x7ff3d5b21220>], [<matplotlib.lines.Line2D object at 0x7ff3d291f400>], [<matplotlib.lines.Line2D object at 0x7ff3d29633a0>], [<matplotlib.lines.Line2D object at 0x7ff3d59e2700>], [<matplotlib.lines.Line2D object at 0x7ff3d59e2e20>], [<matplotlib.lines.Line2D object at 0x7ff3d54ae8b0>], [<matplotlib.lines.Line2D object at 0x7ff3d54aeb80>], [<matplotlib.lines.Line2D object at 0x7ff3d55042b0>], [<matplotlib.lines.Line2D object at 0x7ff3c59f5190>], [<matplotlib.lines.Line2D object at 0x7ff39f7453a0>], [<matplotlib.lines.Line2D object at 0x7ff3d5b08df0>], [<matplotlib.lines.Line2D object at 0x7ff3d2778910>]]

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()

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 0x7ff3d57b39a0>
Data transformation: <pygimli.core._pygimli_.RTrans object at 0x7ff3d54c7540>
Model transformation (cumulative):
         0 <pygimli.core._pygimli_.RTransLogLU object at 0x7ff3d59bd8e0>
min/max (data): 3.5e-04/0.03
min/max (error): 3%/3%
min/max (start model): 2.0e-04/0.002
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
inv.iter 1 ... chi² = 12.31 (dPhi = 91.5%) lam: 20
--------------------------------------------------------------------------------
inv.iter 2 ... chi² = 8.91 (dPhi = 27.36%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 3 ... chi² = 6.65 (dPhi = 24.91%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 4 ... chi² = 5.96 (dPhi = 9.76%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 5 ... chi² = 5.77 (dPhi = 3.1%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 6 ... chi² = 5.3 (dPhi = 7.86%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 7 ... chi² = 4.57 (dPhi = 12.58%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 8 ... chi² = 4.11 (dPhi = 9.0%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 9 ... chi² = 3.92 (dPhi = 4.26%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 10 ... chi² = 3.76 (dPhi = 3.6%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 11 ... chi² = 3.58 (dPhi = 4.07%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 12 ... chi² = 3.52 (dPhi = 1.53%) lam: 20.0
################################################################################
#                 Abort criteria reached: dPhi = 1.53 (< 2.0%)                 #
################################################################################

1090 [904.3153703700611,...,2661.055167733017]

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)

[[<matplotlib.lines.Line2D object at 0x7ff3d58315b0>], [<matplotlib.lines.Line2D object at 0x7ff3d5831cd0>], [<matplotlib.lines.Line2D object at 0x7ff3d583c340>], [<matplotlib.lines.Line2D object at 0x7ff3d583c9a0>], [<matplotlib.lines.Line2D object at 0x7ff3d584a0a0>], [<matplotlib.lines.Line2D object at 0x7ff3d584a6a0>], [<matplotlib.lines.Line2D object at 0x7ff3d584ad00>], [<matplotlib.lines.Line2D object at 0x7ff3d57d43a0>], [<matplotlib.lines.Line2D object at 0x7ff3d57d4a00>], [<matplotlib.lines.Line2D object at 0x7ff3d57e20a0>], [<matplotlib.lines.Line2D object at 0x7ff3d57e2700>], [<matplotlib.lines.Line2D object at 0x7ff3d57e2d60>], [<matplotlib.lines.Line2D object at 0x7ff3d57ec400>], [<matplotlib.lines.Line2D object at 0x7ff3d57eca60>], [<matplotlib.lines.Line2D object at 0x7ff3d57fa100>]]

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

'./20230707-18.10/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

<matplotlib.collections.LineCollection object at 0x7ff38693b430>

You can play around with the gradient starting model (vTop and vBottom arguments) and the regularization strength lam and customize the mesh.