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

@Author: Ettore Biondi - ettore88@stanford.edu

In [ ]:

# Eikonal Tomography module
import sys
pyEikoTomo = "../Scripts/"
sys.path.insert(0,pyEikoTomo)
pySolver = "../Scripts/python-solver/GenericSolver/python/"
sys.path.insert(0,pySolver)

import EikonalTomography2DMod
import pyVector

# Eikonal Tomography module import sys pyEikoTomo = "../Scripts/" sys.path.insert(0,pyEikoTomo) pySolver = "../Scripts/python-solver/GenericSolver/python/" sys.path.insert(0,pySolver) import EikonalTomography2DMod import pyVector

In [ ]:

import numpy as np
import EikonalTomography3DMod
import pyVector
import pyOperator as pyOp

# Plotting
from matplotlib import rcParams
from mpl_toolkits.axes_grid1 import make_axes_locatable
import matplotlib.pyplot as plt
fontsize = 24
rcParams.update({
    'image.cmap'     : 'jet_r',
    'image.aspect'   : 'auto',
    'axes.grid'      : False,
    'figure.figsize' : (20, 10),
    'savefig.dpi'    : 300,
    'axes.labelsize' : 14,
    'axes.titlesize' : 16,
    'font.size'      : 14,
    'legend.fontsize': 14,
    'xtick.labelsize': 14,
    'ytick.labelsize': 14
})

params = {
    'image.interpolation': 'nearest',
    'image.cmap': 'gray',
    'savefig.dpi': 300,  # to adjust notebook inline plot size
    'axes.labelsize': fontsize, # fontsize for x and y labels (was 10)
    'axes.titlesize': fontsize,
    'font.size': fontsize,
    'legend.fontsize': 18,
    'xtick.labelsize': fontsize,
    'ytick.labelsize': fontsize,
    "font.weight": 'bold',
    'text.usetex': False,
    "ytick.major.pad": 6.0,
    "xtick.major.pad": 6.0,
    "axes.labelweight": 'bold',
    'image.cmap'     : 'jet_r',
    'image.aspect'   : 'auto',
    'axes.grid'      : False
}
rcParams.update(params)

import numpy as np import EikonalTomography3DMod import pyVector import pyOperator as pyOp # Plotting from matplotlib import rcParams from mpl_toolkits.axes_grid1 import make_axes_locatable import matplotlib.pyplot as plt fontsize = 24 rcParams.update({ 'image.cmap' : 'jet_r', 'image.aspect' : 'auto', 'axes.grid' : False, 'figure.figsize' : (20, 10), 'savefig.dpi' : 300, 'axes.labelsize' : 14, 'axes.titlesize' : 16, 'font.size' : 14, 'legend.fontsize': 14, 'xtick.labelsize': 14, 'ytick.labelsize': 14 }) params = { 'image.interpolation': 'nearest', 'image.cmap': 'gray', 'savefig.dpi': 300, # to adjust notebook inline plot size 'axes.labelsize': fontsize, # fontsize for x and y labels (was 10) 'axes.titlesize': fontsize, 'font.size': fontsize, 'legend.fontsize': 18, 'xtick.labelsize': fontsize, 'ytick.labelsize': fontsize, "font.weight": 'bold', 'text.usetex': False, "ytick.major.pad": 6.0, "xtick.major.pad": 6.0, "axes.labelweight": 'bold', 'image.cmap' : 'jet_r', 'image.aspect' : 'auto', 'axes.grid' : False } rcParams.update(params)

In [ ]:

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

# Velocity model
vv = pyVector.vectorIC(np.zeros((1,nx,nz),dtype=np.float32))
vv.getNdArray()[:] = np.expand_dims(np.expand_dims(np.expand_dims(1. + 0.1*z, axis=0), axis=0), axis=0)

# Fast-Marching-Method (FMM) dx = dz = 0.1 nx, nz = 201, 101 oy = ox = oz = 0.0 x = np.linspace(ox,(nx-1)*dx, nx) z = np.linspace(oz,(nz-1)*dz, nz) # Velocity model vv = pyVector.vectorIC(np.zeros((1,nx,nz),dtype=np.float32)) vv.getNdArray()[:] = np.expand_dims(np.expand_dims(np.expand_dims(1. + 0.1*z, axis=0), axis=0), axis=0)

