Seismic Wave Tomography via Fast Marching - Demo on switching regularization and L-curve

If you are running this notebook locally, make sure you’ve followed steps here to set up the environment. (This environment.yml file specifies a list of packages required to run the notebooks)

In this notebook, we would like to demonstrate the capability of CoFI to easily switch between different types of regularizations.

We will use cofi to run a seismic wave tomography example, in which the forward calculation is based on the Fast Marching Fortran code by Nick Rawlinson. The Fast Marching code is wrapped in package espresso.

We refer you to fmm_tomography.ipynb for further theretical details.

0. Import modules

# -------------------------------------------------------- #
#                                                          #
#     Uncomment below to set up environment on "colab"     #
#                                                          #
# -------------------------------------------------------- #

# !pip install -U cofi geo-espresso

import numpy as np
import matplotlib.pyplot as plt
import pprint

import cofi
import espresso

Understanding the inference problem

Before we starting working with cofi, let’s get familiar with the problem itself.

Below is a plot of the true model and the paths generated from this model. As you can see, there are two anomalies, one with lower velocity (red, top left) and the other with higher velocity (blue, bottom right).

fmm = espresso.FmmTomography()

fmm.plot_model(fmm.good_model, with_paths=True);

 New data set has:
 10  receivers
 10  sources
 100  travel times
 Range of travel times:  0.008911182496368759 0.0153757024856463
 Mean travel time: 0.01085811731230709

<Axes: xlabel='x (km)', ylabel='y (km)'>

pprint.pprint(fmm.metadata)

{'author_names': ['Nick Rawlinson', 'Malcolm Sambridge'],
 'citations': [('Rawlinson, N., de Kool, M. and Sambridge, M., 2006. Seismic '
                'wavefront tracking in 3-D heterogeneous media: applications '
                'with multiple data classes, Explor. Geophys., 37, 322-330.',
                ''),
               ('Rawlinson, N. and Urvoy, M., 2006. Simultaneous inversion of '
                'active and passive source datasets for 3-D seismic structure '
                'with application to Tasmania, Geophys. Res. Lett., 33 L24313',
                '10.1029/2006GL028105'),
               ('de Kool, M., Rawlinson, N. and Sambridge, M. 2006. A '
                'practical grid based method for tracking multiple refraction '
                'and reflection phases in 3D heterogeneous media, Geophys. J. '
                'Int., 167, 253-270',
                ''),
               ('Saygin, E. 2007. Seismic receiver and noise correlation based '
                'studies in Australia, PhD thesis, Australian National '
                'University.',
                '10.25911/5d7a2d1296f96')],
 'contact_email': 'Malcolm.Sambridge@anu.edu.au',
 'contact_name': 'Malcolm Sambridge',
 'linked_sites': [('Software published on iEarth',
                   'http://iearth.edu.au/codes/FMTOMO/')],
 'problem_short_description': 'The wave front tracker routines solves boundary '
                              'value ray tracing problems into 2D '
                              'heterogeneous wavespeed media, defined by '
                              'continuously varying velocity model calculated '
                              'by 2D cubic B-splines.',
 'problem_title': 'Fast Marching Wave Front Tracking'}

1. Problem setup and utilities

# get problem information from  espresso FmmTomography
model_size = fmm.model_size         # number of model parameters
model_shape = fmm.model_shape       # 2D spatial grids
data_size = fmm.data_size           # number of data points
ref_start_slowness = fmm.starting_model

def objective_func(slowness, reg):
    ttimes = fmm.forward(slowness)
    residual = fmm.data - ttimes
    data_misfit = residual.T @ residual
    model_reg = reg(slowness)
    return data_misfit + model_reg

def gradient(slowness, reg):
    ttimes, A = fmm.forward(slowness, return_jacobian=True)
    data_misfit_grad = -2 * A.T @ (fmm.data - ttimes)
    model_reg_grad = reg.gradient(slowness)
    return data_misfit_grad + model_reg_grad

def hessian(slowness, reg):
    A = fmm.jacobian(slowness)
    data_misfit_hess = 2 * A.T @ A
    model_reg_hess = reg.hessian(slowness)
    return data_misfit_hess + model_reg_hess

2. Invert with quadratic smoothing and damping regularization terms

