Skip to content

WP5: walkaway VSP anisotropy and the gamma -> (eps, delta) calibration - #6

Merged
ingehap merged 1 commit into
mainfrom
claude/blissful-gauss-3y1u5w
Oct 14, 2026
Merged

ingehap merged 1 commit into
mainfrom
claude/blissful-gauss-3y1u5w

Conversation

@ingehap

@ingehap ingehap commented Oct 14, 2026

Copy link
Copy Markdown
Owner

Implements work package 5 of docs/IMPLEMENTATION_PLAN.md: the walkaway-VSP ε and δ of Lewis et al. (2016), following Leaney & Hornby (2007), and the calibration of the sonic γ with them.

What is added

  • iwt.raytrace (exact qP kinematics in VTI, fluids included; Tsvankin, 2012):
    • phase velocity, vertical slowness q(p) and dq/dp, and the polarization;
    • two-point direct rays in layered VTI.
    • Rays steeper than 85° in any layer they cross are dropped. Past that point the first arrival would be a refraction, which is not modelled.
  • iwt.vsp.walkaway:
    • Synthetic survey: several azimuths, source statics and noise.
    • Measurements: p, q and polarization, with their joint covariance propagated from the picks. Overlapping source windows make neighbouring p values correlated.
    • Inversion: (Vp0, Vs0, ε, δ) fitted on V(θ) and ψ(θ) with the exact Christoffel equation, or with Thomsen's weak form for comparison.
    • Error model: Gaussian priors and an effective-variance covariance, frozen within each iteration. The Laplace posterior reports a data_weight per parameter.
    • Checks: a χ² test across azimuths, and a polarization check for lateral heterogeneity.
  • iwt.calibration:
    • upscale_gamma: Backus upscaling, with C66 averaged arithmetically and C44 harmonically.
    • bayes_linear: straight-line fit with an extra scatter term (type-II maximum likelihood), leave-one-out predictions and predictive bands.
  • examples/wp5_walkaway.py and docs/figures/wp5_*.png cover:
    • six levels at 30–60°, inverted five ways;
    • the near-angle trade-off and the consistency checks;
    • heterogeneity inside the sub-array;
    • the calibrated ε and δ curves.
  • Tests: 17 new, including a Monte Carlo check of the covariance, a coverage test, the near-angle trade-off, the weak-anisotropy bias, the Vs0 prior, both checks, and missing picks.
  • README and plan updated; version 0.5.0.

Main findings (synthetic well)

  • Fast stringers limit the direct arrivals to 775–1450 m of offset, which is 23–78° at the levels. This may be why the article had only about 30–60°.
  • Adding polarizations halves the ε error and triples the δ precision.
    • Their cost is an errors-in-variables bias of about 0.5σ, which scales with the noise variance.
    • With polarizations, ε and δ also depend on the Vs0 prior (Grechka & Mateeva, 2007).
    • Fitting slowness only gives calibrated errors, but Vp0 then comes from its prior alone.
  • Heterogeneity inside the 40 m sub-array shifts ε by up to 0.03 while χ² stays small, so the posterior σ is too small. Add a model error of about 0.01.
  • The polarization check found a 5 ms/km source static at only four of six levels. That static biases the slowness-only δ by up to 0.11.
  • The γ → ε calibration fits well inside 640–860 m but fails below. Below 860 m half of the true ε curve falls outside its 2σ band, because the WP4 γ bias changes with depth there.

Validation

  • Python 3.13 with RuntimeWarning and DeprecationWarning as errors: 128 passed.
  • Python 3.9 with numpy 1.21 and scipy 1.7: 124 passed, 1 skipped (plotting).
  • pyflakes is clean, and the README quick start runs.

🤖 Generated with Claude Code

https://claude.ai/code/session_019vaUDUPicoFeQaUhfhnoYq


Generated by Claude Code

- iwt.raytrace: exact qP phase velocity, vertical slowness q(p) and dq/dp,
  polarization (Christoffel eigenvector) in VTI incl. fluids; two-point
  direct rays in layered VTI, excluding rays steeper than max_angle (85 deg)
  in any layer (refraction regime).