In [ ]:

# Source/Receiver positions on actual axis not in grid points!
SouPos = np.array([[x[0],z[nz-1]], [x[int(nx/2)], z[nz-1]], [x[nx-1], z[nz-1]]])
RecPos = np.array([[x[ix], z[0]] for ix in np.arange(0,nx-1)])

# Data vector
tt_data = pyVector.vectorIC(np.zeros((SouPos.shape[0], RecPos.shape[0]),dtype=np.float32))
tt_data.zero()

# Setting Forward non-linear operator
Eik2D_Op = EikonalTomography2DMod.EikonalTT_2D(vv, tt_data, ox, oz,
                                                            dx, dz,
                                                               SouPos, RecPos,
                                                               verbose=False, allocateTT=False)
Eik2D_Op.forward(False, vv, tt_data)

# Source/Receiver positions on actual axis not in grid points! SouPos = np.array([[x[0],z[nz-1]], [x[int(nx/2)], z[nz-1]], [x[nx-1], z[nz-1]]]) RecPos = np.array([[x[ix], z[0]] for ix in np.arange(0,nx-1)]) # Data vector tt_data = pyVector.vectorIC(np.zeros((SouPos.shape[0], RecPos.shape[0]),dtype=np.float32)) tt_data.zero() # Setting Forward non-linear operator Eik2D_Op = EikonalTomography2DMod.EikonalTT_2D(vv, tt_data, ox, oz, dx, dz, SouPos, RecPos, verbose=False, allocateTT=False) Eik2D_Op.forward(False, vv, tt_data)

In [ ]:

# Plotting traveltime vector
fig, ax = plt.subplots()
ax.plot(RecPos[:,0], tt_data.getNdArray()[0,:], lw=2)
ax.plot(RecPos[:,0], tt_data.getNdArray()[1,:], lw=2)
ax.plot(RecPos[:,0], tt_data.getNdArray()[2,:], lw=2)
ax.grid()
plt.ylabel("Traveltime [s]")
ax.set_ylim(0, 25)
ax.set_xlim(RecPos[:,0].min(), RecPos[:,0].max())
plt.tight_layout()
plt.show()

# Plotting traveltime vector fig, ax = plt.subplots() ax.plot(RecPos[:,0], tt_data.getNdArray()[0,:], lw=2) ax.plot(RecPos[:,0], tt_data.getNdArray()[1,:], lw=2) ax.plot(RecPos[:,0], tt_data.getNdArray()[2,:], lw=2) ax.grid() plt.ylabel("Traveltime [s]") ax.set_ylim(0, 25) ax.set_xlim(RecPos[:,0].min(), RecPos[:,0].max()) plt.tight_layout() plt.show()

In [ ]:

# Dot-product test for the linearized Eikonal equation in 1D
Eik2D_Lin_Op = EikonalTomography2DMod.EikonalTT_lin_2D(vv, tt_data, ox, oz,
                                                                dx, dz,
                                                               SouPos, RecPos,
                                                               verbose=False, ginsu=False)
Eik2D_Lin_Op.dotTest(True)

# Dot-product test for the linearized Eikonal equation in 1D Eik2D_Lin_Op = EikonalTomography2DMod.EikonalTT_lin_2D(vv, tt_data, ox, oz, dx, dz, SouPos, RecPos, verbose=False, ginsu=False) Eik2D_Lin_Op.dotTest(True)

Inversion¶

In [ ]:

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

# Velocity model
vv0 = pyVector.vectorIC(np.zeros((1,nx,nz),dtype=np.float32))

# Background Velocity model
vv0.getNdArray()[:] = np.expand_dims(np.expand_dims(1. + 0.1*z, axis=0), axis=0)

# 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.getNdArray()[:] -= np.expand_dims(gauss*0.5,axis=0)

