Toy Problem: Surface Wave Inversion#
Notebook Outline:
Toy Surface-Wave Inversion Workflow (Python-native)#
Model parameterization (toy)#
We use a 1-D model with 3 layers over a halfspace, all horizontally layered.
Parameters:
thicknesses:
h1, h2, h3(meters)shear velocities:
Vs1, Vs2, Vs3, Vs4(m/s), whereVs4is the halfspacedensity is fixed (toy):
rho = 2000 kg/m^3everywhere
This is nonlinear because Love-wave dispersion depends on Vs(z) through a transcendental dispersion relation.
def love_dispersion_residual(c, omega, h, vs, rho):
'''
Residual of Love-wave dispersion for an N-layer stack over a halfspace.
We enforce traction-free surface and evanescent decay in the halfspace.
State vector: [u; tau] where u is SH displacement, tau is shear traction.
Propagate from surface (tau=0) downward to the halfspace interface.
In halfspace: tau = -mu * alpha * u, alpha = sqrt(k^2 - (omega/vs_half)^2).
Returns a real residual; root ~ 0 corresponds to a mode.
'''
h = np.asarray(h, dtype=float)
vs = np.asarray(vs, dtype=float)
N = len(h)
if np.isscalar(rho):
rho = np.full(N+1, float(rho))
else:
rho = np.asarray(rho, dtype=float)
k = omega / c # rad/m
# Surface BC: traction-free
u = 1.0 + 0.0j
tau = 0.0 + 0.0j
# Propagate layers (top -> bottom)
for i in range(N):
beta = vs[i]
mu = rho[i] * beta**2
q = np.sqrt((omega / beta)**2 - k**2 + 0j) # complex-safe
if abs(q) < 1e-12:
q = q + 1e-12j
ch = np.cos(q * h[i])
sh = np.sin(q * h[i])
u_new = u * ch + (tau / (mu * q)) * sh
tau_new = -mu * q * u * sh + tau * ch
u, tau = u_new, tau_new
# Halfspace decay condition
beta_h = vs[-1]
mu_h = rho[-1] * beta_h**2
alpha = np.sqrt(k**2 - (omega / beta_h)**2 + 0j)
if np.real(alpha) < 0:
alpha = -alpha
F = tau + mu_h * alpha * u
return float(np.real(F))
def love_phase_velocity_fundamental(freqs, h, vs, rho=2000.0, ngrid=600):
'''
Compute the *fundamental* Love-wave phase velocity curve c(f)
by bracketing and root finding at each frequency.
'''
freqs = np.asarray(freqs, dtype=float)
h = np.asarray(h, dtype=float)
vs = np.asarray(vs, dtype=float)
c_out = np.full_like(freqs, np.nan, dtype=float)
cmin = 0.95 * np.min(vs[:-1]) # slightly below slowest layer Vs
cmax = 0.999 * vs[-1] # slightly below halfspace Vs (guided)
for j, f in enumerate(freqs):
omega = 2*np.pi*f
cgrid = np.linspace(cmin, cmax, ngrid)
D = np.array([love_dispersion_residual(c, omega, h, vs, rho) for c in cgrid])
s = np.sign(D)
idx = np.where(s[:-1]*s[1:] < 0)[0]
if len(idx) == 0:
c_out[j] = cgrid[np.argmin(np.abs(D))]
continue
roots = []
for i0 in idx:
a, b = cgrid[i0], cgrid[i0+1]
try:
r = brentq(lambda c: love_dispersion_residual(c, omega, h, vs, rho), a, b, maxiter=200)
roots.append(r)
except ValueError:
pass
c_out[j] = np.min(roots) if roots else cgrid[np.argmin(np.abs(D))]
return c_out
Synthetic array seismograms from a dispersion curve (no PDE solve)#
We synthesize a dispersive wave packet with the target phase velocity curve c(f) by building
a complex spectrum at each receiver:
\[
U_x(f) = A(f)\,e^{-i 2\pi f t_0}\,e^{-i k(f) x},\qquad k(f)=\frac{2\pi f}{c(f)}
\]
Then irfft gives u(x,t). This is intentionally “physics-lite” but perfect for showing what an F–K plot does.
def synth_array_from_dispersion(x, dt, nt, f_model, c_model, f0=3.0, bw=1.0, t0=5.0, noise_std=0.02):
'''
Build dispersive array seismograms u(x,t) by frequency-domain synthesis.
x: receiver positions (m)
dt, nt: time sampling
f_model, c_model: dispersion curve for interpolation
'''
x = np.asarray(x, dtype=float)
nx = len(x)
t = np.arange(nt)*dt
f_fft = np.fft.rfftfreq(nt, dt)
c_fft = np.interp(f_fft, f_model, c_model, left=np.nan, right=np.nan)
A = np.exp(-0.5*((f_fft - f0)/bw)**2)
A[np.isnan(c_fft)] = 0.0
k = 2*np.pi*f_fft / np.where(c_fft>0, c_fft, np.nan)
k[np.isnan(k)] = 0.0
S = A * np.exp(-1j*2*np.pi*f_fft*t0)
u = np.zeros((nx, nt), dtype=float)
for i, xi in enumerate(x):
spec = S * np.exp(-1j*k*xi)
u[i, :] = np.fft.irfft(spec, n=nt)
u = u + noise_std*np.std(u)*np.random.randn(*u.shape)
return u, t
def compute_group_velocity(f, c):
'''
Compute group velocity U = dω/dk from phase velocity c(f).
U(f) = c(f) / (1 - (f/c) * dc/df)
'''
f = np.asarray(f)
c = np.asarray(c)
dcdf = np.gradient(c, f)
denom = 1.0 - (f / c) * dcdf
U = c / denom
U = np.clip(U, 0.5*np.min(c), 2.0*np.max(c))
return U
def plot_record_section(u, x, t, title="Record Section", tmin=None, tmax=None,
clip=3.0, show_wiggle=True, group_vel_overlay=None):
'''
Plot array seismograms u(x,t) as record section.
'''
u = np.asarray(u)
x = np.asarray(x)
t = np.asarray(t)
if tmin is None:
tmin = t[0]
if tmax is None:
tmax = t[-1]
tmask = (t >= tmin) & (t <= tmax)
t_plot = t[tmask]
u_plot = u[:, tmask]
std = np.std(u_plot)
u_norm = u_plot / (clip * std)
u_norm = np.clip(u_norm, -1, 1)
fig, ax = plt.subplots(figsize=(11, 6))
extent = [t_plot[0], t_plot[-1], x[0]/1000, x[-1]/1000]
im = ax.imshow(u_norm, aspect='auto', extent=extent, cmap='seismic',
vmin=-1, vmax=1, origin='lower', interpolation='bilinear')
if show_wiggle:
scale = 0.8 * (x[1] - x[0]) / 1000.0 / clip
for i in range(0, len(x), max(1, len(x)//20)):
trace = u_plot[i, :] / std
offset = x[i] / 1000.0
ax.plot(t_plot, offset + scale * trace, 'k-', lw=0.5, alpha=0.6)
if group_vel_overlay is not None:
f_band = group_vel_overlay['f']
U_band = group_vel_overlay['U']
t_start = group_vel_overlay.get('t_start', 0.0)
colors = plt.cm.plasma(np.linspace(0.1, 0.9, len(f_band)))
for i, (f, U) in enumerate(zip(f_band, U_band)):
t_arrival = t_start + x / U
ax.plot(t_arrival, x/1000, '--', color=colors[i], lw=2, alpha=0.8,
label=f'{f:.1f} Hz: U={U:.0f} m/s')
ax.legend(loc='upper left', fontsize=9, framealpha=0.8)
ax.set_xlabel('Time (s)', fontsize=12)
ax.set_ylabel('Distance (km)', fontsize=12)
ax.set_title(title, fontsize=13)
ax.grid(True, alpha=0.3, ls=':')
cbar = plt.colorbar(im, ax=ax, pad=0.01)
cbar.set_label('Normalized amplitude', fontsize=10)
plt.tight_layout()
plt.show()
F–K spectrum + ridge picking#
We compute a 2-D FFT across space and time. With the sign convention (u(x,t)\sim \cos(\omega t - kx)), energy concentrates along a ridge (k(\omega)). Convert to phase velocity:
\[
c(\omega)=\frac{\omega}{k}=\frac{2\pi f}{k}
\]
We then pick the ridge by taking, at each frequency, the wavenumber (k) of maximum amplitude (within a k-window).
def fk_spectrum(u, dx, dt):
'''
Compute |U(k,f)| from u(x,t) using a 2D FFT.
Returns:
f (Hz), k (rad/m), A (|U|) with f>=0 (half-spectrum).
'''
u = np.asarray(u, dtype=float)
nx, nt = u.shape
U = np.fft.fftshift(np.fft.fft2(u), axes=(0,1))
A = np.abs(U)
k = np.fft.fftshift(np.fft.fftfreq(nx, d=dx)) * 2*np.pi # rad/m
f = np.fft.fftshift(np.fft.fftfreq(nt, d=dt)) # Hz
pos = f >= 0
return f[pos], k, A[:, pos]
def pick_dispersion_ridge(f, k, A, fmin=0.5, fmax=10.0, kmin=0.0, kmax=None, smooth_kernel=9):
'''
Ridge pick: for each frequency, pick k at max amplitude in a k-window.
Convert to c = 2*pi*f/k.
'''
f = np.asarray(f)
k = np.asarray(k)
A = np.asarray(A)
if kmax is None:
kmax = np.max(k)
# positive k only
kpos = k > 0
k2 = k[kpos]
A2 = A[kpos, :]
fmask = (f >= fmin) & (f <= fmax)
f2 = f[fmask]
A3 = A2[:, fmask]
kmask = (k2 >= kmin) & (k2 <= kmax)
k3 = k2[kmask]
A4 = A3[kmask, :]
c_pick = np.full_like(f2, np.nan, dtype=float)
for j in range(len(f2)):
i_max = np.argmax(A4[:, j])
k_star = k3[i_max]
c_pick[j] = 2*np.pi*f2[j] / k_star
if smooth_kernel is not None and smooth_kernel >= 3 and smooth_kernel % 2 == 1:
c_pick = medfilt(c_pick, kernel_size=smooth_kernel)
return f2, c_pick
Demo: forward → synth seismograms → F–K → extracted dispersion#
# --- True model (3 layers + halfspace) ---
h_true = np.array([200.0, 500.0, 1000.0]) # m
vs_true = np.array([500.0, 900.0, 1400.0, 2500.0]) # m/s
# Frequencies for dispersion
f_model = np.linspace(0.5, 10.0, 40)
c_true = love_phase_velocity_fundamental(f_model, h_true, vs_true, rho=2000.0)
# Compute group velocity
U_true = compute_group_velocity(f_model, c_true)
# --- Synthesize array data with improved parameters ---
nx = 48
dx = 30.0 # m
x = np.arange(nx)*dx
dt = 0.01 # s
nt = 8192 # longer time window for dispersion
u, t = synth_array_from_dispersion(x, dt, nt, f_model, c_true,
f0=3.0, bw=2.5, t0=10.0, noise_std=0.02)
# --- Plot record section with group velocity overlay ---
f_overlay = np.array([1.5, 2.5, 3.5, 5.0])
U_overlay = np.interp(f_overlay, f_model, U_true)
plot_record_section(u, x, t,
title="Synthetic Dispersed Surface Waves (Love Mode)",
tmin=5, tmax=40, clip=2.5, show_wiggle=True,
group_vel_overlay={'f': f_overlay, 'U': U_overlay, 't_start': 10.0})
# --- F-K spectrum ---
f_fk, k_fk, A_fk = fk_spectrum(u, dx=dx, dt=dt)
# --- Pick ridge ---
f_obs, c_obs = pick_dispersion_ridge(f_fk, k_fk, A_fk, fmin=0.6, fmax=8.0, kmin=0.02, smooth_kernel=11)
# Plot both phase and group velocities
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))
ax1.plot(f_model, c_true, label="true c(f)", lw=2)
ax1.plot(f_obs, c_obs, "o", ms=4, label="picked from F-K")
ax1.set_xlabel("Frequency (Hz)")
ax1.set_ylabel("Phase velocity (m/s)")
ax1.set_title("Phase Velocity")
ax1.legend()
ax1.grid(True)
ax2.plot(f_model, U_true, label="true U(f)", lw=2, color='C1')
ax2.scatter(f_overlay, U_overlay, s=80, marker='*', zorder=10,
label="overlay freqs", edgecolor='k', linewidth=0.5)
ax2.set_xlabel("Frequency (Hz)")
ax2.set_ylabel("Group velocity (m/s)")
ax2.set_title("Group Velocity")
ax2.legend()
ax2.grid(True)
plt.tight_layout()
plt.show()
# Visualize the F-K amplitude (log scale)
plt.figure(figsize=(7,4))
plt.pcolormesh(f_fk, k_fk, np.log10(A_fk + 1e-12), shading="auto")
plt.xlabel("Frequency f (Hz)")
plt.ylabel("Wavenumber k (rad/m)")
plt.ylim(0, 0.35)
plt.title("F-K amplitude (log10)")
plt.colorbar(label="log10 |U|")
plt.show()
Inversion: fit c(f) with a 3-layer Vs(z)#
We define the model vector:
\[
m = [h_1,h_2,h_3, V_{S1},V_{S2},V_{S3},V_{S4}]
\]
Forward operator: (g(m) = c_{\text{pred}}(f))
Data vector: (d = c_{\text{obs}}(f))
Least-squares objective (with optional smoothness regularization):
\[
\phi(m)=\left\|\frac{g(m)-d}{\sigma}\right\|_2^2 + \lambda^2\left\|\Delta V_S\right\|_2^2
\]
We solve with scipy.optimize.least_squares using bounds to keep parameters physical.
def pack_params(h, vs):
return np.hstack([h, vs])
def unpack_params(p):
h = p[:3]
vs = p[3:]
return h, vs
sigma = 30.0 # m/s (toy constant uncertainty)
def residuals(p, freqs, c_obs, lam=0.0):
h, vs = unpack_params(p)
c_pred = love_phase_velocity_fundamental(freqs, h, vs, rho=2000.0, ngrid=500)
r_data = (c_pred - c_obs) / sigma
if lam > 0:
dv = np.diff(vs)
r_reg = lam * dv / 500.0 # scaling for conditioning (toy)
return np.hstack([r_data, r_reg])
return r_data
# Initial guess
h0 = np.array([300.0, 400.0, 800.0])
vs0 = np.array([600.0, 1000.0, 1600.0, 2300.0])
p0 = pack_params(h0, vs0)
lb = pack_params([50.0, 50.0, 50.0], [200.0, 200.0, 200.0, 500.0])
ub = pack_params([2000.0, 2000.0, 3000.0], [2000.0, 3000.0, 4000.0, 5000.0])
lam = 0.8
res = least_squares(residuals, p0, bounds=(lb, ub), args=(f_obs, c_obs, lam), verbose=1, max_nfev=40)
h_est, vs_est = unpack_params(res.x)
c_est = love_phase_velocity_fundamental(f_obs, h_est, vs_est, rho=2000.0)
print("Estimated thicknesses (m):", np.round(h_est, 1))
print("Estimated Vs (m/s):", np.round(vs_est, 1))
plt.figure()
plt.plot(f_obs, c_obs, ".", label="picked data")
plt.plot(f_model, c_true, label="true")
plt.plot(f_obs, c_est, "--", label="fit")
plt.xlabel("Frequency (Hz)")
plt.ylabel("Phase velocity (m/s)")
plt.legend()
plt.grid(True)
plt.show()
/var/folders/js/lzmy975n0l5bjbmr9db291m00000gn/T/ipykernel_74403/2599251142.py:50: RuntimeWarning: overflow encountered in scalar multiply
F = tau + mu_h * alpha * u
/var/folders/js/lzmy975n0l5bjbmr9db291m00000gn/T/ipykernel_74403/2599251142.py:39: RuntimeWarning: overflow encountered in scalar multiply
tau_new = -mu * q * u * sh + tau * ch
/var/folders/js/lzmy975n0l5bjbmr9db291m00000gn/T/ipykernel_74403/2599251142.py:50: RuntimeWarning: overflow encountered in scalar add
F = tau + mu_h * alpha * u
/var/folders/js/lzmy975n0l5bjbmr9db291m00000gn/T/ipykernel_74403/2599251142.py:39: RuntimeWarning: overflow encountered in scalar add
tau_new = -mu * q * u * sh + tau * ch
/var/folders/js/lzmy975n0l5bjbmr9db291m00000gn/T/ipykernel_74403/2599251142.py:38: RuntimeWarning: overflow encountered in scalar multiply
u_new = u * ch + (tau / (mu * q)) * sh
/var/folders/js/lzmy975n0l5bjbmr9db291m00000gn/T/ipykernel_74403/2599251142.py:50: RuntimeWarning: invalid value encountered in scalar multiply
F = tau + mu_h * alpha * u
/var/folders/js/lzmy975n0l5bjbmr9db291m00000gn/T/ipykernel_74403/2599251142.py:38: RuntimeWarning: overflow encountered in scalar add
u_new = u * ch + (tau / (mu * q)) * sh
---------------------------------------------------------------------------
KeyboardInterrupt Traceback (most recent call last)
Cell In[16], line 31
28 ub = pack_params([2000.0, 2000.0, 3000.0], [2000.0, 3000.0, 4000.0, 5000.0])
30 lam = 0.8
---> 31 res = least_squares(residuals, p0, bounds=(lb, ub), args=(f_obs, c_obs, lam), verbose=1, max_nfev=40)
33 h_est, vs_est = unpack_params(res.x)
34 c_est = love_phase_velocity_fundamental(f_obs, h_est, vs_est, rho=2000.0)
File ~/GitHub/ess-412-512-intro2seismology/.pixi/envs/default/lib/python3.11/site-packages/scipy/_lib/_util.py:696, in _workers_wrapper.<locals>.inner(*args, **kwds)
694 with MapWrapper(_workers) as mf:
695 kwargs['workers'] = mf
--> 696 return func(*args, **kwargs)
File ~/GitHub/ess-412-512-intro2seismology/.pixi/envs/default/lib/python3.11/site-packages/scipy/optimize/_lsq/least_squares.py:1019, in least_squares(fun, x0, jac, bounds, method, ftol, xtol, gtol, x_scale, loss, f_scale, diff_step, tr_solver, tr_options, jac_sparsity, max_nfev, verbose, args, kwargs, callback, workers)
1015 result = call_minpack(vector_fun.fun, x0, vector_fun.jac, ftol, xtol, gtol,
1016 max_nfev, x_scale, jac_method=jac)
1018 elif method == 'trf':
-> 1019 result = trf(vector_fun.fun, vector_fun.jac, x0, f0, J0, lb, ub, ftol, xtol,
1020 gtol, max_nfev, x_scale, loss_function, tr_solver,
1021 tr_options.copy(), verbose, callback=callback_wrapped)
1023 elif method == 'dogbox':
1024 if tr_solver == 'lsmr' and 'regularize' in tr_options:
File ~/GitHub/ess-412-512-intro2seismology/.pixi/envs/default/lib/python3.11/site-packages/scipy/optimize/_lsq/trf.py:124, in trf(fun, jac, x0, f0, J0, lb, ub, ftol, xtol, gtol, max_nfev, x_scale, loss_function, tr_solver, tr_options, verbose, callback)
120 return trf_no_bounds(
121 fun, jac, x0, f0, J0, ftol, xtol, gtol, max_nfev, x_scale,
122 loss_function, tr_solver, tr_options, verbose, callback=callback)
123 else:
--> 124 return trf_bounds(
125 fun, jac, x0, f0, J0, lb, ub, ftol, xtol, gtol, max_nfev, x_scale,
126 loss_function, tr_solver, tr_options, verbose, callback=callback)
File ~/GitHub/ess-412-512-intro2seismology/.pixi/envs/default/lib/python3.11/site-packages/scipy/optimize/_lsq/trf.py:376, in trf_bounds(fun, jac, x0, f0, J0, lb, ub, ftol, xtol, gtol, max_nfev, x_scale, loss_function, tr_solver, tr_options, verbose, callback)
372 f_true = f.copy()
374 cost = cost_new
--> 376 J = jac(x)
377 njev += 1
379 if loss_function is not None:
File ~/GitHub/ess-412-512-intro2seismology/.pixi/envs/default/lib/python3.11/site-packages/scipy/optimize/_differentiable_functions.py:741, in VectorFunction.jac(self, x)
739 def jac(self, x):
740 self._update_x(x)
--> 741 self._update_jac()
742 if hasattr(self.J, "astype"):
743 # returns a copy so that downstream can't overwrite the
744 # internal attribute. But one can't copy a LinearOperator
745 return self.J.astype(self.J.dtype)
File ~/GitHub/ess-412-512-intro2seismology/.pixi/envs/default/lib/python3.11/site-packages/scipy/optimize/_differentiable_functions.py:710, in VectorFunction._update_jac(self)
707 else:
708 self._njev += 1
--> 710 self.J = self.jac_wrapped(xp_copy(self.x), f0=self.f)
711 self.J_updated = True
File ~/GitHub/ess-412-512-intro2seismology/.pixi/envs/default/lib/python3.11/site-packages/scipy/optimize/_differentiable_functions.py:445, in _VectorJacWrapper.__call__(self, x, f0, **kwds)
443 self.njev += 1
444 elif self.jac in FD_METHODS:
--> 445 J, dct = approx_derivative(
446 self.fun,
447 x,
448 f0=f0,
449 **self.finite_diff_options,
450 )
451 self.nfev += dct['nfev']
453 if self.sparse_jacobian:
File ~/GitHub/ess-412-512-intro2seismology/.pixi/envs/default/lib/python3.11/site-packages/scipy/optimize/_numdiff.py:593, in approx_derivative(fun, x0, method, rel_step, abs_step, f0, bounds, sparsity, as_linear_operator, args, kwargs, full_output, workers)
591 with MapWrapper(workers) as mf:
592 if sparsity is None:
--> 593 J, _nfev = _dense_difference(fun_wrapped, x0, f0, h,
594 use_one_sided, method,
595 mf)
596 else:
597 if not issparse(sparsity) and len(sparsity) == 2:
File ~/GitHub/ess-412-512-intro2seismology/.pixi/envs/default/lib/python3.11/site-packages/scipy/optimize/_numdiff.py:686, in _dense_difference(fun, x0, f0, h, use_one_sided, method, workers)
684 f_evals = workers(fun, x_generator2(x0, h))
685 dx = [(x0[i] + h[i]) - x0[i] for i in range(n)]
--> 686 df = [f_eval - f0 for f_eval in f_evals]
687 df_dx = [delf / delx for delf, delx in zip(df, dx)]
688 nfev += len(df_dx)
File ~/GitHub/ess-412-512-intro2seismology/.pixi/envs/default/lib/python3.11/site-packages/scipy/optimize/_numdiff.py:686, in <listcomp>(.0)
684 f_evals = workers(fun, x_generator2(x0, h))
685 dx = [(x0[i] + h[i]) - x0[i] for i in range(n)]
--> 686 df = [f_eval - f0 for f_eval in f_evals]
687 df_dx = [delf / delx for delf, delx in zip(df, dx)]
688 nfev += len(df_dx)
File ~/GitHub/ess-412-512-intro2seismology/.pixi/envs/default/lib/python3.11/site-packages/scipy/optimize/_numdiff.py:879, in _Fun_Wrapper.__call__(self, x)
876 if xp.isdtype(x.dtype, "real floating"):
877 x = xp.astype(x, self.x0.dtype)
--> 879 f = np.atleast_1d(self.fun(x, *self.args, **self.kwargs))
880 if f.ndim > 1:
881 raise RuntimeError("`fun` return value has "
882 "more than 1 dimension.")
File ~/GitHub/ess-412-512-intro2seismology/.pixi/envs/default/lib/python3.11/site-packages/scipy/optimize/_differentiable_functions.py:415, in _VectorFunWrapper.__call__(self, x)
413 def __call__(self, x):
414 self.nfev += 1
--> 415 return np.atleast_1d(self.fun(x))
File ~/GitHub/ess-412-512-intro2seismology/.pixi/envs/default/lib/python3.11/site-packages/scipy/optimize/_lsq/least_squares.py:263, in _WrapArgsKwargs.__call__(self, x)
262 def __call__(self, x):
--> 263 return self.f(x, *self.args, **self.kwargs)
Cell In[16], line 13, in residuals(p, freqs, c_obs, lam)
11 def residuals(p, freqs, c_obs, lam=0.0):
12 h, vs = unpack_params(p)
---> 13 c_pred = love_phase_velocity_fundamental(freqs, h, vs, rho=2000.0, ngrid=500)
14 r_data = (c_pred - c_obs) / sigma
16 if lam > 0:
Cell In[10], line 85, in love_phase_velocity_fundamental(freqs, h, vs, rho, ngrid)
83 a, b = cgrid[i0], cgrid[i0+1]
84 try:
---> 85 r = brentq(lambda c: love_dispersion_residual(c, omega, h, vs, rho), a, b, maxiter=200)
86 roots.append(r)
87 except ValueError:
File ~/GitHub/ess-412-512-intro2seismology/.pixi/envs/default/lib/python3.11/site-packages/scipy/optimize/_zeros_py.py:846, in brentq(f, a, b, args, xtol, rtol, maxiter, full_output, disp)
844 raise ValueError(f"rtol too small ({rtol:g} < {_rtol:g})")
845 f = _wrap_nan_raise(f)
--> 846 r = _zeros._brentq(f, a, b, xtol, rtol, maxiter, args, full_output, disp)
847 return results_c(full_output, r, "brentq")
File ~/GitHub/ess-412-512-intro2seismology/.pixi/envs/default/lib/python3.11/site-packages/scipy/optimize/_zeros_py.py:94, in _wrap_nan_raise.<locals>.f_raise(x, *args)
93 def f_raise(x, *args):
---> 94 fx = f(x, *args)
95 f_raise._function_calls += 1
96 if np.isnan(fx):
Cell In[10], line 85, in love_phase_velocity_fundamental.<locals>.<lambda>(c)
83 a, b = cgrid[i0], cgrid[i0+1]
84 try:
---> 85 r = brentq(lambda c: love_dispersion_residual(c, omega, h, vs, rho), a, b, maxiter=200)
86 roots.append(r)
87 except ValueError:
Cell In[10], line 51, in love_dispersion_residual(c, omega, h, vs, rho)
48 alpha = -alpha
50 F = tau + mu_h * alpha * u
---> 51 return float(np.real(F))
File ~/GitHub/ess-412-512-intro2seismology/.pixi/envs/default/lib/python3.11/site-packages/numpy/lib/_type_check_impl.py:80, in _real_dispatcher(val)
76 return 'D'
77 return min(intersection, key=_typecodes_by_elsize.index)
---> 80 def _real_dispatcher(val):
81 return (val,)
84 @array_function_dispatch(_real_dispatcher)
85 def real(val):
KeyboardInterrupt:
def plot_vs_profile(h, vs, label=None):
h = np.asarray(h)
vs = np.asarray(vs)
z_edges = np.concatenate([[0], np.cumsum(h)])
z_max = z_edges[-1] + 2000
z = [0.0]
v = [vs[0]]
z0 = 0.0
for i, hi in enumerate(h):
z1 = z0 + hi
z += [z1, z1]
v += [vs[i], vs[i+1]]
z0 = z1
z += [z_max]
v += [vs[-1]]
plt.plot(v, z, label=label)
plt.figure()
plot_vs_profile(h_true, vs_true, label="true")
plot_vs_profile(h_est, vs_est, label="estimated")
plt.gca().invert_yaxis()
plt.xlabel("Vs (m/s)")
plt.ylabel("Depth (m)")
plt.legend()
plt.grid(True)
plt.title("1-D Vs model (3 layers + halfspace)")
plt.show()
Discussion prompts#
Which frequencies constrain the shallowest layer most strongly? Why?
Increase dx: when do you see spatial aliasing in F–K? How does it corrupt picked c(f)?
Hold
hfixed and invert onlyVs: does the fit improve or worsen? What does that say about trade-offs?Increase regularization
lam: what changes in Vs(z) and in the misfit?
Extensions#
Add a second mode in the synthetic wavefield (two ridges) and discuss mode mixing.
Replace ridge-picking by a weighted centroid in k (soft picks with uncertainties).
Compare phase-velocity inversion to group-velocity inversion (time–frequency methods).