Using wavelet decomposition to retrieve information from seismograms

Point of contact: Petros Bogiatzis


Determination of P- and S-wave arrivals is important for a variety of seismological applications, including earthquake identification and seismic tomography. The arrival times are typically picked by dedicated human analysts through visual inspection of the waveforms. However, with the rapid expansion of seismic instrumentation, availability of robust and reliable automatic techniques that can process a large number of seismic traces is becoming essential. In this work, we present a pair of algorithms for efficient picking of P and S onset times. Both algorithms are based on the continuous wavelet transform (Figure 1) of the seismic signals that, unlike Fourier transform, uses basis functions that are localized in time and frequency. The analysis, therefore, considers data in both the time and frequency domains, and it is suitable for the processing of non-stationary signals. For determining the P-wave arrival, the wavelet coefficients are calculated using the vertical component of the seismogram, and the onset time of the wave is identified. In the case of the S-wave arrival, we take advantage of the polarization of the shear waves, and cross-examine the wavelet coefficients from the two horizontal components. We evaluate the performance of our algorithms using real data from Hi-net and F-net (Okada et al., 2004; Obara et al., 2005). We also show that demonstrate that the algorithms can be fully automated to measure additional properties such as P- and S-wave frequency content, dispersion, and shear-wave splitting.


Figure 1: The Continuous Wavelet Transform (CWT) of a signal f, calculates the cross-correlation of the signal with dilated or compressed versions of the wavelet.




Here is a general work flow of P-picking algorithm (Figure 2)

  1. The vertical component of the seismogram is used without being filtered.
  2. A rough analysis of the arrival type is performed to select the appropriate wavelet; Lines are fit to the envelope of the seismogram for two overlapping windows; one that includes the peak amplitude (red) and another that does not (blue). The comparison of their slopes allows to distinguish between sharp and emerging phase arrivals.
  3. CWT is applied to the seismogram. The arrival of the first compressional wave causes changes of the power, phase and frequency content of the seismogram. These changes can readily be detected in the wavelet domain.
  4. To enhance the non-stationary seismic signal over the background noise, two additional steps are taken. The first is the application of a moving 2-D range-filter to the wavelet coefficients. The second step evaluates a simple detection function to each scale (white curves).
    5. The final pick is a weighted average of the picks for each scale. The weighing is based upon the energy and the SNR for each scale.


Figure 2: General work flow of P-picking algorithm




S-waves arrive within the coda of P-waves and other phases. The polarization implies that the horizontal components should record coherent signal when S-wave arrives. To detect this arrival, we utilize the wavelet cross-spectrum (WCS; Hudgins et al., 1993; Maraun & Kurths, 2004). We calculate the WCS as the dot product between the corresponding CWT coefficients of the horizontal components of the seismogram. In the synthetic example above, we use the Ricker wavelet to represent the S-arrival and a phase difference of π between the two horizontal components. The WCS shows a sharp region of high anti-correlation (Figure 3). The detection is performed in a similar way as the P-wave arrival.



Figure 3: General work flow of S-picking algorithm




By keeping track of the individual picks and their weights over different scales, the frequency content of the final pick can be estimated. This data-driven additional information is useful for calculating the sensitivity kernel of the ray-path in seismic tomography. Moreover, body-wave dispersion can be measured directly from the picks along different scales as shown in Figure 4.

Figure 4 (Left column) Top row shows the vertical component of the Hi-net YMKH station around P-wave arrival time. Middle row shows the CWT of the seismogram, and bottom row shows the filtered coefficients. White curves show the variation of the detection function and they are normalized for each scale. White triangles mark the pick at each scale. The P-onset time (vertical line) is calculated as the weighted average of the individual picks. P-wave dispersion, which is also visible in the seismogram, is well captured from the picks along at different scales. (Right column) Dispersion diagram showing the pick time and the corresponding apparent speed as a function of scale and the associated center frequency. The intensity of the circles shows the importance of each pick that is, dark circles are associated with larger weights used to calculate the arrival time.




Toward automating S-wave splitting measurements