# Fast-Marching-Method (FMM) # Fast-Marching-Method (FMM) dx = dz = 0.1 nx, nz = 201, 101 ox = oz = 0.0 x = np.linspace(ox,(nx-1)*dx, nx) z = np.linspace(oz,(nz-1)*dz, nz) # Velocity model vv0 = pyVector.vectorIC(np.zeros((1,nx,nz),dtype=np.float32)) # Background Velocity model vv0.getNdArray()[:] = np.expand_dims(np.expand_dims(1. + 0.1*z, axis=0), axis=0) # 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.getNdArray()[:] -= np.expand_dims(gauss*0.5,axis=0)

In [ ]:

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

# Data vector
tt_data = pyVector.vectorIC(np.zeros((SouPos.shape[0], RecPos.shape[0]),dtype=np.float32))
tt_data.zero()

# Setting Forward non-linear operator
Eik2D_Op = EikonalTomography2DMod.EikonalTT_2D(vv, tt_data, ox, oz,
                                                             dx, dz,
                                                               SouPos, RecPos,
                                                               verbose=False, allocateTT=False)

Eik2D_Lin_Op = EikonalTomography2DMod.EikonalTT_lin_2D(vv, tt_data, ox, oz,
                                                                dx, dz,
                                                               SouPos, RecPos,
                                                               verbose=False, ginsu=False)

# Source/Receiver positions SouPos = np.array([[x[ix],z[nz-1]] for ix in np.arange(0,nx,20)]) RecPos = np.array([[x[ix],z[0]] for ix in np.arange(0,nx-1)]) # Data vector tt_data = pyVector.vectorIC(np.zeros((SouPos.shape[0], RecPos.shape[0]),dtype=np.float32)) tt_data.zero() # Setting Forward non-linear operator Eik2D_Op = EikonalTomography2DMod.EikonalTT_2D(vv, tt_data, ox, oz, dx, dz, SouPos, RecPos, verbose=False, allocateTT=False) Eik2D_Lin_Op = EikonalTomography2DMod.EikonalTT_lin_2D(vv, tt_data, ox, oz, dx, dz, SouPos, RecPos, verbose=False, ginsu=False)

In [ ]:

fig, ax = plt.subplots()
im = ax.imshow(vv0.getNdArray()[0,:,:].T, extent=[x[0], x[-1], z[-1], z[0]], aspect=0.5)
ax = plt.gca()
ax.grid()
plt.xlabel("x [km]")
plt.ylabel("z [km]")
cbar = plt.colorbar(im, orientation="horizontal",
                    cax=make_axes_locatable(ax).append_axes("bottom", size="5%", pad=-.8))

_=ax.scatter(RecPos[:,0], RecPos[:,1]+0.2, marker="v", c="c", s=600.0, edgecolors= "k")
_=ax.scatter(SouPos[:,0], SouPos[:,1]-0.3, marker="*", c="r", s=1000.0, edgecolors= "k")
cbar.set_label('Velocity [km/s]')
plt.tight_layout()
plt.show()

fig, ax = plt.subplots() im = ax.imshow(vv0.getNdArray()[0,:,:].T, extent=[x[0], x[-1], z[-1], z[0]], aspect=0.5) ax = plt.gca() ax.grid() plt.xlabel("x [km]") plt.ylabel("z [km]") cbar = plt.colorbar(im, orientation="horizontal", cax=make_axes_locatable(ax).append_axes("bottom", size="5%", pad=-.8)) _=ax.scatter(RecPos[:,0], RecPos[:,1]+0.2, marker="v", c="c", s=600.0, edgecolors= "k") _=ax.scatter(SouPos[:,0], SouPos[:,1]-0.3, marker="*", c="r", s=1000.0, edgecolors= "k") cbar.set_label('Velocity [km/s]') plt.tight_layout() plt.show()

In [ ]:

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

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

In [ ]:

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

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

In [ ]:

# Plotting traveltime vector
fig, ax = plt.subplots(figsize=(20,10))
for ids in range(SouPos.shape[0]):
    ax.plot(RecPos[:,0], tt_data.getNdArray()[ids,:], lw=2)
