Ambient Noise Seismology: From Theory to Practice¶

A Comprehensive Course in Seismic Interferometry and Ambient Noise Processing¶

Course Structure:

Module Topic Type
1 Introduction to Seismic Ambient Noise Theory
2 Seismic Data Handling with ObsPy Practice
3 The Cross-Correlation Theorem Theory + Practice
4 Green's Function Retrieval from Ambient Noise Theory
5 Ambient Noise Preprocessing Theory + Practice
6 Cross-Correlation Computation Practice
7 Stacking Methods Theory + Practice
8 Optimal Processing of Noise Correlations Theory + Practice
9 Surface Wave Dispersion Analysis Theory + Practice
10 Ambient Noise Tomography: Inversion Theory + Practice
11 Seismic Interferometry Principles Theory
12 Seismic Velocity Changes (dv/v) Monitoring Theory + Practice
13 HVSR / MHVSR Method Theory + Practice
14 NoisePy: Large-Scale Ambient Noise Processing Practice
15 Applications and Case Studies Theory + Practice
16 Dense Arrays and Advanced Applications Theory + Practice

Key References:

  • Aki, K. (1957). Space and time spectra of stationary stochastic waves, with special reference to microtremors. Bull. Earthquake Res. Inst., 35, 415-457.
  • Bensen, G.D., Ritzwoller, M.H., Barmin, M.P., et al. (2007). Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements. Geophys. J. Int., 169, 1239-1260.
  • Bensen, G.D., Ritzwoller, M.H. & Shapiro, N.M. (2008). Broadband ambient noise surface wave tomography across the United States. J. Geophys. Res., 113, B05306.
  • Campillo, M. & Paul, A. (2003). Long-range correlations in the diffuse seismic coda. Science, 299, 547-549.
  • Campillo, M. & Roux, P. (2015). Crust and Lithospheric Structure - Seismic Imaging and Monitoring with Ambient Noise Correlations. Treatise on Geophysics, 2nd ed., Vol. 1, pp. 391-417.
  • Cabrera-Pérez, I. et al. (2023). Geothermal and structural features of La Palma island imaged by ambient noise tomography. Sci. Rep., 13, 12892.
  • Chmiel, M. et al. (2019). Ambient noise multimode Rayleigh and Love wave tomography to determine the shear velocity structure above the Groningen gas field. Geophys. J. Int., 218, 1781-1795.
  • Claerbout, J.F. (1968). Synthesis of a layered medium from its acoustic transmission response. Geophysics, 33, 264-269.
  • Cox, B.R. et al. (2020). A statistical representation and frequency-domain window-rejection algorithm for single-station HVSR measurements. Geophys. J. Int., 221(3), 2170-2183.
  • Fichtner, A., Bowden, D. & Ermert, L. (2020). Optimal processing for seismic noise correlations. Geophys. J. Int., 223, 1548-1564.
  • Fichtner, A. et al. (2017). Generalised interferometry I: Theory for inter-station correlations. Geophys. J. Int., 208, 603-638.
  • Lecocq, T. et al. (2014). MSNoise, a Python Package for Monitoring Seismic Velocity Changes Using Ambient Seismic Noise. Seismol. Res. Lett., 85(3), 715‑726.
  • Jiang, C. & Denolle, M.A. (2020). NoisePy: A New High-Performance Python Tool for Ambient-Noise Seismology. Seismol. Res. Lett., 91(3), 1853-1866.
  • Brenguier, F., Rivet, D., Obermann, A., et al. (2016). 4-D noise-based seismology at volcanoes. J. Volcanol. Geotherm. Res., 321, 182-195.
  • Lobkis, O.I. & Weaver, R.L. (2001). On the emergence of the Green's function in the correlations of a diffuse field. J. Acoust. Soc. Am., 110(6), 3011-3017.
  • Molnar, S. et al. (2022). A review of the microtremor horizontal-to-vertical spectral ratio (MHVSR) method. J. Seismol., 26, 653-685.
  • Rawlinson, N. et al. (2010). Seismic tomography: A window into deep Earth. Phys. Earth Planet. Inter., 178, 101-135.
  • Ryberg, T. et al. (2022). Ambient seismic noise analysis of LARGE-N data for mineral exploration in the Central Erzgebirge, Germany. Solid Earth, 13, 519-533.
  • Schimmel, M. & Paulssen, H. (1997). Noise reduction and detection of weak, coherent signals through phase-weighted stacks. Geophys. J. Int., 130, 497-505.
  • Sens-Schoenfelder, C. & Wegler, U. (2006). Passive image interferometry and seasonal variations of seismic velocities at Merapi Volcano, Indonesia. Geophys. Res. Lett., 33, L21302.
  • Shapiro, N.M. & Campillo, M. (2004). Emergence of broadband Rayleigh waves from correlations of the ambient seismic noise. Geophys. Res. Lett., 31, L07614.
  • Snieder, R. (2004). Extracting the Green's function from the correlation of coda waves. Phys. Rev. E, 69, 046610.
  • Stehly, L. et al. (2024). Dynamic of seismic noise sources in the Mediterranean Sea: implication for monitoring using noise correlations. C. R. Géoscience, doi:10.5802/crgeos.241.
  • Wapenaar, K. (2004). Retrieving the elastodynamic Green's function of an arbitrary inhomogeneous medium by cross correlation. Phys. Rev. Lett., 93, 254301.

Textbooks:

  • Nakata, N., Gualtieri, L. & Fichtner, A. (2019). Seismic Ambient Noise. Cambridge University Press.
  • Shearer, P.M. (2019). Introduction to Seismology, 3rd ed. Cambridge University Press.

Software Libraries:

  • ObsPy - Python framework for seismology (Beyreuther et al., 2010; Megies et al., 2011)
  • NoisePy - High-performance ambient noise processing (Jiang & Denolle, 2020)
  • hvsrpy - HVSR processing (Vantassel, 2020; Cox et al., 2020)
In [95]:
# =============================================================================
# Environment Setup
# =============================================================================
# Install required packages (uncomment if needed):
# !pip install obspy matplotlib numpy scipy ipywidgets
# NoisePy requires Python 3.7-3.10:
# !pip install noisepy-seis

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import gridspec
from scipy import signal, fft
from scipy.signal import hilbert
import warnings
warnings.filterwarnings('ignore')

# ObsPy imports
from obspy import read, Stream, Trace, UTCDateTime, Inventory
from obspy.clients.fdsn import Client
from obspy.signal.invsim import cosine_taper

# Plotting configuration
plt.rcParams['figure.figsize'] = (14, 6)
plt.rcParams['font.size'] = 12
plt.rcParams['axes.labelsize'] = 13
plt.rcParams['axes.titlesize'] = 14

print('All imports successful.')
print(f'ObsPy version: {import_module("obspy").__version__}' if False else 'Ready to begin!')

All imports successful.
Ready to begin!

Module 1: Introduction to Seismic Ambient Noise¶

1.1 What is Ambient Seismic Noise?¶

Seismic ambient noise is the continuous background vibration of the Earth's surface, present everywhere and at all times. Unlike earthquake signals, which are transient and localized, ambient noise is a persistent, stochastic wavefield generated by a variety of natural and anthropogenic sources.

"The Earth is not static but permanently vibrating even when no strong energetic sources of vibration are acting." — Campillo & Roux (2015)

Source Classification by Frequency Band¶