- iwt.vsp.walkaway: synthetic walkaway (azimuths, source statics, noise),
  slowness/polarization measurement with the joint covariance propagated
  from the picks, (Vp0, Vs0, eps, delta) inversion on V(theta) and
  psi(theta) with the exact Christoffel equation (or Thomsen's weak form),
  Gaussian priors, frozen effective-variance covariance and Laplace
  posterior with data weights; azimuthal chi2 consistency test and a
  polarization check for lateral heterogeneity.
- iwt.calibration: Backus upscaling of the sonic gamma, Bayesian straight
  line with extra scatter (type-II ML), leave-one-out predictions.
- examples/wp5_walkaway.py and docs/figures/wp5_*.png: six levels at
  30-60 degrees, five inversion variants, near-angle trade-off, checks,
  sub-array heterogeneity, and the calibrated eps/delta curves.
- Tests (17 new), README and plan updated; version 0.5.0.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_019vaUDUPicoFeQaUhfhnoYq
@ingehap
ingehap marked this pull request as ready for review October 14, 2026 17:03
Copilot AI balanced review requested due to automatic review settings October 14, 2026 17:03
@ingehap
ingehap merged commit 19171c6 into main Oct 14, 2026
4 checks passed

Copilot AI left a comment

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🟡 Changes recommended

Calibration leakage, uncertainty diagnostics, and kinematic boundary errors remain unresolved.

6 open findings
What changed in this PR

Implements WP5, connecting walkaway VSP anisotropy estimates with the existing sonic calibration workflow.

Changes:

  • Adds exact VTI ray tracing, walkaway measurements, inversion, and consistency checks.
  • Adds gamma upscaling and Bayesian calibration to epsilon and delta.
  • Adds tests, a synthetic example, documentation, and version 0.5.0.
File Description
tests/​test_walkaway.py Tests ray tracing and walkaway inversion.
tests/​test_calibration.py Tests upscaling and linear calibration.
README.md Documents WP5 usage, results, and limitations.
pyproject.toml Updates package version.
iwt/​vsp/​walkaway.py Implements survey processing, inversion, and checks.
iwt/​vsp/​__init__.py Documents the walkaway module.
iwt/​raytrace.py Adds exact qP kinematics and layered rays.
iwt/​calibration.py Adds gamma upscaling and Bayesian regression.
iwt/​__init__.py Exposes ray tracing and updates version.
examples/​wp5_walkaway.py Demonstrates inversion and calibration.
docs/​IMPLEMENTATION_PLAN.md Records WP5 completion and deviations.

🧠 Review effort: Balanced


💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.

Comment thread iwt/calibration.py
Comment on lines +241 to +243
n = max(1, int(round(length / dz)) | 1)
c44m = 1.0 / _running_mean(1.0 / np.asarray(c44, float), n)
c66m = _running_mean(np.asarray(c66, float), n)
Comment thread iwt/calibration.py
Comment on lines +301 to +322
def log_evidence(tau):
var = s ** 2 + tau ** 2
G = np.column_stack([np.ones_like(x), x])
C = np.diag(var) + G @ np.linalg.inv(P0) @ G.T
r = y - G @ m0
sign, logdet = np.linalg.slogdet(C)
return -0.5 * (logdet + r @ np.linalg.solve(C, r))

if scatter is None:
spread = max(float(np.std(y)), float(np.max(s)), 1e-6)
grid = np.concatenate([[0.0], np.geomspace(1e-4 * spread, 3.0 * spread, 200)])
scatter = float(grid[np.argmax([log_evidence(t) for t in grid])])
var = s ** 2 + scatter ** 2
mean, cov = _line_posterior(x, y, var, m0, P0)
loo_m = np.empty(x.size)
loo_s = np.empty(x.size)
for i in range(x.size):
keep = np.arange(x.size) != i
mi, ci = _line_posterior(x[keep], y[keep], var[keep], m0, P0)
g = np.array([1.0, x[i]])
loo_m[i] = g @ mi
loo_s[i] = np.sqrt(g @ ci @ g + scatter ** 2)
Comment thread iwt/raytrace.py
Q_s = np.where(b < 0, 2.0 * c / (-b + disc), (-b - disc) / (2.0 * a))
Q_f = 1.0 / a33 - P
Q = np.where(fluid, Q_f, Q_s)
q = np.where(Q > 0, np.sqrt(np.where(Q > 0, Q, 1.0)), np.nan)
Comment thread iwt/raytrace.py
Comment on lines +104 to +105
ux = (a13 + a44) * p * q
uz = 1.0 - a11 * p * p - a44 * q * q
Comment thread iwt/vsp/walkaway.py
Comment on lines +307 to +308
return np.where(np.isfinite(self.prior_std),
1.0 - np.diag(self.cov) / self.prior_std ** 2, 1.0)
Comment thread iwt/vsp/walkaway.py
e[j] = 1e-6 * max(abs(x[j]), 1.0)
G[:, j] = (_forward(x + e, data.theta)[1] - _forward(x - e, data.theta)[1]) / (2 * e[j])
cov = np.where(np.outer(ok, ok), posterior.cov, 0.0)
R = R + G @ cov @ G.T
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants