There will be a webinar on pyGIMLi hosted by SEG on March 19, 2024 at 4 pm CET. Register for free here.

Source code for pygimli.physics.traveltime.refraction1d

from math import sqrt
import numpy as np
import pygimli as pg


def simulateNlayerRefraction(offsets, thk, vel, muteDirect=False):
    """First arrival of layered medium for one shot at the origin.

    TODO
    ----
    * Example
    * Test

    Parameters
    ----------
    offset : iterable
        Offsets distances between shot and geophone
    thk : iterable
        Thickness of the layers
    vel : iterable
        Velocities for the layers
    muteDirect : bool [False]
        Direct first arivals will be muted.

    Returns
    -------
    tt : np.array
        Traveltimes from origin to the offset positions.
    """
    s = 1. / np.array(vel)
    tt = offsets * s[0]
    
    if muteDirect:
        tt[:] = max(tt)

    nlay = len(vel)

    for i in range(1, nlay):  # i-th refracted
        tn = offsets * s[i]  # slope
        for j in range(i):  # sum over intercepts
            dsq = np.maximum(s[j]**2 - s[i]**2, 0)
            tn += 2 * thk[j] * sqrt(dsq)

        tt = np.minimum(tt, tn)
    return tt


[docs] class RefractionNLayer(pg.core.ModellingBase): """Forward operator for 1D Refraction seismic with layered model."""
[docs] def __init__(self, offset=0, nlay=3, vbase=1100, verbose=True): """Init forward operator. Model are velocity increases. Parameters ---------- offset : iterable vector of offsets between shot and geophone nlay : int number of layers vbase : float base velocity (all values are above) """ super().__init__(verbose=verbose) self.nlay = nlay self.offset = offset self.vbase = vbase mesh = pg.meshtools.createMesh1DBlock(nlay) self.setMesh(mesh)
[docs] def thkVel(self, model): """Return thickness and velocity vectors from model.""" return model(0, self.nlay-1), \ np.cumprod(model(self.nlay-1, self.nlay*2-1)) * self.vbase
[docs] def response(self, model): """Return forward response f(m).""" assert len(model) == self.nlay*2-1 return simulateNlayerRefraction(self.offset, *self.thkVel(model))
[docs] class RefractionNLayerFix1stLayer(pg.core.ModellingBase): """FOP for 1D Refraction seismic with layered model (e.g. water layer)."""
[docs] def __init__(self, offset=0, nlay=3, v0=1465, d0=200, muteDirect=False, verbose=True): """Init forward operator for velocity increases with fixed 1st layer. Parameters ---------- offset : iterable vector of offsets between shot and geophone nlay : int number of layers v0 : float First layer velocity (at the same time base velocity) d0 : float Depth of first layer in meter. muteDirect : bool [False] Mute the direct arrivels fron the first layer. """ super().__init__(verbose=verbose) self.nlay = nlay self.offset = offset self.v0 = v0 self.d0 = d0 self.mDirect = muteDirect mesh = pg.meshtools.createMesh1DBlock(nlay) self.setMesh(mesh)
[docs] def thkVel(self, model): """Return thickness and velocity vectors from model.""" thk = pg.cat(pg.Vector(1, self.d0), model(0, self.nlay-1)) relvel = model(self.nlay-1, self.nlay*2-1) vel = self.v0 * np.cumprod(pg.cat(pg.Vector(1, 1.0), relvel)) return thk, vel
[docs] def response(self, model): """Return forward response f(m).""" assert len(model) == self.nlay*2-1 return simulateNlayerRefraction(self.offset, *self.thkVel(model), muteDirect=self.mDirect)
if __name__ == '__main__': pass