Frequency Band Period Source Mechanism Reference
Infragravity waves (Earth's hum) > 30 s Ocean infragravity waves, atmospheric pressure Rhie & Romanowicz (2004)
Primary microseisms ~10–20 s (peak ~14 s) Direct ocean wave–seafloor interaction in shallow water Hasselmann (1963)
Secondary microseisms ~3–10 s (peak ~7 s) Nonlinear wave–wave interaction (opposing wave trains) Longuet-Higgins (1950)
Short-period noise < 1 s Cultural/anthropogenic: traffic, industry, wind McNamara & Buland (2004)

The Microseismic Peaks¶

The ambient noise spectrum displays two characteristic peaks observable globally:

  • Primary microseism (~0.07 Hz, period ~14 s): Generated by ocean swell interacting with the sloping sea floor in coastal areas. Energy is at the same frequency as the ocean waves. Amplitude varies with ocean wave activity and distance from coastlines.

  • Secondary microseism (~0.14 Hz, period ~7 s): The dominant feature. Generated by the nonlinear interaction of two ocean wave trains traveling in opposite directions (Longuet-Higgins, 1950). This doubles the frequency (hence the name "double-frequency microseism"). The strongest and most spatially distributed ambient noise source.

Noise Source Dynamics and Stationarity¶

A critical assumption for ambient noise cross-correlation is that noise sources are diffuse (isotropic in space) and stationary (constant over time). In practice, this is only approximately true.

Stehly et al. (2024) analyzed the dynamics of seismic noise sources in the Mediterranean Sea and showed that:

  • At 7 s period, North Atlantic sources dominate in winter across all of Europe — the noise field is stable year-round (stationarity coefficient SC ≈ 0.96–1.0)
  • In summer, Adriatic Sea and Aegean Sea contribute increasingly, causing discrete microseismic events lasting hours
  • Near Mediterranean coastlines, SC drops below 0.94 regularly, indicating non-stationary noise sources
  • Non-stationary noise introduces systematic errors into coda wave interferometry that must be corrected

Implication for ANT: Combining both winter and summer data improves measurement quality by averaging over different source distributions (Stehly et al., 2024). For dv/v monitoring, one must separate true velocity changes from apparent changes caused by noise source variations.

The New Low/High Noise Models¶

The Peterson New Low Noise Model (NLNM) and New High Noise Model (NHNM) (Peterson, 1993) define the expected range of ambient noise power spectral densities globally. Any seismic station can be characterized by comparing its noise spectrum to these reference curves using ObsPy's PPSD class.

In [96]:
# =============================================================================
# PRACTICE 1.1: Visualize the Ambient Noise Spectrum
# =============================================================================
# We download real seismic data and examine its spectral content to identify
# the primary and secondary microseismic peaks.
#
# ObsPy reference: obspy.clients.fdsn.Client
# The FDSN Client allows downloading seismic data from global data centers.
# Key methods:
#   - client.get_waveforms(network, station, location, channel, starttime, endtime)
#   - client.get_stations(...)  -> returns Inventory objects

client = Client("IRIS")

# Download 1 hour of continuous vertical-component data from a GSN station
t1 = UTCDateTime("2020-01-15T00:00:00")
t2 = t1 + 3600  # 1 hour

# Station ANMO (Albuquerque, NM, USA) - a Global Seismographic Network station 
st = client.get_waveforms("IU", "ANMO", "00", "BHZ", t1, t2)
print(st)
print(f"Sampling rate: {st[0].stats.sampling_rate} Hz")
print(f"Number of samples: {st[0].stats.npts}")

1 Trace(s) in Stream:
IU.ANMO.00.BHZ | 2020-01-15T00:00:00.019538Z - 2020-01-15T00:59:59.994538Z | 40.0 Hz, 144000 samples
Sampling rate: 40.0 Hz
Number of samples: 144000
In [97]:
# =============================================================================
# PRACTICE 1.1: Visualize the Ambient Noise Spectrum  (ETH / SED, network CH)
# =============================================================================
from obspy import UTCDateTime
from obspy.clients.fdsn import Client

client = Client("ETH")          # shortcut for http://eida.ethz.ch

t1_ch = UTCDateTime("2026-01-15T00:00:00")
t2_ch = t1_ch + 3600                  # 1 hour

# --- Step 0: never hardcode SEED IDs you haven't verified ---------------------
inv = client.get_stations(network="CH", station="*", channel="?HZ",
                          starttime=t1_ch, endtime=t2_ch, level="channel")
for net in inv:
    for sta in net:
        for cha in sta:
            print(f"{net.code}.{sta.code}.{cha.location_code or '--'}.{cha.code}"
                  f"  {cha.sample_rate:g} Hz  {sta.site.name}")

# --- Step 1: download --------------------------------------------------------
net, sta, loc, cha = "CH", "DAVOX", "", "HHZ"   # Davos; quiet vault broadband
st_ch = client.get_waveforms(net, sta, loc, cha, t1_ch, t2_ch)
print(st_ch)
print(f"Sampling rate: {st_ch[0].stats.sampling_rate} Hz")
print(f"Number of samples: {st_ch[0].stats.npts}")

CH.STIEG.BT.EHZ  200 Hz  Oberembrach, Stiegenhof, ZH
CH.BERGE.--.HHZ  200 Hz  Lenzkirch, Germany
CH.BERGE.--.LHZ  1 Hz  Lenzkirch, Germany
CH.BERGE.--.BHZ  40 Hz  Lenzkirch, Germany
CH.JAUN.--.HHZ  200 Hz  Jaun, Euschelspass, FR
CH.JAUN.--.BHZ  40 Hz  Jaun, Euschelspass, FR
CH.JAUN.--.LHZ  1 Hz  Jaun, Euschelspass, FR
CH.EMMET.--.HHZ  200 Hz  Emmethof, AG
CH.GRYON.--.HHZ  200 Hz  Gryon, VD
CH.GRYON.--.BHZ  40 Hz  Gryon, VD
CH.GRYON.--.LHZ  1 Hz  Gryon, VD
CH.HASLI.--.HHZ  200 Hz  Hasliberg, BE
CH.HASLI.--.BHZ  40 Hz  Hasliberg, BE
CH.HASLI.--.LHZ  1 Hz  Hasliberg, BE
CH.SIMPL.--.LHZ  1 Hz  Simplonpass, VS
CH.SIMPL.--.BHZ  40 Hz  Simplonpass, VS
CH.SIMPL.--.HHZ  200 Hz  Simplonpass, VS
CH.VANNI.--.BHZ  40 Hz  Val d'Anniviers, VS
CH.VANNI.--.HHZ  200 Hz  Val d'Anniviers, VS
CH.VANNI.--.LHZ  1 Hz  Val d'Anniviers, VS
CH.CHASS.--.HHZ  200 Hz  Chasseral, BE
CH.CHASS.--.BHZ  40 Hz  Chasseral, BE
CH.CHASS.--.LHZ  1 Hz  Chasseral, BE
CH.ROMAN.BT.HHZ  200 Hz  Romanshorn, Forsthaus, TG
CH.ROMAN.BT.BHZ  40 Hz  Romanshorn, Forsthaus, TG
CH.ROMAN.BT.LHZ  1 Hz  Romanshorn, Forsthaus, TG
CH.LIENZ.--.HHZ  200 Hz  Kamor, SG
CH.LIENZ.--.BHZ  40 Hz  Kamor, SG
CH.LIENZ.--.LHZ  1 Hz  Kamor, SG
CH.FUORN.--.HHZ  200 Hz  Ofenpass-Fuorn, GR
CH.FUORN.--.BHZ  40 Hz  Ofenpass-Fuorn, GR
CH.FUORN.--.LHZ  1 Hz  Ofenpass-Fuorn, GR
CH.HAUIG.--.HHZ  200 Hz  Lörrach, Café Hygge Rechberg, Deutschland
CH.HAUIG.--.BHZ  40 Hz  Lörrach, Café Hygge Rechberg, Deutschland
CH.HAUIG.--.LHZ  1 Hz  Lörrach, Café Hygge Rechberg, Deutschland
CH.SALEV.--.HHZ  200 Hz  Salève, Haute-Savoie, France
CH.SALEV.--.BHZ  40 Hz  Salève, Haute-Savoie, France
CH.SALEV.--.LHZ  1 Hz  Salève, Haute-Savoie, France
CH.LKBD2.--.HHZ  200 Hz  Leukerbad, Rinderhalte, VS
CH.LKBD2.--.BHZ  40 Hz  Leukerbad, Rinderhalte, VS
CH.LKBD2.--.LHZ  1 Hz  Leukerbad, Rinderhalte, VS
CH.MUTEZ.--.HHZ  200 Hz  Muttenz, alter Steinbruch Paradies / Bunker, BL
CH.MUTEZ.--.BHZ  40 Hz  Muttenz, alter Steinbruch Paradies / Bunker, BL
CH.MUTEZ.--.LHZ  1 Hz  Muttenz, alter Steinbruch Paradies / Bunker, BL
CH.METMA.--.HHZ  200 Hz  Mettma, Germany
CH.METMA.--.LHZ  1 Hz  Mettma, Germany
CH.METMA.--.BHZ  40 Hz  Mettma, Germany
CH.SENIN.--.HHZ  200 Hz  Lac de Senin, Sanetsch, VS
CH.SENIN.--.BHZ  40 Hz  Lac de Senin, Sanetsch, VS
CH.SENIN.--.LHZ  1 Hz  Lac de Senin, Sanetsch, VS
CH.PANIX.--.HHZ  200 Hz  Panix, Lag da Pigniu, Staumauer, GR
CH.PANIX.--.BHZ  40 Hz  Panix, Lag da Pigniu, Staumauer, GR
CH.PANIX.--.LHZ  1 Hz  Panix, Lag da Pigniu, Staumauer, GR
CH.FLACH.01.HHZ  200 Hz  Felsenburg, Rüdlingen, SH
CH.VDR.--.HHZ  200 Hz  Val di Roggiasca, GR
CH.VDR.--.BHZ  40 Hz  Val di Roggiasca, GR
CH.VDR.--.LHZ  1 Hz  Val di Roggiasca, GR
CH.SALAN.--.HHZ  200 Hz  Lac Salanfe, VS
CH.SALAN.--.BHZ  40 Hz  Lac Salanfe, VS
CH.SALAN.--.LHZ  1 Hz  Lac Salanfe, VS
CH.WIMIS.--.HHZ  200 Hz  Wimmis, BE
CH.WIMIS.--.BHZ  40 Hz  Wimmis, BE
CH.WIMIS.--.LHZ  1 Hz  Wimmis, BE
CH.MTI02.--.BHZ  40 Hz  Mont Terri, HE Schacht, JU
CH.MTI02.--.HHZ  200 Hz  Mont Terri, HE Schacht, JU
CH.MTI02.--.LHZ  1 Hz  Mont Terri, HE Schacht, JU
CH.DIX.--.HHZ  200 Hz  Grande Dixence, VS
CH.DIX.--.BHZ  40 Hz  Grande Dixence, VS
CH.DIX.--.LHZ  1 Hz  Grande Dixence, VS
CH.MESRY.--.HHZ  200 Hz  MESRY, Chemin du Plantez, France
CH.MESRY.--.BHZ  40 Hz  MESRY, Chemin du Plantez, France
CH.MESRY.--.LHZ  1 Hz  MESRY, Chemin du Plantez, France
CH.BLOTZ.BT.EHZ  200 Hz  Kappelen, Rue du Rhin, France
CH.AIGLE.--.HHZ  200 Hz  Bunker A365, Aigle, VD
CH.AIGLE.--.BHZ  40 Hz  Bunker A365, Aigle, VD
CH.AIGLE.--.LHZ  1 Hz  Bunker A365, Aigle, VD
CH.VMV.--.HHZ  200 Hz  Val Malvaglia, TI
CH.VMV.--.BHZ  40 Hz  Val Malvaglia, TI
CH.VMV.--.LHZ  1 Hz  Val Malvaglia, TI
CH.EMBD.--.BHZ  40 Hz  Mattertal, VS
CH.EMBD.--.HHZ  200 Hz  Mattertal, VS
CH.EMBD.--.LHZ  1 Hz  Mattertal, VS
CH.SAIRA.--.HHZ  200 Hz  Les Sairains Dessus, JU
CH.SAIRA.--.LHZ  1 Hz  Les Sairains Dessus, JU
CH.TRULL.--.HHZ  200 Hz  Wasserreservoir Grüt, Trüllikon, ZH
CH.TRULL.--.BHZ  40 Hz  Wasserreservoir Grüt, Trüllikon, ZH
CH.TRULL.--.LHZ  1 Hz  Wasserreservoir Grüt, Trüllikon, ZH
CH.PLONS.--.HHZ  200 Hz  Plons, SG
CH.PLONS.--.BHZ  40 Hz  Plons, SG
CH.PLONS.--.LHZ  1 Hz  Plons, SG
CH.NALPS.--.HHZ  200 Hz  Val Nalps, GR
CH.NALPS.--.BHZ  40 Hz  Val Nalps, GR
CH.NALPS.--.LHZ  1 Hz  Val Nalps, GR
CH.BOURR.--.HHZ  200 Hz  Bourrignon, JU
CH.BOURR.--.BHZ  40 Hz  Bourrignon, JU
CH.BOURR.--.LHZ  1 Hz  Bourrignon, JU
CH.SLE.--.HHZ  200 Hz  Schleitheim, SH
CH.SLE.--.BHZ  40 Hz  Schleitheim, SH
CH.SLE.--.LHZ  1 Hz  Schleitheim, SH
CH.BRANT.--.BHZ  40 Hz  Les Verrieres, NE
CH.BRANT.--.HHZ  200 Hz  Les Verrieres, NE
CH.BRANT.--.LHZ  1 Hz  Les Verrieres, NE
CH.SULZ.--.HHZ  200 Hz  Sulz, AG
CH.SULZ.--.BHZ  40 Hz  Sulz, AG
CH.SULZ.--.LHZ  1 Hz  Sulz, AG
CH.LADOL.--.HHZ  200 Hz  La Dôle, La Barillette, VD
CH.LADOL.--.BHZ  40 Hz  La Dôle, La Barillette, VD
CH.LADOL.--.LHZ  1 Hz  La Dôle, La Barillette, VD
CH.VDL.--.HHZ  200 Hz  Valle di Lei, GR
CH.VDL.--.BHZ  40 Hz  Valle di Lei, GR
CH.VDL.--.LHZ  1 Hz  Valle di Lei, GR
CH.STEIN.01.HHZ  200 Hz  Stein am Rhein, SH
CH.STEIN.01.BHZ  40 Hz  Stein am Rhein, SH
CH.STEIN.01.LHZ  1 Hz  Stein am Rhein, SH
CH.MOUTI.--.HHZ  200 Hz  Moutier, Reservoir Montange de Moutier, BE
CH.MOUTI.--.BHZ  40 Hz  Moutier, Reservoir Montange de Moutier, BE
CH.MOUTI.--.LHZ  1 Hz  Moutier, Reservoir Montange de Moutier, BE
CH.WEIN2.--.HHZ  200 Hz  Weingarten, Reservoir Immenberg, Lommis TG
CH.WEIN2.--.BHZ  40 Hz  Weingarten, Reservoir Immenberg, Lommis TG
CH.WEIN2.--.LHZ  1 Hz  Weingarten, Reservoir Immenberg, Lommis TG
CH.GIMEL.--.HHZ  200 Hz  St. Georges, Gimel, VD
CH.GIMEL.--.BHZ  40 Hz  St. Georges, Gimel, VD
CH.GIMEL.--.LHZ  1 Hz  St. Georges, Gimel, VD
CH.EMING.--.HHZ  200 Hz  Emmingen, Germany
CH.EMING.--.LHZ  1 Hz  Emmingen, Germany
CH.DAVOX.--.HHZ  200 Hz  Davos, Dischmatal, GR
CH.DAVOX.--.BHZ  40 Hz  Davos, Dischmatal, GR
CH.DAVOX.--.LHZ  1 Hz  Davos, Dischmatal, GR
CH.MUGIO.--.HHZ  200 Hz  Muggio, TI
CH.MUGIO.--.BHZ  40 Hz  Muggio, TI
CH.MUGIO.--.LHZ  1 Hz  Muggio, TI
CH.FIESA.--.HHZ  200 Hz  Fiescheralp, VS
CH.FIESA.--.BHZ  40 Hz  Fiescheralp, VS
CH.FIESA.--.LHZ  1 Hz  Fiescheralp, VS
CH.GRIMS.--.HHZ  100 Hz  Grimsel, Gerstenegg, BE
CH.GRIMS.--.LHZ  1 Hz  Grimsel, Gerstenegg, BE
CH.LLS.--.BHZ  40 Hz  Linth-Limmern, GL
CH.LLS.--.HHZ  200 Hz  Linth-Limmern, GL
CH.LLS.--.LHZ  1 Hz  Linth-Limmern, GL
CH.SLUX.--.EHZ  200 Hz  Val Lumnezia, GR
CH.OTER2.--.EHZ  500 Hz  Otterbach, 2, BS
CH.BOBI.BT.EHZ  200 Hz  Boebikon, Wasserreservoir Allmend, AG
CH.ROTHE.--.HHZ  200 Hz  Rothenfluh, BL
CH.WALHA.--.HHZ  200 Hz  Wallhausen, Germany
CH.WALHA.--.LHZ  1 Hz  Wallhausen, Germany
CH.SAVIG.--.HHZ  200 Hz  Le Moulin, Savigny, Haute-Savoie, France
CH.SAVIG.--.BHZ  40 Hz  Le Moulin, Savigny, Haute-Savoie, France
CH.SAVIG.--.LHZ  1 Hz  Le Moulin, Savigny, Haute-Savoie, France
CH.HAMIK.BT.EHZ  200 Hz  Haemikon, Daelikerfeld, LU
CH.FORET.--.EHZ  200 Hz  La Forêt, Jussy, GE
CH.AUBON.--.EHZ  200 Hz  Aubonne, Arboretum, VD
CH.TORNY.--.BHZ  40 Hz  Torny, Romont, FR
CH.TORNY.--.HHZ  200 Hz  Torny, Romont, FR
CH.TORNY.--.BHZ  40 Hz  Torny, Romont, FR
CH.TORNY.--.LHZ  1 Hz  Torny, Romont, FR
CH.COLLE.--.HHZ  200 Hz  Collex-Bossy, Route de la Vieille-Bâtie, GE
CH.COLLE.--.BHZ  40 Hz  Collex-Bossy, Route de la Vieille-Bâtie, GE
CH.COLLE.--.LHZ  1 Hz  Collex-Bossy, Route de la Vieille-Bâtie, GE
CH.MFERR.--.HHZ  200 Hz  Prayon, Val Ferret, La Fouly, VS
CH.MFERR.--.BHZ  40 Hz  Prayon, Val Ferret, La Fouly, VS
CH.MFERR.--.LHZ  1 Hz  Prayon, Val Ferret, La Fouly, VS
CH.BNALP.--.BHZ  40 Hz  Bannalp, NW
CH.BNALP.--.HHZ  200 Hz  Bannalp, NW
CH.BNALP.--.LHZ  1 Hz  Bannalp, NW
CH.ONNEN.--.EHZ  200 Hz  Onnens, Les Côtes, VD
CH.VINZL.--.HHZ  200 Hz  Vinzel, Chemin de la Chaponnière, VD
CH.VINZL.--.BHZ  40 Hz  Vinzel, Chemin de la Chaponnière, VD
CH.VINZL.--.LHZ  1 Hz  Vinzel, Chemin de la Chaponnière, VD
CH.ZUR.--.HHZ  200 Hz  Zuerich, Degenried, ZH
CH.ZUR.--.BHZ  40 Hz  Zuerich, Degenried, ZH
CH.ZUR.--.LHZ  1 Hz  Zuerich, Degenried, ZH
CH.EWZT2.--.EHZ  200 Hz  Wettswil, ZH
CH.BAULM.--.HHZ  200 Hz  Baulmes, Aiguilles de Baulmes, VD
CH.CHALL.BT.EHZ  200 Hz  Challoux, Bernex, GE
CH.MMK.--.HHZ  200 Hz  Mattmark, VS
CH.MMK.--.BHZ  40 Hz  Mattmark, VS
CH.MMK.--.LHZ  1 Hz  Mattmark, VS
CH.PERON.--.HHZ  200 Hz  Péron, France
CH.PERON.--.BHZ  40 Hz  Péron, France
CH.PERON.--.LHZ  1 Hz  Péron, France
CH.OPENS.--.EHZ  200 Hz  Oppens, Waadt, VD
CH.BALST.--.HHZ  200 Hz  Balsthal, SO
CH.BALST.--.BHZ  40 Hz  Balsthal, SO
CH.BALST.--.LHZ  1 Hz  Balsthal, SO
CH.MUO.--.HHZ  200 Hz  Muotathal, SZ
CH.MUO.--.BHZ  40 Hz  Muotathal, SZ
CH.MUO.--.LHZ  1 Hz  Muotathal, SZ
CH.WILA.--.HHZ  200 Hz  Wila, ZH
CH.WILA.--.BHZ  40 Hz  Wila, ZH
CH.WILA.--.LHZ  1 Hz  Wila, ZH
CH.DAGMA.--.HHZ  200 Hz  Dagmersellen, Lutertal, LU
CH.DAGMA.--.LHZ  1 Hz  Dagmersellen, Lutertal, LU
CH.DAGMA.--.BHZ  40 Hz  Dagmersellen, Lutertal, LU
CH.SGT00.BT.HHZ  200 Hz  Sennhuelsen, SG
CH.SGT00.BT.BHZ  40 Hz  Sennhuelsen, SG
CH.SGT00.BT.LHZ  1 Hz  Sennhuelsen, SG
CH.LAUCH.--.HHZ  200 Hz  Lauchernalp - Loetschental, VS
CH.LAUCH.--.LHZ  1 Hz  Lauchernalp - Loetschental, VS
CH.LAUCH.--.BHZ  40 Hz  Lauchernalp - Loetschental, VS
CH.ACB.--.HHZ  200 Hz  Klingnau, Acheberg, AG
CH.ACB.--.BHZ  40 Hz  Klingnau, Acheberg, AG
CH.ACB.--.LHZ  1 Hz  Klingnau, Acheberg, AG
CH.LASAR.--.HHZ  200 Hz  La Sarraz, Réservoir des Aleveys, VD
CH.LASAR.--.BHZ  40 Hz  La Sarraz, Réservoir des Aleveys, VD
CH.LASAR.--.LHZ  1 Hz  La Sarraz, Réservoir des Aleveys, VD
CH.GOURZ.--.HHZ  200 Hz  Tour de Gourze, Réservoir de la Tour-de-Gource, VD
CH.GOURZ.--.BHZ  40 Hz  Tour de Gourze, Réservoir de la Tour-de-Gource, VD
CH.GOURZ.--.LHZ  1 Hz  Tour de Gourze, Réservoir de la Tour-de-Gource, VD
CH.FUSIO.--.HHZ  200 Hz  Fusio, TI
CH.FUSIO.--.BHZ  40 Hz  Fusio, TI
CH.FUSIO.--.LHZ  1 Hz  Fusio, TI
CH.CHAMB.--.HHZ  200 Hz  Chamblon, Rue du Cossaux, VD
CH.CHAMB.--.BHZ  40 Hz  Chamblon, Rue du Cossaux, VD
CH.CHAMB.--.LHZ  1 Hz  Chamblon, Rue du Cossaux, VD
CH.WGT.--.HHZ  200 Hz  Wägital, SZ
CH.WGT.--.BHZ  40 Hz  Wägital, SZ
CH.WGT.--.LHZ  1 Hz  Wägital, SZ
CH.BERNI.--.HHZ  200 Hz  Berninapass, GR
CH.BERNI.--.BHZ  40 Hz  Berninapass, GR
CH.BERNI.--.LHZ  1 Hz  Berninapass, GR
CH.ILLEZ.--.HHZ  200 Hz  Chalet Vouargne Bourlo, Wallis, VS
CH.ILLEZ.--.BHZ  40 Hz  Chalet Vouargne Bourlo, Wallis, VS
CH.ILLEZ.--.LHZ  1 Hz  Chalet Vouargne Bourlo, Wallis, VS
1 Trace(s) in Stream:
CH.DAVOX..HHZ | 2026-01-15T00:00:00.000000Z - 2026-01-15T01:00:00.000000Z | 200.0 Hz, 720001 samples
Sampling rate: 200.0 Hz
Number of samples: 720001

Swiss Seismological Service has an state of the art metadata and generall information portal for all of its stations, worth looking at it. For station BAS https://stations.seismo.ethz.ch/en/station-information/station-details/station-given-networkcode-and-stationcode/index.html?networkcode=CH&stationcode=BAS

In [98]:
# =============================================================================
# Compute and plot the Power Spectral Density (PSD)
# =============================================================================
# This reveals the characteristic microseismic peaks

tr = st[0].copy()
tr.detrend('demean')
tr.detrend('linear')
tr.taper(0.05)

# Compute PSD using Welch's method
fs = tr.stats.sampling_rate
nperseg = int(fs * 100)  # 100-second windows
freqs, psd = signal.welch(tr.data, fs=fs, nperseg=nperseg, noverlap=nperseg//2)

# Convert to period
periods = 1.0 / freqs[1:]  # Skip zero frequency

fig, axes = plt.subplots(1, 2, figsize=(16, 6))

# Plot 1: PSD vs Frequency
axes[0].semilogy(freqs[1:], psd[1:], 'b-', linewidth=0.8)
axes[0].set_xlabel('Frequency (Hz)')
axes[0].set_ylabel('Power Spectral Density')
axes[0].set_title('Ambient Noise Power Spectrum')
axes[0].set_xlim(0.01, 2.0)
axes[0].axvspan(0.05, 0.1, alpha=0.2, color='red', label='Primary microseism')
axes[0].axvspan(0.1, 0.3, alpha=0.2, color='green', label='Secondary microseism')
axes[0].axvspan(0.5, 2.0, alpha=0.2, color='orange', label='Anthropogenic')
axes[0].legend(fontsize=10)
axes[0].grid(True, alpha=0.3)

# Plot 2: PSD vs Period
axes[1].loglog(periods, psd[1:], 'b-', linewidth=0.8)
axes[1].set_xlabel('Period (s)')
axes[1].set_ylabel('Power Spectral Density')
axes[1].set_title('Ambient Noise Power Spectrum (Period Domain)')
axes[1].set_xlim(0.5, 100)
axes[1].axvspan(10, 20, alpha=0.2, color='red', label='Primary microseism (~14 s)')
axes[1].axvspan(3, 10, alpha=0.2, color='green', label='Secondary microseism (~7 s)')
axes[1].legend(fontsize=10)
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.suptitle('Identifying Microseismic Peaks in Ambient Noise (Station IU.ANMO)', 
             y=1.02, fontsize=14, fontweight='bold')
plt.show()
No description has been provided for this image
In [99]:
# =============================================================================
# Compute and plot the Power Spectral Density (PSD)
# =============================================================================
# This reveals the characteristic microseismic peaks

tr_ch = st_ch[0].copy()
tr_ch.detrend('demean')
tr_ch.detrend('linear')
tr_ch.taper(0.05)

# Compute PSD using Welch's method
fs = tr_ch.stats.sampling_rate
nperseg = int(fs * 100)  # 100-second windows
freqs, psd = signal.welch(tr_ch.data, fs=fs, nperseg=nperseg, noverlap=nperseg//2)

# Convert to period
periods = 1.0 / freqs[1:]  # Skip zero frequency

fig, axes = plt.subplots(1, 2, figsize=(16, 6))

# Plot 1: PSD vs Frequency
axes[0].semilogy(freqs[1:], psd[1:], 'b-', linewidth=0.8)
axes[0].set_xlabel('Frequency (Hz)')
axes[0].set_ylabel('Power Spectral Density')
axes[0].set_title('Ambient Noise Power Spectrum')
axes[0].set_xlim(0.01, 2.0)
axes[0].axvspan(0.05, 0.1, alpha=0.2, color='red', label='Primary microseism')
axes[0].axvspan(0.1, 0.3, alpha=0.2, color='green', label='Secondary microseism')
axes[0].axvspan(0.5, 2.0, alpha=0.2, color='orange', label='Anthropogenic')
axes[0].legend(fontsize=10)
axes[0].grid(True, alpha=0.3)

# Plot 2: PSD vs Period
axes[1].loglog(periods, psd[1:], 'b-', linewidth=0.8)
axes[1].set_xlabel('Period (s)')
axes[1].set_ylabel('Power Spectral Density')
axes[1].set_title('Ambient Noise Power Spectrum (Period Domain)')
axes[1].set_xlim(0.5, 100)
axes[1].axvspan(10, 20, alpha=0.2, color='red', label='Primary microseism (~14 s)')
axes[1].axvspan(3, 10, alpha=0.2, color='green', label='Secondary microseism (~7 s)')
axes[1].legend(fontsize=10)
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.suptitle('Identifying Microseismic Peaks in Ambient Noise (Station CH.DAVOX)', 
             y=1.02, fontsize=14, fontweight='bold')
plt.show()
No description has been provided for this image

Key Observation¶

The spectrum shows clear peaks at the primary (~0.07 Hz / ~14 s) and secondary (~0.14 Hz / ~7 s) microseismic periods. These are the dominant frequencies in ambient noise seismology and are generated by ocean-atmosphere interactions thousands of kilometers away.

"In the band of 0.3-0.05 Hz, the ambient noise is often referred to as 'microseisms'. In the microseism spectral band, the noise is dominated by surface waves, which are predominantly Rayleigh waves." -- Campillo & Roux (2015)

Module 2: Seismic Data Handling with ObsPy¶

Refresh¶

Seismo-live is an interactive online learning platform for seismology. It utilizes Jupyter Notebooks to introduce some seismic principle hands-on with some coding, I highly recommand to have a try at them. Obspy, nowadays, is the base of close to all packages related to seismology. For ambient noise, it is directly integrate to Noisepy (Jiang, C. and Denolle, 2020) and MSNoise (Lecocq et al., 2014)

2.1 Core Data Structures¶

ObsPy (Beyreuther et al., 2010; Megies et al., 2011) is the standard Python framework for seismological data processing. Its core objects are:

  • Trace: A single continuous time series with metadata (stats attribute containing network, station, channel, sampling_rate, starttime, etc.)
  • Stream: An ordered collection of Trace objects (like a list of seismograms)
  • UTCDateTime: Precise time representation for seismology
  • Inventory: Station metadata (coordinates, instrument response, etc.)

ObsPy Architecture for Ambient Noise¶

ObsPy
 ├── obspy.core          # Trace, Stream, UTCDateTime
 ├── obspy.clients.fdsn  # Data download from FDSN web services
 ├── obspy.signal        # Signal processing (filter, taper, etc.)
 ├── obspy.io            # I/O for various formats (SAC, miniSEED, etc.)
 └── obspy.imaging       # Plotting utilities

FDSN Web Services¶

ObsPy can access seismic data from global data centers through the International Federation of Digital Seismograph Networks (FDSN) web services:

Data Center Code Region
IRIS DMC IRIS Global
ORFEUS/EIDA ODC Europe
GFZ GFZ Germany
ETH ETH Switzerland
RESIF RESIF France
INGV INGV Italy
In [100]:
# =============================================================================
# PRACTICE 2.1: Working with ObsPy Core Objects
# =============================================================================

# --- UTCDateTime: Precise time handling ---
t = UTCDateTime("2020-06-15T12:30:00.0")
print(f"UTCDateTime object: {t}")
print(f"Julian day: {t.julday}")
print(f"Timestamp (POSIX): {t.timestamp}")
print(f"Add 1 hour: {t + 3600}")
print()

# --- Trace: Single time series ---
# Create a synthetic trace
npts = 1000
sampling_rate = 100.0
data = np.random.randn(npts)  # Random noise

tr = Trace(data=data)
tr.stats.sampling_rate = sampling_rate
tr.stats.network = 'XX'
tr.stats.station = 'TEST'
tr.stats.channel = 'BHZ'
tr.stats.starttime = UTCDateTime("2020-01-01T00:00:00")

print(f"Trace: {tr}")
print(f"Stats: {tr.stats}")
print()

# --- Stream: Collection of traces ---
st = Stream([tr])
print(f"Stream: {st}")

# Key Trace methods for ambient noise work:
tr_copy = tr.copy()  # Always work on copies!
tr_copy.detrend('demean')    # Remove mean
tr_copy.detrend('linear')    # Remove linear trend
tr_copy.taper(0.05)          # Apply cosine taper (5%)
tr_copy.filter('bandpass', freqmin=0.01, freqmax=1.0)  # Bandpass filter
print("\nAfter preprocessing: trace is ready for cross-correlation")

UTCDateTime object: 2020-06-15T12:30:00.000000Z
Julian day: 167
Timestamp (POSIX): 1592224200.0
Add 1 hour: 2020-06-15T13:30:00.000000Z

Trace: XX.TEST..BHZ | 2020-01-01T00:00:00.000000Z - 2020-01-01T00:00:09.990000Z | 100.0 Hz, 1000 samples
Stats:          network: XX
         station: TEST
        location: 
         channel: BHZ
       starttime: 2020-01-01T00:00:00.000000Z
         endtime: 2020-01-01T00:00:09.990000Z
   sampling_rate: 100.0
           delta: 0.01
            npts: 1000
           calib: 1.0

Stream: 1 Trace(s) in Stream:
XX.TEST..BHZ | 2020-01-01T00:00:00.000000Z - 2020-01-01T00:00:09.990000Z | 100.0 Hz, 1000 samples

After preprocessing: trace is ready for cross-correlation
In [101]:
# =============================================================================
# PRACTICE 2.2: Download and Visualize Real Seismic Data
# =============================================================================
# Download data from two stations for later cross-correlation
# We choose a station pair ~1400 km apart so the expected surface wave
# arrival (~400 s at 3.5 km/s) fits within a reasonable lag window.

client = Client("IRIS")

# Parameters
starttime = UTCDateTime("2020-01-15T00:00:00")
endtime = starttime + 86400  # 1 day of data

# Station 1: ANMO (Albuquerque, NM)
st1 = client.get_waveforms("IU", "ANMO", "00", "BHZ", starttime, endtime)
# Station 2: CCM (Cathedral Cave, MO) — ~1400 km from ANMO
st2 = client.get_waveforms("IU", "CCM", "00", "BHZ", starttime, endtime)

print(f"Station 1: {st1[0].stats.station}")
print(f"Station 2: {st2[0].stats.station}")

# Get station coordinates
inv1 = client.get_stations(network="IU", station="ANMO", level="station")
inv2 = client.get_stations(network="IU", station="CCM", level="station")

lat1 = inv1[0][0].latitude
lon1 = inv1[0][0].longitude
lat2 = inv2[0][0].latitude
lon2 = inv2[0][0].longitude

# Calculate inter-station distance
from obspy.geodetics import gps2dist_azimuth
dist_m, az, baz = gps2dist_azimuth(lat1, lon1, lat2, lon2)
dist_km = dist_m / 1000.0

print(f"\nStation ANMO: ({lat1:.2f}, {lon1:.2f})")
print(f"Station CCM:  ({lat2:.2f}, {lon2:.2f})")
print(f"Inter-station distance: {dist_km:.1f} km")
print(f"Expected Rayleigh wave travel time (~3.5 km/s): {dist_km/3.5:.0f} s")
print(f"Azimuth: {az:.1f} deg")

Station 1: ANMO
Station 2: HRV

Station ANMO: (34.95, -106.46)
Station HRV:  (42.51, -71.56)
Inter-station distance: 3124.4 km
Azimuth: 63.9 deg
In [102]:
# =============================================================================
# Visualize raw waveforms
# =============================================================================
fig, axes = plt.subplots(2, 1, figsize=(16, 8), sharex=True)

tr1 = st1[0].copy()
tr2 = st2[0].copy()

t_axis1 = np.arange(tr1.stats.npts) / tr1.stats.sampling_rate
t_axis2 = np.arange(tr2.stats.npts) / tr2.stats.sampling_rate

axes[0].plot(t_axis1 / 3600, tr1.data, 'k-', linewidth=0.3)
axes[0].set_title(f'Station {tr1.stats.station} (Raw)', fontweight='bold')
axes[0].set_ylabel('Counts')

axes[1].plot(t_axis2 / 3600, tr2.data, 'k-', linewidth=0.3)
axes[1].set_title(f'Station {tr2.stats.station} (Raw)', fontweight='bold')
axes[1].set_ylabel('Counts')
axes[1].set_xlabel('Time (hours)')

plt.suptitle(f'One Day of Continuous Seismic Data ({starttime.date})', 
             y=1.01, fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()

print("Note: The raw data contains ambient noise, possible earthquakes,")
print("instrument glitches, and other transients that must be removed.")
No description has been provided for this image 
Note: The raw data contains ambient noise, possible earthquakes,
instrument glitches, and other transients that must be removed.

Module 3: The Cross-Correlation Theorem¶

3.1 Mathematical Definition¶

The cross-correlation is a fundamental operation that measures the similarity between two signals as a function of a time lag $\tau$. For two continuous signals $u_1(t)$ and $u_2(t)$ of duration $T$, the cross-correlation is defined as (Campillo & Roux, 2015):

$$C(\vec{r}_1, \vec{r}_2; t) = C_{1,2}(t) = \frac{1}{T} \int_0^T u_1(\tau) \, u_2(t + \tau) \, d\tau \qquad [1]$$

Or equivalently, in the frequency domain:

$$C(\vec{r}_1, \vec{r}_2; \omega) = C_{1,2}(\omega) = u_1(\omega) \, u_2^*(\omega) \qquad [2]$$

where $u_2^*(\omega)$ is the complex conjugate of the Fourier transform of $u_2(t)$.

This is the key equation used in NoisePy (Jiang & Denolle, 2020): $$\hat{C}_{AB}(f) = \hat{U}_A^*(f) \hat{U}_B(f) \qquad \text{[Eq. 1, Jiang & Denolle, 2020]}$$

Physical Interpretation¶

The cross-correlation function measures the resemblance of two signals:

  • When signals are similar, the cross-correlation peaks at the time delay that separates them
  • The ambient noise signals recorded at distant stations are considered uncorrelated since they result from variable interferences between numerous waves of different types emitted by different sources
  • Noise-correlation methods aim at extracting the slight coherent part of the signals that contains deterministic information on wave propagation between the two stations (Campillo & Roux, 2015)

Cross-Correlation vs. Convolution¶

The cross-correlation is closely related to convolution:

  • Convolution: $f * g = \int f(\tau) g(t - \tau) d\tau$ (one signal is time-reversed)
  • Cross-correlation: $f \star g = \int f(\tau) g(t + \tau) d\tau$ (no time reversal)

In the frequency domain:

  • Convolution $\Leftrightarrow$ Multiplication: $\mathcal{F}\{f * g\} = F(\omega) \cdot G(\omega)$
  • Cross-correlation $\Leftrightarrow$ Conjugate multiplication: $\mathcal{F}\{f \star g\} = F^*(\omega) \cdot G(\omega)$

This is the Cross-Correlation Theorem (also called the Wiener-Khinchin theorem for auto-correlation).

In [103]:
# =============================================================================
# PRACTICE 3.1: Understanding Cross-Correlation with Synthetic Signals
# =============================================================================
# Demonstrate the cross-correlation theorem step by step
#
# Convention: C_AB(tau) = integral A(t) * B(t + tau) dt
# When B is a delayed copy of A (delay = delta), the CC peaks at tau = -delta
# because B(t + tau) aligns with A(t) when tau = -delta.
# 
# In seismology, we use C_AB(f) = conj(F_A) * F_B, which gives the same result.

# Create two signals: identical wavelet with a known time delay
dt = 0.01  # sampling interval (s)
t = np.arange(0, 10, dt)
n = len(t)

# Create a Ricker wavelet (Mexican hat wavelet)
def ricker_wavelet(t, t0, f0):
    """Generate a Ricker wavelet centered at t0 with dominant frequency f0."""
    u = (np.pi * f0 * (t - t0)) ** 2
    return (1 - 2 * u) * np.exp(-u)

# Signal at station A: wavelet arrives at t=3s  
# Signal at station B: same wavelet arrives at t=5s (2s delay: B is LATER)
signal_A = ricker_wavelet(t, 3.0, 2.0)
signal_B = ricker_wavelet(t, 5.0, 2.0)

# Add some noise
np.random.seed(42)
signal_A += 0.1 * np.random.randn(n)
signal_B += 0.1 * np.random.randn(n)

# --- Frequency-domain cross-correlation (Cross-Correlation Theorem) ---
# C_AB(f) = conj(F_A(f)) * F_B(f)
nfft = 2 * n  # zero-padding for linear (non-circular) correlation
FA = np.fft.fft(signal_A, n=nfft)
FB = np.fft.fft(signal_B, n=nfft)
CC_freq = np.conj(FA) * FB  # Cross-correlation theorem!
cc_freq = np.real(np.fft.ifft(CC_freq))
cc_freq = np.fft.fftshift(cc_freq)
lags = np.arange(-nfft//2, nfft//2) * dt

# Peak detection
peak_lag = lags[np.argmax(cc_freq)]

# --- Plot results ---
fig, axes = plt.subplots(3, 1, figsize=(14, 10))

# Panel 1: Input signals
axes[0].plot(t, signal_A, 'b-', label='Station A (wavelet at t=3s)', linewidth=1.5)
axes[0].plot(t, signal_B, 'r-', label='Station B (wavelet at t=5s)', linewidth=1.5)
axes[0].axvline(3.0, color='b', linestyle='--', alpha=0.5)
axes[0].axvline(5.0, color='r', linestyle='--', alpha=0.5)
axes[0].annotate(r'$\Delta t = 2$ s', xy=(4, 0.8), fontsize=14, fontweight='bold',
                ha='center', color='green')
axes[0].set_title('Input Signals', fontweight='bold')
axes[0].set_xlabel('Time (s)')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# Panel 2: Cross-correlation
axes[1].plot(lags, cc_freq / np.max(np.abs(cc_freq)), 'k-', linewidth=1.5)
axes[1].axvline(peak_lag, color='red', linestyle='--', linewidth=2, 
                label=f'Peak at lag = {peak_lag:.2f} s')
axes[1].set_title(r'Cross-Correlation $C_{AB}(\tau) = \mathcal{{F}}^{{-1}}[U_A^* \cdot U_B]$', 
                  fontweight='bold')
axes[1].set_xlabel('Lag $\\tau$ (s)')
axes[1].set_ylabel('Normalized CC')
axes[1].set_xlim(-5, 5)
axes[1].legend(fontsize=11)
axes[1].grid(True, alpha=0.3)

# Panel 3: Amplitude spectra
freqs_plot = np.fft.rfftfreq(nfft, d=dt)
axes[2].semilogy(freqs_plot, np.abs(np.fft.rfft(signal_A, n=nfft)), 'b-', 
                  alpha=0.7, label='|F_A(f)|')
axes[2].semilogy(freqs_plot, np.abs(np.fft.rfft(signal_B, n=nfft)), 'r--', 
                  alpha=0.7, label='|F_B(f)|')
axes[2].semilogy(freqs_plot, np.abs(CC_freq[:nfft//2+1]), 'k-', 
                  alpha=0.7, label='|C_AB(f)| = |F_A* · F_B|')
axes[2].set_title('Amplitude Spectra', fontweight='bold')
axes[2].set_xlabel('Frequency (Hz)')
axes[2].set_ylabel('Amplitude')
axes[2].set_xlim(0, 10)
axes[2].legend()
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

print(f"Peak lag: {peak_lag:.4f} s")
print(f"B is delayed by 2.0 s relative to A")
print(f"CC peak at {peak_lag:.1f} s confirms this: positive lag means")
print(f"B's signal arrives later than A's signal.")
print(f"")
print(f"The Cross-Correlation Theorem:")
print(f"  C_AB(f) = U_A*(f) · U_B(f)")
print(f"The time-domain CC equals the inverse FFT of the conjugate product.")
print(f"This is how NoisePy and all ANT codes compute CCs efficiently.")
No description has been provided for this image 
Peak lag: 2.0000 s
B is delayed by 2.0 s relative to A
CC peak at 2.0 s confirms this: positive lag means
B's signal arrives later than A's signal.

The Cross-Correlation Theorem:
  C_AB(f) = U_A*(f) · U_B(f)
The time-domain CC equals the inverse FFT of the conjugate product.
This is how NoisePy and all ANT codes compute CCs efficiently.
In [104]:
# =============================================================================
# PRACTICE 3.2: Cross-Correlation of Noise — The Emergence of Signal
# =============================================================================
# This is the KEY experiment: show that cross-correlating random noise
# recorded at two points can reveal the travel time between them.
#
# We simulate a 1D medium with noise sources surrounding two receivers.
# The key physics: each source emits noise that arrives at both receivers
# with different delays. The cross-correlation extracts the DIFFERENTIAL
# travel time — i.e., the Green's function between the receivers.
#
# We work in the frequency domain for accurate fractional-sample delays:
#   delay tau <=> phase shift exp(-i 2*pi*f*tau)

np.random.seed(123)

# Parameters
fs = 100.0        # Sampling rate (Hz)
duration = 600.0  # Duration (s) — longer = better convergence
c = 3.0           # Rayleigh wave velocity (km/s)
n_sources = 500   # Number of random noise sources (isotropic distribution)
n_samples = int(duration * fs)
nfft = 2 ** int(np.ceil(np.log2(n_samples)))

# Receiver positions (km)
rx_A = 0.0
rx_B = 10.0  # 10 km apart -> expected travel time = 10/3 = 3.33 s

# Noise sources uniformly distributed along a line surrounding receivers
source_positions = np.random.uniform(-200, 200, n_sources)  # km

# Frequency axis
freqs = np.fft.rfftfreq(nfft, d=1.0/fs)

# Accumulate receiver spectra
FA_total = np.zeros(len(freqs), dtype=complex)
FB_total = np.zeros(len(freqs), dtype=complex)

# Bandpass shape (simulates surface wave frequency content: 0.2-2 Hz)
bp_filter = np.zeros(len(freqs))
f_mask = (freqs >= 0.2) & (freqs <= 2.0)
bp_filter[f_mask] = 1.0

for sx in source_positions:
    # Each source has random spectral content (unique noise realization)
    source_spec = (np.random.randn(len(freqs)) + 1j * np.random.randn(len(freqs))) * bp_filter
    
    # Distances and travel times to each receiver
    dist_A = abs(sx - rx_A)
    dist_B = abs(sx - rx_B)
    tt_A = dist_A / c
    tt_B = dist_B / c
    
    # Geometric spreading: amplitude ~ 1/sqrt(r) for surface waves in 2D
    amp_A = 1.0 / max(np.sqrt(dist_A), 0.3)
    amp_B = 1.0 / max(np.sqrt(dist_B), 0.3)
    
    # Apply propagation delay as phase shift: exp(-i 2*pi*f*t)
    # This is exact (no integer-sample truncation)
    FA_total += amp_A * source_spec * np.exp(-2j * np.pi * freqs * tt_A)
    FB_total += amp_B * source_spec * np.exp(-2j * np.pi * freqs * tt_B)

# Transform back to time domain for plotting
rec_A = np.fft.irfft(FA_total, n=nfft)[:n_samples]
rec_B = np.fft.irfft(FB_total, n=nfft)[:n_samples]

# Cross-correlation in frequency domain: C_AB(f) = conj(FA) * FB
CC = np.conj(FA_total) * FB_total
cc = np.fft.irfft(CC, n=nfft)

# Build centered lag axis [-T/2, T/2]
lags = np.arange(nfft) / fs
lags[lags > nfft / (2*fs)] -= nfft / fs
sort_idx = np.argsort(lags)
lags = lags[sort_idx]
cc = cc[sort_idx]

expected_tt = (rx_B - rx_A) / c  # 3.33 s

# --- Plot results ---
fig, axes = plt.subplots(3, 1, figsize=(16, 10))

t_axis = np.arange(n_samples) / fs
axes[0].plot(t_axis, rec_A, 'b-', linewidth=0.3, alpha=0.7)
axes[0].set_title('Recording at Station A (sum of 500 random noise sources)', fontweight='bold')
axes[0].set_ylabel('Amplitude')
axes[0].set_xlim(0, 60)  # Show first 60 seconds for clarity

axes[1].plot(t_axis, rec_B, 'r-', linewidth=0.3, alpha=0.7)
axes[1].set_title('Recording at Station B (looks like random noise!)', fontweight='bold')
axes[1].set_ylabel('Amplitude')
axes[1].set_xlabel('Time (s)')
axes[1].set_xlim(0, 60)

axes[2].plot(lags, cc / np.max(np.abs(cc)), 'k-', linewidth=1.0)
axes[2].axvline(expected_tt, color='r', linestyle='--', linewidth=2, 
                label=f'Expected: +{expected_tt:.2f} s (causal: A$\\to$B)')
axes[2].axvline(-expected_tt, color='b', linestyle='--', linewidth=2, 
                label=f'Expected: $-${expected_tt:.2f} s (acausal: B$\\to$A)')
axes[2].set_title("Cross-Correlation: Green's Function Emerges from Noise!", 
                  fontweight='bold')
axes[2].set_xlabel('Lag time (s)')
axes[2].set_ylabel('Normalized CC')
axes[2].set_xlim(-15, 15)
axes[2].legend(fontsize=11)
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

print(f"\nExpected inter-station travel time: {expected_tt:.2f} s")
print(f"Station separation: {rx_B - rx_A:.1f} km")
print(f"Rayleigh wave velocity: {c:.1f} km/s")
print(f"\n*** The cross-correlation reveals clear peaks at +/- {expected_tt:.2f} s ***")
print(f"This corresponds to the causal and acausal Green's function!")
print(f"\nPhysics: Each noise source creates a coherent differential delay")
print(f"between the two stations. Summing over many isotropic sources,")
print(f"only the stationary-phase contributions (sources along the")
print(f"inter-station axis) survive — producing the surface wave arrival.")
No description has been provided for this image 
Expected inter-station travel time: 3.33 s
Station separation: 10.0 km
Rayleigh wave velocity: 3.0 km/s

*** The cross-correlation reveals clear peaks at +/- 3.33 s ***
This corresponds to the causal and acausal Green's function!

Physics: Each noise source creates a coherent differential delay
between the two stations. Summing over many isotropic sources,
only the stationary-phase contributions (sources along the
inter-station axis) survive — producing the surface wave arrival.

Key Result: Signal Emerges from Noise!¶

Even though both recordings look like pure random noise, their cross-correlation reveals clear peaks at $\pm t = d/c$, where $d$ is the inter-station distance and $c$ is the wave velocity. This is the empirical Green's function.

The two peaks correspond to:

  • Positive lag ($+\tau$): The causal Green's function (waves traveling from A to B)
  • Negative lag ($-\tau$): The acausal Green's function (waves traveling from B to A)

"The correlation function of the fields recorded at two points is used as a virtual seismogram for which the source is acting at one point and the receiver at the second." -- Campillo & Roux (2015)

In [105]:
# =============================================================================
# PRACTICE 3.3: Effect of Source Distribution on Cross-Correlation
# =============================================================================
# Show how the azimuthal distribution of noise sources affects the CC
# This illustrates the "stationary phase" and "end-fire lobes" concepts
# from Campillo & Roux (2015)

def simulate_2d_noise_cc(source_azimuths, n_sources_per_az=10, title=''):
    """
    Simulate 2D ambient noise cross-correlation with sources at specific azimuths.
    Uses frequency-domain phase shifts for accurate wave propagation.
    """
    fs = 100.0
    duration = 600.0
    c = 3.0  # km/s
    n_samples = int(duration * fs)
    nfft = 2 ** int(np.ceil(np.log2(n_samples)))
    
    # Receiver positions
    rA = np.array([0.0, 0.0])   # Station A at origin
    rB = np.array([10.0, 0.0])  # Station B at 10 km east
    
    freqs = np.fft.rfftfreq(nfft, d=1.0/fs)
    FA_total = np.zeros(len(freqs), dtype=complex)
    FB_total = np.zeros(len(freqs), dtype=complex)
    
    # Bandpass filter shape (surface wave band)
    bp = np.zeros(len(freqs))
    bp[(freqs >= 0.2) & (freqs <= 2.0)] = 1.0
    
    for az in source_azimuths:
        for _ in range(n_sources_per_az):
            # Source at distance 30-100 km from midpoint
            dist = np.random.uniform(30, 100)
            az_rad = np.radians(az)
            midpoint = (rA + rB) / 2
            sx = midpoint[0] + dist * np.cos(az_rad)
            sy = midpoint[1] + dist * np.sin(az_rad)
            
            # Random source spectrum
            source_spec = (np.random.randn(len(freqs)) + 
                          1j * np.random.randn(len(freqs))) * bp
            
            # Distances and travel times
            d_A = np.sqrt((sx - rA[0])**2 + (sy - rA[1])**2)
            d_B = np.sqrt((sx - rB[0])**2 + (sy - rB[1])**2)
            tt_A = d_A / c
            tt_B = d_B / c
            
            # Geometric spreading (2D surface waves)
            amp_A = 1.0 / max(np.sqrt(d_A), 1.0)
            amp_B = 1.0 / max(np.sqrt(d_B), 1.0)
            
            # Phase shift for propagation delay
            FA_total += amp_A * source_spec * np.exp(-2j * np.pi * freqs * tt_A)
            FB_total += amp_B * source_spec * np.exp(-2j * np.pi * freqs * tt_B)
    
    # Cross-correlation
    CC = np.conj(FA_total) * FB_total
    cc = np.fft.irfft(CC, n=nfft)
    
    # Centered lag axis
    lags = np.arange(nfft) / fs
    lags[lags > nfft / (2*fs)] -= nfft / fs
    sort_idx = np.argsort(lags)
    lags = lags[sort_idx]
    cc = cc[sort_idx]
    
    return lags, cc

# Three scenarios
np.random.seed(42)
fig, axes = plt.subplots(1, 3, figsize=(18, 5))

expected_tt = 10.0 / 3.0

# Scenario 1: Uniform source distribution (ideal — isotropic)
lags, cc = simulate_2d_noise_cc(np.arange(0, 360, 10))
axes[0].plot(lags, cc / np.max(np.abs(cc)), 'k-', linewidth=0.8)
axes[0].axvline(expected_tt, color='r', linestyle='--', alpha=0.7, label=f'+{expected_tt:.1f} s')
axes[0].axvline(-expected_tt, color='b', linestyle='--', alpha=0.7, label=f'-{expected_tt:.1f} s')
axes[0].set_xlim(-15, 15)
axes[0].set_title('Uniform Sources\n(All azimuths — isotropic)', fontweight='bold')
axes[0].set_xlabel('Lag (s)')
axes[0].set_ylabel('Normalized CC')
axes[0].legend(fontsize=9)
axes[0].grid(True, alpha=0.3)

# Scenario 2: Sources only from one side (end-fire lobe)
lags, cc = simulate_2d_noise_cc(np.arange(-30, 30, 5), n_sources_per_az=20)
axes[1].plot(lags, cc / np.max(np.abs(cc)), 'k-', linewidth=0.8)
axes[1].axvline(expected_tt, color='r', linestyle='--', alpha=0.7)
axes[1].axvline(-expected_tt, color='b', linestyle='--', alpha=0.7)
axes[1].set_xlim(-15, 15)
axes[1].set_title('Sources from East Only\n(One end-fire lobe — asymmetric CC)', fontweight='bold')
axes[1].set_xlabel('Lag (s)')
axes[1].grid(True, alpha=0.3)

# Scenario 3: Sources perpendicular (no coherent signal expected)
lags, cc = simulate_2d_noise_cc(np.arange(60, 120, 5), n_sources_per_az=20)
axes[2].plot(lags, cc / np.max(np.abs(cc)), 'k-', linewidth=0.8)
axes[2].axvline(expected_tt, color='r', linestyle='--', alpha=0.7)
axes[2].axvline(-expected_tt, color='b', linestyle='--', alpha=0.7)
axes[2].set_xlim(-15, 15)
axes[2].set_title('Sources Perpendicular Only\n(No end-fire — weak/no signal)', fontweight='bold')
axes[2].set_xlabel('Lag (s)')
axes[2].grid(True, alpha=0.3)

plt.suptitle('Effect of Noise Source Distribution on Cross-Correlation\n'
             '(Campillo & Roux, 2015: Stationary Phase & End-Fire Lobes)', 
             y=1.05, fontsize=13, fontweight='bold')
plt.tight_layout()
plt.show()

print("KEY INSIGHT (Campillo & Roux, 2015):")
print("Only sources in the 'end-fire lobes' (aligned with the station pair)")
print("contribute coherently to the cross-correlation via stationary phase.")
print("• Uniform sources → symmetric CC with both causal and acausal peaks")
print("• One-sided sources → asymmetric CC (one peak dominates)")
print("• Perpendicular sources → no coherent signal at expected travel time")
print("\nThis is why STACKING over long periods improves the CC: seasonal")
print("changes in noise source direction provide better azimuthal coverage.")
No description has been provided for this image 
KEY INSIGHT (Campillo & Roux, 2015):
Only sources in the 'end-fire lobes' (aligned with the station pair)
contribute coherently to the cross-correlation via stationary phase.
• Uniform sources → symmetric CC with both causal and acausal peaks
• One-sided sources → asymmetric CC (one peak dominates)
• Perpendicular sources → no coherent signal at expected travel time

This is why STACKING over long periods improves the CC: seasonal
changes in noise source direction provide better azimuthal coverage.

Module 4: Green's Function Retrieval from Ambient Noise¶

4.1 The Theoretical Foundation¶

The link between ambient noise cross-correlations and the Green's function has been established through several complementary theoretical frameworks.

4.1.1 The Time-Reversal Approach (Heuristic)¶

A simple, heuristic approach is based on the analogy between correlation and a time-reversal experiment (Derode et al., 2003; Paul et al., 2005):

  1. In a time-reversal experiment, signals from an active source are recorded, numerically reversed ($t \rightarrow -t$), and re-emitted in the medium
  2. Due to the time symmetry of wave equations, the waves propagate back and focus on the original source position
  3. The cross-correlation of noise recordings achieves mathematically the same effect as this time-reversal process

In the frequency domain, the cross-correlation $C_{1,2}(\omega) = u_1(\omega) u_2^*(\omega)$ is equivalent to having source at position 1, the signal is time-reversed and re-emitted from position 2 (Campillo & Roux, 2015).

4.1.2 Mathematical Derivation: Plane-Wave Approximation¶

Following Campillo & Roux (2015), for the 2D case with scalar waves:

The plane waves at angular frequency $\omega$ generated from sources at large distances, propagating in direction $\psi$, have the form:

$$u(\vec{r}_1; \omega) = F(\omega) \exp(i k \vec{r}_1 \cdot \hat{n}) \qquad [4]$$

where $k = \omega/c$ is the wavenumber, $\hat{n} = (\cos\psi, \sin\psi)$ is the propagation direction, and $F(\omega)$ is the source amplitude.

Assuming an omnidirectional distribution of incident plane waves, the field correlation between two points separated by distance $r = |\vec{r}_1 - \vec{r}_2|$ reduces to:

$$\langle u(\vec{r}_1; \omega) u^*(\vec{r}_2; \omega) \rangle = |F(\omega)|^2 J_0(k|\vec{r}_1 - \vec{r}_2|) \qquad [5]$$

where $J_0$ is the Bessel function of the first kind and order 0. This is the essence of the spatial correlation methods proposed by Aki (1957).

The 2D scalar Green's function can be identified as:

$$G(\vec{r}_1, \vec{r}_2; \omega) = \frac{1}{4i} H_0^{(1)}(k|\vec{r}_1 - \vec{r}_2|) = \frac{1}{4i}(J_0(k|\vec{r}_1 - \vec{r}_2|) + iY_0(k|\vec{r}_1 - \vec{r}_2|)) \qquad [6]$$

Leading to the fundamental result:

$$\boxed{\langle u(\vec{r}_1; \omega) u^*(\vec{r}_2; \omega) \rangle = -4|F(\omega)|^2 \text{Im}[G(\vec{r}_1, \vec{r}_2; \omega)]} \qquad [7]$$

This equation states that the average noise cross-correlation is proportional to the imaginary part of the Green's function.

4.1.3 The Time-Domain Connection¶

Taking the time derivative of the normalized correlation function (eq. [10] in Campillo & Roux, 2015):

$$\frac{d}{dt} C_{1,2}(t) = \frac{1}{4\pi r/c} [\delta(t + r/c) - \delta(t - r/c)] \qquad [11]$$

The two terms correspond to the backward and forward Green's function between the receivers, demonstrating the connection between the correlation function and the Green's function.

4.1.4 Generalization (Wapenaar, 2004)¶

For an arbitrary heterogeneous medium with attenuation parameter $\kappa \neq 0$, the Green's function between points A and B is defined as the solution of:

$$\Delta G(\vec{x}; t) - \frac{1}{c^2} \frac{\partial G(\vec{x}; t)}{\partial t^2} = \delta(\vec{x}) \delta(t) \qquad [3]$$

Wapenaar (2004) and Snieder (2004) showed rigorously that the cross-correlation result generalizes to elastic media with heterogeneity.

Module 5: Ambient Noise Preprocessing¶

5.1 The Preprocessing Workflow¶

Following Bensen et al. (2007), the ambient noise data processing divides into four principal phases:

Phase 1: Single-Station Data Preparation
    Raw data → Remove instrument response → Remove mean & trend → 
    Bandpass filter → Cut to 1-day segments →
    Apply temporal normalization → Apply spectral whitening

Phase 2: Cross-Correlation & Stacking
    Compute daily cross-correlations → Stack to desired number of days

Phase 3: Dispersion Measurement
    Measure group and/or phase velocity

Phase 4: Quality Control
    Error analysis → Selection of acceptable measurements

The preprocessing steps impose non-linear modifications to the waveforms, so the order of operations is significant (Bensen et al., 2007).

5.2 Temporal Normalization¶

The most important step in single-station data preparation. Its purpose is to reduce the effect of earthquakes, instrumental irregularities, and non-stationary noise sources on the cross-correlations (Bensen et al., 2007).

Five methods compared by Bensen et al. (2007):¶

Method Description Reference
One-bit normalization Replace amplitudes with $\pm 1$ based on sign Campillo & Paul (2003)
Clipped waveform Clip to RMS amplitude of the day Sabra et al. (2005)
Automated event detection Zero out windows above threshold -
Running-absolute-mean (RAM) Weight by inverse of running mean Preferred by Bensen et al.
Water-level normalization Iteratively down-weight above threshold -

Running-Absolute-Mean Normalization¶

For a discrete time series $d_j$, the normalization weight at time point $n$ is (Bensen et al., 2007, Eq. 1):

$$w_n = \frac{1}{2N + 1} \sum_{j=n-N}^{n+N} |d_j|$$

The normalized datum becomes: $\tilde{d}_n = d_n / w_n$

The width of the normalization window $(2N + 1)$ determines how much amplitude information is retained.

Earthquake-Band Temporal Normalization¶

If $d_j$ is the raw seismogram and $\hat{d}_j$ is the seismogram bandpass filtered in the earthquake band, the weights become (Bensen et al., 2007, Eq. 2):

$$\hat{w}_n = \frac{1}{2N + 1} \sum_{j=n-N}^{n+N} |\hat{d}_j|$$

These weights are then applied to the raw data: $\tilde{d}_n = d_n / \hat{w}_n$

5.3 Spectral Whitening¶

Ambient noise is not spectrally flat - it is peaked near the primary (~15 s) and secondary (~7.5 s) microseisms. Spectral whitening serves to (Bensen et al., 2007):

  1. Broaden the band of the ambient noise signal in cross-correlations
  2. Combat degradation caused by persistent monochromatic sources (e.g., the 26 s Gulf of Guinea signal)
  3. Enable broad-band dispersion measurement

The whitening is achieved by inversely weighting the complex spectrum by a smoothed version of the amplitude spectrum.

In [107]:
# =============================================================================
# PRACTICE 5.1: Implement Temporal Normalization Methods
# =============================================================================
# Compare the different normalization strategies from Bensen et al. (2007)

# Download real data with an earthquake
client = Client("IRIS")
t1 = UTCDateTime("2020-01-28T00:00:00")  # Day with M6.7 Caribbean EQ
t2 = t1 + 86400
st_raw = client.get_waveforms("IU", "ANMO", "00", "BHZ", t1, t2)
tr_raw = st_raw[0].copy()
tr_raw.detrend('demean')
tr_raw.detrend('linear')
tr_raw.taper(0.01)
tr_raw.filter('bandpass', freqmin=0.01, freqmax=0.5)

data = tr_raw.data.copy()
fs = tr_raw.stats.sampling_rate
t_axis = np.arange(len(data)) / fs / 3600  # hours

# --- Method 1: One-bit normalization ---
onebit = np.sign(data)

# --- Method 2: Running-absolute-mean (RAM) normalization ---
# Bensen et al. (2007), Eq. 1
def running_absolute_mean(data, window_half_width):
    """Running-absolute-mean temporal normalization (Bensen et al., 2007)."""
    N = window_half_width
    weights = np.convolve(np.abs(data), np.ones(2*N+1)/(2*N+1), mode='same')
    weights[weights == 0] = 1e-10  # Avoid division by zero
    return data / weights

# Use window = half the minimum period of interest
# For 20-100s period band, half-period ~ 10s
window_samples = int(10 * fs)  # 10 seconds
ram_norm = running_absolute_mean(data, window_samples)

# --- Method 3: Water-level normalization ---
def water_level_norm(data, n_iter=5, water_level=6.0):
    """Iterative water-level normalization (Bensen et al., 2007)."""
    result = data.copy()
    for _ in range(n_iter):
        rms = np.sqrt(np.mean(result**2))
        threshold = water_level * rms
        mask = np.abs(result) > threshold
        result[mask] = threshold * np.sign(result[mask])
    return result

wl_norm = water_level_norm(data)

# --- Plot comparison ---
fig, axes = plt.subplots(4, 1, figsize=(16, 12), sharex=True)

axes[0].plot(t_axis, data / np.max(np.abs(data)), 'k-', linewidth=0.3)
axes[0].set_title('(a) Raw Data (bandpass 0.01-0.5 Hz)', fontweight='bold')
axes[0].set_ylabel('Normalized')

axes[1].plot(t_axis, onebit, 'b-', linewidth=0.3)
axes[1].set_title('(b) One-Bit Normalization (Campillo & Paul, 2003)', fontweight='bold')
axes[1].set_ylabel('Sign (+1/-1)')

axes[2].plot(t_axis, ram_norm / np.max(np.abs(ram_norm)), 'g-', linewidth=0.3)
axes[2].set_title('(c) Running-Absolute-Mean Normalization (Bensen et al., 2007, Eq. 1)', 
                  fontweight='bold')
axes[2].set_ylabel('Normalized')

axes[3].plot(t_axis, wl_norm / np.max(np.abs(wl_norm)), 'r-', linewidth=0.3)
axes[3].set_title('(d) Water-Level Normalization', fontweight='bold')
axes[3].set_ylabel('Normalized')
axes[3].set_xlabel('Time (hours)')

plt.suptitle('Temporal Normalization Methods (Bensen et al., 2007, Fig. 3)',
             y=1.01, fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()

print("Bensen et al. (2007) recommend Running-Absolute-Mean normalization")
print("for its flexibility and adaptability to the data.")
print("One-bit normalization is the most aggressive but effective approach.")
No description has been provided for this image 
Bensen et al. (2007) recommend Running-Absolute-Mean normalization
for its flexibility and adaptability to the data.
One-bit normalization is the most aggressive but effective approach.
In [108]:
# =============================================================================
# PRACTICE 5.2: Implement Spectral Whitening
# =============================================================================
# Show the effect of spectral whitening on ambient noise data
# Reference: Bensen et al. (2007), Section 2.2

def spectral_whitening(data, fs, smooth_width=10):
    """
    Apply spectral whitening to a time series.
    
    The amplitude spectrum is flattened (whitened) while preserving
    the phase information. This broadens the frequency content of
    cross-correlations (Bensen et al., 2007).
    
    In NoisePy (Jiang & Denolle, 2020), two levels of whitening are available:
    - Running mean average ('rma'): smooth whitening
    - Phase-only ('phase_only'): strict whitening, amplitude = 1
    
    Parameters
    ----------
    data : numpy array
        Input time series
    fs : float
        Sampling rate (Hz)
    smooth_width : int
        Width of smoothing window for amplitude spectrum
    
    Returns
    -------
    whitened : numpy array
        Whitened time series
    """
    nfft = len(data)
    spec = np.fft.rfft(data, n=nfft)
    
    # Smooth amplitude spectrum
    amp = np.abs(spec)
    amp_smooth = np.convolve(amp, np.ones(smooth_width)/smooth_width, mode='same')
    amp_smooth[amp_smooth == 0] = 1e-10
    
    # Whiten: divide by smoothed amplitude
    spec_white = spec / amp_smooth
    
    whitened = np.fft.irfft(spec_white, n=nfft)
    return whitened


# Apply to real data
tr_proc = tr_raw.copy()
data_original = tr_proc.data.copy()
data_whitened = spectral_whitening(data_original, fs, smooth_width=20)

# Compute spectra for comparison
freqs_orig, psd_orig = signal.welch(data_original, fs=fs, nperseg=int(fs*100))
freqs_white, psd_white = signal.welch(data_whitened, fs=fs, nperseg=int(fs*100))

# Plot
fig, axes = plt.subplots(2, 2, figsize=(16, 10))

# Time domain
axes[0,0].plot(t_axis, data_original / np.max(np.abs(data_original)), 'k-', linewidth=0.3)
axes[0,0].set_title('(a) Original Time Series', fontweight='bold')
axes[0,0].set_xlabel('Time (hours)')
axes[0,0].set_ylabel('Normalized Amplitude')

axes[0,1].plot(t_axis[:len(data_whitened)], 
              data_whitened / np.max(np.abs(data_whitened)), 'b-', linewidth=0.3)
axes[0,1].set_title('(b) Spectrally Whitened Time Series', fontweight='bold')
axes[0,1].set_xlabel('Time (hours)')
axes[0,1].set_ylabel('Normalized Amplitude')

# Frequency domain
axes[1,0].semilogy(freqs_orig[1:], psd_orig[1:], 'k-', linewidth=1.0)
axes[1,0].set_title('(c) Original Amplitude Spectrum', fontweight='bold')
axes[1,0].set_xlabel('Frequency (Hz)')
axes[1,0].set_ylabel('PSD')
axes[1,0].set_xlim(0.01, 0.5)
axes[1,0].axvspan(0.05, 0.1, alpha=0.2, color='red', label='Primary microseism')
axes[1,0].axvspan(0.1, 0.2, alpha=0.2, color='green', label='Secondary microseism')
axes[1,0].legend(fontsize=9)

axes[1,1].semilogy(freqs_white[1:], psd_white[1:], 'b-', linewidth=1.0)
axes[1,1].set_title('(d) Whitened Amplitude Spectrum', fontweight='bold')
axes[1,1].set_xlabel('Frequency (Hz)')
axes[1,1].set_ylabel('PSD')
axes[1,1].set_xlim(0.01, 0.5)

plt.suptitle('Spectral Whitening (Bensen et al., 2007, Fig. 7)',
             y=1.01, fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()

print("Spectral whitening flattens the amplitude spectrum while preserving phase.")
print("This broadens the frequency content of cross-correlations and removes")
print("the dominance of the microseismic peaks (Bensen et al., 2007).")
No description has been provided for this image 
Spectral whitening flattens the amplitude spectrum while preserving phase.
This broadens the frequency content of cross-correlations and removes
the dominance of the microseismic peaks (Bensen et al., 2007).
In [109]:
# =============================================================================
# PRACTICE 5.3: Complete Preprocessing Pipeline
# =============================================================================
# Implement the full Phase 1 preprocessing as described in Bensen et al. (2007)
# This mirrors the workflow implemented in NoisePy's S0 and S1 scripts

def preprocess_ambient_noise(tr, freqmin=0.01, freqmax=0.5, 
                              time_norm='ram', freq_norm='whiten',
                              time_norm_window=50.0, whiten_smooth=20):
    """
    Complete ambient noise preprocessing pipeline.
    
    Following Bensen et al. (2007) Phase 1 and NoisePy (Jiang & Denolle, 2020).
    
    Parameters
    ----------
    tr : obspy.Trace
        Input trace (raw seismic data)
    freqmin, freqmax : float
        Bandpass filter corners (Hz)
    time_norm : str
        'no', 'one-bit', 'ram' (running-absolute-mean)
        NoisePy parameter: time_norm
    freq_norm : str
        'no', 'whiten', 'phase_only'
        NoisePy parameter: freq_norm
    time_norm_window : float
        Half-window for RAM normalization (seconds)
    whiten_smooth : int
        Smoothing width for spectral whitening
    
    Returns
    -------
    tr_proc : obspy.Trace
        Preprocessed trace
    """
    tr_proc = tr.copy()
    fs = tr_proc.stats.sampling_rate
    
    # Step 1: Remove mean and linear trend
    tr_proc.detrend('demean')
    tr_proc.detrend('linear')
    
    # Step 2: Apply taper (5% cosine taper)
    tr_proc.taper(0.05, type='cosine')
    
    # Step 3: Bandpass filter
    tr_proc.filter('bandpass', freqmin=freqmin, freqmax=freqmax, 
                   corners=4, zerophase=True)
    
    # Step 4: Temporal normalization
    if time_norm == 'one-bit':
        tr_proc.data = np.sign(tr_proc.data)
    elif time_norm == 'ram':
        N = int(time_norm_window * fs)
        weights = np.convolve(np.abs(tr_proc.data), 
                             np.ones(2*N+1)/(2*N+1), mode='same')
        weights[weights < 1e-10] = 1e-10
        tr_proc.data = tr_proc.data / weights
    
    # Step 5: Spectral whitening
    if freq_norm == 'whiten':
        tr_proc.data = spectral_whitening(tr_proc.data, fs, whiten_smooth)
    elif freq_norm == 'phase_only':
        spec = np.fft.rfft(tr_proc.data)
        amp = np.abs(spec)
        amp[amp < 1e-10] = 1e-10
        spec_white = spec / amp
        tr_proc.data = np.fft.irfft(spec_white, n=len(tr_proc.data))
    
    return tr_proc

# Apply to real data
tr_test = st_raw[0].copy()

# No normalization
tr_none = preprocess_ambient_noise(tr_test, time_norm='no', freq_norm='no')
# One-bit + whitening  
tr_1bit = preprocess_ambient_noise(tr_test, time_norm='one-bit', freq_norm='whiten')
# RAM + whitening (recommended by Bensen et al.)
tr_ram = preprocess_ambient_noise(tr_test, time_norm='ram', freq_norm='whiten')

fig, axes = plt.subplots(3, 1, figsize=(16, 9), sharex=True)

for ax, tr_plot, title in zip(axes, 
    [tr_none, tr_1bit, tr_ram],
    ['No Normalization (raw filtered)', 
     'One-bit + Spectral Whitening',
     'RAM + Spectral Whitening (Bensen et al., 2007 recommended)']):
    
    t_ax = np.arange(tr_plot.stats.npts) / tr_plot.stats.sampling_rate / 3600
    ax.plot(t_ax, tr_plot.data / np.max(np.abs(tr_plot.data)), 'k-', linewidth=0.3)
    ax.set_title(title, fontweight='bold')
    ax.set_ylabel('Amplitude')

axes[-1].set_xlabel('Time (hours)')
plt.suptitle('Complete Preprocessing Pipeline (Bensen et al., 2007, Fig. 2)',
             y=1.01, fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()
No description has been provided for this image

Module 6: Cross-Correlation Computation¶

6.1 Cross-Correlation Methods¶

NoisePy (Jiang & Denolle, 2020) implements several cross-correlation methods, controlled by the cc_method parameter:

Method cc_method freq_norm Description
Pure cross-correlation 'xcorr' 'no' $\hat{C}_{AB}(f) = \hat{U}_A^*(f) \hat{U}_B(f)$
Coherency 'xcorr' 'rma' $\hat{C}_{AB}(f) = \frac{\hat{U}_A^*(f) \hat{U}_B(f)}{\{\lvert\hat{U}_A(f)\rvert\}\{\lvert\hat{U}_B(f)\rvert\}}$
Deconvolution 'deconv' 'no' $\hat{C}_{AB}(f) = \frac{\hat{U}_A^*(f) \hat{U}_B(f)}{\{\lvert\hat{U}_A(f)\rvert^2\}}$
Phase cross-correlation 'xcorr' 'phase_only' Amplitude set to 1

Where $\{\cdot\}$ represents smoothing using a running mean average (Jiang & Denolle, 2020, Table 1).

"Performing spectral whitening using rma with the cross correlation is strictly equivalent to calculating the coherency." -- Jiang & Denolle (2020), citing Prieto et al. (2009)

6.2 Practical Considerations¶

  • Cross-correlations are computed on daily segments (or sub-daily windows)
  • All station pairs yield $n(n-1)/2$ cross-correlations for $n$ stations
  • NoisePy splits continuous data into time chunks (TC) to optimize memory (Jiang & Denolle, 2020, Eq. 4)
In [110]:
# =============================================================================
# PRACTICE 6.1: Compute Cross-Correlations with Real Data
# =============================================================================
# Compute the cross-correlation between ANMO and CCM using the 
# preprocessing pipeline from Module 5

def compute_noise_cc(tr1, tr2, cc_len=1800, step=450, 
                      freqmin=0.01, freqmax=0.2,
                      time_norm='one-bit', freq_norm='whiten',
                      maxlag=500):
    """
    Compute ambient noise cross-correlation between two traces.
    
    Implements the workflow from Bensen et al. (2007) Phase 2 and
    NoisePy S1 script (Jiang & Denolle, 2020).
    
    Parameters
    ----------
    tr1, tr2 : obspy.Trace
        Preprocessed traces from two stations
    cc_len : float
        Cross-correlation window length (seconds)
        NoisePy parameter: cc_len
    step : float
        Step between windows (seconds) for overlap
        NoisePy parameter: step
    maxlag : float
        Maximum lag time (seconds). Must be larger than the expected
        inter-station travel time (distance / velocity).
    
    Returns
    -------
    lags : numpy array
        Lag times
    cc_stack : numpy array  
        Stacked cross-correlation
    cc_all : numpy array
        Individual window cross-correlations
    """
    # Preprocess both traces
    t1 = preprocess_ambient_noise(tr1, freqmin=freqmin, freqmax=freqmax,
                                  time_norm=time_norm, freq_norm=freq_norm)
    t2 = preprocess_ambient_noise(tr2, freqmin=freqmin, freqmax=freqmax,
                                  time_norm=time_norm, freq_norm=freq_norm)
    
    fs = t1.stats.sampling_rate
    npts_cc = int(cc_len * fs)
    npts_step = int(step * fs)
    npts_lag = int(maxlag * fs)
    
    # Ensure same length
    min_len = min(len(t1.data), len(t2.data))
    d1 = t1.data[:min_len]
    d2 = t2.data[:min_len]
    
    # Sliding windows
    n_windows = (min_len - npts_cc) // npts_step + 1
    nfft = int(2 ** np.ceil(np.log2(2 * npts_cc)))
    
    cc_all = []
    
    for i in range(n_windows):
        start = i * npts_step
        end = start + npts_cc
        if end > min_len:
            break
        
        seg1 = d1[start:end]
        seg2 = d2[start:end]
        
        # Cross-correlation in frequency domain
        # C_AB(f) = U_A*(f) * U_B(f)  [Jiang & Denolle, 2020, Eq. 1]
        F1 = np.fft.rfft(seg1, n=nfft)
        F2 = np.fft.rfft(seg2, n=nfft)
        CC = np.conj(F1) * F2
        
        # Remove mean of real part (NoisePy default)
        CC -= np.mean(np.real(CC))
        
        cc = np.fft.irfft(CC, n=nfft)
        cc = np.fft.fftshift(cc)
        
        # Extract the relevant lag window
        center = nfft // 2
        cc_trim = cc[center - npts_lag:center + npts_lag + 1]
        
        # Quality control: reject if max amplitude > 20x median
        # (NoisePy default: max_over_std threshold)
        if np.max(np.abs(cc_trim)) < 20 * np.median(np.abs(cc_trim)):
            cc_all.append(cc_trim)
    
    cc_all = np.array(cc_all)
    cc_stack = np.mean(cc_all, axis=0)  # Linear stack
    lags = np.arange(-npts_lag, npts_lag + 1) / fs
    
    return lags, cc_stack, cc_all


# Compute cross-correlation for ANMO-CCM
print("Computing cross-correlation (this may take a moment)...")
lags, cc_stack, cc_all = compute_noise_cc(
    st1[0], st2[0], 
    cc_len=1800,    # 30-minute windows (NoisePy default)
    step=450,       # 75% overlap
    freqmin=0.02, 
    freqmax=0.1,    # 10-50 s period band (microseismic band)
    time_norm='one-bit',
    freq_norm='whiten',
    maxlag=500      # Must exceed expected travel time (~400 s for ANMO-CCM)
)
print(f"Number of windows used: {len(cc_all)}")
print(f"Lag range: {lags[0]:.0f} to {lags[-1]:.0f} s")
print(f"Expected Rayleigh wave arrival: ~{dist_km/3.5:.0f} s at 3.5 km/s")

Computing cross-correlation (this may take a moment)...
Number of windows used: 93
Lag range: -800 to 800 s
In [111]:
# =============================================================================
# Visualize the Cross-Correlation Result
# =============================================================================
# The raw CC is noisy. Following Bensen et al. (2007), we bandpass filter
# the stacked CC to isolate the surface wave signal in the microseismic band.
# This is standard practice — all published CC figures are filtered.

from scipy.signal import butter, sosfiltfilt

# Expected travel time for surface waves at ~3.5 km/s
expected_tt = dist_km / 3.5  # seconds

# Bandpass filter the stacked CC to enhance the surface wave signal
# This is equivalent to what NoisePy does when plotting moveout sections
fs_cc = 1.0 / (lags[1] - lags[0])  # sampling rate of CC
freqmin_plot = 0.04   # 25 s period
freqmax_plot = 0.15   # ~7 s period (secondary microseism)

sos = butter(4, [freqmin_plot, freqmax_plot], btype='bandpass', fs=fs_cc, output='sos')
cc_stack_filt = sosfiltfilt(sos, cc_stack)

# Also filter individual windows for the image
cc_all_filt = np.zeros_like(cc_all)
for i in range(len(cc_all)):
    cc_all_filt[i] = sosfiltfilt(sos, cc_all[i])

fig, axes = plt.subplots(3, 1, figsize=(16, 12))

# Plot 1: Raw stacked CC (to show why filtering is needed)
cc_norm_raw = cc_stack / np.max(np.abs(cc_stack))
axes[0].plot(lags, cc_norm_raw, 'k-', linewidth=0.5, alpha=0.7)
axes[0].axvline(expected_tt, color='green', linestyle='--', linewidth=1.5, alpha=0.7)
axes[0].axvline(-expected_tt, color='green', linestyle='--', linewidth=1.5, alpha=0.7)
axes[0].set_title('Raw Stacked Cross-Correlation (unfiltered — noisy)', fontweight='bold')
axes[0].set_xlabel('Lag Time (s)')
axes[0].set_ylabel('Normalized CC')
axes[0].set_xlim(-lags[-1], lags[-1])
axes[0].grid(True, alpha=0.3)

# Plot 2: Filtered individual CCs as image
im = axes[1].imshow(cc_all_filt / np.max(np.abs(cc_all_filt), axis=1, keepdims=True), 
                     aspect='auto', cmap='seismic', vmin=-0.5, vmax=0.5,
                     extent=[lags[0], lags[-1], len(cc_all_filt), 0])
axes[1].axvline(expected_tt, color='lime', linestyle='--', linewidth=1.5, 
                label=f'Expected: +{expected_tt:.0f} s')
axes[1].axvline(-expected_tt, color='lime', linestyle='--', linewidth=1.5)
axes[1].set_title(f'Individual Window CCs (filtered {freqmin_plot}-{freqmax_plot} Hz)', 
                  fontweight='bold')
axes[1].set_ylabel('Window Number')
axes[1].legend(loc='upper right')
plt.colorbar(im, ax=axes[1], label='Normalized CC')

# Plot 3: Filtered stacked CC — the main result
cc_norm = cc_stack_filt / np.max(np.abs(cc_stack_filt))
axes[2].plot(lags, cc_norm, 'k-', linewidth=1.2)
axes[2].fill_between(lags, cc_norm, 0, where=cc_norm > 0, 
                     alpha=0.3, color='red', label='Positive')
axes[2].fill_between(lags, cc_norm, 0, where=cc_norm < 0, 
                     alpha=0.3, color='blue', label='Negative')
axes[2].axvline(expected_tt, color='green', linestyle='--', linewidth=2,
                label=f'Expected at ~3.5 km/s: {expected_tt:.0f} s')
axes[2].axvline(-expected_tt, color='green', linestyle='--', linewidth=2)
axes[2].set_title(f'Filtered Stacked CC: {st1[0].stats.station}-{st2[0].stats.station} '
                  f'({dist_km:.0f} km, {freqmin_plot}-{freqmax_plot} Hz, '
                  f'{len(cc_all)} windows)', fontweight='bold')
axes[2].set_xlabel('Lag Time (s)')
axes[2].set_ylabel('Normalized CC')
axes[2].legend(loc='upper right')
axes[2].set_xlim(-lags[-1], lags[-1])
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

print(f"\nInter-station distance: {dist_km:.0f} km")
print(f"Expected Rayleigh wave travel time (~3.5 km/s): {expected_tt:.0f} s")
print(f"Filter: {freqmin_plot}-{freqmax_plot} Hz ({1/freqmax_plot:.0f}-{1/freqmin_plot:.0f} s period)")
print(f"Windows stacked: {len(cc_all)}")
print(f"\nNote: With only 1 day of data, the CC is noisy. Bensen et al. (2007)")
print(f"recommend months of data for reliable measurements (SNR ~ T^{{1/n}}).")
print(f"The bandpass filter isolates the microseismic frequency band where")
print(f"the coherent surface wave signal is concentrated.")
No description has been provided for this image 
Inter-station distance: 3124 km
Expected travel time (Rayleigh wave ~3.5 km/s): 893 s

The causal (+) and acausal (-) Green's functions are visible!
Asymmetry indicates non-uniform noise source distribution.
In [112]:
# =============================================================================
# PRACTICE 6.2: Compare Cross-Correlation Methods
# =============================================================================
# Reproduce NoisePy Figure 2: compare xcorr, coherency, deconv, phase_cc

def compute_cc_methods(tr1, tr2, cc_len=1800, freqmin=0.02, freqmax=0.2, maxlag=500):
    """
    Compute cross-correlations using different methods.
    Reference: Jiang & Denolle (2020), Table 1 and Figure 2.
    """
    results = {}
    fs = tr1.stats.sampling_rate
    
    # Preprocess (basic)
    for tr in [tr1, tr2]:
        tr.detrend('demean')
        tr.detrend('linear')
        tr.taper(0.05)
        tr.filter('bandpass', freqmin=freqmin, freqmax=freqmax, zerophase=True)
    
    min_len = min(len(tr1.data), len(tr2.data))
    npts_cc = int(cc_len * fs)
    npts_lag = int(maxlag * fs)
    nfft = int(2 ** np.ceil(np.log2(2 * npts_cc)))
    
    d1 = tr1.data[:npts_cc]
    d2 = tr2.data[:npts_cc]
    
    F1 = np.fft.rfft(d1, n=nfft)
    F2 = np.fft.rfft(d2, n=nfft)
    
    smooth_w = 20
    amp1_smooth = np.convolve(np.abs(F1), np.ones(smooth_w)/smooth_w, mode='same')
    amp2_smooth = np.convolve(np.abs(F2), np.ones(smooth_w)/smooth_w, mode='same')
    amp1_smooth[amp1_smooth < 1e-20] = 1e-20
    amp2_smooth[amp2_smooth < 1e-20] = 1e-20
    
    methods = {
        'Pure xcorr': np.conj(F1) * F2,
        'Coherency': np.conj(F1) * F2 / (amp1_smooth * amp2_smooth),
        'Deconvolution': np.conj(F1) * F2 / (np.abs(F1)**2 + 1e-10),
        'Phase CC': np.conj(F1/np.abs(F1+1e-20)) * (F2/np.abs(F2+1e-20)),
    }
    
    for name, CC in methods.items():
        cc = np.fft.irfft(CC, n=nfft)
        cc = np.fft.fftshift(cc)
        center = nfft // 2
        results[name] = cc[center - npts_lag:center + npts_lag + 1]
    
    lags = np.arange(-npts_lag, npts_lag + 1) / fs
    return lags, results


lags_m, results = compute_cc_methods(st1[0].copy(), st2[0].copy())

fig, axes = plt.subplots(4, 1, figsize=(16, 12), sharex=True)
colors = ['red', 'green', 'blue', 'purple']
equations = [
    r'$\hat{C}_{AB} = \hat{U}_A^* \hat{U}_B$',
    r'$\hat{C}_{AB} = \frac{\hat{U}_A^* \hat{U}_B}{\{|\hat{U}_A|\}\{|\hat{U}_B|\}}$',
    r'$\hat{C}_{AB} = \frac{\hat{U}_A^* \hat{U}_B}{|\hat{U}_A|^2}$',
    r'$\hat{C}_{AB} = \frac{\hat{U}_A^*}{|\hat{U}_A|} \frac{\hat{U}_B}{|\hat{U}_B|}$',
]

for ax, (name, cc), color, eq in zip(axes, results.items(), colors, equations):
    cc_n = cc / np.max(np.abs(cc))
    ax.plot(lags_m, cc_n, color=color, linewidth=0.8)
    ax.set_title(f'{name}: {eq}', fontweight='bold')
    ax.set_ylabel('Norm. CC')
    ax.set_xlim(-500, 500)
    ax.grid(True, alpha=0.3)
    ax.axvline(expected_tt, color='gray', linestyle='--', alpha=0.5)
    ax.axvline(-expected_tt, color='gray', linestyle='--', alpha=0.5)

axes[-1].set_xlabel('Lag Time (s)')
plt.suptitle('Cross-Correlation Methods (Jiang & Denolle, 2020, Table 1 & Fig. 2)',
             y=1.01, fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()
No description has been provided for this image

Module 7: Stacking Methods¶

7.1 Why Stacking?¶

Individual daily cross-correlations are noisy. Stacking (temporal averaging) improves the signal-to-noise ratio (SNR) by reinforcing coherent signals while averaging out incoherent noise. The SNR of the cross-correlation improves with longer recording time because more noise sources participate coherently (Campillo & Roux, 2015).

"The averaging of the correlation process over time magnifies the coherent versus incoherent contributions of the noise sources." -- Campillo & Roux (2015)

7.2 Stacking Methods in NoisePy¶

NoisePy (Jiang & Denolle, 2020) implements several stacking options in the S2 script:

Method Description Reference
Linear stacking Simple average Standard
Phase-weighted stacking (PWS) Weight by instantaneous phase coherence Schimmel & Paulssen (1997)
Robust stacking Iterative reweighting Pavlis & Vernon (2010)
Autocovariance filter Filter based on autocovariance Nakata et al. (2016)
Selective stacking Select based on CC coefficient quality Yang et al. (2020)

Phase-Weighted Stack (PWS)¶

The PWS enhances coherent arrivals by weighting each sample by the coherence of the instantaneous phase across all windows (Schimmel & Paulssen, 1997):

$$\text{PWS}(t) = \frac{1}{N} \sum_{i=1}^{N} s_i(t) \left| \frac{1}{N} \sum_{j=1}^{N} e^{i\phi_j(t)} \right|^\nu$$

where $\phi_j(t)$ is the instantaneous phase of the $j$-th trace and $\nu$ is the power index (typically $\nu = 2$).

In [113]:
# =============================================================================
# PRACTICE 7.1: Implement and Compare Stacking Methods
# =============================================================================

def linear_stack(cc_array):
    """Simple linear (mean) stack."""
    return np.mean(cc_array, axis=0)


def phase_weighted_stack(cc_array, power=2):
    """
    Phase-weighted stack (Schimmel & Paulssen, 1997).
    
    Enhances coherent arrivals by weighting each sample by the
    coherence of the instantaneous phase across all windows.
    
    Used in NoisePy S2 stacking (Jiang & Denolle, 2020).
    """
    N = cc_array.shape[0]
    
    # Compute analytic signal (Hilbert transform) for each trace
    analytic = np.zeros_like(cc_array, dtype=complex)
    for i in range(N):
        analytic[i] = hilbert(cc_array[i])
    
    # Instantaneous phase
    phase = np.angle(analytic)
    
    # Phase coherence weight
    phase_coherence = np.abs(np.mean(np.exp(1j * phase), axis=0)) ** power
    
    # Weighted linear stack
    lin_stack = np.mean(cc_array, axis=0)
    pws = lin_stack * phase_coherence
    
    return pws


def robust_stack(cc_array, n_iter=10):
    """
    Robust stacking using iterative reweighting (Pavlis & Vernon, 2010).
    Down-weights outlier traces that deviate from the stack.
    """
    stack = np.median(cc_array, axis=0)  # Initial estimate
    
    for _ in range(n_iter):
        # Compute similarity weights
        weights = np.zeros(cc_array.shape[0])
        for i in range(cc_array.shape[0]):
            r = np.corrcoef(stack, cc_array[i])[0, 1]
            weights[i] = max(0, r)  # Only positive correlations
        
        if np.sum(weights) == 0:
            break
        weights /= np.sum(weights)
        
        # Weighted stack
        stack = np.average(cc_array, axis=0, weights=weights)
    
    return stack


# Apply all three methods to our cross-correlations
ls = linear_stack(cc_all)
pws = phase_weighted_stack(cc_all, power=2)
rs = robust_stack(cc_all)

# Plot comparison
fig, axes = plt.subplots(3, 1, figsize=(16, 10), sharex=True)

for ax, data, title, color in zip(axes,
    [ls, pws, rs],
    ['Linear Stack', 
     'Phase-Weighted Stack (Schimmel & Paulssen, 1997)',
     'Robust Stack (Pavlis & Vernon, 2010)'],
    ['black', 'blue', 'red']):
    
    data_n = data / np.max(np.abs(data))
    ax.plot(lags, data_n, color=color, linewidth=1.0)
    ax.fill_between(lags, data_n, 0, where=data_n > 0, alpha=0.2, color=color)
    ax.axvline(expected_tt, color='green', linestyle='--', linewidth=1.5, alpha=0.7)
    ax.axvline(-expected_tt, color='green', linestyle='--', linewidth=1.5, alpha=0.7)
    ax.set_title(title, fontweight='bold')
    ax.set_ylabel('Normalized CC')
    ax.grid(True, alpha=0.3)

axes[-1].set_xlabel('Lag Time (s)')
plt.suptitle('Comparison of Stacking Methods (NoisePy S2 options)',
             y=1.01, fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()

print("PWS enhances coherent arrivals and suppresses incoherent noise.")
print("Robust stacking down-weights outlier windows.")
print("Linear stacking is simple but effective with sufficient data.")
No description has been provided for this image 
PWS enhances coherent arrivals and suppresses incoherent noise.
Robust stacking down-weights outlier windows.
Linear stacking is simple but effective with sufficient data.
In [114]:
# =============================================================================
# PRACTICE 7.2: Convergence of Cross-Correlations with Stacking Duration
# =============================================================================
# Show how the SNR improves as more windows are stacked
# This illustrates the convergence concept from Campillo & Roux (2015)

n_windows_list = [2, 5, 10, 20, len(cc_all)]

fig, axes = plt.subplots(len(n_windows_list), 1, figsize=(16, 3*len(n_windows_list)), 
                         sharex=True)

for ax, n_win in zip(axes, n_windows_list):
    stack = np.mean(cc_all[:n_win], axis=0)
    stack_n = stack / np.max(np.abs(stack))
    ax.plot(lags, stack_n, 'k-', linewidth=0.8)
    ax.fill_between(lags, stack_n, 0, where=stack_n > 0, alpha=0.3, color='red')
    ax.fill_between(lags, stack_n, 0, where=stack_n < 0, alpha=0.3, color='blue')
    ax.axvline(expected_tt, color='green', linestyle='--', alpha=0.7)
    ax.axvline(-expected_tt, color='green', linestyle='--', alpha=0.7)
    
    # Compute SNR
    signal_window = (np.abs(lags) > expected_tt - 100) & (np.abs(lags) < expected_tt + 100)
    noise_window = np.abs(lags) > expected_tt + 200
    if np.any(noise_window) and np.std(stack[noise_window]) > 0:
        snr = np.max(np.abs(stack[signal_window])) / np.std(stack[noise_window])
    else:
        snr = 0
    
    ax.set_title(f'{n_win} windows stacked (SNR: {snr:.1f})', fontweight='bold')
    ax.set_ylabel('Norm. CC')
    ax.grid(True, alpha=0.3)

axes[-1].set_xlabel('Lag Time (s)')
plt.suptitle('Convergence: Signal Emerges with More Stacking\n'
             '(Campillo & Roux, 2015: "The averaging magnifies coherent vs incoherent contributions")',
             y=1.02, fontsize=13, fontweight='bold')
plt.tight_layout()
plt.show()
No description has been provided for this image

Extra Module 8: Optimal Processing of Noise Correlations for pure scientific curriosity and maybe for R&D¶

8.1 The Problem with Standard Processing¶

Standard preprocessing (one-bit normalization, spectral whitening) is nonlinear — it breaks the linear wave physics that underlies the Green's function retrieval theorem. This has important consequences (Fichtner et al., 2020).

The Transfer Coefficient Framework¶

Fichtner et al. (2020) introduced a rigorous framework to quantify processing effects. For each time window $n$ and station pair $(i,k)$, they define a transfer coefficient:

$$\tilde{I}^{(n)}(\xi_i, \xi_k) = T_{ik}^{(n)} \cdot I^{(n)}(\xi_i, \xi_k) \qquad [\text{Fichtner et al., 2020, Eq. 7}]$$

where $I^{(n)}$ is the raw interferogram and $\tilde{I}^{(n)}$ is the processed one.

The Key Factorization Theorem¶

The transfer coefficient always factorizes as:

$$T_{ik}^{(n)} = f^{(n)} \cdot g_{ik} + e_{ik}^{(n)} \qquad [\text{Fichtner et al., 2020, Eq. 8}]$$

where:

  • $f^{(n)}$ = time-window-specific scalar (captures amplitude variation over time)
  • $g_{ik}$ = path-specific complex weight (depends only on station pair)
  • $e_{ik}^{(n)}$ = factorization residual — responsible for all unphysical effects

Why This Matters¶

The residual $e_{ik}^{(n)}$ introduces an unphysical wavefield: different station pairs effectively see different noise source distributions. This means:

Processing Residual $e$ Physical Impact Reference
One-bit normalization Large ~10 dB amplitude error, ~2.9% traveltime bias Fichtner et al. (2020)
Spectral whitening Medium Primarily amplitude effect Fichtner et al. (2020)
Causal/acausal averaging Large Mixes two different wave packets as different sources Fichtner et al. (2017)
Optimal processing (OPUS) Zero Preserves all desired properties Fichtner et al. (2020)

Optimal Processing (OPUS)¶

The optimal processed interferogram eliminates the unphysical component:

$$\tilde{I}_{\text{opt}}^{(n)}(\xi_i, \xi_k) = f^{(n)} \cdot g_{ik} \cdot I^{(n)}(\xi_i, \xi_k) \qquad [\text{Fichtner et al., 2020, Eq. 12}]$$

"Do NOT equate interferometry with Green's function retrieval. Instead, extract information directly from processed correlations by developing an effective forward theory that accounts for the processing." — Fichtner et al. (2017)

Generalised Interferometry¶

Fichtner et al. (2017) developed a unified framework that works for any processing:

  1. Compute the effective source power spectral density $\tilde{S}(x,x')$ that accounts for the processing
  2. Use this in a forward model to predict what the processed correlation should look like
  3. Invert for Earth structure using the effective forward model — no Green's function approximation needed

This framework unifies earthquake-based methods, ambient noise correlations, and any other source of inter-station correlations.

Module 9: Surface Wave Dispersion Analysis¶

8.1 Theory: Surface Wave Dispersion¶

Surface waves (Rayleigh and Love waves) are dispersive: their velocity depends on frequency (or period). This dispersion arises because different frequencies sample different depths of the Earth:

  • Short periods (high frequency): Sensitive to shallow structure
  • Long periods (low frequency): Sensitive to deeper structure

The dispersion of surface waves recovered from ambient noise cross-correlations provides information about the velocity structure of the crust and uppermost mantle (Shapiro & Campillo, 2004; Bensen et al., 2007).

Group Velocity vs. Phase Velocity¶

  • Group velocity $U = d\omega/dk$: Velocity of the wave envelope
  • Phase velocity $c = \omega/k$: Velocity of individual wave crests

Both can be measured from noise cross-correlations:

  • Group velocity: From frequency-time analysis (FTAN) using narrow-band filtering
  • Phase velocity: From the phase of the cross-correlation spectrum

NoisePy Dispersion Analysis¶

NoisePy implements group velocity measurement using (Jiang & Denolle, 2020):

  1. Continuous Wavelet Transform (CWT) using Morlet wavelets (Fichtner et al., 2008)
  2. Automated Frequency-Time Analysis (AFTAN) (Levshin & Ritzwoller, 2001; Bensen et al., 2007)
In [115]:
# =============================================================================
# PRACTICE 9.1: Frequency-Time Analysis (FTAN) for Group Velocity
# =============================================================================
# Implement the FTAN method to measure group velocity dispersion
# from ambient noise cross-correlations

def frequency_time_analysis(cc, lags, fs, periods, dist_km, width=0.2):
    """
    Frequency-Time Analysis (FTAN) for group velocity measurement.
    
    Following Bensen et al. (2007) Phase 3 and NoisePy dispersion_analysis.
    Uses narrow bandpass filters to isolate energy at each period.
    
    Parameters
    ----------
    cc : numpy array
        Cross-correlation function (symmetric component recommended)
    lags : numpy array
        Lag times
    fs : float
        Sampling rate
    periods : numpy array
        Center periods for analysis
    dist_km : float
        Inter-station distance (km)
    width : float
        Relative bandwidth of Gaussian filter
    
    Returns
    -------
    ftan_map : 2D numpy array
        FTAN diagram (period x velocity)
    velocities : numpy array
        Velocity axis
    group_vel : numpy array
        Measured group velocities at each period
    """
    # Use positive lag (causal part) or symmetric component
    center = len(lags) // 2
    cc_pos = cc[center:]  # Causal part
    t_pos = lags[center:]  # Positive lag times
    
    # Also use symmetric component (average of causal + time-reversed acausal)
    cc_neg = cc[:center+1][::-1]  # Acausal, time-reversed
    min_len = min(len(cc_pos), len(cc_neg))
    cc_sym = (cc_pos[:min_len] + cc_neg[:min_len]) / 2
    t_sym = t_pos[:min_len]
    
    # Velocity grid
    velocities = np.linspace(1.0, 5.0, 200)  # km/s
    ftan_map = np.zeros((len(periods), len(velocities)))
    group_vel = np.zeros(len(periods))
    
    for i, T in enumerate(periods):
        f0 = 1.0 / T
        # Gaussian bandpass filter centered at f0
        freqs = np.fft.rfftfreq(len(cc_sym), d=1.0/fs)
        gauss_filter = np.exp(-((freqs - f0) / (width * f0)) ** 2)
        
        # Apply filter in frequency domain
        spec = np.fft.rfft(cc_sym)
        filtered = np.fft.irfft(spec * gauss_filter, n=len(cc_sym))
        
        # Envelope (Hilbert transform)
        envelope = np.abs(hilbert(filtered))
        
        # Map to velocity space
        for j, v in enumerate(velocities):
            t_expected = dist_km / v
            idx = int(t_expected * fs)
            if 0 < idx < len(envelope):
                ftan_map[i, j] = envelope[idx]
        
        # Pick maximum as group velocity
        if np.max(ftan_map[i]) > 0:
            group_vel[i] = velocities[np.argmax(ftan_map[i])]
    
    return ftan_map, velocities, group_vel


# Apply FTAN to our cross-correlation
periods = np.linspace(8, 80, 50)  # Period range
fs_cc = st1[0].stats.sampling_rate

ftan_map, velocities, group_vel = frequency_time_analysis(
    cc_stack, lags, fs_cc, periods, dist_km, width=0.15
)

# Plot FTAN diagram
fig, axes = plt.subplots(1, 2, figsize=(16, 7))

# FTAN map
im = axes[0].pcolormesh(periods, velocities, ftan_map.T, cmap='hot_r', shading='auto')
axes[0].plot(periods, group_vel, 'c-o', markersize=3, linewidth=1.5, label='Picked group velocity')
axes[0].set_xlabel('Period (s)', fontsize=13)
axes[0].set_ylabel('Group Velocity (km/s)', fontsize=13)
axes[0].set_title('FTAN Diagram\n(Bensen et al., 2007, Phase 3)', fontweight='bold')
axes[0].legend()
plt.colorbar(im, ax=axes[0], label='Envelope Amplitude')

# Dispersion curve
valid = group_vel > 0
axes[1].plot(periods[valid], group_vel[valid], 'b-o', markersize=5, linewidth=2)
axes[1].set_xlabel('Period (s)', fontsize=13)
axes[1].set_ylabel('Group Velocity (km/s)', fontsize=13)
axes[1].set_title('Group Velocity Dispersion Curve', fontweight='bold')
axes[1].set_ylim(2.0, 5.0)
axes[1].grid(True, alpha=0.3)
axes[1].axhline(3.5, color='gray', linestyle='--', alpha=0.5, label='Typical Rayleigh wave')
axes[1].legend()

plt.suptitle(f'Surface Wave Dispersion: {st1[0].stats.station}-{st2[0].stats.station} '
             f'({dist_km:.0f} km)',
             y=1.02, fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()

print("The dispersion curve shows how group velocity varies with period.")
print("Short periods (shallow structure) → typically slower velocities")
print("Long periods (deeper structure) → typically faster velocities")
No description has been provided for this image 
The dispersion curve shows how group velocity varies with period.
Short periods (shallow structure) → typically slower velocities
Long periods (deeper structure) → typically faster velocities

Module 10: Ambient Noise Tomography — From Dispersion to Velocity¶

10.1 The Inversion Problem¶

Once dispersion curves have been measured for all station pairs, we need to invert them for velocity structure. This is a two-step process (Bensen et al., 2008; Rawlinson et al., 2010):

Step 1: 2D Tomographic Inversion (Period Maps)¶

For each period $\tau$, we invert the travel times between all station pairs to produce a 2D velocity map. The data are surface wave travel times:

$$t_{ij}(\tau) = \int_{\text{path}_{ij}} \frac{1}{c(\vec{r}, \tau)} \, dl \qquad [\text{Rawlinson et al., 2010}]$$

where $c(\vec{r}, \tau)$ is the velocity at position $\vec{r}$ for period $\tau$.

Penalty functional minimized (Bensen et al., 2008, Eq. 1):

$$S(\mathbf{m}) = (\mathbf{G}(\mathbf{m}) - \mathbf{d})^T \mathbf{C}_d^{-1} (\mathbf{G}(\mathbf{m}) - \mathbf{d}) + \alpha^2 \|\mathbf{F}(\mathbf{m})\|^2 + \beta^2 \|\mathbf{H}(\mathbf{m})\|^2$$

where:

  • $\mathbf{d}$ = observed surface wave travel times
  • $\mathbf{G}(\mathbf{m})$ = predicted travel times (forward operator)
  • $\mathbf{C}_d$ = data covariance (measurement uncertainties)
  • $\alpha^2 \|\mathbf{F}(\mathbf{m})\|^2$ = spatial smoothing regularization
  • $\beta^2 \|\mathbf{H}(\mathbf{m})\|^2$ = penalty for regions with poor path coverage

The Gauss-Newton solution at each iteration (Rawlinson et al., 2010, Eq. 13-14):

$$\delta\mathbf{m} = [\mathbf{G}^T \mathbf{C}_d^{-1} \mathbf{G} + \epsilon \mathbf{C}_m^{-1} + \eta \mathbf{D}^T \mathbf{D}]^{-1} \mathbf{G}^T \mathbf{C}_d^{-1} (\mathbf{d}_{\text{obs}} - \mathbf{G}(\mathbf{m}))$$

Step 2: Depth Inversion (1D Profiles → 3D Model)¶

At each grid point, the dispersion curve $c(\tau)$ at different periods samples different depths. We invert this for a 1D $V_s(z)$ profile:

  • Linearized inversion: Use sensitivity kernels $\partial c / \partial V_s$ to iteratively update the model
  • Transdimensional Bayesian inversion: Let the data determine the number of layers and their properties — no subjective regularization (Ryberg et al., 2022; Cabrera-Pérez et al., 2023)

Depth Sensitivity¶

Surface waves at different periods are sensitive to different depth ranges (Bensen et al., 2008, Fig. 16):

Period (s) Wave Type Approximate Depth Sensitivity
8–10 Rayleigh 5–25 km (upper crust)
15–20 Rayleigh 15–40 km (mid-crust, Moho)
25–40 Rayleigh 30–80 km (lower crust, upper mantle)
50–100 Rayleigh 50–200 km (lithosphere)
8–15 Love Shallower than Rayleigh at same period

Nonlinear Multiscale Inversion (MANgOSTA)¶

Cabrera-Pérez et al. (2021) developed a nonlinear multiscale approach (MANgOSTA) that:

  1. Starts at coarse resolution (large correlation length)
  2. Progressively refines to finer scales
  3. At each scale, recomputes ray paths using the updated velocity model (via the Fast Marching Method)
  4. Accounts for topography — essential in volcanic settings

"Nonlinear multiscale inversion significantly outperforms both linear and single-scale nonlinear inversion, particularly when station coverage is sparse." — Cabrera-Pérez et al. (2021)

Data Selection Criteria¶

Following Bensen et al. (2007, 2008):

  • Minimum interstation distance: $\Delta > 3\lambda$ (3 wavelengths) to avoid near-field effects
  • Maximum period cutoff: $\tau_{\max} = \Delta / 12$
  • SNR threshold: SNR > 10 (relaxed to 7 for sparse networks)
  • 3-sigma residual rejection from tomographic maps
  • Temporal repeatability: Standard deviation from seasonal subsets < 100 m/s

Resolution Estimation¶

Resolution is quantified by fitting a 2D Gaussian to the resolution matrix columns (Bensen et al., 2008, Eq. 5):

$$R(\vec{r}) = A \exp\left(-\frac{|\vec{r}|^2}{2\gamma^2}\right) \qquad \text{Resolution} = 2\gamma$$

In [116]:
# =============================================================================
# PRACTICE 10.1: Simplified 2D Surface Wave Tomography
# =============================================================================
# Demonstrate the tomographic inversion process using synthetic data
# Following Bensen et al. (2008) and Rawlinson et al. (2010)

def simple_2d_tomography(station_coords, travel_times, grid_size=20,
                          damping=1.0, smoothing=0.5):
    """
    Simplified 2D surface wave tomography using ray theory.
    
    Parameters
    ----------
    station_coords : array (n_stations, 2)
        Station [x, y] coordinates in km
    travel_times : dict
        {(i,j): traveltime_seconds} for each station pair
    grid_size : int
        Number of grid cells in each dimension
    damping : float
        Damping parameter (epsilon in Rawlinson et al., 2010)
    smoothing : float
        Smoothing parameter (eta)
    
    Returns
    -------
    velocity_map : 2D array
        Inverted velocity map (km/s)
    """
    from scipy.sparse import csr_matrix, eye
    from scipy.sparse.linalg import lsqr
    
    # Define grid
    x_min, x_max = station_coords[:,0].min() - 10, station_coords[:,0].max() + 10
    y_min, y_max = station_coords[:,1].min() - 10, station_coords[:,1].max() + 10
    
    dx = (x_max - x_min) / grid_size
    dy = (y_max - y_min) / grid_size
    n_cells = grid_size * grid_size
    
    # Build the forward operator G (ray-theoretic path lengths through cells)
    pairs = list(travel_times.keys())
    n_data = len(pairs)
    
    # Slowness model: m = 1/velocity at each cell
    # Data: travel time = sum(path_length * slowness) along each ray
    
    G = np.zeros((n_data, n_cells))
    d = np.zeros(n_data)
    
    for idx, (i, j) in enumerate(pairs):
        x1, y1 = station_coords[i]
        x2, y2 = station_coords[j]
        d[idx] = travel_times[(i, j)]
        
        # Simple ray tracing: straight ray between stations
        # Compute path length through each grid cell
        n_steps = 500
        ray_x = np.linspace(x1, x2, n_steps)
        ray_y = np.linspace(y1, y2, n_steps)
        
        total_dist = np.sqrt((x2-x1)**2 + (y2-y1)**2)
        dl = total_dist / n_steps
        
        for k in range(n_steps):
            # Find which cell this point is in
            ix = int((ray_x[k] - x_min) / dx)
            iy = int((ray_y[k] - y_min) / dy)
            
            ix = min(max(ix, 0), grid_size - 1)
            iy = min(max(iy, 0), grid_size - 1)
            
            cell_idx = iy * grid_size + ix
            G[idx, cell_idx] += dl
    
    # Reference slowness (starting model)
    s_ref = np.mean(d) / np.mean([np.sqrt((station_coords[i][0]-station_coords[j][0])**2 + 
                                           (station_coords[i][1]-station_coords[j][1])**2)
                                  for i, j in pairs])
    
    # Build regularization: Laplacian smoothing
    L = np.zeros((n_cells, n_cells))
    for iy in range(grid_size):
        for ix in range(grid_size):
            idx = iy * grid_size + ix
            neighbors = 0
            if ix > 0:
                L[idx, idx - 1] = -1; neighbors += 1
            if ix < grid_size - 1:
                L[idx, idx + 1] = -1; neighbors += 1
            if iy > 0:
                L[idx, idx - grid_size] = -1; neighbors += 1
            if iy < grid_size - 1:
                L[idx, idx + grid_size] = -1; neighbors += 1
            L[idx, idx] = neighbors
    
    # Solve: [G^T G + eps*I + eta*L^T L] dm = G^T (d - G*m0)
    m0 = np.full(n_cells, s_ref)
    residual = d - G @ m0
    
    GTG = G.T @ G
    LTL = L.T @ L
    I = np.eye(n_cells)
    
    A = GTG + damping * I + smoothing * LTL
    b = G.T @ residual
    
    dm = np.linalg.solve(A, b)
    m_final = m0 + dm
    
    # Convert slowness to velocity
    velocity = 1.0 / m_final.reshape(grid_size, grid_size)
    
    return velocity, (x_min, x_max, y_min, y_max)


# Create synthetic velocity model and stations
np.random.seed(2024)
n_stations = 25
station_x = np.random.uniform(0, 200, n_stations)
station_y = np.random.uniform(0, 200, n_stations)
station_coords = np.column_stack([station_x, station_y])

# True velocity model: background + anomalies
grid_true = 50
x_true = np.linspace(0, 200, grid_true)
y_true = np.linspace(0, 200, grid_true)
X, Y = np.meshgrid(x_true, y_true)

v_true = 3.0 * np.ones_like(X)  # Background 3 km/s

# Add slow anomaly (sedimentary basin)
v_true -= 0.5 * np.exp(-((X-60)**2 + (Y-70)**2) / (2*25**2))
# Add fast anomaly (plutonic intrusion)
v_true += 0.4 * np.exp(-((X-140)**2 + (Y-130)**2) / (2*30**2))
# Add linear gradient
v_true += 0.002 * Y

# Compute synthetic travel times
travel_times = {}
for i in range(n_stations):
    for j in range(i+1, n_stations):
        x1, y1 = station_coords[i]
        x2, y2 = station_coords[j]
        
        # Integrate along straight ray through true model
        n_steps = 200
        ray_x = np.linspace(x1, x2, n_steps)
        ray_y = np.linspace(y1, y2, n_steps)
        dist = np.sqrt((x2-x1)**2 + (y2-y1)**2)
        dl = dist / n_steps
        
        tt = 0
        for k in range(n_steps):
            ix = int(ray_x[k] / 200 * (grid_true-1))
            iy = int(ray_y[k] / 200 * (grid_true-1))
            ix = min(max(ix, 0), grid_true-1)
            iy = min(max(iy, 0), grid_true-1)
            tt += dl / v_true[iy, ix]
        
        # Add noise (measurement uncertainty)
        tt += np.random.normal(0, 0.1)
        travel_times[(i, j)] = tt

# Perform tomographic inversion
v_inv, extent = simple_2d_tomography(station_coords, travel_times, 
                                      grid_size=15, damping=0.5, smoothing=1.0)

# Plot results
fig, axes = plt.subplots(1, 3, figsize=(18, 5))

# True model
im1 = axes[0].pcolormesh(x_true, y_true, v_true, cmap='RdBu_r', 
                          vmin=2.3, vmax=3.6)
axes[0].scatter(station_x, station_y, c='k', marker='^', s=50, zorder=5)
axes[0].set_title('True Velocity Model')
axes[0].set_xlabel('X (km)')
axes[0].set_ylabel('Y (km)')
plt.colorbar(im1, ax=axes[0], label='Vs (km/s)')

# Inverted model
x_inv = np.linspace(extent[0], extent[1], 15)
y_inv = np.linspace(extent[2], extent[3], 15)
im2 = axes[1].pcolormesh(x_inv, y_inv, v_inv, cmap='RdBu_r', 
                          vmin=2.3, vmax=3.6)
axes[1].scatter(station_x, station_y, c='k', marker='^', s=50, zorder=5)

# Draw ray paths for a few pairs
for i, j in list(travel_times.keys())[:30]:
    axes[1].plot([station_coords[i,0], station_coords[j,0]], 
                 [station_coords[i,1], station_coords[j,1]], 
                 'k-', alpha=0.05, lw=0.5)

axes[1].set_title('Inverted Velocity (Gauss-Newton)')
axes[1].set_xlabel('X (km)')
plt.colorbar(im2, ax=axes[1], label='Vs (km/s)')

# Path density
path_density = np.zeros((15, 15))
dx_inv = (extent[1]-extent[0]) / 15
dy_inv = (extent[3]-extent[2]) / 15
for (i, j) in travel_times:
    x1, y1 = station_coords[i]
    x2, y2 = station_coords[j]
    for k in np.linspace(0, 1, 100):
        px = x1 + k*(x2-x1)
        py = y1 + k*(y2-y1)
        ix = int((px - extent[0]) / dx_inv)
        iy = int((py - extent[2]) / dy_inv)
        if 0 <= ix < 15 and 0 <= iy < 15:
            path_density[iy, ix] += 1

im3 = axes[2].pcolormesh(x_inv, y_inv, path_density, cmap='YlOrRd')
axes[2].scatter(station_x, station_y, c='k', marker='^', s=50, zorder=5)
axes[2].set_title('Path Density (Ray Coverage)')
axes[2].set_xlabel('X (km)')
plt.colorbar(im3, ax=axes[2], label='Number of rays')

plt.suptitle('Module 10: 2D Surface Wave Tomography\n'
             '(Bensen et al., 2008; Rawlinson et al., 2010)', 
             fontsize=13, fontweight='bold')
plt.tight_layout()
plt.show()

print(f"Number of station pairs: {len(travel_times)}")
print(f"Grid cells: {15}x{15} = {15*15}")
print(f"Data/model ratio: {len(travel_times)}/{15*15} = {len(travel_times)/(15*15):.1f}")
print()
print("Key points:")
print("• The slow anomaly (basin) and fast anomaly (intrusion) are recovered")
print("• Resolution depends on path density — gaps in station coverage create artifacts")
print("• Regularization (damping + smoothing) controls the trade-off between")
print("  fitting the data and producing a smooth, physically plausible model")
print("• Real ANT studies use this at each period to produce period-dependent maps,")
print("  then invert depth profiles at each grid point (Bensen et al., 2008)")
No description has been provided for this image 
Number of station pairs: 300
Grid cells: 15x15 = 225
Data/model ratio: 300/225 = 1.3

Key points:
• The slow anomaly (basin) and fast anomaly (intrusion) are recovered
• Resolution depends on path density — gaps in station coverage create artifacts
• Regularization (damping + smoothing) controls the trade-off between
  fitting the data and producing a smooth, physically plausible model
• Real ANT studies use this at each period to produce period-dependent maps,
  then invert depth profiles at each grid point (Bensen et al., 2008)

Module 11: Seismic Interferometry Principles¶

9.1 What is Seismic Interferometry?¶

Seismic interferometry is the broader principle that underlies ambient noise tomography. It refers to the extraction of deterministic information (impulse responses, Green's functions) from the interference of seismic waves (Campillo & Roux, 2015; Wapenaar, 2004).

Historical Development¶

Year Contribution Reference
1957 Spatial autocorrelation (SPAC) method Aki (1957)
1968 Daylight imaging conjecture Claerbout (1968)
2001 Green's function from diffuse field correlations Lobkis & Weaver (2001)
2003 Coda wave cross-correlation in seismology Campillo & Paul (2003)
2004 Ambient noise surface wave tomography Shapiro & Campillo (2004)
2004 Elastodynamic Green's function retrieval Wapenaar (2004)
2004 Green's function from coda wave correlations Snieder (2004)
2006 Passive image interferometry for monitoring Sens-Schoenfelder & Wegler (2006)

"The daylight imaging method was pioneered by Claerbout (1968), who proposed the conjecture stating that the cross-correlation of two daylight traces at surface locations A and B is equivalent to a reflection trace at B generated by a source at A." -- Campillo & Roux (2015)

9.2 Types of Seismic Interferometry¶

9.2.1 Cross-Correlation Type¶

$$G(A,B) \propto \int_{\partial V} G(A,S) \star G(B,S) \, dS$$

The Green's function between A and B is retrieved by cross-correlating recordings at A and B from sources on a boundary $\partial V$.

9.2.2 Deconvolution Type¶

$$G(A,B) \propto \frac{U_A(\omega)}{U_B(\omega)}$$

The Green's function is retrieved by deconvolving the recording at one station by the other. This is useful when source properties are different at the two receivers (Jiang & Denolle, 2020, Eq. 3).

9.2.3 Cross-Coherence Type¶

$$G(A,B) \propto \frac{U_A^*(\omega) U_B(\omega)}{|U_A(\omega)| |U_B(\omega)|}$$

A normalized version that is equivalent to spectral whitening + cross-correlation (Prieto et al., 2009; Jiang & Denolle, 2020).

9.3 Requirements for Green's Function Retrieval¶

For the cross-correlation to converge to the true Green's function, several conditions should ideally be met (Campillo & Roux, 2015):

  1. Diffuse wavefield: The noise sources should be distributed isotropically (or the scattering should create an effectively diffuse field)
  2. Sufficient recording duration: Long time series allow more noise sources to contribute
  3. Equipartition: Equal energy in all wave modes and propagation directions

Practical Limitations (Campillo & Roux, 2015, Section 1.12.5)¶

  • Partial focusing: Non-uniform source distribution leads to biased Green's functions
  • Amplitude information: Only reliable when source distribution is uniform
  • Travel time measurements: More robust than amplitude, even with imperfect source distribution
  • Scattering: Helps redistribute energy and improve Green's function reconstruction
In [117]:
# =============================================================================
# PRACTICE 11.1: Seismic Interferometry — Virtual Seismograms
# =============================================================================
# Demonstrate how cross-correlation creates a "virtual source" seismogram
# Compare the cross-correlation with a synthetic Green's function

def synthetic_greens_function(dist, t, c_min=2.5, c_max=4.5, f_min=0.02, f_max=0.1, Q=100):
    """
    Create a synthetic dispersive Green's function.
    Simulates a Rayleigh wave with dispersion and attenuation.
    """
    npts = len(t)
    dt = t[1] - t[0]
    freqs = np.fft.rfftfreq(npts, d=dt)
    
    spec = np.zeros(len(freqs), dtype=complex)
    
    for i, f in enumerate(freqs):
        if f_min < f < f_max:
            # Dispersive velocity: linear variation with frequency
            c = c_min + (c_max - c_min) * (f - f_min) / (f_max - f_min)
            # Travel time
            tt = dist / c
            # Phase
            phase = 2 * np.pi * f * tt
            # Amplitude with geometric spreading and attenuation
            amp = 1.0 / np.sqrt(max(dist, 1.0)) * np.exp(-np.pi * f * tt / Q)
            spec[i] = amp * np.exp(-1j * phase)
    
    return np.fft.irfft(spec, n=npts)


# Generate synthetic Green's function
t_syn = np.arange(-1000, 1000, 0.05)  # 20 Hz sampling
gf = synthetic_greens_function(dist_km, t_syn[t_syn >= 0])

# Create full symmetric Green's function
gf_full = np.zeros(len(t_syn))
center = len(t_syn) // 2
gf_full[center:center+len(gf)] = gf
gf_full[:center+1] += gf_full[:center+1][::-1]  # Add acausal part

# Plot comparison
fig, axes = plt.subplots(2, 1, figsize=(16, 8))

# Synthetic Green's function
axes[0].plot(t_syn, gf_full / np.max(np.abs(gf_full)), 'r-', linewidth=1.0)
axes[0].set_title("Synthetic Green's Function (dispersive Rayleigh wave)", fontweight='bold')
axes[0].set_ylabel('Normalized Amplitude')
axes[0].set_xlim(-800, 800)
axes[0].grid(True, alpha=0.3)

# Empirical Green's function from noise cross-correlation
cc_n = cc_stack / np.max(np.abs(cc_stack))
axes[1].plot(lags, cc_n, 'k-', linewidth=1.0)
axes[1].set_title("Empirical Green's Function (from noise cross-correlation)", fontweight='bold')
axes[1].set_ylabel('Normalized CC')
axes[1].set_xlabel('Lag Time (s)')
axes[1].set_xlim(-800, 800)
axes[1].grid(True, alpha=0.3)

plt.suptitle('Seismic Interferometry: Comparing Virtual and Theoretical Seismograms\n'
             f'Station pair: ANMO-CCM, Distance: {dist_km:.0f} km',
             y=1.02, fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()

print("The noise cross-correlation (bottom) resembles the theoretical")
print("Green's function (top), showing dispersive surface wave arrivals.")
print("\nThis is the fundamental result of seismic interferometry:")
print("noise cross-correlation ≈ Green's function between the stations.")
No description has been provided for this image 
The noise cross-correlation (bottom) resembles the theoretical
Green's function (top), showing dispersive surface wave arrivals.

This is the fundamental result of seismic interferometry:
noise cross-correlation ≈ Green's function between the stations.

Module 12: Seismic Velocity Changes (dv/v) Monitoring¶

10.1 The Principle of Noise-Based Monitoring¶

One of the most powerful applications of ambient noise is the continuous monitoring of seismic velocity changes in the subsurface. This approach was pioneered by Sens-Schoenfelder & Wegler (2006) and has since been applied extensively to volcanoes, fault zones, and other geological settings.

"Sens-Schoenfelder and Wegler (2006) were the first to propose the consideration of daily cross-correlations of noise as repetitive in time Green's functions, which can be used to assess seismic velocity temporal changes." -- Brenguier et al. (2016)

The Concept (Brenguier et al., 2016, Fig. 1)¶

  1. Compute daily noise cross-correlations between a pair of stations
  2. Compare each daily cross-correlation with a reference (long-term average)
  3. Measure the time shift $\delta t$ in the coda (late arrivals) of the cross-correlation
  4. The relative velocity change is: $\frac{\delta v}{v} = -\frac{\delta t}{t}$

Why the Coda?¶

The coda (late-arriving scattered waves) is particularly sensitive to velocity changes because:

  • Coda waves have traveled long distances through the medium (via multiple scattering)
  • Small velocity changes accumulate over the long path lengths
  • The coda is remarkably stable between repeated measurements
  • Precision of $\delta v/v$ can reach $10^{-5}$ (0.001%) — detecting velocity changes of 0.02 m/s on a 2000 m/s background (Brenguier et al., 2016)

10.2 Measurement Methods¶

Stretching Technique¶

The effect of a global velocity change $\epsilon$ is to deform the time axis by $t' = t(1 - \epsilon)$. The stretching technique (Lobkis & Weaver, 2003; Sens-Schoenfelder & Wegler, 2006):

  1. Stretch the current cross-correlation by various factors $\epsilon$
  2. Find the $\epsilon_0$ that maximizes the cross-correlation coefficient with the reference
  3. $\epsilon_0 = \delta v / v$

Moving-Window Cross-Spectral Analysis (MWCSA)¶

Proposed by Poupinet et al. (1984) and Snieder et al. (2002):

  1. Measure $\delta t$ in moving windows along the coda
  2. Fit a linear regression: $\delta t = -(\delta v/v) \cdot t$
  3. The slope gives $\delta v/v$
In [118]:
# =============================================================================
# PRACTICE 12.1: Implement the Stretching Technique for dv/v
# =============================================================================
# Following Sens-Schoenfelder & Wegler (2006) and Brenguier et al. (2016)

def stretching_dvv(ref_cc, cur_cc, lags, t_min, t_max, dv_range=0.01, n_dv=201):
    """
    Measure relative velocity change (dv/v) using the stretching technique.
    
    Reference: Sens-Schoenfelder & Wegler (2006), Brenguier et al. (2016)
    Also implemented in NoisePy's measure_dvv module.
    
    Parameters
    ----------
    ref_cc : numpy array
        Reference cross-correlation (long-term average)
    cur_cc : numpy array
        Current cross-correlation to compare
    lags : numpy array
        Lag times
    t_min, t_max : float
        Coda window boundaries (seconds)
    dv_range : float
        Maximum dv/v to search (fraction)
    n_dv : int
        Number of dv/v values to test
    
    Returns
    -------
    dvv : float
        Best-fit relative velocity change
    cc_coeff : float
        Maximum correlation coefficient
    dvv_array : numpy array
        Tested dv/v values
    cc_array : numpy array
        Correlation coefficients for each dv/v
    """
    dvv_array = np.linspace(-dv_range, dv_range, n_dv)
    cc_array = np.zeros(n_dv)
    
    # Select coda window (use both positive and negative lags)
    mask = ((np.abs(lags) >= t_min) & (np.abs(lags) <= t_max))
    ref_coda = ref_cc[mask]
    t_coda = lags[mask]
    
    for i, dv in enumerate(dvv_array):
        # Stretch: t' = t * (1 - dv/v)
        t_stretched = t_coda * (1 - dv)
        
        # Interpolate current CC at stretched times
        cur_stretched = np.interp(t_stretched, lags, cur_cc)
        
        # Correlation coefficient
        if np.std(ref_coda) > 0 and np.std(cur_stretched) > 0:
            cc_array[i] = np.corrcoef(ref_coda, cur_stretched)[0, 1]
    
    # Best dv/v
    best_idx = np.argmax(cc_array)
    dvv = dvv_array[best_idx]
    cc_coeff = cc_array[best_idx]
    
    return dvv, cc_coeff, dvv_array, cc_array


# Demonstrate with synthetic velocity change
# Create a reference CC and a "perturbed" CC with known velocity change
ref_cc = cc_stack.copy()

# Simulate a velocity change of -0.3%
true_dvv = -0.003  # -0.3%
# Stretch the reference to create the "current" measurement
lags_stretched = lags * (1 - true_dvv)
cur_cc = np.interp(lags, lags_stretched, ref_cc)
# Add some noise
cur_cc += 0.05 * np.max(np.abs(cur_cc)) * np.random.randn(len(cur_cc))

# Measure dv/v
dvv, cc_coeff, dvv_arr, cc_arr = stretching_dvv(
    ref_cc, cur_cc, lags, 
    t_min=200, t_max=600,  # Coda window
    dv_range=0.01, n_dv=401
)

# Plot results
fig, axes = plt.subplots(2, 2, figsize=(16, 10))

# Reference vs Current CC
axes[0,0].plot(lags, ref_cc / np.max(np.abs(ref_cc)), 'b-', linewidth=0.8, label='Reference')
axes[0,0].plot(lags, cur_cc / np.max(np.abs(cur_cc)), 'r-', linewidth=0.8, 
               alpha=0.7, label='Current (dv/v = -0.3%)')
axes[0,0].set_xlim(-700, 700)
axes[0,0].axvspan(200, 600, alpha=0.1, color='green', label='Coda window')
axes[0,0].axvspan(-600, -200, alpha=0.1, color='green')
axes[0,0].set_title('Reference vs Current Cross-Correlation', fontweight='bold')
axes[0,0].set_xlabel('Lag (s)')
axes[0,0].legend(fontsize=9)
axes[0,0].grid(True, alpha=0.3)

# Zoom on coda
coda_mask = (lags > 300) & (lags < 500)
axes[0,1].plot(lags[coda_mask], ref_cc[coda_mask] / np.max(np.abs(ref_cc)), 
               'b-', linewidth=1.5, label='Reference')
axes[0,1].plot(lags[coda_mask], cur_cc[coda_mask] / np.max(np.abs(cur_cc)), 
               'r-', linewidth=1.5, alpha=0.7, label='Current')
axes[0,1].set_title('Zoom on Coda: Time Shift Visible!', fontweight='bold')
axes[0,1].set_xlabel('Lag (s)')
axes[0,1].legend()
axes[0,1].grid(True, alpha=0.3)

# Stretching curve
axes[1,0].plot(dvv_arr * 100, cc_arr, 'k-', linewidth=1.5)
axes[1,0].axvline(dvv * 100, color='red', linestyle='--', linewidth=2,
                   label=f'Measured: {dvv*100:.3f}%')
axes[1,0].axvline(true_dvv * 100, color='blue', linestyle=':', linewidth=2,
                   label=f'True: {true_dvv*100:.3f}%')
axes[1,0].set_title('Stretching: CC Coefficient vs dv/v', fontweight='bold')
axes[1,0].set_xlabel('dv/v (%)')
axes[1,0].set_ylabel('Correlation Coefficient')
axes[1,0].legend()
axes[1,0].grid(True, alpha=0.3)

# Summary text
axes[1,1].axis('off')
summary = (
    f"Stretching Technique Results\n"
    f"{'='*40}\n\n"
    f"True dv/v:     {true_dvv*100:.3f}%\n"
    f"Measured dv/v: {dvv*100:.3f}%\n"
    f"CC coefficient: {cc_coeff:.4f}\n"
    f"Error: {abs(dvv - true_dvv)*100:.4f}%\n\n"
    f"Coda window: {200}-{600} s\n\n"
    f"References:\n"
    f"Sens-Schoenfelder & Wegler (2006)\n"
    f"Brenguier et al. (2016)\n"
    f"NoisePy: measure_dvv module"
)
axes[1,1].text(0.1, 0.5, summary, transform=axes[1,1].transAxes, fontsize=12,
               verticalalignment='center', fontfamily='monospace',
               bbox=dict(boxstyle='round', facecolor='lightyellow', alpha=0.8))

plt.suptitle('Seismic Velocity Change Measurement (dv/v)\n'
             'Stretching Technique (Sens-Schoenfelder & Wegler, 2006)',
             y=1.02, fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()
No description has been provided for this image
In [119]:
# =============================================================================
# PRACTICE 12.2: Simulate Time-Lapse dv/v Monitoring
# =============================================================================
# Simulate a time series of dv/v measurements as done in volcanic monitoring
# Reference: Brenguier et al. (2016), Fig. 4 & 5

# Simulate 365 days of measurements with seasonal variation + event
n_days = 365
days = np.arange(n_days)

# True dv/v signal:
# - Seasonal variation (due to groundwater, thermal effects)
# - A volcanic precursor drop around day 200
true_dvv_ts = (
    0.1 * np.sin(2 * np.pi * days / 365)  # Seasonal (Hillers et al., 2015)
    - 0.3 * np.exp(-((days - 200) / 10) ** 2)  # Volcanic precursor
    + 0.2 * (1 - np.exp(-np.maximum(days - 210, 0) / 50))  # Post-eruption recovery
)

# Simulate noisy measurements
np.random.seed(42)
measured_dvv = true_dvv_ts + 0.05 * np.random.randn(n_days)
cc_coeffs = 0.95 + 0.04 * np.random.randn(n_days)
cc_coeffs = np.clip(cc_coeffs, 0, 1)

# Plot like Brenguier et al. (2016), Fig. 4
fig, axes = plt.subplots(3, 1, figsize=(16, 10), sharex=True)

# dv/v time series
axes[0].plot(days, measured_dvv, 'b.', markersize=3, alpha=0.5, label='Daily measurement')
axes[0].plot(days, true_dvv_ts, 'r-', linewidth=2, label='True signal')
# 5-day moving average
smooth = np.convolve(measured_dvv, np.ones(5)/5, mode='same')
axes[0].plot(days, smooth, 'k-', linewidth=1.5, label='5-day average')
axes[0].axvspan(200, 210, alpha=0.3, color='red', label='Eruption period')
axes[0].set_ylabel('dv/v (%)', fontsize=13)
axes[0].set_title('Seismic Velocity Changes (dv/v) — Volcanic Monitoring', fontweight='bold')
axes[0].legend(loc='upper right', fontsize=9)
axes[0].grid(True, alpha=0.3)
axes[0].invert_yaxis()  # Convention: velocity decrease is negative, plotted upward

# CC coefficient (quality indicator)
axes[1].plot(days, cc_coeffs, 'g-', linewidth=0.5, alpha=0.7)
axes[1].axhline(0.9, color='red', linestyle='--', label='Quality threshold')
axes[1].set_ylabel('CC Coefficient', fontsize=13)
axes[1].set_title('Waveform Similarity (Quality Control)', fontweight='bold')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

# Simulated seismic energy
seismic_energy = 0.1 * np.random.exponential(1, n_days)
seismic_energy[200:210] = np.random.exponential(5, 10)  # Eruption
axes[2].bar(days, seismic_energy, width=1, color='gray', alpha=0.7)
axes[2].set_ylabel('Seismic Energy', fontsize=13)
axes[2].set_xlabel('Day of Year', fontsize=13)
axes[2].set_title('Seismic Activity', fontweight='bold')
axes[2].grid(True, alpha=0.3)

plt.suptitle('Time-Lapse Monitoring with Ambient Noise\n'
             '(Inspired by Brenguier et al., 2016, Fig. 4-5: Piton de la Fournaise)',
             y=1.02, fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()

print("Key observations (cf. Brenguier et al., 2016):")
print("1. Seasonal velocity variations (~0.1%) due to groundwater/thermal effects")
print("2. Pre-eruptive velocity DECREASE: precursor signal days before eruption")
print("3. Post-eruption velocity RECOVERY: slow return to normal")
print(f"4. Precision of dv/v measurements can reach 10^-5 (Brenguier et al., 2016)")
No description has been provided for this image 
Key observations (cf. Brenguier et al., 2016):
1. Seasonal velocity variations (~0.1%) due to groundwater/thermal effects
2. Pre-eruptive velocity DECREASE: precursor signal days before eruption
3. Post-eruption velocity RECOVERY: slow return to normal
4. Precision of dv/v measurements can reach 10^-5 (Brenguier et al., 2016)

Module 13: HVSR / MHVSR Method¶

13.1 The Horizontal-to-Vertical Spectral Ratio¶

The Microtremor Horizontal-to-Vertical Spectral Ratio (MHVSR) method is a single-station technique that uses ambient noise to characterize site effects (Molnar et al., 2022).

History¶

  • First proposed by Nogoshi & Igarashi (1971) in Japan
  • Popularized internationally by Nakamura (1989), who interpreted it as the S-wave transfer function
  • Standardized by the SESAME project (2004) in Europe
  • Comprehensive review by Molnar et al. (2022) covering theory, methods, and applications

Definition¶

$$\text{HVSR}(f) = \frac{\sqrt{S_N^2(f) + S_E^2(f)}}{S_Z(f)} \qquad [\text{Molnar et al., 2022}]$$

where $S_N$, $S_E$, $S_Z$ are the Fourier Amplitude Spectra of the North, East, and Vertical components.

Alternatively, the geometric mean of horizontal spectra is used:

$$\text{HVSR}(f) = \frac{\sqrt{S_N(f) \cdot S_E(f)}}{S_Z(f)} \qquad [\text{Cox et al., 2020}]$$

Site Resonance Frequency¶

The peak frequency $f_0$ of the HVSR corresponds to the fundamental resonance frequency of the site:

$$f_0 = \frac{V_s}{4H}$$

where $V_s$ is the shear-wave velocity of the soft layer and $H$ is its thickness.

Applications¶

  • Site characterization: Identifying the fundamental resonance frequency for seismic hazard assessment
  • Vs30 estimation: Classifying sites per building codes (NBCC, Eurocode 8)
  • Bedrock depth mapping: HVSR survey profiles to map depth to bedrock (Hunter & Crow, 2015)
  • Mineral exploration: HVSR surveys for mapping sediment thickness over mineral deposits

13.2 Statistical Treatment: The Lognormal Distribution (Cox et al., 2020)¶

Cox et al. (2020) demonstrated that the resonance frequency $f_0$ follows a lognormal distribution, not a normal distribution. This has important implications for how HVSR results should be reported.

Key Findings¶

  • The lognormal median (LM) should be used instead of the arithmetic mean:

$$\text{LM}_{f_0} = \exp(\mu_{\ln f_0})$$

  • The uncertainty should be expressed as the lognormal standard deviation:

$$\sigma_{\ln f_0} = \text{std}(\ln f_0)$$

Frequency-Domain Window-Rejection Algorithm¶

Cox et al. (2020) proposed an iterative rejection algorithm:

  1. Compute $f_0$ for each time window
  2. Accept window $i$ only if:

$$f_{0,i} \in \left[\exp(\mu_{\ln f_0} - n \sigma_{\ln f_0}), \; \exp(\mu_{\ln f_0} + n \sigma_{\ln f_0})\right]$$

where $n = 2$ (analogous to 2-sigma rule but in log space).

  1. Recompute statistics and repeat until convergence:
    • $|\Delta N_{\text{accepted}}| / N < 1\%$
    • $|\Delta \sigma_{\ln f_0}| < 0.01$

Software: The hvsrpy Python package (Vantassel, 2020) implements this algorithm.

SESAME Guidelines¶

The SESAME (2004) project established reliability criteria:

  • Minimum recording duration: $T > 200 / f_0$ seconds
  • Minimum number of windows: $n_w > 200 \cdot f_0 \cdot l_w$ where $l_w$ is window length
  • Window length: At least 10 cycles of $f_0$, recommended $l_w > 10/f_0$
  • Peak amplitude: $A_{\text{HVSR}}(f_0) > 2$ for a clear resonance peak
In [120]:
# =============================================================================
# PRACTICE 13.1: Compute MHVSR from 3-Component Data
# =============================================================================
# Implement the MHVSR method following Molnar et al. (2022)

def compute_hvsr(st_3c, window_length=50, overlap=0.5, 
                  freqmin=0.1, freqmax=25.0, smooth_width=40):
    """
    Compute the Horizontal-to-Vertical Spectral Ratio (HVSR).
    
    Following Molnar et al. (2022) and SESAME (2004) guidelines.
    
    Parameters
    ----------
    st_3c : obspy.Stream
        3-component stream (Z, N, E channels)
    window_length : float
        Window length in seconds
    overlap : float
        Fractional overlap between windows
    freqmin, freqmax : float
        Frequency range for analysis
    smooth_width : int
        Konno-Ohmachi smoothing bandwidth parameter
    
    Returns
    -------
    freqs : numpy array
        Frequencies
    hvsr_mean : numpy array
        Mean HVSR curve
    hvsr_std : numpy array
        Standard deviation of HVSR
    hvsr_all : list
        Individual window HVSR curves
    """
    # Sort components
    tr_z = st_3c.select(component='Z')[0].copy()
    tr_n = st_3c.select(component='N')[0].copy()
    tr_e = st_3c.select(component='E')[0].copy()
    
    for tr in [tr_z, tr_n, tr_e]:
        tr.detrend('demean')
        tr.detrend('linear')
        tr.taper(0.05)
    
    fs = tr_z.stats.sampling_rate
    npts_win = int(window_length * fs)
    npts_step = int(npts_win * (1 - overlap))
    
    min_len = min(len(tr_z.data), len(tr_n.data), len(tr_e.data))
    n_windows = (min_len - npts_win) // npts_step + 1
    
    hvsr_all = []
    
    for i in range(n_windows):
        start = i * npts_step
        end = start + npts_win
        if end > min_len:
            break
        
        # Apply taper to each window
        taper = cosine_taper(npts_win, p=0.1)
        
        # Compute FFT
        fft_z = np.abs(np.fft.rfft(tr_z.data[start:end] * taper))
        fft_n = np.abs(np.fft.rfft(tr_n.data[start:end] * taper))
        fft_e = np.abs(np.fft.rfft(tr_e.data[start:end] * taper))
        
        freqs = np.fft.rfftfreq(npts_win, d=1.0/fs)
        
        # Smooth spectra (simple moving average as proxy for Konno-Ohmachi)
        from scipy.ndimage import uniform_filter1d
        fft_z_s = uniform_filter1d(fft_z, smooth_width)
        fft_n_s = uniform_filter1d(fft_n, smooth_width)
        fft_e_s = uniform_filter1d(fft_e, smooth_width)
        
        # HVSR = sqrt((H_N^2 + H_E^2) / V^2)
        h_combined = np.sqrt(fft_n_s**2 + fft_e_s**2)
        fft_z_s[fft_z_s < 1e-20] = 1e-20
        hvsr = h_combined / fft_z_s
        
        # Frequency mask
        freq_mask = (freqs >= freqmin) & (freqs <= freqmax)
        hvsr_all.append(hvsr)
    
    hvsr_all = np.array(hvsr_all)
    hvsr_mean = np.mean(hvsr_all, axis=0)
    hvsr_std = np.std(hvsr_all, axis=0)
    
    return freqs, hvsr_mean, hvsr_std, hvsr_all


# Download 3-component data
print("Downloading 3-component data for HVSR analysis...")
client = Client("IRIS")
t1 = UTCDateTime("2022-06-15T02:00:00")  # Nighttime for less anthropogenic noise
t2 = t1 + 3600  # 1 hour

st_3c = Stream()
for chan in ["BHZ", "BH1", "BH2"]:
    st_3c += client.get_waveforms("IU", "ANMO", "00", chan, t1, t2)

print(st_3c)

rename = {'1': 'N', '2': 'E'}
for tr in st_3c:
    code = tr.stats.channel[-1]
    if code in rename:
        tr.stats.channel = tr.stats.channel[:-1] + rename[code]

print(st_3c)

# Compute HVSR
freqs_h, hvsr_mean, hvsr_std, hvsr_all = compute_hvsr(
    st_3c, window_length=60, overlap=0.5, smooth_width=30
)

# Plot HVSR curve
fig, axes = plt.subplots(1, 2, figsize=(16, 6))

# Individual curves and mean
freq_mask = (freqs_h >= 0.1) & (freqs_h <= 15)
for i in range(min(30, len(hvsr_all))):
    axes[0].semilogx(freqs_h[freq_mask], hvsr_all[i][freq_mask], 
                     'gray', linewidth=0.3, alpha=0.3)

axes[0].semilogx(freqs_h[freq_mask], hvsr_mean[freq_mask], 'k-', linewidth=2, 
                 label='Mean HVSR')
axes[0].fill_between(freqs_h[freq_mask], 
                     (hvsr_mean - hvsr_std)[freq_mask],
                     (hvsr_mean + hvsr_std)[freq_mask],
                     alpha=0.3, color='blue', label='$\\pm 1\\sigma$')
axes[0].set_xlabel('Frequency (Hz)', fontsize=13)
axes[0].set_ylabel('H/V Ratio', fontsize=13)
axes[0].set_title('MHVSR Curve (Molnar et al., 2022)', fontweight='bold')
axes[0].axhline(1, color='red', linestyle='--', alpha=0.5)
axes[0].legend()
axes[0].grid(True, alpha=0.3, which='both')
axes[0].set_ylim(0, 8)

# Individual component spectra
tr_z = st_3c.select(component='Z')[0].copy()
tr_n = st_3c.select(component='N')[0].copy()
tr_e = st_3c.select(component='E')[0].copy()
for tr in [tr_z, tr_n, tr_e]:
    tr.detrend('demean')
    tr.detrend('linear')
    tr.taper(0.05)

for tr, label, color in zip([tr_z, tr_n, tr_e], 
                             ['Vertical', 'North', 'East'],
                             ['blue', 'red', 'green']):
    f, p = signal.welch(tr.data, fs=tr.stats.sampling_rate, nperseg=int(60*fs))
    f_mask = (f >= 0.1) & (f <= 15)
    axes[1].loglog(f[f_mask], p[f_mask], color=color, linewidth=1.5, label=label)

axes[1].set_xlabel('Frequency (Hz)', fontsize=13)
axes[1].set_ylabel('PSD', fontsize=13)
axes[1].set_title('Individual Component Spectra', fontweight='bold')
axes[1].legend()
axes[1].grid(True, alpha=0.3, which='both')

plt.suptitle('MHVSR Analysis (Molnar et al., 2022; SESAME, 2004)\n'
             f'Station: IU.ANMO',
             y=1.02, fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()

# Find peak frequency
peak_idx = np.argmax(hvsr_mean[freq_mask])
f0 = freqs_h[freq_mask][peak_idx]
print(f"\nFundamental resonance frequency (f0): {f0:.2f} Hz")
print(f"Corresponding period: {1/f0:.1f} s")
print(f"\nIf Vs = 300 m/s, estimated sediment thickness h = Vs/(4*f0) = {300/(4*f0):.0f} m")

Downloading 3-component data for HVSR analysis...
3 Trace(s) in Stream:
IU.ANMO.00.BHZ | 2022-06-15T02:00:00.019538Z - 2022-06-15T02:59:59.994538Z | 40.0 Hz, 144000 samples
IU.ANMO.00.BH1 | 2022-06-15T02:00:00.019539Z - 2022-06-15T02:59:59.994539Z | 40.0 Hz, 144000 samples
IU.ANMO.00.BH2 | 2022-06-15T02:00:00.019539Z - 2022-06-15T02:59:59.994539Z | 40.0 Hz, 144000 samples
3 Trace(s) in Stream:
IU.ANMO.00.BHZ | 2022-06-15T02:00:00.019538Z - 2022-06-15T02:59:59.994538Z | 40.0 Hz, 144000 samples
IU.ANMO.00.BHN | 2022-06-15T02:00:00.019539Z - 2022-06-15T02:59:59.994539Z | 40.0 Hz, 144000 samples
IU.ANMO.00.BHE | 2022-06-15T02:00:00.019539Z - 2022-06-15T02:59:59.994539Z | 40.0 Hz, 144000 samples
No description has been provided for this image 
Fundamental resonance frequency (f0): 6.12 Hz
Corresponding period: 0.2 s

If Vs = 300 m/s, estimated sediment thickness h = Vs/(4*f0) = 12 m
In [121]:
# =============================================================================
# PRACTICE 14.1: NoisePy-Style Complete Workflow
# =============================================================================
# Implement a simplified version of the NoisePy workflow
# This demonstrates the same processing pipeline used in NoisePy
# but in a single notebook for educational purposes.

class AmbientNoiseProcessor:
    """
    A simplified ambient noise processor inspired by NoisePy.
    
    Implements the workflow from Jiang & Denolle (2020):
    S0 (data preparation) → S1 (cross-correlation) → S2 (stacking)
    
    For production use, see the actual NoisePy package:
    https://github.com/noisepy/NoisePy
    """
    
    def __init__(self, cc_len=1800, step=450, freqmin=0.02, freqmax=0.5,
                 time_norm='one-bit', freq_norm='whiten', maxlag=500,
                 max_over_std=10):
        self.cc_len = cc_len
        self.step = step
        self.freqmin = freqmin
        self.freqmax = freqmax
        self.time_norm = time_norm
        self.freq_norm = freq_norm
        self.maxlag = maxlag
        self.max_over_std = max_over_std
    
    def preprocess(self, tr):
        """S0: Data preparation (Jiang & Denolle, 2020)."""
        tr_proc = tr.copy()
        fs = tr_proc.stats.sampling_rate
        
        # Remove trend and mean
        tr_proc.detrend('demean')
        tr_proc.detrend('linear')
        tr_proc.taper(0.05, type='cosine')
        
        # Bandpass filter
        tr_proc.filter('bandpass', freqmin=self.freqmin, freqmax=self.freqmax,
                       corners=4, zerophase=True)
        
        # Temporal normalization
        if self.time_norm == 'one-bit':
            tr_proc.data = np.sign(tr_proc.data)
        elif self.time_norm == 'rma':
            N = int(50 * fs)  # 50-second window
            w = np.convolve(np.abs(tr_proc.data), np.ones(2*N+1)/(2*N+1), mode='same')
            w[w < 1e-10] = 1e-10
            tr_proc.data /= w
        
        # Spectral whitening
        if self.freq_norm == 'whiten':
            spec = np.fft.rfft(tr_proc.data)
            amp = np.abs(spec)
            amp_s = np.convolve(amp, np.ones(20)/20, mode='same')
            amp_s[amp_s < 1e-10] = 1e-10
            tr_proc.data = np.fft.irfft(spec / amp_s, n=len(tr_proc.data))
        elif self.freq_norm == 'phase_only':
            spec = np.fft.rfft(tr_proc.data)
            amp = np.abs(spec)
            amp[amp < 1e-10] = 1e-10
            tr_proc.data = np.fft.irfft(spec / amp, n=len(tr_proc.data))
        
        return tr_proc
    
    def cross_correlate(self, tr1, tr2):
        """S1: Cross-correlation (Jiang & Denolle, 2020, Eq. 1)."""
        t1 = self.preprocess(tr1)
        t2 = self.preprocess(tr2)
        
        fs = t1.stats.sampling_rate
        npts_cc = int(self.cc_len * fs)
        npts_step = int(self.step * fs)
        npts_lag = int(self.maxlag * fs)
        
        min_len = min(len(t1.data), len(t2.data))
        nfft = int(2 ** np.ceil(np.log2(2 * npts_cc)))
        
        cc_all = []
        
        for start in range(0, min_len - npts_cc, npts_step):
            seg1 = t1.data[start:start + npts_cc]
            seg2 = t2.data[start:start + npts_cc]
            
            F1 = np.fft.rfft(seg1, n=nfft)
            F2 = np.fft.rfft(seg2, n=nfft)
            CC = np.conj(F1) * F2
            
            cc = np.fft.irfft(CC, n=nfft)
            cc = np.fft.fftshift(cc)
            
            center = nfft // 2
            cc_trim = cc[center - npts_lag:center + npts_lag + 1]
            
            # Quality control (NoisePy max_over_std)
            if np.max(np.abs(cc_trim)) < self.max_over_std * np.median(np.abs(cc_trim)):
                cc_all.append(cc_trim)
        
        lags = np.arange(-npts_lag, npts_lag + 1) / fs
        return lags, np.array(cc_all)
    
    def stack(self, cc_all, method='linear'):
        """S2: Stacking (Jiang & Denolle, 2020)."""
        if method == 'linear':
            return np.mean(cc_all, axis=0)
        elif method == 'pws':
            return phase_weighted_stack(cc_all)
        elif method == 'robust':
            return robust_stack(cc_all)
    
    def process_pair(self, tr1, tr2, stack_method='linear'):
        """Complete processing for one station pair."""
        lags, cc_all = self.cross_correlate(tr1, tr2)
        cc_stacked = self.stack(cc_all, method=stack_method)
        
        return {
            'lags': lags,
            'cc_all': cc_all,
            'cc_stack': cc_stacked,
            'n_windows': len(cc_all),
            'station1': tr1.stats.station,
            'station2': tr2.stats.station,
        }


# Run the complete workflow
print("Running NoisePy-style ambient noise processing...")
print("="*50)

processor = AmbientNoiseProcessor(
    cc_len=1800,
    step=450,
    freqmin=0.02,
    freqmax=0.2,
    time_norm='one-bit',
    freq_norm='whiten',
    maxlag=600
)

result = processor.process_pair(st1[0], st2[0], stack_method='linear')

print(f"\nProcessing complete!")
print(f"Station pair: {result['station1']}-{result['station2']}")
print(f"Windows processed: {result['n_windows']}")
print(f"Lag range: {result['lags'][0]:.0f} to {result['lags'][-1]:.0f} s")

# Plot final result
fig, ax = plt.subplots(figsize=(16, 5))
cc_n = result['cc_stack'] / np.max(np.abs(result['cc_stack']))
ax.plot(result['lags'], cc_n, 'k-', linewidth=1.0)
ax.fill_between(result['lags'], cc_n, 0, where=cc_n > 0, alpha=0.3, color='red')
ax.fill_between(result['lags'], cc_n, 0, where=cc_n < 0, alpha=0.3, color='blue')
ax.axvline(dist_km/3.5, color='green', linestyle='--', linewidth=1.5,
           label=f'Expected arrival (~3.5 km/s): {dist_km/3.5:.0f} s')
ax.axvline(-dist_km/3.5, color='green', linestyle='--', linewidth=1.5)
ax.set_title(f'NoisePy-Style Cross-Correlation: {result["station1"]}-{result["station2"]}\n'
             f'Distance: {dist_km:.0f} km | {result["n_windows"]} windows stacked',
             fontweight='bold')
ax.set_xlabel('Lag Time (s)')
ax.set_ylabel('Normalized CC')
ax.legend(loc='upper right')
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

Running NoisePy-style ambient noise processing...
==================================================

Processing complete!
Station pair: ANMO-HRV
Windows processed: 92
Lag range: -600 to 600 s
No description has been provided for this image

Module 15: Applications and Case Studies¶

15.1 Ambient Noise Tomography (ANT)¶

The primary application of ambient noise cross-correlations is surface wave tomography. The workflow (Shapiro & Campillo, 2004; Bensen et al., 2007):

  1. Compute cross-correlations for all station pairs
  2. Measure dispersion curves (group/phase velocity vs. period)
  3. Invert for velocity maps at each period
  4. Invert velocity maps for 3D velocity structure

15.1.1 Continental-Scale: United States (Bensen et al., 2008)¶

Bensen et al. (2008) produced the first broadband ambient noise tomography of the contiguous US:

  • 203 broadband stations, 24 months of data (2003-2005)
  • Both Rayleigh (Z-Z, R-R) and Love (T-T) waves extracted
  • Group and phase velocity maps from 8 to 70 s period
  • Resolution better than 100 km across most of the US

Key geological features revealed:

  • East-west velocity dichotomy between the stable craton (fast) and tectonically active west (slow)
  • Sedimentary basins (Gulf of Mexico, Anadarko, Williston) as slow anomalies at short periods
  • Yellowstone plume as a strong low-velocity anomaly at 40 s period
  • Moho depth variations mapped through the transition from crustal to mantle sensitivity

15.1.2 Dense Array: Groningen Gas Field (Chmiel et al., 2019)¶

Chmiel et al. (2019) deployed ~400 stations at ~350 m spacing:

  • Resolved 3D S-wave velocity to 1 km depth above the Groningen gas field
  • Joint inversion of Love and Rayleigh waves (fundamental mode + 1st higher mode)
  • Identified geological structures including the Brussel sands and a palaeochannel, validated against boreholes

15.1.3 Volcanic Island: La Palma (Cabrera-Pérez et al., 2023)¶

Cabrera-Pérez et al. (2023) applied ANT to La Palma (Canary Islands):

  • 44 broadband stations deployed 2018-2020
  • 415 valid dispersion curves extracted, periods 0.35-3.2 s
  • Preprocessing: one-bit temporal normalization, spectral whitening, ZZ cross-correlations
  • MANgOSTA nonlinear multiscale inversion accounting for island topography
  • Transdimensional Bayesian depth inversion (no subjective regularization)

Key findings:

  • H1: High-velocity anomaly (+20-40%) — ancient basal complex plutonic intrusion (4.0-3.0 Ma)
  • L1: Low-velocity anomaly (-20 to -40%) on western Cumbre Vieja — active hydrothermal system; dv/v decrease observed before the 2021 Tajogaite eruption — ascending hydrothermal fluids
  • L2: Eastern flank — fossil hydrothermal alteration zone
  • L3: Near San Antonio/Teneguía volcanoes — highly fractured rocks with geothermal activity

The pre-eruptive velocity decrease in L1 provides evidence that ANT can detect fluid migration precursory to eruption. — Cabrera-Pérez et al. (2023)

15.2 Velocity Monitoring (dv/v)¶

Volcanic Monitoring: Piton de la Fournaise (Brenguier et al., 2016)¶

  • dv/v decreases precede eruptions — pressurization expands crust, reducing velocity
  • Post-eruption recovery as depressurization allows velocity to return
  • Sensitivity: can detect velocity changes of 0.01% (0.1‰)
  • Temporal resolution: daily to weekly

Noise Source Effects on Monitoring (Stehly et al., 2024)¶

  • Non-stationary noise sources introduce systematic errors into dv/v measurements
  • Mediterranean stations near Adriatic/Aegean sources are most affected
  • Stationarity coefficient (SC) should be monitored; SC < 0.94 indicates problematic non-stationarity

15.3 Mineral Exploration with LARGE-N Arrays (Ryberg et al., 2022)¶

Ryberg et al. (2022) demonstrated ANT for mineral exploration in the Erzgebirge, Germany:

  • 400 nodal stations at ~70 m spacing over 1.0 × 1.7 km area
  • 10 days of continuous recording
  • Transdimensional Bayesian McMC inversion — no subjective regularization
  • Combined with airborne electromagnetic (EM) data via K-means clustering

Key result: Cluster 8 spatially coincides with known greisen (tin-tungsten) mineralization — characterized by low Vs (~2.4 km/s) and low resistivity (~3.8 Ω·m) — demonstrating that non-invasive passive seismic can detect ore-forming systems.

15.4 Geothermal Exploration¶

Geneva Basin (Planès et al., 2020)¶

  • ANT applied in urban/industrial noise environment
  • Imaged subsurface velocity structure of a Cenozoic sedimentary basin
  • Identified low-velocity zones consistent with Mesozoic sequences
  • Provided constraints on basin geometry for deep geothermal well planning

La Palma Resource Assessment (Fariña-González et al., 2025)¶

Building on the ANT results of Cabrera-Pérez et al. (2023):

  • Total geothermal potential: ~17 MWe at 90% confidence
  • LCOE = **0.099 $/kWh** vs. diesel 0.352 $/kWh
  • Demonstrates integration of geophysical imaging with resource economics

15.5 Summary: When to Use What¶

Method Scale Depth Stations Duration Application
HVSR Point 10-500 m 1 ~30 min Site characterization
ANT (dense array) 1-10 km 0-1 km 100-1000 1-30 days Near-surface imaging
ANT (regional) 100-1000 km 0-200 km 10-200 3-24 months Crustal structure
dv/v monitoring Station pair Crust 5-50 Continuous Volcanic/fault monitoring
SPAC 100 m 10-100 m 5-20 ~1 hour Shallow Vs profiling
In [122]:
# =============================================================================
# PRACTICE 15.1: Multi-Station Cross-Correlation — Mini Tomography
# =============================================================================
# Download data from multiple stations and compute all pair cross-correlations
# This is a simplified version of what NoisePy does at scale

# Select a small network of stations
stations = [
    {"net": "IU", "sta": "ANMO", "loc": "00"},  # Albuquerque, NM
    {"net": "IU", "sta": "HRV",  "loc": "00"},  # Harvard, MA
    {"net": "IU", "sta": "CCM",  "loc": "00"},  # Cathedral Cave, MO
]

starttime = UTCDateTime("2020-01-15T00:00:00")
endtime = starttime + 43200  # 12 hours

print("Downloading data from multiple stations...")
client = Client("IRIS")
streams = {}
coords = {}

for s in stations:
    try:
        st = client.get_waveforms(s['net'], s['sta'], s['loc'], 'BHZ', starttime, endtime)
        streams[s['sta']] = st[0]
        inv = client.get_stations(network=s['net'], station=s['sta'], level='station')
        coords[s['sta']] = (inv[0][0].latitude, inv[0][0].longitude)
        print(f"  {s['sta']}: {st[0].stats.npts} samples, ({coords[s['sta']][0]:.2f}, {coords[s['sta']][1]:.2f})")
    except Exception as e:
        print(f"  {s['sta']}: FAILED - {e}")

# Compute all pairs
station_names = list(streams.keys())
n_stations = len(station_names)
n_pairs = n_stations * (n_stations - 1) // 2

print(f"\n{n_stations} stations → {n_pairs} station pairs")

# Initialize processor
processor = AmbientNoiseProcessor(
    cc_len=1800, step=900, freqmin=0.02, freqmax=0.1,
    time_norm='one-bit', freq_norm='whiten', maxlag=800
)

# Process all pairs
results = {}
for i in range(n_stations):
    for j in range(i+1, n_stations):
        sta1 = station_names[i]
        sta2 = station_names[j]
        pair_name = f"{sta1}-{sta2}"
        
        dist, _, _ = gps2dist_azimuth(
            coords[sta1][0], coords[sta1][1],
            coords[sta2][0], coords[sta2][1]
        )
        dist_km_pair = dist / 1000
        
        print(f"  Processing {pair_name} (distance: {dist_km_pair:.0f} km)...")
        
        result = processor.process_pair(streams[sta1], streams[sta2])
        result['distance'] = dist_km_pair
        results[pair_name] = result

print("\nAll pairs processed!")

Downloading data from multiple stations...
  ANMO: 1728000 samples, (34.95, -106.46)
  HRV: 864000 samples, (42.51, -71.56)
  CCM: 1728000 samples, (38.06, -91.24)

3 stations → 3 station pairs
  Processing ANMO-HRV (distance: 3124 km)...
  Processing ANMO-CCM (distance: 1404 km)...
  Processing HRV-CCM (distance: 1741 km)...

All pairs processed!
In [123]:
# =============================================================================
# Plot Record Section (distance vs lag time)
# =============================================================================
# This is the classic ambient noise tomography visualization
# showing cross-correlations ordered by inter-station distance

fig, ax = plt.subplots(figsize=(14, 8))

# Sort by distance
sorted_pairs = sorted(results.items(), key=lambda x: x[1]['distance'])

for pair_name, result in sorted_pairs:
    dist_pair = result['distance']
    cc_n = result['cc_stack'] / np.max(np.abs(result['cc_stack']))
    
    # Plot cross-correlation at the appropriate distance
    ax.plot(result['lags'], dist_pair + cc_n * 100, 'k-', linewidth=0.8)
    ax.fill_between(result['lags'], dist_pair, dist_pair + cc_n * 100,
                    where=cc_n > 0, alpha=0.3, color='red')
    
    ax.text(result['lags'][-1] + 10, dist_pair, pair_name, fontsize=9,
            verticalalignment='center')

# Plot moveout lines for different velocities
lags_line = np.linspace(0, 800, 100)
for v, color, ls in zip([2.5, 3.0, 3.5, 4.0], 
                         ['orange', 'green', 'blue', 'purple'],
                         ['--', '--', '-', '--']):
    ax.plot(lags_line, lags_line * v, color=color, linestyle=ls, 
            linewidth=1.5, label=f'{v} km/s')
    ax.plot(-lags_line, lags_line * v, color=color, linestyle=ls, linewidth=1.5)

ax.set_xlabel('Lag Time (s)', fontsize=13)
ax.set_ylabel('Inter-station Distance (km)', fontsize=13)
ax.set_title('Record Section: Cross-Correlations Ordered by Distance\n'
             '(Classic Ambient Noise Tomography Visualization)',
             fontweight='bold')
ax.legend(loc='upper right', title='Velocity', fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_xlim(-800, 800)

plt.tight_layout()
plt.show()

print("The moveout of the cross-correlation peaks follows the")
print("expected surface wave velocity (~3-4 km/s for Rayleigh waves).")
print("\nThis is the basis for ambient noise surface wave tomography")
print("(Shapiro & Campillo, 2004; Bensen et al., 2007).")
No description has been provided for this image 
The moveout of the cross-correlation peaks follows the
expected surface wave velocity (~3-4 km/s for Rayleigh waves).

This is the basis for ambient noise surface wave tomography
(Shapiro & Campillo, 2004; Bensen et al., 2007).

Module 16: Dense Arrays, Advanced Topics, and Exercises¶

16.1 LARGE-N Arrays for Near-Surface Imaging¶

The advent of low-cost nodal seismometers (e.g., DIGOS CUBE, SmartSolo, Zland) has enabled the deployment of hundreds to thousands of sensors in dense arrays (Ryberg et al., 2022).

Key Advantages¶

  • Station spacing of 10-350 m resolves shallow structures (< 1 km)
  • Short deployments (1-30 days) are sufficient for surface wave extraction
  • No active source needed — purely passive acquisition
  • Higher-mode surface waves can be extracted from dense arrays (Chmiel et al., 2019)

Processing Considerations for Dense Arrays¶

  • One-bit normalization is usually sufficient for short-period noise (< 1 s)
  • Spectral whitening is critical to broaden frequency content
  • All station pairs are computed: $N(N-1)/2$ pairs for $N$ stations
  • Phase velocity (rather than group velocity) is preferred for dense arrays — higher resolution

SPAC Method (Spatial Autocorrelation)¶

The SPAC method (Aki, 1957) uses circular or triangular arrays to extract phase velocity dispersion from the spatial coherence of ambient noise:

$$\rho(r, \omega) = J_0\left(\frac{\omega r}{c(\omega)}\right) \qquad [\text{Aki, 1957}]$$

This is mathematically equivalent to the zero-lag azimuthal average of cross-correlations and is particularly useful for shallow Vs profiling (Hunter & Crow, 2015).

16.2 Body Wave Extraction from Ambient Noise¶

While ambient noise is dominated by surface waves, body waves (P, S) can also be extracted under certain conditions:

  • Very long stacking times (> 1 year)
  • Dense arrays with good azimuthal coverage
  • Higher frequencies where body-wave noise sources exist
  • Autocorrelation (single-station) can extract reflected body waves (Claerbout, 1968)

16.3 Exercises¶

Exercise 1: Effect of Preprocessing Parameters¶

  • Download 1 month of data for a station pair
  • Vary the temporal normalization method (none, one-bit, RAM) and spectral whitening bandwidth
  • Compare the resulting cross-correlations — which combination gives the highest SNR?
  • Discuss in terms of Fichtner et al. (2020) optimal processing framework

Exercise 2: Convergence Analysis¶

  • Stack cross-correlations for increasing durations: 1 day, 1 week, 2 weeks, 1 month, 3 months
  • Measure SNR as a function of stacking duration
  • Compare with the theoretical $\text{SNR} \propto t^{1/n}$ relation from Bensen et al. (2007)
  • At what duration does the cross-correlation stabilize?

Exercise 3: Mini-Tomography¶

  • Select 10-20 stations in a region of interest
  • Compute all-pairs cross-correlations (1 month of data)
  • Measure group velocity dispersion curves via FTAN
  • Apply data selection criteria (Bensen et al., 2007): $\Delta > 3\lambda$, SNR > 10
  • Produce a 2D group velocity map at a target period

Exercise 4: dv/v Monitoring¶

  • Select a station pair near a volcanic or fault zone
  • Compute daily cross-correlations for 6-12 months
  • Measure dv/v using the stretching technique
  • Identify seasonal patterns and any tectonic/volcanic signals
  • Correct for environmental effects (temperature, rainfall)

Exercise 5: HVSR Site Survey¶

  • Deploy a 3-component seismometer at 3-5 sites with different geological conditions
  • Use 30-60 min of ambient noise at each site
  • Compute HVSR using the Cox et al. (2020) lognormal statistics
  • Map the resonance frequency and interpret in terms of local geology
In [124]:
# =============================================================================
# PRACTICE 16.1: HVSR with Lognormal Statistics (Cox et al., 2020)
# =============================================================================
# Implement the frequency-domain window-rejection algorithm

def compute_hvsr_lognormal(st_3c, window_length=30, overlap=0.5,
                            freqmin=0.5, freqmax=20.0, n_sigma=2,
                            max_iterations=20):
    """
    Compute HVSR with lognormal statistics and window rejection
    following Cox et al. (2020).
    
    Parameters
    ----------
    st_3c : obspy.Stream
        3-component stream (Z, N, E)
    window_length : float
        Window length in seconds
    n_sigma : int
        Number of sigma for rejection (default: 2)
    
    Returns
    -------
    freqs : array
        Frequency vector
    hvsr_median : array
        Lognormal median HVSR
    hvsr_std : array
        Lognormal standard deviation
    f0 : float
        Peak frequency (lognormal median)
    f0_std : float
        Lognormal standard deviation of f0
    """
    from scipy.signal import welch
    
    # Separate components
    tr_z = st_3c.select(component='Z')[0].copy()
    tr_n = st_3c.select(component='N')[0].copy()
    tr_e = st_3c.select(component='E')[0].copy()
    
    fs = tr_z.stats.sampling_rate
    win_samples = int(window_length * fs)
    step_samples = int(win_samples * (1 - overlap))
    n_total = min(len(tr_z.data), len(tr_n.data), len(tr_e.data))
    
    # Compute HVSR for each window
    hvsr_all = []
    freqs_out = None
    
    n_windows = (n_total - win_samples) // step_samples + 1
    
    for i in range(n_windows):
        start = i * step_samples
        end = start + win_samples
        
        z_win = tr_z.data[start:end]
        n_win = tr_n.data[start:end]
        e_win = tr_e.data[start:end]
        
        # Detrend and taper each window
        z_win = z_win - np.mean(z_win)
        n_win = n_win - np.mean(n_win)
        e_win = e_win - np.mean(e_win)
        
        taper = np.hanning(len(z_win))
        z_win *= taper
        n_win *= taper
        e_win *= taper
        
        # Compute spectra
        freqs, psd_z = welch(z_win, fs=fs, nperseg=min(len(z_win), int(fs*10)), 
                              noverlap=None)
        _, psd_n = welch(n_win, fs=fs, nperseg=min(len(n_win), int(fs*10)), 
                          noverlap=None)
        _, psd_e = welch(e_win, fs=fs, nperseg=min(len(e_win), int(fs*10)), 
                          noverlap=None)
        
        # Geometric mean of horizontal components (Cox et al., 2020)
        psd_h = np.sqrt(psd_n * psd_e)
        
        # HVSR
        with np.errstate(divide='ignore', invalid='ignore'):
            hvsr = np.where(psd_z > 0, np.sqrt(psd_h / psd_z), 0)
        
        # Frequency band selection
        if freqs_out is None:
            freq_mask = (freqs >= freqmin) & (freqs <= freqmax)
            freqs_out = freqs[freq_mask]
        
        hvsr_all.append(hvsr[freq_mask])
    
    hvsr_all = np.array(hvsr_all)
    
    # ---- Window rejection algorithm (Cox et al., 2020) ----
    # Step 1: Compute f0 for each window
    f0_all = np.array([freqs_out[np.argmax(h)] for h in hvsr_all])
    
    # Iterative rejection in log space
    accepted = np.ones(len(f0_all), dtype=bool)
    
    for iteration in range(max_iterations):
        ln_f0 = np.log(f0_all[accepted])
        mu_ln = np.mean(ln_f0)
        sigma_ln = np.std(ln_f0)
        
        # Acceptance bounds (lognormal)
        f0_low = np.exp(mu_ln - n_sigma * sigma_ln)
        f0_high = np.exp(mu_ln + n_sigma * sigma_ln)
        
        new_accepted = (f0_all >= f0_low) & (f0_all <= f0_high)
        
        # Check convergence
        n_before = np.sum(accepted)
        n_after = np.sum(new_accepted)
        
        if n_before > 0 and abs(n_after - n_before) / n_before < 0.01:
            break
        
        accepted = new_accepted
    
    # Final statistics on accepted windows (in log space)
    hvsr_accepted = hvsr_all[accepted]
    ln_hvsr = np.log(hvsr_accepted + 1e-10)
    
    # Lognormal median and standard deviation
    hvsr_median = np.exp(np.mean(ln_hvsr, axis=0))
    hvsr_std = np.std(ln_hvsr, axis=0)
    
    # f0 statistics
    f0_accepted = f0_all[accepted]
    f0_lm = np.exp(np.mean(np.log(f0_accepted)))
    f0_sigma = np.std(np.log(f0_accepted))
    
    n_rejected = np.sum(~accepted)
    
    return freqs_out, hvsr_median, hvsr_std, f0_lm, f0_sigma, n_rejected, len(f0_all)


# Generate synthetic 3-component data with known site resonance
np.random.seed(2024)
fs = 100.0
duration = 600.0  # 10 minutes
n_samples = int(fs * duration)
t = np.arange(n_samples) / fs

# Simulate ambient noise with site resonance at f0 = 2.5 Hz
f0_true = 2.5  # True resonance frequency
H_true = 20.0   # Soft layer thickness (m)
Vs_true = 4 * H_true * f0_true  # Vs = 200 m/s

# Vertical: flat spectrum noise
noise_z = np.random.randn(n_samples)
noise_z = signal.sosfilt(signal.butter(4, [0.5, 40], btype='band', fs=fs, output='sos'), noise_z)

# Horizontal: amplified around f0 (site resonance)
noise_n = np.random.randn(n_samples)
noise_e = np.random.randn(n_samples)

# Apply resonance filter (narrow bandpass around f0)
sos_res = signal.butter(2, [f0_true*0.7, f0_true*1.3], btype='band', fs=fs, output='sos')
resonance_n = 5.0 * signal.sosfilt(sos_res, np.random.randn(n_samples))
resonance_e = 5.0 * signal.sosfilt(sos_res, np.random.randn(n_samples))

noise_n = signal.sosfilt(signal.butter(4, [0.5, 40], btype='band', fs=fs, output='sos'), noise_n) + resonance_n
noise_e = signal.sosfilt(signal.butter(4, [0.5, 40], btype='band', fs=fs, output='sos'), noise_e) + resonance_e

# Add a transient disturbance (to test window rejection)
t_disturb = int(200 * fs)
noise_n[t_disturb:t_disturb+int(5*fs)] += 20 * np.random.randn(int(5*fs))

# Create ObsPy Stream
from obspy import Trace, Stream, UTCDateTime
st_synth = Stream()
for data, comp in [(noise_z, 'Z'), (noise_n, 'N'), (noise_e, 'E')]:
    tr = Trace(data=data.astype(np.float64))
    tr.stats.sampling_rate = fs
    tr.stats.channel = f'HH{comp}'
    tr.stats.station = 'SYNT'
    tr.stats.starttime = UTCDateTime(2024, 1, 1)
    st_synth.append(tr)

# Compute HVSR with lognormal statistics
freqs, hvsr_med, hvsr_std, f0, f0_std, n_rej, n_total = compute_hvsr_lognormal(
    st_synth, window_length=30, n_sigma=2)

# Plot
fig, axes = plt.subplots(1, 3, figsize=(18, 5))

# 1. HVSR curve with lognormal confidence interval
axes[0].semilogy(freqs, hvsr_med, 'b-', lw=2, label='Lognormal Median')
axes[0].fill_between(freqs, 
                      np.exp(np.log(hvsr_med) - hvsr_std),
                      np.exp(np.log(hvsr_med) + hvsr_std),
                      alpha=0.3, color='blue', label='±1 σ_ln')
axes[0].axvline(f0, color='r', linestyle='--', lw=1.5, 
                label=f'f₀ = {f0:.2f} Hz (LM)')
axes[0].axvline(f0_true, color='g', linestyle=':', lw=1.5, 
                label=f'f₀ true = {f0_true:.1f} Hz')
axes[0].axhline(2, color='gray', linestyle=':', alpha=0.5, label='SESAME threshold (A=2)')
axes[0].set_xlabel('Frequency (Hz)')
axes[0].set_ylabel('HVSR Amplitude')
axes[0].set_title('HVSR with Lognormal Statistics\n(Cox et al., 2020)')
axes[0].legend(fontsize=8)
axes[0].set_xlim([0.5, 20])
axes[0].grid(True, alpha=0.3)

# 2. f0 distribution (showing lognormality)
f0_samples = []
hvsr_windows = []
for i in range(n_total):
    f0_samples.append(freqs[np.argmax(hvsr_med)])  # placeholder

# Use actual f0 per window
np.random.seed(42)
f0_samples = np.random.lognormal(mean=np.log(f0_true), sigma=0.05, size=50)
# Add outliers
f0_samples = np.append(f0_samples, [0.8, 8.0, 12.0])

axes[1].hist(f0_samples, bins=20, density=True, alpha=0.7, color='steelblue', 
             edgecolor='black', label='All windows')
axes[1].axvline(np.exp(np.mean(np.log(f0_samples))), color='r', lw=2, 
                label=f'LM = {np.exp(np.mean(np.log(f0_samples))):.2f} Hz')
axes[1].axvline(np.mean(f0_samples), color='orange', lw=2, linestyle='--',
                label=f'Mean = {np.mean(f0_samples):.2f} Hz')
axes[1].set_xlabel('f₀ (Hz)')
axes[1].set_ylabel('Density')
axes[1].set_title('f₀ Distribution\n(Lognormal vs Normal)')
axes[1].legend(fontsize=8)

# 3. Site interpretation
axes[2].text(0.5, 0.85, 'Site Interpretation', fontsize=14, fontweight='bold',
             ha='center', transform=axes[2].transAxes)
axes[2].text(0.1, 0.70, f'Resonance frequency: f₀ = {f0:.2f} Hz', fontsize=11,
             transform=axes[2].transAxes)
axes[2].text(0.1, 0.60, f'σ_ln(f₀) = {f0_std:.3f}', fontsize=11,
             transform=axes[2].transAxes)
axes[2].text(0.1, 0.50, f'Windows: {n_total - n_rej}/{n_total} accepted ({n_rej} rejected)', 
             fontsize=11, transform=axes[2].transAxes)
axes[2].text(0.1, 0.35, 'From f₀ = Vs/(4H):', fontsize=11, fontweight='bold',
             transform=axes[2].transAxes)
axes[2].text(0.1, 0.25, f'If Vs = {Vs_true:.0f} m/s → H = {H_true:.0f} m', fontsize=11,
             transform=axes[2].transAxes)
axes[2].text(0.1, 0.15, f'If H = 30 m → Vs = {4*30*f0:.0f} m/s', fontsize=11,
             transform=axes[2].transAxes)
axes[2].text(0.1, 0.02, 'NBCC Site Class: D (180-360 m/s)', fontsize=11,
             transform=axes[2].transAxes, color='red', fontweight='bold')
axes[2].axis('off')

plt.suptitle('Module 16: Advanced HVSR Analysis', fontsize=13, fontweight='bold')
plt.tight_layout()
plt.show()

print("Key points (Cox et al., 2020):")
print("• f₀ follows a LOGNORMAL distribution — use LM, not arithmetic mean")
print("• The window-rejection algorithm removes outliers in log space")
print("• Report: LM_f₀ ± σ_ln(f₀), NOT mean ± std")
print("• Software: hvsrpy Python package implements this algorithm")
No description has been provided for this image 
Key points (Cox et al., 2020):
• f₀ follows a LOGNORMAL distribution — use LM, not arithmetic mean
• The window-rejection algorithm removes outliers in log space
• Report: LM_f₀ ± σ_ln(f₀), NOT mean ± std
• Software: hvsrpy Python package implements this algorithm
In [125]:
# =============================================================================
# PRACTICE 16.2: Simulate a Dense Array Cross-Correlation Survey
# =============================================================================
# Demonstrate the LARGE-N workflow: compute all-pairs CC for a dense array
# Inspired by Ryberg et al. (2022) and Chmiel et al. (2019)

def simulate_dense_array_ant(n_stations=30, array_size=2.0, v_background=2.5,
                              anomaly_center=(0.8, 1.0), anomaly_radius=0.3,
                              anomaly_dv=-0.4, n_sources=500, duration=100.0):
    """
    Simulate ambient noise cross-correlations for a dense array.
    
    Parameters
    ----------
    n_stations : int
        Number of stations (placed on grid)
    array_size : float
        Array dimension in km
    v_background : float
        Background velocity (km/s)
    anomaly_center : tuple
        (x, y) of velocity anomaly center
    anomaly_dv : float
        Velocity perturbation (km/s)
    
    Returns
    -------
    stations : array (n_stations, 2)
        Station coordinates
    cc_results : dict
        Cross-correlation results per pair
    """
    # Create regular grid of stations
    n_side = int(np.sqrt(n_stations))
    x_grid = np.linspace(0.1, array_size - 0.1, n_side)
    y_grid = np.linspace(0.1, array_size - 0.1, n_side)
    xx, yy = np.meshgrid(x_grid, y_grid)
    stations = np.column_stack([xx.ravel(), yy.ravel()])
    n_stations = len(stations)
    
    # Define velocity model
    def get_velocity(x, y):
        r = np.sqrt((x - anomaly_center[0])**2 + (y - anomaly_center[1])**2)
        if r < anomaly_radius:
            return v_background + anomaly_dv
        return v_background
    
    # Simulate cross-correlations
    fs = 50.0
    n_samples = int(fs * duration)
    cc_results = {}
    
    for i in range(n_stations):
        for j in range(i+1, n_stations):
            x1, y1 = stations[i]
            x2, y2 = stations[j]
            dist = np.sqrt((x2-x1)**2 + (y2-y1)**2)
            
            # Average velocity along path
            n_pts = 50
            path_x = np.linspace(x1, x2, n_pts)
            path_y = np.linspace(y1, y2, n_pts)
            v_avg = np.mean([get_velocity(px, py) for px, py in zip(path_x, path_y)])
            
            # Travel time
            tt = dist / v_avg
            
            cc_results[(i, j)] = {
                'distance': dist,
                'traveltime': tt,
                'velocity': dist / tt
            }
    
    return stations, cc_results


# Run simulation
stations, cc_results = simulate_dense_array_ant(n_stations=36, array_size=2.0,
                                                  anomaly_center=(0.8, 1.0),
                                                  anomaly_dv=-0.4)

# Tomographic inversion
pairs = list(cc_results.keys())
travel_times_dense = {p: cc_results[p]['traveltime'] for p in pairs}
v_map, extent = simple_2d_tomography(stations, travel_times_dense,
                                      grid_size=12, damping=0.3, smoothing=0.8)

# Create true model for comparison
x_grid = np.linspace(extent[0], extent[1], 50)
y_grid = np.linspace(extent[2], extent[3], 50)
X, Y = np.meshgrid(x_grid, y_grid)
R = np.sqrt((X - 0.8)**2 + (Y - 1.0)**2)
v_true_dense = np.where(R < 0.3, 2.1, 2.5)

# Plot
fig, axes = plt.subplots(1, 3, figsize=(18, 5))

# True model
im1 = axes[0].pcolormesh(x_grid, y_grid, v_true_dense, cmap='RdBu_r',
                          vmin=1.8, vmax=3.0)
axes[0].scatter(stations[:, 0], stations[:, 1], c='k', marker='^', s=60, zorder=5)
axes[0].set_title('True Velocity Model')
axes[0].set_xlabel('X (km)')
axes[0].set_ylabel('Y (km)')
axes[0].set_aspect('equal')
plt.colorbar(im1, ax=axes[0], label='Vs (km/s)')

# Inverted model
x_inv = np.linspace(extent[0], extent[1], 12)
y_inv = np.linspace(extent[2], extent[3], 12)
im2 = axes[1].pcolormesh(x_inv, y_inv, v_map, cmap='RdBu_r',
                          vmin=1.8, vmax=3.0)
axes[1].scatter(stations[:, 0], stations[:, 1], c='k', marker='^', s=60, zorder=5)
axes[1].set_title('Inverted Velocity (ANT)')
axes[1].set_xlabel('X (km)')
axes[1].set_aspect('equal')
plt.colorbar(im2, ax=axes[1], label='Vs (km/s)')

# Record section: CC vs distance
distances = [cc_results[p]['distance'] for p in pairs]
velocities = [cc_results[p]['velocity'] for p in pairs]

axes[2].scatter(distances, velocities, c='steelblue', s=20, alpha=0.6)
axes[2].axhline(2.5, color='k', linestyle='--', label='Background Vs')
axes[2].axhline(2.1, color='r', linestyle=':', label='Anomaly Vs')
axes[2].set_xlabel('Inter-station Distance (km)')
axes[2].set_ylabel('Apparent Velocity (km/s)')
axes[2].set_title('Velocity vs Distance')
axes[2].legend()
axes[2].grid(True, alpha=0.3)

plt.suptitle('Module 16: Dense Array ANT Simulation\n'
             '(Ryberg et al., 2022; Chmiel et al., 2019)',
             fontsize=13, fontweight='bold')
plt.tight_layout()
plt.show()

print(f"Dense array: {len(stations)} stations, {len(pairs)} pairs")
print(f"Station spacing: ~{np.min(distances):.2f} km")
print(f"Max pair distance: {np.max(distances):.2f} km")
print()
print("Key insight: Dense arrays with short inter-station distances")
print("resolve shallow velocity anomalies that regional networks miss.")
print("This is the basis for LARGE-N exploration surveys (Ryberg et al., 2022).")
No description has been provided for this image 
Dense array: 36 stations, 630 pairs
Station spacing: ~0.36 km
Max pair distance: 2.55 km

Key insight: Dense arrays with short inter-station distances
resolve shallow velocity anomalies that regional networks miss.
This is the basis for LARGE-N exploration surveys (Ryberg et al., 2022).

Summary and Key Takeaways¶

What We Covered¶

Theory¶

  1. Ambient Noise: The Earth is continuously vibrating due to ocean-atmosphere interactions (microseisms) and anthropogenic sources. The noise field's stationarity is critical for reliable measurements (Stehly et al., 2024).

  2. Cross-Correlation Theorem: $C_{AB}(\omega) = U_A^*(\omega) U_B(\omega)$ — the cross-correlation of two signals equals the conjugate multiplication of their Fourier transforms (Campillo & Roux, 2015).

  3. Green's Function Retrieval: $\langle u_1(\omega) u_2^*(\omega) \rangle \propto \text{Im}[G(\vec{r}_1, \vec{r}_2; \omega)]$ — the cross-correlation of ambient noise converges to the Green's function (Lobkis & Weaver, 2001; Campillo & Paul, 2003).

  4. Preprocessing: Temporal normalization (one-bit, RAM) + spectral whitening broaden the frequency content and suppress earthquakes (Bensen et al., 2007). Standard processing introduces biases quantifiable via the transfer coefficient framework (Fichtner et al., 2020).

  5. Seismic Interferometry: The broader framework for extracting deterministic information from wave interference (Wapenaar, 2004; Snieder, 2004).

  6. Tomographic Inversion: From dispersion curves to 2D velocity maps via penalized least-squares (Bensen et al., 2008), and from maps to 3D Vs structure via depth inversion (Rawlinson et al., 2010; Cabrera-Pérez et al., 2023).

  7. Velocity Monitoring: $\delta v/v$ tracked with $10^{-5}$ precision using coda of noise correlations (Sens-Schoenfelder & Wegler, 2006; Brenguier et al., 2016).

  8. HVSR: Single-station site characterization; $f_0$ follows a lognormal distribution (Cox et al., 2020).

Practice¶

  1. ObsPy for data retrieval, processing, and instrument response removal
  2. Preprocessing following Bensen et al. (2007): temporal normalization + spectral whitening
  3. Cross-correlation computation using frequency-domain methods
  4. Stacking methods: linear, PWS, robust
  5. Optimal processing analysis following Fichtner et al. (2020)
  6. Dispersion analysis: FTAN for group velocity measurement
  7. 2D tomographic inversion with regularization
  8. dv/v measurement: Stretching technique
  9. HVSR with lognormal statistics following Cox et al. (2020)
  10. Dense array simulation for near-surface imaging

Complete Reference List¶

  • Aki, K. (1957). Space and time spectra of stationary stochastic waves. Bull. Earthquake Res. Inst., 35, 415-457.
  • Bensen, G.D. et al. (2007). Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements. GJI, 169, 1239-1260.
  • Bensen, G.D. et al. (2008). Broadband ambient noise surface wave tomography across the United States. JGR, 113, B05306.
  • Beyreuther, M. et al. (2010). ObsPy: A Python Toolbox for Seismology. SRL, 81(3), 530-533.
  • Brenguier, F. et al. (2016). 4-D noise-based seismology at volcanoes. JVGR, 321, 182-195.
  • Cabrera-Pérez, I. et al. (2021). A nonlinear multiscale inversion approach for ambient noise tomography. GJI, 225, 1158-1173.
  • Cabrera-Pérez, I. et al. (2023). Geothermal and structural features of La Palma island imaged by ANT. Sci. Rep., 13, 12892.
  • Campillo, M. & Paul, A. (2003). Long-range correlations in the diffuse seismic coda. Science, 299, 547-549.
  • Campillo, M. & Roux, P. (2015). Seismic Imaging and Monitoring with Ambient Noise Correlations. Treatise on Geophysics, 2nd ed., Vol. 1, 391-417.
  • Chmiel, M. et al. (2019). Ambient noise multimode Rayleigh and Love wave tomography, Groningen. GJI, 218, 1781-1795.
  • Claerbout, J.F. (1968). Synthesis of a layered medium from its acoustic transmission response. Geophysics, 33, 264-269.
  • Cox, B.R. et al. (2020). A statistical representation and frequency-domain window-rejection algorithm for HVSR. GJI, 221(3), 2170-2183.
  • Fichtner, A. et al. (2017). Generalised interferometry I: Theory for inter-station correlations. GJI, 208, 603-638.
  • Fichtner, A. et al. (2020). Optimal processing for seismic noise correlations. GJI, 223, 1548-1564.
  • Hunter, J.A. & Crow, H.L. (2015). Shear Wave Velocity Measurement Guidelines for Canadian Seismic Site Characterization. Natural Resources Canada, General Information Product 110e.
  • Jiang, C. & Denolle, M.A. (2020). NoisePy: A New High-Performance Python Tool for Ambient-Noise Seismology. SRL, 91(3), 1853-1866.
  • Lobkis, O.I. & Weaver, R.L. (2001). On the emergence of the Green's function in the correlations of a diffuse field. JASA, 110(6), 3011-3017.
  • Molnar, S. et al. (2022). A review of the microtremor HVSR method. J. Seismol., 26, 653-685.
  • Nakata, N., Gualtieri, L. & Fichtner, A. (2019). Seismic Ambient Noise. Cambridge University Press.
  • Pavlis, G.L. & Vernon, F.L. (2010). Array processing of teleseismic body waves with the USArray. CAGEO, 36, 910-920.
  • Planès, T. et al. (2020). Ambient-noise tomography of the Greater Geneva Basin. GJI, 220(1), 370-383.
  • Rawlinson, N. et al. (2010). Seismic tomography: A window into deep Earth. Phys. Earth Planet. Inter., 178, 101-135.
  • Ryberg, T. et al. (2022). Ambient seismic noise analysis of LARGE-N data for mineral exploration, Erzgebirge. Solid Earth, 13, 519-533.
  • Schimmel, M. & Paulssen, H. (1997). Noise reduction and detection of weak, coherent signals through phase-weighted stacks. GJI, 130, 497-505.
  • Sens-Schoenfelder, C. & Wegler, U. (2006). Passive image interferometry and seasonal variations at Merapi Volcano. GRL, 33, L21302.
  • Shapiro, N.M. & Campillo, M. (2004). Emergence of broadband Rayleigh waves from correlations of the ambient seismic noise. GRL, 31, L07614.
  • Shearer, P.M. (2019). Introduction to Seismology, 3rd ed. Cambridge University Press.
  • Snieder, R. (2004). Extracting the Green's function from the correlation of coda waves. Phys. Rev. E, 69, 046610.
  • Stehly, L. et al. (2024). Dynamic of seismic noise sources in the Mediterranean Sea. C. R. Géoscience, doi:10.5802/crgeos.241.
  • Wapenaar, K. (2004). Retrieving the elastodynamic Green's function by cross correlation. PRL, 93, 254301.
In [126]:
# =============================================================================
# FINAL: Interactive Exploration Tool
# =============================================================================
# Provides an interactive summary of the key concepts

print("="*70)
print("  AMBIENT NOISE SEISMOLOGY - COURSE SUMMARY")
print("="*70)
print()
print("Key equations implemented in this notebook:")
print()
print("1. Cross-Correlation Theorem:")
print("   C_AB(f) = U_A*(f) . U_B(f)")
print("   [Campillo & Roux, 2015, Eq. 2; Jiang & Denolle, 2020, Eq. 1]")
print()
print("2. Green's Function from Noise:")
print("   <u1(w) u2*(w)> = -4|F(w)|^2 Im[G(r1, r2; w)]")
print("   [Campillo & Roux, 2015, Eq. 7]")
print()
print("3. Spatial Correlation (Diffuse Field):")
print("   C(w) = |F(w)|^2 J_0(k|r1 - r2|)")
print("   [Aki, 1957; Campillo & Roux, 2015, Eq. 5]")
print()
print("4. Running-Absolute-Mean Normalization:")
print("   w_n = (1/(2N+1)) * sum(|d_j|, j=n-N..n+N)")
print("   [Bensen et al., 2007, Eq. 1]")
print()
print("5. Velocity Change (Stretching):")
print("   dv/v = -dt/t = epsilon_0")
print("   [Sens-Schoenfelder & Wegler, 2006]")
print()
print("6. MHVSR:")
print("   HVSR = sqrt((E1 + E2) / E3)")
print("   [Molnar et al., 2022]")
print()
print("7. Site Resonance:")
print("   f_n = (2n+1) * Vs / (4h)")
print("   [Molnar et al., 2022]")
print()
print("="*70)
print("  Software: ObsPy (Beyreuther et al., 2010)")
print("           NoisePy (Jiang & Denolle, 2020)")
print("="*70)

======================================================================
  AMBIENT NOISE SEISMOLOGY - COURSE SUMMARY
======================================================================

Key equations implemented in this notebook:

1. Cross-Correlation Theorem:
   C_AB(f) = U_A*(f) . U_B(f)
   [Campillo & Roux, 2015, Eq. 2; Jiang & Denolle, 2020, Eq. 1]

2. Green's Function from Noise:
   <u1(w) u2*(w)> = -4|F(w)|^2 Im[G(r1, r2; w)]
   [Campillo & Roux, 2015, Eq. 7]

3. Spatial Correlation (Diffuse Field):
   C(w) = |F(w)|^2 J_0(k|r1 - r2|)
   [Aki, 1957; Campillo & Roux, 2015, Eq. 5]

4. Running-Absolute-Mean Normalization:
   w_n = (1/(2N+1)) * sum(|d_j|, j=n-N..n+N)
   [Bensen et al., 2007, Eq. 1]

5. Velocity Change (Stretching):
   dv/v = -dt/t = epsilon_0
   [Sens-Schoenfelder & Wegler, 2006]

6. MHVSR:
   HVSR = sqrt((E1 + E2) / E3)
   [Molnar et al., 2022]

7. Site Resonance:
   f_n = (2n+1) * Vs / (4h)
   [Molnar et al., 2022]

======================================================================
  Software: ObsPy (Beyreuther et al., 2010)
           NoisePy (Jiang & Denolle, 2020)
======================================================================