2.1 Define BaseProblem

# define CoFI BaseProblem
fmm_problem_quadratic_reg = cofi.BaseProblem()
fmm_problem_quadratic_reg.set_initial_model(ref_start_slowness)

# add regularization: flattening + smoothing
smoothing_factor = 0.001
reg_smoothing = smoothing_factor * cofi.utils.QuadraticReg(
    model_shape=model_shape,
    weighting_matrix="smoothing"
)

fmm_problem_quadratic_reg.set_objective(objective_func, args=[reg_smoothing])
fmm_problem_quadratic_reg.set_gradient(gradient, args=[reg_smoothing])
fmm_problem_quadratic_reg.set_hessian(hessian, args=[reg_smoothing])

2.2 Define InversionOptions

my_options = cofi.InversionOptions()

my_options.set_tool("cofi.simple_newton")
my_options.set_params(
    num_iterations=15,
    step_length=1,
    obj_tol=1e-16,
    verbose=True,
    hessian_is_symmetric=True
)

2.3 Start an inversion

inv = cofi.Inversion(fmm_problem_quadratic_reg, my_options)
inv_result_quadratic_reg = inv.run()
inv_result_quadratic_reg.summary()

Iteration #0, updated objective function value: 1.7333503598417198e-07
Iteration #1, updated objective function value: 2.4531149751843008e-09
Iteration #2, updated objective function value: 1.666821376148691e-10
Iteration #3, updated objective function value: 3.40997782806854e-11
Iteration #4, updated objective function value: 3.8205147744763477e-11
Iteration #5, updated objective function value: 1.9124945329717177e-11
Iteration #6, updated objective function value: 3.215129076867847e-11
Iteration #7, updated objective function value: 1.603032618496361e-11
Iteration #8, updated objective function value: 2.891380116342641e-11
Iteration #9, updated objective function value: 1.506331785126499e-11
Iteration #10, updated objective function value: 2.7052434154813148e-11
Iteration #11, updated objective function value: 1.4201001080880201e-11
Iteration #12, updated objective function value: 2.6062305943114613e-11
Iteration #13, updated objective function value: 1.3931639754656179e-11
Iteration #14, updated objective function value: 2.5358538408239953e-11
============================
Summary for inversion result
============================
SUCCESS
----------------------------
model: [0.00050485 0.00049715 0.00048991 ... 0.00050558 0.00050198 0.00049801]
num_iterations: 14
objective_val: 2.5358538408239953e-11
n_obj_evaluations: 16
n_grad_evaluations: 15
n_hess_evaluations: 15

2.4 Plotting

clim = (1/np.max(fmm.good_model)-1, 1/np.min(fmm.good_model)+1)

fmm.plot_model(inv_result_quadratic_reg.model, clim=clim);            # inverted model
fmm.plot_model(fmm.good_model);       # true model

• 
• 

<Axes: xlabel='x (km)', ylabel='y (km)'>

---

3. Invert with Gaussian prior as regularization term

Instead of using a smoothing and damping regularization, in this section, we use a model covariance matrix and prior model.

\(\chi_{P}^{2}=\left(\mathbf{y} -\mathbf{f}(\mathbf{m})\right)^T C_d^{-1} \left(\mathbf{y} -\mathbf{f}(\mathbf{m})\right) + \left( \mathbf{m} - \mathbf{m}_p \right)^T C_p^{-1} \left( \mathbf{m} - \mathbf{m}_p \right)\)

\(\Delta \mathbf{m}= ({J}^T {C}_d^{-1} {J}+{C}_p^{-1})^{-1} ({J}^T{C}_d^{-1} (\mathbf{y}-\mathbf{f}(\mathbf{m}))+{C}_p^{-1}(\mathbf{m}_p-\mathbf{m}))\)

We can use CoFI’s utility module to help us generate a the Gaussian prior term.

3.1 Define BaseProblem

# define CoFI BaseProblem
fmm_problem_gaussian_prior = cofi.BaseProblem()
fmm_problem_gaussian_prior.set_initial_model(ref_start_slowness)

