Comparing approaches to measuring seismic phase variations in the time, frequency, and wavelet domains

Introduction

The shallow subsurface is a dynamic environment, responding to tectonic stresses, changes in hydrologic and magmatic conditions, thermoelastic and tidal stresses, and anthropogenic activities. Seismic wave velocities respond to changes in physical properties in the subsurface, so monitoring velocity changes provides a powerful technique we can use to characterize the forces acting on the shallow Earth. In particular, the coda waves of seismic signals broadly sample the medium and are thus sensitive to very small changes in subsurface properties. Assuming a similar source-receiver path, changes in subsurface properties correspond to changes in the phase of the coda waves. By measuring the phase changes in coda waves through time, we can monitor the physical properties of the shallow subsurface. The reliability and accuracy of this monitoring are determined by our ability to measure phase changes in these signals. A detailed description of the following experiments and results can be found in our manuscript:

Yuan, C., Bryan, J., and Denolle, M. Comparing appraoches to measuring phase variations in time, frequency, wavelet domains, submitted to Geophysical Journal International.

To understand how changes in seismic velocities can be measured from phase shifts in seismic waves, consider a ray that travels a distance L between two seismometers in time t at velocity v. These quantities are related via

For a homogeneous velocity change, the ray path does not change (dL=0). Taking the differential of L, we find

This yields a relation between the velocity change, dv, and the change in arrival time, dt

We find that the change in arrival time scales linearly with lag time, resulting in a stretching or compression of the waveform following a decrease or an increase, respectively, in seismic velocity.

The simple linear relationship between change in arrival time and change in velocity may become more complicated when we consider a heterogeneous velocity perturbation. The observed phase changes then depend on the spatial variation of the velocity perturbation (lateral and depth localization). In order to localize the velocity perturbation in 3D, it is useful to measure the phase changes for many station pairs and at many frequencies ("spectrum of dv/v").

There are two main ways we can measure waveform distortion, and thus the change in subsurface velocity.

1. Global stretching factor: linearly stretch/compress two coda waveforms such that you maximize their correlation coefficient.

2. Local phase shifts: extract the time shifts at first for either local windows or individual phase lags and then approximate the global stretching/compression with the linear regression of these time shifts with lag time.

We can categorize existing methods for measuring velocity perturbations by the domain in which they are applied.

• Time-Domain: Windowed Cross-Correlation (WCC), Trace-Stretching (TS), and Dynamic Time-Warping (DTW).

• Frequency-Domain: Moving Window Cross-Spectrum (MWCS)

• Wavelet-Domain: Wavelet Cross-Spectrum (WXS)

We add two new wavelet-domain methods to this collection.

• Wavelet-Transform Stretching (WTS)

• Wavelet-Transform Dynamic Time-Warping (WTDTW)

For the rest of this tutorial, we will focus on comparing the performance of these methods through a variety of modeling exercises. Our roadmap for the tutorial is organized as follows:

1. Half-space model: We apply a homogeneous velocity change. The goal of this exercise is to validate our assumption that a homogeneous velocity change corresponds to a linear stretch of the coda waves.

2. Layer over half-space model: We apply three different velocity perturbations with different spatial localization. The goals of these tests are to explore the role of the background velocity structure, and to see the first-order effects of depth-localized velocity changes.

3. Sliding layer model: We perturb the velocity in a thin layer, which we slide from the surface to specific depths. The goal of this test is to explore whether the dv/v spectrum systematically responds to the depth of a localized velocity perturbation.

Environment Setup

Before we can measure anything, we must first install the necessary packages

using Pkg

Pkg.rm("AbstractPlotting");

Pkg.rm("WGLMakie");

Pkg.update()

# dependencies

Pkg.add("GLM");

Pkg.add("FFTW");

Pkg.add("JLD2");

Pkg.add("DSP");

Pkg.add("Random");

Pkg.add("Dierckx");

Pkg.add("DataFrames")

# seismic processing packages

Pkg.add("SeisIO");

Pkg.add("SeisNoise");

# dv/v package

Pkg.add(PackageSpec(url="https://github.com/jaredbryan881/SeisDvv.jl.git"));

We then import the packages

using Plots, JLD2, DSP, FFTW

using SeisIO, SeisNoise

using SeisDvv

