Ambient Noise Cross-Correlation — Exercises¶

In [ ]:

Copied!

!pip install obspy cartopy -q
!pip install obspy cartopy -q

In [ ]:

Copied!





import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from obspy import UTCDateTime
from obspy.clients.fdsn import Client
from obspy.signal.cross_correlation import correlate
from obspy.geodetics import locations2degrees, degrees2kilometers
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from obspy import UTCDateTime
from obspy.clients.fdsn import Client
from obspy.signal.cross_correlation import correlate
from obspy.geodetics import locations2degrees, degrees2kilometers

Setup: Load Stations and Helper Function¶

Run the cells below to set up the Taiwan station network and the helper function from the lecture.

In [ ]:

Copied!





# Connect to FDSN data center and query Taiwan stations
client = Client("IRIS")

inventory = client.get_stations(
    minlatitude=22.1, maxlatitude=26.0,
    minlongitude=119.6, maxlongitude=121.8,
    starttime=UTCDateTime("2020-01-01"),
    endtime=UTCDateTime("2020-01-02"),
    channel="BH*",
    level="station",
)

stations = {}
for network in inventory:
    for station in network:
        name = f"{network.code}.{station.code}"
        stations[name] = (station.latitude, station.longitude)

print(f"Loaded {len(stations)} stations")
# Connect to FDSN data center and query Taiwan stations
client = Client("IRIS")

inventory = client.get_stations(
    minlatitude=22.1, maxlatitude=26.0,
    minlongitude=119.6, maxlongitude=121.8,
    starttime=UTCDateTime("2020-01-01"),
    endtime=UTCDateTime("2020-01-02"),
    channel="BH*",
    level="station",
)

stations = {}
for network in inventory:
    for station in network:
        name = f"{network.code}.{station.code}"
        stations[name] = (station.latitude, station.longitude)

print(f"Loaded {len(stations)} stations")

In [ ]:

Copied!





def download_and_correlate(sta1, sta2, starttime, endtime, max_lag=300):
    """Download noise for two stations, preprocess, and cross-correlate."""
    net1, name1 = sta1.split(".")
    net2, name2 = sta2.split(".")

    st1 = client.get_waveforms(net1, name1, "", "BHZ", starttime, endtime)
    st2 = client.get_waveforms(net2, name2, "", "BHZ", starttime, endtime)

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

    for tr in [tr1, tr2]:
        tr.detrend("demean")
        tr.detrend("linear")
        tr.taper(0.05)

    if tr1.stats.sampling_rate != tr2.stats.sampling_rate:
        tr2.resample(tr1.stats.sampling_rate)

    t1 = max(tr1.stats.starttime, tr2.stats.starttime)
    t2 = min(tr1.stats.endtime, tr2.stats.endtime)
    tr1.trim(t1, t2)
    tr2.trim(t1, t2)

    shift = int(tr1.stats.sampling_rate * max_lag)
    ccf = correlate(tr1.data, tr2.data, shift=shift, normalize='naive')
    lags = np.arange(-shift, shift + 1) / tr1.stats.sampling_rate

    lat1, lon1 = stations[sta1]
    lat2, lon2 = stations[sta2]
    dist = degrees2kilometers(locations2degrees(lat1, lon1, lat2, lon2))

    return lags, ccf, dist
def download_and_correlate(sta1, sta2, starttime, endtime, max_lag=300):
    """Download noise for two stations, preprocess, and cross-correlate."""
    net1, name1 = sta1.split(".")
    net2, name2 = sta2.split(".")

    st1 = client.get_waveforms(net1, name1, "", "BHZ", starttime, endtime)
    st2 = client.get_waveforms(net2, name2, "", "BHZ", starttime, endtime)

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

    for tr in [tr1, tr2]:
        tr.detrend("demean")
        tr.detrend("linear")
        tr.taper(0.05)

    if tr1.stats.sampling_rate != tr2.stats.sampling_rate:
        tr2.resample(tr1.stats.sampling_rate)

    t1 = max(tr1.stats.starttime, tr2.stats.starttime)
    t2 = min(tr1.stats.endtime, tr2.stats.endtime)
    tr1.trim(t1, t2)
    tr2.trim(t1, t2)

    shift = int(tr1.stats.sampling_rate * max_lag)
    ccf = correlate(tr1.data, tr2.data, shift=shift, normalize='naive')
    lags = np.arange(-shift, shift + 1) / tr1.stats.sampling_rate

    lat1, lon1 = stations[sta1]
    lat2, lon2 = stations[sta2]
    dist = degrees2kilometers(locations2degrees(lat1, lon1, lat2, lon2))

    return lags, ccf, dist

