Sparsity-promoting approach to polarization analysis of seismic signals in the time-frequency domain

In this paper, we present a new approach to the TF-domain PA methods. More precisely, we provide an in-detailed discussion on rearranging the eigenvalue decomposition polarization analysis (EDPA) formalism in the frequency domain to obtain the frequency-dependent polarization properties from the Fourier coeﬃcients owing to the Fourier space orthogonality. Then, by extending the formulation to the TF-domain and incorporating sparsity-promoting time-frequency representation (SP-TFR), we alleviate the limited resolution when estimating the TFdomain polarization parameters. The ﬁnal details of the technique are to apply an adaptive sparsity-promoting time-frequency ﬁltering (SP-TFF) to extract and ﬁlter diﬀerent phases of the seismic wave. By processing earthquake waveforms, we show that by combining amplitude, directivity, and rectilinearity attributes on the sparse TF-domain polarization map of the signal, we are able to extract or ﬁlter diﬀerent phases of seismic waves.


Sparsity-promoting approach to polarization analysis of seismic signals in the time-frequency domain
Hamzeh Mohammadigheymasi, Paul Crocker, Maryam Fathi, Eduardo Almeida, Graça Silveira, Ali Gholami, and Martin Schimmel Abstract-Time-frequency (TF)-domain polarization analysis (PA) methods are widely used as a processing tool to decompose multi-component seismic signals.However, as a drawback, they are unable to obtain sufficient resolution to discriminate between overlapping seismic phases, as they generally rely on a low-resolution time-frequency representation (TFR) method.In this paper, we present a new approach to the TF-domain PA methods.More precisely, we provide an in-detailed discussion on rearranging the eigenvalue decomposition polarization analysis (EDPA) formalism in the frequency domain to obtain the frequency-dependent polarization properties from the Fourier coefficients owing to the Fourier space orthogonality.Then, by extending the formulation to the TF-domain and incorporating sparsity-promoting time-frequency representation (SP-TFR), we alleviate the limited resolution when estimating the TFdomain polarization parameters.The final details of the technique are to apply an adaptive sparsity-promoting time-frequency filtering (SP-TFF) to extract and filter different phases of the seismic wave.By processing earthquake waveforms, we show that by combining amplitude, directivity, and rectilinearity attributes on the sparse TF-domain polarization map of the signal, we are able to extract or filter different phases of seismic waves.The SP-TFF method is evaluated on synthetic and real data associated with the source mechanism of the Mw = 8.2 earthquake that occurred in the southsouthwest of Tres Picos, Mexico.A detailed discussion on the results of these experiments is given, approving the efficiency of the technique in separating not only the Rayleigh from the Love waves but also to discriminate them from the body and coda waves.

