Seismic Time-Frequency Analysis Based on Time-Reassigned Synchrosqueezing Transform

In this paper, we propose a novel seismic time-frequency analysis method via the time-reassigned synchrosqueezing transform (TSST), in which the time-frequency coefficients are reassigned in the time direction rather than in the frequency direction as the short-time Fourier-based synchrosqueezing transform (FSST) does. Such a technique can not only produce a highly concentrated time-frequency representation (TFR) for a wide variety of strongly frequency modulated signal, but also allow for the reconstruction of the modes with a high accuracy. Numerical experiments on synthetic signal and field data demonstrate the effectiveness of this new method, and show that the proposed method is more suitable for extracting seismic time-frequency feature and identifying the thin layers compared with the traditional FSST, which offers the potential in highlighting subtle geological structures.


I. INTRODUCTION
Extracting useful information from large amounts of recorded data is important for a large number of real-world applications [1]- [3]. Time-frequency analysis (TFA) method can be applied to characterizing the non-stationary and nonlinear seismic signal [4]- [8]. Traditional TFA approaches can be roughly divided into two categories. The first one is the linear transform based methods such as short-time Fourier transform (STFT) [9], [10] and continuous wavelet transform (CWT) [11], [12]. Although both transforms are invertible and allow one to cope with the non-stationary signal, they suffer from the Heisenberg uncertain principle, which limits the adaptivity of the method itself [13]. The second type is based on quadratic transform, e.g. Wigner-Ville distribution (WVD) [14], [15]. Such method is characterized by high time-frequency resolution, however, the strong interference reduces the readability of TFR and the WVD is not invertible.
To solve those problems, many efforts have been made. Auger and Flandrin (1995) provided a means of improving the TFR readability, termed as reassignment method (RM) [16], which transforms the time-frequency coefficients The associate editor coordinating the review of this manuscript and approving it for publication was Wei-Wen Hu .
to their corresponding center of gravity along both the time axis and the frequency one. Unfortunately, the reassigned transform is no longer invertible and does not allow for mode reconstruction. Matching pursuit (MP) improves the time-frequency resolution greatly by decomposing a seismic trace into a series of wavelets which belong to a comprehensive dictionary of functions [17], [18]. However, the high resolution is achieved at the cost of increased computational burden. Han and van der Baan (2013) illustrated the suitability of EMD-based algorithms for seismic time-frequency analysis [19]. However, a lack of solid mathematical foundation still limits the wide application of such methods. Daubechies and Maes [20] developed a phase-based technique, called synchrosqueezing transform (SST) [20], [21]. It aims at improving the TFR by combining CWT and reassignment technique, and allowing for mode retrieval [22]. Meanwhile, Thakur and Wu (2011) proposed an extension of SST to the STFT setting, i.e., so-called FSST [23], which yields an excellent TFR and mode reconstruction, and is proven to be robust to noise [24]. Subsequently, several similar synchrosqueezing transforms based on different frameworks are emerging, such as synchrosqueezing S-transform (SSST) [25], [26] and synchrosqueezing generalized S-transform (SSGST) [27], [28], which have been successfully used to analyze seismic data. However, it is worth noting that the aforementioned techniques mostly reassign the time-frequency energy in the frequency direction, so that the strongly modulated multicomponent signal with fast varying instantaneous frequency (IF) may not be satisfactorily dealt with [29], [30].
In this paper, we proposed a new seismic TFA method based on the TSST, in which the reassignment operation is implemented in the time direction [31], rather than frequency direction as the SST and FSST do. This results in the perfect energy concentration on the TFR for the signal with fast varying IF. Meanwhile, the proposed method allows for reconstruction of each component making up of the input signal and the computation complexity remains the same as the FSST. In addition, it is noteworthy that our motivation is different from that presented in [32], where the key idea is to only hold the time-frequency information most related to time-varying features of the signal and to remove most smeared time-frequency energy. We mainly focus on the time-frequency coefficients reassignment in the time direction. The contributions of the paper can be summarized as below: (1) we first investigate the time-reassigned property of TSST and achieve seismic time-frequency feature extraction via the TSST, (2) our method shows the excellent potential in subtle geological structures characterization compared to the conventional STFT and FSST approaches.
The outline of this paper is as follows: in Section II, some fundamental definitions on STFT and FSST are firstly recalled. Then, we introduce the TSST and Rényi entropy, and describe the practical implementation of the proposed method. Section III delivers numerical results on both synthetic and real examples, comparing the proposed method with standard STFT and FSST. Finally, a discussion of the results is given in section IV, and conclusions are drawn in Section V.

