Fast Marching Method test using a two-layer modelΒΆ

This example shows how the FMM implementation of pyGIMLi works. FMM is an alternative to the Shortest path method and utilized by the Refraction manager. We illustrate the effectiveness and compare the solutions.

The first-arrival travel-time \(t\) from a shot point at the origin is described by a Hamilton-Jacobi (HJ) type equation, known as Eikonal equation [PL91]

\[(\nabla t)^2 = s^2\]

where \(s\) is the spatially distributed slowness (i.e. 1 / velocity). Under the assumption of piece-wise slowness (for individual cells) one can solve this equation by the so-called Fast-Marching method as first developed for 2D grids by Podvon & Lecompte (1991). Starting from the source location, a front is successively increased and filled by the travel times.

This example demonstrates the core features of an FMM Python implementation and compares the travel time for a two-layer model with the analytical solution (known for straight or sloped layers from the text books) and the shortest-path algorithm (Moser, 1991) implemented in the pyGIMLi core. [Mos91]

import time
from math import asin, tan

import numpy as np
import matplotlib.pyplot as plt

import pygimli as pg
from pygimli.mplviewer import drawMesh, drawField, drawStreamLines
from pygimli.physics.traveltime import fastMarch

First we provide the analytical solution for a given offset vector x.

def analyticalSolution2Layer(x, zlay, v1, v2):
    """Analytical solution: minimum of direct and critically refracted wave."""
    tdirect = np.abs(x) / v1  # direct wave
    alfa = asin(v1 / v2)  # critically refracted wave angle
    xreflec = tan(alfa) * zlay * 2.  # first critically refracted
    trefrac = (x - xreflec) / v2 + xreflec * v2 / v1**2
    return np.minimum(tdirect, trefrac)

The model consists of two boxes, in the first is the source (1)

0--1--------2
|           |
6-----------3
|           |
5-----------4

We create a PLC (piece-wise linear complex) and insert the nodes.

xmin, xmax, zlay = -20., 150., 20.  # model dimensions
plc = pg.Mesh(2)
nodes = []
nodes.append(plc.createNode(xmin, 0., 0.))  # 0
nodes.append(plc.createNode(0.0, 0., 0.))  # 1
nodes.append(plc.createNode(xmax, 0., 0.))  # 2
nodes.append(plc.createNode(xmax, -zlay, 0.))  # 3
nodes.append(plc.createNode(xmax, -zlay * 2, 0.))  # 4
nodes.append(plc.createNode(xmin, -zlay * 2, 0.))  # 5
nodes.append(plc.createNode(xmin, -zlay, 0.))  # 6

The nodes are connected from from 0 to 6 and back to 0. An additional edge is drawn from 6 to 3. Node/edge markers do not matter.

for i in range(6):
    plc.createEdge(nodes[i], nodes[i + 1])

plc.createEdge(nodes[6], nodes[0])
plc.createEdge(nodes[6], nodes[3])

We insert region markers (0 and 1) into the two layers and generate the mesh.

tri = pg.TriangleWrapper(plc)
plc.addRegionMarker(pg.RVector3(0., -zlay + .1), 0, 3.)  # 10m^2 max area
plc.addRegionMarker(pg.RVector3(0., -zlay - .1), 1, 10.)
tri.setSwitches('-pzeAfaq34.6')
mesh = pg.Mesh(2)
tri.generate(mesh)
mesh.createNeighbourInfos()
print(mesh)

Out:

Mesh: Nodes: 1781 Cells: 3405 Boundaries: 5185

Next we generate a velocity model from the markers by using a map. The values are associated to the markers and stored as attributes.

v = [1000., 3000.]
slomap = pg.stdMapF_F()  # mapping markers to real slowness values
for i, vi in enumerate(v):
    slomap.insert(i, 1. / vi)

mesh.mapCellAttributes(slomap)  # map values to attributes using map

We initialize the source position and the travel time vector initialize sets and tags and define the initial condition.

source = pg.RVector3(0., 0.)  # does not have to be a node!
times = pg.RVector(mesh.nodeCount(), 0.)
upwind, downwind = set(), set()
upTags, downTags = np.zeros(mesh.nodeCount()), np.zeros(mesh.nodeCount())
cell = mesh.findCell(source)
for i, n in enumerate(cell.nodes()):
    times[n.id()] = cell.attribute() * n.pos().distance(source)
    upTags[n.id()] = 1