I. Introduction
A Seismic wavefield recorded as a seismogram is a superposition of overlapping direct, reflected, refracted, converted, and scattered body and surface waves, contaminated by various background sources and the signal generated noise [1].It is well known that surface and body waves can carry considerable information about the subsurface structure.Depending on the research scope, any of these seismic phases can be studied, while their detection and extraction require advanced processing and analysis tools.Accordingly, multicomponent processing techniques have been developed to analyze the nonlinear and time-varying processes behind the seismic sources and the propagating environment [Refer to [2] as a rigorous survey].Among these techniques, polarization analysis methods have attracted significant attention.
Generally, the polarization analysis methods can be divided into three broad categories: time, frequency, and time-frequency (TF) domain methods.As a pioneer, Flinn [3] introduced the eigenvalue decomposition polarization analysis in the time-domain.Likewise and in the realm of Hilbert transform, Vidale [4] introduced the analytic signal polarization technique.Although these methods are equipped by time-windowing to extract the non-stationary signal properties, they cannot discriminate between overlapping events with different frequencies.Comparatively, studies were carried out on achieving polarization properties of the seismic time-series in the frequency domain.Primarily developed by [5], the propagation direction of surface waves was obtained using the amplitude and phase spectrum of seismic waves.Likewise, Samson and Olson [6] proposed a technique to provoke the polarization state based on eigenvalue decomposition of the spectral matrix in different frequency narrow bands.However, these techniques are incapable of analyzing non-stationary signals in the time domain.
Due to the non-stationary nature of seismic signals with overlapping phases in time and frequency, pure time-or frequency-domain methods are often difficult to discriminate between the non-stationary seismic phases adequately.To alleviate this, Jurkevics [7] proposed filtering the signals into a series of narrow frequency bands, applying short sliding time windows, and then estimating the polarization ellipse from the covariance matrix in each window at each band.Despite providing a TF-domain insight to polarization analysis, the resolution was still limited in the frequency domain due to the restricted functionality of the bandpass filtering.Therefore, by applying different time-frequency representation (TFR) methods [8,9,10] and using various polarization estimation criteria, several TF-based polarization analysis methods were proposed to analyze the non-stationary seismic signals.These TF-domain polarization estimation tools range from eigenvalue decomposition of the covariance matrix [11,12] and complex trace analysis [4,13,14], to fitting the particle motion to a parametric polarization ellipse [15].Despite attaining a TF-domain estimation these methods are not yet able to resolve closely-spaced overlapping seismic events in both time and frequency domains.
Obtaining a high-resolution TFR has always been a challenge for the scientific community and a variety of methods have been introduced, including synchrosqueezing transform (SST) [16], sparse transforms [17], SP-TFR [18,19,20], to name but a few.Remarkably, by formulating the TFR as an inverse problem, and taking advantages of sparsity-promoting (SP) regularization as a promising tool to obtain a high-resolution solution, several methods have been developed to improve the TFR resolution [18,19,20].SP-TFR found many applications in seismology, including the denoising of microseismic data [20] and attenuation of seismic ground roll noise [19].
On the other hand, filtering or extracting different phases of seismic waves is a concerning challenge in the seismic community.Although in a laterally homogeneous structure, the Love wave appears mainly on the transverse and the Rayleigh wave on the vertical and radial components, due to the lateral heterogeneity and anisotropy, it is seldom the case in real data [4].Likewise, this assumption fails in miss-oriented horizontal sensor components, which is a global problem even for the best-installed seismic networks [21].Hence, it is always challenging to discriminate between the body and Rayleigh waves on the vertical and radial components as well as SH and Love waves on the transverse component.Accordingly, polarization filtering techniques have been developed to extract or filter a specific phase from other phases using the analyzed polarization information.As pioneers, Flinn [3], Montalbetti and Kanasewich [22], and Vidale [4] utilized time-domain rectilinearity and directivity attributes to amplify body wave phases teleseismic data.In like manner, Samson and Olson [6] applied these criteria to filter ultra low-frequency magnetic field signal fluctuations in the frequency domain.Similarly, Schimmel et al. [14] designed a TF domain filter based on the degree of polarization (DOP) extracted from the semi-major and semi-minor axis of the elliptical motion of the three-component ambient noise data to filter the elliptical particle motion.In an intuitive method, Pinnegar [15] utilized semi-major, semi-minor, inclination, and azimuth parameters to discriminate between the circular and linear polarization to filter the Rayleigh waves.Although the introduced filtering scheme is efficient in filtering the elliptically polarized phases like Rayleigh waves, it still faces challenges in filtering the linear phases because the method attains a null value of azimuth and inclination angles analyzing linear particle motions.
This article begins with an elaborative review of the eigenvalue decomposition polarization analysis EDPA; we formulate EDPA as a function of frequency.On this basis, we obtain high-resolution TF-domain polarization parameters by extending the formulation from frequency to the TF-domain and combining with the SP-TFR.Afterward, by extending the polarization filtering method of Pinnegar [15] to incorporate high-resolution TF-domain information of directivity, rectilinearity, and amplitude properties, we design suitably defined filters to accept (or reject) linear and elliptical seismic phases, including Rayleigh and Love, making it possible to separate them from body and coda waves.The main focus of the paper is to discriminate between the Love and Rayleigh from the body and coda waves.
This paper is organized in the following manner.First, in section II we elaborate the theory behind the EDPA in the time, frequency, and TF domains.We show that EDPA can be applied independently for every single frequency by taking advantage of the orthogonality of the Fourier domain.This property can be extended to the TF-domain, making it possible to obtain TF-domain polarization parameters.Then, by reviewing the SP-TFR and combining with EDPA, we obtain a high-resolution TF-domain polarization map of the signal.Afterward, by implementing TFdomain rectilinearity, directivity, and amplitude attribute and combining them with SP-TFR we introduce sparsitypromoting time-frequency filtering (SP-TFF) method to be used for filtering different phases of seismic signal.Next, in section III by conducting numerical experiments on synthetic and real earthquake data examples, we show that SP-TFF can efficiently extract and filter Rayleigh and Love waves and discriminate between linearly polarized seismic phases like Love and body and coda waves.Finally, in sections IV we discuss the results and conclude the paper.

II. Theory
To keep this paper self-contained, we rearrange the EDPA formalism in the frequency and TF domain to be combined by SP-TFR.Furthermore, to be consistent with the numerical algorithms, we exclusively present a discrete version of the mathematical equations.

A. Polarization analysis using eigenvalue decomposition
Suppose that is a seismic time-series recorded with three components aligned with the base vectors of a right-handed coordinate system {e E , e N , e Z } (east-north-vertical), {e T , e R , e Z } (transverse-radial-vertical), or {e L , e Q , e T } (tangent-normal-binormal or Frenet wave).The latter coordinate systems are obtained from the ordinary eastnorth-vertical system by rotation [see [23] for more details].Each component, x i ∈ R L×1 (i = 1, 2, 3), samples the wavefield arriving to the sensor along the time axis on t k = kδt, with δt being sampling interval and k = 0, 1, ..., 2n being the time index assuming an odd length L = 2n + 1 of seismogram.Then, the polarization properties of X can be extracted by the eigenvalue decomposition of the covariance matrix [ 22,3,7,24].The elements of the symmetric and realvalued matrix V in (2) are auto-and cross-variances of the components of the time-series defined as In (3), µ is the mean or expected value of components and is defined as [22].Assuming a weakly stationary condition for all the components of the time series, Ψ{x i } ∼ = 0, i = 1, 2, 3, (3) can be rewritten as which simplifies the definition of elements in (2) to autoand cross-correlations, or equally, where (.) T in (7) denotes the transposition operator.
The weakly stationary assumption in ( 5) is satisfied by applying DC removal or detrending.
The quadratic matrix of correlation coefficient, C, fits the particle motion ellipsoid in a least-squares sense; the parameters of this ellipsoid are obtained by solving the system of equations [7].Geometrically, solutions to (8) give directions (eigenvectors, (u 1 , u 2 , u 3 )) that the linear transformation by operator C merely elongates or shrinks; the ratio of elongation/shrinkage is given by eigenvalues, (λ 1 , λ 2 , λ 3 ).Indeed, the eigenvectors direct principal axes of the polarization motion ellipsoid, and the eigenvalues are the size of those axes; eigenvalues are sorted such that λ j ≥ λ k for j < k.

1) Eigenvalue decomposition in the frequency domain:
The first implementations of the eigenvalue decomposition on the frequency domain were proposed by [25,6]; the spectral matrix corresponding to a perturbation around the central frequency [ω − δω, ω + δω] was decomposed using eigenvalue decomposition to provoke the polarization states of the signal.This decomposition scheme is more compatible with natural signals, which are rarely composed of single-frequency polarized elements.However, by decomposing to a strictly polarized single frequency state, one can benefit from the orthogonality of the Fourier transform to extract frequency-dependent polarization properties.Here, we review and simplify the process.
The discrete Fourier domain counterpart of (1), , is obtained by modulating x i , i = 1, 2, 3, with pure sinusoids having discrete frequencies as follows Here, l = −n, ..., n, is the frequency index giving frequency content of the signal on discrete frequencies ω l = l (2n+1)δt , and  is the imaginary unit.The original signal is reconstructed by applying the inverse Fourier transform, Accordingly, the frequency-domain counterpart of ( 5) is defined as in which, (.) * , and • are complex conjugate and Hadamard or element-wise product operator, respectively, and ZL{.} is a function that picks the zero-lag element.By taking advantage of the orthogonality property of the Fourier domain, (11) can be written as Hence, the auto-and cross-correlation terms for every single frequency can be decomposed via (13).Instead of applying operation in (13) to obtain gives the frequency-dependent auto-and crosscorrelation elements.Solving system of equations ( 8) for the frequency-dependent elements gives the eigenvectors, (u 1 (l), u 2 (l), u 3 (l)) and eigenvalues, (λ 1 (l), λ 2 (l), λ 3 (l)), as a function of frequency.
2) Eigenvalue decomposition in the time-frequency domain: The process described in section (II-A1) effectively decomposes stationary signals into it's frequencydependent polarization components (eigenvectors and eigenvalues).Having to deal with the non-stationary nature of seismic data, the same definition is extended to give polarization components that depend on time and frequency by substituting a TFR in place of the ordinary Fourier transform.Much research has been devoted to finding efficient TF analysis methods and many powerful methods have been developed in the past decades, including the wavelet transform [26], the Wigner-Ville distribution [27], short time fourier transform (STFT) [8], Stockwell transform (ST) [9], etc.As an instance, the STFT representation of signal x i in ( 1) is obtained by with being a Gaussian window with standard deviation σ, centered on the time index k.The definition in ( 16) can be extended to the ST [9] by applying a time-frequency spectral localization using a window function scalable with frequency as with k and l are defined the same as (16).The parameter σ in ( 17) and ( 18) controls the resolution of the transform in the time and frequency domain; higher values of σ attains higher frequency resolution, while lower values improves the time resolution.By obtaining TF of the 3-components of the signal, the TF-domain auto-and cross-correlation terms is obtained as for time and frequency indexes , k = 0, 1, ...., 2n + 1, l = 0, 1, ...., n.Consequently, the TF-dependent eigenvectors, (u giving a TF map of polarization state of signal.This decomposition process is similar to the method introduced by Pinnegar [15].However, as we will discuss in the following sections, it can be used to filter the linearly polarized seismic phases, which is not able to be done with the Pinnegar [15] method.

