Part II: Rayleigh and Love
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
LOVE_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)
rayleigh = pd(PERIODS, mode=0, wave="rayleigh").velocity
love = pd(PERIODS, mode=0, wave="love").velocity
rayleigh_obs = rayleigh + np.random.normal(0, RAYLEIGH_STD, rayleigh.size)
love_obs = love + np.random.normal(0, LOVE_STD, love.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, rayleigh, 'r--', label='Rayleigh (true Vs)', lw=0.5)
ax2.plot(PERIODS, love, 'b--', label='Love (true Vs)', lw=0.5)
ax2.plot(PERIODS, rayleigh_obs, 'ro', label='Noisy Rayleigh')
ax2.plot(PERIODS, love_obs, 'bo', label='Noisy Love')
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, wave='rayleigh', mode=0):
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=mode, wave=wave).velocity
return d_pred
from functools import partial
forward_rayleigh = partial(forward_sw, wave='rayleigh', mode=0)
forward_love = partial(forward_sw, wave='love', mode=0)
target_rayleigh = bb.Target("rayleigh",
rayleigh_obs,
std_min=0.001,
std_max=0.1,
std_perturb_std=0.002)
target_love = bb.Target("love",
love_obs,
std_min=0.001,
std_max=0.1,
std_perturb_std=0.002)
log_likelihood = bb.LogLikelihood(
targets=[target_rayleigh, target_love],
fwd_functions=[forward_rayleigh, forward_love]
)
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: 134242/350000 (38.35 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(['rayleigh', 'love']): 38580/49924 (77.28%)
ParamSpacePerturbation(param_space_name=voronoi):
BirthPerturbation(voronoi): 9071.00/49862.00 (18.19%)
DeathPerturbation(voronoi): 9069.00/50147.00 (18.08%)
ParamPerturbation(voronoi.discretization): 21900.00/49816.00 (43.96%)
ParamPerturbation(voronoi.vs): 55622.00/150251.00 (37.02%)
Chain ID: 1
TEMPERATURE: 1
EXPLORED MODELS: 350000
ACCEPTANCE RATE: 132036/350000 (37.72 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(['rayleigh', 'love']): 38667/50386 (76.74%)
ParamSpacePerturbation(param_space_name=voronoi):
BirthPerturbation(voronoi): 8864.00/49617.00 (17.86%)
DeathPerturbation(voronoi): 8853.00/49918.00 (17.74%)
ParamPerturbation(voronoi.discretization): 21730.00/50150.00 (43.33%)
ParamPerturbation(voronoi.vs): 53922.00/149929.00 (35.97%)
Chain ID: 2
TEMPERATURE: 1
EXPLORED MODELS: 350000
ACCEPTANCE RATE: 135593/350000 (38.74 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(['rayleigh', 'love']): 38776/50310 (77.07%)
ParamSpacePerturbation(param_space_name=voronoi):
BirthPerturbation(voronoi): 9249.00/50001.00 (18.50%)
DeathPerturbation(voronoi): 9244.00/49932.00 (18.51%)
ParamPerturbation(voronoi.discretization): 22447.00/50041.00 (44.86%)
ParamPerturbation(voronoi.vs): 55877.00/149716.00 (37.32%)
Chain ID: 3
TEMPERATURE: 1
EXPLORED MODELS: 350000
ACCEPTANCE RATE: 132019/350000 (37.72 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(['rayleigh', 'love']): 38382/49786 (77.09%)
ParamSpacePerturbation(param_space_name=voronoi):
BirthPerturbation(voronoi): 8783.00/50290.00 (17.46%)
DeathPerturbation(voronoi): 8779.00/49773.00 (17.64%)
ParamPerturbation(voronoi.discretization): 21902.00/50140.00 (43.68%)
ParamPerturbation(voronoi.vs): 54173.00/150011.00 (36.11%)
Chain ID: 4
TEMPERATURE: 1
EXPLORED MODELS: 350000
ACCEPTANCE RATE: 131004/350000 (37.43 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(['rayleigh', 'love']): 38689/50136 (77.17%)
ParamSpacePerturbation(param_space_name=voronoi):
BirthPerturbation(voronoi): 8820.00/49951.00 (17.66%)
DeathPerturbation(voronoi): 8809.00/49748.00 (17.71%)
ParamPerturbation(voronoi.discretization): 21806.00/49879.00 (43.72%)
ParamPerturbation(voronoi.vs): 52880.00/150286.00 (35.19%)
Chain ID: 5
TEMPERATURE: 1
EXPLORED MODELS: 350000
ACCEPTANCE RATE: 133402/350000 (38.11 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(['rayleigh', 'love']): 38876/50183 (77.47%)
ParamSpacePerturbation(param_space_name=voronoi):
BirthPerturbation(voronoi): 9127.00/50181.00 (18.19%)
DeathPerturbation(voronoi): 9121.00/50091.00 (18.21%)
ParamPerturbation(voronoi.discretization): 21883.00/49821.00 (43.92%)
ParamPerturbation(voronoi.vs): 54395.00/149724.00 (36.33%)
Chain ID: 6
TEMPERATURE: 1
EXPLORED MODELS: 350000
ACCEPTANCE RATE: 133889/350000 (38.25 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(['rayleigh', 'love']): 38602/49946 (77.29%)
ParamSpacePerturbation(param_space_name=voronoi):
BirthPerturbation(voronoi): 8960.00/50129.00 (17.87%)
DeathPerturbation(voronoi): 8954.00/50192.00 (17.84%)
ParamPerturbation(voronoi.discretization): 22200.00/49905.00 (44.48%)
ParamPerturbation(voronoi.vs): 55173.00/149828.00 (36.82%)
Chain ID: 7
TEMPERATURE: 1
EXPLORED MODELS: 350000
ACCEPTANCE RATE: 131188/350000 (37.48 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(['rayleigh', 'love']): 38543/50008 (77.07%)
ParamSpacePerturbation(param_space_name=voronoi):
BirthPerturbation(voronoi): 8760.00/49726.00 (17.62%)
DeathPerturbation(voronoi): 8757.00/49932.00 (17.54%)
ParamPerturbation(voronoi.discretization): 21798.00/49817.00 (43.76%)
ParamPerturbation(voronoi.vs): 53330.00/150517.00 (35.43%)
Chain ID: 8
TEMPERATURE: 1
EXPLORED MODELS: 350000
ACCEPTANCE RATE: 130058/350000 (37.16 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(['rayleigh', 'love']): 38619/49875 (77.43%)
ParamSpacePerturbation(param_space_name=voronoi):
BirthPerturbation(voronoi): 8854.00/49913.00 (17.74%)
DeathPerturbation(voronoi): 8845.00/49954.00 (17.71%)
ParamPerturbation(voronoi.discretization): 21914.00/49992.00 (43.84%)
ParamPerturbation(voronoi.vs): 51826.00/150266.00 (34.49%)
Chain ID: 9
TEMPERATURE: 1
EXPLORED MODELS: 350000
ACCEPTANCE RATE: 134029/350000 (38.29 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(['rayleigh', 'love']): 38844/50274 (77.26%)
ParamSpacePerturbation(param_space_name=voronoi):
BirthPerturbation(voronoi): 8986.00/49990.00 (17.98%)
DeathPerturbation(voronoi): 8984.00/49579.00 (18.12%)
ParamPerturbation(voronoi.discretization): 21971.00/50002.00 (43.94%)
ParamPerturbation(voronoi.vs): 55244.00/150155.00 (36.79%)
Chain ID: 10
TEMPERATURE: 1
EXPLORED MODELS: 350000
ACCEPTANCE RATE: 135972/350000 (38.85 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(['rayleigh', 'love']): 38148/49598 (76.91%)
ParamSpacePerturbation(param_space_name=voronoi):
BirthPerturbation(voronoi): 9363.00/50333.00 (18.60%)
DeathPerturbation(voronoi): 9364.00/49854.00 (18.78%)
ParamPerturbation(voronoi.discretization): 22596.00/50179.00 (45.03%)
ParamPerturbation(voronoi.vs): 56501.00/150036.00 (37.66%)
Chain ID: 11
TEMPERATURE: 1
EXPLORED MODELS: 350000
ACCEPTANCE RATE: 134231/350000 (38.35 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(['rayleigh', 'love']): 38495/49803 (77.29%)
ParamSpacePerturbation(param_space_name=voronoi):
BirthPerturbation(voronoi): 9103.00/50009.00 (18.20%)
DeathPerturbation(voronoi): 9100.00/49855.00 (18.25%)
ParamPerturbation(voronoi.discretization): 22342.00/50259.00 (44.45%)
ParamPerturbation(voronoi.vs): 55191.00/150074.00 (36.78%)
Chain ID: 12
TEMPERATURE: 1
EXPLORED MODELS: 350000
ACCEPTANCE RATE: 129166/350000 (36.90 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(['rayleigh', 'love']): 38785/50142 (77.35%)
ParamSpacePerturbation(param_space_name=voronoi):
BirthPerturbation(voronoi): 8643.00/50184.00 (17.22%)
DeathPerturbation(voronoi): 8644.00/50383.00 (17.16%)
ParamPerturbation(voronoi.discretization): 21510.00/49804.00 (43.19%)
ParamPerturbation(voronoi.vs): 51584.00/149487.00 (34.51%)
Chain ID: 13
TEMPERATURE: 1
EXPLORED MODELS: 350000
ACCEPTANCE RATE: 132291/350000 (37.80 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(['rayleigh', 'love']): 38469/49841 (77.18%)
ParamSpacePerturbation(param_space_name=voronoi):
BirthPerturbation(voronoi): 9016.00/49906.00 (18.07%)
DeathPerturbation(voronoi): 9013.00/50131.00 (17.98%)
ParamPerturbation(voronoi.discretization): 22135.00/50316.00 (43.99%)
ParamPerturbation(voronoi.vs): 53658.00/149806.00 (35.82%)
Chain ID: 14
TEMPERATURE: 1
EXPLORED MODELS: 350000
ACCEPTANCE RATE: 131594/350000 (37.60 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(['rayleigh', 'love']): 38394/49982 (76.82%)
ParamSpacePerturbation(param_space_name=voronoi):
BirthPerturbation(voronoi): 8844.00/49837.00 (17.75%)
DeathPerturbation(voronoi): 8842.00/50096.00 (17.65%)
ParamPerturbation(voronoi.discretization): 22021.00/50297.00 (43.78%)
ParamPerturbation(voronoi.vs): 53493.00/149788.00 (35.71%)
Chain ID: 15
TEMPERATURE: 1
EXPLORED MODELS: 350000
ACCEPTANCE RATE: 135836/350000 (38.81 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(['rayleigh', 'love']): 38797/50198 (77.29%)
ParamSpacePerturbation(param_space_name=voronoi):
BirthPerturbation(voronoi): 9045.00/49969.00 (18.10%)
DeathPerturbation(voronoi): 9042.00/49742.00 (18.18%)
ParamPerturbation(voronoi.discretization): 22442.00/50090.00 (44.80%)
ParamPerturbation(voronoi.vs): 56510.00/150001.00 (37.67%)
Chain ID: 16
TEMPERATURE: 1
EXPLORED MODELS: 350000
ACCEPTANCE RATE: 135742/350000 (38.78 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(['rayleigh', 'love']): 38718/50132 (77.23%)
ParamSpacePerturbation(param_space_name=voronoi):
BirthPerturbation(voronoi): 9187.00/50131.00 (18.33%)
DeathPerturbation(voronoi): 9184.00/49771.00 (18.45%)
ParamPerturbation(voronoi.discretization): 22473.00/49996.00 (44.95%)
ParamPerturbation(voronoi.vs): 56180.00/149970.00 (37.46%)
Chain ID: 17
TEMPERATURE: 1
EXPLORED MODELS: 350000
ACCEPTANCE RATE: 131447/350000 (37.56 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(['rayleigh', 'love']): 38702/50280 (76.97%)
ParamSpacePerturbation(param_space_name=voronoi):
BirthPerturbation(voronoi): 8960.00/49775.00 (18.00%)
DeathPerturbation(voronoi): 8956.00/49831.00 (17.97%)
ParamPerturbation(voronoi.discretization): 21340.00/49892.00 (42.77%)
ParamPerturbation(voronoi.vs): 53489.00/150222.00 (35.61%)
Chain ID: 18
TEMPERATURE: 1
EXPLORED MODELS: 350000
ACCEPTANCE RATE: 137481/350000 (39.28 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(['rayleigh', 'love']): 38762/49887 (77.70%)
ParamSpacePerturbation(param_space_name=voronoi):
BirthPerturbation(voronoi): 9335.00/49754.00 (18.76%)
DeathPerturbation(voronoi): 9330.00/49959.00 (18.68%)
ParamPerturbation(voronoi.discretization): 22314.00/49715.00 (44.88%)
ParamPerturbation(voronoi.vs): 57740.00/150685.00 (38.32%)
Chain ID: 19
TEMPERATURE: 1
EXPLORED MODELS: 350000
ACCEPTANCE RATE: 134230/350000 (38.35 %)
PARTIAL ACCEPTANCE RATES:
NoisePerturbation(['rayleigh', 'love']): 38929/50259 (77.46%)
ParamSpacePerturbation(param_space_name=voronoi):
BirthPerturbation(voronoi): 8983.00/50231.00 (17.88%)
DeathPerturbation(voronoi): 8978.00/50190.00 (17.89%)
ParamPerturbation(voronoi.discretization): 21945.00/49872.00 (44.00%)
ParamPerturbation(voronoi.vs): 55395.00/149448.00 (37.07%)
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=0.05)
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()
rayleigh_pred = np.array(results['rayleigh.dpred'])
love_pred = np.array(results['love.dpred'])
rayleigh_percentiles = np.percentile(rayleigh_pred, (10, 90), axis=0)
love_percentiles = np.percentile(love_pred, (10, 90), axis=0)
# Predicted dispersion curves
fig, ax = plt.subplots(1, 1, figsize=(7, 4))
ax.plot(PERIODS, rayleigh, 'r--', label='Rayleigh (True Vs)', lw=0.5)
ax.plot(PERIODS, rayleigh_obs, 'ro', label='Rayleigh Obs.')
ax.fill_between(PERIODS,
rayleigh_percentiles[0],
rayleigh_percentiles[1],
color='r', alpha=0.2,
label='Inferred Rayleigh (10th-90th perc.)',
zorder=100)
ax.plot(PERIODS, love, 'b--', label='Love (True Vs)', lw=0.5)
ax.plot(PERIODS, love_obs, 'bo', label='Love Obs.')
ax.fill_between(PERIODS,
love_percentiles[0],
love_percentiles[1],
color='b',
alpha=0.2,
label='Inferred Love (10th-90th perc.)',
zorder=100)
ax.set_xlabel('Period [s]')
ax.set_ylabel('Phase velocity [km/s]')
ax.grid()
ax.legend()
plt.show()
# Histograms of inferred data noise
std_min = target_rayleigh.std_min
std_max = target_rayleigh.std_max
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(6, 5))
ax1.axvline(x=RAYLEIGH_STD, color='r', lw=3, alpha=0.3, label='True data noise', zorder=0)
pdf, bins, _ = ax1.hist(results['rayleigh.std'], density=True, bins=20, ec='w', label='Posterior')
ax1.fill_between([std_min, std_max], 1 / (std_max - std_min), alpha=0.2, label='Prior')
ax1.set_xlabel('Noise standard deviation')
ax1.ticklabel_format(style='sci', axis='x', scilimits=(0,0))
ax1.legend(framealpha=0.9, loc='upper right')
ax1.set_title('Rayleigh')
ax2.axvline(x=LOVE_STD, color='r', lw=3, alpha=0.3, label='True data noise', zorder=0)
pdf, bins, _ = ax2.hist(results['love.std'], density=True, bins=20, ec='w', label='Posterior')
ax2.fill_between([std_min, std_max], 1 / (std_max - 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, loc='upper right')
ax2.set_title('Love')
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()