II. THEORY A. SHORT-TIME FOURIER TRANSFORM
The (modified) STFT of a given signal f is defined as follows: where g is a real-valued and symmetrical window function, and V f (t, η) 2 is called the spectrogram of signal f . Furthermore, the original signal f can be reconstructed from its STFT by the following equation: B. STFT-BASED SYNCHROSQUEEZING TRANSFORM (FSST) The purpose of FSST is to sharpen the blurry STFT representation by reassigning the time-frequency coefficients of V f (t, η) from the (t, η) to t, ∧ ω f (t, η) , so that the energy-concentrated TFR is obtained.
The local instantaneous frequency estimation at time t and frequency η based on the (modified) STFT, where ∂ t is the partial derivative with respect to t, R {•} denotes the real part of a complex number. The FSST using the synchrosqueezing operator is defined as follows: where δ is the Dirac distribution, ω is the frequency variable, and γ is some threshold.
The ith mode can be approximately retrieved by: where is the phase function of the ith mode of an AM-FM signal defined by: where A i (t) and ϕ i (t) are the instantaneous amplitude and instantaneous phase of the ith mode, respectively. It's important to note that the FSST squeezes the time-frequency coefficients along the frequency direction, which is suitable for a class of weak frequency modulated signals. However, for strong frequency modulated ones, it usually does not give a satisfactory result. Thus, a more advanced TFA method is required.

C. TIME-REASSIGNED SYNCHROSQUEEZING TRANSFORM (TSST)
Unlike FSST, the TSST reassigns the time-frequency energy in the time direction, and preserves the invertible properties [31]. Now, we consider an impulse signal defined as follows: The standard STFT of signal f δ (τ ) with the real symmetric window g (τ ) is described as: (8) where u and η are the time variable and the frequency variable, respectively.
Then we have Thus, . (10) Note that Eq. (10) is different from the conventional group delay (GD), in which the (modified) STFT is used. Herein, we define the (generalized) GD estimator of the impulse signal f δ (τ ) is defined as: The TSST is defined as: Eq. (12) shows that the TSST reassigns the information from the (u, η) plane Finally, the signal reconstruction from TSST is achievable by:

D. Rényi ENTROPY
In the paper, we employ the Rényi entropy to assess the distribution concentration of a TFR [33]. The α-Rényi entroy with regard to a nonzero function s is described as follows: where s α := |s (x)| α dx 1 / α and α > 0. Generally speaking, α > 2 is recommended for TFR measure [34], and we choose α = 3 in this paper. A lower Rényi entropy means the more concentrated TFR. The measure of distribution concentration is defined as: and, where u denotes the time instant, d is the size of the neighborhood, and I u : Herein, the parameter d is set to 0.12 by several trials because of its insensitivity to the final result.   Figure 1 illustrates the proposed seismic time-frequency feature extraction workflow, and the detailed implementation steps are summarized as follows:

E. PROCEDURE OF THE PROPOSED METHOD
Step 1: Determine the window duration parameter and the threshold applied on window value.
Step 2: Perform TSST on the input signal to calculate timefrequency coefficients.
Step 4: Select the characteristic frequency within the frequency band of seismic signal.
Step 5: Repeat steps (1) through (4) for all traces in the seismic data.
In the TSST algorithm, the window duration parameter (σ ) and threshold (γ ) are two key parameters. The σ is often related with the time-frequency energy concentration, while the γ controls the computation accuracy with mode  retrieval. It could result in a blurred TFR and the inaccurate signal reconstruction if the parameters are not appropriately selected. In practice, the σ and γ can be obtained by several trials.