B. Regularized sparsity-promoting TF decomposition
The system of equations for STFT and ST linear which allows us to define the TF coefficients as a solution of a linear system equations where α is a vectorized rearrangement of TF coefficients in (16) [see [20,19] for more details about the structure of the forward operator G].Since the linear system in ( 21) in under-determined, there exists an infinite number of TF maps for representing the signal.The desired TF map can be obtained by using some form of a priori information under the frame of regularization techniques [20,19].A sparsity-promoting regularization enables selecting a TF model with a minimum number of non-zero coefficient by solving a constrained optimization problem where α p = ( i |α(i)| p ) 1/p is the p norm of a vector x and µ > 0 is the sparsity parameter [20,19].By choosing a proper µ, one can control the resolution of the TF map and allows us to being able to discriminate between closely spaced events in time and frequency, while reconstructing the data.The optimization problem ( 22) can be solved by a variety of methods such as the split Bregman method [19] or fast iterative soft thresholding algorithm (FISTA) [28].In this study, we utilized the FISTA method to solve (22).
The obtained SP-TFR through (22) is used to design an adaptive filtering for extracting (or filtering) different phases of seismic waves.In the next section, we briefly review the adaptive filtering approach in the TF domain.

C. Adaptive filtering in the TF domain
Adaptive filtering has been extensively applied in seismology [22,4,7,14].In an intuitive scheme, Pinnegar [15] utilized a combination of inclination, azimuth, and rectilinearity attributes in the TF domain to filter the Rayleigh waves.However, his method is not able to scrutinize the pure linear polarization because of an undefined inclination angle for linear particle motions [see section (III) for more details].As a result, the method is not applicable to filter the Love wave or any other seismic phase with a linear polarity.Here, we extend his methodology by combining rectilinearity, directivity, and amplitude attributes [7,22] with the TF-domain polarization parameters obtained from SP-TFR of 3-components of the signal to introduce SP-TFF method.
1) Rectilinearity attribute: Rectilinearity is a critical parameter for discriminating between the elliptical and linear particle motion states.The purely rectilinear ground motion is modeled by one nonzero eigenvalue in the TF plane Nevertheless, due to the presence of contaminating noise, out-of-plane energy, and scattering distortions, it is seldom the case for the real data [29].To circumvent this, a degree of rectilinearity is defined as a rectilinearity measure to discriminate between the rectilinear motion of Love and body waves and elliptical motion of Rayleigh waves [7], [13], [30].Accordingly, a rectilinearity filter is designed in the TF domain as while to avoid the Gibbs phenomenon caused by abrupt frequency cut-off, the accept or reject regions are cosine tapered in (23) by incorporating suitably defined adjusting parameters α and β.The proposed filter is similar to the TF filter introduced by Pinnegar [15], whereas it can be designed to filter the linear polarization particle motion.
2) Directivity attribute: Directivity is another crucial parameters to discriminate between different seismic phases based on the direction of particle motion.More precisely, a directivity measure is defined as the absolute value of the dot product of the first eigenvector by the base vectors then normalizing the measure in the TF plane.Correspondingly, a directivity filter is designed in the TF domain as The adjusting parameters γ and λ have the role of both cosine tapering to avoid the Gibbs phenomenon and threshholding as a percentage of the maximum measure.
In the next section we present the combination of these attribute to filter seismic data.
3) Amplitude attribute: Although generally surface waves manifest themselves with higher amplitude than the body and coda waves [31], it is challenging in practice to work with amplitude attributes to discriminate them.Having a TF-domain insight to analyze the signal, enforced by SP-TFR to present it with a few sparse coefficients, accentuates the amplitude difference between the surface and coda waves.In other words, the energy of surface waves is extracted locally with a few sparse coefficients, while the body and coda waves are distributed to a broader range of coefficients due to having a broader frequency content.Hence, an amplitude attribute can be more efficiently used to discriminate surface waves from body and coda waves.Specifically, it can be employed as a tool to separate Love and SH waves, which have the same type of polarization directivity and rectilinearity.Defining an amplitude attribute as a corresponding amplitude filter is designed in the TF domain as normalized in the whole TF plane.The adjusting parameters ζ and η acts as a measure to pass (or reject) the coefficient in the TF plane, while applying the cosine tapering.In the next section we present the combination of these attribute to filter seismic data.