Exercise 1: Cross-Correlate a New Station Pair¶

In the lecture, we used TW.NACB and TW.YHNB. Now try a different pair: TW.SSLB and TW.YULB.

Complete the steps below to download data, preprocess, cross-correlate, and check whether the surface wave arrival matches the expected time.

1a. Download and preprocess¶

Fill in the three preprocessing steps. Recall from the lecture:

"demean": removes the DC offset
"linear": removes slow instrumental drift
taper(0.05): smooths edges to prevent spectral artifacts

In [ ]:

Copied!





# Download one day of data
sta1, sta2 = "TW.SSLB", "TW.YULB"
starttime = UTCDateTime("2020-01-01")
endtime = UTCDateTime("2020-01-02")

net1, name1 = sta1.split(".")
net2, name2 = sta2.split(".")
st1 = client.get_waveforms(net1, name1, "", "BHZ", starttime, endtime)
st2 = client.get_waveforms(net2, name2, "", "BHZ", starttime, endtime)

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

# Preprocessing: fill in the three steps
for tr in [tr1, tr2]:
    tr.detrend(???)    # step 1: remove DC offset
    tr.detrend(???)    # step 2: remove linear drift
    tr.taper(???)      # step 3: taper edges

# Trim to common time window
t_start = max(tr1.stats.starttime, tr2.stats.starttime)
t_end = min(tr1.stats.endtime, tr2.stats.endtime)
tr1.trim(t_start, t_end)
tr2.trim(t_start, t_end)

print(f"Traces ready: {tr1.stats.npts} samples at {tr1.stats.sampling_rate} Hz")
# Download one day of data
sta1, sta2 = "TW.SSLB", "TW.YULB"
starttime = UTCDateTime("2020-01-01")
endtime = UTCDateTime("2020-01-02")

net1, name1 = sta1.split(".")
net2, name2 = sta2.split(".")
st1 = client.get_waveforms(net1, name1, "", "BHZ", starttime, endtime)
st2 = client.get_waveforms(net2, name2, "", "BHZ", starttime, endtime)

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

# Preprocessing: fill in the three steps
for tr in [tr1, tr2]:
    tr.detrend(???)    # step 1: remove DC offset
    tr.detrend(???)    # step 2: remove linear drift
    tr.taper(???)      # step 3: taper edges

# Trim to common time window
t_start = max(tr1.stats.starttime, tr2.stats.starttime)
t_end = min(tr1.stats.endtime, tr2.stats.endtime)
tr1.trim(t_start, t_end)
tr2.trim(t_start, t_end)

print(f"Traces ready: {tr1.stats.npts} samples at {tr1.stats.sampling_rate} Hz")

1b. Cross-correlate and compute expected arrival¶

Fill in the ??? to:

Compute the maximum shift in samples from the max lag in seconds
Compute the expected arrival time using $\tau = D / c$ with $c = 3$ km/s

In [ ]:

Copied!





# Cross-correlate
max_lag = 300  # seconds
shift = ???  # convert max_lag to samples (hint: multiply by sampling rate)
ccf = correlate(tr1.data, tr2.data, shift=shift, normalize='naive')
lags = np.arange(-shift, shift + 1) / tr1.stats.sampling_rate

# Inter-station distance
lat1, lon1 = stations[sta1]
lat2, lon2 = stations[sta2]
dist_km = degrees2kilometers(locations2degrees(lat1, lon1, lat2, lon2))

