Large-N cross-correlation functions of ambient seismic noise recorded in the Vienna basin, Austria.

This repository contains the C₁cross-correlation functions of the ambient seismic field and station metadata required to reproduce the results described in the manuscript "Source effects in higher-order ambient seismic field correlations" by Schippkus et al., 2025.

We provide a complete set of Jupyter notebooks that implements all processing described in the manuscript for full reproducibility. They make use of the data provided in this repository. The notebooks are hosted on GitHub at https://github.com/schipp/higher_order_correlations_c2.

The files hosted here are:

<strong>stations_receivers.csv</strong>: List of 1990 stations used as receiver stations. Station OMV.GDT is used as the master station in the manuscript. Station names are arbitrary.
<strong>stations_auxiliary.csv</strong>: List of 304 stations used as auxiliary stations.
<strong>correlations_for_c1_data.pt</strong>: C₁ cross-correlations of all 1990 receiver stations with the master station OMV.GDT in the center. Saved as a PyTorch torch.tensor with shape [1990, 3001]. Sampling rate is 5 Hz, correlation functions are limited to 300 seconds of anti-causal and causal lapse time (=&gt; 3001 samples). First dimension (the receiver stations) is sorted alphabetically by station name.
<strong>correlations_for_c2_data.pt</strong>: C₁ cross-correlations of all 1990 receiver stations with the 304 auxiliary stations surrounding them. Saved as a PyTorch torch.tensor with shape [1990, 304, 3001]. Sampling rate is 5 Hz, correlation functions are limited to 300 seconds of anti-causal and causal lapse time (=&gt; 3001 samples). First dimension (the receiver stations) and second dimension (the auxiliary stations) are sorted alphabetically by station name.

The cross-correlations are computed as described in the manuscript: ~4 weeks of continuous recordings are cut into 1-hr windows and spectrally whitened. All windows are cross-correlated and linearly stacked. No additional processing is applied. For more details on these data, please see the manuscript.

We provide a minimal Python code snippet to load the C₁ correlation functions with the master station ("OMV.GDT"), read station locations, filter the correlation functions with a narrowband-filter around 0.3 Hz, and plot the focal spot of the C₁ correlation wavefield in space:

<pre><code class="language-python">import torch
import polars as pl
import matplotlib.pyplot as plt
from scipy.signal import butter, filtfilt

# load correlation functions
sampling_rate = 5
lapse_times = torch.arange(-300, 300 + 1 / sampling_rate, 1 / sampling_rate)
corrs = torch.load("correlations_for_c1_data.pt", weights_only=False)

# load metadata
stations = pl.read_csv("stations_receivers.csv")
names = stations.select("station").to_series().to_list()
x, y = stations.select(["X", "Y"]).to_numpy().T
master_station = "OMV.GDT"

# apply acausal narrowband filter
frequency_band = (0.29, 0.31)
a, b = butter(4, frequency_band, btype="bandpass", fs=sampling_rate)
filtered_corrs = torch.tensor(filtfilt(a, b, corrs, axis=-1).copy())

# extract focal spot and set master station (=auto-correlation) to zero
focal_spot = filtered_corrs[:, torch.argmin(torch.abs(lapse_times))]
focal_spot[names.index(master_station)] = 0

# plot focal spot
fig, ax = plt.subplots()
ax.scatter(
x,
y,
c=focal_spot,
cmap="RdBu_r",
vmin=-focal_spot.abs().max(),
vmax=focal_spot.abs().max(),
)
ax.set(
xlabel="Distance (km)",
ylabel="Distance (km)",
title="Focal spot",
aspect="equal"
)</code></pre>

本仓库包含Schippkus等人2025年发表的论文《高阶环境地震场互相关中的源效应》中所述研究结果复现所需的环境地震场C₁互相关函数，以及台站元数据。

我们提供了完整的Jupyter Notebook集合，实现了论文中所述的全部处理流程以确保结果可完全复现。这些Notebook依托本仓库提供的数据运行，托管于GitHub平台：https://github.com/schipp/higher_order_correlations_c2。

本仓库托管的文件包括：
**stations_receivers.csv**：作为接收台站的1990个台站列表。论文中将OMV.GDT台站作为主台站，台站名称为任意标识。
**stations_auxiliary.csv**：作为辅助台站的304个台站列表。
**correlations_for_c1_data.pt**：以主台站OMV.GDT为中心的1990个接收台站与该主台站之间的C₁互相关函数。保存为PyTorch张量（torch.tensor），形状为[1990, 3001]。采样率为5 Hz，相关函数覆盖300秒的反因果与因果走时范围，对应3001个采样点。第一维度（接收台站）按台站名称字母顺序排序。
**correlations_for_c2_data.pt**：1990个接收台站与其周围304个辅助台站之间的C₁互相关函数。保存为PyTorch张量，形状为[1990, 304, 3001]。采样率为5 Hz，相关函数覆盖300秒的反因果与因果走时范围，对应3001个采样点。第一维度（接收台站）与第二维度（辅助台站）均按台站名称字母顺序排序。

本互相关函数的计算流程如论文所述：将约4周的连续记录切割为1小时时长的窗口，并进行频谱白化处理。对所有窗口执行互相关后进行线性叠加，未施加额外处理流程。如需了解这些数据的更多细节，请参阅该论文。

我们提供了一段极简Python代码片段，用于加载包含主台站("OMV.GDT")的C₁互相关函数、读取台站位置信息、以0.3 Hz为中心的窄带滤波器对相关函数进行滤波，并绘制C₁相关波场的空间焦斑：

<pre><code class="language-python">import torch
import polars as pl
import matplotlib.pyplot as plt
from scipy.signal import butter, filtfilt

# load correlation functions
sampling_rate = 5
lapse_times = torch.arange(-300, 300 + 1 / sampling_rate, 1 / sampling_rate)
corrs = torch.load("correlations_for_c1_data.pt", weights_only=False)

# load metadata
stations = pl.read_csv("stations_receivers.csv")
names = stations.select("station").to_series().to_list()
x, y = stations.select(["X", "Y"]).to_numpy().T
master_station = "OMV.GDT"

# apply acausal narrowband filter
frequency_band = (0.29, 0.31)
a, b = butter(4, frequency_band, btype="bandpass", fs=sampling_rate)
filtered_corrs = torch.tensor(filtfilt(a, b, corrs, axis=-1).copy())

# extract focal spot and set master station (=auto-correlation) to zero
focal_spot = filtered_corrs[:, torch.argmin(torch.abs(lapse_times))]
focal_spot[names.index(master_station)] = 0

# plot focal spot
fig, ax = plt.subplots()
ax.scatter(
x,
y,
c=focal_spot,
cmap="RdBu_r",
vmin=-focal_spot.abs().max(),
vmax=focal_spot.abs().max(),
)
ax.set(
xlabel="Distance (km)",
ylabel="Distance (km)",
title="Focal spot",
aspect="equal"
)</code></pre>