If a complex wavelet is used, then both CWT and WCS coefficients are complex. From the amplitude of the coefficients, we can select the high-correlation region that corresponds to the S-wave window. From the phase, we retrieve the time delay between the components at the window selected using the amplitude. We can use this property to automate S-wave splitting measurements (i.e., fast splitting direction and time-delay). To demonstrate this technique we use an example with an earthquake that occurred at the subducting Philippine's slab at the depth of 140 km and the F-net station AMM that is located right above the hypocenter of the event (M4.3, 02/03/2002; Figure 5). This configuration allows S waves to sample the mantle wedge between the subducting and the overriding plate.

Figure 5 Left: Map of Japan. The star shows the epicenter of the selected earthquake and the triangle the location of the F-net station AMM. Right: (a) Horizontal components of F-net AMM station around S-wave arrival time. (b) Amplitude of the WCS coefficients. The coherent region constrains the S-wave window. (c) Phase of the WCS coefficients show the phase-difference between the two signals. (d) Horizontal components after time delay correction inferred from the phase difference at the peak amplitude. (e) WCS amplitude for the corrected components. (f) WCS phase for the corrected components. Note that the two signals are now in phase for scales greater than ~10.




We rotate the horizontal components at various angles, and estimate and correct the delay time, δt (Figures 5 and 6). The fast splitting direction φ, (shaded region) is found at the maximum amplitude that corresponds to 0 phase difference. Our results suggest trench parallel splitting, with φ=48-50o and δt=0.68 sec, which are in good agreement with previous studies for the same station (e.g., Long & van der Hilst, 2005).




Figure 6: Determination of the fast-splitting direction and delay time of S waves recorded at the AMM station. The two vertical lines indicate the maximum coefficient values that correspond to phase difference of 0 (solid line) and π (dashed line), respectively. (a) Peak amplitude before (dashed curve) and after (solid curve) the delay-time correction as a function of rotation angles of the horizontal components. (b) The phase difference at the peak amplitude before (dashed curve) and after (solid curve) the delay-time correction. (c) The delay time correction inferred from the phase difference at various rotation angles.





For more information check

  • Bogiatzis, P. & Ishii, M., 2015. Continuous Wavelet Decomposition Algorithms for Automatic Detection of Compressional- and Shear-Wave Arrival Times Bull. Seismol. Soc. Am. 105 1628-1641, doi:10.1785/0120140267.


    • If you are interested in obtaining the source code, please contact Petros Bogiatzis.



    References

    1. Grossman, A., Morlet, J., 1984. Decomposition of Hardy functions into square integrable wavelets of constant shape, SIAM J. Math. Anal., 15, 723-736
    2. Hudgins, L., Friebe, C., and Mayer, M., 1993. Wavelet Transforms and Atmospheric Turbulence, Phys. Rev. Lett., 71, 3279–3282.
    3. Long M.D., van der Hilst R.D., 2005. Upper mantle anisotropy beneath Japan from shear wave splitting, Phys. Earth Planet. Inter., 151, 206-222.
    4. Maraun, D., Kurths, J., 2004. Cross wavelet analysis: significance testing and pitfalls, Nonlin. Processes Geophys., 11, 505-514.
    5. Okada, Y., Kasahara, K., Hori, S., Obara, K., Sekiguchi, S., Fujiwara, H., Yamamoto, A., 2004. Recent progress of seismic observation networks in Japan Hi-net, F-net, K-NET and KiK-net, Earth, Planets 6. and Space, 56, xv-xviii.
    6. Obara, K., Kasahara, K., Hori, S., Okada, Y., 2005. A densely distributed high-sensitivity seismograph network in Japan: Hi-net by National Research Institute for Earth Science and Disaster Prevention, Review ofScientific Instruments, 76, 021301.

    Back to other research projects.

    Department of Earth and Planetary Sciences / Harvard University / 20 Oxford Street / Cambridge / MA 02138 / U.S.A. / Telephone: +1 617 495 2350 / Fax: +1 617 496 1907 / Email: reilly@eps.hartvard.edu