ax.grid()
plt.xlabel('x [km]')
plt.ylabel("Traveltime [s]")
ax.set_ylim(5.0, 15.0)
ax.set_xlim(RecPos[:,0].min(), RecPos[:,0].max())
ax.invert_yaxis()
plt.tight_layout()
plt.show()

# Plotting traveltime vector fig, ax = plt.subplots(figsize=(20,10)) for ids in range(SouPos.shape[0]): ax.plot(RecPos[:,0], tt_data.getNdArray()[ids,:], lw=2) ax.grid() plt.xlabel('x [km]') plt.ylabel("Traveltime [s]") ax.set_ylim(5.0, 15.0) ax.set_xlim(RecPos[:,0].min(), RecPos[:,0].max()) ax.invert_yaxis() plt.tight_layout() plt.show()

In [ ]:

# instantiate solver
from pyNonLinearSolver import LBFGSsolver as LBFGS
from pyNpOperator import GaussianFilter
from pyProblem import ProblemL2NonLinear
from pyStopper import BasicStopper
BFGSsolver = LBFGS(BasicStopper(niter=50, tolg_proj=1e-32), m_steps=30)
# Saving objective function
BFGSsolver.setDefaults(save_obj=True)

# Creating problem object using Smoothing filter
G = GaussianFilter(vv0, 2.0)

Eik2D_Inv_NlOp = pyOp.NonLinearOperator(Eik2D_Op, Eik2D_Lin_Op * G, Eik2D_Lin_Op.set_vel)

L2_tt_prob = ProblemL2NonLinear(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

# instantiate solver from pyNonLinearSolver import LBFGSsolver as LBFGS from pyNpOperator import GaussianFilter from pyProblem import ProblemL2NonLinear from pyStopper import BasicStopper BFGSsolver = LBFGS(BasicStopper(niter=50, tolg_proj=1e-32), m_steps=30) # Saving objective function BFGSsolver.setDefaults(save_obj=True) # Creating problem object using Smoothing filter G = GaussianFilter(vv0, 2.0) Eik2D_Inv_NlOp = pyOp.NonLinearOperator(Eik2D_Op, Eik2D_Lin_Op * G, Eik2D_Lin_Op.set_vel) L2_tt_prob = ProblemL2NonLinear(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

In [ ]:

BFGSsolver.run(L2_tt_prob, verbose=True)

BFGSsolver.run(L2_tt_prob, verbose=True)

In [ ]:

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

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

In [ ]:

objVal = BFGSsolver.obj
fig = plt.figure(figsize=(20,6))
ax = plt.subplot(121)
plt.plot(objVal/objVal[0],lw=3)
ax.set_xlabel("Iteration number")
ax.set_ylabel("Scaled $\phi$")
ax.set_ylim([0,1])
ax.grid()
ax.autoscale(enable=True, axis='x', tight=True)
ax = plt.subplot(122)
ax.plot(np.log10(objVal/objVal[0]), lw=3)
ax.set_xlabel("Iteration number")
ax.set_ylabel("Log10 scaled $\phi$")
ax.set_ylim([None,0])
ax.grid()
ax.autoscale(enable=True, axis='x', tight=True)

objVal = BFGSsolver.obj fig = plt.figure(figsize=(20,6)) ax = plt.subplot(121) plt.plot(objVal/objVal[0],lw=3) ax.set_xlabel("Iteration number") ax.set_ylabel("Scaled $\phi$") ax.set_ylim([0,1]) ax.grid() ax.autoscale(enable=True, axis='x', tight=True) ax = plt.subplot(122) ax.plot(np.log10(objVal/objVal[0]), lw=3) ax.set_xlabel("Iteration number") ax.set_ylabel("Log10 scaled $\phi$") ax.set_ylim([None,0]) ax.grid() ax.autoscale(enable=True, axis='x', tight=True)

Try to compute invert the observed traveltimes with different level of smoothing and number of iterations to assess how the velocity anomaly is retrieved.

In [ ]:

### Your code ###

### Your code ###