liaohaiyang1534/cross-wavelet-transform

This function provides travel-time shifts with Morlet cross-wavelet transformation

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 ([email protected]) and Aurélien Mordret ([email protected])

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 cwt 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 = 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

Plotting Example

Using the cwt 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)

Coming soon

Implementation of the cone of influence for each frequency bands

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, In Review.

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