for i, n in enumerate(cell.nodes()):
    tmpNodes = pg.commonNodes(n.cellSet())
    for nn in tmpNodes:
        if not upTags[nn.id()] and not downTags[nn.id()]:
            downwind.add(nn)
            downTags[nn.id()] = 1

Then we start marching until all fields are filled.

tic = time.time()
while len(downwind) > 0:
    fastMarch(mesh, downwind, times, upTags, downTags)

print(time.time() - tic, "s")

Out:

25.379729986190796 s

First, we plot the traveltime field and streamlines

fig, ax = plt.subplots(figsize=(10, 5))
drawMesh(ax, mesh)
ax.set_xlabel('x [m]')
ax.set_ylabel('y [m]')
drawField(ax, mesh, times, cmap='Spectral', fillContour=True)
drawStreamLines(ax, mesh, -times, nx=50, ny=50)
../../_images/sphx_glr_plot_2_FastMarching_2layer_001.png

We compare the result with the analytical solution along the x axis

x = np.arange(0., 140., 0.5)
tFMM = pg.interpolate(mesh, times, x, x * 0., x * 0.)
tAna = analyticalSolution2Layer(x, zlay, v[0], v[1])
print("min(dt)={} ms  max(dt)={} ms".format(min(tFMM - tAna) * 1000,
                                            max(tFMM - tAna) * 1000))

Out:

min(dt)=-0.0033639721106065723 ms  max(dt)=1.6065909149150792 ms

In order to use the Dijkstra, we extract the surface positions >0

mx = pg.x(mesh.positions())
my = pg.y(mesh.positions())
fi = pg.find((my == 0.0) & (mx >= 0))
px = np.sort(mx(fi))

A data container with index arrays named s (shot) and g (geophones) is created and filled with the positions and shot/geophone indices.

data = pg.DataContainer()
data.registerSensorIndex('s')
data.registerSensorIndex('g')
for pxi in px:
    data.createSensor(pg.RVector3(pxi, 0.0))

ndata = len(px) - 1
data.resize(ndata)
data.set('s', pg.RVector(ndata, 0))  # only one shot at first sensor
data.set('g', pg.utils.grange(1, ndata, 1))  # all others and geophones
fop = pg.TravelTimeDijkstraModelling(mesh, data)
tDijkstra = fop.response(mesh.cellAttributes())

We plot the calculated and measured travel times and relative differences

fig, ax = plt.subplots()
ax.plot(x, tAna*1000, 'r-', label='analytical')
ax.plot(x, tFMM*1000, 'b+-', label='FMM')
ax.plot(px[1:], tDijkstra*1000, 'gx-', label='Dijkstra')
ax.set_xlabel('x [m]')
ax.set_ylabel('t [s]')
ax.grid(True)
ax.legend()
ax2 = ax.twinx()
dtFMM = (tFMM - tAna) * 1000
tAnaD = analyticalSolution2Layer(px[1:], zlay, v[0], v[1])
dtDijkstra = (tDijkstra - tAnaD) * 1000
ax2.set_ylabel('dt [ms]')
if 1:  # relative differences in percent
    ax2.set_ylabel('dt [%]')
    dtFMM /= (tAna * 10)
    dtDijkstra /= (tAnaD * 10)
ax2.plot(x, dtFMM, 'b.-')
ax2.plot(px[1:], dtDijkstra, 'g.-')
pg.wait()
../../_images/sphx_glr_plot_2_FastMarching_2layer_002.png

Note that the Fast Marching Method is implemented in a modelling class (right now Python and very slow but to be replaced by fast C++ soon) FMModelling that can be more easily used with the Refraction manager by ra.useFMM(True)

Total running time of the script: ( 0 minutes 27.665 seconds)

Generated by Sphinx-Gallery

News (Aug 2017):
pyGIMLi Paper out now! Download PDF here.


© 2017 - GIMLi Development Team
Created using Bootstrap and Sphinx. Last updated on Nov 21, 2017.