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

Seismic Travel time 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 tomography example.

Theoretical background#

# display theory on travel time tomography
from IPython.display import display, Markdown

with open("../../theory/geo_travel_time_tomography.md", "r") as f:
    content = f.read()

display(Markdown(content))

<IPython.core.display.Markdown object>

For forward modelling, a fast marching wave front tracker is used, utilizing the Fast Marching Fortran code within the package `FMTOMO <http://iearth.edu.au/codes/FMTOMO/>`__ by Nick Rawlinson. The Fast Marching code is wrapped in package pyfm2d. Further details can be found in:

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.

0. Import modules#

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

# !pip install -U cofi pyfm2d

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

import cofi
# NB You will need to separately install pyfm2d in your python environment with `pip install pyfm2d'
import pyfm2d as wt # import fmm package

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).

# read in problem data
loaded_dict = np.load('../../data/travel_time_tomography/nonlinear_tomo_example.npz')
nonlinear_tomo_example = dict(loaded_dict)
loaded_dict.close()

# set up problem
good_model = nonlinear_tomo_example["_mtrue"]
extent = nonlinear_tomo_example["extent"]
sources = nonlinear_tomo_example["sources"]
receivers = nonlinear_tomo_example["receivers"]
obstimes = nonlinear_tomo_example["_data"]
print(' New data set have:\n',len(receivers),' receivers\n',len(sources),' sources\n',len(obstimes),' travel times\n',
'Range of travel times: ',np.min(obstimes),'to',np.max(obstimes),'\n Mean travel time:',np.mean(obstimes))

New data set have:
10  receivers
10  sources
100  travel times
Range of travel times:  0.009490006 to 0.01558705
Mean travel time: 0.011210016954999999

# display true model and raypaths
options = wt.WaveTrackerOptions(paths=True,cartesian=True) # set wavetracker options
result = wt.calc_wavefronts(good_model,receivers,sources,extent=extent, options=options) # track wavefronts
wt.display_model(good_model,paths=result.paths,extent=extent,line=0.3,alpha=0.82)

1. Problem setup and utilities#

# get problem information
model_size = good_model.size                           # number of model parameters
model_shape = good_model.shape                         # 2D spatial grid shape
data_size = data_size = len(obstimes)                  # number of data
ref_start_slowness = nonlinear_tomo_example["_sstart"] # use the starting guess supplied by the nonlinear example

def objective_func(slowness, reg, sigma, reduce_data=None):  # reduce_data=(idx_from, idx_to)
    if reduce_data is None: idx_from, idx_to = (0, data_size)
    else: idx_from, idx_to = reduce_data
    if(True):
        options = wt.WaveTrackerOptions(
            cartesian=True,
        )
        result = wt.calc_wavefronts(1./slowness.reshape(model_shape),receivers,sources,extent=extent,options=options) # track wavefronts
        ttimes = result.ttimes
    residual = obstimes[idx_from:idx_to] - ttimes[idx_from:idx_to]
    data_misfit = residual.T @ residual / sigma**2
    model_reg = reg(slowness)
    return  data_misfit + model_reg

def gradient(slowness, reg, sigma, reduce_data=None):       # reduce_data=(idx_from, idx_to)
    if reduce_data is None: idx_from, idx_to = (0, data_size)
    else: idx_from, idx_to = reduce_data
    if(True):
        options = wt.WaveTrackerOptions(
                    paths=True,
                    frechet=True,
                    cartesian=True,
                    )
        result = wt.calc_wavefronts(1./slowness.reshape(model_shape),receivers,sources,extent=extent,options=options) # track wavefronts
        ttimes = result.ttimes
        A = result.frechet.toarray()
    ttimes = ttimes[idx_from:idx_to]
    A = A[idx_from:idx_to]
    data_misfit_grad = -2 * A.T @ (obstimes[idx_from:idx_to] - ttimes) / sigma**2
    model_reg_grad = reg.gradient(slowness)
    return  data_misfit_grad + model_reg_grad

def hessian(slowness, reg, sigma, reduce_data=None):        # reduce_data=(idx_from, idx_to)
    if reduce_data is None: idx_from, idx_to = (0, data_size)
    else: idx_from, idx_to = reduce_data
    if(True):
        options = wt.WaveTrackerOptions(
                    paths=True,
                    frechet=True,
                    cartesian=True,
                    )
        result = wt.calc_wavefronts(1./slowness.reshape(model_shape),receivers,sources,extent=extent,options=options)
        A = result.frechet.toarray()
    A = A[idx_from:idx_to]
    data_misfit_hess = 2 * A.T @ A / sigma**2
    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.flatten())

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

sigma = 0.000008          # data standard deviation of noise
fmm_problem_quadratic_reg.set_objective(objective_func, args=[reg, sigma, None])
fmm_problem_quadratic_reg.set_gradient(gradient, args=[reg, sigma, None])
fmm_problem_quadratic_reg.set_hessian(hessian, args=[reg, sigma, None])

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: 37346.04444673342
Iteration #1, updated objective function value: 1496.5260409442917
Iteration #2, updated objective function value: 186.6540128296109
Iteration #3, updated objective function value: 49.72207315667848
Iteration #4, updated objective function value: 44.59552267851312
Iteration #5, updated objective function value: 87.36430757337139
Iteration #6, updated objective function value: 25.36158241794411
Iteration #7, updated objective function value: 36.87090297015299
Iteration #8, updated objective function value: 12.617214770325097
Iteration #9, updated objective function value: 13.857067302487787
Iteration #10, updated objective function value: 7.8961362057078714
Iteration #11, updated objective function value: 10.282043752198193
Iteration #12, updated objective function value: 14.08368705414384
Iteration #13, updated objective function value: 10.513448351701289
Iteration #14, updated objective function value: 5.3676081138584415
============================
Summary for inversion result
============================
SUCCESS
----------------------------
model: [0.00037221 0.00037399 0.00037577 ... 0.00061149 0.00061299 0.00061453]
num_iterations: 14
objective_val: 5.3676081138584415
n_obj_evaluations: 16
n_grad_evaluations: 15
n_hess_evaluations: 15