4) Regularized Sparsity-promoting Time Frequency Filtering :
To combine the properties of different attributes, a similar methodology to [15] can be followed.More precisely, the total TF reject filter to reject a phase is obtained by combining the rectilinearity, directivity, and amplitude filters, as Similarly, a special seismic phase can be extracted by defining an extract filter as Ψ E = 1 − Ψ R .Finally, the filtering process is applied by element-wise multiplication of Ψ with the SP-TFR of the three components; then, the filtered signal is reconstructed in the time domain by applying (21) giving the SP-TFF.
Besides all the defined criteria, subjective information can be incorporated as a constraint to control the filtering domain while rejecting or extracting seismic phases.As an example, in the case of filtering or extracting the surface waves, an approximate initial time of the Love wave can be introduced to the algorithm to limit the domain of filtering in (30).It can be picked visually on the transverse component or estimated by using the standard global dispersion curves.
In the next section, we examine the application of the proposed filtering method to filter different seismic wave phases.

III. Numerical examples
To evaluate the SP-TFF method, we test it with synthetic and real data examples by extracting and filtering the Love and Rayleigh waves.The implementation results are compared with those of the method introduced in Pinnegar [15].

A. Synthetic examples
The synthetic data corresponds to the source mechanism of the M w = 8.  mechanism and the source-receiver geometry were chosen such that the amplitude of body and coda waves is almost comparable to the surface waves making separation more challenging.
A three-dimensional synthetic seismic data was generated through the 1D ak135f earth model [32] with spectralelement method assuming 3D (an-)elastic, anisotropic and acoustic wave propagation in spherical domains.The simulation was run by using the AxiSEM library through the IRIS Synthetics Engine (Syngine) client of ObsPy software [33,34].
In the simulation, the seismic wavefield is recorded in the College Outpost, Alaska, USA [see Table .I for more details], at the azimuth of 61.56 • to the epicenter.
The generated data were preprocessed; detrended and decimated by a factor of 8 to attain a data set with a sampling rate of 2 sec.Then, the traces were rotated to the transverse-radial-vertical coordinate system.To make the simulation more realistic, the data was contaminated by a Gaussian noises, n ∈ R L×1 (bandpass-filtered in the range of [0.02, 0.5] Hz) to give a signal to noise ratio (SNR) of 10.The Transverse, Radial, and Vertical components of the total motion are shown in gray color in Fig. 1.
To evaluate the efficiency of SP-TFF, the results are compared with those obtained by Pinnegar [15].In this method, the ordinary ST is used as the TFR, and the TFdomain polarization parameters are obtained by fitting the particle motion to a parametric ellipse by incorporating the TFR of 3-components.More precisely, a set of the The TFR of the transverse, radial and vertical components of the synthetic data obtained by the ordinary ST (applied by Pinnegar [15]) are shown in the panels (a), (c), and (e) of Fig. 2; the corresponding TFR for the SP-TFR method are shown in the panels (b), (d), and (f).The SP-TFR attains a highly compact TFR with a maximum amplitude higher than the ST, while ST distributed the energy in the TF plane in a wider area.The up-chirp characteristics of surface waves are obvious in both TFRs.
The TF domain SM and Sm axes of the particle motion obtained by Pinnegar [15] method are depicted in (a) and (c) panels of Fig. 3; the corresponding TF domain SM and Sm axes for EDPA using SP-TFR are shown in in (b) and (d) panels.The SM and Sm axes for EDPA is obtained as with k and l are defined similar to (16).By considering the left panels of Figs. 2 and 3 it is evident that by applying the ordinary ST, the Rayleigh and Love waves are inseparably overlapping both in time and frequency.It can e deduced either from the TFR maps or from the SM and Sm maps.In contrast, for the SP-TFR (Right panels), the high-resolution TFR successfully separated the wavefields, giving the possibility of discriminating between different seismic wave phases.It is also seen in the SM and Sm maps in Fig. 3. Another interesting feature is the SH and Love wave pattern in the TFR.The Love wave has been concentrated around the dispersion curve; however, the SH wave has been spread in a wider frequency bandwidth around the time 800 secs.It makes using amplitude attribute more efficient in separating them.Subsequently, by incorporating the obtained TF domain polarization parameters, we process the data to filter different seismic phases.

