Stack

Moving stack computation.

Reads per-day or per-window CCF NetCDF files written by the CC step and produces moving-window stacked CCFs for every (pair, component, mov_stack) combination.

The stacking window is defined by |stack.mov_stack| as a list of (window, step) tuples. For example ('7D', '1D') produces a rolling 7-day mean, sampled daily. Each output is written as a NetCDF file under the lineage output path:

OUTPUT/.../stack_N/_output/<mov_stack>/<component>/<sta1>_<sta2>.nc

When |stack.wienerfilt| is Y, a Wiener filter is applied to the CCF time series before the rolling mean, attenuating isolated spikes and filling short gaps.

This step is upstream of refstack: once Done, propagate_downstream creates both a REF sentinel job (triggering s04_stack_refstack) and direct MWCS/stretching/wavelet T jobs for the new days.

To run this step:

$ msnoise cc stack

Parallel processing:

$ msnoise -t 4 cc stack

Configuration Parameters

• stack.mov_stack : A list of two parameters: the time to “roll” over (default 1 day) and the granularity (step) of the resulting stacked CCFs (default 1 day) to stack for the Moving-window stacks. This can be a list of tuples, e.g. ((‘1d’,’1d’),(‘2d’,’1d’)) corresponds to the MSNoise 1.6 “1,2” before. Time deltas can be anything pandas can interpret (“d”, “min”, “sec”, etc). (default=((‘1D’,’1D’)))

• stack.wienerfilt : Apply wiener filter before stacking? Y[N] (default=N)

• stack.wiener_mlen : Smoothing along date axis (time delta) (default=24h)

• stack.wiener_nlen : Smoothing along lag time axis (time delta) (default=0.5s)

• cc.keep_all : Keep all cross-corr (length: corr_duration) (default=Y)

• cc.corr_duration : Data windows to correlate (in seconds) (default=1800.0)

• cc.cc_sampling_rate : Sampling Rate for the CrossCorrelation (in Hz) (default=20.0)

• cc.components_to_compute : List (comma separated) of components to compute between two different stations (default=ZZ)

• cc.components_to_compute_single_station : List (comma separated) of components within a single station. ZZ would be the autocorrelation of Z component, while ZE or ZN are the cross-components. Defaults to [], no single-station computations are done. (default=)

• global.hpc : Is MSNoise going to run on an HPC? (default=N)

Stacking Methods

MSNoise supports three stacking methods, selected via the stack_method configuration key (valid for both the moving stack and the reference stack). All three reduce an array of \(N\) daily CCF traces \(\{d_j(t)\}_{j=1}^{N}\) to a single representative trace.

Linear stack

The arithmetic mean:

\[s(t) = \frac{1}{N} \sum_{j=1}^{N} d_j(t)\]

Incoherent noise cancels at a rate of \(1/\sqrt{N}\). This is the fastest and most transparent option. It provides no suppression of high-amplitude transients (earthquakes, instrument glitches) that survive pre-processing.

Set stack_method = linear in the stack or refstack configset.

Phase-weighted stack (pws)

Introduced by Schimmel and Paulssen[1]. Each sample is weighted by the instantaneous phase coherence \(c(t)\) of the analytic signal across all traces:

\[c(t) = \frac{1}{N} \left| \sum_{j=1}^{N} e^{i\,\phi_j(t)} \right|, \qquad c(t) \in [0,\,1]\]

where \(\phi_j(t) = \arg\!\bigl(d_j(t) + i\,\mathcal{H}\{d_j\}(t)\bigr)\) is the instantaneous phase of trace \(j\) obtained via the Hilbert transform \(\mathcal{H}\). The coherence is smoothed with a boxcar window of pws_timegate seconds and then raised to the power \(v\) = pws_power:

\[s(t) = \frac{1}{N} \sum_{j=1}^{N} d_j(t) \cdot c(t)^v\]

Incoherent transients produce random instantaneous phases across traces, so \(c(t) \approx 0\) at those times; truly coherent arrivals yield \(c(t) \approx 1\) and are preserved.

