diff --git a/literature/figs/demo_14_network_pairs.png b/literature/figs/demo_14_network_pairs.png new file mode 100644 index 0000000..89c0256 Binary files /dev/null and b/literature/figs/demo_14_network_pairs.png differ diff --git a/literature/figs/demo_15_window_envelope.png b/literature/figs/demo_15_window_envelope.png new file mode 100644 index 0000000..b84a961 Binary files /dev/null and b/literature/figs/demo_15_window_envelope.png differ diff --git a/literature/figs/demo_1_methods.png b/literature/figs/demo_1_methods.png index 5661e6e..90ee1a3 100644 Binary files a/literature/figs/demo_1_methods.png and b/literature/figs/demo_1_methods.png differ diff --git a/literature/figs/demo_2_aggregation.png b/literature/figs/demo_2_aggregation.png index 99b2e4c..35bb134 100644 Binary files a/literature/figs/demo_2_aggregation.png and b/literature/figs/demo_2_aggregation.png differ diff --git a/paper/manuscript_marine.qmd b/paper/manuscript_marine.qmd index 1b17b03..40bde71 100644 --- a/paper/manuscript_marine.qmd +++ b/paper/manuscript_marine.qmd @@ -80,10 +80,10 @@ The elevated sensitivity comes at the price a long series of processing choices and at almost every step the analyst makes a choice. Among these are choices of estimators between windowed phase measurements or strething [@Mikesell2015, @Mao2020, @Yuan2021], frequency band and coda window (which together set the sampled depth; @Obermann2013 [@Obermann2016]), - reference window [@Brenguier2014, @Ermert2023, @Okubo24], how much to substack to increase + reference window [@Brenguier2014, @Ermert2023, @Okubo2024], how much to substack to increase coherence among windows (Hadzianou Celine), and how to aggregate and weight the many cross-component and -station-pair measurements that make up a single reported \dvv\ time series (e.g., @hobiger14). +station-pair measurements that make up a single reported \dvv\ time series (e.g., @Hobiger2012). These choices are typically made by habit, justified briefly if at all , and rarely reported in enough detail to reproduce. The community has long flagged *individual* pitfalls (spurious changes from non-stationary noise, @Zhan2013; measurement-error @@ -159,7 +159,6 @@ we do not model these effects using full waveforms. Because the imposed $\dvv(t) a recovered series from it is an artefact of the processing, not of the data. - We utilize seven \dvv\ estimators that were implemented in ``noisepy`` in the `monitoring_methods` module [@Jiang2020]: trace stretching (TS, @lobkis01), windowed cross-correlation (WCC, @poupinet84), dynamic time warping (DTW, @Mikesell2015), @@ -227,62 +226,113 @@ Clock / timing & causal vs acausal branch handling & A clock error fabricates \end{table} ``` -# Measurement: how the choices move the velocity change and its error {#sec:results} +# parameter-dependent \dvv\ and its errors {#sec:results} -The literature agrees on the components of a well-posed \dvv\ measurement --- a +The literature agrees on the components of a well-posed \dvv\ measurement: a stretching-family estimator for robustness at low SNR and large change [@Mikesell2015; @Yuan2021], a coherence-based error model [@Clarke2011; @Weaver2011], a long stable reference [@Wang2017], and -cross-validation against a second estimator [@Obermann2019] --- yet +cross-validation against a second estimator [@Obermann2019]. Yet studies do not always report the same set, and the uncertainty convention is rarely, if ever, quantified (Appendix~\ref{app:survey}). The sections below address each -component in turn and quantify, against a known truth, how far a reasonable -deviation moves the recovered value and its stated error. The deliverable of this +component in turn and quantify, against a known truth, how far a parameter choice +impacts the recovered \dvv\ and its error. The deliverable of this section is a measurement covariance, which every subsequent step in the inference chain (Sections~\ref{sec:depth}--\ref{sec:stress}) consumes. -## Estimator family {#sec:methods-fig} +Each of the parametric component is implemented in a single comprehensive package ``codameter`` +that borrows from ``msnoise`` [@Lecocq2014] and ``noisepy`` (@Jiang2020) to extract the +measurement of \dvv\ from the ambient noise monitoring workflow and generalize it +to any coda wave from repeated source-receiver paths. + +Table~\ref{tab:results-synthesis} previews the RMS error against the known +synthetic ground truth for the best- and worst-case option on each axis +covered in this section, each derived in its own dedicated synthetic (detailed +in the corresponding subsection below); it is a synthesis of *this section's* +per-choice numbers, distinct from Table~\ref{tab:bp-measure}'s one-at-a-time +sweep on a single shared scenario in Section~\ref{sec:multiverse}. + +```{=latex} +\begin{table} +\footnotesize +\caption{Synthesis of Section~\ref{sec:results}: RMS error against the known +synthetic truth for the best- and worst-case option on each axis, each from +its own dedicated synthetic scenario (see the cross-referenced subsection).} +\label{tab:results-synthesis} +\begin{tabularx}{\textwidth}{@{}>{\raggedright\arraybackslash}p{3.2cm} L L L@{}} +\toprule +\textbf{Axis} & \textbf{Best-case RMS} & \textbf{Worst-case RMS} & \textbf{Section} \\ +\midrule +Estimator (family split) & $<\!0.01\,\%$ (TS/WTS, up to 5\,\% true \dvv) & +cycle-skip $>\!0.5\,\%$ past $\sim\!1.4\,\%$ true \dvv\ (MWCS) & +Section~\ref{sec:methods-fig} \\ \hline +Cross-component aggregation & $\sim\!0.03\,\%$ (Approach B, averaged images) & +$\sim\!0.31\,\%$ (Approach A, unweighted) & Section~\ref{sec:aggregation} \\ \hline +Network aggregation (per-pair spread) & $\sim\!0.005\,\%$ (network SE) & +$\sim\!0.05\,\%$ (individual-pair range) & Section~\ref{sec:uncertainty} \\ \hline +Frequency band & $\sim\!0.003$--$0.03\,\%$ (matched to depth) & +$\sim\!0.10\,\%$ (mismatched) & Section~\ref{sec:params} \\ \hline +Coda window & $\sim\!0.01\,\%$ (adapted to band) & +$\sim\!3.8\,\%$ (fixed, wrong band) & Section~\ref{sec:window} \\ \hline +Reference scheme & $\sim\!0.03$--$0.04\,\%$ (fixed stack / joint inversion) & +$\sim\!0.16\,\%$ (60-day moving) & Section~\ref{sec:params} \\ \hline +Stack length & $\sim\!0.018\,\%$ (10-day) & +$\sim\!0.044\,\%$ (1-day, noisy) & Section~\ref{sec:params} \\ +\bottomrule +\end{tabularx} +\end{table} +``` -## TO REVIEW +## Estimator family {#sec:methods-fig} On small, clean \dvv\, all seven estimators agree, demonstrating robustness of the methods -(Fig.~\ref{fig:methods}a). The estimator choice becomes consequential at large -\dvv\ (Fig.~\ref{fig:methods}b), where the methods -split by family. The stretching family (TS, WTS) and WCC match the whole dilated -coda and remain accurate; the phase methods read a wrapped phase, so MWCS -cycle-skips, while the *same* cross-wavelet phase, once unwrapped in 2-D, -recovers the change --- success or failure turning on the unwrapping sub-choice -alone. The warping methods (DTW, WTDTW) track but under-shoot the largest -strains. No estimator is simply "right"; an undocumented estimator choice +(Fig.~\ref{fig:methods}a). Sweeping the same clean recovery out to +$\pm 5\,\%$ (Fig.~\ref{fig:methods}b) shows exactly where and how each family +first departs from the 1:1 line, and the estimator choice becomes consequential +at large, noisy \dvv\ (Fig.~\ref{fig:methods}c), where the effect of the methods +is split according to their phase measurement approaches. The stretching family (TS, WTS) and WCC match the whole dilated +coda and remain accurate for high SNR coda waves; the phase methods (MWCS) read a wrapped phase and may +cycle-skips, while the *same* cross-wavelet phase (WCS), once unwrapped in 2-D, +recovers the change. The warping methods (DTW, WTDTW) track but under-shoot the largest +strains. No estimator is simply "right" and choices in the estimators can alters the \dvv\ measurements for larger strain changes. -The wavelet methods use a -dependency-free Morlet continuous wavelet transform; WCS optionally unwraps the -cross-wavelet phase in two dimensions (along lapse, anchored at $\tau\!\to\!0$, -then along frequency) following @Mao2020. - \begin{figure} \centering \includegraphics[width=\textwidth]{demo_1_methods.png} \caption{Estimator choice across the seven NoisePy methods. (a) Clean, - small \dvv: all agree. (b) Large \dvv\ (a pre-failure landslide signal): MWCS - cycle-skips, 2-D-unwrapped WCS and the stretching family track, the warping - methods under-shoot.} + small \dvv: all agree. (b) The same clean recovery swept over $\pm 5\,\%$ + true \dvv: MWCS cycle-skips almost immediately past $\pm 1$--$2\,\%$, DTW/WTDTW + break past $\pm 3$--$4\,\%$, WCS degrades smoothly, while TS/WTS/WCC track the + 1:1 line throughout. (c) Large, noisy \dvv\ (a pre-failure landslide signal): + MWCS cycle-skips, 2-D-unwrapped WCS and the stretching family track, the + warping methods under-shoot.} \label{fig:methods} \end{figure} -## Aggregating cross-components {#sec:aggregation} +On a clean, noiseless sweep of true \dvv\ from 0 to 5\,\% (Fig.~\ref{fig:methods}b), +the phase-wrapped MWCS estimator is the first to break: its error stays below +$0.1\,\%$ up to a true \dvv\ of $\sim\!1.2\,\%$, then jumps discontinuously past +$0.5\,\%$ error by $\sim\!1.4\,\%$ true \dvv\ --- the cycle-skip. The stretching +family (TS, WTS) stays below $0.01\,\%$ error out to the full 5\,\% tested, and +WCC stays below $0.1\,\%$; the warping methods (DTW, WTDTW) are accurate at small +\dvv\ but develop intermittent, large ($>\!1\,\%$) errors from $\sim\!3\,\%$ +true \dvv\ onward as the warp path becomes ill-conditioned, while WCS degrades +smoothly rather than catastrophically, crossing $1\,\%$ error only beyond +$\sim\!4\,\%$ true \dvv. + +## Aggregating cross-component results {#sec:aggregation} Each three-component seismic station (e.g., Z, N, E) carries 6 cross-componet correlations (ZZ, NN, EE, ZE, ZN, NE), whether they are calculated at the single station or a inter-station pair. Each carry a signature of the changes in velocity, components -may be dominated by Love or Rayleigh waves (REF), but there scattering and non-straight -ray path induce cross-component leakage between modes (REF), and thus it is often +may be dominated by Love or Rayleigh waves [@Lin2008; @Stehly2006], but there scattering and non-straight +ray path induce cross-component leakage between modes [@Hennino2001; @Margerin2019], and thus it is often assumed in practice that coda waves of cross-components with multi-scattering characteristixcs -(e.g., no clearly separated phases) are compoed of "surface waves" with strong S-wave sensitivity. -Combining them together requires parameter choices, such as averaging them directly (REF), -or weighted (e.g., using coherence-based weighting @hobiger14). Combining is another workflow +(e.g., no clearly separated phases) are composed of "surface waves" with strong S-wave sensitivity. +Combining them together requires parameter choices, such as averaging them directly [@Liu2014], +or weighted (e.g., using coherence-based weighting @Hobiger2012, @DePlaen2016). Combining is another workflow choice that can change both the value and the uncertainty (Fig.~\ref{fig:aggregation}). One may peak-pick each component's correlation-coefficient curve $\mathrm{CC}(\varepsilon,t)$ and then average the @@ -292,26 +342,33 @@ $\mathrm{CC}(\varepsilon,t)$ *images* across components first and peak-pick once (Approach B). All three conventions appear in the literature; on the same pair they give visibly different time series, and they propagate uncertainty along incompatible pathways (the ensemble spread of the per-component picks for A, the -width of the averaged correlation peak for B). +width of the averaged correlation peak for B). On our six-component synthetic (three +good, three poor SNR), the RMS error against the known truth is $\sim\!0.31\,\%$ for +unweighted Approach A, $\sim\!0.08\,\%$ for coherence-weighted Approach A (the poor +components suppressed, $\sim\!4\times$ better), and $\sim\!0.03\,\%$ for Approach B +(averaging the images before peak-picking, $\sim\!11\times$ better than the +unweighted mean and $\sim\!3\times$ better than the weighted one). \begin{figure} \centering \includegraphics[width=\textwidth]{demo_2_aggregation.png} - \caption{Cross-component aggregation for one station pair. (a) Unweighted - \dvv-averaging is biased by the poor components, while coherence-weighting and - image-averaging track the truth. (b) The averaged $\mathrm{CC}(\dvv,t)$ image - of Approach B with its peak ridge.} + \caption{Cross-component aggregation for one station pair, both panels on the + same \dvv\ axis. (a) Unweighted \dvv-averaging (Approach A) is biased by the + poor components (grey), while coherence-weighting and image-averaging + (Approach B) track the truth; the shaded band is B's local peak-width + uncertainty. (b) The averaged $\mathrm{CC}(\dvv,t)$ image of Approach B + (dark = high coherence) with its peak ridge (white) tracking the truth + (black, dashed).} \label{fig:aggregation} \end{figure} -## Network aggregation and the reported error bar {#sec:uncertainty} +## Aggregating across station pairs {#sec:uncertainty} -Adding the next layer up --- combining many station pairs into a network series +Adding the next layer up --- combining many station pairs into a unified network series --- exposes the most consequential and least-reported choice of all: how to summarize the uncertainty. Three conventions are common: a coherence-weighted -standard error, an unweighted standard error ($\sigma=\mathrm{std}/\sqrt{N}$), -and the between-pair standard deviation (the scatter, *not* divided by -$\sqrt{N}$). On the same synthetic network the recovered means nearly coincide, +standard error (e.g., @Clarke2011), an unweighted standard error ($\sigma=\mathrm{std}/\sqrt{N}$; e.g., @Brenguier2008), +and the between-pair standard deviation (e.g., @Clements2018). On the same synthetic network the recovered means nearly coincide, but the reported $1\sigma$ spans a factor of $\sim\!\sqrt{N}$ (Fig.~\ref{fig:uncertainty}). A velocity change that is "$3\sigma$ significant" under the tightest convention is "$1\sigma$, not significant" under the most @@ -329,19 +386,57 @@ and the standard-error-versus-standard-deviation convention are all stated. \label{fig:uncertainty} \end{figure} -## Frequency band, coda window, reference and stacking {#sec:params} +Figure~\ref{fig:uncertainty} plots only the *network-aggregate* series, which +hides how much the individual pairs actually disagree. Fig.~\ref{fig:network-pairs} +plots the same nine-pair network's *individual* \dvv(t) curves, styled after a +basin-scale, urban ambient-noise deployment such as the San Gabriel Valley +groundwater network [@Clements2018] --- an illustrative geometry rather than a +literal reproduction of that network's exact station spacing. The individual +pairs range in quality from a coherence-weighted SNR of $\sim\!2.5$ to +$\sim\!11$, and their spread at any given day (median range $\sim\!0.053\,\%$) +is nearly $10\times$ wider than the coherence-weighted network standard error +(median $\sim\!0.005\,\%$) and more than $3\times$ wider than the more +conservative between-pair standard deviation (median $\sim\!0.017\,\%$). A +network-level error bar, however it is computed, describes the precision of +the *mean*, not the *dispersion* of what individual pairs actually report --- +the two are routinely conflated when a single station-pair result is compared +against a published network value. -The remaining choices are no less consequential. The frequency band sets the -sampled depth: in a two-layer medium, banding high recovers a shallow seasonal -signal and banding low recovers a deep trend --- the band chooses *what is -measured* (Fig.~\ref{fig:params}a). The coda window does not transfer across -bands: because the coda fades first at high frequency, a late window that is full -of signal at low frequency is pure noise at high frequency, and reusing it -recovers only noise (Fig.~\ref{fig:params}b). The reference defines what survives: a -moving reference re-baselines continuously and erases slow trends that a fixed -reference or a joint inversion [@Brenguier2014] preserve -(Fig.~\ref{fig:params}c). And the stacking length trades noise against temporal -resolution, rounding off and delaying a coseismic step (Fig.~\ref{fig:params}d). +\begin{figure} + \centering + \includegraphics[width=0.75\textwidth]{demo_14_network_pairs.png} + \caption{Individual station-pair \dvv(t) for the same nine-pair network as + Fig.~\ref{fig:uncertainty}, coloured by pair SNR. The pair-to-pair spread + (grey band) is far wider than any of the three network-level $1\sigma$ + conventions in Fig.~\ref{fig:uncertainty}b.} + \label{fig:network-pairs} +\end{figure} + +## Frequency band, reference and stacking {#sec:params} + +The remaining choices are no less consequential. The **frequency band** sets the +sampled depth: in a two-layer medium, high frequencies recover shallow, often seasonal +signals and low frequencies recover a deeper, maybe more tectonic, signal +(Fig.~\ref{fig:params}a). On our synthetic two-layer groundwater scenario, a band +matched to each layer (1.5--6\,Hz for the shallow, seasonal signal; 0.2--0.8\,Hz +for the deep, drought-trend signal) recovers each with RMS error +$\sim\!0.003$--$0.03\,\%$; using the wrong band for a given depth inflates the +error to $\sim\!0.10\,\%$ for both --- a $\sim\!4$--30$\times$ degradation. + +The **reference** defines what survives: a moving reference re-baselines +continuously and erases slow trends that a fixed reference or a joint inversion +[@Brenguier2014] preserve (Fig.~\ref{fig:params}c). On the volcano synthetic, a +fixed total-stack reference gives RMS $\sim\!0.03\,\%$ and a Brenguier-style +joint inversion $\sim\!0.04\,\%$ (both preserve the pre-eruptive trend), while a +60-day moving reference gives RMS $\sim\!0.16\,\%$ --- roughly 4--5$\times$ worse +--- because it re-baselines away the very trend being measured. + +The **stacking length** trades noise against temporal resolution, rounding off +and delaying a coseismic step (Fig.~\ref{fig:params}d). On the earthquake +synthetic, a 1-day stack (noisy) gives RMS $\sim\!0.044\,\%$, a 10-day stack +$\sim\!0.018\,\%$ (the best of the three), and an over-long 45-day stack +$\sim\!0.031\,\%$ --- worse than the 10-day stack despite averaging down more +noise, because it smears the step itself. \begin{figure} \centering @@ -356,9 +451,81 @@ resolution, rounding off and delaying a coseismic step (Fig.~\ref{fig:params}d). \label{fig:params} \end{figure} +The accepted range for each of these parameters across the published literature +is catalogued directly, study by study, in the survey of +Appendix~\ref{app:survey} (Table~\ref{tab:survey}): the frequency band, coda +window, estimator, and uncertainty treatment actually reported by 103 ambient-noise +\dvv\ studies. + +The **reference** choice also has a construction axis beyond fixed-versus-moving: +how much of the record a "fixed" reference actually stacks. On the volcano +synthetic, a reference built from the whole pre-eruptive record gives RMS +$\sim\!0.03\,\%$, one built from only the earliest 15\,\% of that record (an +"early stack," noisier for having fewer days) gives $\sim\!0.04\,\%$ +($\sim\!1.5\times$ worse), a 60-day moving stack gives $\sim\!0.16\,\%$ +(erasing the trend, as above), and the Brenguier-style joint inversion +[@Brenguier2014] gives $\sim\!0.04\,\%$ while additionally preserving the trend +--- so among the reference-construction choices, whole-period and joint +inversion are comparably good, early-stacking is a modest ($\sim\!1.5\times$) +noise penalty, and only the moving reference actively destroys signal. + +## Coda window and its covariation with frequency band {#sec:window} + +The coda window is not independent of the frequency band --- it deserves its own +treatment because the two covary strongly, and getting this wrong is one of the +larger, more avoidable sources of error in \S\ref{sec:results}. Because intrinsic +and scattering attenuation both grow with frequency, high-frequency coda energy +falls into the noise floor much sooner than low-frequency coda: a coda window +that is well past the direct arrival for a $\sim\!0.3$--$0.8\,$Hz band is, at +$\sim\!3$--$6\,$Hz, sampling almost pure noise (Fig.~\ref{fig:params}b). On our +synthetic, a fixed 20--40\,s window at the high band gives RMS error +$\sim\!3.8\,\%$ (the noise floor, not the signal), while a window hand-adapted +to the band (3--12\,s) recovers the truth at $\sim\!0.01\,\%$ --- nearly a +400$\times$ difference from this one choice alone. This is itself an instance +of a literature-documented tension: band-matched windowing is recommended, but +Table~\ref{tab:survey} shows many surveyed studies instead reuse one fixed +window across bands. + +Hand-adapting the window per band, as above, requires knowing the band in +advance and re-tuning per deployment. A more principled alternative --- used in +our group --- is to track the coda envelope directly and stop the window where +it *flattens* onto the noise floor, rather than pre-specifying a window from a +rule of thumb. We implement this as `coda_window_from_envelope()`: band-pass a +long-term reference stack, smooth its envelope, estimate the noise floor from a +common late-lapse window, and take the window end as the first lapse time past +a short onset where the envelope stays within a factor of that floor for a +sustained interval (not a single noisy dip). Applied blind (without being told +which band it is) to three bands spanning low, mid, and high frequency +(Fig.~\ref{fig:window-envelope}a), the detector recovers windows of +$\sim\!(3,29)\,$s, $\sim\!(3,37)\,$s, and $\sim\!(3,14)\,$s respectively --- +correctly shrinking at the high band, though the low-versus-mid ordering is not +perfectly monotonic on this synthetic (an artifact of how the fixed additive +noise floor interacts with each band's filter, not a claim that the detector is +exact). Recovering \dvv\ with each band's own detected window instead of one +universal fixed (10--30\,s) window (Fig.~\ref{fig:window-envelope}b) gives RMS +$\sim\!0.032\,\%$ vs.\ $\sim\!0.036\,\%$ at the low band (a modest, +$\sim\!1.1\times$ gain), $\sim\!0.017\,\%$ vs.\ $\sim\!0.035\,\%$ at the mid +band ($\sim\!2\times$), and $\sim\!0.019\,\%$ vs.\ $\sim\!1.8\,\%$ at the high +band ($\sim\!93\times$) --- the fixed window is adequate at low frequency and +catastrophic at high frequency, while the envelope-derived window is close to +the best achievable at every band without ever being told what band it is +measuring. + +\begin{figure} + \centering + \includegraphics[width=\textwidth]{demo_15_window_envelope.png} + \caption{Coda window / frequency-band covariation. (a) Smoothed coda + envelopes at three bands (log scale), shaded by each band's + envelope-detected window --- shrinking automatically at higher frequency. + (b) RMS error against the known truth for a single universal fixed window + versus each band's own envelope-derived window: comparable at low + frequency, $\sim\!93\times$ better at high frequency.} + \label{fig:window-envelope} +\end{figure} + ## Choices that create spurious \texorpdfstring{\dvv}{dv/v} {#sec:artifacts} -Some choices do not merely bias the result --- they create spurious signal. A station +Some choices may create spurious signal. A station clock error delays the whole correlation by a lapse-independent shift, producing an apparent \dvv\ that appears with opposite sign on the causal and acausal branches; measuring the two branches separately is the diagnostic @@ -367,6 +534,9 @@ late coda, so a late measurement window reports a coherent spurious *seasonal* \dvv\ many times the real signal while an earlier window stays clean [the waveform-level version of @Zhan2013] (Fig.~\ref{fig:artifacts}b). +Biases from spurious arrivals could be quantified but mostly we should +decontaminate our workflow from these artefacts or not interpret the results. + \begin{figure} \centering \includegraphics[width=\textwidth]{demo_8_artifacts.png} @@ -376,11 +546,16 @@ late coda, so a late measurement window reports a coherent spurious *seasonal* \label{fig:artifacts} \end{figure} -## Causal and acausal branches: a further aggregation choice {#sec:branches} +## Causal and acausal branches {#sec:branches} -The clock-error diagnostic above uses the two branches of the correlation, and -how they are *combined* is itself an under-reported choice. The causal -(positive-lag) and acausal (negative-lag) branches sample opposite-direction +In a symmetric cross-correlation, both sides of the coda (positive or negative lags) +should exhibit the same \dvv\. Due to the directionality of the wavefield +recorded at the two stations, the correlated wavefield in the coda may differ [@Stehly2006]. +While the interpretation of such coda in terms of Earth's structure effect is difficult [@Snieder2002], +the stability of the wavefield excited in the coda is the main requirement for +stable \dvv\ measurements [@Hadziioannou2009]. Given the challenge in interpreting both sides independently, +researchers typically would measure \dvv\ on each lag and then report its average [@Kidiwela2026]. + The causal (positive-lag) and acausal (negative-lag) branches sample opposite-direction paths with different source-side illumination, and in a 3D medium their sensitivity kernels sample partly different volumes, so the two can report genuinely different \dvv\ without either being wrong. Two regimes bound the @@ -388,15 +563,15 @@ choice (Fig.~\ref{fig:branches}). When a change is *localized* to the volume one branch samples, symmetrizing or averaging the branches --- the common default --- dilutes it toward zero, while the branch that carries the change recovers it (Fig.~\ref{fig:branches}a); here preferring the branch of greatest change is -recovering a real quantity. When instead both branches share the *same* change, +an researchers' judgement. +When instead both branches share the *same* change, their difference is measurement noise, and selecting the branch of greatest change over-reports it --- a max-of-two-estimators selection bias that grows as -SNR falls (Fig.~\ref{fig:branches}b). The rule "take the side of the greatest -change" therefore helps only after a structural check: the two branches must -change with the *same sign* (opposite signs are a clock error, -Fig.~\ref{fig:artifacts}a) and be coherent in time. The defensible procedure is -to measure both branches, gate on that consistency, select by the coherence of -the change rather than its amplitude (a criterion independent of the answer, so +SNR falls (Fig.~\ref{fig:branches}b). When both side exhibit the same change +(sign, coherence) but with different magnitude, it is reasonable to use the +\dvv\ of greatest change given the already low sensitivity in coda waves [@Obermann2013]. +An acceptable workflow is to measure both branches, evaluate on that consistency, +select by the coherence of the change rather than its amplitude (a criterion independent of the answer, so it carries no selection bias), and carry the between-branch difference as an explicit term of the measurement covariance $C_d$ (Section~\ref{sec:bayes}) rather than discarding it by averaging. @@ -420,18 +595,31 @@ the error-bar change it induces. # The multiverse of processing choices {#sec:multiverse} -The figures so far each isolate one choice. We now ask the cumulative question -directly. Starting from a single **best-practice baseline** (trace stretching, a +The scenario is a single representative station pair monitoring a shallow volcanic +edifice, in the style of the permanent broadband deployments used at effusive/dome +volcanoes such as Piton de la Fournaise [@Brenguier2008]: a coda band matched to +shallow depths (0.4--1.0\,Hz), a coda window past the direct arrival (10--30\,s), +and daily correlations sampled every 3 days over 2.5 years at a per-day +correlation-coefficient SNR of 7, typical of a continuously operating station. +The synthetic ground truth combines an annual seasonal \dvv\ cycle (as from +near-surface thermoelastic/hydrologic effects), a slow pre-eruptive inflation +ramp, and a sharp co-eruptive velocity drop with partial recovery. This single-pair +scenario isolates the measurement-step choices from the network-aggregation choices +already covered in Sections~\ref{sec:aggregation}--\ref{sec:uncertainty}. + +The previous sections only identified single choices, but the overall research workflow +involves them all. We now estimate the combined effects of these parametic choices. + Starting from a single **best-practice baseline** (trace stretching, a band matched to the target depth, a coda window well past the direct arrival, a 10-day stack, a long stable reference, coherence gating; the cross-cutting rules -of @Brenguier2014 [@Weaver2011; @Clarke2011] as distilled in our survey), we flip -*one* choice at a time to a documented deviation and measure the resulting error against -the known truth (Fig.~\ref{fig:deviations}). The ranking is unambiguous: relative +of @Brenguier2014 [@Weaver2011; @Clarke2011] as distilled in our survey), we change +every parameter at a time to a deviation from best practice documented in the literature + and measure the resulting error against the known truth (Fig.~\ref{fig:deviations}). The ranking is unambiguous: relative to a best-practice RMS error of $\sim\!0.02\,\%$, an unwrapped phase estimator that cycle-skips is catastrophic, a moving reference and a wrapped-phase MWCS each inflate the error roughly an order of magnitude, an over-long stack and a late low-SNR window distort the co-eruptive drop, while the frequency band (which mostly sets -precision here) and the gating are second-order. The Brenguier joint inversion is +precision here) and the gating are second-order. The joint inversion is nearly as good as the fixed reference and preserves the trend. \begin{figure} @@ -445,20 +633,25 @@ nearly as good as the fixed reference and preserves the trend. \label{fig:deviations} \end{figure} -The choices do not act in isolation, so we run the **full factorial**: every -combination of three estimators, two bands, three coda windows, three stack -lengths and two reference schemes --- 108 individually reasonable pipelines --- -on one volcano dataset (Fig.~\ref{fig:multiverse}a). The recovered \dvv\ fans out -across more than an order of magnitude in RMS error, and the fan is widest -exactly at the sharp co-eruptive drop. Attributing the variance of the outcome to -each axis with a first-order (main-effect) sensitivity index +We test the compounding effects of these choices through 108 reasonable workflows selecting + three estimators, two frequency bands, three coda windows, three stacking +lengths and two reference schemes on a synthetic "volcano" dv/v time series +that includes a seasonal oscillation, a slow pre-eruptive inflation ramp, and a +sharp co-eruptive drop with partial exponential recovery (Fig.~\ref{fig:multiverse}a). +The per-day standard deviation across the 108 pipelines varies by a factor of +$\sim\!4$--5 over the time series (from $\sim\!0.4\,\%$ to $\sim\!1.7\,\%$), and is +widest exactly at the sharp co-eruptive drop; the RMS error against the known +truth spans over two orders of magnitude across pipelines ($\sim\!0.02$--$2.8\,\%$). + +Attributing the variance of the outcome to each axis with a first-order (main-effect) sensitivity index (Fig.~\ref{fig:multiverse}b) shows that, for this dataset, the **coda window and the stack length** control the RMS error most, the **stack length** dominates the recovered drop amplitude, and the estimator is third; the first-order indices sum to well under one, so a large part of the spread is *interaction* between choices ---- compounding, not additive. The appropriate object is therefore not a single curve +compounding. The appropriate object is therefore not a single curve but a distribution of \dvv\ over the processing choices. + \begin{figure} \centering \includegraphics[width=\textwidth]{demo_11_multiverse.png} @@ -471,12 +664,12 @@ but a distribution of \dvv\ over the processing choices. \label{fig:multiverse} \end{figure} -Table~\ref{tab:bp-measure} distils this into the measurement-step best practice +Table~\ref{tab:bp-measure} summarizes this into the measurement-step best practice and the documented deviation for each choice, with the consequence the synthetic makes concrete. The baseline and deviation sets are the ones our survey extracts from the literature (the cross-cutting rules of @Snieder2002 [@Clarke2011; -@Weaver2011; @Brenguier2014; @Wang2017; @Obermann2019]) and that -`codameter.deviations` runs. +@Weaver2011; @Brenguier2014; @Wang2017; @Obermann2019]). are implemented in +`codameter.deviations`. ```{=latex} \begin{table} diff --git a/paper/manuscript_marine.tex b/paper/manuscript_marine.tex index 2c0b767..9720ac4 100644 --- a/paper/manuscript_marine.tex +++ b/paper/manuscript_marine.tex @@ -246,7 +246,7 @@ \title{The reproducibility cost of ad-hoc processing choices in ambient-noise seismic velocity-change monitoring} \author{M. A. Denolle} -\date{2026-07-22} +\date{2026-07-24} \begin{document} \maketitle \begin{abstract} @@ -321,11 +321,11 @@ \section{Introduction}\label{sec:intro} frequency band and coda window (which together set the sampled depth; \citeauthor{Obermann2013} \citetext{\citeyear{Obermann2013}; \citealp{Obermann2016}}), reference -window \citep[\citet{Ermert2023}, \citet{Okubo24}]{Brenguier2014}, how +window \citep[\citet{Ermert2023}, \citet{Okubo2024}]{Brenguier2014}, how much to substack to increase coherence among windows (Hadzianou Celine), and how to aggregate and weight the many cross-component and station-pair measurements that make up a single reported \dvv~time -series (e.g., \citet{hobiger14}). These choices are typically made by +series (e.g., \citet{Hobiger2012}). These choices are typically made by habit, justified briefly if at all , and rarely reported in enough detail to reproduce. The community has long flagged \emph{individual} pitfalls (spurious changes from non-stationary noise, \citet{Zhan2013}; @@ -487,70 +487,128 @@ \section{Synthetic framework}\label{sec:methods} \end{tabularx} \end{table} -\section{Measurement: how the choices move the velocity change and its -error}\label{sec:results} +\section{\texorpdfstring{parameter-dependent \dvv~and its +errors}{parameter-dependent ~and its errors}}\label{sec:results} -The literature agrees on the components of a well-posed \dvv~measurement ---- a stretching-family estimator for robustness at low SNR and large -change \citep{Mikesell2015, Yuan2021}, a coherence-based error model -\citep{Clarke2011, Weaver2011}, a long stable reference +The literature agrees on the components of a well-posed +\dvv~measurement: a stretching-family estimator for robustness at low +SNR and large change \citep{Mikesell2015, Yuan2021}, a coherence-based +error model \citep{Clarke2011, Weaver2011}, a long stable reference \citep{Wang2017}, and cross-validation against a second estimator -\citep{Obermann2019} --- yet studies do not always report the same set, -and the uncertainty convention is rarely, if ever, quantified +\citep{Obermann2019}. Yet studies do not always report the same set, and +the uncertainty convention is rarely, if ever, quantified (Appendix\textasciitilde{}\ref{app:survey}). The sections below address each component in turn and quantify, against a known truth, how far a -reasonable deviation moves the recovered value and its stated error. The +parameter choice impacts the recovered \dvv~and its error. The deliverable of this section is a measurement covariance, which every subsequent step in the inference chain (Sections\textasciitilde{}\ref{sec:depth}--\ref{sec:stress}) consumes. -\subsection{Estimator family}\label{sec:methods-fig} +Each of the parametric component is implemented in a single +comprehensive package \texttt{codameter} that borrows from +\texttt{msnoise} \citep{Lecocq2014} and \texttt{noisepy} +(\citet{Jiang2020}) to extract the measurement of \dvv~from the ambient +noise monitoring workflow and generalize it to any coda wave from +repeated source-receiver paths. + +Table\textasciitilde{}\ref{tab:results-synthesis} previews the RMS error +against the known synthetic ground truth for the best- and worst-case +option on each axis covered in this section, each derived in its own +dedicated synthetic (detailed in the corresponding subsection below); it +is a synthesis of \emph{this section's} per-choice numbers, distinct +from Table\textasciitilde{}\ref{tab:bp-measure}'s one-at-a-time sweep on +a single shared scenario in +Section\textasciitilde{}\ref{sec:multiverse}. -\subsection{TO REVIEW}\label{to-review} +\begin{table} +\footnotesize +\caption{Synthesis of Section~\ref{sec:results}: RMS error against the known +synthetic truth for the best- and worst-case option on each axis, each from +its own dedicated synthetic scenario (see the cross-referenced subsection).} +\label{tab:results-synthesis} +\begin{tabularx}{\textwidth}{@{}>{\raggedright\arraybackslash}p{3.2cm} L L L@{}} +\toprule +\textbf{Axis} & \textbf{Best-case RMS} & \textbf{Worst-case RMS} & \textbf{Section} \\ +\midrule +Estimator (family split) & $<\!0.01\,\%$ (TS/WTS, up to 5\,\% true \dvv) & +cycle-skip $>\!0.5\,\%$ past $\sim\!1.4\,\%$ true \dvv\ (MWCS) & +Section~\ref{sec:methods-fig} \\ \hline +Cross-component aggregation & $\sim\!0.03\,\%$ (Approach B, averaged images) & +$\sim\!0.31\,\%$ (Approach A, unweighted) & Section~\ref{sec:aggregation} \\ \hline +Network aggregation (per-pair spread) & $\sim\!0.005\,\%$ (network SE) & +$\sim\!0.05\,\%$ (individual-pair range) & Section~\ref{sec:uncertainty} \\ \hline +Frequency band & $\sim\!0.003$--$0.03\,\%$ (matched to depth) & +$\sim\!0.10\,\%$ (mismatched) & Section~\ref{sec:params} \\ \hline +Coda window & $\sim\!0.01\,\%$ (adapted to band) & +$\sim\!3.8\,\%$ (fixed, wrong band) & Section~\ref{sec:window} \\ \hline +Reference scheme & $\sim\!0.03$--$0.04\,\%$ (fixed stack / joint inversion) & +$\sim\!0.16\,\%$ (60-day moving) & Section~\ref{sec:params} \\ \hline +Stack length & $\sim\!0.018\,\%$ (10-day) & +$\sim\!0.044\,\%$ (1-day, noisy) & Section~\ref{sec:params} \\ +\bottomrule +\end{tabularx} +\end{table} + +\subsection{Estimator family}\label{sec:methods-fig} On small, clean \dvv, all seven estimators agree, demonstrating -robustness of the methods (Fig.\textasciitilde{}\ref{fig:methods}a). The -estimator choice becomes consequential at large -\dvv~(Fig.\textasciitilde{}\ref{fig:methods}b), where the methods split -by family. The stretching family (TS, WTS) and WCC match the whole -dilated coda and remain accurate; the phase methods read a wrapped -phase, so MWCS cycle-skips, while the \emph{same} cross-wavelet phase, -once unwrapped in 2-D, recovers the change --- success or failure -turning on the unwrapping sub-choice alone. The warping methods (DTW, -WTDTW) track but under-shoot the largest strains. No estimator is simply -``right''; an undocumented estimator choice can alters the -\dvv~measurements for larger strain changes. - -The wavelet methods use a dependency-free Morlet continuous wavelet -transform; WCS optionally unwraps the cross-wavelet phase in two -dimensions (along lapse, anchored at \(\tau\!\to\!0\), then along -frequency) following \citet{Mao2020}. +robustness of the methods (Fig.\textasciitilde{}\ref{fig:methods}a). +Sweeping the same clean recovery out to \(\pm 5\,\%\) +(Fig.\textasciitilde{}\ref{fig:methods}b) shows exactly where and how +each family first departs from the 1:1 line, and the estimator choice +becomes consequential at large, noisy +\dvv~(Fig.\textasciitilde{}\ref{fig:methods}c), where the effect of the +methods is split according to their phase measurement approaches. The +stretching family (TS, WTS) and WCC match the whole dilated coda and +remain accurate for high SNR coda waves; the phase methods (MWCS) read a +wrapped phase and may cycle-skips, while the \emph{same} cross-wavelet +phase (WCS), once unwrapped in 2-D, recovers the change. The warping +methods (DTW, WTDTW) track but under-shoot the largest strains. No +estimator is simply ``right'' and choices in the estimators can alters +the \dvv~measurements for larger strain changes. \begin{figure} \centering \includegraphics[width=\textwidth]{demo_1_methods.png} \caption{Estimator choice across the seven NoisePy methods. (a) Clean, - small \dvv: all agree. (b) Large \dvv\ (a pre-failure landslide signal): MWCS - cycle-skips, 2-D-unwrapped WCS and the stretching family track, the warping - methods under-shoot.} + small \dvv: all agree. (b) The same clean recovery swept over $\pm 5\,\%$ + true \dvv: MWCS cycle-skips almost immediately past $\pm 1$--$2\,\%$, DTW/WTDTW + break past $\pm 3$--$4\,\%$, WCS degrades smoothly, while TS/WTS/WCC track the + 1:1 line throughout. (c) Large, noisy \dvv\ (a pre-failure landslide signal): + MWCS cycle-skips, 2-D-unwrapped WCS and the stretching family track, the + warping methods under-shoot.} \label{fig:methods} \end{figure} -\subsection{Aggregating cross-components}\label{sec:aggregation} +On a clean, noiseless sweep of true \dvv~from 0 to 5,\% +(Fig.\textasciitilde{}\ref{fig:methods}b), the phase-wrapped MWCS +estimator is the first to break: its error stays below \(0.1\,\%\) up to +a true \dvv~of \(\sim\!1.2\,\%\), then jumps discontinuously past +\(0.5\,\%\) error by \(\sim\!1.4\,\%\) true \dvv~--- the cycle-skip. The +stretching family (TS, WTS) stays below \(0.01\,\%\) error out to the +full 5,\% tested, and WCC stays below \(0.1\,\%\); the warping methods +(DTW, WTDTW) are accurate at small \dvv~but develop intermittent, large +(\(>\!1\,\%\)) errors from \(\sim\!3\,\%\) true \dvv~onward as the warp +path becomes ill-conditioned, while WCS degrades smoothly rather than +catastrophically, crossing \(1\,\%\) error only beyond \(\sim\!4\,\%\) +true \dvv. + +\subsection{Aggregating cross-component results}\label{sec:aggregation} Each three-component seismic station (e.g., Z, N, E) carries 6 cross-componet correlations (ZZ, NN, EE, ZE, ZN, NE), whether they are calculated at the single station or a inter-station pair. Each carry a signature of the changes in velocity, components may be dominated by -Love or Rayleigh waves (REF), but there scattering and non-straight ray -path induce cross-component leakage between modes (REF), and thus it is -often assumed in practice that coda waves of cross-components with -multi-scattering characteristixcs (e.g., no clearly separated phases) -are compoed of ``surface waves'' with strong S-wave sensitivity. -Combining them together requires parameter choices, such as averaging -them directly (REF), or weighted (e.g., using coherence-based weighting -\citet{hobiger14}). Combining is another workflow choice that can change -both the value and the uncertainty +Love or Rayleigh waves \citep{Lin2008, Stehly2006}, but there scattering +and non-straight ray path induce cross-component leakage between modes +\citep{Hennino2001, Margerin2019}, and thus it is often assumed in +practice that coda waves of cross-components with multi-scattering +characteristixcs (e.g., no clearly separated phases) are composed of +``surface waves'' with strong S-wave sensitivity. Combining them +together requires parameter choices, such as averaging them directly +\citep{Liu2014}, or weighted (e.g., using coherence-based weighting +\citet{Hobiger2012}, \citet{DePlaen2016}). Combining is another workflow +choice that can change both the value and the uncertainty (Fig.\textasciitilde{}\ref{fig:aggregation}). One may peak-pick each component's correlation-coefficient curve \(\mathrm{CC}(\varepsilon,t)\) and then average the per-component \dvv~(Approach A) --- unweighted, a @@ -560,35 +618,45 @@ \subsection{Aggregating cross-components}\label{sec:aggregation} All three conventions appear in the literature; on the same pair they give visibly different time series, and they propagate uncertainty along incompatible pathways (the ensemble spread of the per-component picks -for A, the width of the averaged correlation peak for B). +for A, the width of the averaged correlation peak for B). On our +six-component synthetic (three good, three poor SNR), the RMS error +against the known truth is \(\sim\!0.31\,\%\) for unweighted Approach A, +\(\sim\!0.08\,\%\) for coherence-weighted Approach A (the poor +components suppressed, \(\sim\!4\times\) better), and \(\sim\!0.03\,\%\) +for Approach B (averaging the images before peak-picking, +\(\sim\!11\times\) better than the unweighted mean and \(\sim\!3\times\) +better than the weighted one). \begin{figure} \centering \includegraphics[width=\textwidth]{demo_2_aggregation.png} - \caption{Cross-component aggregation for one station pair. (a) Unweighted - \dvv-averaging is biased by the poor components, while coherence-weighting and - image-averaging track the truth. (b) The averaged $\mathrm{CC}(\dvv,t)$ image - of Approach B with its peak ridge.} + \caption{Cross-component aggregation for one station pair, both panels on the + same \dvv\ axis. (a) Unweighted \dvv-averaging (Approach A) is biased by the + poor components (grey), while coherence-weighting and image-averaging + (Approach B) track the truth; the shaded band is B's local peak-width + uncertainty. (b) The averaged $\mathrm{CC}(\dvv,t)$ image of Approach B + (dark = high coherence) with its peak ridge (white) tracking the truth + (black, dashed).} \label{fig:aggregation} \end{figure} -\subsection{Network aggregation and the reported error -bar}\label{sec:uncertainty} - -Adding the next layer up --- combining many station pairs into a network -series --- exposes the most consequential and least-reported choice of -all: how to summarize the uncertainty. Three conventions are common: a -coherence-weighted standard error, an unweighted standard error -(\(\sigma=\mathrm{std}/\sqrt{N}\)), and the between-pair standard -deviation (the scatter, \emph{not} divided by \(\sqrt{N}\)). On the same -synthetic network the recovered means nearly coincide, but the reported -\(1\sigma\) spans a factor of \(\sim\!\sqrt{N}\) -(Fig.\textasciitilde{}\ref{fig:uncertainty}). A velocity change that is -``\(3\sigma\) significant'' under the tightest convention is -``\(1\sigma\), not significant'' under the most conservative one --- -from identical data. Error bars on published \dvv~are therefore not -comparable across studies unless the aggregation, the weighting, and the -standard-error-versus-standard-deviation convention are all stated. +\subsection{Aggregating across station pairs}\label{sec:uncertainty} + +Adding the next layer up --- combining many station pairs into a unified +network series --- exposes the most consequential and least-reported +choice of all: how to summarize the uncertainty. Three conventions are +common: a coherence-weighted standard error (e.g., \citet{Clarke2011}), +an unweighted standard error (\(\sigma=\mathrm{std}/\sqrt{N}\); e.g., +\citet{Brenguier2008}), and the between-pair standard deviation (e.g., +\citet{Clements2018}). On the same synthetic network the recovered means +nearly coincide, but the reported \(1\sigma\) spans a factor of +\(\sim\!\sqrt{N}\) (Fig.\textasciitilde{}\ref{fig:uncertainty}). A +velocity change that is ``\(3\sigma\) significant'' under the tightest +convention is ``\(1\sigma\), not significant'' under the most +conservative one --- from identical data. Error bars on published +\dvv~are therefore not comparable across studies unless the aggregation, +the weighting, and the standard-error-versus-standard-deviation +convention are all stated. \begin{figure} \centering @@ -600,23 +668,65 @@ \subsection{Network aggregation and the reported error \label{fig:uncertainty} \end{figure} -\subsection{Frequency band, coda window, reference and -stacking}\label{sec:params} - -The remaining choices are no less consequential. The frequency band sets -the sampled depth: in a two-layer medium, banding high recovers a -shallow seasonal signal and banding low recovers a deep trend --- the -band chooses \emph{what is measured} -(Fig.\textasciitilde{}\ref{fig:params}a). The coda window does not -transfer across bands: because the coda fades first at high frequency, a -late window that is full of signal at low frequency is pure noise at -high frequency, and reusing it recovers only noise -(Fig.\textasciitilde{}\ref{fig:params}b). The reference defines what -survives: a moving reference re-baselines continuously and erases slow -trends that a fixed reference or a joint inversion \citep{Brenguier2014} -preserve (Fig.\textasciitilde{}\ref{fig:params}c). And the stacking -length trades noise against temporal resolution, rounding off and -delaying a coseismic step (Fig.\textasciitilde{}\ref{fig:params}d). +Figure\textasciitilde{}\ref{fig:uncertainty} plots only the +\emph{network-aggregate} series, which hides how much the individual +pairs actually disagree. Fig.\textasciitilde{}\ref{fig:network-pairs} +plots the same nine-pair network's \emph{individual} \dvv(t) curves, +styled after a basin-scale, urban ambient-noise deployment such as the +San Gabriel Valley groundwater network \citep{Clements2018} --- an +illustrative geometry rather than a literal reproduction of that +network's exact station spacing. The individual pairs range in quality +from a coherence-weighted SNR of \(\sim\!2.5\) to \(\sim\!11\), and +their spread at any given day (median range \(\sim\!0.053\,\%\)) is +nearly \(10\times\) wider than the coherence-weighted network standard +error (median \(\sim\!0.005\,\%\)) and more than \(3\times\) wider than +the more conservative between-pair standard deviation (median +\(\sim\!0.017\,\%\)). A network-level error bar, however it is computed, +describes the precision of the \emph{mean}, not the \emph{dispersion} of +what individual pairs actually report --- the two are routinely +conflated when a single station-pair result is compared against a +published network value. + +\begin{figure} + \centering + \includegraphics[width=0.75\textwidth]{demo_14_network_pairs.png} + \caption{Individual station-pair \dvv(t) for the same nine-pair network as + Fig.~\ref{fig:uncertainty}, coloured by pair SNR. The pair-to-pair spread + (grey band) is far wider than any of the three network-level $1\sigma$ + conventions in Fig.~\ref{fig:uncertainty}b.} + \label{fig:network-pairs} +\end{figure} + +\subsection{Frequency band, reference and stacking}\label{sec:params} + +The remaining choices are no less consequential. The \textbf{frequency +band} sets the sampled depth: in a two-layer medium, high frequencies +recover shallow, often seasonal signals and low frequencies recover a +deeper, maybe more tectonic, signal +(Fig.\textasciitilde{}\ref{fig:params}a). On our synthetic two-layer +groundwater scenario, a band matched to each layer (1.5--6,Hz for the +shallow, seasonal signal; 0.2--0.8,Hz for the deep, drought-trend +signal) recovers each with RMS error \(\sim\!0.003\)--\(0.03\,\%\); +using the wrong band for a given depth inflates the error to +\(\sim\!0.10\,\%\) for both --- a \(\sim\!4\)--30\(\times\) degradation. + +The \textbf{reference} defines what survives: a moving reference +re-baselines continuously and erases slow trends that a fixed reference +or a joint inversion \citep{Brenguier2014} preserve +(Fig.\textasciitilde{}\ref{fig:params}c). On the volcano synthetic, a +fixed total-stack reference gives RMS \(\sim\!0.03\,\%\) and a +Brenguier-style joint inversion \(\sim\!0.04\,\%\) (both preserve the +pre-eruptive trend), while a 60-day moving reference gives RMS +\(\sim\!0.16\,\%\) --- roughly 4--5\(\times\) worse --- because it +re-baselines away the very trend being measured. + +The \textbf{stacking length} trades noise against temporal resolution, +rounding off and delaying a coseismic step +(Fig.\textasciitilde{}\ref{fig:params}d). On the earthquake synthetic, a +1-day stack (noisy) gives RMS \(\sim\!0.044\,\%\), a 10-day stack +\(\sim\!0.018\,\%\) (the best of the three), and an over-long 45-day +stack \(\sim\!0.031\,\%\) --- worse than the 10-day stack despite +averaging down more noise, because it smears the step itself. \begin{figure} \centering @@ -631,20 +741,105 @@ \subsection{Frequency band, coda window, reference and \label{fig:params} \end{figure} +The accepted range for each of these parameters across the published +literature is catalogued directly, study by study, in the survey of +Appendix\textasciitilde{}\ref{app:survey} +(Table\textasciitilde{}\ref{tab:survey}): the frequency band, coda +window, estimator, and uncertainty treatment actually reported by 103 +ambient-noise \dvv~studies. + +The \textbf{reference} choice also has a construction axis beyond +fixed-versus-moving: how much of the record a ``fixed'' reference +actually stacks. On the volcano synthetic, a reference built from the +whole pre-eruptive record gives RMS \(\sim\!0.03\,\%\), one built from +only the earliest 15,\% of that record (an ``early stack,'' noisier for +having fewer days) gives \(\sim\!0.04\,\%\) (\(\sim\!1.5\times\) worse), +a 60-day moving stack gives \(\sim\!0.16\,\%\) (erasing the trend, as +above), and the Brenguier-style joint inversion \citep{Brenguier2014} +gives \(\sim\!0.04\,\%\) while additionally preserving the trend --- so +among the reference-construction choices, whole-period and joint +inversion are comparably good, early-stacking is a modest +(\(\sim\!1.5\times\)) noise penalty, and only the moving reference +actively destroys signal. + +\subsection{Coda window and its covariation with frequency +band}\label{sec:window} + +The coda window is not independent of the frequency band --- it deserves +its own treatment because the two covary strongly, and getting this +wrong is one of the larger, more avoidable sources of error in +\S\ref{sec:results}. Because intrinsic and scattering attenuation both +grow with frequency, high-frequency coda energy falls into the noise +floor much sooner than low-frequency coda: a coda window that is well +past the direct arrival for a \(\sim\!0.3\)--\(0.8\,\)Hz band is, at +\(\sim\!3\)--\(6\,\)Hz, sampling almost pure noise +(Fig.\textasciitilde{}\ref{fig:params}b). On our synthetic, a fixed +20--40,s window at the high band gives RMS error \(\sim\!3.8\,\%\) (the +noise floor, not the signal), while a window hand-adapted to the band +(3--12,s) recovers the truth at \(\sim\!0.01\,\%\) --- nearly a +400\(\times\) difference from this one choice alone. This is itself an +instance of a literature-documented tension: band-matched windowing is +recommended, but Table\textasciitilde{}\ref{tab:survey} shows many +surveyed studies instead reuse one fixed window across bands. + +Hand-adapting the window per band, as above, requires knowing the band +in advance and re-tuning per deployment. A more principled alternative +--- used in our group --- is to track the coda envelope directly and +stop the window where it \emph{flattens} onto the noise floor, rather +than pre-specifying a window from a rule of thumb. We implement this as +\texttt{coda\_window\_from\_envelope()}: band-pass a long-term reference +stack, smooth its envelope, estimate the noise floor from a common +late-lapse window, and take the window end as the first lapse time past +a short onset where the envelope stays within a factor of that floor for +a sustained interval (not a single noisy dip). Applied blind (without +being told which band it is) to three bands spanning low, mid, and high +frequency (Fig.\textasciitilde{}\ref{fig:window-envelope}a), the +detector recovers windows of \(\sim\!(3,29)\,\)s, \(\sim\!(3,37)\,\)s, +and \(\sim\!(3,14)\,\)s respectively --- correctly shrinking at the high +band, though the low-versus-mid ordering is not perfectly monotonic on +this synthetic (an artifact of how the fixed additive noise floor +interacts with each band's filter, not a claim that the detector is +exact). Recovering \dvv~with each band's own detected window instead of +one universal fixed (10--30,s) window +(Fig.\textasciitilde{}\ref{fig:window-envelope}b) gives RMS +\(\sim\!0.032\,\%\) vs.~\(\sim\!0.036\,\%\) at the low band (a modest, +\(\sim\!1.1\times\) gain), \(\sim\!0.017\,\%\) vs.~\(\sim\!0.035\,\%\) +at the mid band (\(\sim\!2\times\)), and \(\sim\!0.019\,\%\) +vs.~\(\sim\!1.8\,\%\) at the high band (\(\sim\!93\times\)) --- the +fixed window is adequate at low frequency and catastrophic at high +frequency, while the envelope-derived window is close to the best +achievable at every band without ever being told what band it is +measuring. + +\begin{figure} + \centering + \includegraphics[width=\textwidth]{demo_15_window_envelope.png} + \caption{Coda window / frequency-band covariation. (a) Smoothed coda + envelopes at three bands (log scale), shaded by each band's + envelope-detected window --- shrinking automatically at higher frequency. + (b) RMS error against the known truth for a single universal fixed window + versus each band's own envelope-derived window: comparable at low + frequency, $\sim\!93\times$ better at high frequency.} + \label{fig:window-envelope} +\end{figure} + \subsection{\texorpdfstring{Choices that create spurious \texorpdfstring{\dvv}{dv/v}}{Choices that create spurious }}\label{sec:artifacts} -Some choices do not merely bias the result --- they create spurious -signal. A station clock error delays the whole correlation by a -lapse-independent shift, producing an apparent \dvv~that appears with -opposite sign on the causal and acausal branches; measuring the two -branches separately is the diagnostic +Some choices may create spurious signal. A station clock error delays +the whole correlation by a lapse-independent shift, producing an +apparent \dvv~that appears with opposite sign on the causal and acausal +branches; measuring the two branches separately is the diagnostic (Fig.\textasciitilde{}\ref{fig:artifacts}a). Seasonally varying noise sources warp the low-SNR late coda, so a late measurement window reports a coherent spurious \emph{seasonal} \dvv~many times the real signal while an earlier window stays clean \citep[the waveform-level version of][]{Zhan2013} (Fig.\textasciitilde{}\ref{fig:artifacts}b). +Biases from spurious arrivals could be quantified but mostly we should +decontaminate our workflow from these artefacts or not interpret the +results. + \begin{figure} \centering \includegraphics[width=\textwidth]{demo_8_artifacts.png} @@ -654,36 +849,41 @@ \subsection{\texorpdfstring{Choices that create spurious \label{fig:artifacts} \end{figure} -\subsection{Causal and acausal branches: a further aggregation -choice}\label{sec:branches} - -The clock-error diagnostic above uses the two branches of the -correlation, and how they are \emph{combined} is itself an -under-reported choice. The causal (positive-lag) and acausal -(negative-lag) branches sample opposite-direction paths with different -source-side illumination, and in a 3D medium their sensitivity kernels -sample partly different volumes, so the two can report genuinely +\subsection{Causal and acausal branches}\label{sec:branches} + +In a symmetric cross-correlation, both sides of the coda (positive or +negative lags) should exhibit the same \dvv. Due to the directionality +of the wavefield recorded at the two stations, the correlated wavefield +in the coda may differ \citep{Stehly2006}. While the interpretation of +such coda in terms of Earth's structure effect is difficult +\citep{Snieder2002}, the stability of the wavefield excited in the coda +is the main requirement for stable \dvv~measurements +\citep{Hadziioannou2009}. Given the challenge in interpreting both sides +independently, researchers typically would measure \dvv~on each lag and +then report its average \citep{Kidiwela2026}. The causal (positive-lag) +and acausal (negative-lag) branches sample opposite-direction paths with +different source-side illumination, and in a 3D medium their sensitivity +kernels sample partly different volumes, so the two can report genuinely different \dvv~without either being wrong. Two regimes bound the choice (Fig.\textasciitilde{}\ref{fig:branches}). When a change is \emph{localized} to the volume one branch samples, symmetrizing or averaging the branches --- the common default --- dilutes it toward zero, while the branch that carries the change recovers it (Fig.\textasciitilde{}\ref{fig:branches}a); here preferring the branch -of greatest change is recovering a real quantity. When instead both +of greatest change is an researchers' judgement. When instead both branches share the \emph{same} change, their difference is measurement noise, and selecting the branch of greatest change over-reports it --- a max-of-two-estimators selection bias that grows as SNR falls -(Fig.\textasciitilde{}\ref{fig:branches}b). The rule ``take the side of -the greatest change'' therefore helps only after a structural check: the -two branches must change with the \emph{same sign} (opposite signs are a -clock error, Fig.\textasciitilde{}\ref{fig:artifacts}a) and be coherent -in time. The defensible procedure is to measure both branches, gate on -that consistency, select by the coherence of the change rather than its -amplitude (a criterion independent of the answer, so it carries no -selection bias), and carry the between-branch difference as an explicit -term of the measurement covariance \(C_d\) -(Section\textasciitilde{}\ref{sec:bayes}) rather than discarding it by -averaging. +(Fig.\textasciitilde{}\ref{fig:branches}b). When both side exhibit the +same change (sign, coherence) but with different magnitude, it is +reasonable to use the \dvv~of greatest change given the already low +sensitivity in coda waves \citep{Obermann2013}. An acceptable workflow +is to measure both branches, evaluate on that consistency, select by the +coherence of the change rather than its amplitude (a criterion +independent of the answer, so it carries no selection bias), and carry +the between-branch difference as an explicit term of the measurement +covariance \(C_d\) (Section\textasciitilde{}\ref{sec:bayes}) rather than +discarding it by averaging. \begin{figure} \centering @@ -705,24 +905,40 @@ \subsection{Causal and acausal branches: a further aggregation \section{The multiverse of processing choices}\label{sec:multiverse} -The figures so far each isolate one choice. We now ask the cumulative -question directly. Starting from a single \textbf{best-practice -baseline} (trace stretching, a band matched to the target depth, a coda -window well past the direct arrival, a 10-day stack, a long stable -reference, coherence gating; the cross-cutting rules of -\citeauthor{Brenguier2014} +The scenario is a single representative station pair monitoring a +shallow volcanic edifice, in the style of the permanent broadband +deployments used at effusive/dome volcanoes such as Piton de la +Fournaise \citep{Brenguier2008}: a coda band matched to shallow depths +(0.4--1.0,Hz), a coda window past the direct arrival (10--30,s), and +daily correlations sampled every 3 days over 2.5 years at a per-day +correlation-coefficient SNR of 7, typical of a continuously operating +station. The synthetic ground truth combines an annual seasonal +\dvv~cycle (as from near-surface thermoelastic/hydrologic effects), a +slow pre-eruptive inflation ramp, and a sharp co-eruptive velocity drop +with partial recovery. This single-pair scenario isolates the +measurement-step choices from the network-aggregation choices already +covered in +Sections\textasciitilde{}\ref{sec:aggregation}--\ref{sec:uncertainty}. + +The previous sections only identified single choices, but the overall +research workflow involves them all. We now estimate the combined +effects of these parametic choices. Starting from a single +\textbf{best-practice baseline} (trace stretching, a band matched to the +target depth, a coda window well past the direct arrival, a 10-day +stack, a long stable reference, coherence gating; the cross-cutting +rules of \citeauthor{Brenguier2014} \citetext{\citeyear{Brenguier2014}; \citealp{Weaver2011}; \citealp{Clarke2011}} -as distilled in our survey), we flip \emph{one} choice at a time to a -documented deviation and measure the resulting error against the known -truth (Fig.\textasciitilde{}\ref{fig:deviations}). The ranking is -unambiguous: relative to a best-practice RMS error of -\(\sim\!0.02\,\%\), an unwrapped phase estimator that cycle-skips is -catastrophic, a moving reference and a wrapped-phase MWCS each inflate -the error roughly an order of magnitude, an over-long stack and a late -low-SNR window distort the co-eruptive drop, while the frequency band -(which mostly sets precision here) and the gating are second-order. The -Brenguier joint inversion is nearly as good as the fixed reference and -preserves the trend. +as distilled in our survey), we change every parameter at a time to a +deviation from best practice documented in the literature and measure +the resulting error against the known truth +(Fig.\textasciitilde{}\ref{fig:deviations}). The ranking is unambiguous: +relative to a best-practice RMS error of \(\sim\!0.02\,\%\), an +unwrapped phase estimator that cycle-skips is catastrophic, a moving +reference and a wrapped-phase MWCS each inflate the error roughly an +order of magnitude, an over-long stack and a late low-SNR window distort +the co-eruptive drop, while the frequency band (which mostly sets +precision here) and the gating are second-order. The joint inversion is +nearly as good as the fixed reference and preserves the trend. \begin{figure} \centering @@ -735,22 +951,28 @@ \section{The multiverse of processing choices}\label{sec:multiverse} \label{fig:deviations} \end{figure} -The choices do not act in isolation, so we run the \textbf{full -factorial}: every combination of three estimators, two bands, three coda -windows, three stack lengths and two reference schemes --- 108 -individually reasonable pipelines --- on one volcano dataset -(Fig.\textasciitilde{}\ref{fig:multiverse}a). The recovered \dvv~fans -out across more than an order of magnitude in RMS error, and the fan is -widest exactly at the sharp co-eruptive drop. Attributing the variance -of the outcome to each axis with a first-order (main-effect) sensitivity -index (Fig.\textasciitilde{}\ref{fig:multiverse}b) shows that, for this +We test the compounding effects of these choices through 108 reasonable +workflows selecting three estimators, two frequency bands, three coda +windows, three stacking lengths and two reference schemes on a synthetic +``volcano'' dv/v time series that includes a seasonal oscillation, a +slow pre-eruptive inflation ramp, and a sharp co-eruptive drop with +partial exponential recovery +(Fig.\textasciitilde{}\ref{fig:multiverse}a). The per-day standard +deviation across the 108 pipelines varies by a factor of \(\sim\!4\)--5 +over the time series (from \(\sim\!0.4\,\%\) to \(\sim\!1.7\,\%\)), and +is widest exactly at the sharp co-eruptive drop; the RMS error against +the known truth spans over two orders of magnitude across pipelines +(\(\sim\!0.02\)--\(2.8\,\%\)). + +Attributing the variance of the outcome to each axis with a first-order +(main-effect) sensitivity index +(Fig.\textasciitilde{}\ref{fig:multiverse}b) shows that, for this dataset, the \textbf{coda window and the stack length} control the RMS error most, the \textbf{stack length} dominates the recovered drop amplitude, and the estimator is third; the first-order indices sum to well under one, so a large part of the spread is \emph{interaction} -between choices --- compounding, not additive. The appropriate object is -therefore not a single curve but a distribution of \dvv~over the -processing choices. +between choices compounding. The appropriate object is therefore not a +single curve but a distribution of \dvv~over the processing choices. \begin{figure} \centering @@ -764,13 +986,13 @@ \section{The multiverse of processing choices}\label{sec:multiverse} \label{fig:multiverse} \end{figure} -Table\textasciitilde{}\ref{tab:bp-measure} distils this into the +Table\textasciitilde{}\ref{tab:bp-measure} summarizes this into the measurement-step best practice and the documented deviation for each choice, with the consequence the synthetic makes concrete. The baseline and deviation sets are the ones our survey extracts from the literature (the cross-cutting rules of \citeauthor{Snieder2002} -\citetext{\citeyear{Snieder2002}; \citealp{Clarke2011}; \citealp{Weaver2011}; \citealp{Brenguier2014}; \citealp{Wang2017}; \citealp{Obermann2019}}) -and that \texttt{codameter.deviations} runs. +\citetext{\citeyear{Snieder2002}; \citealp{Clarke2011}; \citealp{Weaver2011}; \citealp{Brenguier2014}; \citealp{Wang2017}; \citealp{Obermann2019}}). +are implemented in \texttt{codameter.deviations}. \begin{table} \footnotesize diff --git a/paper/references.bib b/paper/references.bib index 496e7f5..434a7da 100644 --- a/paper/references.bib +++ b/paper/references.bib @@ -132,6 +132,34 @@ @article{Stehly2015 doi = {10.1093/gji/ggv110} } +@article{Lin2008, + author = {Lin, Fan-Chi and Moschetti, Morgan P. and Ritzwoller, Michael H.}, + title = {Surface wave tomography of the western {United States} from ambient seismic noise: {Rayleigh} and {Love} wave phase velocity maps}, + journal = {Geophysical Journal International}, volume = {173}, number = {1}, pages = {281--298}, year = {2008}, + doi = {10.1111/j.1365-246X.2008.03720.x} +} + +@article{Stehly2006, + author = {Stehly, Laurent and Campillo, Michel and Shapiro, Nikolai M.}, + title = {A study of the seismic noise from its long-range correlation properties}, + journal = {Journal of Geophysical Research}, volume = {111}, pages = {B10306}, year = {2006}, + doi = {10.1029/2005JB004237} +} + +@article{Hennino2001, + author = {Hennino, R. and Tr{\'e}gour{\`e}s, N. and Shapiro, N. M. and Margerin, L. and Campillo, M. and van Tiggelen, B. A. and Weaver, R. L.}, + title = {Observation of Equipartition of Seismic Waves}, + journal = {Physical Review Letters}, volume = {86}, pages = {3447--3450}, year = {2001}, + doi = {10.1103/PhysRevLett.86.3447} +} + +@article{Margerin2019, + author = {Margerin, Ludovic and Bajaras, Andres and Campillo, Michel}, + title = {A scalar radiative transfer model including the coupling between surface and body waves}, + journal = {Geophysical Journal International}, volume = {219}, number = {2}, pages = {1092--1108}, year = {2019}, + doi = {10.1093/gji/ggz348} +} + @article{Steegen2016, author = {Steegen, Sara and Tuerlinckx, Francis and Gelman, Andrew and Vanpaemel, Wolf}, title = {Increasing transparency through a multiverse analysis}, @@ -175,6 +203,13 @@ @article{Okubo2024 doi = {10.1029/2023JB028084} } +@article{Kidiwela2026, + author = {Kidiwela, Maleen and Denolle, Marine A. and Wilcock, William S. D. and Feng, K. F.}, + title = {Active protothrusts and fluid highways: Seismic noise reveals hidden subduction dynamics in Cascadia}, + journal = {Science Advances}, volume = {12}, number = {9}, pages = {eaea3684}, year = {2026}, + doi = {10.1126/sciadv.aea3684} +} + % Rock-physics references for the depth/stress framework sections % (density/partial-saturation and the acoustoelastic stress conversion). @article{Gassmann1951, diff --git a/src/codameter/synthetic_demo.py b/src/codameter/synthetic_demo.py index bf307e7..e860f02 100644 --- a/src/codameter/synthetic_demo.py +++ b/src/codameter/synthetic_demo.py @@ -1070,6 +1070,54 @@ def _envelope(x: np.ndarray, fs: float, smooth_s: float = 2.0) -> np.ndarray: return np.convolve(np.abs(x), np.ones(n) / n, mode="same") +def coda_window_from_envelope( + t: np.ndarray, + ref: np.ndarray, + fs: float, + band: tuple[float, float], + *, + t_start: float = 3.0, + floor_factor: float = 1.8, + floor_window: tuple[float, float] = (40.0, 50.0), + persist_s: float = 3.0, + smooth_s: float = 1.5, +) -> tuple[float, float]: + """Pick the coda window automatically: track the envelope, stop where it flattens. + + Band-passes ``ref`` (ideally a long-term, low-noise reference stack, not a + single noisy daily CCF), smooths its envelope, and estimates a noise floor + from the late lapse-time window ``floor_window`` (assumed, for a maxlag long + enough, to already be dominated by noise regardless of band). The window end + is the first lapse time past ``t_start`` where the envelope stays within + ``floor_factor`` of that floor for at least ``persist_s`` seconds (a + sustained flattening, not a single noisy dip) -- the same practice as + tracking the coda envelope by eye and stopping where it visibly flattens, + made automatic and reproducible. + + Frequency-dependent intrinsic attenuation means high-frequency coda energy + falls into the noise floor much sooner than low-frequency coda (see + :func:`make_freqdep_coda`), so this returns a *shorter* window at high + frequency and a *longer* one at low frequency without being told the band + in advance -- it discovers the covariation from the data. + """ + bp = bandpass(ref, fs, *band) + env = _envelope(bp, fs, smooth_s=smooth_s) + pos = t >= 0 + tp, envp = t[pos], env[pos] + fmask = (tp >= floor_window[0]) & (tp <= floor_window[1]) + floor = np.median(envp[fmask]) if fmask.any() else envp[-1] + flat = envp <= floor_factor * floor + persist_n = max(1, int(persist_s * fs)) + sustained = ( + np.convolve(flat.astype(float), np.ones(persist_n), mode="valid") + >= persist_n - 0.5 + ) + tsus = tp[: len(sustained)] + candidates = np.where((tsus >= t_start) & sustained)[0] + t_end = float(tsus[candidates[0]]) if len(candidates) else float(tp[-1]) + return t_start, t_end + + # One colour + line style per NoisePy estimator, grouped by family: # time-domain warp (solid), phase (dashed), wavelet (dash-dot). _MSTYLE = { @@ -1096,7 +1144,15 @@ def fig_methods(seed: int = 11): recs = { m: measure(m, cur, s.ref, s.t, band=band, fs=s.fs, window=win) for m in METHODS } - # (b) a large, smoothly varying change (landslide pre-failure); decimate days + # (b) the same clean recovery, swept over the full +/-5 % range, to show + # exactly where and how each estimator family breaks from the 1:1 line. + trues_wide = np.linspace(-0.05, 0.05, 41) + cur_wide = np.stack([impose_dvv(s.ref, s.t, x) for x in trues_wide]) + recs_wide = { + m: measure(m, cur_wide, s.ref, s.t, band=band, fs=s.fs, window=win) + for m in METHODS + } + # (c) a large, smoothly varying change (landslide pre-failure); decimate days # so the per-day DTW/WTDTW warps stay fast. days = _days(3.0)[::3] truth = landslide_truth(days) @@ -1110,7 +1166,7 @@ def fig_methods(seed: int = 11): kw.update(sub) recL[m] = measure(m, ccfs, s.ref, s.t, **kw) - fig, (axA, axB) = plt.subplots(1, 2, figsize=(6.9, 3.6)) + fig, (axA, axB, axC) = plt.subplots(1, 3, figsize=(9.8, 3.6)) axA.plot([-0.5, 0.5], [-0.5, 0.5], color="0.6", lw=1, ls=":", label="1:1 (truth)") for m in METHODS: col, ls = _MSTYLE[m] @@ -1130,18 +1186,36 @@ def fig_methods(seed: int = 11): title="(a) clean, small dv/v", ) axA.legend(loc="upper left", fontsize=8, ncol=2) - axB.plot(_yrs(days), truth * PCT, color=C["truth"], lw=2.6, label="truth") + axB.plot([-5, 5], [-5, 5], color="0.6", lw=1, ls=":", label="1:1 (truth)") for m in METHODS: col, ls = _MSTYLE[m] axB.plot( - _yrs(days), recL[m] * PCT, ls=ls, color=col, lw=1.2, alpha=0.9, label=m + trues_wide * PCT, + recs_wide[m] * PCT, + ls=ls, + lw=1.3, + color=col, + alpha=0.9, + label=m, ) axB.set( + xlabel="true dv/v (%)", + ylabel="recovered dv/v (%)", + title="(b) clean, $\\pm 5\\,\\%$ sweep", + ) + axB.legend(loc="upper left", fontsize=7.5, ncol=2) + axC.plot(_yrs(days), truth * PCT, color=C["truth"], lw=2.6, label="truth") + for m in METHODS: + col, ls = _MSTYLE[m] + axC.plot( + _yrs(days), recL[m] * PCT, ls=ls, color=col, lw=1.2, alpha=0.9, label=m + ) + axC.set( xlabel="time (years)", ylabel="dv/v (%)", - title="(b) large dv/v — MWCS cycle-skips", + title="(c) large dv/v — MWCS cycle-skips", ) - axB.legend(loc="lower left", fontsize=8, ncol=2) + axC.legend(loc="lower left", fontsize=8, ncol=2) fig.tight_layout() return fig @@ -1175,13 +1249,22 @@ def fig_aggregation(seed: int = 88): A_unw = dvv_c.mean(axis=0) A_wt = (cc_c * dvv_c).sum(axis=0) / (cc_c.sum(axis=0) + 1e-12) # Approach B — average the CC images, then peak-pick once. Its uncertainty is - # the *width of the averaged CC peak* (treating CC as a likelihood over eps), - # a different statistical object from A's ensemble spread. + # the *local width of the averaged CC peak* (a half-max/FWHM-style width + # around the peak, not a moment over the full search range -- the latter is + # dominated by the width of the epsilon search window itself, not by how + # sharp the peak actually is, and barely varies day to day). mean_img = images.mean(axis=0) B, _ = peak_dvv(es, mean_img) - w = np.clip(mean_img, 0, None) - mu = (w * es).sum(1) / w.sum(1) - sig_B = np.sqrt((w * (es - mu[:, None]) ** 2).sum(1) / w.sum(1)) + peak_val = mean_img.max(axis=1) + w_local = np.clip(mean_img - (peak_val / 2.0)[:, None], 0, None) + mu_local = (w_local * es).sum(1) / w_local.sum(1) + sig_B = np.sqrt((w_local * (es - mu_local[:, None]) ** 2).sum(1) / w_local.sum(1)) + + # Shared y-limits so (a) and (b) are directly comparable, sized to fit + # Approach A's real (unclipped) excursions -- the per-component grey lines + # go further still (poor components swing to +/-6 %) but are background + # context, not the point, so they are allowed to clip at the edges. + ylim = (-1.4, 1.4) fig, (axA, axB) = plt.subplots(1, 2, figsize=(6.9, 3.6)) for d in dvv_c: @@ -1222,21 +1305,28 @@ def fig_aggregation(seed: int = 88): axA.set( xlabel="time (years)", ylabel="dv/v (%)", - ylim=(-0.4, 0.4), - title="(a) 3 defensible recipes, 3 answers", + ylim=ylim, + title="(a) Component aggregation: three recipes", ) leg = axA.legend(loc="lower left", fontsize=7.5, frameon=True) leg.get_frame().set(facecolor="white", alpha=0.9, edgecolor="0.7") extent = [_yrs(days)[0], _yrs(days)[-1], es[0] * PCT, es[-1] * PCT] - axB.imshow( - mean_img.T, aspect="auto", origin="lower", extent=extent, cmap="magma", vmin=0 + im = axB.imshow( + mean_img.T, + aspect="auto", + origin="lower", + extent=extent, + cmap="magma_r", # reversed: dark = high CC, matching the paper's convention + vmin=0, ) - axB.plot(_yrs(days), B * PCT, color="white", lw=0.8, label="peak of mean CC (B)") - axB.plot(_yrs(days), truth * PCT, color="cyan", lw=1.0, ls="--", label="truth") + cbar = fig.colorbar(im, ax=axB, pad=0.02) + cbar.set_label("coherence CC (dark = high)") + axB.plot(_yrs(days), B * PCT, color="white", lw=1.2, label="peak of mean CC (B)") + axB.plot(_yrs(days), truth * PCT, color="black", lw=1.0, ls="--", label="truth") axB.set( xlabel="time (years)", ylabel="dv/v candidate (%)", - ylim=(-0.35, 0.3), + ylim=ylim, title="(b) averaged CC(dv/v, t) image (B)", ) axB.legend(loc="upper right", fontsize=8) @@ -1302,6 +1392,8 @@ def network_dvv( return { "truth": truth, "es": es, + "pair_dvv": pair_dvv, # [n_pairs, ndays] -- the individual pair curves + "pair_snr": pair_snr, "weighted_se": {"dvv": W, "sigma": W_se}, "unweighted_se": {"dvv": U, "sigma": sd / np.sqrt(n_pairs)}, "unweighted_sd": {"dvv": U, "sigma": sd}, @@ -1355,6 +1447,59 @@ def fig_uncertainty(seed: int = 123): return fig +def fig_network_pairs(seed: int = 123): + """The per-pair spread that Fig. fig:uncertainty's network-level view hides. + + A basin-scale, urban ambient-noise deployment (in the style of the San + Gabriel Valley groundwater network of Clements & Denolle 2018 -- an + illustrative, not a literal, reproduction of that network's exact station + geometry) has station pairs of heterogeneous quality: some pairs sit on + thick, well-coupled sediment with high SNR, others are noisier. Plotting + the individual per-pair dv/v(t) curves (not just the network-aggregate + mean and its error bars, as in Fig. fig:uncertainty) shows that the true + pair-to-pair spread is wider than any of the three network conventions' + error bars communicate on their own. + """ + import matplotlib.pyplot as plt + + days = _days(3.0) + truth = _seasonal(days, 0.0012, 60) - 0.0008 * (days >= int(1.5 * YEAR_D)) + R = network_dvv(truth, seed=seed) + pair_dvv, pair_snr = R["pair_dvv"], R["pair_snr"] + yr = _yrs(days) + + fig, ax = plt.subplots(figsize=(6.0, 4.0)) + order = np.argsort(-pair_snr) # best SNR first, for a readable legend/colour ramp + cmap = plt.get_cmap("viridis") + for rank, p in enumerate(order): + ax.plot( + yr, + pair_dvv[p] * PCT, + color=cmap(rank / max(1, len(order) - 1)), + lw=0.9, + alpha=0.85, + ) + ax.plot(yr, truth * PCT, color=C["truth"], lw=2.4, label="network truth") + lo = np.nanmin(pair_dvv, axis=0) * PCT + hi = np.nanmax(pair_dvv, axis=0) * PCT + ax.fill_between( + yr, lo, hi, color="0.5", alpha=0.15, lw=0, label="individual-pair range" + ) + sm = plt.cm.ScalarMappable( + cmap=cmap, norm=plt.Normalize(vmin=pair_snr.min(), vmax=pair_snr.max()) + ) + cbar = fig.colorbar(sm, ax=ax, pad=0.02) + cbar.set_label("pair SNR") + ax.set( + xlabel="time (years)", + ylabel="dv/v (%)", + title="Individual station-pair dv/v -- wider than the network error bar", + ) + ax.legend(loc="lower left", fontsize=8.5) + fig.tight_layout() + return fig + + def fig_window_band(seed: int = 66): """Coda window must scale with frequency band: a fixed late window is full of signal at low frequency but pure noise at high frequency.""" @@ -1428,6 +1573,89 @@ def fig_window_band(seed: int = 66): return fig +def fig_window_envelope(seed: int = 71): + """The coda window / frequency-band covariation, and an automatic fix. + + Three bands, each with its own :func:`coda_window_from_envelope` detection: + (a) the smoothed envelopes with the detected window marked per band, + showing the window shrinking automatically as the band moves to higher + frequency. (b) dv/v RMS against the known truth for a single universal + fixed window (chosen for the low band) versus each band's own + envelope-derived window -- the fixed window is fine at low frequency but + catastrophic at high frequency, where it samples almost pure noise; the + envelope-derived window recovers the truth at every band without being + told the band in advance. + """ + import matplotlib.pyplot as plt + + s = Synth() + fs = s.fs + tf, cf = make_freqdep_coda(fs=fs, seed=2) + bands = [(0.3, 0.8), (1.0, 2.5), (3.0, 6.0)] + band_labels = ["low 0.3–0.8 Hz", "mid 1.0–2.5 Hz", "high 3.0–6.0 Hz"] + band_cols = [C["alt"], C["landslide"], C["groundwater"]] + fixed_window = (10.0, 30.0) # a single universal window, ignoring the band + + ref_stack = daily_ccfs( + tf, [cf], [np.zeros(60)], fs=fs, snr=8.0, gen_band=(0.2, 8.0), seed=5 + ).mean(axis=0) + days = _days(1.5)[::3] + truth = _seasonal(days, 0.0015, 60) + + windows, rms_fixed, rms_adapt = [], [], [] + fig, (axA, axB) = plt.subplots(1, 2, figsize=(6.9, 3.6)) + m = tf >= 0 + for band, lab, col in zip(bands, band_labels, band_cols, strict=True): + t1, t2 = coda_window_from_envelope(tf, ref_stack, fs, band) + windows.append((t1, t2)) + env = _envelope(bandpass(ref_stack, fs, *band), fs, smooth_s=1.5) + norm = env[m].max() + axA.semilogy(tf[m], env[m] / norm, color=col, lw=1.4, label=lab) + axA.axvspan(t1, t2, color=col, alpha=0.12, lw=0) + + ccfs = daily_ccfs( + tf, [cf], [truth], fs=fs, snr=8.0, gen_band=(0.2, 8.0), seed=seed + ) + rf, _ = measure_stretching(ccfs, cf, tf, band=band, fs=fs, window=fixed_window) + ra, _ = measure_stretching(ccfs, cf, tf, band=band, fs=fs, window=(t1, t2)) + v = np.isfinite(rf) + rms_fixed.append(float(np.sqrt(np.mean((rf[v] - truth[v]) ** 2))) * PCT) + v = np.isfinite(ra) + rms_adapt.append(float(np.sqrt(np.mean((ra[v] - truth[v]) ** 2))) * PCT) + + axA.set( + xlabel="lapse time (s)", + ylabel="coda envelope (norm.)", + ylim=(1e-3, 2), + title="(a) window shrinks with frequency", + ) + axA.legend(loc="upper right", fontsize=8) + + x = np.arange(len(bands)) + w = 0.35 + axB.bar( + x - w / 2, rms_fixed, width=w, color=C["bad"], label="universal fixed window" + ) + axB.bar( + x + w / 2, + rms_adapt, + width=w, + color=C["groundwater"], + label="envelope-adaptive window", + ) + axB.set_yscale("log") + axB.set( + xticks=x, + xticklabels=band_labels, + ylabel="RMS error vs truth (dv/v, %, log)", + title="(b) fixed window fails at high band", + ) + axB.tick_params(axis="x", labelsize=7.5, rotation=15) + axB.legend(loc="upper left", fontsize=8) + fig.tight_layout() + return fig + + def fig_stacking(seed: int = 22): """Earthquake: stack length trades noise against coseismic-step sharpness.""" import matplotlib.pyplot as plt @@ -1837,8 +2065,10 @@ def fig_branch_asymmetry(seed: int = 131): "demo_1_methods": fig_methods, "demo_2_aggregation": fig_aggregation, "demo_3_uncertainty": fig_uncertainty, + "demo_14_network_pairs": fig_network_pairs, "demo_4_frequency_depth": fig_frequency_depth, "demo_5_window_band": fig_window_band, + "demo_15_window_envelope": fig_window_envelope, "demo_6_stacking": fig_stacking, "demo_7_reference": fig_reference, "demo_8_artifacts": fig_artifacts,