And finally we open the file containing waveform data

# open up the file to get synthetic waveforms

f = JLD2.jldopen(
dvv_validation.jld2);

Half-Space Model

The first test we will look at is for a homogeneous velocity change in a half-space model. The purpose of this section is to evaluate whether a homogeneous velocity change corresponds to a linear stretch of the coda waves.

We read a base waveform and a perturbed waveform, which is generated by applying a velocity perturbation of magnitude 0.1% everywhere in the medium.

# read waveforms from the first model configuration: half-space medium with homogeneous velocity perturbation

halfspace = f["halfspace"]

time = halfspace["t"]; # time axis

S1 = halfspace["base"]; # unperturbed velocity structure

S2 = halfspace["perturbed"]; # homogeneously perturbed velocity structure

# find the indices corresponding to 5-35 s

window = findall(t->t<=35.0 && t>=5.0, time);

# true dv/v

dvv = 0.1;

We can visualize the effect of this 0.1% velocity perturbation in the time domain. Make note of the increased separation between the "reference" and "current" waveforms with increasing lag time.

# slide a window across the waveform we use to show the change in phase with increasing lag

mag_len = 400

anim = @animate for i∈1:length(time[window])-mag_len

  mag_win=i:i+mag_len

  # Plot coda section of the signals

  p1=Plots.plot(time[window], S1[window], color=:blue, xaxis="Time [s]", labels="Reference")

  Plots.plot!(time[window], S2[window], color=:red, labels="Current")

  # Plot the magnification window

  Plots.vspan!(p1, [time[window[mag_win[1]]], time[window[mag_win[end]]]], color=:red, alpha=0.1, labels="Magnified Section")

  # Plot the magnified section

  p2 = Plots.plot(time[window[mag_win]], S1[window[mag_win]], color=:blue, xaxis="Time [s]", labels="Reference")

  Plots.plot!(time[window[mag_win]], S2[window[mag_win]], color=:red, labels="Current", ylims=(-0.00075, 0.00075))

  # Plot

  p=Plots.plot(p1, p2, layout=(2,1), legend=:outertopright)

  end every 5

gif(anim, fps=15)

We can also look at the change in the signals in the frequency domain (via their power spectra). The scaling theorem of the Fourier transform tells us that a linear stretch of a signal in the time domain corresponds to a linear compression of its Fourier transform.

"""

Compute the power spectrum of a given signal

# Arguments

- `A::AbstractArray`: Input signal

# Returns

-  `power::AbstractArray`: Power spectrum

"""

function power(A::AbstractArray, fs::AbstractFloat)

  power=abs.(rfft(A)).^2

  freqs=FFTW.rfftfreq(length(A), fs)

	return power, freqs

end

# Find the power spectra of S1 and S2

P1,freqs = power(S1[window], 100.0);

P2,freqs = power(S2[window], 100.0);

# plot against frequencies in fft

p=Plots.plot(freqs[1:100],P1[1:100], lw=2, color=:blue, labels="Reference", xaxis="Frequency [Hz]");

Plots.plot!(freqs[1:100],P2[1:100], lw=2, color=:red, labels="Current");

display(p)

Time-Domain Methods (WCC, TS, DTW)

Seeing the difference between the waveforms is one thing, but now we can quantify it with our first dv/v measurements! Let's start with the time-domain methods.

# parameters for dv/v algorithms

fs=100.0     # sample frequency

win_len=5.0  # sliding window length

win_step=2.5 # sliding window step

fmin=0.5     # frequency minimum

fmax=3.0     # frequency maximum

# windowed cross-correlation (WCC)

println("Calculating dv/v with WCC")

dvv_wcc, dvv_err_wcc, int_wcc, int_err_wcc, dvv0_wcc, dvv0_err_wcc = wcc(S1, S2, time, window, fs, win_len, win_step);

# trace stretching (TS)

println("Calculating dv/v with TS")

dvv_ts, cc_ts, cdp_Ts, eps_ts, err_ts, allC_ts = SeisDvv.stretching(S1, S2, time, window, fmin, fmax, dvmin=-0.005, dvmax=0.005, ntrial=100);

println("Correlation coefficient: $(cc_ts)")

# dynamic time warping (DTW)

println("Calculating dv/v with DTW")

