Earthquake body wave extraction using sparsity-promoting polarization filtering in the time-frequency domain

TL;DR

This paper advances seismic signal processing by developing an enhanced sparsity-promoting polarization filtering method in the time-frequency domain to isolate earthquake body waves from surface wave interference.

Contribution

It introduces a tailored filter set based on predicted seismic ray angles, improving automated extraction of body waves in complex seismic data.

Findings

Effective separation of body and surface waves demonstrated on synthetic data.

Successful application to real earthquake data from Guerrero, Mexico.

Enhanced polarization filtering improves phase discrimination in seismic analysis.

Abstract

Seismic waves generated by earthquakes consist of multiple phases that carry critical information about Earth's internal structure as they propagate through heterogeneous media. These phases provide constraints from different regions of the Earth, such as the crust, mantle, and even the cores. The choice of phase depends on the study target and scientific objective: surface waves are suited for imaging shallow structures, whereas body waves yield higher-resolution information at depth. A key challenge in body-wave studies is that the low-amplitude P and S arrivals are often masked by surface waves overlapping in both time and frequency. Although body waves typically contain higher-frequency content, their spectral overlap with surface waves limits the effectiveness of conventional filtering approaches. Addressing this issue requires advanced signal-processing techniques. One such method, Sparsity-Promoting Time-Frequency Filtering (SP-TFF, Mohammadigheymasi et al., 2022), exploits high-resolution polarization information in the time-frequency domain. SP-TFF integrates amplitude, directivity, and rectilinearity to enhance phase discrimination. Here, we further develop SP-TFF by designing a filter set tailored to isolate body-wave arrivals otherwise masked by high-amplitude surface waves. The directivity filters are constructed from predicted seismic ray incidence angles, enabling focused extraction of body-wave energy and suppression of interfering phases. We evaluate the method on both synthetic tests and waveform data from the Mw 7.0 Guerrero, Mexico, earthquake of September 8, 2021, recorded by the United States National Seismic Network (USNSN). Our results show that SP-TFF provides a robust computational framework for automated body-wave extraction, integrating polarization-informed filtering into seismological data-processing pipelines.

Equations

$$\mathbf{X}(t)=\left[\begin{array}{c}{x}_1(t)\\ x_2(t)\\ x_3(t)\end{array}\right],t\in \mathbb{R},$$

$$\mathbf{X}(t,f)=\mathcal{T}\mathcal{F}\left\{\mathbf{X}(t)\right\}=\left[\begin{array}{c}{X}_1(t,f)\\ X_2(t,f)\\ X_3(t,f)\end{array}\right],(t,f)\in \mathbb{R}\times \mathbb{R},$$

$$X_i^{\text{STFT}}(t,f)={\int }_{-\infty }^{\infty }{x}_i(\tau )w(\tau -t)e^{-ȷ2\pi f\tau }d\tau ,$$

$${\mathbf{X}}^{\text{STFT}}(t,f)=\left[\begin{array}{c}{X}_1^{\text{STFT}}(t,f)\\ X_2^{\text{STFT}}(t,f)\\ X_3^{\text{STFT}}(t,f)\end{array}\right],$$

$$w(\tau -t)=\frac{1}{\sigma \sqrt{2\pi }}\exp \left[-\frac{(\tau -t{)}^2}{2{\sigma }^2}\right]$$

$$\big{(}\mathbf{C}(t,f)-\lambda_{i}(t,f)\,\mathbf{I}\big{)}\,\mathbf{u}_{i}(t,f)=0.$$

$$\mathbf{C}(t,f)=\big{[}C_{ij}(t,f)\big{]}_{i,j=1}^{3}\in\mathbb{C}^{3\times 3}.$$

$$C_{ij}(t,f)=X_i(t,f)X_j^{\ast }(t,f),$$

$$\mathbf{x}=\mathbf{G}\mathbf{\alpha },\mathbf{G}\in {\mathbb{C}}^{L\times L^2},\mathbf{\alpha }\in {\mathbb{C}}^{L^2\times 1},$$

$$\mathbf{\alpha }=\arg \underset{\mathbf{\alpha }}{\min }\frac{1}{2}\Vert \mathbf{G}\mathbf{\alpha }-\mathbf{x}{\Vert }_2^2+\mu \Vert \mathbf{\alpha }{\Vert }_1,$$