1) Love and Rayleigh wave filtering using SP-TFF:
The TF domain polarization parameters obtained form SP-TFR, are used to define and adaptive filter according to section (II-C) to filter Love and Rayleigh waves.
To extract and filter the Love waves, we design a directivity filter by defining the directivity measure with respect to the transverse axis e T in (26), and a set of adjusting parameters γ = 0.13 and λ = 0.16 for amplitude threshholding and cosine tapering of the directivity measure.An amplitude filter by the set of parameters ζ = 0.26 and η = 0.23 was combined to define a Love-reject filter accroding to (30).The results of applying the filter on the SP-TFRs of the transverse, radial and vertical components are shown in panels (a), (c), and (d) of Fig. 4. As is shown, the energy corresponds to the Love wave in the TF plane has been significantly removed, and only scattered energy remains, which corresponds to the body and coda waves, and noise (top panel).The SP-TFRs of the radial (panel (c)) and vertical (panel (e)) components have not been affected by filtering.The black color waveform in (a), (c), and (e) panels of Fig. 1 depict reconstructed transverse, radial and vertical components after filtering in the time domain; the Love wave is almost entirely removed in the time domain, while the other phases, including the body and coda wave, and also the noise has remained in the seismogram.It is a promising feature of the SP-TFF algorithm compared to the Pinnegar [15] method, which attains a null value for the inclination and azimuth parameters corresponding to a linear particle motion.The filtering process has not affected other phases in the radial and vertical components, except a minor effect on the Rayleigh phase around the time 700s, in which the SP-TFR of the Rayleigh and Love phases fully overlaps.An interesting result of applying SP-TFF in this example is that a SH phase around 800 sec masked by high-amplitude Love waves has been recovered after filtering.To filter the Rayleigh phase, the directivity measure is computed with respect to the radial-vertical plane computed as The adjusting parameters are set to γ = 0.25 and λ = 0.3.Furthermore, a rectilinearity filter is defined by setting the parameters α = 0.1 and β = 0.12.The results of applying the filter on the SP-TFRs of 3-components are shown in the panels (b), (d), and (f) of Fig. 1.Similar to the Love wave filtering, the filtered SP-TFR only contains scattered energy of the noise, body, and coda waves in the radial and vertical components.The reconstructed filtered components are shown in the right panels shown in the right panel of Fig. 1.As shown, SP-TFF successfully filtered the Rayleigh wave without affecting the body and coda waves in the radial and vertical components and substantially affected the other phases in the transverse components.The same as for the Love wave filtering, around 800s, the Love wave has slightly been filtered.The results obtained from the SP-TFF are superior to those from Pinnegar [15] by having a very high-resolution TFR enable to separate the Rayleigh and Love waves, while in the ordinary ST the TF resolution is limited.As a final test to assess the SP-TFF method, the Love and Rayleigh phases are extracted by applying (30)    processing methods like dispersion curve inversion.In the following subsection, we examine the performance of the method on real data.

