2D Travel-time tomography using the Eikonal equation

Authors

Affiliations

Francesco Picetti

Image and Sound Processing Lab, Politecnico di Milano

Ettore Biondi

SeismoLab, California Institute of Technology

April 26, 2022

Traveltime tomography using the Eikonal equation (2D case)¶

@Author: Ettore Biondi - ebiondi@caltech.edu

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

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

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,
})

WARNING! DATAPATH not found. The folder /tmp will be used to write binary files

/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/dask_jobqueue/core.py:20: FutureWarning: tmpfile is deprecated and will be removed in a future release. Please use dask.utils.tmpfile instead.
  from distributed.utils import tmpfile

# Fast-Marching-Method (FMM)
dx = dz = 0.1
nx, nz = 201, 101
x = np.linspace(0,(nx-1)*dx, nx)
z = np.linspace(0,(nz-1)*dz, nz)

# Velocity model
vv = o.VectorNumpy((nx,nz), ax_info=[o.AxInfo(nx, 0, dx, "x [km]"), o.AxInfo(nz, 0, dz, "z [km]")])
vv[:] = 1. + 0.1*z

# Source/Receiver positions
SouPos = np.array([[0,nz-1], [int(nx/2), nz-1], [nx-1, nz-1]])
RecPos = np.array([[ix,0] for ix in np.arange(0,nx)])

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

# Setting Forward non-linear operator
Eik2D_Op = k.EikonalTT_2D(vv, tt_data, SouPos, RecPos)
Eik2D_Op.forward(False, vv, tt_data)

# Plotting traveltime vector
fig, ax = plt.subplots()
ax.plot(tt_data.ax_info[1].plot(), tt_data[0], lw=2)
ax.plot(tt_data.ax_info[1].plot(), tt_data[1], lw=2)
ax.plot(tt_data.ax_info[1].plot(), tt_data[2], lw=2)
ax.grid()
plt.xlabel(tt_data.ax_info[1].l)
plt.ylabel("Traveltime [s]")
ax.set_ylim(0, 25)
ax.set_xlim(tt_data.ax_info[1].o, tt_data.ax_info[1].last)
plt.tight_layout()
plt.show()

# Dot-product test for the linearized Eikonal equation in 1D
Eik2D_Lin_Op = k.EikonalTT_lin_2D(vv, tt_data, SouPos, RecPos)
Eik2D_Lin_Op.dotTest(True)

Dot-product tests of forward and adjoint operators
--------------------------------------------------
Applying forward operator add=False
 Runs in: 2.6184580326080322 seconds
Applying adjoint operator add=False
 Runs in: 3.965895175933838 seconds
Dot products add=False: domain=1.542013e+00 range=1.542013e+00 
Absolute error: 8.881784e-16
Relative error: 5.759863e-16 

Applying forward operator add=True
 Runs in: 0.11392402648925781 seconds
Applying adjoint operator add=True
 Runs in: 0.514909029006958 seconds
Dot products add=True: domain=3.084026e+00 range=3.084026e+00 
Absolute error: 8.881784e-16
Relative error: 2.879932e-16 

-------------------------------------------------

Inversion¶

# Fast-Marching-Method (FMM)
# dx = dz = 0.05
dx = dz = 0.1
nx, nz = 201, 101
x = np.linspace(0,(nx-1)*dx, nx)
z = np.linspace(0,(nz-1)*dz, nz)

# Background Velocity model
vv0 = o.VectorNumpy((nx,nz), ax_info=[o.AxInfo(nx, 0, dx, "x [km]"), o.AxInfo(nz, 0, dz, "z [km]")])
vv0[:] = 1.0 + z*0.1

# Gaussian anomaly
zz, xx = np.meshgrid(z, x)
dst = np.sqrt(xx*xx+zz*zz)
sigma = 1.0
xloc = 101*dx
zloc = 51*dz
gauss = np.exp(-( ((xx-xloc)**2 + (zz-zloc)**2) / (2.0*sigma**2)))
# Constructing true model
vv = vv0.clone()
vv[:] -= gauss*0.5

