diff --git a/README.md b/README.md index 9067fa9..317f7ba 100644 --- a/README.md +++ b/README.md @@ -18,8 +18,9 @@ The roadmap, assumptions, critical assessment and full reference list are in **Status:** work packages 1 (foundations), 2 (forward modelling), 3 (ZO-VSP processing, Q and the sonic calibration), 4 (sonic anisotropy), 5 (walkaway anisotropy and the γ → (ε, δ) -calibration) and 6 (ray tracing and Q-compensated depth imaging, the Fig. 3 experiment) are -implemented. WP7–WP8 (the integrated loop, real-data I/O) are still to do. +calibration), 6 (ray tracing and Q-compensated depth imaging, the Fig. 3 experiment) and 7 (the +integrated loop with its discrepancy register, multiples, the Fig. 2, 4 and 5 analogues) are +implemented. WP8 (real-data I/O and a command line) is still to do. ![Synthetic well](docs/figures/wp1_well.png) @@ -135,6 +136,25 @@ gather = np.zeros((x.size, 1000)) # data[offset, image = imaging.migrate_cmp(gather, 0.001, x, mig, tab.depth, 6.0, 8.0, table=tab, q_compensation="amplitude", gain_limit=10.0, high_cut=80.0) fit = avo.intercept_gradient(tab.theta[0], np.ones(x.size), (20.0, 44.0)) # A + B sin^2(theta) + +# the integrated loop: model versions, discrepancy register, tie metrics, multiples (WP7) +from iwt import multiples, welltie, workflow +from iwt.synthetic import SEA_WATER + +reg = workflow.Register() +v0 = reg.new_version("logs", workflow.model_from_logs(logs, 417.5, SEA_WATER.vp0, SEA_WATER.rho), + "sonic, dipole and density as logged") +v1 = reg.new_version("ZO-VSP calibrated", # sonic moved to 30 Hz, gap from VSP + workflow.model_from_logs(logs, 417.5, SEA_WATER.vp0, SEA_WATER.rho, qp=60.0, + qs=30.0, gap_vp=2130.0), "Q, dispersion, gap") +reg.log("synthetic late against the corridor stack", "blind gap above the logs", + "fill the gap from the VSP", "shift 1 ms", "adopt", version=v1.number) +reg.check("corridor stack vs synthetic, correlation", 0.93, 0.8, ">=") +print(reg.to_markdown()) # versions, discrepancies, criteria +syn = normal_incidence(v1.model, 0.002, 600, wav)[1] +m = welltie.tie_metrics(np.roll(syn, 2), syn, 0.002) # correlation, shift, phase, scale, PEP +g = multiples.gathers(v1.model.subset(0.0, 780.0), [300.0], 6.0, 8.0, 0.002, 600, wav, f_max=90.0) +g.keys() # full, primaries, no mode conversion, ..., Born primaries ``` ## What WP1 provides @@ -223,6 +243,19 @@ about 1400 slownesses) with `workers=4`. | AVA | `iwt.avo.pick_event`, `intercept_gradient` | event tracking across a migrated gather (follows residual moveout); A + B sin²θ fit | | figures | `examples/wp6_imaging.py` | `docs/figures/wp6_*.png` | +## What WP7 provides + +| Plan item | Code | Notes | +|---|---|---| +| model from logs | `iwt.workflow.model_from_logs` | 0.5 m elastic model: sonic and dipole moved to the seismic frequency with Q (Kjartansson, 1979), ε, δ, γ and Q curves, the sea, a gap layer above the logs (velocity from the VSP if given) | +| loop and register | `iwt.workflow.Register` | numbered model versions with their changes; discrepancy entries (observation, hypothesis, test, result, decision, version); acceptance criteria fixed before tuning; Markdown output | +| well-tie metrics (Fig. 2) | `iwt.welltie.tie_metrics`, `estimate_wavelet`, `rotate_phase` | zero-lag and best correlation; time shift and constant phase from the complex cross-correlation of analytic signals; scale; PEP. Least-squares multi-trace wavelet (Walden & White, 1998) | +| multiples and mode conversion | `iwt.multiples.gathers`, `PHYSICS`, `relative_rms` | the same CMP gather with physics switched off step by step (full; primaries with their mode conversions; no mode conversion; P-wave primaries; Born primaries) | +| identifiability | `iwt.avo.linear_resolution` | Gaussian posterior of linearized parameters: posterior std, data weight, posterior correlation | +| Fig. 5 section | `iwt.synthetic.dipping_reservoir_models` | columns with a wedge-shaped anisotropic shale, a dipping reservoir and a flat oil–water contact; an all-brine twin for the fluid signal | +| exact PP in slowness | `iwt.modeling.reflection.exact_pp_slowness` | the VTI coefficient as a function of the horizontal slowness (complex beyond critical) | +| the loop | `examples/reproduce_paper.py` | `docs/discrepancy_register.md`, `docs/figures/wp7_*.png` | + ## The synthetic well The Hanssen data are proprietary, and the article gives no velocities, densities, anisotropy or @@ -773,6 +806,176 @@ V_nmo ≈ 1820 m/s, with errors of at most 0.71 ms out to offset/depth 2. ![AVA after migration](docs/figures/wp6_ava.png) +## WP7 results (synthetic well; `python examples/reproduce_paper.py`, about 6 min with 4 workers) + +**The loop.** One script runs the article's workflow (its Fig. 1) using estimated models only: +logs, then the ZO-VSP, then the walkaway, then the surface seismic. The true earth is used for +scoring. In the entries marked *synthetic only*, it also attributes misfits that field data +could not attribute. Every model version, discrepancy and acceptance criterion is written to +[`docs/discrepancy_register.md`](docs/discrepancy_register.md). + +| version | from | change | +|---|---|---| +| v0 | logs | sonic, dipole and density as logged (12 kHz, 2 kHz); isotropic; no Q | +| v1 | + ZO-VSP | Vp moved to 30 Hz with the VSP intrinsic Q; Qs = Qp/2; the blind gap above the logs (417.5–470 m) filled from the VSP | +| v2 | + walkaway | ε and δ from the sonic γ through the walkaway calibration | +| (v3) | + "rock-physics" Vs | dipole Vs moved to 30 Hz with Qs = Qp/2: tested and rejected | + +**1. ZO-VSP (Fig. 2 analogue, left).** + +- **Intrinsic Q per block:** 46, 60, 123 and 72 (true 46, 59, 80 and 70). The thin Fruholmen + block (60 m) is the poorly constrained one. +- **Drift (VSP minus integrated sonic):** 11.9 ms at 1650 m, 7.5 ms RMS. After the physical + calibration the residual is 0.22 ms RMS. +- **Corridor stack against the v1 synthetic:** with the gap above the logs filled from the first + 10 m of logs, the correlation is −0.06 and the shift 8 ms. With the VSP interval velocity in + the gap (2129 m/s; true 2067 m/s), the correlation is 0.93 (0.97 after the best shift and + rotation), the shift +1.3 ms and the phase +4°. +- **The 52 m blind interval between the seabed and the first log controls the tie.** In this + synthetic, a model built from logs alone does not tie a reservoir 250 m below the seabed. + +**2. Anisotropy.** + +- **Walkaway ε** (true 40 m value) at 650, 690, 720, 770, 810 and 850 m: 0.145 (0.145), + 0.060 (0.038), 0.082 (0.084), 0.049 (0.053), 0.103 (0.093) and 0.101 (0.102). δ is off by up + to 0.03. +- **Calibration:** ε = 0.039 + 0.66 γ and δ = −0.007 + 0.49 γ. +- **Unit means of the calibrated curves** (true values in brackets): + - δ is 0.035–0.04 too high in the shallow shales: Lower Cretaceous 0.117 (0.077) and + Hekkingen 0.135 (0.100); + - Hekkingen ε is too low: 0.230 (0.260). + + Carrying the six-level calibration up into the overburden is an extrapolation, and this is + where it fails (plan, critical point 3). + +**3. Surface seismic (Fig. 2 right and Fig. 4 analogues).** + +- **The measurement** is the full-physics reflectivity response of the true earth, including + internal multiples, mode conversion and Q. Offsets run from 100 to 1000 m. There is no free + surface: the data are assumed deghosted and free-surface demultipled. +- **Processing:** the measurement is migrated with v2 (20 m Backus, VTI, VSP Q compensation, + gain limit 10, 80 Hz high cut). It is compared with v2's elastic synthetic, processed the + same way without Q compensation. + +| criterion (fixed before tuning) | value | threshold | passed | +|---|---|---|---| +| VSP first-break residual after calibration | 0.22 ms | ≤ 1 ms | yes | +| corridor stack vs synthetic, correlation | 0.93 | ≥ 0.8 | yes | +| reservoir event time shift (near-offset tie) | 0.02 ms | ≤ 2 ms | yes | +| Stø-top AVA gradient, contamination-free, \|predicted/measured − 1\| | 0.40 | ≤ 0.10 | **no** | +| Stø-top AVA gradient, as acquired, \|predicted/measured − 1\| | 0.73 | ≤ 0.10 | **no** | + +- **The tie closes; the AVA does not.** + - **Near offsets (≤ 300 m):** correlation 0.87 (0.98 after the best rotation), no depth + shift, scale 0.96 and a residual phase rotation of +27°, which was not investigated + further. + - **As acquired**, the reservoir AVA is erratic beyond 25° (Fig. 4 analogue, top row). The + Stø misfit is 0.62 and the predicted/measured gradient 0.27; at Fruholmen they are 1.28 and + 3.59. This reproduces the article's observation. + +![Well ties](docs/figures/wp7_tie.png) + +**4. Discrepancies.** + +- **Contamination: elastic, and growing with angle.** + - After identical migration, full − Born is 0.93 × Born in the Stø window and 1.53 × + at Fruholmen. At offsets up to 300 m it is only 0.26 and 0.40. + - With mode conversion switched off, the P-wave interbed multiples are only 0.14 and + 0.23 × Born. The large far-offset contamination appears only with mode conversion on, + so an acoustic (P-wave) demultiple would not repair this AVA. + - The switches cannot split it further. The primaries with their mode conversions differ + from Born by 1.35 × and 2.34 ×, and the multiples cancel most of that again. This + is the elastic counterpart of the O'Doherty & Anstey (1971) cancellation. + - "Converted waves" and "interbed multiples" are therefore not separate culprits here. The + article's interbed multiples are, in this synthetic, elastic multiples with converted legs. + - The remaining hypotheses are tested on contamination-free (Born) data, which is what an + ideal demultiple would deliver. +- **"Rock-physics" Vs: rejected.** Moving the dipole Vs to 30 Hz (Qs = Qp/2) raises the + contamination-free Stø misfit from 0.141 to 0.164 and lowers Fruholmen's from 0.062 to + 0.051, so it gives no overall improvement. + - The identifiability study below fits the Stø residual best with Vs +3 %, the opposite + sign to the physical dispersion correction. + - Fitting Vs to the AVA would therefore have produced a wrong "rock-physics Vs" that + absorbs errors made elsewhere. +- **Model or processing? (synthetic only)** + - The true earth's own elastic prediction, processed like v2's, lies 2–14 % above the + measurement over 20–44° (mean misfit 0.088). + - Re-migrating the measurement with the true Q above the logs lowers the misfit to 0.056; + with the true Q everywhere it falls to 0.040. That is the floor set by the first-order Q + operator and the shear Q. + - The cause: the VSP has no receivers between the seabed and 450 m, so the Quaternary (true + Q 30) takes the first block's Q of 46. This leaves about 3 % of amplitude error at the + reservoir that no reservoir parameter should absorb. +- **Can the model error be split one part at a time? (synthetic only) No.** + - Swapping parts of v2 into the true earth gives a Stø misfit of 0.40 with v2's ε and δ, and + 0.23 with v2's Vp, Vs and ρ. With both parts (that is, v2 itself) it is 0.14. + - The cause: the 5 m-smoothed γ gives the isotropic calcite stringers the shale's + anisotropy, ε 0.13–0.15 and δ 0.06–0.07 instead of 0. + - On the sharp 4300 m/s stringers of the true earth, this lowers their critical angle from + about 40° to about 35°, and the predicted Stø amplitude collapses beyond about 33°. + - In v2 the stringers are also smoothed, by the log apertures (3800 m/s instead of 4300 m/s + for the 1 m stringer), so the two errors partly cancel. +- **Identifiability: the reservoir AVA alone cannot choose between the article's + explanations.** + - The contamination-free Stø residual (measured/predicted − 1) is +13 % to +20 % over + 20–38° and −12 % at 44°. + - Each candidate alone explains only part of it, leaving 56–94 % of its RMS: Vs below + +3.1 %, Vs above −4.5 %, ρ below −3.1 %, δ below −0.078, ε below −0.073, or a 20 % change + of the overburden Q. + - Fitted jointly, with 2 % noise per angle, the data constrain ρ below (data weight 0.89) + and the overburden Q (0.60). They barely constrain the Vs on either side (0.38 and 0.17) + and not δ (0.05). + - Vs below and the overburden Q correlate +0.59 in the posterior. The overburden Q, which + the true-earth test identifies as a real culprit here, trades off against the reservoir Vs. + +![AVA](docs/figures/wp7_ava.png) + +**5. Fig. 5 analogue.** + +- **The section:** zero-offset columns every 20 m over 1 km. The Hekkingen shale thickens at + 1.2°, so the reservoir dips under a flat oil–water contact at 694.5 m and the oil column + thins from 22.5 to 1.5 m. +- **The fluid signal** is the image of an oil section minus that of an all-brine twin (5 % + noise). Both are migrated with the 20 m Backus averages of the true + oil-section models. + +| imaging | fluid-signal depth (m) | half-max thickness (m) | peak / noise RMS | +|---|---|---|---| +| VTI + Q compensation | 683.9 | 43.2 | 13.6 | +| VTI, no Q compensation | 681.3 | 44.6 | 16.7 | +| isotropic (NMO velocity) + Q compensation | 700.8 | 45.2 | 13.8 | +| isotropic (NMO velocity), no Q compensation | 698.0 | 46.5 | 16.7 | + +- **Anisotropy positions the fluid signal.** VTI places it inside the oil column (top Stø at + 672–693 m, contact at 694.5 m). The isotropic NMO-velocity image puts it 17 m deeper, below + the contact. +- **Q compensation hardly sharpens it here.** It thins the signal by only 3 %, and it costs + signal-to-noise (16.7 → 13.6) in 1D at zero offset. +- **What this 1D test lacks.** The article's clearer contact (its Fig. 5) needs 2D/3D + migration, where restored high frequencies help focus a flat event against dipping ones. +- **At 30 Hz the 2–22 m oil column is below tuning thickness everywhere,** so the contact + cannot be seen in the oil section itself (second panel below). + +![Fig. 5 analogue](docs/figures/wp7_section.png) + +**What WP7 says about the article.** + +1. **The time tie closes with the ZO-VSP alone,** once the blind interval above the logs is + filled from it. Anisotropy matters for depth positioning, here by 17 m. +2. **The AVA cannot close without removing elastic contamination.** P-wave interbed multiples + alone are small here; with mode conversion the contamination reaches 1–2 × Born beyond + 300 m. Elastic modelling and subtraction, or elastic inversion of the reservoir + interval, would be needed; an acoustic demultiple would not suffice. +3. **A reservoir Vs change can stand in for an overburden Q error** (posterior correlation + +0.59). A "rock-physics Vs" is justified only by independent evidence, such as + dipole-dispersion analysis or a Vs from the walkaway shear arrivals, not by the AVA fit. +4. **The parts of the model must share one resolution.** A better calibration would fit the + intrinsic γ → (ε, δ) relation so that the Backus average of the log model at the walkaway + scale (40 m) reproduces the walkaway ε and δ. Anisotropy would then be zeroed in log-flagged + stringers. The direct fit used here is the article's method. +5. **The shallowest Q** (seabed to the first VSP receiver) needs its own constraint: the + source signature or the seabed multiples. + ## Assumptions 1. **"LFP article" means Lewis et al. (2016)**, the only article in [`References/`](References/). @@ -856,14 +1059,36 @@ V_nmo ≈ 1820 m/s, with errors of at most 0.71 ms out to offset/depth 2. the whole ray. The reflectivity modelling also lets the shear Q act on oblique P waves. The difference is at most 2–3 % in amplitude at 35°. 25. **The "observed" data are noise-free Born primaries without free surface** (WP6): no ghosts, - no multiples, no converted waves. Interbed multiples and the free surface are WP7's subject. - The migration models are 20 m Backus averages of the true model; WP7 will use the - estimated ones (ZO-VSP Q, walkaway ε and δ). + no multiples, no converted waves. The migration models are 20 m Backus averages of the true + model. WP7 adds internal multiples and mode conversion and uses the estimated models (ZO-VSP + Q, walkaway ε and δ). 26. **AVA conventions** (WP6): amplitudes are trough picks on the migrated gather, followed across offsets. The angle is that of the migration model at the picked depth. The fits use 20–44°: the 3D data in the article start at 20°, and beyond 44° the seabed reflection approaches the reservoir. +27. **The WP7 surface "measurement" has no free surface and no noise** (deghosted and + free-surface demultipled). Noise is added only to the ZO-VSP, the walkaway and the Fig. 5 + section. In 417.5 m of water the first seabed multiple arrives near 1.1 s, after the + reservoir, but peg-legs from the shallow section would add to the contamination. +28. **Every estimated model uses Qs = Qp/2,** because the VSP gives only the P-wave Q. The model + takes the Q of the first VSP block between the seabed and the first receiver (450 m). +29. **"Contamination-free" means the Born primaries:** an ideal demultiple that also removes + the transmission and conversion losses. A real demultiple keeps those losses, so its + output would also need amplitude balancing. +30. **The identifiability Jacobian is linearized** around 8 m unit means on either side of the + Stø top, with the exact VTI interface coefficient. It ignores the tuning by the stringers. + The overburden-Q column is −π f t\*(θ) per relative change of Q, with f = 35 Hz and t\* + from the migration table. The prior standard deviations are assumed: Vs 3 %, ρ 2 %, ε and + δ 0.03, and Q 30 %. +31. **The AVA acceptance criterion uses 20–44°, not the plan's 30–50°,** because the seabed + reflection crosses the reservoir at about 54° (WP6). The criteria were fixed before any + tuning and were not changed after the results. +32. **The Fig. 5 section is imaged column by column at zero offset** with the 20 m Backus + averages of the true oil-section models. The isotropic variant uses the NMO velocity Vp0 √(1 + 2δ). The + contact is scored by the fluid signal (oil image minus brine image, without the noise + image), not by picking: a 2–22 m column at 30 Hz gives no separate contact event to pick. + ## Notes found while implementing - **The plan defined interbed multiples as "full minus primaries-only" (WP7). That is misleading @@ -927,6 +1152,29 @@ V_nmo ≈ 1820 m/s, with errors of at most 0.71 ms out to offset/depth 2. suite found counter-examples, in about 10 % of random stacks. The correct statements are γ ≥ 0 and ε ≥ δ (η ≥ 0); ε can be negative when Vp contrasts are small and Vs contrasts large (Berryman et al., 1999). The plan and the tests now say so. +- **Ray-Born synthetics fail for thin fast stringers beyond about 27°** (WP7). A ray-Born + gather in a smoothed background loses the cancellation between the tops and bases of thin + beds once they go post-critical, and the reservoir AVA was wrong by tens of per cent. The + predictions therefore use the full reflectivity method. +- **A real cross-correlation mistakes a phase rotation for a time shift.** For a band-limited + wavelet a constant rotation φ looks like a shift of about φ/(2πf), about 3 ms per 30° at + 30 Hz. The tie metrics use the + complex cross-correlation of the analytic signals: its envelope gives the shift and its + phase the rotation. +- **The blind interval above the first log controls a shallow tie.** Filling the 52 m gap with + log velocities made the synthetic 8 ms late (correlation −0.06). The VSP interval velocity + fixed it (0.93). +- **The contamination cannot be split into "multiples" and "converted waves" by switching + physics off** (WP7). Without mode conversion, P-wave multiples are small (0.11–0.34 × Born); + with it, the converted primaries and the multiples are each larger than their sum and + largely cancel. Switching mode conversion off is not a physical medium either: at far + offsets the P-only response loses the converted energy and gets none of it back. The + contamination is reported as full − Born, near and far. +- **Swapping the truth in one part at a time does not attribute the misfit when the parts have + different resolutions** (synthetic only). The smoothed anisotropy curves give the stringers + shale anisotropy; that error is harmless on smoothed velocities and severe on sharp ones. +- **The physical Vs correction and the AVA-fitted Vs have opposite signs** on this well. A Vs + tuned to the AVA would have absorbed the overburden Q error. ## Repository layout @@ -941,7 +1189,10 @@ iwt/ sonic.py C33, C44, C66 and gamma; tube wave with tool; exact TI Stoneley mode; uncertainty raytrace.py exact qP phase and group velocity, slowness, polarization; direct and reflected rays imaging.py constant-Q operators, gain-limited inverse Q; Q-compensated CMP depth migration - avo.py event tracking on migrated gathers; intercept-gradient fits + avo.py event tracking on migrated gathers; intercept-gradient fits; linear resolution + welltie.py wavelet estimation, phase rotation, tie metrics (correlation, shift, phase, PEP) + workflow.py model from logs; model versions, discrepancy register, acceptance criteria + multiples.py gathers with physics switched off step by step; relative RMS plotting.py well, anisotropy, AVO, gather, depth-profile and trace displays (matplotlib) modeling/ elastic.py VTI plane waves with constant Q: slownesses, eigenvectors, interface scattering @@ -958,8 +1209,9 @@ examples/wp3_zovsp.py ZO-VSP processing, Q, drift calibration; doc examples/wp4_sonic_anisotropy.py Stoneley dispersion, gamma processing ladder, uncertainty; wp4_*.png examples/wp5_walkaway.py walkaway eps and delta, checks, gamma calibration; wp5_*.png examples/wp6_imaging.py moveout, Q-compensated migration variants (Fig. 3), gain limit; wp6_*.png -tests/ pytest suite (143 tests) -docs/ implementation plan and figures +examples/reproduce_paper.py the integrated loop (WP7); register, Figs. 2, 4, 5 analogues; wp7_*.png +tests/ pytest suite (154 tests) +docs/ implementation plan, discrepancy register and figures References/ the article (PDF) ``` @@ -1008,6 +1260,7 @@ Run the tests with `pytest`. The plotting tests are skipped when matplotlib is m * Thomsen, L. (1986). Weak elastic anisotropy. *Geophysics*, 51(10), 1954–1966. https://doi.org/10.1190/1.1442051 * Tsvankin, I. (2012). *Seismic Signatures and Analysis of Reflection Data in Anisotropic Media* (3rd ed.). SEG. https://doi.org/10.1190/1.9781560803003 * Valero, H.-P., Djikpesse, H. A., & Sinha, B. K. (2009). Estimation of borehole fluid slowness using sonic array waveforms. *SEG Technical Program Expanded Abstracts 2009*, 361–365. https://doi.org/10.1190/1.3255617 +* Walden, A. T., & White, R. E. (1998). Seismic wavelet estimation: A frequency domain solution to a geophysical noisy input-output problem. *IEEE Transactions on Geoscience and Remote Sensing*, 36(1), 287–297. https://doi.org/10.1109/36.655337 * Walsh, J., Sinha, B., Plona, T., & Ammerman, M. (2007). Derivation of anisotropy parameters in a shale using borehole sonic data. *SEG Technical Program Expanded Abstracts 2007*, 323–327. https://doi.org/10.1190/1.2792435 * Wang, Y. (2002). A stable and efficient approach of inverse Q filtering. *Geophysics*, 67(2), 657–663. https://doi.org/10.1190/1.1468627 * White, J. E. (1983). *Underground Sound: Application of Seismic Waves*. Elsevier. @@ -1017,4 +1270,4 @@ Run the tests with `pytest`. The plotting tests are skipped when matplotlib is m * OMV (2014). OMV finds hydrocarbons at Hanssen well in Norway (press release). https://www.omv.com/en/media/press-releases/2014/omv-finds-hydrocarbons-at-hanssen-well-in-norway * Sokkeldirektoratet. Wellbore 7324/7-2 factpage. https://factpages.sodir.no/en/wellbore/pageview/exploration/all/7450 -DOIs were checked against the Scite index; the WP3 and WP4 references were also checked against publisher records (the White & Simm page range could not be confirmed and is omitted). The exact tool-term notation of Norris (1990) could not be checked against its full text; the form used here was derived from the annulus compliance. The WP5 references Grechka & Mateeva (2007) and Grechka et al. (2024) were checked against the Scite index; Gaiser (1990), Miller & Spencer (1994) and Tsvankin (2012) by their DOIs. The WP6 references (Alkhalifah & Tsvankin, 1995; Wang, 2002; Cavalca et al., 2013) are those of the plan, checked there by DOI. The plan lists further references for WP7–WP8. +DOIs were checked against the Scite index; the WP3 and WP4 references were also checked against publisher records (the White & Simm page range could not be confirmed and is omitted). The exact tool-term notation of Norris (1990) could not be checked against its full text; the form used here was derived from the annulus compliance. The WP5 references Grechka & Mateeva (2007) and Grechka et al. (2024) were checked against the Scite index; Gaiser (1990), Miller & Spencer (1994) and Tsvankin (2012) by their DOIs. The WP6 references (Alkhalifah & Tsvankin, 1995; Wang, 2002; Cavalca et al., 2013) and the WP7 references (Walden & White, 1998; White & Simm, 2003; O'Doherty & Anstey, 1971; Kjartansson, 1979) are those of the plan, checked there by DOI. The plan lists further references for WP8. diff --git a/docs/IMPLEMENTATION_PLAN.md b/docs/IMPLEMENTATION_PLAN.md index 66461d7..281e900 100644 --- a/docs/IMPLEMENTATION_PLAN.md +++ b/docs/IMPLEMENTATION_PLAN.md @@ -1,6 +1,6 @@ # Implementation plan: integrated borehole and surface seismic well tie (Lewis et al., 2016) -Status: **WP1–WP6 implemented** (see the README); WP7–WP8 to do. This document says what will be built, in what +Status: **WP1–WP7 implemented** (see the README); WP8 to do. This document says what will be built, in what order, how each part will be checked, and where the method can be improved. It follows the conventions of the earlier repositories in this series (for example `ES_PHYS_P010`): a Python package, a `paper.py` that records where each published number comes from, figures redrawn on @@ -375,7 +375,39 @@ reservoir AVA gradient against the elastic synthetic, scored numerically. accuracy; - Q compensation of a forward-attenuated spike recovers it up to the gain limit. -### WP7: Integration, AVO, multiples and the figures +### WP7: Integration, AVO, multiples and the figures (done) + +*As built:* `iwt.welltie`, `iwt.workflow`, `iwt.multiples`, `iwt.avo.linear_resolution`, +`iwt.synthetic.dipping_reservoir_models` and `examples/reproduce_paper.py`, which runs the whole +loop (logs → ZO-VSP → walkaway → surface seismic) on estimated models only and writes +`docs/discrepancy_register.md`. Changes from the plan: + +- **The predictions are full-waveform reflectivity, not ray-Born synthetics.** A ray-Born + gather in a smoothed background was built and dropped: beyond about 27° it loses the + cancellation between the tops and bases of the thin fast stringers, and the reservoir AVA + went wrong by tens of per cent. +- **Tie metrics use the complex cross-correlation of the analytic signals.** A real + cross-correlation trades a constant phase rotation for a time shift. The envelope peak gives + the shift; its phase gives the rotation (White & Simm, 2003). +- **The "observed" surface data are the full-physics response of the true earth** (internal + multiples, mode conversion, Q), without free surface (assumed deghosted and free-surface + demultipled). The contamination is measured after identical migration against the Born + primaries. It grows with offset and needs mode conversion: without it, the P-wave interbed + multiples are small. With it, the converted primaries and the multiples largely cancel, so + switching physics off does not split the contamination further. +- **The Vs scenarios (Fig. 4) are tested on contamination-free (Born) data**, as an ideal + demultiple would provide. On synthetic data the truth can also be swapped into the model one + part at a time, which attributes the remaining misfit; field data cannot do this, so the + identifiability study (section 5) is run on the same residual. +- **Acceptance criteria:** the AVA criterion uses 20–44°, not 30–50°, because the seabed + reflection crosses the reservoir at about 54° (WP6). +- **The Fig. 5 analogue is a zero-offset section of 1D columns** (a dipping wedge and a flat + contact), imaged column by column. Its fluid signal (oil section minus brine section) is + scored by depth and thickness. 2D migration is not modelled, so Q compensation can only + sharpen, not focus. +- **Not built:** hold-out validation across data families and Monte Carlo ensembles over the + WP4/WP5 posteriors (both still recommended), and free-surface multiples. + **Builds:** diff --git a/docs/discrepancy_register.md b/docs/discrepancy_register.md new file mode 100644 index 0000000..043531d --- /dev/null +++ b/docs/discrepancy_register.md @@ -0,0 +1,34 @@ +# Discrepancy register of `examples/reproduce_paper.py` + +Generated by the script on the synthetic well; do not edit by hand. + +## Model versions + +| # | name | changes | +|---|---|---| +| v0 | logs | sonic, dipole and density as logged (12 kHz / 2 kHz); isotropic; no Q | +| v1 | ZO-VSP calibrated | Vp at 30 Hz with VSP intrinsic Q; Qp from the VSP, Qs = Qp/2; gap from VSP | +| v2 | anisotropic | eps, delta from the calibrated sonic gamma | + +## Discrepancies + +| # | observation | hypothesis | test | result | decision | version | +|---|---|---|---|---|---|---| +| 1 | Integrated sonic time short of the VSP time by 12 ms at 1650 m | sonic dispersion (12 kHz vs 30 Hz) and layering | move the sonic to 30 Hz with the VSP intrinsic Q, Backus over 40 m | residual drift 0.22 ms RMS | adopt | v1 | +| 2 | Synthetic 8 ms late against the corridor stack (zero-lag correlation -0.06) | the 418-470 m gap above the logs (Quaternary) is filled with log velocities | fill the gap with the VSP interval velocity | correlation 0.93, shift +1.3 ms | adopt | v1 | +| 3 | Isotropic v1 cannot predict non-hyperbolic moveout or walkaway slownesses | VTI overburden (shales) | slowness-polarization inversion at six levels, gamma -> (eps, delta) calibration | eps slope 0.66, delta slope 0.49 per gamma | adopt eps, delta, gamma curves | v2 | +| 4 | Reservoir AVA beyond 20 deg erratic and not predicted (misfit 0.62 at the Stø top) | interbed multiples and mode conversion from the overburden stringers | compare physics variants of the true-earth gathers after identical processing | full - Born = 0.9 x Born (Stø), 1.5 x (Fruholmen), growing with offset (0.26 at <= 300 m); without mode conversion the P-wave multiples are only 0.14, 0.23 x Born; converted primaries and multiples each exceed the total and largely cancel; contamination-free misfit 0.14 | elastic (converted-wave) contamination: an acoustic demultiple would not repair the AVA; test the remaining hypotheses on contamination-free data | – | +| 5 | Contamination-free reservoir AVA still off: misfit 0.14 (Stø), 0.06 (Fruholmen) | sonic and seismic see Vs differently: dipole Vs is dispersed (2 kHz vs 30 Hz) | Vs moved to 30 Hz with Qs = Qp/2 (the article's rock-physics Vs analogue) | misfit 0.16 (Stø), 0.05 (Fruholmen) | reject: no overall improvement | – | +| 6 | (synthetic only) The true earth itself predicts +10% amplitude at all angles | Q compensation too weak: no VSP receivers between the seabed and 450 m, so the Quaternary (true Q 30) takes the first block's Q (46) | re-migrate the measurement with the true Q above the logs, then everywhere | mean misfit 0.088 -> 0.056 -> 0.040 | explained; constrain the shallow t* (VSP source signature, seabed multiples) on field data | – | +| 7 | (synthetic only) v2's Stø AVA falls with angle against the true earth's | the calibrated eps and delta, or the log resolution of the stringers | swap one part of v2 into the true earth | Stø misfit 0.40 with v2's eps and delta, 0.23 with v2's Vp, Vs and rho; 0.14 with both | not separable: the 5 m-smoothed gamma gives the isotropic stringers shale anisotropy, which v2's equally smoothed velocities partly cancel; keep one resolution | – | +| 8 | Flat oil-water contact hard to see under a dipping reservoir (Fig. 5) | Q-PSDM sharpens it; anisotropy positions it | zero-offset images of an oil and a brine section, VTI vs isotropic, with and without Q | VTI places the fluid signal at 684 m, isotropic at 701 m; Q compensation thins it from 44.6 to 43.2 m | anisotropy is what matters in 1D; the image improvement needs 2D/3D imaging | – | + +## Acceptance criteria (fixed before tuning) + +| criterion | value | threshold | passed | +|---|---|---|---| +| VSP first-break residual after calibration | 0.221 ms | <= 1 ms | yes | +| corridor stack vs synthetic, correlation | 0.934 | >= 0.8 | yes | +| reservoir event time shift (near-offset tie) | 0.0181 ms | <= 2 ms | yes | +| Stø-top AVA gradient, contamination-free, \|predicted/measured - 1\| | 0.396 | <= 0.1 | **no** | +| Stø-top AVA gradient, as acquired, \|predicted/measured - 1\| | 0.727 | <= 0.1 | **no** | diff --git a/docs/figures/wp7_ava.png b/docs/figures/wp7_ava.png new file mode 100644 index 0000000..a5dd76e Binary files /dev/null and b/docs/figures/wp7_ava.png differ diff --git a/docs/figures/wp7_section.png b/docs/figures/wp7_section.png new file mode 100644 index 0000000..f6247aa Binary files /dev/null and b/docs/figures/wp7_section.png differ diff --git a/docs/figures/wp7_tie.png b/docs/figures/wp7_tie.png new file mode 100644 index 0000000..1d9da7d Binary files /dev/null and b/docs/figures/wp7_tie.png differ diff --git a/examples/reproduce_paper.py b/examples/reproduce_paper.py new file mode 100644 index 0000000..d75644c --- /dev/null +++ b/examples/reproduce_paper.py @@ -0,0 +1,685 @@ +"""Work package 7: the integrated loop of Lewis et al. (2016) on the synthetic well. + +The article refines one earth model with the logs, the ZO-VSP, the walkaway VSP and +the surface seismic "until the predictions match the measurements, or until the +remaining discrepancies are explained". This script runs that loop with the +methods of WP1-WP6 and records it in a register (``iwt.workflow``): every model +version, every discrepancy (observation, hypothesis, test, result, decision) and +the acceptance criteria fixed in the plan before any tuning. + +0. v0: the model from the logs as logged (isotropic, no Q). +1. ZO-VSP (WP3): time-depth, interval Q, drift; v1 = sonic moved to 30 Hz with the + VSP intrinsic Q, the blind gap above the logs filled from the VSP. + Fig. 2 analogue: corridor stack against the v1 synthetic. +2. Anisotropy (WP4 + WP5): sonic gamma, walkaway eps and delta at six levels + (priors from this loop's ZO-VSP and dipole), the gamma -> (eps, delta) + calibration; v2 = v1 with eps, delta and gamma curves. +3. Surface seismic (WP6): the "observed" CMP gather is the full-physics + reflectivity response of the true earth (internal multiples, mode conversion, + Q; no free surface: assumed deghosted). It is migrated with v2 (VTI, Q + compensation) and compared with elastic synthetics of v2, as in the article + (Figs. 2 and 4): the near-offset tie, and the AVA of the Stø and Fruholmen + tops. +4. Discrepancies: contamination by multiples and mode conversion (full against + Born gathers), then, on contamination-free (Born) data, the reservoir Vs + (dipole against dispersion-corrected: the article's "rock-physics Vs"), and an + identifiability study. On synthetic data the truth is known: the true earth's + own prediction separates model from processing errors (the measurement is + re-migrated with the true Q), and swapping one part of v2 into the true earth + tests whether the model errors can be attributed one part at a time. +5. Fig. 5 analogue: a dipping reservoir with a flat oil-water contact, imaged at + zero offset with and without anisotropy and Q compensation. + +Writes ``docs/figures/wp7_*.png`` and ``docs/discrepancy_register.md`` and prints +the numbers quoted in the README (about 6 minutes with 4 workers). + + python examples/reproduce_paper.py [--workers 4] +""" + +import argparse +import sys +import time +from pathlib import Path + +import numpy as np + +sys.path.insert(0, str(Path(__file__).resolve().parent)) + +import wp3_zovsp as wp3 # noqa: E402 +import wp5_walkaway as wp5 # noqa: E402 + +from iwt import avo, paper, sonic as S, welltie as T, workflow as WF # noqa: E402 +from iwt import calibration as C # noqa: E402 +from iwt import imaging as I # noqa: E402 +from iwt import multiples as M # noqa: E402 +from iwt.calibration import bayes_linear # noqa: E402 +from iwt.modeling.convolution import normal_incidence, ricker # noqa: E402 +from iwt.modeling.reflection import exact_pp # noqa: E402 +from iwt.qmodel import phase_velocity # noqa: E402 +from iwt.raytrace import reflection_table # noqa: E402 +from iwt.synthetic import (OWC_TVDSS, SEA_WATER, dipping_reservoir_models, # noqa: E402 + wisting_like_well) +from iwt.vsp import attenuation as A # noqa: E402 +from iwt.vsp import walkaway as W # noqa: E402 +from iwt.vsp import zovsp as Z # noqa: E402 + +ROOT = Path(__file__).resolve().parents[1] +FIG = ROOT / "docs" / "figures" +OFFSETS = np.arange(100.0, 1001.0, 50.0) +DT, NT = 0.001, 1150 +Z_S, Z_R = 6.0, 8.0 +X_AVA = 850.0 +FIT = (20.0, 44.0) +EVENTS = {"Stø top": paper.STO_TOP_TVDSS, "Fruholmen top": 740.0} +BORN = dict(internal_multiples=False, converted_waves=False, transmission_loss=False) +OBSERVED = {"full": {}, "primaries": dict(internal_multiples=False), + "P only": dict(converted_waves=False), + "P primaries": dict(internal_multiples=False, converted_waves=False), "Born": BORN} +BOTTOM = 780.0 + + +def unit_mean(m, a, b, name): + s = m.split_at([a, b]) + i = (s.mid > a) & (s.mid < b) + return float(np.sum(getattr(s, name)[i] * s.thickness[i]) / np.sum(s.thickness[i])) + + +# -- 1. ZO-VSP -------------------------------------------------------------------------------- + + +def zovsp_step(model, logs, reg, workers): + z = logs.depth + v0 = reg.new_version("logs", WF.model_from_logs(logs, paper.WATER_DEPTH, SEA_WATER.vp0, + SEA_WATER.rho, disperse_vs=False), + "sonic, dipole and density as logged (12 kHz / 2 kHz); isotropic; no Q") + kw = dict(f_max=wp3.F_MAX, workers=workers) + clean = Z.synthetic_zovsp(model, wp3.DEPTHS, wp3.SOURCE_DEPTH, wp3.DT, wp3.NT, wp3.WAVELET, **kw) + vsp = clean.with_data(Z.add_noise(clean.data, wp3.NOISE, wp3.DT, wp3.F_MAX, seed=3)) + strat = Z.synthetic_zovsp(v0.model, wp3.DEPTHS, wp3.SOURCE_DEPTH, wp3.DT, wp3.NT, wp3.WAVELET, + attenuation=False, **kw) + t = vsp.t + r, rs = wp3.process(t, vsp.data), wp3.process(t, strat.data) + tv = Z.vertical_times(r["phase"], wp3.DEPTHS, wp3.SOURCE_DEPTH, datum_velocity=SEA_WATER.vp0) + q = A.interval_q(wp3.DEPTHS, tv, r["cs"] - rs["cs"], wp3.Q_BLOCKS) + qz = q.at(z, q.q[0]) + d0 = C.drift(z, logs.vp, wp3.DEPTHS, tv) + vcal = C.physical_calibration(z, logs.vp, logs.rho, qz, wp3.F_SONIC, wp3.F_SEIS, + backus_length=wp3.BACKUS) + d1 = C.drift(z, vcal, wp3.DEPTHS, tv) + edges = np.concatenate([[wp3.DEPTHS[0]], wp3.Q_BLOCKS, [wp3.DEPTHS[-1]]]) + q_true = 1.0 / wp3.true_inv_q(model, edges) + print("1. ZO-VSP: intrinsic Q per block (centroid shift, strat. from the log model) | true: " + + " ".join(f"{a:.0f} | {b:.0f}" for a, b in zip(q.q, q_true))) + print(f" drift (VSP - integrated sonic): {1e3 * d0[-1]:.1f} ms at {wp3.DEPTHS[-1]:.0f} m, " + f"RMS {1e3 * np.sqrt(np.nanmean(d0 ** 2)):.1f} ms; after the physical calibration " + f"RMS {1e3 * np.sqrt(np.nanmean(d1 ** 2)):.2f} ms") + reg.log("Integrated sonic time short of the VSP time by " + f"{1e3 * d0[-1]:.0f} ms at {wp3.DEPTHS[-1]:.0f} m", + "sonic dispersion (12 kHz vs 30 Hz) and layering", + "move the sonic to 30 Hz with the VSP intrinsic Q, Backus over 40 m", + f"residual drift {1e3 * np.sqrt(np.nanmean(d1 ** 2)):.2f} ms RMS", "adopt", 1) + # the corridor stack against the synthetic, before and after filling the blind gap + target = ricker(30.0, wp3.DT) + dec = Z.deconvolve(t, r["up"], r["down"], tv, target) + corridor, _ = Z.corridor_stack(t, dec, tv, length=0.1) + w = (t > 2 * tv.min() + 0.02) & (t < 2 * tv.max()) + base = dict(qp=qz, qs=0.5 * qz, disperse_vs=False) + m_a = WF.model_from_logs(logs, paper.WATER_DEPTH, SEA_WATER.vp0, SEA_WATER.rho, **base) + tie_a = T.tie_metrics(corridor[w], normal_incidence(m_a, wp3.DT, wp3.NT, target)[1][w], wp3.DT) + t_sb = paper.WATER_DEPTH / SEA_WATER.vp0 + v_gap = (z[0] - paper.WATER_DEPTH) / (np.interp(z[0], wp3.DEPTHS, tv) - t_sb) + m_b = WF.model_from_logs(logs, paper.WATER_DEPTH, SEA_WATER.vp0, SEA_WATER.rho, gap_vp=v_gap, + **base) + syn = normal_incidence(m_b, wp3.DT, wp3.NT, target)[1] + tie_b = T.tie_metrics(corridor[w], syn[w], wp3.DT) + true_gap = (z[0] - paper.WATER_DEPTH) / float(model.vertical_time(z[0]) - model.vertical_time(paper.WATER_DEPTH)) + print(f" corridor stack vs synthetic: gap from the first logs: correlation {tie_a.correlation:.2f}, " + f"shift {1e3 * tie_a.shift:+.1f} ms; gap from the VSP ({v_gap:.0f} m/s, true {true_gap:.0f}): " + f"correlation {tie_b.correlation:.2f} (best {tie_b.max_correlation:.2f}), shift " + f"{1e3 * tie_b.shift:+.1f} ms, phase {np.degrees(tie_b.phase):+.0f} deg") + reg.log(f"Synthetic {1e3 * tie_a.shift:.0f} ms late against the corridor stack " + f"(zero-lag correlation {tie_a.correlation:.2f})", + "the 418-470 m gap above the logs (Quaternary) is filled with log velocities", + "fill the gap with the VSP interval velocity", + f"correlation {tie_b.correlation:.2f}, shift {1e3 * tie_b.shift:+.1f} ms", "adopt", 1) + v1 = reg.new_version("ZO-VSP calibrated", m_b, + "Vp at 30 Hz with VSP intrinsic Q; Qp from the VSP, Qs = Qp/2; gap from VSP") + reg.check("VSP first-break residual after calibration", 1e3 * np.sqrt(np.nanmean(d1 ** 2)), 1.0, + "<=", "ms") + reg.check("corridor stack vs synthetic, correlation", tie_b.correlation, 0.8, ">=") + return dict(tv=tv, qz=qz, v_gap=v_gap, t=t, corridor=corridor, syn=syn, w=w, v1=v1, tie=tie_b) + + +# -- 2. anisotropy ---------------------------------------------------------------------------- + + +def anisotropy_step(model, logs, reg, zo): + z = logs.depth + g_sonic, g_up = wp5.sonic_gamma(model, logs) + earth40 = model.backus(wp5.SUBARRAY, respect=[model.interfaces[1]]) + survey = W.synthetic_walkaway(earth40, wp5.OFFSETS, wp5.DEPTHS, time_noise=wp5.SIGMA_T, + psi_noise=wp5.SIGMA_PSI, seed=5) + est, rows = [], [] + for lv in wp5.LEVELS: + d = W.measure_slowness(survey, lv, time_noise=wp5.SIGMA_T, psi_noise=wp5.SIGMA_PSI) + th = np.degrees(d.theta) + d = d.select((th >= wp5.ANGLES[0]) & (th <= wp5.ANGLES[1])) + a, b = lv - wp5.SUBARRAY / 2, lv + wp5.SUBARRAY / 2 + vp_prior = (b - a) / (np.interp(b, wp3.DEPTHS, zo["tv"]) - np.interp(a, wp3.DEPTHS, zo["tv"])) + win = (z >= a) & (z <= b) + c44 = 1.0 / np.mean(1.0 / (logs.rho[win] * logs.vs[win] ** 2)) + vs_b = np.sqrt(c44 / np.mean(logs.rho[win])) + vs_prior = float(phase_velocity(vs_b, 0.5 * np.interp(lv, z, zo["qz"]), 30.0, + logs.meta["f_sonic_s"])) + post = W.invert_level(d, {"vp0": (vp_prior, 0.01 * vp_prior), + "vs0": (vs_prior, 0.05 * vs_prior)}) + est.append((post["eps"], post["delta"])) + t = earth40.at(lv) + rows.append(f"{lv:.0f} m: eps {post['eps'][0]:.3f} ({float(t['eps']):.3f}), delta " + f"{post['delta'][0]:+.3f} ({float(t['delta']):+.3f})") + gx = np.array([np.interp(lv, z, g_up) for lv in wp5.LEVELS]) + cal = {} + for k, n in enumerate(("eps", "delta")): + y = np.array([e[k][0] for e in est]) + s = np.array([e[k][1] for e in est]) + cal[n] = bayes_linear(gx, y, s, prior_mean=(0.0, 0.5), prior_std=(0.1, 2.0)) + ok = np.isfinite(g_sonic) + g5 = S.running_average(z, np.interp(z, z[ok], g_sonic[ok]), 5.0) + eps_c, delta_c = cal["eps"].predict(g5)[0], cal["delta"].predict(g5)[0] + print("\n2. Walkaway (priors: Vp0 from this ZO-VSP, Vs0 from the dipole at 30 Hz), estimate (true):") + print(" " + "; ".join(rows)) + print(f" eps = {cal['eps'].coef[0]:+.3f} + {cal['eps'].coef[1]:.2f} gamma, delta = " + f"{cal['delta'].coef[0]:+.3f} + {cal['delta'].coef[1]:.2f} gamma") + ti = model.layer_index(z) + units = ((470, 600, "L. Cret."), (600, 631, "Hekkingen"), (631, 672, "Fuglen"), + (672, 702, "Stø"), (702, 740, "Nordmela"), (740, 800, "Fruholmen")) + print(" unit means, calibrated (true): " + "; ".join( + f"{n} eps {eps_c[(z >= a) & (z < b)].mean():.3f} ({model.eps[ti][(z >= a) & (z < b)].mean():.3f})" + f" delta {delta_c[(z >= a) & (z < b)].mean():+.3f} ({model.delta[ti][(z >= a) & (z < b)].mean():+.3f})" + for a, b, n in units)) + reg.log("Isotropic v1 cannot predict non-hyperbolic moveout or walkaway slownesses", + "VTI overburden (shales)", "slowness-polarization inversion at six levels, " + "gamma -> (eps, delta) calibration", + f"eps slope {cal['eps'].coef[1]:.2f}, delta slope {cal['delta'].coef[1]:.2f} per gamma", + "adopt eps, delta, gamma curves", 2) + v1 = zo["v1"].model + v2m = WF.model_from_logs(logs, paper.WATER_DEPTH, SEA_WATER.vp0, SEA_WATER.rho, qp=zo["qz"], + qs=0.5 * zo["qz"], eps=eps_c, delta=delta_c, gamma=g5, + gap_vp=zo["v_gap"], disperse_vs=False) + assert np.allclose(v2m.vp0, v1.vp0) + v2 = reg.new_version("anisotropic", v2m, "eps, delta from the calibrated sonic gamma") + return dict(v2=v2, g5=g5, eps=eps_c, delta=delta_c, cal=cal, qz=zo["qz"]) + + +# -- 3. surface seismic ----------------------------------------------------------------------- + + +def process(data, mig, q_comp, sel, z, tab): + t = np.arange(data.shape[1]) * DT + t_sb = reflection_table(mig, [mig.interfaces[1] - 0.01], OFFSETS, Z_S, Z_R).time[0] + w = np.clip((t[None, :] - (t_sb[:, None] + 0.035)) / 0.01, 0.0, 1.0) + mute = 0.5 - 0.5 * np.cos(np.pi * w) + return I.migrate_cmp((data * mute)[sel], DT, OFFSETS[sel], mig, z, Z_S, Z_R, + q_compensation="amplitude" if q_comp else None, gain_limit=10.0, + high_cut=80.0, table=tab) + + +def ava(img, z0): + return avo.pick_event(img.depth, img.data, img.theta, z0, 6.0) + + +def misfit(ref, other, scale=1.0): + g = np.radians(np.arange(FIT[0], FIT[1] + 0.01, 1.0)) + a = np.interp(g, ref.theta, ref.amplitude) + b = np.interp(g, other.theta, other.amplitude * scale) + return float(np.sqrt(np.mean((a - b) ** 2) / np.mean(a ** 2))) + + +def surface_step(model, reg, an, workers): + v2 = an["v2"].model + wav = ricker(30.0, DT) + run = dict(dt=DT, nt=NT, wavelet=wav, f_max=100.0, workers=workers) + t0 = time.time() + obs = M.gathers(model.subset(0.0, BOTTOM), OFFSETS, Z_S, Z_R, physics=OBSERVED, + attenuation=True, **run) + pred = M.gathers(v2.subset(0.0, BOTTOM), OFFSETS, Z_S, Z_R, physics={"full": {}, "Born": BORN}, + attenuation=False, **run) + print(f"\n3. Surface CMP gathers (reflectivity, 19 offsets to 1000 m, 7 physics variants): " + f"{time.time() - t0:.0f} s") + mig = v2.backus(20.0, respect=[v2.interfaces[1]]) + sel = OFFSETS <= X_AVA + z = np.arange(600.0, BOTTOM, 0.5) + tab = reflection_table(mig, z, OFFSETS[sel], Z_S, Z_R) + meas = {k: process(v, mig, True, sel, z, tab) for k, v in obs.items()} + pr = {k: process(v, mig, False, sel, z, tab) for k, v in pred.items()} + near = OFFSETS[sel] <= 300.0 + wz = (z >= 620.0) & (z <= 770.0) + stack = lambda im: np.nanmean(im.data[:, near], axis=1) # noqa: E731 + tie = T.tie_metrics(stack(meas["full"])[wz], stack(pr["full"])[wz], 0.5, max_shift=10.0) + v_res = unit_mean(mig, 640.0, 760.0, "vp0") + print(f" near-offset (<= 300 m) migrated stack vs elastic synthetic of v2, 620-770 m: " + f"correlation {tie.correlation:.2f} (best {tie.max_correlation:.2f}), shift {tie.shift:+.1f} m " + f"({2e3 * tie.shift / v_res:+.1f} ms), phase {np.degrees(tie.phase):+.0f} deg, scale " + f"{tie.scale:.2f}") + reg.check("reservoir event time shift (near-offset tie)", abs(2e3 * tie.shift / v_res), 2.0, + "<=", "ms") + out = dict(meas=meas, pred=pr, z=z, sel=sel, tab=tab, mig=mig, tie=tie, stack=stack, wz=wz, + obs=obs) + rows = {} + for name, z0 in EVENTS.items(): + z0m = z0 - 2.0 + pm, pp = ava(meas["full"], z0m), ava(pr["full"], z0m) + pmb, ppb = ava(meas["Born"], z0m), ava(pr["Born"], z0m) + rows[name] = dict(full=(pm, pp), born=(pmb, ppb)) + fm, fp = (avo.intercept_gradient(p.theta, p.amplitude, FIT) for p in (pm, pp)) + print(f" {name}: as acquired (full physics), predicted/measured gradient " + f"{fp.gradient / fm.gradient:+.2f}, AVA misfit {misfit(pm, pp):.2f}; " + f"contamination-free (Born) misfit {misfit(pmb, ppb):.2f}") + out["rows"] = rows + return out + + +def contamination(reg, sf): + """RMS differences between physics variants in the migrated reservoir windows, relative to + the Born primaries: all offsets, near (<= 300 m) and far.""" + meas, z = sf["meas"], sf["z"] + near = OFFSETS[sf["sel"]] <= 300.0 + out = {} + for name, (lo, hi) in (("Stø", (650.0, 690.0)), ("Fruholmen", (720.0, 760.0))): + wz = (z >= lo) & (z <= hi) + for part, cols in (("all", slice(None)), ("near", near), ("far", ~near)): + d = {k: im.data[wz][:, cols] for k, im in meas.items()} + born = np.sqrt(np.nanmean(d["Born"] ** 2)) + rel = lambda a: float(np.sqrt(np.nanmean(a ** 2)) / born) # noqa: E731 + out[name, part] = {"full - Born": rel(d["full"] - d["Born"]), + "P multiples": rel(d["P only"] - d["P primaries"]), + "P primaries": rel(d["P primaries"]), + "primaries - Born": rel(d["primaries"] - d["Born"]), + "full - primaries": rel(d["full"] - d["primaries"])} + print(" contamination after migration, rms relative to the Born primaries " + "(all offsets / <= 300 m / beyond):") + for key in out["Stø", "all"]: + print(f" {key:17s} " + "; ".join( + f"{n} " + " / ".join(f"{out[n, part][key]:.2f}" for part in ("all", "near", "far")) + for n in ("Stø", "Fruholmen"))) + st = sf["rows"]["Stø top"] + a, b = out["Stø", "all"], out["Fruholmen", "all"] + reg.log("Reservoir AVA beyond 20 deg erratic and not predicted (misfit " + f"{misfit(*st['full']):.2f} at the Stø top)", + "interbed multiples and mode conversion from the overburden stringers", + "compare physics variants of the true-earth gathers after identical processing", + f"full - Born = {a['full - Born']:.1f} x Born (Stø), {b['full - Born']:.1f} x (Fruholmen), " + f"growing with offset ({out['Stø', 'near']['full - Born']:.2f} at <= 300 m); " + f"without mode conversion the P-wave multiples are only {a['P multiples']:.2f}, " + f"{b['P multiples']:.2f} x Born; converted primaries and multiples each exceed the " + "total and largely cancel; contamination-free misfit " + f"{misfit(*st['born']):.2f}", + "elastic (converted-wave) contamination: an acoustic demultiple would not repair " + "the AVA; test the remaining hypotheses on contamination-free data") + return out + + +VS = "dipole Vs at 30 Hz (Qs = Qp/2)" +TRUE = "true earth" +TRUE_ANISO = "true earth, calibrated eps and delta" +TRUE_ELASTIC = "true earth, log Vp, Vs and rho" + + +def scenarios(model, logs, reg, an, sf, workers): + """Born-domain predictions: the Vs hypothesis, then (synthetic only) the truth, alone and + with one part of v2 swapped in, and the measurement re-migrated with the true Q.""" + v2 = an["v2"].model + v3m = WF.model_from_logs(logs, paper.WATER_DEPTH, SEA_WATER.vp0, SEA_WATER.rho, + qp=an["qz"], qs=0.5 * an["qz"], eps=an["eps"], delta=an["delta"], + gamma=an["g5"], gap_vp=float(v2.vp0[1]), disperse_vs=True) + true = model.subset(0.0, BOTTOM) + below = true.mid > logs.depth[0] + est = v2.at(true.mid[below]) + t_aniso, t_elastic = true.copy(), true.copy() + for k in ("eps", "delta"): + getattr(t_aniso, k)[below] = est[k] + for k in ("vp0", "vs0", "rho"): + getattr(t_elastic, k)[below] = est[k] + models = {VS: v3m.subset(0.0, BOTTOM), TRUE: true, TRUE_ANISO: t_aniso, TRUE_ELASTIC: t_elastic} + run = dict(dt=DT, nt=NT, wavelet=ricker(30.0, DT), f_max=100.0, workers=workers) + t0 = time.time() + preds = {k: M.gathers(m, OFFSETS, Z_S, Z_R, physics={"Born": BORN}, attenuation=False, + **run)["Born"] for k, m in models.items()} + print(f"\n4. Scenario predictions (Born, elastic): {time.time() - t0:.0f} s") + out = {"v2 (dipole Vs as logged)": sf["pred"]["Born"]} + for k, d in preds.items(): + out[k] = process(d, sf["mig"], False, sf["sel"], sf["z"], sf["tab"]) + table, ratio = {}, {} + for name, z0 in EVENTS.items(): + pm = sf["rows"][name]["born"][0] + row = {} + for k, im in out.items(): + p = ava(im, z0 - 2.0) + row[k] = misfit(pm, p) + ratio[name, k] = np.interp(pm.theta, p.theta, p.amplitude) / pm.amplitude + table[name] = row + print(f" {name}: AVA misfit against the contamination-free measurement: " + + "; ".join(f"{k} {v:.3f}" for k, v in row.items())) + v3 = None + gain = {k: np.mean([table[e][k] for e in EVENTS]) for k in (VS, "v2 (dipole Vs as logged)")} + adopt = gain[VS] < gain["v2 (dipole Vs as logged)"] + st, fr = table["Stø top"], table["Fruholmen top"] + reg.log(f"Contamination-free reservoir AVA still off: misfit {st['v2 (dipole Vs as logged)']:.2f} " + f"(Stø), {fr['v2 (dipole Vs as logged)']:.2f} (Fruholmen)", + "sonic and seismic see Vs differently: dipole Vs is dispersed (2 kHz vs 30 Hz)", + "Vs moved to 30 Hz with Qs = Qp/2 (the article's rock-physics Vs analogue)", + f"misfit {st[VS]:.2f} (Stø), {fr[VS]:.2f} (Fruholmen)", + "adopt" if adopt else "reject: no overall improvement", 3 if adopt else None) + if adopt: + v3 = reg.new_version("seismic-frequency Vs", v3m, "dipole Vs moved to 30 Hz (Qs = Qp/2)") + # model or processing? The true earth's own prediction (synthetic only) + th = sf["rows"]["Stø top"]["born"][0].theta + r_true = ratio["Stø top", TRUE][(th >= np.radians(FIT[0])) & (th <= np.radians(FIT[1]))] + sel = sf["sel"] + tq = model.backus(20.0, respect=[model.interfaces[1]]) + q_fix = {} + for label, upper in (("above the logs", logs.depth[0]), ("everywhere", np.inf)): + mq = sf["mig"].copy() + k = mq.mid < upper + mq.qp[k] = tq.at(mq.mid[k])["qp"] + mq.qp[0] = np.inf + tab = reflection_table(mq, sf["z"], OFFSETS[sel], Z_S, Z_R) + mb = process(sf["obs"]["Born"], mq, True, sel, sf["z"], tab) + q_fix[label] = np.mean([misfit(ava(mb, z0 - 2.0), ava(out[TRUE], z0 - 2.0)) + for z0 in EVENTS.values()]) + print(f" true-earth prediction / measurement at the Stø top, 20-44 deg: " + f"{r_true.min():.2f}-{r_true.max():.2f}; " + f"mean misfit {np.mean([table[e][TRUE] for e in EVENTS]):.3f}, with the true Q above the " + f"logs {q_fix['above the logs']:.3f}, everywhere {q_fix['everywhere']:.3f}") + reg.log(f"(synthetic only) The true earth itself predicts {100 * (np.median(r_true) - 1):+.0f}% " + "amplitude at all angles", + "Q compensation too weak: no VSP receivers between the seabed and 450 m, so the " + f"Quaternary (true Q 30) takes the first block's Q ({v2.qp[1]:.0f})", + "re-migrate the measurement with the true Q above the logs, then everywhere", + f"mean misfit {np.mean([table[e][TRUE] for e in EVENTS]):.3f} -> " + f"{q_fix['above the logs']:.3f} -> {q_fix['everywhere']:.3f}", + "explained; constrain the shallow t* (VSP source signature, seabed multiples) on " + "field data") + reg.log("(synthetic only) v2's Stø AVA falls with angle against the true earth's", + "the calibrated eps and delta, or the log resolution of the stringers", + "swap one part of v2 into the true earth", + f"Stø misfit {st[TRUE_ANISO]:.2f} with v2's eps and delta, {st[TRUE_ELASTIC]:.2f} with " + f"v2's Vp, Vs and rho; {st['v2 (dipole Vs as logged)']:.2f} with both", + "not separable: the 5 m-smoothed gamma gives the isotropic stringers shale anisotropy, " + "which v2's equally smoothed velocities partly cancel; keep one resolution") + return dict(table=table, images=out, v3=v3, ratio=ratio, q_fix=q_fix) + + +def identifiability(sf, an): + """Which single change explains the contamination-free Stø-top AVA residual, and what + the residual can resolve jointly (linearized, exact interface coefficient + Q).""" + v2 = an["v2"].model + pm, pp = sf["rows"]["Stø top"]["born"] + th = np.radians(np.arange(FIT[0], FIT[1] + 0.01, 2.0)) + resid = np.interp(th, pm.theta, pm.amplitude) / np.interp(th, pp.theta, pp.amplitude) - 1.0 + z0 = paper.STO_TOP_TVDSS + up = [unit_mean(v2, z0 - 8.0, z0, k) for k in ("vp0", "vs0", "rho", "eps", "delta")] + lo = [unit_mean(v2, z0, z0 + 8.0, k) for k in ("vp0", "vs0", "rho", "eps", "delta")] + R0 = exact_pp(th, tuple(up), tuple(lo)).real + t_star = np.interp(th, sf["tab"].theta[np.searchsorted(sf["z"], z0) - 1], + sf["tab"].tstar[np.searchsorted(sf["z"], z0) - 1]) + + def rel(side, k, step): + u, lw = list(up), list(lo) + (u if side == "above" else lw)[k] += step + return (exact_pp(th, tuple(u), tuple(lw)).real / R0 - 1.0) / step + + params = { # name: (sensitivity of the relative amplitude, prior std, unit) + "Vs below (Stø)": (rel("below", 1, 10.0) * lo[1] / 100.0, 3.0, "%"), + "Vs above (Fuglen)": (rel("above", 1, 10.0) * up[1] / 100.0, 3.0, "%"), + "rho below": (rel("below", 2, 5.0) * lo[2] / 100.0, 2.0, "%"), + "delta below": (rel("below", 4, 0.005), 0.03, ""), + "eps below": (rel("below", 3, 0.005), 0.03, ""), + "overburden Q": (-np.pi * 35.0 * t_star / 100.0, 30.0, "%"), # 1 % lower Q in the model + } + print("\n identifiability: the contamination-free Stø-top residual (measured/predicted - 1) is " + + ", ".join(f"{100 * r:+.0f}%" for r in resid[::3]) + f" at {', '.join(f'{d:.0f}' for d in np.degrees(th[::3]))} deg") + sig = 0.02 + for k, (j, _, unit) in params.items(): + dp = float(j @ resid / (j @ j)) + chi = np.sqrt(np.mean((resid - j * dp) ** 2)) / np.sqrt(np.mean(resid ** 2)) + print(f" {k:18s} alone: change {dp:+.3g}{unit}, explains the residual to " + f"{100 * chi:.0f}% of its RMS") + J = np.column_stack([v[0] for v in params.values()]) + res = avo.linear_resolution(J, [v[1] for v in params.values()], sig ** 2 * np.eye(th.size)) + names = list(params) + print(" jointly (2 % data noise per angle): data weight " + ", ".join( + f"{n} {w:.2f}" for n, w in zip(names, res.data_weight))) + i, j_ = names.index("Vs below (Stø)"), names.index("delta below") + print(f" posterior correlation Vs below / delta below {res.correlation[i, j_]:+.2f}, " + f"Vs below / overburden Q {res.correlation[i, names.index('overburden Q')]:+.2f}") + return dict(resid=resid, th=th, params=names, resolution=res) + + +# -- 5. the Fig. 5 analogue ------------------------------------------------------------------- + + +def section_step(): + xs = np.arange(0.0, 1001.0, 20.0) + wav = ricker(30.0, DT) + oil, brine = dipping_reservoir_models(xs), dipping_reservoir_models(xs, oil=False) + nt = 1000 + sec = {k: np.array([normal_incidence(m, DT, nt, wav, attenuation=True, f_max=100.0)[1] + for m in ms]) for k, ms in (("oil", oil), ("brine", brine))} + rng = np.random.default_rng(4) + t = np.arange(nt) * DT + noise = rng.standard_normal(sec["oil"].shape) + f = np.fft.rfftfreq(nt, DT) + spec = np.fft.rfft(noise) + spec[..., f > 100.0] = 0.0 + noise = np.fft.irfft(spec, n=nt) + win = (t > 0.7) & (t < 0.85) + noise *= 0.05 * np.sqrt(np.mean(sec["oil"][:, win] ** 2)) / noise.std() + z = np.arange(600.0, 760.0, 0.5) + from scipy.signal import hilbert + + def image(data, iso, qc): + out = np.empty((xs.size, z.size)) + for j, m in enumerate(oil): + mm = m.backus(20.0, respect=[m.interfaces[1]]) + if iso: + mm = mm.copy() + mm.vp0 = mm.vp0 * np.sqrt(1.0 + 2.0 * mm.delta) + mm.eps[:] = 0.0 + mm.delta[:] = 0.0 + out[j] = I.migrate_cmp(data[j][None, :], DT, [0.0], mm, z, 0.0, 0.0, + q_compensation="amplitude" if qc else None, gain_limit=10.0, + high_cut=80.0, true_amplitude=False).data[:, 0] + return out + + col = OWC_TVDSS - (paper.STO_TOP_TVDSS + xs * np.tan(np.radians(1.2))) + variants = {"VTI + Q compensation": (False, True), "VTI, no Q compensation": (False, False), + "isotropic (NMO v) + Q compensation": (True, True), + "isotropic (NMO v), no Q compensation": (True, False)} + print("\n5. Dipping reservoir, flat oil-water contact at 694.5 m (oil column 22 -> 2 m over 1 km), " + "zero-offset images, 5 % noise") + images, fluid, metrics = {}, {}, {} + for name, (iso, qc) in variants.items(): + so, sb, sn = image(sec["oil"] + noise, iso, qc), image(sec["brine"], iso, qc), image(noise, iso, qc) + diff = so - sb - sn # the noise-free fluid signal + env = np.abs(hilbert(diff, axis=1)) + thick = np.flatnonzero(col >= 12.0) + pk = [z[np.argmax(env[j])] for j in thick] + width = [] + for j in thick: + e = env[j] / env[j].max() + width.append(np.ptp(z[e > 0.5])) + snr = float(np.max(np.abs(diff)) / np.sqrt(np.mean(sn ** 2))) + metrics[name] = (np.mean(pk), np.ptp(pk), np.mean(width), snr) + images[name], fluid[name] = so, diff + print(f" {name:38s} fluid signal at {np.mean(pk):.1f} m (spread {np.ptp(pk):.1f} m), " + f"half-max thickness {np.mean(width):.1f} m, peak / noise RMS {snr:.1f}") + return dict(xs=xs, z=z, images=images, fluid=fluid, metrics=metrics, + time=sec["oil"] + noise, t=t) + + +# -- figures ---------------------------------------------------------------------------------- + + +def figures(zo, sf, sc, sec): + try: + import matplotlib + + matplotlib.use("Agg") + import matplotlib.pyplot as plt + + from iwt.plotting import PALETTE + except ImportError: + print("\nmatplotlib not installed: no figures") + return + FIG.mkdir(parents=True, exist_ok=True) + col = PALETTE["series"] + + def style(ax): + ax.set_facecolor(PALETTE["surface"]) + ax.grid(color=PALETTE["grid"], lw=0.5) + + # Fig. 2 analogue + fig, axes = plt.subplots(1, 2, figsize=(9.5, 6.0)) + fig.patch.set_facecolor(PALETTE["surface"]) + t, w = zo["t"], zo["w"] + ax = axes[0] + sc_ = np.max(np.abs(zo["corridor"][w])) + ax.plot(zo["corridor"] / sc_, t, color=col[0], lw=1.1, label="VSP corridor stack") + ax.plot(zo["syn"] / np.max(np.abs(zo["syn"][w])) - 1.2, t, color=col[1], lw=1.1, + label="synthetic, v1 (shifted left)") + ax.set_ylim(0.95, 0.6) + ax.set_xlabel("amplitude") + ax.set_ylabel("two-way time (s)") + ax.set_title(f"ZO-VSP tie: correlation {zo['tie'].correlation:.2f}, shift " + f"{1e3 * zo['tie'].shift:+.1f} ms", fontsize=9, loc="left") + ax.legend(fontsize=7, frameon=False, loc="lower left") + ax = axes[1] + z, wz = sf["z"], sf["wz"] + m_, p_ = sf["stack"](sf["meas"]["full"]), sf["stack"](sf["pred"]["full"]) + s = np.nanmax(np.abs(m_[wz])) + ax.plot(m_ / s, z, color=col[0], lw=1.1, label="surface, near offsets (Q-PSDM, v2)") + ax.plot(p_ * sf["tie"].scale / s - 1.2, z, color=col[1], lw=1.1, + label="elastic synthetic of v2 (shifted left)") + for name, z0 in EVENTS.items(): + ax.axhline(z0, color=PALETTE["ink2"], lw=0.6, ls="--") + ax.annotate(name, (0.62, z0 - 1.5), fontsize=7, color=PALETTE["ink2"]) + ax.set_ylim(770.0, 620.0) + ax.set_xlabel("amplitude") + ax.set_ylabel("depth (m)") + ax.set_title(f"Surface tie: correlation {sf['tie'].correlation:.2f}, shift " + f"{round(sf['tie'].shift, 1) + 0.0:+.1f} m", fontsize=9, loc="left") + ax.legend(fontsize=7, frameon=False, loc="lower left") + for ax in axes: + style(ax) + fig.tight_layout() + fig.savefig(FIG / "wp7_tie.png", dpi=150) + + # Fig. 4 analogue + fig, axes = plt.subplots(2, 2, figsize=(11.0, 7.5)) + fig.patch.set_facecolor(PALETTE["surface"]) + for c, (name, z0) in enumerate(EVENTS.items()): + pm, pp = sf["rows"][name]["full"] + ax = axes[0, c] + ax.plot(np.degrees(pm.theta), pm.amplitude, "o-", color=PALETTE["ink"], ms=3, lw=1.6, + label="measured (as acquired)") + ax.plot(np.degrees(pp.theta), pp.amplitude, "s--", color=col[1], ms=3, lw=1.2, + label="elastic synthetic, v2 (full physics)") + ax.set_title(f"{name}: full physics (multiples, mode conversion)", fontsize=9, loc="left") + ax = axes[1, c] + pm = sf["rows"][name]["born"][0] + ax.plot(np.degrees(pm.theta), pm.amplitude, "o-", color=PALETTE["ink"], ms=3, lw=1.6, + label="measured, contamination-free") + for k, (key, im) in enumerate(sc["images"].items()): + p = ava(im, z0 - 2.0) + ax.plot(np.degrees(p.theta), p.amplitude, ("s--", "^-.", "v:", "d--", "x-")[k], + color=(tuple(col) + (PALETTE["muted"],))[k], + ms=3, lw=1.1, label=key) + ax.set_title(f"{name}: contamination-free (Born)", fontsize=9, loc="left") + for ax in axes.ravel(): + ax.axvspan(*FIT, color=PALETTE["grid"], alpha=0.5, zorder=0) + ax.set_xlabel("reflection angle (deg)") + ax.set_ylabel("migrated amplitude") + ax.legend(fontsize=6.5, frameon=False) + style(ax) + fig.tight_layout() + fig.savefig(FIG / "wp7_ava.png", dpi=150) + + # Fig. 5 analogue + panels = {"oil section, time,\nno Q compensation": None, + "oil section,\nVTI + Q compensation": sec["images"]["VTI + Q compensation"]} + panels.update({f"fluid signal (oil - brine),\n{k}": sec["fluid"][k] + for k in ("VTI + Q compensation", "isotropic (NMO v) + Q compensation")}) + fig, axes = plt.subplots(1, 4, figsize=(15.0, 4.8)) + fig.patch.set_facecolor(PALETTE["surface"]) + from matplotlib.colors import LinearSegmentedColormap + + cmap = LinearSegmentedColormap.from_list("d", ("#1c5cab", "#f0efec", "#c23b3a")) + xs, zz = sec["xs"], sec["z"] + for ax, (name, im) in zip(axes, panels.items()): + if im is None: + tsel = (sec["t"] > 0.72) & (sec["t"] < 0.88) + d = sec["time"][:, tsel] + clip = np.percentile(np.abs(d), 99) + ax.imshow(d.T, aspect="auto", cmap=cmap, vmin=-clip, vmax=clip, + extent=[xs[0], xs[-1], sec["t"][tsel][-1], sec["t"][tsel][0]]) + ax.set_ylabel("two-way time (s)") + else: + clip = np.max(np.abs(im)) + ax.imshow(im.T, aspect="auto", cmap=cmap, vmin=-clip, vmax=clip, + extent=[xs[0], xs[-1], zz[-1], zz[0]]) + ax.axhline(OWC_TVDSS, color=PALETTE["ink2"], lw=0.7, ls="--") + ax.set_ylim(740.0, 640.0) + ax.set_ylabel("depth (m)") + ax.set_title(name, fontsize=8.5, loc="left") + ax.set_xlabel("distance (m)") + fig.tight_layout() + fig.savefig(FIG / "wp7_section.png", dpi=150) + print("\nfigures written to", FIG) + + +def main(): + ap = argparse.ArgumentParser() + ap.add_argument("--workers", type=int, default=4) + args = ap.parse_args() + t_start = time.time() + model, logs = wisting_like_well() + reg = WF.Register() + zo = zovsp_step(model, logs, reg, args.workers) + an = anisotropy_step(model, logs, reg, zo) + sf = surface_step(model, reg, an, args.workers) + contamination(reg, sf) + sc = scenarios(model, logs, reg, an, sf, args.workers) + identifiability(sf, an) + sec = section_step() + pm, pp = sf["rows"]["Stø top"]["born"] + fm, fp = (avo.intercept_gradient(p.theta, p.amplitude, FIT) for p in (pm, pp)) + reg.check("Stø-top AVA gradient, contamination-free, |predicted/measured - 1|", + abs(fp.gradient / fm.gradient - 1.0), 0.10, "<=") + pm, pp = sf["rows"]["Stø top"]["full"] + fm, fp = (avo.intercept_gradient(p.theta, p.amplitude, FIT) for p in (pm, pp)) + reg.check("Stø-top AVA gradient, as acquired, |predicted/measured - 1|", + abs(fp.gradient / fm.gradient - 1.0), 0.10, "<=") + m = sec["metrics"] + reg.log("Flat oil-water contact hard to see under a dipping reservoir (Fig. 5)", + "Q-PSDM sharpens it; anisotropy positions it", + "zero-offset images of an oil and a brine section, VTI vs isotropic, with and without Q", + f"VTI places the fluid signal at {m['VTI + Q compensation'][0]:.0f} m, isotropic at " + f"{m['isotropic (NMO v) + Q compensation'][0]:.0f} m; Q compensation thins it from " + f"{m['VTI, no Q compensation'][2]:.1f} to {m['VTI + Q compensation'][2]:.1f} m", + "anisotropy is what matters in 1D; the image improvement needs 2D/3D imaging") + text = ("# Discrepancy register of `examples/reproduce_paper.py`\n\n" + "Generated by the script on the synthetic well; do not edit by hand.\n\n" + + reg.to_markdown()) + (ROOT / "docs" / "discrepancy_register.md").write_text(text) + print("\nAcceptance criteria:") + for c in reg.criteria: + print(f" {c.name}: {c.value:.3g} {c.unit} ({c.sense} {c.threshold:g}) " + f"{'pass' if c.passed else 'FAIL'}") + print(f"\nregister written to docs/discrepancy_register.md ({len(reg.versions)} versions, " + f"{len(reg.discrepancies)} discrepancies); total {time.time() - t_start:.0f} s") + figures(zo, sf, sc, sec) + + +if __name__ == "__main__": + main() diff --git a/iwt/__init__.py b/iwt/__init__.py index 1a01c72..902746a 100644 --- a/iwt/__init__.py +++ b/iwt/__init__.py @@ -17,11 +17,15 @@ package 5 exact qP ray tracing (:mod:`iwt.raytrace`), the walkaway VSP inversion (:mod:`iwt.vsp.walkaway`) and the gamma -> (eps, delta) calibration, and work package 6 Q-compensated depth migration of CMP gathers (:mod:`iwt.imaging`) with -AVA extraction (:mod:`iwt.avo`). +AVA extraction (:mod:`iwt.avo`), and work package 7 the integrated loop: well-tie +metrics (:mod:`iwt.welltie`), versioned models and the discrepancy register +(:mod:`iwt.workflow`) and the contamination by multiples and mode conversion +(:mod:`iwt.multiples`). Plotting (:mod:`iwt.plotting`) needs matplotlib. """ -from . import avo, calibration, imaging, paper, qmodel, raytrace, rockphysics, sonic +from . import (avo, calibration, imaging, multiples, paper, qmodel, raytrace, rockphysics, + sonic, welltie, workflow) from .model import ( ElasticModel1D, Stiffness, @@ -34,17 +38,20 @@ ) from .synthetic import WellLogs, synthetic_logs, wisting_like_model, wisting_like_well -__version__ = "0.6.0" +__version__ = "0.7.0" __all__ = [ "avo", "calibration", "imaging", + "multiples", "paper", "qmodel", "raytrace", "rockphysics", "sonic", + "welltie", + "workflow", "ElasticModel1D", "Stiffness", "backus_average", diff --git a/iwt/avo.py b/iwt/avo.py index 277eacd..c8f868d 100644 --- a/iwt/avo.py +++ b/iwt/avo.py @@ -4,7 +4,9 @@ offset), taking the extreme of a fixed polarity in a window that moves with the event, so residual moveout from a wrong migration model is tracked rather than cut. :func:`intercept_gradient` fits the two-term form ``R = A + B sin^2(theta)`` -(Shuey, 1985) by least squares over an angle range. +(Shuey, 1985) by least squares over an angle range. :func:`linear_resolution` +says which model changes such data can tell apart (the identifiability study of +the plan): Gaussian linear inversion of the sensitivities with prior widths. References ---------- @@ -18,7 +20,8 @@ import numpy as np -__all__ = ["Pick", "pick_event", "AVAFit", "intercept_gradient"] +__all__ = ["Pick", "pick_event", "AVAFit", "intercept_gradient", "Resolution", + "linear_resolution"] class Pick(NamedTuple): @@ -109,3 +112,28 @@ def intercept_gradient(theta, amplitude, theta_range: Sequence[float] = (0.0, 90 scale = 1.0 if sigma is not None else (r @ r) / max(n - 2, 1) cov = np.linalg.inv((G / s[:, None]).T @ (G / s[:, None])) * scale return AVAFit(float(coef[0]), float(coef[1]), cov, float(np.sqrt(np.mean(r * r))), n) + + +class Resolution(NamedTuple): + """Linear (Gaussian) resolution of parameters by data: posterior std, data weight, correlation.""" + + posterior_std: np.ndarray + data_weight: np.ndarray + correlation: np.ndarray + + +def linear_resolution(J, prior_std, data_cov) -> Resolution: + """What linearized data ``d = J m`` with covariance ``data_cov`` resolve about ``m``. + + ``prior_std`` are the prior standard deviations of the parameters. The + data weight ``1 - posterior variance / prior variance`` is near 0 for a + parameter the data cannot see; a correlation near +-1 between two + parameters means the data can explain a change of one by the other. + """ + J = np.atleast_2d(np.asarray(J, dtype=float)) + s = np.asarray(prior_std, dtype=float) + Cd = np.atleast_2d(np.asarray(data_cov, dtype=float)) + H = J.T @ np.linalg.solve(Cd, J) + np.diag(1.0 / s ** 2) + C = np.linalg.inv(H) + sd = np.sqrt(np.diag(C)) + return Resolution(sd, 1.0 - sd ** 2 / s ** 2, C / np.outer(sd, sd)) diff --git a/iwt/modeling/reflection.py b/iwt/modeling/reflection.py index 29fde66..b59fc5e 100644 --- a/iwt/modeling/reflection.py +++ b/iwt/modeling/reflection.py @@ -34,7 +34,8 @@ from ..model import phase_velocities, thomsen_to_cij from .elastic import eigenvectors, interface_scattering, vertical_slowness -__all__ = ["zoeppritz_pp", "exact_pp", "aki_richards", "shuey", "ruger", "slowness_from_angle"] +__all__ = ["zoeppritz_pp", "exact_pp", "exact_pp_slowness", "aki_richards", "shuey", "ruger", + "slowness_from_angle"] def zoeppritz_pp(theta, vp1, vs1, rho1, vp2, vs2, rho2): @@ -70,7 +71,18 @@ def exact_pp(theta, upper, lower): the phase angle of the incident qP wave. Both media must be solids. """ theta = np.atleast_1d(np.asarray(theta, dtype=float)) - p = slowness_from_angle(theta, upper[0], upper[1], upper[3], upper[4]) + 0j + return exact_pp_slowness(slowness_from_angle(theta, upper[0], upper[1], upper[3], upper[4]), + upper, lower) + + +def exact_pp_slowness(p, upper, lower): + """:func:`exact_pp` as a function of the horizontal slowness ``p`` (s/m). + + Beyond the critical slowness of either medium the coefficient is complex + (``exp(-i omega t)`` convention); beyond the upper medium's own horizontal + slowness the incident wave is evanescent and the value is formal. + """ + p = np.atleast_1d(np.asarray(p, dtype=float)) + 0j Ds = [] for vp0, vs0, rho, eps, delta in (upper, lower): c = thomsen_to_cij(vp0, vs0, rho, eps, delta, 0.0) diff --git a/iwt/multiples.py b/iwt/multiples.py new file mode 100644 index 0000000..f6ba37b --- /dev/null +++ b/iwt/multiples.py @@ -0,0 +1,74 @@ +"""Interbed-multiple and mode-conversion contamination of reservoir reflections. + +The article found the shallow Wisting reservoir reflections "contaminated by +interbed multiples originating in the overburden" by full-waveform modelling. +Here the contamination is measured by switching physics on and off in the +reflectivity method (:mod:`iwt.modeling.reflectivity`) and comparing the gathers +after identical processing (mute, migration, AVA picking): + +* ``full``: internal multiples and P-S mode conversion; +* ``primaries``: single reflections, with every mode conversion on the way + (P down, S up and back to P at the seabed, for example) and the transmission + and conversion losses; +* ``no mode conversion``: multiples of P waves only (exact PP coefficients, the + converted energy dropped); +* ``primaries with transmission loss``: P-wave primaries only; +* ``Born primaries``: single scattering, unit transmission. + +The contamination is "full minus Born primaries". Differences between the +intermediate variants do not isolate mechanisms in finely layered media. +Short-period multiples restore much of the amplitude that transmission losses +remove (O'Doherty & Anstey, 1971), so "full minus primaries" mostly measures +that restoration; with mode conversion on, the converted primaries and the +multiples are each larger than their sum and largely cancel. Dropping mode +conversion is not a physical medium: at oblique incidence the P-only stack loses +the converted energy at every thin layer without getting any of it back. It +still bounds the part of the contamination that P-wave interbed multiples alone +could cause. + +References +---------- +O'Doherty, R. F., & Anstey, N. A., 1971, Reflections on amplitudes: Geophysical +Prospecting, 19(3), 430-458. https://doi.org/10.1111/j.1365-2478.1971.tb00610.x +""" + +from __future__ import annotations + +from typing import Dict, Optional + +import numpy as np + +from .model import ElasticModel1D +from .modeling.reflectivity import Reflectivity + +__all__ = ["PHYSICS", "gathers", "relative_rms"] + +PHYSICS: Dict[str, dict] = { + "full": {}, + "primaries": dict(internal_multiples=False), + "no mode conversion": dict(converted_waves=False), + "primaries with transmission loss": dict(internal_multiples=False, converted_waves=False), + "Born primaries": dict(internal_multiples=False, converted_waves=False, + transmission_loss=False), +} + + +def gathers(model: ElasticModel1D, offsets, z_source: float, z_receiver: float, dt: float, + nt: int, wavelet, physics: Optional[Dict[str, dict]] = None, + attenuation: bool = True, f_max: Optional[float] = None, + workers: int = 1) -> Dict[str, np.ndarray]: + """Hydrophone CMP gathers ``data[offset, time]`` for each physics variant (no free surface).""" + out = {} + for name, kw in (PHYSICS if physics is None else physics).items(): + r = Reflectivity(model, attenuation=attenuation, free_surface=False, **kw) + out[name] = r.record(offsets, z_source, [z_receiver], dt, nt, wavelet, "pressure", + f_max=f_max, workers=workers).data[0] + return out + + +def relative_rms(a, b, mask=None) -> float: + """``rms(a - b) / rms(b)`` over ``mask`` (NaN ignored).""" + a, b = np.asarray(a, dtype=float), np.asarray(b, dtype=float) + if mask is not None: + a, b = a[mask], b[mask] + return float(np.sqrt(np.nanmean((a - b) ** 2) / np.nanmean(b * b))) diff --git a/iwt/raytrace.py b/iwt/raytrace.py index e9b2480..32461f0 100644 --- a/iwt/raytrace.py +++ b/iwt/raytrace.py @@ -269,6 +269,17 @@ def reflection_table(model: ElasticModel1D, depths, offsets, z_source: float = 0 p_lim = np.where(legs > 0, np.sin(th_max) / qp_phase_velocity(th_max, *A), np.inf) k_out = np.searchsorted(m.interfaces, z_out) # layers [0, k) lie above z_out n_lay = int(k_out.max()) + if not np.any(x > 0.0): # zero offset: vertical sums, exact + q, _, d2q = qp_vertical_slowness(0.0, *A[:, :n_lay], derivative=2) + w = legs[:n_lay] * h[:n_lay] + T = np.concatenate([[0.0], np.cumsum(w * q)])[k_out] + ts = np.concatenate([[0.0], np.cumsum(w * q * qinv[:n_lay])])[k_out] + dX = np.concatenate([[0.0], np.cumsum(-w * d2q)])[k_out] + i_s = int(np.searchsorted(m.interfaces, zs, side="right")) - 1 + q_s = qp_vertical_slowness(0.0, *A[:, i_s]) + one = np.ones((1, x.size)) + return ReflectionTable(z_out, x, T[:, None] * one, 0.0 * T[:, None] * one, + 0.0 * T[:, None] * one, ts[:, None] * one, (q_s * dX)[:, None] * one) p_cap = np.minimum.accumulate(p_lim[:n_lay]) # smallest limit above each depth i_src = int(np.searchsorted(m.interfaces, zs, side="right")) - 1 diff --git a/iwt/synthetic.py b/iwt/synthetic.py index f028a54..40f8b5e 100644 --- a/iwt/synthetic.py +++ b/iwt/synthetic.py @@ -63,6 +63,7 @@ "OWC_TVDSS", "STO_PHI", "wisting_like_model", + "dipping_reservoir_models", "synthetic_logs", "wisting_like_well", "WASHOUT", @@ -258,6 +259,49 @@ def wisting_like_model(dz: float = paper.MODEL_SAMPLING, td: float = paper.TD_TV return ElasticModel1D(interfaces, label=labels, f_ref=f_ref, **full).validate() +def dipping_reservoir_models(x, dip_deg: float = 1.2, wedge: str = "Hekkingen", + td: float = 800.0, dz: float = paper.MODEL_SAMPLING, + seed: int = 2016, f_ref: float = 30.0, oil: bool = True): + """Columns of a 2D section with a dipping reservoir and a flat oil-water contact. + + The ``wedge`` unit thickens by ``x tan(dip_deg)`` (x in m), so everything + below it, stringers included, dips; the oil-water contact stays at + :data:`OWC_TVDSS`, so the oil column thins with x and disappears where top + Stø reaches the contact. The wedge is an anisotropic shale, so an + isotropic velocity model mispositions the deeper events by an amount that + grows with x. Returns one :class:`~iwt.model.ElasticModel1D` per x (down to + ``td``), with the same fine-scale fluctuation seed. ``oil=False`` fills the + whole reservoir with brine (the reference for the fluid signal). + """ + names = [u.name for u in WISTING_LIKE_UNITS] + k = names.index(wedge) + wedge_top = WISTING_LIKE_UNITS[k].top + out = [] + for shift in np.asarray(x, dtype=float) * np.tan(np.radians(dip_deg)): + units = [] + for j, u in enumerate(WISTING_LIKE_UNITS): + if j <= k: + units.append(u) + elif u is _STO_WATER: + water = replace(u, top=(max(OWC_TVDSS, _STO_OIL.top + shift) if oil + else _STO_OIL.top + shift)) + units.append(water) + elif u is _STO_OIL: + if oil and u.top + shift < OWC_TVDSS - dz: + units.append(replace(u, top=u.top + shift)) + else: + units.append(replace(u, top=u.top + shift)) + units.sort(key=lambda u: u.top) + below = units[units.index(water) + 1:] + if below and below[0].top <= water.top: + raise ValueError("the dip moves the base of the reservoir above the contact") + stringers = [replace(st, top=st.top + shift) if st.top > wedge_top else st + for st in STRINGERS] + out.append(wisting_like_model(dz=dz, td=td, units=units, stringers=stringers, + seed=seed, f_ref=f_ref)) + return out + + @dataclass class WellLogs: """Wireline logs in SI units, indexed by depth (m TVDSS). diff --git a/iwt/welltie.py b/iwt/welltie.py new file mode 100644 index 0000000..2bb142d --- /dev/null +++ b/iwt/welltie.py @@ -0,0 +1,115 @@ +"""Well tie: wavelet estimation and tie metrics. + +:func:`estimate_wavelet` finds the wavelet that, convolved with a reflectivity +series (or with a synthetic made with a unit spike), best predicts the seismic +trace(s), by damped least squares over a centred window of lags. It is the +time-domain form of the matched-filter (Wiener) estimate; Walden & White (1998) +give the frequency-domain treatment of the same noisy input-output problem. +Several traces (angles, offsets, nearby traces) can be used together. + +:func:`tie_metrics` reports what a tie is judged by (White & Simm, 2003): the +correlation at zero lag and its maximum over lags (the time shift, refined to a +fraction of a sample), the constant phase rotation that maximizes the +correlation, the least-squares amplitude scale, and the proportion of the trace +energy the scaled synthetic predicts (PEP). + +References +---------- +Walden, A. T., & White, R. E., 1998, Seismic wavelet estimation: A frequency +domain solution to a geophysical noisy input-output problem: IEEE Transactions on +Geoscience and Remote Sensing, 36(1), 287-297. https://doi.org/10.1109/36.655337 +White, R. E., & Simm, R., 2003, Tutorial: Good practice in well ties: First Break, +21(10). https://doi.org/10.3997/1365-2397.21.10.25640 +""" + +from __future__ import annotations + +from typing import NamedTuple + +import numpy as np + +__all__ = ["estimate_wavelet", "TieMetrics", "tie_metrics", "rotate_phase"] + + +def estimate_wavelet(trace, reflectivity, dt: float, half_length: float, + damping: float = 1e-3): + """Least-squares wavelet ``(t, w)`` on lags ``-half_length .. +half_length`` (s). + + ``trace`` and ``reflectivity`` are ``[n_traces, nt]`` (or 1D) on the same + time axis; the wavelet minimizes ``|trace - reflectivity * w|^2 + + damping * |R|^2 |w|^2`` summed over traces, with ``|R|^2`` the reflectivity + energy (so the damping is relative). + """ + d = np.atleast_2d(np.asarray(trace, dtype=float)) + r = np.atleast_2d(np.asarray(reflectivity, dtype=float)) + L = int(round(half_length / dt)) + lags = np.arange(-L, L + 1) + nt = d.shape[1] + rows = [] + for k in lags: # column k: reflectivity shifted by k samples + sh = np.zeros_like(r) + if k >= 0: + sh[:, k:] = r[:, :nt - k] + else: + sh[:, :k] = r[:, -k:] + rows.append(sh.ravel()) + G = np.array(rows).T + A = G.T @ G + A = A + damping * np.trace(A) / A.shape[0] * np.eye(A.shape[0]) + w = np.linalg.solve(A, G.T @ d.ravel()) + return lags * dt, w + + +def rotate_phase(trace, phi: float): + """Constant phase rotation by ``phi`` (rad) via the analytic signal.""" + from scipy.signal import hilbert + + a = hilbert(np.asarray(trace, dtype=float)) + return np.real(a * np.exp(-1j * phi)) + + +class TieMetrics(NamedTuple): + """Correlation at zero lag, envelope-peak correlation, its lag (s), phase (rad), scale, PEP.""" + + correlation: float + max_correlation: float + shift: float + phase: float + scale: float + pep: float + + +def tie_metrics(observed, synthetic, dt: float, max_shift: float = 0.02) -> TieMetrics: + """Tie metrics of ``synthetic`` against ``observed`` (same time axis and window). + + For a band-limited wavelet a constant phase rotation looks much like a time + shift, so both come from the complex cross-correlation of the analytic + signals: its envelope peaks at the time shift (positive when the synthetic + must be delayed), and its phase there is the constant rotation that turns + the shifted synthetic into the observed trace (:func:`rotate_phase`). + ``max_correlation`` is the envelope value there (the correlation after the + best shift and rotation). ``scale`` and ``pep`` refer to the unshifted, + unrotated synthetic. + """ + from scipy.signal import hilbert + + o = np.asarray(observed, dtype=float) + s = np.asarray(synthetic, dtype=float) + no, ns = np.linalg.norm(o), np.linalg.norm(s) + oa, sa = hilbert(o), hilbert(s) + L = int(round(max_shift / dt)) + lags = np.arange(-L, L + 1) + # |analytic| norms are twice the real ones (equal energy in the quadrature part) + cc = np.array([np.vdot(np.roll(sa, k), oa) for k in lags]) / (2.0 * no * ns) + env = np.abs(cc) + i = int(np.argmax(env)) + off = 0.0 + if 0 < i < lags.size - 1: + a, b, c = env[i - 1], env[i], env[i + 1] + den = a - 2 * b + c + off = 0.5 * (a - c) / den if den != 0 else 0.0 + phi = float(-np.angle(cc[i])) + scale = float(o @ s / (ns * ns)) + pep = float(1.0 - np.sum((o - scale * s) ** 2) / np.sum(o * o)) + return TieMetrics(float(o @ s / (no * ns)), float(env[i]), float((lags[i] + off) * dt), + phi, scale, pep) diff --git a/iwt/workflow.py b/iwt/workflow.py new file mode 100644 index 0000000..fd924e2 --- /dev/null +++ b/iwt/workflow.py @@ -0,0 +1,179 @@ +"""The integrated loop: versioned earth models, acceptance criteria and the discrepancy register. + +The article's workflow (Fig. 1) refines one earth model with every data set "until +the predictions match the measurements, or until the remaining discrepancies are +explained". That stopping rule is qualitative. Here the loop is made auditable: + +* :class:`Register` keeps every model version with what changed and why; +* each discrepancy is an entry with the observation, the hypothesis, the test, + its result and the decision, and the model version it produced (if any); +* acceptance criteria are fixed in advance (plan, section 5) and evaluated on + the final version. + +:func:`model_from_logs` builds the 0.5 m elastic model from wireline logs, the +starting point of the loop: sonic and dipole velocities moved from the tool +frequencies to the seismic frequency with a given Q (Kjartansson, 1979), optional +anisotropy and Q curves, the sea, and a gap layer between the seabed and the top +of the logs. + +References +---------- +Kjartansson, E., 1979, Constant Q-wave propagation and attenuation: JGR, 84(B9), +4737-4748. https://doi.org/10.1029/JB084iB09p04737 +Lewis, O. J., et al., 2016, Integrated processing of borehole and surface seismic +measurements for enhanced well tie and AVO analysis: 78th EAGE Conference & +Exhibition. https://doi.org/10.3997/2214-4609.201601553 +""" + +from __future__ import annotations + +from dataclasses import dataclass, field +from typing import List, Optional + +import numpy as np + +from .model import ElasticModel1D +from .qmodel import phase_velocity + +__all__ = ["model_from_logs", "ModelVersion", "Discrepancy", "Criterion", "Register"] + + +def _at(values, n): + v = np.asarray(values, dtype=float) + return np.broadcast_to(v, (n,)).copy() if v.ndim == 0 else v + + +def model_from_logs(logs, water_depth: float, water_velocity: float, water_density: float, + qp=np.inf, qs=np.inf, eps=0.0, delta=0.0, gamma=0.0, + f_seismic: float = 30.0, disperse_vs: bool = True, + gap_length: float = 10.0, gap_vp: Optional[float] = None) -> ElasticModel1D: + """Elastic model at the log sampling from sonic, dipole and density logs. + + ``qp``, ``qs``, ``eps``, ``delta`` and ``gamma`` are scalars or curves at + the log depths. The sonic Vp is moved from the monopole frequency and the + dipole Vs from the dipole frequency (``logs.meta``) to ``f_seismic`` with + the given Q (no change where Q is infinite); ``disperse_vs=False`` keeps the + dipole Vs as logged. Layers lie between log samples (slowness and density + averaged). The gap between the seabed and the top log takes the mean of + the first ``gap_length`` m of logs, isotropic, with the first Q values; + ``gap_vp`` replaces its P velocity (for example from a VSP, the logs being + blind there). + """ + d = np.asarray(logs.depth, dtype=float) + n = d.size + qp, qs = _at(qp, n), _at(qs, n) + vp = phase_velocity(np.asarray(logs.vp, float), qp, f_seismic, logs.meta["f_sonic_p"]) + vs = np.asarray(logs.vs, float) + if disperse_vs: + vs = phase_velocity(vs, qs, f_seismic, logs.meta["f_sonic_s"]) + mid = lambda x: 0.5 * (x[1:] + x[:-1]) # noqa: E731 + s_p, s_s, rho = mid(1.0 / vp), mid(1.0 / vs), mid(np.asarray(logs.rho, float)) + curves = {k: mid(_at(v, n)) for k, v in (("eps", eps), ("delta", delta), ("gamma", gamma), + ("qp", qp), ("qs", qs))} + top = d < d[0] + gap_length + gap = (np.mean(vp[top]) if gap_vp is None else float(gap_vp), np.mean(vs[top]), + np.mean(np.asarray(logs.rho)[top])) + return ElasticModel1D( + np.concatenate([[0.0, water_depth], d]), + vp0=np.concatenate([[water_velocity, gap[0]], 1.0 / s_p]), + vs0=np.concatenate([[0.0, gap[1]], 1.0 / s_s]), + rho=np.concatenate([[water_density, gap[2]], rho]), + eps=np.concatenate([[0.0, 0.0], curves["eps"]]), + delta=np.concatenate([[0.0, 0.0], curves["delta"]]), + gamma=np.concatenate([[0.0, 0.0], curves["gamma"]]), + qp=np.concatenate([[np.inf, qp[0]], curves["qp"]]), + qs=np.concatenate([[np.inf, qs[0]], curves["qs"]]), + f_ref=f_seismic) + + +@dataclass +class ModelVersion: + """One state of the earth model and the reason for it.""" + + number: int + name: str + model: ElasticModel1D + changes: str + + +@dataclass +class Discrepancy: + """An entry of the discrepancy register.""" + + observation: str + hypothesis: str + test: str + result: str + decision: str + version: Optional[int] = None + + +@dataclass +class Criterion: + """An acceptance criterion fixed before tuning: ``value`` must satisfy ``sense threshold``.""" + + name: str + value: float + threshold: float + sense: str = "<=" + unit: str = "" + + @property + def passed(self) -> bool: + if self.sense == "<=": + return bool(self.value <= self.threshold) + if self.sense == ">=": + return bool(self.value >= self.threshold) + raise ValueError("sense must be '<=' or '>='") + + +@dataclass +class Register: + """Model versions, discrepancies and acceptance criteria of one run of the loop.""" + + versions: List[ModelVersion] = field(default_factory=list) + discrepancies: List[Discrepancy] = field(default_factory=list) + criteria: List[Criterion] = field(default_factory=list) + + def new_version(self, name: str, model: ElasticModel1D, changes: str) -> ModelVersion: + v = ModelVersion(len(self.versions), name, model, changes) + self.versions.append(v) + return v + + @property + def current(self) -> ModelVersion: + return self.versions[-1] + + def log(self, observation: str, hypothesis: str, test: str, result: str, decision: str, + version: Optional[int] = None) -> Discrepancy: + entry = Discrepancy(observation, hypothesis, test, result, decision, version) + self.discrepancies.append(entry) + return entry + + def check(self, name: str, value: float, threshold: float, sense: str = "<=", + unit: str = "") -> Criterion: + c = Criterion(name, float(value), float(threshold), sense, unit) + c.passed # validates ``sense`` + self.criteria.append(c) + return c + + def to_markdown(self) -> str: + """The register as Markdown: versions, discrepancies, acceptance criteria.""" + def row(*cells): + return "| " + " | ".join(str(c).replace("|", "\\|") for c in cells) + " |" + + out = ["## Model versions", "", "| # | name | changes |", "|---|---|---|"] + out += [row(f"v{v.number}", v.name, v.changes) for v in self.versions] + out += ["", "## Discrepancies", "", + "| # | observation | hypothesis | test | result | decision | version |", + "|---|---|---|---|---|---|---|"] + for i, d in enumerate(self.discrepancies, start=1): + out.append(row(i, d.observation, d.hypothesis, d.test, d.result, d.decision, + "–" if d.version is None else f"v{d.version}")) + out += ["", "## Acceptance criteria (fixed before tuning)", "", + "| criterion | value | threshold | passed |", "|---|---|---|---|"] + for c in self.criteria: + out.append(row(c.name, f"{c.value:.3g} {c.unit}".strip(), + f"{c.sense} {c.threshold:g} {c.unit}".strip(), + "yes" if c.passed else "**no**")) + return "\n".join(out) + "\n" diff --git a/pyproject.toml b/pyproject.toml index b9d3665..6c278f0 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta" [project] name = "iwt" -version = "0.6.0" +version = "0.7.0" description = "Integrated borehole and surface seismic well tie (Lewis et al., 2016, 78th EAGE) in Python" readme = "README.md" requires-python = ">=3.9" diff --git a/tests/test_workflow.py b/tests/test_workflow.py new file mode 100644 index 0000000..273f36f --- /dev/null +++ b/tests/test_workflow.py @@ -0,0 +1,250 @@ +"""Tests of the WP7 tools: well tie (iwt.welltie), the integrated loop (iwt.workflow), +contamination by multiples (iwt.multiples), linear resolution (iwt.avo) and the +dipping reservoir section (iwt.synthetic).""" + +import numpy as np +import pytest + +from iwt import avo +from iwt import multiples as M +from iwt import raytrace as RT +from iwt import welltie as T +from iwt import workflow as WF +from iwt.model import ElasticModel1D, phase_velocities +from iwt.modeling.convolution import ricker +from iwt.modeling.reflection import exact_pp, exact_pp_slowness +from iwt.qmodel import phase_velocity +from iwt.synthetic import (OWC_TVDSS, SEA_WATER, WellLogs, dipping_reservoir_models, + wisting_like_well) + +DT = 0.002 + + +def _reflectivity(n, seed=1): + rng = np.random.default_rng(seed) + r = np.zeros(n) + k = rng.choice(np.arange(30, n - 30), 25, replace=False) + r[k] = rng.normal(0.0, 0.1, k.size) + return r + + +def _convolve(r, w): + L = (w.size - 1) // 2 + return np.convolve(r, w)[L:L + r.size] + + +# -- well tie --------------------------------------------------------------------------- + + +def test_estimate_wavelet_recovers_a_rotated_wavelet(): + t, w = ricker(25.0, DT) + w = T.rotate_phase(np.pad(w, 20), np.radians(40.0))[20:-20] # padded: no wrap-around + rs = np.array([_reflectivity(500, s) for s in (1, 2)]) + d = np.array([_convolve(r, w) for r in rs]) + tw, we = T.estimate_wavelet(d, rs, DT, t[-1], damping=1e-8) + np.testing.assert_allclose(tw, t, atol=1e-12) + np.testing.assert_allclose(we, w, atol=1e-6 * np.abs(w).max()) + # noise: the estimate degrades gracefully + rng = np.random.default_rng(3) + dn = d + 0.05 * d.std() * rng.standard_normal(d.shape) + _, wn = T.estimate_wavelet(dn, rs, DT, t[-1]) + assert np.linalg.norm(wn - w) / np.linalg.norm(w) < 0.1 + + +def test_tie_metrics_separates_shift_and_phase(): + _, w = ricker(30.0, DT) + s = _convolve(_reflectivity(600), w) + for k, phi in ((0, 0.0), (4, np.radians(35.0)), (-3, np.radians(-60.0))): + o = 2.0 * T.rotate_phase(np.roll(s, k), phi) + m = T.tie_metrics(o, s, DT, max_shift=0.02) + assert m.shift == pytest.approx(k * DT, abs=0.05 * DT) + assert m.phase == pytest.approx(phi, abs=np.radians(1.0)) + assert m.max_correlation == pytest.approx(1.0, abs=1e-3) + if k == 0 and phi == 0.0: + assert m.correlation == pytest.approx(1.0) + assert m.scale == pytest.approx(2.0) + assert m.pep == pytest.approx(1.0) + else: + assert m.correlation < 0.9 + # a pure time shift is not mistaken for a phase rotation (and vice versa) + m = T.tie_metrics(np.roll(s, 5), s, DT) + assert abs(m.phase) < np.radians(1.0) + m = T.tie_metrics(T.rotate_phase(s, np.radians(90.0)), s, DT) + assert abs(m.shift) < 0.05 * DT + + +# -- model from logs -------------------------------------------------------------------- + + +def _logs(): + z = np.arange(470.0, 520.01, 0.5) + rng = np.random.default_rng(0) + vp = 2400.0 * (1.0 + 0.03 * rng.standard_normal(z.size)) + vs = 1000.0 * (1.0 + 0.03 * rng.standard_normal(z.size)) + rho = 2200.0 + 20.0 * rng.standard_normal(z.size) + return WellLogs(z, vp, vs, rho, meta=dict(f_sonic_p=12e3, f_sonic_s=2e3)) + + +def test_model_from_logs_preserves_sonic_time_and_disperses(): + logs = _logs() + z = logs.depth + m = WF.model_from_logs(logs, 418.0, SEA_WATER.vp0, SEA_WATER.rho) + np.testing.assert_allclose(m.interfaces, np.concatenate([[0.0, 418.0], z])) + # layers between samples: the vertical time is the trapezoidal integral of slowness + t = m.vertical_time(z[-1]) - m.vertical_time(z[0]) + assert float(t) == pytest.approx(np.sum(0.5 * (1.0 / logs.vp[1:] + 1.0 / logs.vp[:-1]) * np.diff(z)), rel=1e-12) + # the gap: means of the first 10 m, or the given velocity + top = z < z[0] + 10.0 + assert m.vp0[1] == pytest.approx(np.mean(logs.vp[top])) + assert m.rho[1] == pytest.approx(np.mean(logs.rho[top])) + assert m.qp[0] == np.inf and m.vs0[0] == 0.0 + mg = WF.model_from_logs(logs, 418.0, SEA_WATER.vp0, SEA_WATER.rho, gap_vp=2100.0) + assert mg.vp0[1] == 2100.0 + np.testing.assert_array_equal(mg.vp0[2:], m.vp0[2:]) + # Q: Vp moved from 12 kHz to 30 Hz (slower), Vs too unless disperse_vs=False + q = WF.model_from_logs(logs, 418.0, SEA_WATER.vp0, SEA_WATER.rho, qp=50.0, qs=25.0, + eps=0.1, delta=np.full(z.size, 0.05)) + vp30 = phase_velocity(logs.vp[:2], 50.0, 30.0, 12e3) + assert q.vp0[2] == pytest.approx(2.0 / np.sum(1.0 / vp30)) + assert np.all(q.vp0[2:] < m.vp0[2:]) and np.all(q.vs0[2:] < m.vs0[2:]) + nv = WF.model_from_logs(logs, 418.0, SEA_WATER.vp0, SEA_WATER.rho, qp=50.0, qs=25.0, + disperse_vs=False) + np.testing.assert_allclose(nv.vs0[2:], m.vs0[2:]) + assert np.all(q.qp[1:] == 50.0) and np.all(q.qs[1:] == 25.0) + assert np.all(q.eps[2:] == 0.1) and np.all(q.delta[2:] == 0.05) and q.eps[1] == 0.0 + + +def test_register_and_criteria(): + reg = WF.Register() + m = ElasticModel1D([0.0, 100.0, 200.0], vp0=[1500.0, 2000.0], vs0=[0.0, 800.0], + rho=[1000.0, 2000.0]) + v0 = reg.new_version("logs", m, "as logged") + v1 = reg.new_version("calibrated", m.copy(), "Q | dispersion") + assert (v0.number, v1.number) == (0, 1) and reg.current is v1 + reg.log("late synthetic", "gap velocity", "VSP gap", "shift 1 ms", "adopt", 1) + reg.log("AVA erratic", "multiples", "Born vs full", "0.9 x Born", "explain") + assert reg.check("tie", 0.93, 0.8, ">=").passed + assert not reg.check("drift", 1.5, 1.0, "<=", "ms").passed + with pytest.raises(ValueError): + reg.check("bad", 1.0, 1.0, "==") + assert len(reg.criteria) == 2 + md = reg.to_markdown() + assert "| v1 | calibrated |" in md + assert "| 1 | late synthetic |" in md and "| adopt | v1 |" in md and "| – |" in md + assert "| tie | 0.93 | >= 0.8 | yes |" in md and "| drift | 1.5 ms | <= 1 ms | **no** |" in md + reg.check("|ratio - 1|", 0.5, 0.1) + assert "| \\|ratio - 1\\| | 0.5 |" in reg.to_markdown() # pipes escaped + + +# -- contamination ---------------------------------------------------------------------- + + +def test_relative_rms(): + b = np.sin(np.linspace(0.0, 10.0, 200)) + assert M.relative_rms(1.1 * b, b) == pytest.approx(0.1) + a = b.copy() + a[:100] += 1.0 + assert M.relative_rms(a, b, mask=np.arange(200) >= 100) == pytest.approx(0.0) + a = 1.2 * b + a[5] = np.nan + assert M.relative_rms(a, b) == pytest.approx(0.2, rel=0.02) + + +def test_gathers_switch_physics(): + m = ElasticModel1D([0.0, 200.0, 260.0, 600.0], vp0=[1500.0, 2200.0, 3000.0], + vs0=[0.0, 900.0, 1500.0], rho=[1000.0, 2100.0, 2400.0]) + dt, nt = 0.002, 400 + g = M.gathers(m, [50.0], 5.0, 5.0, dt, nt, ricker(25.0, dt), attenuation=False, f_max=60.0) + assert list(g) == list(M.PHYSICS) + t = np.arange(nt) * dt + seabed = (t > 0.24) & (t < 0.295) # 0.262 s: identical in every variant + base = (t > 0.30) & (t < 0.33) # 0.316 s: below the seabed transmission + for k in M.PHYSICS: + assert M.relative_rms(g[k][0, seabed], g["full"][0, seabed]) < 0.03 + # unit transmission (Born) strengthens the deeper primary by about 1 / (1 - R^2) + r = (2100.0 * 2200.0 - 1000.0 * 1500.0) / (2100.0 * 2200.0 + 1000.0 * 1500.0) + ratio = (np.abs(g["Born primaries"][0, base]).max() + / np.abs(g["primaries with transmission loss"][0, base]).max()) + assert ratio == pytest.approx(1.0 / (1.0 - r * r), rel=0.15) + + +# -- resolution ------------------------------------------------------------------------- + + +def test_linear_resolution(): + J = np.eye(3) + r = avo.linear_resolution(J, [1.0, 2.0, 3.0], 0.5 ** 2 * np.eye(3)) + np.testing.assert_allclose(r.posterior_std ** 2, 1.0 / (4.0 + 1.0 / np.array([1, 4, 9]))) + np.testing.assert_allclose(r.data_weight, 1.0 - r.posterior_std ** 2 / np.array([1, 4, 9])) + np.testing.assert_allclose(r.correlation, np.eye(3), atol=1e-12) + # a parameter the data cannot see keeps its prior; two the data confuse trade off + x = np.linspace(-1.0, 1.0, 20) + J = np.column_stack([x, 1.001 * x, np.zeros_like(x)]) + r = avo.linear_resolution(J, [1.0, 1.0, 1.0], 0.01 ** 2 * np.eye(x.size)) + assert r.data_weight[2] == pytest.approx(0.0) + assert r.correlation[0, 1] < -0.99 + + +# -- kinematic and coefficient helpers -------------------------------------------------- + + +def test_reflection_table_zero_offset_fast_path(): + m = ElasticModel1D([0.0, 300.0, 340.0, 420.0, 600.0], vp0=[1500.0, 2400.0, 2900.0, 2600.0], + vs0=[0.0, 1000.0, 1400.0, 1200.0], rho=[1000.0, 2200.0, 2350.0, 2300.0], + eps=[0.0, 0.15, 0.05, 0.1], delta=[0.0, 0.05, 0.02, 0.04], + qp=[np.inf, 40.0, 80.0, 60.0]) + z = np.array([320.0, 340.0, 400.0, 500.0]) + fast = RT.reflection_table(m, z, [0.0], 6.0, 8.0) + general = RT.reflection_table(m, z, [0.0, 200.0], 6.0, 8.0) + for a, b in zip(fast[2:], general[2:]): + np.testing.assert_allclose(a[:, 0], b[:, 0], rtol=1e-9, atol=1e-15) + + +def test_exact_pp_slowness_matches_angle_form(): + up, lo = (2400.0, 1000.0, 2200.0, 0.15, 0.05), (2900.0, 1500.0, 2350.0, 0.05, 0.02) + th = np.radians(np.linspace(0.0, 40.0, 9)) + V = phase_velocities(th, up[0], up[1], up[3], up[4], 0.0)[0] + np.testing.assert_allclose(exact_pp_slowness(np.sin(th) / V, up, lo), exact_pp(th, up, lo), + rtol=1e-12) + # beyond the critical slowness of the lower medium the coefficient is complex + pc = 1.0 / (lo[0] * np.sqrt(1.0 + 2.0 * lo[3])) + r = exact_pp_slowness([0.5 * pc, 1.05 * pc], up, lo) + assert abs(r[0].imag) < 1e-12 and abs(r[1].imag) > 1e-3 + + +# -- the Fig. 5 analogue ---------------------------------------------------------------- + + +def test_dipping_reservoir_models(): + xs = np.array([0.0, 500.0, 1000.0]) + oil = dipping_reservoir_models(xs, dz=0.5, td=760.0) + brine = dipping_reservoir_models(xs, dz=0.5, td=760.0, oil=False) + shift = xs * np.tan(np.radians(1.2)) + for m, b, s in zip(oil, brine, shift): + zc = np.arange(560.0, 759.0, 0.5) + 0.25 + po, pb = m.at(zc), b.at(zc) + fluid = (zc > 672.0 + s) & (zc < OWC_TVDSS) + # only the oil column differs; the hydrocarbon lowers Vp and density + assert np.all(po["vp0"][fluid] < pb["vp0"][fluid]) + assert np.all(po["rho"][fluid] < pb["rho"][fluid]) + np.testing.assert_array_equal(po["vp0"][~fluid], pb["vp0"][~fluid]) + # the reservoir top dips with the wedge, the contact stays flat + assert m.interfaces[np.searchsorted(m.interfaces, 672.0 + s - 0.1)] == pytest.approx( + 672.0 + s, abs=0.5) + assert np.any(np.isclose(m.interfaces, OWC_TVDSS)) + # units above the wedge are untouched + a, b = oil[0].at([500.0, 590.0]), oil[2].at([500.0, 590.0]) + np.testing.assert_array_equal(a["vp0"], b["vp0"]) + # the oil column (22.5 m at x = 0) thins to 1.5 m at 1 km + assert OWC_TVDSS - 672.0 == pytest.approx(22.5) + assert OWC_TVDSS - (672.0 + shift[-1]) == pytest.approx(1.5, abs=0.1) + + +def test_model_from_wisting_logs_matches_truth_time(): + model, logs = wisting_like_well() + m = WF.model_from_logs(logs, 418.0, SEA_WATER.vp0, SEA_WATER.rho, disperse_vs=False) + z = logs.depth + t_true = model.vertical_time(z[-1]) - model.vertical_time(z[0]) + t_logs = m.vertical_time(z[-1]) - m.vertical_time(z[0]) + # logged at 12 kHz: faster than the 30 Hz truth (sonic drift of order 10 ms) + assert 0.0 < float(t_true - t_logs) < 0.02