Part I: Rayleigh
import bayesbay as bb
from bayesbay.discretization import Voronoi1D
import numpy as np
import matplotlib.pyplot as plt
from disba import PhaseDispersion
np.random.seed(30)
THICKNESS = np.array([10, 10, 15, 20, 20, 20, 20, 20, 0])
VS = np.array([3.38, 3.44, 3.66, 4.25, 4.35, 4.32, 4.315, 4.38, 4.5])
VP_VS = 1.77
VP = VS * VP_VS
RHO = 0.32 * VP + 0.77
PERIODS = np.geomspace(3, 80, 20)
RAYLEIGH_STD = 0.02
def initialize_vs(param, positions=None):
vmin, vmax = param.get_vmin_vmax(positions)
sorted_vals = np.sort(np.random.uniform(vmin, vmax, positions.size))
return sorted_vals
vs = bb.prior.UniformPrior(name="vs",
vmin=[2.2, 2.8, 3.3, 4],
vmax=[3.9, 4.6, 4.8, 5],
position=[0, 20, 60, 150],
perturb_std=0.15)
vs.set_custom_initialize(initialize_vs)
pd = PhaseDispersion(THICKNESS, VP, VS, RHO)
phase_vel = pd(PERIODS, mode=0, wave="rayleigh").velocity
d_obs = phase_vel + np.random.normal(0, RAYLEIGH_STD, phase_vel.size)
# For plotting the prior
depth_prior = np.linspace(0, 200)
vmin_prior, vmax_prior = vs.get_vmin_vmax(depth_prior)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(9, 4), gridspec_kw={'width_ratios': [1, 2.5]})
ax1.fill_betweenx(depth_prior, vmin_prior, vmax_prior, alpha=0.2, label='Prior')
Voronoi1D.plot_tessellation(THICKNESS,
VS,
label='True Vs',
ax=ax1,
color='k',
lw=2,
input_type='extents')
ax1.set_xlabel('Vs [km/s]')
ax1.set_ylabel('Depth [km]')
ax1.set_ylim(np.cumsum(THICKNESS)[-1] + max(THICKNESS), 0)
ax1.grid()
ax1.legend()
ax2.plot(PERIODS, phase_vel, 'r--', label='Rayleigh (true Vs)', lw=0.5)
ax2.plot(PERIODS, d_obs, 'ro', label='Rayleigh obs.')
ax2.set_xlabel('Period [s]')
ax2.set_ylabel('Phase velocity [km/s]')
ax2.grid()
ax2.legend()
plt.tight_layout(w_pad=2)
plt.show()
voronoi = Voronoi1D(
name="voronoi",
vmin=0,
vmax=150,
perturb_std=10,
n_dimensions=None,
n_dimensions_min=4,
n_dimensions_max=15,
parameters=[vs],
birth_from='neighbour'
)
parameterization = bb.parameterization.Parameterization(voronoi)
def forward_sw(state):
voronoi = state["voronoi"]
voronoi_sites = voronoi["discretization"]
thickness = Voronoi1D.compute_cell_extents(voronoi_sites)
vs = voronoi["vs"]
vp = vs * VP_VS
rho = 0.32 * vp + 0.77
pd = PhaseDispersion(thickness, vp, vs, rho)
d_pred = pd(PERIODS, mode=0, wave="rayleigh").velocity
return d_pred
target = bb.Target("rayleigh",
d_obs,
std_min=0.001,
std_max=0.1,
std_perturb_std=0.002)
log_likelihood = bb.LogLikelihood(targets=target, fwd_functions=forward_sw)
inversion = bb.BayesianInversion(
parameterization=parameterization,
log_likelihood=log_likelihood,
n_chains=20
)
inversion.run(
sampler=None,
n_iterations=350_000,
burnin_iterations=75_000,
save_every=150,
verbose=False
)
for chain in inversion.chains:
chain.print_statistics()
Chain ID: 0
TEMPERATURE: 1
EXPLORED MODELS: 350000
ACCEPTANCE RATE: 157575/350000 (45.02 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(rayleigh): 42611/49927 (85.35%)
ParamSpacePerturbation(param_space_name=voronoi):
BirthPerturbation(voronoi): 10045.00/49686.00 (20.22%)
DeathPerturbation(voronoi): 10038.00/50036.00 (20.06%)
ParamPerturbation(voronoi.discretization): 23767.00/49894.00 (47.63%)
ParamPerturbation(voronoi.vs): 71114.00/150457.00 (47.27%)
Chain ID: 1
TEMPERATURE: 1
EXPLORED MODELS: 350000
ACCEPTANCE RATE: 158312/350000 (45.23 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(rayleigh): 42762/49912 (85.67%)
ParamSpacePerturbation(param_space_name=voronoi):
BirthPerturbation(voronoi): 10253.00/49649.00 (20.65%)
DeathPerturbation(voronoi): 10245.00/50535.00 (20.27%)
ParamPerturbation(voronoi.discretization): 23743.00/49854.00 (47.63%)
ParamPerturbation(voronoi.vs): 71309.00/150050.00 (47.52%)
Chain ID: 2
TEMPERATURE: 1
EXPLORED MODELS: 350000
ACCEPTANCE RATE: 156074/350000 (44.59 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(rayleigh): 42757/50123 (85.30%)
ParamSpacePerturbation(param_space_name=voronoi):
BirthPerturbation(voronoi): 10095.00/49908.00 (20.23%)
DeathPerturbation(voronoi): 10091.00/50041.00 (20.17%)
ParamPerturbation(voronoi.discretization): 23725.00/50175.00 (47.28%)
ParamPerturbation(voronoi.vs): 69406.00/149753.00 (46.35%)
Chain ID: 3
TEMPERATURE: 1
EXPLORED MODELS: 350000
ACCEPTANCE RATE: 157139/350000 (44.90 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(rayleigh): 42713/49971 (85.48%)
ParamSpacePerturbation(param_space_name=voronoi):
BirthPerturbation(voronoi): 9963.00/49944.00 (19.95%)
DeathPerturbation(voronoi): 9960.00/49988.00 (19.92%)
ParamPerturbation(voronoi.discretization): 23490.00/50377.00 (46.63%)
ParamPerturbation(voronoi.vs): 71013.00/149720.00 (47.43%)
Chain ID: 4
TEMPERATURE: 1
EXPLORED MODELS: 350000
ACCEPTANCE RATE: 157877/350000 (45.11 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(rayleigh): 42717/50089 (85.28%)
ParamSpacePerturbation(param_space_name=voronoi):
BirthPerturbation(voronoi): 9938.00/49792.00 (19.96%)
DeathPerturbation(voronoi): 9929.00/49667.00 (19.99%)
ParamPerturbation(voronoi.discretization): 23948.00/50342.00 (47.57%)
ParamPerturbation(voronoi.vs): 71345.00/150110.00 (47.53%)
Chain ID: 5
TEMPERATURE: 1
EXPLORED MODELS: 350000
ACCEPTANCE RATE: 157169/350000 (44.91 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(rayleigh): 42668/50059 (85.24%)
ParamSpacePerturbation(param_space_name=voronoi):
BirthPerturbation(voronoi): 10315.00/50251.00 (20.53%)
DeathPerturbation(voronoi): 10310.00/49852.00 (20.68%)
ParamPerturbation(voronoi.discretization): 23940.00/50049.00 (47.83%)
ParamPerturbation(voronoi.vs): 69936.00/149789.00 (46.69%)
Chain ID: 6
TEMPERATURE: 1
EXPLORED MODELS: 350000
ACCEPTANCE RATE: 156444/350000 (44.70 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(rayleigh): 42234/49696 (84.98%)
ParamSpacePerturbation(param_space_name=voronoi):
BirthPerturbation(voronoi): 10014.00/49954.00 (20.05%)
DeathPerturbation(voronoi): 10016.00/50122.00 (19.98%)
ParamPerturbation(voronoi.discretization): 23593.00/50362.00 (46.85%)
ParamPerturbation(voronoi.vs): 70587.00/149866.00 (47.10%)
Chain ID: 7
TEMPERATURE: 1
EXPLORED MODELS: 350000
ACCEPTANCE RATE: 154587/350000 (44.17 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(rayleigh): 42702/50076 (85.27%)
ParamSpacePerturbation(param_space_name=voronoi):
BirthPerturbation(voronoi): 9982.00/49995.00 (19.97%)
DeathPerturbation(voronoi): 9977.00/49966.00 (19.97%)
ParamPerturbation(voronoi.discretization): 23187.00/50014.00 (46.36%)
ParamPerturbation(voronoi.vs): 68739.00/149949.00 (45.84%)
Chain ID: 8
TEMPERATURE: 1
EXPLORED MODELS: 350000
ACCEPTANCE RATE: 156155/350000 (44.62 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(rayleigh): 42842/50128 (85.47%)
ParamSpacePerturbation(param_space_name=voronoi):
BirthPerturbation(voronoi): 9893.00/49898.00 (19.83%)
DeathPerturbation(voronoi): 9884.00/49822.00 (19.84%)
ParamPerturbation(voronoi.discretization): 23483.00/50149.00 (46.83%)
ParamPerturbation(voronoi.vs): 70053.00/150003.00 (46.70%)
Chain ID: 9
TEMPERATURE: 1
EXPLORED MODELS: 350000
ACCEPTANCE RATE: 158024/350000 (45.15 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(rayleigh): 42519/49801 (85.38%)
ParamSpacePerturbation(param_space_name=voronoi):
BirthPerturbation(voronoi): 10198.00/49930.00 (20.42%)
DeathPerturbation(voronoi): 10198.00/50138.00 (20.34%)
ParamPerturbation(voronoi.discretization): 23329.00/50250.00 (46.43%)
ParamPerturbation(voronoi.vs): 71780.00/149881.00 (47.89%)
Chain ID: 10
TEMPERATURE: 1
EXPLORED MODELS: 350000
ACCEPTANCE RATE: 156907/350000 (44.83 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(rayleigh): 42610/50036 (85.16%)
ParamSpacePerturbation(param_space_name=voronoi):
BirthPerturbation(voronoi): 10195.00/49979.00 (20.40%)
DeathPerturbation(voronoi): 10187.00/49482.00 (20.59%)
ParamPerturbation(voronoi.discretization): 23983.00/50121.00 (47.85%)
ParamPerturbation(voronoi.vs): 69932.00/150382.00 (46.50%)
Chain ID: 11
TEMPERATURE: 1
EXPLORED MODELS: 350000
ACCEPTANCE RATE: 156273/350000 (44.65 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(rayleigh): 43059/50345 (85.53%)
ParamSpacePerturbation(param_space_name=voronoi):
BirthPerturbation(voronoi): 10246.00/50172.00 (20.42%)
DeathPerturbation(voronoi): 10241.00/50219.00 (20.39%)
ParamPerturbation(voronoi.discretization): 23557.00/49873.00 (47.23%)
ParamPerturbation(voronoi.vs): 69170.00/149391.00 (46.30%)
Chain ID: 12
TEMPERATURE: 1
EXPLORED MODELS: 350000
ACCEPTANCE RATE: 157690/350000 (45.05 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(rayleigh): 42874/50089 (85.60%)
ParamSpacePerturbation(param_space_name=voronoi):
BirthPerturbation(voronoi): 10175.00/50286.00 (20.23%)
DeathPerturbation(voronoi): 10165.00/49920.00 (20.36%)
ParamPerturbation(voronoi.discretization): 23344.00/49883.00 (46.80%)
ParamPerturbation(voronoi.vs): 71132.00/149822.00 (47.48%)
Chain ID: 13
TEMPERATURE: 1
EXPLORED MODELS: 350000
ACCEPTANCE RATE: 156574/350000 (44.74 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(rayleigh): 42724/50039 (85.38%)
ParamSpacePerturbation(param_space_name=voronoi):
BirthPerturbation(voronoi): 10157.00/50381.00 (20.16%)
DeathPerturbation(voronoi): 10150.00/49985.00 (20.31%)
ParamPerturbation(voronoi.discretization): 23531.00/49806.00 (47.25%)
ParamPerturbation(voronoi.vs): 70012.00/149789.00 (46.74%)
Chain ID: 14
TEMPERATURE: 1
EXPLORED MODELS: 350000
ACCEPTANCE RATE: 156945/350000 (44.84 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(rayleigh): 42858/50200 (85.37%)
ParamSpacePerturbation(param_space_name=voronoi):
BirthPerturbation(voronoi): 9993.00/50107.00 (19.94%)
DeathPerturbation(voronoi): 9985.00/49906.00 (20.01%)
ParamPerturbation(voronoi.discretization): 23571.00/49990.00 (47.15%)
ParamPerturbation(voronoi.vs): 70538.00/149797.00 (47.09%)
Chain ID: 15
TEMPERATURE: 1
EXPLORED MODELS: 350000
ACCEPTANCE RATE: 155274/350000 (44.36 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(rayleigh): 42797/50264 (85.14%)
ParamSpacePerturbation(param_space_name=voronoi):
BirthPerturbation(voronoi): 10126.00/49942.00 (20.28%)
DeathPerturbation(voronoi): 10127.00/50141.00 (20.20%)
ParamPerturbation(voronoi.discretization): 23597.00/49780.00 (47.40%)
ParamPerturbation(voronoi.vs): 68627.00/149873.00 (45.79%)
Chain ID: 16
TEMPERATURE: 1
EXPLORED MODELS: 350000
ACCEPTANCE RATE: 157937/350000 (45.12 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(rayleigh): 42835/50035 (85.61%)
ParamSpacePerturbation(param_space_name=voronoi):
BirthPerturbation(voronoi): 10291.00/50209.00 (20.50%)
DeathPerturbation(voronoi): 10288.00/49674.00 (20.71%)
ParamPerturbation(voronoi.discretization): 23801.00/50011.00 (47.59%)
ParamPerturbation(voronoi.vs): 70722.00/150071.00 (47.13%)
Chain ID: 17
TEMPERATURE: 1
EXPLORED MODELS: 350000
ACCEPTANCE RATE: 160715/350000 (45.92 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(rayleigh): 42967/50294 (85.43%)
ParamSpacePerturbation(param_space_name=voronoi):
BirthPerturbation(voronoi): 10320.00/50209.00 (20.55%)
DeathPerturbation(voronoi): 10318.00/49934.00 (20.66%)
ParamPerturbation(voronoi.discretization): 23699.00/49813.00 (47.58%)
ParamPerturbation(voronoi.vs): 73411.00/149750.00 (49.02%)
Chain ID: 18
TEMPERATURE: 1
EXPLORED MODELS: 350000
ACCEPTANCE RATE: 156671/350000 (44.76 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(rayleigh): 42509/49964 (85.08%)
ParamSpacePerturbation(param_space_name=voronoi):
BirthPerturbation(voronoi): 10131.00/50394.00 (20.10%)
DeathPerturbation(voronoi): 10123.00/49778.00 (20.34%)
ParamPerturbation(voronoi.discretization): 23119.00/49653.00 (46.56%)
ParamPerturbation(voronoi.vs): 70789.00/150211.00 (47.13%)
Chain ID: 19
TEMPERATURE: 1
EXPLORED MODELS: 350000
ACCEPTANCE RATE: 156817/350000 (44.80 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(rayleigh): 42613/49979 (85.26%)
ParamSpacePerturbation(param_space_name=voronoi):
BirthPerturbation(voronoi): 9960.00/49922.00 (19.95%)
DeathPerturbation(voronoi): 9963.00/50123.00 (19.88%)
ParamPerturbation(voronoi.discretization): 23339.00/49955.00 (46.72%)
ParamPerturbation(voronoi.vs): 70942.00/150021.00 (47.29%)
results = inversion.get_results(concatenate_chains=True)
sampled_voronoi_nuclei = results['voronoi.discretization']
sampled_thickness = [Voronoi1D.compute_cell_extents(n) for n in sampled_voronoi_nuclei]
sampled_vs = results['voronoi.vs']
interp_depths = np.linspace(0, 160, 160)
statistics_vs = Voronoi1D.get_tessellation_statistics(
sampled_thickness, sampled_vs, interp_depths, input_type='extents'
)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(6, 8), gridspec_kw={'width_ratios': [2.5, 1]})
ax1, cbar = Voronoi1D.plot_tessellation_density(sampled_thickness,
sampled_vs,
input_type='extents',
ax=ax1,
cmap='binary',
vmin=0,
vmax=2000)
Voronoi1D.plot_tessellation(THICKNESS,
VS,
ax=ax1,
color='yellow',
lw=3,
label='True Vs',
input_type='extents')
ax1.plot(statistics_vs['mean'], interp_depths, 'b', lw=2, label='Vs Ensemble Mean')
ax1.plot(statistics_vs['median'], interp_depths, 'r', lw=2, label='Vs Ensemble Median')
ax1.set_xlabel("Vs [km/s]")
ax1.set_ylabel("Depth [km]")
ax1.legend()
Voronoi1D.plot_interface_hist(sampled_voronoi_nuclei, bins=75, ec='w', ax=ax2)
ax2.tick_params(labelleft=False)
ax2.set_ylabel('')
ax2.set_ylim(*ax1.get_ylim())
plt.tight_layout()
plt.show()
d_pred = np.array(results['rayleigh.dpred'])
percentiles = np.percentile(d_pred, (10, 90), axis=0)
# Predicted dispersion curves
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(6, 6))
ax1.plot(PERIODS, phase_vel, 'r--', label='Rayleigh (true Vs)', lw=0.5)
ax1.plot(PERIODS, d_obs, 'ro', label='Observations')
ax1.fill_between(PERIODS, percentiles[0], percentiles[1], color='r', alpha=0.2, label='Inferred Data (10th-90th percentiles)', zorder=100)
ax1.set_xlabel('Period [s]')
ax1.set_ylabel('Phase velocity [km/s]')
ax1.grid()
ax1.legend()
# Histogram of inferred data noise
ax2.axvline(x=RAYLEIGH_STD, color='r', lw=3, alpha=0.3, label='True data noise')
pdf, bins, _ = ax2.hist(results['rayleigh.std'], density=True, bins=20, ec='w', zorder=100, label='Posterior')
ax2.fill_between([target.std_min, target.std_max], 1 / (target.std_max - target.std_min), alpha=0.2, label='Prior')
ax2.set_xlabel('Noise standard deviation')
ax2.ticklabel_format(style='sci', axis='x', scilimits=(0,0))
ax2.legend(framealpha=0.9)
plt.tight_layout()
plt.show()
def get_subplot_layout(n_subplots):
rows = int(np.sqrt(n_subplots))
cols = int(np.ceil(n_subplots / rows))
return rows, cols
rows, cols = get_subplot_layout(len(inversion.chains))
fig, axes = plt.subplots(rows, cols, figsize=(15, 15))
for ipanel, (ax, chain) in enumerate(zip(np.ravel(axes), inversion.chains)):
saved_states = chain.saved_states
saved_nuclei = saved_states["voronoi.discretization"]
saved_vs = saved_states['voronoi.vs']
Voronoi1D.plot_tessellations(saved_nuclei,
saved_vs,
ax=ax,
linewidth=0.1,
color="k",
bounds=(0, 150))
Voronoi1D.plot_tessellation(THICKNESS,
VS,
color='yellow',
lw=2,
ax=ax,
input_type='extents')
Voronoi1D.plot_tessellation_statistics(
saved_nuclei, saved_vs, interp_depths, ax=ax
)
ax.set_title(f'Chain {chain.id}')
ax.tick_params(direction='in', labelleft=False, labelbottom=False)
ax.set_xlabel('')
ax.set_ylabel('')
if not ipanel % cols:
ax.set_ylabel('Depth [km]')
ax.tick_params(labelleft=True)
if ipanel >= (rows-1) * cols:
ax.set_xlabel('Vs [km/s]')
ax.tick_params(labelbottom=True)
plt.tight_layout()
plt.show()