fig, ax = plt.subplots()
im = ax.imshow(vv0.plot().T, extent=[x[0], x[-1], z[-1], z[0]], aspect=0.5, cmap="jet_r")
ax = plt.gca()
plt.xlabel(vv0.ax_info[0].l)
plt.ylabel(vv0.ax_info[1].l)
cbar = plt.colorbar(im, orientation="horizontal",
                    cax=make_axes_locatable(ax).append_axes("bottom", size="5%", pad=0))
cbar.set_label('Velocity [km/s]')
plt.tight_layout()
plt.show()
ax.grid()

fig, ax = plt.subplots()
im = ax.imshow(vv.plot().T, extent=[x[0], x[-1], z[-1], z[0]], aspect=0.5, cmap="jet_r")
ax = plt.gca()
plt.xlabel(vv0.ax_info[0].l)
plt.ylabel(vv0.ax_info[1].l)
cbar = plt.colorbar(im, orientation="horizontal",
                    cax=make_axes_locatable(ax).append_axes("bottom", size="5%", pad=0))
cbar.set_label('Velocity [km/s]')
plt.tight_layout()
plt.show()
ax.grid()

# Source/Receiver positions
SouPos = np.array([[ix,nz-1] for ix in np.arange(0,nx,20)])
RecPos = np.array([[ix,0] for ix in np.arange(0,nx)])

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

# Instantiating non-linear operator
Eik2D_Op = k.EikonalTT_2D(vv, tt_data, SouPos, RecPos)
Eik2D_Lin_Op = k.EikonalTT_lin_2D(vv, tt_data, SouPos, RecPos, tt_maps=Eik2D_Op.tt_maps)
Eik2D_NlOp = o.NonlinearOperator(Eik2D_Op, Eik2D_Lin_Op, Eik2D_Lin_Op.set_vel)

# Creating observed data
Eik2D_Op.forward(False, vv, tt_data)
tt_data_obs = tt_data.clone()

# Plotting traveltime vector
fig, ax = plt.subplots(figsize=(10,5))
for ids in range(SouPos.shape[0]):
    ax.plot(tt_data.ax_info[1].plot(), tt_data[ids], lw=2)
ax.grid()
plt.xlabel(tt_data.ax_info[1].l)
plt.ylabel("Traveltime [s]")
ax.set_ylim(5.0, 15.0)
ax.set_xlim(tt_data.ax_info[1].o, tt_data.ax_info[1].last)
ax.invert_yaxis()
plt.tight_layout()
plt.show()

# instantiate solver
BFGSBsolver = o.LBFGS(o.BasicStopper(niter=25, tolg_proj=1e-32), m_steps=30)

# Creating problem object using Smoothing filter
G = o.GaussianFilter(vv0, 1.5)

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

L2_tt_prob = o.NonlinearLeastSquares(vv0.clone(), tt_data_obs, Eik2D_Inv_NlOp, 
                                     minBound=vv0.clone().set(1.),  # min velocity: 1 km/s
                                     maxBound=vv0.clone().set(2.5)) # max velocity: 2.5 km/s

BFGSBsolver.run(L2_tt_prob, verbose=True)

##########################################################################################
Limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) algorithm 
Maximum number of steps to be used for Hessian inverse estimation: 30 
Restart folder: /tmp/restart_2022-04-26T08-05-57.869068/
##########################################################################################