$$\arg \underset{\mathbf{\alpha }}{\min }\Vert \mathbf{\alpha }{\Vert }_1,\text{subject to}\frac{1}{2}\Vert \mathbf{G}\mathbf{\alpha }-\mathbf{x}{\Vert }_2^2\le ϵ.$$

$$Re(t,f)=\frac{{\lambda }_1(t,f)-{\lambda }_2(t,f)-{\lambda }_3(t,f)}{{\lambda }_1(t,f)}$$

$${\Psi }_{Re}(Re(t,f))=\{\begin{array}{ll}1,& -1<Re(t,f)<\alpha ,\\ \cos \left(\frac{\pi (Re(t,f)-\alpha )}{2(\beta -\alpha )}\right),& \alpha <Re(t,f)<\beta ,\\ 0,& \beta <Re(t,f)<1.\end{array}$$

$$D_i(t,f)$$

$$\displaystyle\Psi_{D}\big{(}D_{i}(t,f)\big{)}$$

$$i\in \{T,R,Z\}.$$

$$A(t,f)=\sqrt{2{\lambda }_1(t,f)}$$

$$\displaystyle\Psi_{A}\big{(}A(t,f)\big{)}$$

$$\boldsymbol{\Phi}_{R}(t,f)=1-\Big{\{}1-\boldsymbol{\Phi}_{Re}(t,f)\Big{\}}\circ\Big{\{}1-\boldsymbol{\Phi}_{D}(t,f)\Big{\}}\circ\Big{\{}1-\boldsymbol{\Phi}_{A}(t,f)\Big{\}},$$

$${\mathbf{\Phi }}_E(t,f)=1-{\mathbf{\Phi }}_R(t,f).$$

$$SM(k,l)=\sqrt{2{\lambda }_1(k,l)}$$

$$Sm(k,l)=\sqrt{2{\lambda }_2(k,l)}$$

$$D(k,l)=\sqrt{D_R(k,l{)}^2+D_Z(k,l{)}^2}$$

Full text

1 Introduction

Seismic waves are elastic disturbances generated by earthquakes, volcanic activity, or artificial sources, propagating through the interior of the Earth as vibrations of mechanical particles [1]. The propagation of the seismic wavefield is governed by the elastodynamic wave equation, with the original source function modified by interactions with the Earth’s heterogeneous three-dimensional structure. These interactions generate different seismic phases and impart additional characteristics such as transmission, reflection, diffraction, dispersion, and attenuation. Ultimately, the full wavefield carries critical information about the elastic and structural properties of the subsurface, which can be retrieved once it is recorded on the Earth’s surface by seismometers, highly sensitive instruments designed to measure ground motion [2].

In practice, this physical foundation has motivated the development of a wide range of seismological methods that often focus exclusively on specific seismic phases, relying on first-order estimates of wave propagation patterns and arrival times [3]. This has further underscored the importance of developing robust techniques for extracting and separating individual seismic phases from the fully recorded wavefield, including applications to the separation of the upward and downward wavefield [4], shear-waves extraction passed from the lithosphere [5], the extraction of surface wave phases from ambient seismic noise for time-lapse monitoring [6], and body and coda wave discrimination [7, 8], In particular, the extraction of high-frequency body waves has become an important objective in seismology due to their strong sensitivity to fine-scale velocity perturbations in relatively deeper structures [9]. Extracting such phases, including direct, refracted, and reflected body waves, from seismograms is challenging, as their signals are typically masked by the dominant energy of surface waves, which decay more slowly with distance and therefore dominate body wave amplitudes [10]. Despite these difficulties, advances in phase extraction techniques have demonstrated that body wave phases can be separated from the continuous seismic wavefield [11].

[12] proposed a processing scheme to extract Moho-reflected PpPp phases emerging in the teleseismic P coda using a deconvolution approach, in which the source wavelet of teleseismic P waves was estimated from vertical component records of a seismic array through non-linear waveform analysis. Building on passive approaches, [13] applied seismic interferometry to ground motions from clusters of regional earthquakes, retrieving direct and reflected plane waves that provide structural information. Subsequently, [14] demonstrated the retrieval of diving P and refracted waves from ambient noise records of a dense urban network in Long Beach, California. Similarly, [11] tackled the extraction of PcP—a weak phase often hidden in the P coda—by employing a data-independent strategy based on the slant-stacklet transform and applying it to USArray data. Extending interferometric concepts, [15] constructed stacked autocorrelograms of ambient seismic noise to image body wave reflections, while [16] reported clear observations of triplicated PKP branches (df, bc, ab) from stacked empirical Green’s functions obtained by correlating noise records across dense seismic arrays in South America and China. Brenguier et al. (2016) and Nakata et al. (2016) further demonstrated the temporal stability of direct virtual body waves between arrays at Piton de la Fournaise volcano, paving the way for continuous, passive monitoring of ballistic phases [17]. More recently, [18] retrieved the P-wave reflectivity response of a thick sedimentary basin beneath Jakarta, Indonesia, from phase-weighted stacks of ambient-noise autocorrelations recorded on vertical components.

[19] developed and applied a frequency–wavenumber (f–k) filter to enhance data quality and detect small phases, such as PP precursors, in vespagrams. [9] demonstrated that vehicle traffic, one of the most energetic and permanent anthropogenic seismic sources, can be exploited to reconstruct high-frequency body waves propagating across the San Jacinto Fault in Southern California to depths of several kilometers. Extending this concept, [20] identified highway and railroad traffic as primary sources of high-frequency body-wave energy and showed that selective stacking of cross-correlation functions during vehicle and train passages yields clear body-wave arrivals. Similarly, [21] applied ambient-noise cross-correlation to dense seismic arrays in northeast China and successfully retrieved body-wave signals. Working with data from a dense array in the Dehdasht region of southwestern Iran, [22] simultaneously recovered high-frequency body waves and surface waves from the correlated noise field. To address synthetic passive signals, [23] proposed and tested a Radon correlation approach for body-wave extraction. In parallel, [24] evaluated the performance of a publicly accessible CNN-based P-phase picker across multiple seismic networks of varying scale and instrumentation deployed to monitor long wall coal mining activity.

Improving the resolution of body-wave phase extraction remains an important and active area of research in seismology, requiring the development of increasingly sophisticated methods. In this study, we implement a high-resolution time–frequency filter for extracting earthquake body-wave phases using the Sparsity-Promoting Time–Frequency Filtering (SP-TFF) technique recently introduced by [7]. This approach builds upon the Sparsity-Promoting Time–Frequency Representation (SP-TFR) framework, providing a refined decomposition of seismic wavefields in the time–frequency domain based on particle-motion characteristics, thereby enabling the separation of distinct seismic phases. This decomposition further allows the application of rectilinearity, directivity, and amplitude attributes in the time–frequency domain for either extracting or suppressing different seismic phases. Focusing on the propagation pattern and ray path of body waves, we employ ray tracing to estimate the incidence angle of the desired phases and to design directivity filters in the time–frequency domain. These filters capture the energy along the target direction and maximize the extraction of the wavefield associated with those phases.

Our methodology builds on ideas first introduced at the EGU General Assembly [25]. We briefly review the methodology and provide supplementary analysis of the underlying assumptions. We then schematically illustrate the directivity filter design. Finally, we present an optimal iterative strategy for solving the SP-TFR, based on the Fast Iterative Shrinkage–Thresholding Algorithm (FISTA). We organize the article as follows: first, we review the sparsity-promoting time-frequency filtering (SP-TFF) methodology in Sec. (1.1), then, we discuss setting up filer parameters for body waves extraction in Sec. (2). Afterward, in Sec. (2) numerical examples with high-quality synthetic data and real data are presented to show the proposed method’s performance on body wave extraction. Finally, we provide a discussion on the results and conclude the paper.

1.1 Polarization analysis in the TF-domain

Noteworthy, we use a continuous notation in this section. The Time-Frequency (TF)-domain polarization parameters of particle motion, recorded by three-component seismic data, are defined for a continuous three-component signal

[TABLE]

where $x_{i}(t)$ denotes the $i$th component of ground motion.

The TF coefficients of the three components are obtained via a continuous TF operator $\mathcal{TF}\{\cdot\}$:

[TABLE]

where $X_{i}(t,f)$ is the TF representation of $x_{i}(t)$. Several well-established methods exist for computing (2), including the continuous wavelet transform [26], the Wigner–Ville distribution [27], the Short-Time Fourier Transform ( short time fourier transform (STFT)) [28], and the Stockwell Transform ( Stockwell transform (ST)) [29].

Taking the STFT as an example, the continuous TF representation of $x_{i}(t)$ is given by

[TABLE]

where $w(\cdot)$ is an analysis window, $f$ is frequency (Hz), and $\jmath^{2}=-1$.

Accordingly, the three-component TF representation can be expressed as

[TABLE]

where

[TABLE]

represents a Gaussian window with standard deviation $\sigma$, centered at time $t$. In this continuous formulation, $t\in\mathbb{R}$ denotes time, and $f\in\mathbb{R}$ frequency.

Returning to our main topic, by integrating the TF coefficients of the three components specified in (2), we can calculate the TF-dependent polarization characteristics of 3D motion. This is done in terms of the eigenvectors $\mathbf{u}_{1}(t,f)$, $\mathbf{u}_{2}(t,f)$, $\mathbf{u}_{3}(t,f)$ and eigenvalues $\lambda_{1}(t,f)$, $\lambda_{2}(t,f)$, $\lambda_{3}(t,f)$, which are obtained by solving

[TABLE]

The cross-spectral (covariance) matrix $\mathbf{C}(t,f)$ is defined at each time-frequency location:

[TABLE]

Here, $i,j=1,2,3$ index the different component pairs of the signal. The individual TF-domain cross-spectral terms are computed as

[TABLE]

where ∗ denotes the complex conjugate. In practice, suitable normalization factors are applied depending on the chosen TF decomposition (e.g., STFT, wavelet, or Stockwell transform).

In the computation of (6) using (7), a weakly stationary assumption is often made for the three components of the time series, implying $\mu_{i}\approx 0$ for $i=1,2,3$. This assumption essentially reduces the definition of the covariance matrix elements to auto- and cross-correlations. While this simplifies calculations, the resulting $\mathbf{C}(t,f)$ may not fully capture the variability and dynamics of the original signal, particularly at very low and very high frequencies. To illustrate the effect of this approximation, we perform a numerical simulation of a monochromatic wave of frequency $f$, propagating in 3D space and recorded on two components of a sensor in the $X$ and $Y$ directions. The simulation covers a frequency range from 0 to the Nyquist frequency (100 Hz). The exact cross-covariance of the two components versus the detrended (demeaned) cross-correlation is shown in Fig. 1.

As shown, the assumption provides a reasonable estimate of the true cross-covariance value, except in the very low-frequency range (<1 Hz) and at very high frequencies near the Nyquist frequency. In practice, seismic data naturally have a limited frequency band, and under this condition the assumption appears satisfactory in terms of accuracy.

1.2 Enhancing the resolution by Sparsity-promoting TF-decomposition

The STFT representation (3), and similarly the ST transform, can be expressed in matrix form as a linear system:

[TABLE]

where $\boldsymbol{x}$ is the vector of observed time-domain samples, $\boldsymbol{G}$ is the forward operator, and $\boldsymbol{\alpha}$ is a vectorized form of the TF coefficients defined in (3). Further details about the construction of $\boldsymbol{G}$ can be found in [30, 31]. Since (9) is an underdetermined system, infinitely many TF maps can represent the signal. The optimal map can be obtained through the inclusion of a priori information within a regularization framework [30, 31].

A sparsity-promoting regularization selects a TF model with the fewest non-zero coefficients by solving

[TABLE]

where the $\ell_{2}$- and $\ell_{1}$-norms represent the fit-to-data and regularization terms, respectively, and $\mu>0$ controls the resolution of the TF map. When properly tuned, $\mu$ allows discrimination between closely spaced events in time and frequency while ensuring accurate reconstruction [30, 31]. Importantly, the use of the $\ell_{1}$ norm retains convexity of the cost function, enabling efficient optimization and avoiding local minima associated with non-convex penalties.

This problem is equivalently written in constrained form, following the Karush–Kuhn–Tucker (KKT) conditions [32]:

[TABLE]

Here, $\epsilon$ controls the tolerance on the reconstruction error, preventing overfitting to Gaussian noise and improving robustness. Several algorithms can solve (10), including Iteratively Reweighted Least Squares (IRLS) [33, 34], Split-Bregman iterations [31, 35], and the Fast Iterative Shrinkage-Thresholding Algorithm ( fast iterative soft thresholding algorithm (FISTA)) [36]. In this study, we adopt the FISTA approach, which provided the most efficient solution. A pseudocode implementation is shown in Algorithm (1).

Subsequently, the derived sparsity-promoting time-frequency representation (SP-TFR) from (10) was applied to formulate an adaptive filter capable of extracting different seismic wave phases. A brief overview of the adaptive filtering procedure in the TF domain is presented in Subsection. 1.3.

Obtaining a high-resolution representation in the time-frequency domain has always been a big challenge. On this subject, combining SP-TFR with formulation of the EDPA responds favorably to this need of the scientific community. In the last decade, sparsity has been proved as a promising tool to solve inverse problems for a high-resolution solution which mathematically may be nonunique. It has successfully been applied to many areas of data processing and inversion, including deconvolution/deblurring, migration, tomography, interpolation, and Radon transform, just to name a few [31]. Nevertheless, TFR is considered as an inverse problem in a way that the TF coefficients defined as a solution of a linear equation. Since the desired linear system is considered under-determined, there are infinitely many TF maps to represent the signal. An SP regularization enables selecting a TF model with minimum nonzero coefficients by solving a constrained optimization problem.

1.3 Filter components and filter design

1. Rectilinearity Attribute: Rectilinearity is one of the three parameters used to distinguish between linear and elliptical particle motion. The degree of rectilinearity is

[TABLE]

defined as a rectilinearity measure to discriminate between the rectilinear motion of Love and body waves and the elliptical motion of Rayleigh waves. Hence, the rectilinearity filter is designed in the TF-domain as

[TABLE]

2. Directivity Attribute: Another parameter to discriminate between different seismic phases is the directivity attribute. This parameter is based on the direction of particle motion. A directivity measure is defined as the absolute value of the dot product of the first eigenvector with the base vectors

[TABLE]

Then normalizing the measure in the TF plane. Accordingly, the directivity filter is designed in the TF-domain as

[TABLE]

3. Amplitude Attribute: Given that the amplitudes of body waves are much smaller than those of surface waves, separating these waves from surface waves is a challenge worldwide. Thus, the amplitude parameter can be recognized as the most critical parameter to extract body waves from surface waves, specifically at low frequencies. The amplitude attribute is defined as

[TABLE]

and the amplitude filter is designed in the TF-domain as

[TABLE]

Integrating all these filters the total time–frequency (TF) reject filter for suppressing a seismic phase is constructed by combining the rectilinearity, directivity, and amplitude filters as

[TABLE]

where $\circ$ denotes the Hadamard (element-wise) product.

Similarly, an extract filter for isolating a specific seismic phase can be defined as

[TABLE]

The filtering process is performed by element-wise multiplication of the chosen filter $\boldsymbol{\Phi}(t,f)$ with the SP-TFR of the three components. The filtered signal in the time domain is then reconstructed via the inverse transform formulation in (9), yielding the sparsity-promoting time–frequency filter.

2 Results and Discussion

Body wave extrcation from earthquake waveforms

In this section, we demonstrate the performance of the proposed method in extracting body-wave phases from earthquake waveforms, using both synthetic and real data. For the synthetic tests, we model an earthquake of moment magnitude $M_{w}=7.0$, occurred at a depth of 21.8 km in Guerrero, Mexico, on Wednesday, September 8, 2021, at 01:47:43 UTC, associated with a reverse-thrust fault. The analysis focuses on recordings from the TZNT and VBMS stations of the United States National Seismic Network (USNSN), evaluating the extraction of body waves from overlapping Rayleigh waves on the radial and vertical components (Fig. (3)). For the real data case, we analyze recordings from the OGNE station of the same network.

Synthetic examples

To generate synthetic data, we integrated the seismic network geometry and hypocenteral parameters and focal mechanism information of the 2021 Guerrero earthquake. The 3-component synthetic waveforms of particle velocity are enerated using using the AxiSEM library through the IRIS Synthetics Engine (Syngine) within the ObsPy software [37]. An spectral-element method generated wave propagation in a three-dimensional, non-elastic, anisotropic spherical model in the f135ak D-1 Earth model. We synthesize recordings for both station VBMS to evaluate the effectiveness of the SP-TFF method in separating body waves from Rayleigh waves. The figure below shows the locations of the earthquake and station VBMS.

The three-dimensional synthetic earthquake data were generated The simulation was performed using the AxiSEM library through the IRIS Synthetics Engine (Syngine) within the ObsPy software. The generated data were processed as follows: first, detrending was applied; then, time decimation by a factor of 8 was performed to obtain a dataset with a sampling interval of 2 seconds. Next, the seismograms were rotated into the transverse–radial–vertical coordinate system. To make the simulation more realistic, Gaussian noise $n\in\mathbb{R}^{L\times 1}$ was added (bandpass filtered in the range 0.02–0.5 Hz) to achieve a signal-to-noise ratio (SNR) of 10. The waveforms recorded on the transverse, radial, and vertical components for the synthetic earthquake at station VBMS are shown in panels (a)–(c) of the figure below.

To ensure the efficiency of the SP-TFF method, its results are compared with those obtained using the Pingar method. In this approach, the Stockwell Transform (ST) is used to compute the time–frequency representation (TFR), and the time–frequency domain polarization parameters are obtained by fitting particle motion to a parametric ellipse using the combined TFR of the three components. More precisely, the set of polarization parameters, including $a(k,l)$ (length of the SM axis of the parametric ellipse), $b(k,l)$ (length of the Sm axis of the parametric ellipse), $I(k,l)$ (tilt angle of the ellipse relative to the horizontal plane at $(k,l)$), $\Omega(k,l)$ (azimuth of the major axis), and $\phi(k,l)$ (phase measured at the time of maximum displacement), are determined. Here, $l=0,\dots,n$ and $k=0,1,\dots,L$ denote the time and non-negative frequency indices, respectively.

The time–frequency representations of the radial, transverse, and vertical components of the synthetic earthquake data obtained by the ST method can be seen in panels (a), (c), and (e) of the figure below. The corresponding time–frequency representations for the SP-TFF method are shown in panels (b), (d), and (f).

It is clearly observable that the results obtained using the SP-TFF method are significantly more accurate than those obtained with the ST method proposed by Pingar (2006). It is also noted that the maximum amplitude in SP-TFF is larger than that in ST, while the energy dispersion in ST is considerably higher than in SP-TFF. As seen on the left side of the figure, with the ST method, waves are inseparably overlapping and cannot be distinguished. In contrast, the SP-TFF method (right side) separates seismic waves with high clarity and precision, allowing the distinct identification of seismic phases. The overall patterns of surface and body waves are discernible in both TFRs.

The major and minor axes of the particle-motion ellipses are denoted as $SM$ and $Sm$, respectively. In panels (a) and (c), the major and minor axes obtained using Pingar’s method (2006) are shown, whereas panels (b) and (d) display the corresponding axes derived using Eigenvalue Decomposition Polarization Analysis (EDPA) applied to SP-TFF.

[TABLE]

[TABLE]

As previously defined, $k$ and $l$ correspond to the frequency and time indices, respectively.

Next, by combining the time–frequency domain polarization parameters obtained from SP-TFF, the data are processed to filter out the Rayleigh wave phases. The polarization parameters in the time–frequency domain obtained from SP-TFF are used to define an adaptive filter for separating S and P waves as well as removing Love and Rayleigh waves.

To filter the Rayleigh wave phases, the directivity measure relative to the radial–vertical plane is calculated as:

[TABLE]

Additionally, a rectilinearity filter is designed and applied to remove these waves. Accordingly, the SP-TFF components of the Rayleigh-wave-filtered data are shown in panels (b), (d), and (f) of Figure 5-4. As clearly observable, the Rayleigh waves on the radial and vertical components are removed with high accuracy, while the scattered energy of the body waves in the SP-TFF-filtered data remains on the radial and vertical components.

The transverse, radial, and vertical components of the Rayleigh-wave-filtered data, reconstructed in the time domain, are shown in panels (b), (d), and (f) of the figure below, respectively.

As shown, the SP-TFF method effectively filters the Rayleigh waves without affecting the body waves on the radial and vertical components, and has only a minor effect on other phases on the transverse component. As a final assessment, by applying Equation (3-33) to the SP-TFR, the Rayleigh-wave phases are extracted on all three components. Panels (d) and (f) of Figure 7-4 show the radial and vertical seismograms corresponding to the separated Rayleigh waves. It is noteworthy that the Rayleigh phases have been extracted from the complete seismograms with high accuracy.

Another significant challenge is the identification and separation of Love and SH waves, both of which are recorded on the transverse component and often overlap. Considering that body waves have a broader frequency range than surface waves, SH waves typically occupy a wider bandwidth than Love waves. This characteristic can be very helpful for distinguishing SH waves from Love waves in the time–frequency representation. To demonstrate the effectiveness of the SP-TFF method in addressing this challenge, we analyze the data recorded at station TZTN. The figure below shows the locations of the earthquake and this station.

The waveforms recorded on the transverse, radial, and vertical components for the arrival of the synthetic earthquake at station TZTN are shown in panels (a)–(c) of the figure below.

The time–frequency representations of the radial, transverse, and vertical components of the synthetic earthquake data at station TZTN obtained using the ST method are shown in panels (a), (c), and (e) of the figure below. The corresponding time–frequency representations obtained using the SP-TFF method are shown in panels (b), (d), and (f).

In Figure 10-4, panel (b), the SH wave can be observed around 480 seconds, overlapping with the Love wave; we will address the separation of these two waves in the following section.

The figure below shows the time–frequency representations of the major and minor axes of the particle-motion ellipses. Panels (a) and (c) correspond to the $SM$ and $Sm$ axes obtained using Pingar’s method (2006), while panels (b) and (d) show the $SM$ and $Sm$ axes derived from applying EDPA to the SP-TFF method.

In this section, by combining the time–frequency domain polarization parameters obtained, the data are processed to filter out the Love-wave phases. For Love-wave filtering, a directivity filter is designed, which defines the directivity measure relative to the transverse axis based on a combination of adaptive parameters. An amplitude filter is also constructed to remove the Love waves according to Equation (3-32). By visual inspection, the approximate arrival time of the Love wave can be determined. This procedure is performed to limit the filtering region and affects both the high-amplitude SH waves and the Love waves.

The results of applying the Love-wave removal filter to the SP-TFF for the transverse, radial, and vertical components are shown in panels (a), (c), and (e) of Figure 12-4. As observed, the energy corresponding to the Love wave in the time–frequency plane is significantly removed, leaving only the scattered energy related to the body waves. The SP-TFFs of the radial (see Figure 12-4(c)) and vertical (see Figure 12-4(e)) components remain unaffected by the filter.

The black waveforms in panels (a), (c), and (e) of Figure 13-4 represent the transverse, radial, and vertical components reconstructed in the time domain after filtering. The Love wave is almost completely removed from the transverse component in the time domain, while the SH waves remain in the seismograms. It is observed that the Love-wave filtering process has no effect on other phases in the radial and vertical components. A promising result of time–frequency filtering using the SP-TFF method in this example is that, as shown in panel (a) of the figure below, an SH phase around 480 seconds, which was previously masked by high-amplitude Love waves, is successfully recovered after filtering.

As a final assessment, by applying Equation (3-33) to the SP-TFR, the Love-wave phases are extracted on all three components. In Figure 14-4(a), the separated Love wave is clearly visible. Notably, the Love-wave phase has been extracted from the complete seismograms with high accuracy.

As shown, the SP-TFF method effectively filters the Love wave without affecting the body waves on the radial and vertical components, and has only a minor effect on other phases on the transverse component.

Real data examples

To validate the SP-TFF method, we use the real seismograms of the Guerrero 2021 earthquake recorded at OGNE station. The waveforms were first processed by applying detrending, time-domain resampling at 2-second intervals, deconvolution of the instrument response, and conversion to velocity. The figure below shows the locations of the earthquake and station OGNE, which is analyzed in this study. The processed seismograms of station OGNE are shown in Figure 16-4. These seismograms are represented in the transverse–radial–vertical coordinate system. The time-frequency representations of these components using both Pingar’s method and SP-TFR are shown in the figure below. Panels (a), (c), and (e) correspond to the transverse, radial, and vertical components obtained using Pingar’s method (2006), while panels (b), (d), and (f) show the time–frequency representations of the same components obtained using the SP-TFR method. As observed for the real data, the time–frequency representation obtained using SP-TFR is significantly more accurate than that obtained using Pingar’s method (2006).

The results are similar to those observed in the synthetic dataset, where SP-TFR (right column) provides a more accurate time–frequency representation compared to the ST method implemented by Pingar (2006) (left column). Although the surface and body waves in the real data do not exhibit a clear pattern as in the synthetic data, they can still be partially distinguished. The high-amplitude regions in panels (b), (d), and (f), whose frequency increases over time, correspond to surface waves, while the body waves exhibit a broader frequency band. In each of these figures, examples of surface and body waves have been delineated using blue and green lines, respectively. To better visualize the separation, we focused on body waves that overlap with surface waves.

The distinct polarization patterns between Love and Rayleigh waves are more clearly observed in the $SM$ and $Sm$ axes in the time–frequency domain shown above. The elliptical particle motion of Rayleigh waves can be distinguished from the linear motion of Love and body waves in the $SM$ and $Sm$ axes based on differences in their amplitudes. Filtering Body Waves Using SP-TFF

The filtering process at this stage is carried out in the same manner as described in the synthetic dataset analysis section. The results are presented below.

The results of applying the Love-wave removal filter to the SP-TFR for the transverse, radial, and vertical components are shown in panels (a), (c), and (e) of Figure 19-4. The SP-TFRs of the radial (see Figure 19-4(c)) and vertical (see Figure 19-4(e)) components remain unaffected by the filter.

The transverse component in the figure above corresponds to the separation of body waves from the Love wave. As can be seen by comparing panels (a) and (b), the high-amplitude Love waves around 500 seconds visible in (b) are effectively removed in (a), leaving only the scattered energy of the body waves. To observe the effect of the filter on Rayleigh-wave removal, the radial and vertical components are also examined. Panels (d) and (f) show that, compared to (c) and (e), the Rayleigh waves are well removed while the body waves remain.

The black waveforms in panels (a), (c), and (e) of the figure below display the transverse, radial, and vertical components of the Love-wave-filtered data in the time domain. These results confirm that the method used in this study specifically removes the Love wave from the transverse component in the time domain without significantly affecting other phases.

Finally, the body-wave phases of the earthquake for the real dataset are separated from the surface waves by applying Equation (3-33) to the SP-TFR of all three components. The separated Love wave and body waves on the transverse component are shown in panel (a) above. Similarly, Figure 21-4(b) and (c) display the radial and vertical components of the Rayleigh and body waves, which have been successfully separated.

3 Conclusion and future works

This study demonstrates that Sparsity-Promoting Time–Frequency Filtering (SP-TFF) provides a powerful framework for the reliable extraction of earthquake body-wave phases that are otherwise obscured by surface-wave energy. By incorporating polarization-based attributes—amplitude, directivity, and rectilinearity—into a tailored filtering strategy, the method achieves focused isolation of body-wave arrivals and suppression of interfering phases. Applications to both synthetic and observational data from the 2021 $M_{w}7.0$ Guerrero earthquake confirm that SP-TFF is effective, scalable, and well-suited for automated integration into modern seismological workflows. Beyond improving the resolution of phase extraction, this approach offers a flexible computational tool that can enhance earthquake characterization, support large-scale waveform analyses, and ultimately advance Earth structure imaging. Open access to the codes ensures transparency, reproducibility, and broader adoption within the seismological community. The next step in our study will focus on extracting body waves from Empirical Green’s Functions derived from cross-correlations of ambient seismic noise at local scales.

WWSSN World Wide Standard Seismograph Network GSD Geological Survey Department AR accuracy rate WAC West African Craton MAR Mid-Atlantic Ridge SP-TFF Sparsity Promoting Time-Frequency Filtering KKT Karush–Kuhn–Tucker IRLS Iterated Reweighted Least Squares PA polarization analysis HR high-resolution ROSES Remote Online Sessions for Emerging Seismologists TFR time-frequency representation SP-TFR sparsity-promoting time-frequency representation SP-TFF sparsity-promoting time-frequency filtering DOP degree of polarization FISTA fast iterative soft thresholding algorithm STFT short time fourier transform ST Stockwell transform WT Wavelet transform TSR time-scale representation SST synchrosqueezing transform SM Semi-major Sm Semi-minor EDPA eigenvalue decomposition polarization analysis SNR signal to noise ratio SEM Spectral Element Method TF Time-Frequency AFZ Akwapim Fault Zone CBF Coastal Boundary Fault RFZ Romanche Fractured Zone