# add regularization: Gaussian prior
corrx = 3.0
corry = 3.0
sigma_slowness = 0.5**2
gaussian_prior = cofi.utils.GaussianPrior(
    model_covariance_inv=((corrx, corry), sigma_slowness),
    mean_model=ref_start_slowness.reshape(model_shape)
)

fmm_problem_gaussian_prior.set_objective(objective_func, args=[gaussian_prior])
fmm_problem_gaussian_prior.set_gradient(gradient, args=[gaussian_prior])
fmm_problem_gaussian_prior.set_hessian(hessian, args=[gaussian_prior])

3.2 Start an inversion

# reuse the previously defined InversionOptions object
inv = cofi.Inversion(fmm_problem_gaussian_prior, my_options)
inv_result_gaussian_prior = inv.run()
inv_result_gaussian_prior.summary()

Iteration #0, updated objective function value: 3.633219726912795e-07
Iteration #1, updated objective function value: 2.350718903302729e-07
Iteration #2, updated objective function value: 2.3157837862886104e-07
Iteration #3, updated objective function value: 2.3147517910024575e-07
Iteration #4, updated objective function value: 2.3140410377459331e-07
Change in model parameters below tolerance, stopping.
============================
Summary for inversion result
============================
SUCCESS
----------------------------
model: [0.00049703 0.00049595 0.00049451 ... 0.000503   0.00050221 0.00050161]
num_iterations: 4
objective_val: 2.3140410377459331e-07
n_obj_evaluations: 6
n_grad_evaluations: 5
n_hess_evaluations: 5

3.3 Plotting

fmm.plot_model(inv_result_gaussian_prior.model, clim=clim);            # gaussian prior
fmm.plot_model(fmm.good_model);       # true model

• 
• 

<Axes: xlabel='x (km)', ylabel='y (km)'>

4. L-curve

Now we plot an L-curve for the smoothing regularization case.

lambdas = np.logspace(-4, 4, 15)

my_lcurve_problems = []
for lamb in lambdas:
    my_reg = lamb * reg_smoothing
    my_problem = cofi.BaseProblem()
    my_problem.set_objective(objective_func, args=[my_reg])
    my_problem.set_gradient(gradient, args=[my_reg])
    my_problem.set_hessian(hessian, args=[my_reg])
    my_problem.set_initial_model(ref_start_slowness)
    my_lcurve_problems.append(my_problem)

my_options.set_params(verbose=False)

def my_callback(inv_result, i):
    m = inv_result.model
    res_norm = np.linalg.norm(fmm.forward(m) - fmm.data)
    reg_norm = np.sqrt(reg_smoothing(m))
    print(f"Finished inversion with lambda={lambdas[i]}: {res_norm}, {reg_norm}")
    return res_norm, reg_norm

my_inversion_pool = cofi.utils.InversionPool(
    my_lcurve_problems,
    my_options,
    my_callback,
    True
)
all_res, all_cb_returns = my_inversion_pool.run()

l_curve_points = list(zip(*all_cb_returns))

# plot the L-curve
res_norm, reg_norm = l_curve_points
plt.plot(reg_norm, res_norm, '.-')
plt.xlabel(r'Norm of regularization term $||Wm||_2$')
plt.ylabel(r'Norm of residual $||g(m)-d||_2$')
for i in range(len(lambdas)):
    plt.annotate(f'{lambdas[i]:.1e}', (reg_norm[i], res_norm[i]), fontsize=8)

# plot the previously solved model
my_inverted_model = inv_result_quadratic_reg.model
my_reg_norm = np.sqrt(reg_smoothing(my_inverted_model))
my_residual_norm = np.linalg.norm(fmm.forward(my_inverted_model) - fmm.data)
plt.plot(my_reg_norm, my_residual_norm, "x")
plt.annotate(f"{smoothing_factor:.1e}", (my_reg_norm, my_residual_norm), fontsize=8);

Text(2.7726666314754742e-06, 4.203671985174676e-06, '1.0e-03')

---

Watermark

watermark_list = ["cofi", "espresso", "numpy", "matplotlib"]
for pkg in watermark_list:
    pkg_var = __import__(pkg)
    print(pkg, getattr(pkg_var, "__version__"))

cofi 0.2.7
espresso 0.3.13
numpy 1.24.4
matplotlib 3.8.3

sphinx_gallery_thumbnail_number = -1