iter = 00, obj = 2.48052e+01, resnorm = 7.04e+00, gradnorm = 8.73e+00, feval = 1, geval = 1
iter = 01, obj = 8.19519e+00, resnorm = 4.05e+00, gradnorm = 5.07e+00, feval = 2, geval = 2
iter = 02, obj = 2.58445e+00, resnorm = 2.27e+00, gradnorm = 1.52e+00, feval = 3, geval = 3
iter = 03, obj = 1.23840e+00, resnorm = 1.57e+00, gradnorm = 9.19e-01, feval = 4, geval = 4
iter = 04, obj = 3.15355e-01, resnorm = 7.94e-01, gradnorm = 4.07e-01, feval = 5, geval = 5
iter = 05, obj = 2.03645e-01, resnorm = 6.38e-01, gradnorm = 2.46e-01, feval = 6, geval = 6
iter = 06, obj = 9.79590e-02, resnorm = 4.43e-01, gradnorm = 1.75e-01, feval = 7, geval = 7
iter = 07, obj = 4.01127e-02, resnorm = 2.83e-01, gradnorm = 1.04e-01, feval = 8, geval = 8
iter = 08, obj = 2.53056e-02, resnorm = 2.25e-01, gradnorm = 7.52e-02, feval = 9, geval = 9
iter = 09, obj = 1.84552e-02, resnorm = 1.92e-01, gradnorm = 5.70e-02, feval = 10, geval = 10
iter = 10, obj = 1.35375e-02, resnorm = 1.65e-01, gradnorm = 3.87e-02, feval = 11, geval = 11
iter = 11, obj = 1.13803e-02, resnorm = 1.51e-01, gradnorm = 3.05e-02, feval = 12, geval = 12
iter = 12, obj = 9.72617e-03, resnorm = 1.39e-01, gradnorm = 2.81e-02, feval = 13, geval = 13
iter = 13, obj = 7.87789e-03, resnorm = 1.26e-01, gradnorm = 2.74e-02, feval = 14, geval = 14
iter = 14, obj = 5.84876e-03, resnorm = 1.08e-01, gradnorm = 2.24e-02, feval = 15, geval = 15
iter = 15, obj = 4.64104e-03, resnorm = 9.63e-02, gradnorm = 1.48e-02, feval = 16, geval = 16
iter = 16, obj = 4.15353e-03, resnorm = 9.11e-02, gradnorm = 1.24e-02, feval = 17, geval = 17
iter = 17, obj = 3.77232e-03, resnorm = 8.69e-02, gradnorm = 1.23e-02, feval = 18, geval = 18
iter = 18, obj = 3.22827e-03, resnorm = 8.04e-02, gradnorm = 1.25e-02, feval = 19, geval = 19
iter = 19, obj = 2.76429e-03, resnorm = 7.44e-02, gradnorm = 1.17e-02, feval = 20, geval = 20
iter = 20, obj = 2.43497e-03, resnorm = 6.98e-02, gradnorm = 1.11e-02, feval = 21, geval = 21
iter = 21, obj = 2.16737e-03, resnorm = 6.58e-02, gradnorm = 1.12e-02, feval = 22, geval = 22
iter = 22, obj = 1.89961e-03, resnorm = 6.16e-02, gradnorm = 1.12e-02, feval = 23, geval = 23
iter = 23, obj = 1.66281e-03, resnorm = 5.77e-02, gradnorm = 9.80e-03, feval = 24, geval = 24
iter = 24, obj = 1.49744e-03, resnorm = 5.47e-02, gradnorm = 8.04e-03, feval = 25, geval = 25
iter = 25, obj = 1.38691e-03, resnorm = 5.27e-02, gradnorm = 7.12e-03, feval = 26, geval = 26
Terminate: maximum number of iterations reached

##########################################################################################
Limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) algorithm end
##########################################################################################

fig, ax = plt.subplots()
im = ax.imshow(L2_tt_prob.model.plot().T, extent=[x[0], x[-1], z[-1], z[0]], aspect=0.5, cmap="jet_r")
ax = plt.gca()
plt.xlabel(vv0.ax_info[0].l)
plt.ylabel(vv0.ax_info[1].l)
cbar = plt.colorbar(im, orientation="horizontal",
                    cax=make_axes_locatable(ax).append_axes("bottom", size="5%", pad=0))
cbar.set_label('Velocity [km/s]')
plt.tight_layout()
plt.show()
ax.grid()

Notebooks

LS-RTM with devito: a use case

Notebooks

Rosenbrock: a nonlinear problem survey