dvv_dtw, dvv_err_dtw, int_dtw, int_err_dtw, dvv0_dtw, dvv0_err_dtw = dtw(S1, S2, time, window, fs);

Linearly stretching the waveform gives a very high correlation coefficient! Now if stretching also recovers the true dv/v value, we can be confident that a homogeneous velocity change corresponds to linearly stretching the waveform.

p=Plots.scatter([0, 1, 2], [dvv0_wcc, dvv_ts, dvv0_dtw], xticks=(collect(0:3), ["WCC", "TS", "DTW"]), label=nothing, color=:black, markersize=5)

Plots.plot!([0,2], [dvv, dvv], color=:black, line=:dash, ylims=(0,0.15), label="True dv/v", lw=2., xaxis="Method", yaxis="dv/v [%]")

display(p)

Both TS and DTW accurately recover the true velocity perturbation of 0.1%, and WCC only does slightly worse. TS and DTW are both global measurements that estimate a stretching factor. Since we assumed (and validated) that a homogeneous velocity change corresponds to a linear stretch of the waveforms, we should not be surprised t hat TS and DTW perform very well.

Frequency Bands (MWCS, WTS, WTDTW, WXS)

We can also apply the dv/v methods in a set of frequency bands. This is desirable because coda waves of different frequencies will be sensitive to different depths. In cases where the velocity perturbation is localized in depth, frequency dependence of the dv/v measurements can provide valuable information about the depth of the velocity change. In this case of homogeneous velocity change, we should expect to see the perturbation at all frequencies.

fbands = [[0.6 1.0]; [0.8 1.2]; [1.0 1.6]; [1.2 2.4]; [2.0 3.0]]

(nbands,_) = size(fbands)

# MWCS

println("Calculating dv/v with MWCS")

dvv_mwcs = zeros(nbands)

dvv0_mwcs = zeros(nbands)

for iband=1:nbands

  t_axis_mwcs, dt_mwcs, error_mwcs, mcoh_mwcs = SeisDvv.mwcs(S1, S2, fbands[iband,1], fbands[iband,2], fs, time[window[1]], win_len, win_step, 0);

  dvv_mwcs[iband], dvv_err_mwcs, int_mwcs, int_err_mwcs, dvv0_mwcs[iband], dvv0_err_mwcs = SeisDvv.mwcs_dvv(t_axis_mwcs, dt_mwcs, error_mwcs, mcoh_mwcs, "static", 0.0, 0.0, time[window[1]], time[window[end]]-time[window[1]], "right");

end

# Wavelet transform stretching (WTS)

println("Calculating dv/v with WTS")

freqbands, dvv_wts, cc, cdp, ε, err, allC = SeisDvv.stretching(S1, S2, time, window, fbands, dvmin=-0.005, dvmax=0.005, ntrial=100);

# Wavelet transform dynamic time warping (WTDTW)

println("Calculating dv/v with WTDTW")

freqbands, dvv_wtdtw, dvv_err_wtdtw, int_wtdtw, int_err_wtdtw, dvv0_wtdtw, dvv0_err_wtdtw = dtw(S1, S2, time, window, fs, fbands);

# Wavelet cross-spectrum (WXS)

println("Calculating dv/v with WXS")

freqbands, dvv_wxs, dvv_err_wxs, int_wxs, int_err_wxs, dvv0_wxs, dvv0_err_wxs = wxs(S1, S2, time, window, fs, fbands);

p=Plots.scatter(collect(1:nbands), -dvv_mwcs*100, label="MWCS")

Plots.scatter!(collect(1:nbands), dvv_wts, label="WTS")

Plots.scatter!(collect(1:nbands), dvv0_wtdtw, label="WTDTW")

Plots.scatter!(collect(1:nbands), dvv0_wxs, label="WXS")

Plots.plot!([1, 5], [dvv, dvv], label="True dv/v", xaxis="Frequency Band", yaxis="dv/v [%]", legend=:bottomright)

display(p)

Here we see that all methods perform fairly well in the first four frequency bands, but WXS performs poorly in the final frequency band, where there is the least power. WTS and WTDTW remain accurate over all frequencies, and MWCS is fairly accurate, though tends to overestimate the dv/v, for all frequencies.

All Frequencies (WTS, WTDTW, WXS)

One feature of the wavelet transform is its high frequency resolution. We can measure dv/v at a range of discrete frequencies, corresponding to a perfectly narrow bandpass filter at that frequency. This is the "spectrum of dv/v."

fwindow = [0.5, 3.0]

# WTS

println("Calculating dv/v with WTS")

(freqbands, dvv_wts, err_wts) = SeisDvv.stretching(S1, S2, time, window, fwindow, ntrial=1000);

# WTDTW

println("Calculating dv/v with WTDTW")

(freqbands, dvv_wtdtw, err_wtdtw, a, ea, m0, em0) = dtw(S1, S2, time, window, fs, fwindow);

# WXS

println("Calculating dv/v with WXS")

(freqbands, dvv_wxs, err_wxs) = wxs(S1, S2, time, window, fs, fwindow);

p=Plots.scatter(freqbands[:,1], dvv_wts, label="WTS")

Plots.scatter!(freqbands[:,1], dvv_wtdtw, label="WTDTW")

Plots.scatter!(freqbands[:,1], dvv_wxs, label="WXS")

Plots.plot!([0.6, 3], [dvv, dvv], line=:dash, color=:black, label="True dv/v", xaxis="Frequency [Hz]", yaxis="dv/v [%]")

display(p)

The patterns from the previous test with frequency bands are accentuated by this group of methods. WXS is unstable for low-energy signals, resulting in under-prediction with increasing frequency. We can also see inaccuracies in WTS and WTDTW due to variations in signal strength, with WTDTW having up to 0.005% positive bias at the spectral peak at 1 Hz.

Layer Over Half-Space Model

We can now move to a more complicated and realistic velocity model. The goal of this section is to examine whether the background velocity structure has an effect on the dv/v results. For direct comparison, we will first apply a homogeneous velocity perturbation to the layer over half-space model and calculate the resulting dv/v. We will then perturb the upper layer only, and then perturb only the lower half-space. This will allow us to observe the first-order effects of depth-localized velocity changes. For each of these tests, we will also look at two source-receiver configurations: colocated source and receiver (approximation for the ambient noise auto-correlation and yields the reflectivity response) and a distant source-receiver pair (approximation for the ambient noise cross-correlation and yields the interstation Green's function).

Let's begin by reading the data for each layer perturbation (upper, lower, both) and each source-receiver configuration (auto-correlation, cross-correlation). To reduce the number of plots, we will look only at the "spectrum of dv/v," but dv/v could be computed in frequency bands or with no frequency decomposition as in the first example.

lohs = f["layeroverhalfspace"];

time_ac = lohs["acorr/t"];

time_xc = lohs["xcorr/t"];

# unperturbed velocities

S1_ac = lohs["acorr/base"];

S1_xc = lohs["xcorr/base"];

# perturbed velocities

# upper layer perturbed

S2_ac_up = lohs["acorr/upper"];

S2_xc_up = lohs["xcorr/upper"];

# lower layer perturbed

S2_ac_low = lohs["acorr/lower"];

S2_xc_low = lohs["xcorr/lower"];

# whole space perturbed

S2_ac_both = lohs["acorr/both"];

S2_xc_both = lohs["xcorr/both"];

We will use different time and frequency limits for the auto-correlation and cross-correlation cases. The larger separation in the cross-correlation case results in longer travel times and higher attenuation of high-frequency content.

# time window

window_ac = findall(t->t<=60.0 && t>=20.0, time_ac);

window_xc = findall(t->t<=240.0 && t>=180.0, time_xc);

# frequency window

fwindow_ac = [0.5, 5.0];

fwindow_xc = [0.01, 1.0];

# true dv/v

dvv = 0.1;

Auto-correlation

We now compute dv/v for each layer perturbation for the auto-correlation source-receiver configuration.

# WTS

println("Calculating dv/v with WTS")

(freqbands, dvv_wts_up, err_wts_up) = SeisDvv.stretching(S1_ac, S2_ac_up, time_ac, window_ac, fwindow_ac, dvmin=-0.005, dvmax=0.005, ntrial=100);

(freqbands, dvv_wts_low, err_wts_low) = SeisDvv.stretching(S1_ac, S2_ac_low, time_ac, window_ac, fwindow_ac, dvmin=-0.005, dvmax=0.005, ntrial=100);

(freqbands, dvv_wts_both, err_wts_both) = SeisDvv.stretching(S1_ac, S2_ac_both, time_ac, window_ac, fwindow_ac, dvmin=-0.005, dvmax=0.005, ntrial=100);

# WTDTW

println("Calculating dv/v with WTDTW")

(freqbands, dvv_wtdtw_up, err_wtdtw_up, a, ea, m0, em0) = dtw(S1_ac, S2_ac_up, time_ac, window_ac, fs, fwindow_ac);

(freqbands, dvv_wtdtw_low, err_wtdtw_low, a, ea, m0, em0) = dtw(S1_ac, S2_ac_low, time_ac, window_ac, fs, fwindow_ac);

(freqbands, dvv_wtdtw_both, err_wtdtw_both, a, ea, m0, em0) = dtw(S1_ac, S2_ac_both, time_ac, window_ac, fs, fwindow_ac);

# WXS

println("Calculating dv/v with WXS")

(freqbands, dvv_wxs_up, err_wxs_up) = wxs(S1_ac, S2_ac_up, time_ac, window_ac, fs, fwindow_ac);

(freqbands, dvv_wxs_low, err_wxs_low) = wxs(S1_ac, S2_ac_low, time_ac, window_ac, fs, fwindow_ac);

(freqbands, dvv_wxs_both, err_wxs_both) = wxs(S1_ac, S2_ac_both, time_ac, window_ac, fs, fwindow_ac);

Let's visualize these results:

p_both=Plots.scatter(freqbands[:,1], dvv_wts_both, label="WTS")

Plots.scatter!(freqbands[:,1], dvv_wtdtw_both, label="WTDTW")

Plots.scatter!(freqbands[:,1], dvv_wxs_both, label="WXS")

Plots.plot!([0.5, 5], [dvv, dvv], label="True dv/v", yaxis="dv/v [%]", ylims=(0.0,0.15), title="Whole Space", color=:black)

p_up=Plots.scatter(freqbands[:,1], dvv_wts_up, label="WTS")

Plots.scatter!(freqbands[:,1], dvv_wtdtw_up, label="WTDTW")

Plots.scatter!(freqbands[:,1], dvv_wxs_up, label="WXS")

Plots.plot!([0.5, 5], [dvv, dvv], label="True dv/v", yaxis="dv/v [%]", ylims=(0.0,0.15),title="Upper Layer", color=:black)

p_low=Plots.scatter(freqbands[:,1], dvv_wts_low, label="WTS")

Plots.scatter!(freqbands[:,1], dvv_wtdtw_low, label="WTDTW")

Plots.scatter!(freqbands[:,1], dvv_wxs_low, label="WXS")

Plots.plot!([0.5, 5], [dvv, dvv], label="True dv/v", xaxis="Frequency [Hz]", yaxis="dv/v [%]", ylims=(0.0,0.15), title="Lower Layer", color=:black)

p=Plots.plot(p_both, p_up, p_low, layout=(3,1))

When both layers are perturbed (homogeneous velocity change), we see an identical pattern as in the first example. WTS and WTDTW accurately recover the true dv/v for all frequencies, and WXS under-predicts dv/v at increasing frequency due to loss of signal strength.

We now see the interesting contrast between dv/v spectra with shallow vs deep perturbations. The depth of the velocity perturbation has a strong effect on the shape of the dv/v spectrum! When the shallow layer is perturbed, the dv/v values are under-predicted at low frequencies, with WTS and WTDTW approaching the true dv/v value as frequency increases. When the deep layer is perturbed, the dv/v values are closes to the true value for low frequencies, but the under-prediction worsens as frequency increases. In the view that coda waves are composed of surface waves, this is easily attributed to the frequency-dependent depth sensitivity. Shallow perturbations will be seen by high frequencies, and deep perturbations will be seen by low frequencies.

It is also interesting to note that in contexts where the dv/v is localized in depth, we essentially always underpredict the true value of dv/v. This is likely because the scattered waves that comprise the coda sample both the perturbed and unperturbed medium. The competing influence of the unperturbed and perturbed layers result in dv/v estimates that lie between 0% and the true dv/v value.

Cross-Correlation

We now move to the case of a distant source and receiver, which approximates the ambient noise cross-correlation.

# WTS

println("Calculating dv/v with WTS")

(freqbands, dvv_wts_up, err_wts_up) = SeisDvv.stretching(S1_xc, S2_xc_up, time_xc, window_xc, fwindow_xc, dvmin=-0.005, dvmax=0.005, ntrial=100);

(freqbands, dvv_wts_low, err_wts_low) = SeisDvv.stretching(S1_xc, S2_xc_low, time_xc, window_xc, fwindow_xc, dvmin=-0.005, dvmax=0.005, ntrial=100);

(freqbands, dvv_wts_both, err_wts_both) = SeisDvv.stretching(S1_xc, S2_xc_both, time_xc, window_xc, fwindow_xc, dvmin=-0.005, dvmax=0.005, ntrial=100);

# WTDTW

println("Calculating dv/v with WTDTW")

(freqbands, dvv_wtdtw_up, err_wtdtw_up, a, ea, m0, em0) = dtw(S1_xc, S2_xc_up, time_xc, window_xc, fs, fwindow_xc);

(freqbands, dvv_wtdtw_low, err_wtdtw_low, a, ea, m0, em0) = dtw(S1_xc, S2_xc_low, time_xc, window_xc, fs, fwindow_xc);

(freqbands, dvv_wtdtw_both, err_wtdtw_both, a, ea, m0, em0) = dtw(S1_xc, S2_xc_both, time_xc, window_xc, fs, fwindow_xc);

# WXS

println("Calculating dv/v with WXS")

(freqbands, dvv_wxs_up, err_wxs_up) = wxs(S1_xc, S2_xc_up, time_xc, window_xc, fs, fwindow_xc);

(freqbands, dvv_wxs_low, err_wxs_low) = wxs(S1_xc, S2_xc_low, time_xc, window_xc, fs, fwindow_xc);

(freqbands, dvv_wxs_both, err_wxs_both) = wxs(S1_xc, S2_xc_both, time_xc, window_xc, fs, fwindow_xc);

We can visualize these dv/v results for each layer perturbation:

p_both=Plots.scatter(freqbands[:,1], dvv_wts_both, label="WTS")

Plots.scatter!(freqbands[:,1], dvv_wtdtw_both, label="WTDTW")

Plots.scatter!(freqbands[:,1], dvv_wxs_both, label="WXS")

Plots.plot!([0.01, 1.0], [dvv, dvv], label="True dv/v", yaxis="dv/v [%]", ylims=(0.0,0.25), title="Whole Space", color=:black)

p_up=Plots.scatter(freqbands[:,1], dvv_wts_up, label="WTS")

Plots.scatter!(freqbands[:,1], dvv_wtdtw_up, label="WTDTW")

Plots.scatter!(freqbands[:,1], dvv_wxs_up, label="WXS")

Plots.plot!([0.01, 1.0], [dvv, dvv], label="True dv/v", yaxis="dv/v [%]", ylims=(0.0,0.25),title="Upper Layer", color=:black)

p_low=Plots.scatter(freqbands[:,1], dvv_wts_low, label="WTS")

Plots.scatter!(freqbands[:,1], dvv_wtdtw_low, label="WTDTW")

Plots.scatter!(freqbands[:,1], dvv_wxs_low, label="WXS")

Plots.plot!([0.01, 1.0], [dvv, dvv], label="True dv/v", xaxis="Frequency [Hz]", yaxis="dv/v [%]", ylims=(0.0,0.25), title="Lower Layer", color=:black)

p=Plots.plot(p_both, p_up, p_low, layout=(3,1))

We find essentially the same pattern in the distant source-receiver configuration. The methods accurately recover the dv/v for almost all frequencies in the case of homogeneous velocity perturbation, and show clear frequency-dependence for the cases with localized velocity perturbations. The upper layer perturbation is most accurately recovered at high frequencies, and the lower layer perturbation is best recovered at low frequencies.

Sliding Layer Model

The previous tests demonstrate that depth-localization has a strong effect on the shape of the spectrum of dv/v. The goal of this section is to systematically explore this effect. We perturb a 200m-thick layer starting at the surface, and we increase the depth to the perturbed layer by 100m until we reach 1 km.

Because we are creating dv/v spectra for many different models, we will focus only on WTS in order to prevent the figures from getting too cluttered. You could repeat the same analysis with WXS or WTDTW.

We begin by reading the unperturbed waveforms.

ds = f["indepth"];

time_ac = ds["acorr/t"];

time_xc = ds["xcorr/t"];

# unperturbed velocities

S1_ac = ds["acorr/base"];

S1_xc = ds["xcorr/base"];

As with the previous exercise, we use different time and frequency limits due to the travel time and attenuation differences between the auto-correlation and cross-correlation exercises.

# time window

window_ac = findall(t->t<=6.5 && t>=2.5, time_ac);

window_xc = findall(t->t<=225.0 && t>=185.0, time_xc);

# frequency window

fwindow_ac = [0.5, 30.0];

fwindow_xc = [0.05, 1.5];

Here we define the depth of the perturbations. This gives 11 distinct simulations.

# depth array

depths_ac = collect(0.0:0.1:1.0); # just take the first 1 km since all methods lose sensitivity deeper

depths_xc = collect(0.0:1.0:10.0); # just take the first 10 km since all methods lose sensitivity deeper

Auto-Correlation

Before we look at the spectrum of dv/v for each layer, we calculate dv/v with the raw (broadband) signals.

# initialize dv/v array before iterating over depths

dvv_ts_ac = zeros(length(depths_ac));

for i=1:length(depths_ac)

  println("Computing dv/v spectrum at $(depths_ac[i]) km")

  # get the waveforms corresponding to current depth perturbation

  S2_ac = ds["acorr/$(i-1)"];

  # compute dv/v

  dvv_ts_ac[i], cc_ts, cdp_ts, eps_ts, err_ts, allC_ts = SeisDvv.stretching(S1_ac, S2_ac, time_ac, window_ac, fwindow_ac[1], fwindow_ac[end], dvmin=-0.005, dvmax=0.005, ntrial=100);

end

We can visualize these dv/v results, coloring the points for use in the dv/v spectra figures latter.

colors = range(HSL(colorant"red"), stop=HSL(colorant"green"), length=length(depths_ac));

p=Plots.scatter(depths_ac, dvv_ts_ac, color=colors, ms=6)

Plots.plot!([0.0, 1.0], [0.1, 0.1], label="True dv/v", xaxis="Depth [km]", yaxis="dv/v [%]", ylims=(0.0,0.12), color=:black, legend=false)

display(p)

The recovered velocity perturbation always under-predicts the true value, and the sensitivity to the perturbed layer rapidly decreases with depth. This rapid loss of sensitivity with depth suggests that dv/v measurements made with broadband coda waves are likely inaccurate when the perturbation is localized in depth.

Let's now find the spectrum of dv/v as a function of depth.

# initialize dv/v array of dimension [n_depths, n_frequencies] before iterating over depths

dvv_ts_ac = zeros(length(depths_ac), 71);

freqbands = zeros(71)

for i=1:length(depths_ac)

  println("Computing dv/v spectrum at $(depths_ac[i]) km")

  # get the waveforms corresponding to current depth perturbation

  S2_ac = ds["acorr/$(i-1)"];

  # compute dv/v

  fbands, dvv_ts_ac[i,:], cc_ts_ac, cdp_ts_ac, eps_ts_ac, err_ts_ac, allC_ts_ac = SeisDvv.stretching(S1_ac, S2_ac, time_ac, window_ac, fwindow_ac, dvmin=-0.005, dvmax=0.005, ntrial=100);

  freqbands[:] = fbands[:,1]

end

We plot the results with the same coloring as used in the broadband dv/v measurements

colors = range(HSL(colorant"red"), stop=HSL(colorant"green"), length=length(depths_ac));

p=Plots.plot([fwindow_ac[1], fwindow_ac[end]], [0.1, 0.1], label="True dv/v", xaxis="Frequency [Hz]", yaxis="dv/v [%]", ylims=(0.0,0.12), color=:black, legend=false, xscale=:log10);

# plot dv/v spectrum for each perturbation depth

for i=1:length(depths_ac)

  Plots.scatter!(freqbands, dvv_ts_ac[i,:], color=colors[i], xscale=:log10)

end

display(p)

We see that the dv/v measurements systematically vary with the depth of the perturbation. For shallow velocity perturbations, the spectrum behaves similarly to the "upper layer perturbation" models of the previous section such that the dv/v estimates approach the true value as frequency increases.

As the depth of the perturbation increases, the dv/v spectra have global maxima at decreasing frequency. In addition, as the depth to the perturbation increases, the maximum in the dv/v spectra decreases in amplitude (under-prediction worsens). The loss of sensitivity to deep perturbations at high frequencies results from the fact that high frequencies are most sensitive to the shallow, unperturbed layer. Loss of sensitivity at low frequencies is likely due to large averaging over the unperturbed medium.

We also observe a "limit frequency" below which all dv/v spectra converge. Tests of layer thickness, dominant source frequency, and lapse time show that this limit frequency (1.38 Hz) depends on the lapse time selection. This is likely due to the increasing contribution of body waves to the coda as lapse time increases. Body waves sample deeper structures in the medium, resulting in decreases limit frequency.

Cross-Correlation

We now repeat the previous analysis for the source-receiver configuration that approximates the ambient noise cross-correlation. In this case, we perturb a 4 km thick layer, which we slide in depth with a 1 km increment.

We begin by measuring dv/v with the broadband signals.

# initialize dv/v array before iterating over depths

dvv_ts_xc = zeros(length(depths_xc));

for i=1:length(depths_xc)

  println("Computing dv/v spectrum at $(depths_xc[i]) km")

  # get the waveforms corresponding to current depth perturbation

  S2_xc = ds["xcorr/$(i-1)"];

  # compute dv/v

  dvv_ts_xc[i], cc_ts, cdp_ts, eps_ts, err_ts, allC_ts = SeisDvv.stretching(S1_xc, S2_xc, time_xc, window_xc, fwindow_xc[1], fwindow_xc[end], ntrial=1000);

end

We plot and color the points as in the previous example.

p=Plots.scatter(depths_xc, dvv_ts_xc, color=colors, ms=6)

Plots.plot!([0.0, 10.0], [dvv, dvv], label="True dv/v", xaxis="Depth [km]", yaxis="dv/v [%]", ylims=(0.0,0.12), color=:black, legend=false)

display(p)

The results for this test are similar to those from the auto-correlation tests. The dv/v estimates are closest to the true value for shallow perturbations, and the sensitivity to the perturbation rapidly decreases with increasing perturbation depth.

Now let's find the spectrum of dv/v as a function of perturbation depth.

# initialize dv/v array of dimension [n_depths, n_frequencies] before iterating over depths

dvv_ts_xc = zeros(length(depths_xc), 59);

freqbands = zeros(59)

for i=1:length(depths_xc)

  println("Computing dv/v spectrum at $(depths_xc[i]) km")

  # get the waveforms corresponding to current depth perturbation

  S2_xc = ds["xcorr/$(i-1)"];

  # compute dv/v

  fbands, dvv_ts_xc[i,:], dvv_err_ts_xc = SeisDvv.stretching(S1_xc, S2_xc, time_xc, window_xc, fwindow_xc, dvmin=-0.005, dvmax=0.005, ntrial=100);

  freqbands[:] = fbands[:,1]

end

Plotting with the colors from the broadband dv/v:

p=Plots.plot([fwindow_xc[1], fwindow_xc[end]], [dvv, dvv], label="True dv/v", xaxis="Frequency [Hz]", yaxis="dv/v [%]", ylims=(0.0,0.12), color=:black, legend=false)

# plot dv/v spectrum for each perturbation depth

for i=1:length(depths_xc)

  Plots.scatter!(freqbands, dvv_ts_xc[i, :], color=colors[i])

end

display(p)

Similar to the results from the auto-correlation configuration, the spectra obtained from shallow perturbations approach the true value as frequency increases. For increasingly deep perturbations, the global maximum of the spectra decrease in amplitude and occur at decreasing frequency.

Conclusions

• The depth of the velocity perturbation imparts clear structure on the dv/v spectrum.

• The measured dv/v always, often largely, under-predicts the true velocity perturbation.

• The maximum value and corresponding frequency of the global maximum of the dv/v spectrum decrease with increasing perturbation depth.

• WTS is our recommended method for measurement of dv/v spectra.