CODE for MEASURING SEISMIC TRAVEL-TIME CHANGES with the WAVELET METHOD

This python program is based on the work of Mao and Mordret, that you can find here : https://github.com/shujuanmao/dt-wavelet

The original code was in MATLAB, and here is the Python equivalent

Contact : Higueret Quentin (quentin.higueret@univ-grenoble-alpes.fr) and Aurélien Mordret (aurelien.mordret@univ-grenoble-alpes.fr)

This package contains codes and test data for measuring seismic travel-time shifts in the time-frequency domain using the wavelet cross-spectrum analysis.

Requirements

All the files must be in the same folder when you launch the program

Python 3 and the PyCWT packages from https://github.com/regeirk/pycwt are needed to run the codes

You can use pip to install this package :

$ pip install pycwt

Contents in this package

• xwt.py
  
  The core function to calculate dt in the time-frequency domain by wavelet cross-spectrum analysis.

• plotting_example.py
  
  An example of using the xwt function on synthetic data. Plots come with one click.

• ori_waveform.npy
  new_waveform.npy
  
  Two synthetic waveforms for testing the codes, ori_waveform and new_waveform
  The synthetic seismograms are generated using velocity models by a homogeneous background superimposed with random heterogeneities.
  The perturbation between the current and reference velocity models is a 0.05% homogeneous dv/v throughout the medium. (If interested, see Section 3.1 in the following reference for more details.)

• time.npy
  fs.npy
  
  The time vector and sampling frequency associated with the synthetic waveforms

Usage

WXamp, WXspec, WXangle, Wcoh, WXdt, freqs, coi = xwt(trace_ref,trace_current,fs,ns,nt,vpo,freqmin,freqmax,nptsfreq)

• Input
  
  trace_ref,trace_current : Two vectors, reference and current time series.
  fs : Sampling Frequency // Positive scalar, sampling frequency.
  ns : NumScalesToSmooth // Positive integer, indicating the length of boxcar window.
  nt : DegTimeToSmooth // Positive scalar, indicating the length of the Gaussian window.
  vpo : VoicesPerOctave // Even integer from 4 to 48, indicates how fine the frequency is discretized. Should be no less than 10.
  freqmin : The starting value of the frequency vector, in Hz.
  freqmax : The ending value of the frequency vector, in Hz.
  nptsfreq : Number of frequency samples to generate between the starting and ending value.

• Output
  
  WXamp : Matrix of amplitude product of two CWT in time-frequency domain.
  WXspec : Complex-valued matrix, the wavelet cross-spectrum.
  WXangle : Matrix of the angle of the complex argument in WXspec.
  Wcoh: Matrix of wavelet coherence.
  WXdt : Matrix of time difference and phase difference, respectively between the two input time series in time-frequency domain.
  freqs : Vector of frequencies used in CWT, in Hz.
  coi : Cone of influce. Vector of the maximum period of useful information at each arrival time.
  Note that here coi means MAX PERIODS, in Matlab codes it refers to MIN FREQUENCY.

Plotting Example

Using the xwt function and the synthetic data, we can perform a cross-wavelet transform.
This is the plot we obtain using the plotting_example.py program :

(If interested, compare the image obtain here with the Figure 3 in this article Geophysical Journal International, Volume 221, Issue 1, April 2020, Pages 550–568, https://doi.org/10.1093/gji/ggz495)

Reference

S.Mao, A.Mordret, M.Campillo, H.Fang, R.D.van der Hilst, (2019), On the Measurement of Seismic Travel-Time Changes in the Time-Frequency Domain with Wavelet Cross-Spectrum Analysis, GJI.

Torrence, C. and Compo, G. P.. A Practical Guide to Wavelet Analysis.