2.4 Plotting#

vmodel_inverted = 1./inv_result_quadratic_reg.model.reshape(model_shape)
wt.display_model(vmodel_inverted,extent=extent) # inverted model
wt.display_model(good_model,extent=extent) # true model

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.flatten())

# add regularization: Gaussian prior
corrx = 3.0
corry = 3.0
sigma_slowness = 0.5**2
sigma_slowness = 2.5E-6
gaussian_prior = 0.01 * 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, sigma])
fmm_problem_gaussian_prior.set_gradient(gradient, args=[gaussian_prior, sigma])
fmm_problem_gaussian_prior.set_hessian(hessian, args=[gaussian_prior, sigma])

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: 2248.1364063962765
Iteration #1, updated objective function value: 92.79030072251801
Iteration #2, updated objective function value: 28.63098462661843
Iteration #3, updated objective function value: 26.08983730748967
Iteration #4, updated objective function value: 26.279799379783867
Iteration #5, updated objective function value: 26.071267115197394
Iteration #6, updated objective function value: 26.298470633668593
Iteration #7, updated objective function value: 26.079080995248482
Iteration #8, updated objective function value: 26.309599744775284
Iteration #9, updated objective function value: 26.07521806350919
Iteration #10, updated objective function value: 26.310385039780268
Iteration #11, updated objective function value: 26.074474718579328
Iteration #12, updated objective function value: 26.31244771548592
Iteration #13, updated objective function value: 26.075721652167612
Iteration #14, updated objective function value: 26.31302076319902
============================
Summary for inversion result
============================
SUCCESS
----------------------------
model: [0.00049594 0.00049449 0.00049261 ... 0.00050232 0.00050181 0.00050139]
num_iterations: 14
objective_val: 26.31302076319902
n_obj_evaluations: 16
n_grad_evaluations: 15
n_hess_evaluations: 15

3.3 Plotting#

vmodel_inverted = 1./inv_result_gaussian_prior.model.reshape(model_shape)
wt.display_model(vmodel_inverted,extent=extent) # inverted model
wt.display_model(good_model,extent=extent) # true model

4. L-curve#

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

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

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, sigma])
    my_problem.set_gradient(gradient, args=[my_reg, sigma])
    my_problem.set_hessian(hessian, args=[my_reg, sigma])
    my_problem.set_initial_model(ref_start_slowness.flatten())
    my_lcurve_problems.append(my_problem)

my_options.set_params(verbose=False)

def my_callback(inv_result, i):
    m = inv_result.model
    slowness=m
    options = wt.WaveTrackerOptions(
            cartesian=True,
            )
    result = wt.calc_wavefronts(1./slowness.reshape(model_shape),receivers,sources,extent=extent,options=options) # track wavefronts
    ttimes = result.ttimes
    res_norm = np.linalg.norm(ttimes - obstimes)/sigma**2
    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,
    False
)
all_res, all_cb_returns = my_inversion_pool.run()

l_curve_points = list(zip(*all_cb_returns))

Finished inversion with lambda=0.0001: 211500.44428776175, 0.13483254428280728
Finished inversion with lambda=0.000774263682681127: 260731.1744337502, 0.13496174203337225
Finished inversion with lambda=0.005994842503189409: 271386.6695804191, 0.13502504598278992
Finished inversion with lambda=0.046415888336127774: 451370.3900748823, 0.1327427726565774
Finished inversion with lambda=0.3593813663804626: 453312.5962048154, 0.13102864469172437
Finished inversion with lambda=2.782559402207126: 264044.6349749823, 0.13441418768920918
Finished inversion with lambda=21.54434690031882: 618549.3299736021, 0.13133502515213039
Finished inversion with lambda=166.81005372000558: 723609.8417236926, 0.12674068778075587
Finished inversion with lambda=1291.5496650148827: 316054.01961665374, 0.11347506334192359
Finished inversion with lambda=10000.0: 503616.0297332419, 0.09505699526967784

# 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))
slowness=my_inverted_model
options = wt.WaveTrackerOptions(cartesian=True)
result = wt.calc_wavefronts(1./slowness.reshape(model_shape),receivers,sources,extent=extent,options=options) # track wavefronts
ttimes = result.ttimes
my_residual_norm = np.linalg.norm(ttimes - obstimes)/sigma**2
plt.plot(my_reg_norm, my_residual_norm, "x")
plt.annotate(f"{smoothing_factor:.1e}", (my_reg_norm, my_residual_norm), fontsize=8);

Text(0.1341853290376085, 289115.0941701654, '5.0e+06')

Watermark#

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

cofi 0.2.11+71.gb28b5b0
numpy 2.2.6
matplotlib 3.10.9

sphinx_gallery_thumbnail_number = -1

Total running time of the script: (2 minutes 55.467 seconds)

Gallery generated by Sphinx-Gallery