# Expected surface wave arrival time
v_surface = 3.0  # km/s
t_expected = ???  # hint: distance / velocity

print(f"Distance: {dist_km:.1f} km")
print(f"Expected arrival: {t_expected:.1f} s")

# Plot
plt.figure(figsize=(10, 4))
plt.plot(lags, ccf, 'k', lw=0.8)
plt.axvline(t_expected, color='r', ls='--', alpha=0.7, label=f'Expected ({t_expected:.1f} s)')
plt.axvline(-t_expected, color='r', ls='--', alpha=0.7)
plt.title(f'Cross-Correlation: {sta1} \u2013 {sta2} ({dist_km:.0f} km)')
plt.xlabel('Lag time (seconds)')
plt.ylabel('Correlation coefficient')
plt.xlim(-100, 100)
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# Cross-correlate
max_lag = 300  # seconds
shift = ???  # convert max_lag to samples (hint: multiply by sampling rate)
ccf = correlate(tr1.data, tr2.data, shift=shift, normalize='naive')
lags = np.arange(-shift, shift + 1) / tr1.stats.sampling_rate

# Inter-station distance
lat1, lon1 = stations[sta1]
lat2, lon2 = stations[sta2]
dist_km = degrees2kilometers(locations2degrees(lat1, lon1, lat2, lon2))

# Expected surface wave arrival time
v_surface = 3.0  # km/s
t_expected = ???  # hint: distance / velocity

print(f"Distance: {dist_km:.1f} km")
print(f"Expected arrival: {t_expected:.1f} s")

# Plot
plt.figure(figsize=(10, 4))
plt.plot(lags, ccf, 'k', lw=0.8)
plt.axvline(t_expected, color='r', ls='--', alpha=0.7, label=f'Expected ({t_expected:.1f} s)')
plt.axvline(-t_expected, color='r', ls='--', alpha=0.7)
plt.title(f'Cross-Correlation: {sta1} \u2013 {sta2} ({dist_km:.0f} km)')
plt.xlabel('Lag time (seconds)')
plt.ylabel('Correlation coefficient')
plt.xlim(-100, 100)
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

1c. Measure the actual velocity¶

Find the lag time of the strongest peak in the CCF (at positive lags), then compute the surface wave velocity.

In [ ]:

Copied!





# Find the peak in the positive-lag side of the CCF
positive_mask = lags > 0
positive_lags = lags[positive_mask]
positive_ccf = ccf[positive_mask]

peak_lag = ???  # hint: positive_lags at the index where |positive_ccf| is largest
v_measured = ???  # hint: distance / peak_lag

print(f"Peak lag time: {peak_lag:.1f} s")
print(f"Measured velocity: {v_measured:.2f} km/s")
print(f"Reference velocity: {v_surface:.2f} km/s")
# Find the peak in the positive-lag side of the CCF
positive_mask = lags > 0
positive_lags = lags[positive_mask]
positive_ccf = ccf[positive_mask]

peak_lag = ???  # hint: positive_lags at the index where |positive_ccf| is largest
v_measured = ???  # hint: distance / peak_lag

print(f"Peak lag time: {peak_lag:.1f} s")
print(f"Measured velocity: {v_measured:.2f} km/s")
print(f"Reference velocity: {v_surface:.2f} km/s")

Exercise 2: Interpretation Questions¶

Answer the following questions based on what you learned in the lecture and Exercise 1.

2a. Asymmetry¶

In the lecture moveout plot (Section 5), some station pairs show a stronger signal on one side (positive or negative lag) than the other. What causes this asymmetry? What does it tell you about where the noise is coming from?

Your answer here:

2b. Stacking¶

The lecture showed that individual 3-hour CCFs are noisy but the surface wave signal is still visible. Explain why averaging (stacking) many short-window CCFs improves the signal-to-noise ratio.

Your answer here:

2c. Connection to tomography¶

Suppose you measured the surface wave travel time between every pair of stations in the Taiwan network. How could you use those measurements to create a map of seismic velocity across Taiwan? What determines the spatial resolution of such a map?

Your answer here: