# Traveltime tomography for Seismic Exploration

Authors
Affiliations
Francesco Picetti
Image and Sound Processing Lab, Politecnico di Milano
Ettore Biondi
SeismoLab, California Institute of Technology

# Eikonal traveltime tomography for seismic exploration scenarios¶

@Author: Ettore Biondi - ebiondi@caltech.edu

To use this notebook, we need pykonal installed. In your env, you can install it by:

1. conda install cython
2. pip install git+https://github.com/malcolmw/pykonal@0.2.3b3

In this notebook, we will apply the Eikonal tomography operators on a more realistic example using the Marmousi2 model. We will test two scenarios:

import numpy as np
import occamypy as o
import eikonal2d as k

# Plotting
from matplotlib import rcParams
from mpl_toolkits.axes_grid1 import make_axes_locatable
import matplotlib.pyplot as plt
rcParams.update({
'image.cmap'     : 'jet',
'image.aspect'   : 'equal',
'axes.grid'      : True,
'figure.figsize' : (19, 8),
'savefig.dpi'    : 300,
'axes.labelsize' : 14,
'axes.titlesize' : 16,
'font.size'      : 14,
'legend.fontsize': 14,
'xtick.labelsize': 14,
'ytick.labelsize': 14,
})

## Cross-well tomography¶

Let's plot the true model. Along with the geometry for the cross-well tomography example using a section of the Marmousi.

# Reading velocity model
nx = 1700
nz = 350
dx = 0.01
dz = 0.01
ox = 0.0
oz = 0.0
MarmTrue = np.reshape(np.fromfile("data/Marmousi2Vp.bin", dtype=">f"), (nx,nz)).astype(np.float32)

MarmTrue = MarmTrue[600:1200,:]

desampling=4
dx *= desampling
dz *= desampling
MarmTrue = MarmTrue[::desampling,::desampling]
nx,nz = MarmTrue.shape

# Velocity model
vvTrue = o.VectorNumpy(MarmTrue, ax_info=[o.AxInfo(nx, ox, dx, "x [km]"),
o.AxInfo(nz, oz, dz, "z [km]")])

extent = [vvTrue.ax_info[0].o, vvTrue.ax_info[0].last, vvTrue.ax_info[1].last, vvTrue.ax_info[1].o]
# Source/Receiver positions
SouPos = np.array([[int(nx*0.5),iz] for iz in np.arange(0,nz,2)])
RecPos = np.array([[0,iz] for iz in np.arange(0,nz,2)]+[[nx-1,iz] for iz in np.arange(0,nz,2)])
fig, ax = plt.subplots()
im = ax.imshow(vvTrue.plot().T, extent=extent)

ax.scatter(SouPos[:,0]*dx, SouPos[:,1]*dz,color="white", marker="*", s=100)
ax.scatter(RecPos[:,0]*dx, RecPos[:,1]*dz, color="tab:pink", marker="v",s=100)

plt.xlabel(vvTrue.ax_info[0].l)
plt.ylabel(vvTrue.ax_info[1].l)