III. EXAMPLES
A. SYNTHETIC DATA First, we investigate the performance of the proposed method using a synthetic signal (Figure 2(a)), which is composed of two components: We calculate the time-frequency maps using the STFT, FSST and TSST, respectively. We select the window duration parameter for the proposed method to σ = 18, and a threshold parameter γ is set as 0.00001. The results are shown in Figure 3. It can be observed that the STFT suffers from a poor time-frequency resolution due to its fixed window (see Figure 3(a)), the FSST further improves the time-frequency energy concentration of the STFT, but the time-frequency coefficients still spread out along the IF trajectories, which result in a blurred TFR (see Figure 3(b)). By contrast, the TSST achieves the best TFR result (see Figure 3(c)). This is because the TSST reassigns the time-frequency energy in the time direction, which is more beneficial to characterizing the signal made up of strongly modulated modes. For a better view, we enlarge a local area from Figure 3 (see the pink rectangle), which is plotted in Figure 4. As illustrated in Figure 4, the time-frequency energy is perfectly concentrated by the use of TSST in comparison with STFT and FSST, and the TSST can better delineate each component constituting the synthetic signal.
To demonstrate the TSST's robustness with respect to noise, a series of noisy signals with the different signal-tonoise ratios (SNRs) ranging from 0 to 30 dB are generated by adding the white noise into the synthetic signal. Figure 2(b) shows a result with a SNR of 2 dB, the time-frequency maps obtained by using STFT, FSST, and TSST are present in Figure 5, and the corresponding enlarged versions are displayed in Figure 6. Our method is run with the same parameters as that in Figure 3. It can be clearly seen that the noise has a strong influence on time-frequency feature in the STFT result (see Figures 5(a) and 6(a)). The FSST produces a blurry TFR (see Figures 5(b) and 6(b)), which is not good for IF estimation. However, the TSST still exhibits the better time-frequency energy concentration (see Figures 5(c) and 6(c)), in which the local feature can be characterized well. In this sense, the TSST has better noise robustness.
Then, we employ the Rényi entropy to evaluate the performance of different TFA methods quantitatively. A lower Renyi entropy value means a more concentrated TFR. As reported in Figure 7, the TSST result has the lowest Rényi entropy among all mentioned TFA methods, which indicates that it can produce the most energy-concentrated TFR. VOLUME 9, 2021 FIGURE 5. Time-frequency maps of the noisy signal (Figure 2(b)). (a) STFT, (b) FSST, and (c) TSST. In the results of STFT and FSST, the noise has a greater impact on time-frequency feature compared to that from TSST.

B. FIELD DATA
In this section, a field data (Figure 8(a)) is employed to test the proposed method and compare to the STFT and FSST approaches. This data is composed of 150 traces with 500 samples per trace and a sampling interval of 2 ms. We first take the seismic trace 75, which is plotted in Figure 8 (b), and apply STFT, FSST, and TSST on the data, the results are shown in Figure 9. For this example, the parameters σ and γ are selected as 30 and 0.00005 in the proposed method, respectively. As can be clearly seen, all TFRs exhibit some similar features including the strong spectral energy at 0.20s, 0.33s and 0.73s, and a decrease in spectral energy with time. More specifically, the STFT produces a poorer TFR owing to its intrinsic shortcoming. The FSST and TSST present more detailed information than the STFT due to the higher time-frequency resolution of both methods. However, the TSST is obviously superior to the FSST in time localization, which facilitates the accurate characterization of transient spectral anomalies. The key reason is that the reassignment operation is implemented along the time direction in the algorithm of TSST, rather than the frequency direction as the FSST does. In addition, we also test the reconstruction accuracy with regard to TSST, which is shown in Figure 10. As illustrated, there is nearly no observable difference between the original and reconstructed signals. The TSST has a nearly perfect reconstruction with an overall   reconstruction error (b). The TSST provides a good estimation. There is no obvious difference between the recovered signal and the original one, and the reconstruction error is negligible.
Finally, we apply the three methods to all traces in the field data and extract the 40Hz and 55Hz frequency slices. Figure 11 shows the STFT, FSST, and TSST results. It can be clearly seen that both of STFT and FSST show a set of relatively blurred frequency slices ( Figures. 11(a), (b), (d), and (e)), and cannot effectively highlight the structural and geological information. The STFT brings out a whole spectral content with the poorer time-frequency resolution due to the Heisenberg uncertainty principle (Figure 11(a)), while The squeezing along the frequency axis heavily reduces the time resolution of the FSST, so that it barely identifies the thin layers between 0.2s and 0.35s (Figure 11(b)). By contrast, The TSST depicts a cleaner spectral representation with higher time-frequency resolution and better geological structures. Furthermore, the main reflectors located at 0.1s, 0.2s and 0.35s (indicated by the arrows in Figures 11(c) and (f)) are well identified, which is very crucial in seismic interpretation.
In addition, we also evaluate the computation costs of the aforementioned methods operating on synthetic and real examples, which are listed in Table I. As reported in Table I, the STFT takes the shortest time, the FSST and TSST are almost equivalent in time while the latter is a little shortern, which means the TSST remains the same level of computation complexity as the FSST. All experiments are done on a PC station equipped with an Intel Pentium Duo Core CPU clocked at 2.53 GHz and 2 GB of RAM.

