Note

Go to the end to download the full example code

1D Rayleigh wave phase velocity inversion#

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)

What we do in this notebook#

Here we look at applying CoFI to the inversion of Rayleigh wave phase velocities for a 1D layered earth.

Learning outcomes

A demonstration of CoFI’s ability to switch between parameter estimation and ensemble methods.
A comparison between different McMC samplers that is fixed-d and trans-d samplers
An application of CoFI to field data

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

# !pip install -U cofi git+https://github.com/inlab-geo/pysurf96.git
# !git clone https://github.com/inlab-geo/cofi-examples.git
# %cd cofi-examples/examples/sw_rf_joint

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

from pysurf96 import surf96
import bayesbay
import cofi

Problem description#

Here we illustrate the range of inversion methods made avaialbe by CoFI. That is we first define a base problem and then explore the use of an iterative non linear apporach to find the MAP solution and then employ a range of Markov Chain Monte Carlo strategies to recover the posterior distribution. The forward problem is solved using pysurf 96 (miili/pysurf96) and the field data example is taken from (https://www.eas.slu.edu/eqc/eqc_cps/TUTORIAL/STRUCT/index.html) and we will be inverting observed rayleigh wave phase velocities

Inference problem

# display theory on the inference problem
from IPython.display import display, Markdown

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

display(Markdown(content))

<IPython.core.display.Markdown object>

Solving methods

# display theory on the optimisation approach
with open("../../theory/inv_optimisation.md", "r") as f:
    content = f.read()

display(Markdown(content))

<IPython.core.display.Markdown object>

# display theory on the optimisation approach
with open("../../theory/inv_mcmc.md", "r") as f:
    content = f.read()

display(Markdown(content))

<IPython.core.display.Markdown object>

Further reading

https://en.wikipedia.org/wiki/Surface_wave_inversion

Utilities#

1D model paramterisation#

# display theory on the 1D model parameterisation
with open("../../theory/misc_1d_model_parameterisation.md", "r") as f:
    content = f.read()

display(Markdown(content))

<IPython.core.display.Markdown object>

# layercake model utilities
def form_layercake_model(thicknesses, vs):
    model = np.zeros((len(vs)*2-1,))
    model[1::2] = thicknesses
    model[::2] = vs
    return model

def split_layercake_model(model):
    thicknesses = model[1::2]
    vs = model[::2]
    return thicknesses, vs

# voronoi model utilities
def form_voronoi_model(voronoi_sites, vs):
    return np.hstack((vs, voronoi_sites))

def split_voronoi_model(model):
    voronoi_sites = model[len(model)//2:]
    vs = model[:len(model)//2]
    return voronoi_sites, vs

def voronoi_to_layercake(voronoi_vector: np.ndarray) -> np.ndarray:
    n_layers = len(voronoi_vector) // 2
    velocities = voronoi_vector[:n_layers]
    voronoi_sites = voronoi_vector[n_layers:]
    depths = (voronoi_sites[:-1] + voronoi_sites[1:]) / 2
    thicknesses = depths - np.insert(depths[:-1], 0, 0)
    layercake_vector = np.zeros((2*n_layers-1,))
    layercake_vector[::2] = velocities
    layercake_vector[1::2] = thicknesses
    return layercake_vector

def layercake_to_voronoi(layercake_vector: np.ndarray, first_voronoi_site: float = 0.0) -> np.ndarray:
    n_layers = len(layercake_vector) // 2 + 1
    thicknesses = layercake_vector[1::2]
    velocities = layercake_vector[::2]
    depths = np.cumsum(thicknesses)
    voronoi_sites = np.zeros((n_layers,))
    for i in range(1,len(voronoi_sites)):
        voronoi_sites[i] = 2 * depths[i-1] - voronoi_sites[i-1]
    voronoi_vector = np.hstack((velocities, voronoi_sites))
    return voronoi_vector

Forward solver#

# display theory on the using the forward solver
with open("../../theory/geo_surface_wave_dispersion2.md", "r") as f:
    content = f.read()

display(Markdown(content))

<IPython.core.display.Markdown object>

# Constants
VP_VS = 1.77
RHO_VP_K = 0.32
RHO_VP_B = 0.77

# forward through pysurf96
def forward_sw(model, periods):
    thicknesses, vs = split_layercake_model(model)
    thicknesses = np.append(thicknesses, 10)
    vp = vs * VP_VS
    rho = RHO_VP_K * vp + RHO_VP_B
    return surf96(
        thicknesses,
        vp,
        vs,
        rho,
        periods,
        wave="rayleigh",
        mode=1,
        velocity="phase",
        flat_earth=False,
    )

# numerical jacobian
def jacobian_sw(model, periods, fwd=forward_sw, relative_step=0.01):
    jacobian = np.zeros((len(periods), len(model)))
    original_dpred = fwd(model, periods)
    for i in range(len(model)):
        perturbed_model = model.copy()
        step = relative_step * model[i]
        perturbed_model[i] += step
        perturbed_dpred = fwd(perturbed_model, periods)
        derivative = (perturbed_dpred - original_dpred) / abs(step)
        jacobian[:, i] = derivative
    return jacobian

Visualisation#

For conveninece we also implement two functions to plot the data here the Rayleigh wave phase velocity and a model given in the layer based parametrisation.

def plot_model(model, ax=None, title="model", **kwargs):
    # process data
    thicknesses = np.append(model[1::2], max(model[1::2]))
    velocities = model[::2]
    y = np.insert(np.cumsum(thicknesses), 0, 0)
    x = np.insert(velocities, 0, velocities[0])

    # plot depth profile
    if ax is None:
        _, ax = plt.subplots()
    plotting_style = {
        "linewidth": kwargs.pop("linewidth", kwargs.pop("lw", 0.5)),
        "alpha": 0.2,
        "color": kwargs.pop("color", kwargs.pop("c", "blue")),
    }
    plotting_style.update(kwargs)
    ax.step(x, y, where="post", **plotting_style)
    if ax.get_ylim()[0] < ax.get_ylim()[1]:
        ax.invert_yaxis()
    ax.set_xlabel("Vs (km/s)")
    ax.set_ylabel("Depth (km)")
    ax.set_title(title)
    return ax

def plot_data(rayleigh_phase_velocities, periods, ax=None, scatter=False,
              title="data", **kwargs):
    if ax is None:
        _, ax = plt.subplots()
    plotting_style = {
        "linewidth": kwargs.pop("linewidth", kwargs.pop("lw", 1)),
        "alpha": 1,
        "color": kwargs.pop("color", kwargs.pop("c", "blue")),
    }
    plotting_style.update(**kwargs)
    if scatter:
        ax.scatter(periods, rayleigh_phase_velocities, **plotting_style)
    else:
        ax.plot(periods, rayleigh_phase_velocities, **plotting_style)
    ax.set_xlabel("Periods (s)")
    ax.set_ylabel("Rayleigh phase velocities (km/s)")
    ax.set_title(title)
    return ax

def plot_model_and_data(model, label_model, color_model,
                        forward_func, periods, label_d_pred, color_d_pred,
                        axes=None, light_thin=False):
    if axes is None:
        _, axes = plt.subplots(1, 2, figsize=(10, 4), gridspec_kw={"width_ratios": [1, 2.5]})
    ax1, ax2 = axes
    if light_thin:
        plot_model(model, ax=ax1, color=color_model, alpha=0.2, lw=0.5, label=label_model)
        plot_data(forward_func(model, periods), periods, ax=ax2, color=color_d_pred, alpha=0.2, lw=0.5, label=label_d_pred)
    else:
        plot_model(model, ax=ax1, color=color_model, alpha=1, lw=1, label=label_model)
        plot_data(forward_func(model, periods), periods, ax=ax2, color=color_d_pred, label=label_d_pred)
    ax1.legend()
    ax2.legend()
    return ax1, ax2

Watermark#

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

cofi 0.2.8
numpy 1.26.4
matplotlib 3.8.3
scipy 1.12.0
emcee 3.1.5
bayesbay 0.3.0

sphinx_gallery_thumbnail_number = -1

Total running time of the script: (3 minutes 13.920 seconds)

Gallery generated by Sphinx-Gallery

1D Rayleigh wave phase velocity inversion

Contents

1D Rayleigh wave phase velocity inversion#

What we do in this notebook#

Problem description#

Utilities#

1D model paramterisation#

Forward solver#

Visualisation#

Synthetic example#

Generate synthetic data#

Optimisation#

Optimisation with no damping#

Optimal damping#

Optimisation with damping#

Fixed-dimensional sampling#

Sample the prior#

Sample the posterior#

Trans-dimensional sampling#

Field data example#

Read data#

Modified forward utility#

Optimisation#

Optimisation with no damping#

Optimal damping#

Optimisation with damping#

Fixed-dimensional sampling#

Sample the posterior#

Trans-dimensional sampling#

Watermark#