cbar = plt.colorbar(im, orientation="horizontal", label="velocity [km/s]",

plt.tight_layout()
plt.show()

Now, we can proceed to generate the observed travel times.

# Data vector
tt_obs = o.VectorNumpy((SouPos.shape[0], RecPos.shape[0]),
ax_info=[o.AxInfo(SouPos.shape[0], 0, 1, "sou_id"),
o.AxInfo(RecPos.shape[0], 0, 1, "rec_id")]).zero()

# Setting Forward non-linear operator
Eik2D_Op = k.EikonalTT_2D(vvTrue, tt_obs, SouPos, RecPos)
Eik2D_Op.forward(False, vvTrue, tt_obs)
sou_id = -1

fig, ax = plt.subplots()

im = ax.imshow(Eik2D_Op.tt_maps[sou_id].T, extent=extent)

ax.scatter(SouPos[sou_id,0]*dx,SouPos[sou_id,1]*dz,color="white", marker="*", s=100)
ax.scatter(RecPos[:,0]*dx, RecPos[:,1]*dz, color="tab:pink", marker="v",s=100)

plt.xlabel(vvTrue.ax_info[0].l)
plt.ylabel(vvTrue.ax_info[1].l)

cbar = plt.colorbar(im, orientation="horizontal", label="traveltime [s]",

plt.tight_layout()
plt.show()
# Plotting traveltime vector
fig, ax = plt.subplots(figsize=(10,4))
ax.plot(tt_obs.ax_info[1].plot(), tt_obs[0,:], lw=2)
ax.plot(tt_obs.ax_info[1].plot(), tt_obs[15,:], lw=2)
ax.plot(tt_obs.ax_info[1].plot(), tt_obs[-1,:], lw=2)

plt.xlabel(tt_obs.ax_info[1].l)
plt.ylabel("Traveltime [s]")

ax.set_ylim(tt_obs.max()+1.0, 0)
ax.set_xlim(tt_obs.ax_info[1].o, tt_obs.ax_info[1].last)

plt.tight_layout()
plt.show()

The initial model will be constructed by smoothing the true model.

G = o.GaussianFilter(vvTrue, 2.5)
vvInit = G * vvTrue
fig, ax = plt.subplots()

im = ax.imshow(vvInit.plot().T, extent=extent, vmin=vvTrue.min(), vmax=vvTrue.max())

ax.scatter(SouPos[:,0]*dx, SouPos[:,1]*dz,color="white", marker="*", s=100)
ax.scatter(RecPos[:,0]*dx, RecPos[:,1]*dz, color="tab:pink", marker="v",s=100)

plt.xlabel(vvTrue.ax_info[0].l)
plt.ylabel(vvTrue.ax_info[1].l)

cbar = plt.colorbar(im, orientation="horizontal", label="velocity [km/s]",

plt.tight_layout()
plt.show()

Let's set up the inversion problem and solve it.

# linear operator
Eik2D_Lin_Op = k.EikonalTT_lin_2D(vvInit, tt_obs, SouPos, RecPos, tt_maps=Eik2D_Op.tt_maps)

# nonlinear operator
Eik2D_NlOp = o.NonlinearOperator(Eik2D_Op, Eik2D_Lin_Op, Eik2D_Lin_Op.set_vel)
# instantiate solver
BFGSsolver = o.LBFGS(o.BasicStopper(niter=20, tolg_proj=1e-32), m_steps=30)

G = o.GaussianFilter(vvInit, 3.5)

Eik2D_Inv_NlOp = o.NonlinearOperator(Eik2D_Op, Eik2D_Lin_Op * G, Eik2D_Lin_Op.set_vel)

L2_tt_prob = o.NonlinearLeastSquares(vvInit.clone(), tt_obs, Eik2D_Inv_NlOp,
minBound=vvInit.clone().set(1.0),  # min velocity: 1.0 km/s
maxBound=vvInit.clone().set(5.0))  # max velocity: 5.0 km/s
BFGSsolver.run(L2_tt_prob, verbose=True)
fig, ax = plt.subplots()

im = ax.imshow(L2_tt_prob.model.plot().T, extent=extent, vmin=vvTrue.min(), vmax=vvTrue.max())

ax.scatter(SouPos[:,0]*dx, SouPos[:,1]*dz,color="white", marker="*", s=100)
ax.scatter(RecPos[:,0]*dx, RecPos[:,1]*dz, color="tab:pink", marker="v",s=100)

plt.xlabel(vvTrue.ax_info[0].l)
plt.ylabel(vvTrue.ax_info[1].l)

cbar = plt.colorbar(im, orientation="horizontal", label="velocity [km/s]",

plt.tight_layout()
plt.show()

## First-break surface tomography¶

In this test we will be estimating the subsurface velocity using a surface acquisition.

# Reading velocity model
nx = 1700
nz = 350
dx = 0.01
dz = 0.01
ox = 0.0
oz = 0.0
MarmTrue = np.reshape(np.fromfile("data/Marmousi2Vp.bin", dtype=">f"), (nx,nz)).astype(np.float32)

desampling=4
dx *= desampling
dz *= desampling
MarmTrue = MarmTrue[::desampling,::desampling]
nx,nz = MarmTrue.shape

# Velocity model
vvTrue = o.VectorNumpy(MarmTrue, ax_info=[o.AxInfo(nx, 0, dx, "x [km]"),
o.AxInfo(nz, 0, dz, "z [km]")])

extent = [vvTrue.ax_info[0].o, vvTrue.ax_info[0].last, vvTrue.ax_info[1].last, vvTrue.ax_info[1].o]
# Source/Receiver positions
SouPos = np.array([[ix,0] for ix in np.arange(0,nx,10)])
RecPos = np.array([[ix,0] for ix in np.arange(0,nx,2)])
fig, ax = plt.subplots()

im = ax.imshow(vvTrue.plot().T, extent=extent)

ax.scatter(RecPos[:,0]*dx, RecPos[:,1]*dz, color="tab:pink", marker="v",s=100)
ax.scatter(SouPos[:,0]*dx, SouPos[:,1]*dz,color="white", marker="*", s=100)

plt.xlabel(vvTrue.ax_info[0].l)
plt.ylabel(vvTrue.ax_info[1].l)

cbar = plt.colorbar(im, orientation="horizontal", label="velocity [km/s]",

plt.tight_layout()
plt.show()
# Data vector
tt_obs = o.VectorNumpy((SouPos.shape[0], RecPos.shape[0]),
ax_info=[o.AxInfo(SouPos.shape[0], 0, 1, "sou_id"),
o.AxInfo(RecPos.shape[0], 0, 1, "rec_id")]).zero()

# Setting Forward non-linear operator
Eik2D_Op = k.EikonalTT_2D(vvTrue, tt_obs, SouPos, RecPos)
Eik2D_Op.forward(False, vvTrue, tt_obs)
sou_id = 20

fig, ax = plt.subplots()

im = ax.imshow(Eik2D_Op.tt_maps[sou_id].T, extent=extent)

ax.scatter(RecPos[:,0]*dx, RecPos[:,1]*dz, color="tab:pink", marker="v",s=100)
ax.scatter(SouPos[sou_id,0]*dx, SouPos[sou_id,1]*dz, color="white", marker="*", s=100)

plt.xlabel(vvTrue.ax_info[0].l)
plt.ylabel(vvTrue.ax_info[1].l)

cbar = plt.colorbar(im, orientation="horizontal", label="traveltime [s]",

plt.tight_layout()
plt.show()
# Plotting traveltime vector
fig, ax = plt.subplots(figsize=(10,4))
ax.plot(tt_obs.ax_info[1].plot(), tt_obs[0,:], lw=2)
ax.plot(tt_obs.ax_info[1].plot(), tt_obs[15,:], lw=2)
ax.plot(tt_obs.ax_info[1].plot(), tt_obs[-1,:], lw=2)

plt.xlabel(tt_obs.ax_info[1].l)
plt.ylabel("Traveltime [s]")

ax.set_ylim(tt_obs.max()+1.0, 0)
ax.set_xlim(tt_obs.ax_info[1].o, tt_obs.ax_info[1].last)

plt.tight_layout()
plt.show()

Again, let's create a initial smooth model

# Create Gaussian Filter
G = o.GaussianFilter(vvTrue, 2.5)
# Smoothing true model
vvInit = G * vvTrue
fig, ax = plt.subplots()

im = ax.imshow(vvInit.plot().T, extent=extent, vmin=vvTrue.min(), vmax=vvTrue.max())

ax.scatter(RecPos[:,0]*dx, RecPos[:,1]*dz, color="tab:pink", marker="v",s=100)
ax.scatter(SouPos[:,0]*dx, SouPos[:,1]*dz,color="white", marker="*", s=100)

plt.xlabel(vvTrue.ax_info[0].l)
plt.ylabel(vvTrue.ax_info[1].l)

cbar = plt.colorbar(im, orientation="horizontal", label="velocity [km/s]",

plt.tight_layout()
plt.show()
# Linear operator
Eik2D_Lin_Op = k.EikonalTT_lin_2D(vvInit, tt_obs, SouPos, RecPos, tt_maps=Eik2D_Op.tt_maps)
# Nonlinear operator
Eik2D_NlOp = o.NonlinearOperator(Eik2D_Op, Eik2D_Lin_Op, Eik2D_Lin_Op.set_vel)
# Instantiate solver
BFGSsolver = o.LBFGS(o.BasicStopper(niter=20, tolg_proj=1e-32))
# Create non-linear operator
Eik2D_Inv_NlOp = o.NonlinearOperator(Eik2D_Op, Eik2D_Lin_Op * G, Eik2D_Lin_Op.set_vel)

L2_tt_prob = o.NonlinearLeastSquares(vvInit.clone(), tt_obs, Eik2D_Inv_NlOp,
minBound=vvInit.clone().set(1.0),  # min velocity: 1.0 km/s
maxBound=vvInit.clone().set(5.0))  # max velocity: 5.0 km/s
BFGSsolver.run(L2_tt_prob, verbose=True)
fig, ax = plt.subplots()

im = ax.imshow(L2_tt_prob.model.plot().T, extent=extent, vmin=vvTrue.min(), vmax=vvTrue.max())

ax.scatter(SouPos[:,0]*dx, SouPos[:,1]*dz,color="white", marker="*", s=100)
ax.scatter(RecPos[:,0]*dx, RecPos[:,1]*dz, color="tab:pink", marker="v",s=100)

plt.xlabel(vvTrue.ax_info[0].l)
plt.ylabel(vvTrue.ax_info[1].l)

cbar = plt.colorbar(im, orientation="horizontal", label="velocity [km/s]",
plt.show()