B. Real data example
The real seismogram corresponds to the real data recorded for the same earthquake and station of the synthetic model.The waveform was pre-processed including detrending, decimation, deconvolving the instrument response, and converting to velocity with a sampling rate of 2 sec.The traces were mapped to the transverse-radialvertical coordinate system.The pre-processed transverse, radial, and vertical components are shown in gray color in the top, middle, and bottom panels of Fig. 6.
The obtained TFR by applying [15] and SP-TFR methods are shown in the left and right panels of Fig. 7; similar to the synthetic example, the SP-TFR (right panels) presents a highly compact TFR comparing to the ordinary ST implemented by Pinnegar [15] (left panels).Although there is no sharp up-chirp pattern for the surface waves like in the synthetic data, there are still two separate energy panels in the SP-TFR of different components showing an increasing value of frequency by time.These two panels marked by dash-dot and continuous line ellipsis correspond to Love and Rayleigh waves, respectively.Two other panels shown by dashed ellipsis contain mostly body and coda waves.On the other hand, the TFR obtained by ST and shown in the left panel of Fig. 7 depicts a mixed and inseparable pattern of Love and Rayleigh waves.The distinct polarization pattern between the Love and Rayleigh waves is better visible in the TF domain SM and Sm axes of particle motion as shown in the right panel of Fig. 8.The elliptical particle motion of Rayleigh waves is separable from the Linear particle motion of Love, body, and coda waves in the SM and Sm axes figures.Contrarily, they have been mixed in time and frequency in the results obtained by Pinnegar [15] method, as shown in the left panel of Fig. 8.In the following subsection, we perform the adaptive filtering method to reject or extract Love and Rayleigh waves.
1) Love and Rayleigh wave filtering using SP-TFF: The TF domain polarization parameters obtained from SP-TFR are used to design an adaptive filter to filter Love and Rayleigh waves.For the real data, similar adjusting parameters of γ = 0.1, λ = 0.13, α = 0.03, and β = 0.04 to the ones defined for the synthetic example were set to design directivity and rectilinearity attributes.We only slightly changes the for amplitude filtering parameters by setting ζ = 0.16 and η = 0.19.A visual assessment chose an initial time of 630 seconds as a constraint to limit the filtering region.It affects to discriminate between the high amplitude of SH waves and Love waves.domain while rejecting or extracting seismic phases The (a), (c), and (e) panels of Fig. 9 show the filtered SP-TFRs of the transverse, radial, and vertical components of the data, respectively.As it can be seen, the focused energy corresponding to the Love waves [shown by blue oval in   To filter the Rayleigh phase, the adjusting parameters

IV. Discussion and Conclusions
We presented a SP-TFF method by combining SP-TFR and EDPA methods as a robust seismic processing tool to separate different phases of seismic waves according to their polarization state.Taking advantage of SP-TFR, high-resolution polarization information is attained to be analyzed by TF-domain polarization attributes for resolving closely spaced seismic events in time and frequency.
Conducting numerical examples on synthetic and real earthquake data, we showed that the SP-TFF can be used as a sophisticated tool to filter elliptical and linear particle motion by designing suitably defined directivity, rectilinearity, and amplitude attributes.Remarkably, not only SP-TFF is efficient to filter Love and Rayleigh waves from the other seismic phases, but it also can handle a more challenging problem of discrimination between the Love and SH and coda waves.This is a promising feature of this method.SP-TFF can find application in various seismic processing methods like anisotropy parameters estimation using the Shear Wave Splitting method [35], surface waves extraction for dispersion curve inversion [36] and sensor miss-orientation test [21], and elimination of the surface waves to extract the coda waves.
The highest computational cost of the algorithm is due to solving the regularized (22).As a result, considerable memory and high computation time are required for a massive input data set.Furthermore, the weakly stationary condition assumption in (5) violates for very low frequencies.Notwithstanding, these are drawbacks of SP-TFF methods.
This paper intended to present the methodology of SP-TFF; only processing earthquake waveforms evaluated the efficiency.The research is ongoing with studying other critical issues, including (a) stability of SP-TFF at different SNR levels, (b) evaluation of the efficiency of SP-TFF for processing ambient noise data, (c) application of SP-TFF for extraction of low amplitude seismic phases, and (d) extraction of surface wave dispersion curves using SP-TFF.

V. Code and data availability
The numerical results from the synthetic and real data examples presented in this paper are reproducible by running a set of computer codes available at the Github account ("Will be inserted when the manuscript is accepted").

Figure 1 .
Figure 1.The background gray color waveforms in top, middle, and bottom panels corresponds to transverse, radial, and vertical components of a 3D synthetic seismogram generated for the source mechanism of Mw = 8.2 south-southwest of Tres Picos, Mexico [see text and Table.I].The foreground black color diagram in panels (a), (c), and (e) are Love wave filtered transverse, radial, and vertical components by using the SP-TFF, and the panels (b), (d), and (f) are Rayleigh wave filtered waveforms of the corresponding components.

Figure 2 .
Figure 2. The TFR of the transverse, radial and vertical components of the synthetic data obtained by the ordinary ST (applied by Pinnegar [15]) are shown in the panels (a), (c), and (e).The corresponding TFR for the SP-TFR method are shown in the panels (b), (d), and (f) [Refer to the text for more explanations].

Figure 3 .
Figure 3. Panels (a) and (c): TFR of measure of SM and Sm axis of particle motion obtained by using the Pinnegar [15] method.Panels (b) and (f): The corresponding TFR of measure of SM and Sm axis of particle motion obtained by by implementing EDPA on the SP-TFR.[Refer to the text for more explanations].
on the SP-TFR of three components.The extracted Love wave is shown in the (a) panel of Fig. 5. Similarly, panels (b) and (c) of this figure show the radial and vertical components of Rayleigh waves.Both the Love and Rayleigh phases have been cleanly extracted from the entire waveform.The extracted surface waves can be used as an input to other

Figure 4 .
Figure 4. Left panel: Adaptively Filtered SP-TFRs of the transverse, radial, and vertical components of synthetic data to eliminate the Love wave.Right panel: Adaptively Filtered SP-TFRs of components to eliminate the Rayleigh wave.

Figure 6 .
Figure 6.The background gray color waveforms in top, middle, and bottom panels corresponds to transverse, radial, and vertical components of Mw = 8.2 earthquake [see text and Table.I].The foreground black color diagram in panels (a), (c), and (e) are Love wave filtered transverse, radial, and vertical components by using the SP-TFF, and the panels (b), (d), and (f) are corresponding Rayleigh wave filtered components.

Figure 7 .
Figure 7. Left panel: The TFR of the transverse, radial, and vertical compenents of Mw = 8.2 earthquake obtained by applying conventional ST.Right panel: The SP-TFR of the components.

Figs. 8
Figs. 8 and 7] are highly damped as a result of filtering, and only a scattered signal corresponding to body and coda waves remained in the transverse component.The SP-TFRs of the radial (panel (c)) and vertical ( panel (e)) components have not significantly affected by filtering.The black color waveforms at the (a), (e), and (e) panels of Fig. 6 depict reconstructed signal for the transverse, radial, and vertical components in the time domain.The results confirm that the SP-TFF filtering significantly canceled the Love wave in the time domain without affecting other phases, including the body and coda waves.To filter the Rayleigh phase, the adjusting parameters

Figure 8 .
Figure 8. Left panel: The TFR of the length of the SM and Sm axis particle motion of the transverse, radial, and vertical components of Mw = 8.2 earthquake obtained by using the Pinnegar [15] method.Right: SM and Sm axes by implementing EDPA on the SP-TFR.

Figure 9 .
Figure 9. Left panel: Adaptive Filtered SP-TFRs of the transverse, radial, and vertical components of the Mw = 8.2 earthquake to eliminate the Love wave.Right panel: Adaptive Filtered SP-TFRs of components to eliminate the Rayleigh wave.

Figure 10 .
Figure 10.Panel (a): extracted Love wave by applying SP-TFF method on the synthetic data.Panel (b) and (c): the radial and vertical component of extracted Rayleigh wave by applying SP-TFF method.
2 earthquake occurred in the 101km SSW of Tres Picos, Mexico, on September 8th, 2017, 04:49:19 (UTC), as a result of normal faulting at an intermediate depth of 47.4 km [see Table.I].The source

Table I
Information of synthetic and real data example correspond to Mw = 8.2 earthquake occurred near Coast Of Chiapas, Mexico recorded at COLA station, IU network Alaska, USA.