Practical notes:

• pws_timegate (default 10 s) controls the width of the boxcar smoothing applied to \(c(t)\) before the power is taken. Shorter gates resolve rapid phase changes; longer gates produce a smoother weight and are more robust with few traces.

• pws_power (default 2): higher values increase selectivity at the cost of amplitude fidelity on weak coherent arrivals.

Time-frequency phase-weighted stack (tfpws)

The TF extension of PWS introduced by Schimmel and Gallart[2]. Phase coherence is evaluated independently at each scale of a continuous wavelet transform (CWT), yielding a time–frequency coherence map \(c(a, t)\) that is both scale- and time-dependent.

A complex Morlet wavelet is used:

\[\psi_{s}(t) = \frac{1}{\pi^{1/4}\sqrt{s}}\, e^{i\omega_0 t/s}\, e^{-t^2/(2s^2)}, \qquad \omega_0 = 5\]

Scales \(\{s_k\}\) are log-spaced between freqmin and freqmax (inherited from the parent filter configset), with the scale–frequency relationship \(s_k = \omega_0 f_s / (2\pi f_k)\).

The phase coherence at each (scale, lag) point is:

\[c(s_k, t) = \frac{1}{N} \left| \sum_{j=1}^{N} e^{i\,\arg\mathcal{W}_j(s_k,\,t)} \right|\]

where \(\mathcal{W}_j(s_k, t)\) is the CWT coefficient of trace \(j\). Averaging over the \(A\) = tfpws_nscales scales and raising to the power \(v\) yields the per-lag weight:

\[w(t) = \left[ \frac{1}{A} \sum_{k=1}^{A} c(s_k,\,t) \right]^{v}\]

The final stack is:

\[s(t) = \frac{1}{N} \sum_{j=1}^{N} d_j(t) \cdot w(t)\]

Because coherence is assessed independently at each frequency, tf-PWS is more sensitive to narrow-band coherent arrivals than time-domain PWS. Romero and Schimmel[3] demonstrate this advantage for noise autocorrelations targeting shallow P-wave reflections (3–18 Hz), where amplitude-based transients from Pyrenean seismicity would otherwise corrupt the CCGN result.

Practical notes:

• tfpws_nscales (default 20): more scales improve frequency resolution but increase memory and compute time proportionally. Memory scales as \(O(N \times A \times T)\) complex128.

• pws_power is shared with the pws method (default 2).

• freqmin / freqmax are taken automatically from the parent filter configset — no extra configuration is needed.

Comparison

Method	Amplitude sensitivity	Frequency selectivity	Compute cost	Best for
linear	High	None	Lowest	Dense arrays, clean data
pws	Low	None (time-domain)	Low	General use, impulsive noise
tfpws	Low	High	Moderate	Narrow-band targets, autocorrelations

[1]

M. Schimmel and H. Paulssen. Noise reduction and detection of weak, coherent signals through phase-weighted stacks. Geophysical Journal International, 130(2):497–505, 1997. doi:10.1111/j.1365-246X.1997.tb05664.x.

[2]

M. Schimmel and J. Gallart. Frequency-dependent phase coherence for noise suppression in seismic array data. Journal of Geophysical Research, 112:B04303, 2007. doi:10.1029/2006JB004680.

[3]

Paula Romero and Martin Schimmel. Mapping the basement of the Ebro basin in Spain with seismic ambient noise autocorrelations. Journal of Geophysical Research: Solid Earth, 123:5052–5067, 2018. doi:10.1029/2018JB015498.

See also

Reading these results in Python — use MSNoiseResult:

from msnoise.results import MSNoiseResult
from msnoise.core.db import connect
db = connect()
r = MSNoiseResult.from_ids(db, ...)  # include the steps you need
# then call r.get_ccf(...) or r.get_ref(...)

See Reading outputs with MSNoiseResult for the full guide and all available methods.

See Reading outputs with MSNoiseResult for the full guide and all available methods.