IV. DISCUSSION
TFA is one of widely used tools for characterizing the time-varying feature of a non-stationary signal. The STFT usually fails to produce a satisfactory result due to the Heisenberg uncertainty principle. The original SST and its extended version such as FSST reassign the time-frequency coefficients in the frequency direction, thus, they are more suitable for analyzing the weak frequency modulated signal. In the TSST method, the reassignment operation is implemented in the time direction and the GD estimator is employed to replace the instantaneous frequency estimator in the original SST, which enable it to generate a highly concentrated TFR for a wide variety of strongly time-varying signal. From the comparison with STFT and FSST methods, it can be clearly seen that the TSST can better achieve the feature extraction based on the temporal localization. In addition, it is worth noting that the signal reconstruction is the main advantage of TSST over RM, which means that we can make the most of the invertible property of TSST for signal denoising and mode retrieval.
In the algorithm of TSST, two parameters that need to be predefined, window duration parameter σ and threshold γ , usually play an important role. The window duration parameter determines the energy concentration of a TFR, the threshold parameter controls the signal reconstruction accuracy. Now, we take the synthetic signal in Figure 2(a) for an example to test the above-mentioned parameters. In Figure 12, the effect of the window duration parameter variations is shown on time-frequency energy concentration. As can be observed from Figure 12, a larger value will result in inaccurate IF estimation (see Figure 12(a)) while a smaller value may reduce the energy concentration of a TFR (see Figure 12(b)). Thus, a suitable compromise solution is often   required by several trials in practice. Then, we calculate the MSE for various values of the threshold parameter. The obtained result is shown in Figure 13. As can be clearly seen, the MSE has a slightly increasing against the large variations of the threshold parameter, hence, the TSST is not very sensitive to the threshold parameter.
In terms of computational cost, experimental results indicate that the TSST and the FSST are almost equivalent and slightly higher than the STFT. Although the TSST shows great advantage in characterizing strongly time-varying signal, it does not seriously reduce the computational efficiency. Therefore, the TSST can be applied to seismic processing in the real world.

V. CONCLUSION
In this paper, we present a new approach for seismic time-frequency analysis method based on the TSST. Such a method allows us to better handle a wide variety of strongly time-varying signal. The advantage of the proposed method is demonstrated by numerical experiments both for synthetic and real data. It successfully produces a highly concentrated TFR by reassigning the time-frequency coefficients in the time direction compared with other methods such as the STFT and FSST, while allows for better inevitability of the TFR. Thus, the proposed method is more suitable for characterizing the non-stationary seismic signals with strongly time-varying feature, which could probably help to better understand the subtle geological structures. Future works include the behavior analysis of this method when applied to seismic data with complex noise, and developing new applications. YANG LIU received the B.E. degree in safety engineering from Beijing University of Chemical Technology, in 2020, where he is currently pursuing the M.E. degree. His research interests include seismic signal analysis, deep learning, and intelligent fault diagnosis. VOLUME 9, 2021