An Adaptive Generalized S-Transform Algorithm for Seismic Signal Analysis

As a powerful method in signal processing, time-frequency analysis shows the characteristics of signals in the form of a joint domain distribution of time and frequency. The S-transform is one of the most effective and frequently used algorithms for time-frequency analysis of seismic signals. To address the problems of fixed time window, such as single time-frequency resolution and poor energy aggregation, an adaptive generalized S-transform algorithm is proposed. In this method, we designed a new generalized Gaussian window function with four variable parameters to obtain the optimal time-frequency distribution. The proposed window function uses time-frequency aggregation as the objective function for parameter optimization. Combined with the calculation results of the time-frequency aggregation and time-frequency distribution, we proved that the proposed method has a satisfactory analysis effect and anti-noise performance. In addition, a simulation experiment of the seismic signal is carried out using this algorithm, and the results show that the proposed adaptive generalized S-transform enhances the time-frequency resolution and energy aggregation in seismic signal analysis.


I. INTRODUCTION
Ground penetrating radar (GPR) is often used in underground exploration [1], [2]. GPR transmits high-frequency electromagnetic pulses to the ground, and then the information of the underground medium is resolved through signal analysis of the reflected wave. However, owing to the uneven distribution and the variety of underground media, the signal propagation path and vibration amplitude vary. In addition, there are various noise effects in the seismic wave acquisition, resulting in time-varying nonstationary signals.
Therefore, it is necessary to analyze the seismic signals. The traditional Fourier transform is a global signal analysis method that cannot analyze the local details of the signal. Because of the complexity of GPR data, this method cannot meet actual needs. Time-frequency (TF) analysis methods [3] represent the TF characteristics of signals, and they are used extensively in seismic signal processing [4], [5]. This The associate editor coordinating the review of this manuscript and approving it for publication was Jiafeng Xie. algorithm overcomes the shortcomings of the time-frequency domain separation of Fourier transform. The approach can more precisely analyze the frequency components contained in each time bin and the changes in each frequency component with time. TF analysis of seismic signals and the extraction of feature information [6], [7], [8] are extremely important in geological identification and exploration.
Existing and mature TF analysis algorithms include the short-time Fourier transform (STFT) [9], continuous wavelet transform (CWT) [10], and Wigner-Ville distribution (WVD) [11]. The STFT analysis signal was converted from a 1-dimensional time domain to a 2-dimensional TF joint distribution. The fast Fourier transform intercepts signals through a window function for analysis. Then, the window is shifted, and the process is repeated until the spectrum of the signal is obtained along the time axis. The STFT can adequately describe the local characteristics of seismic signals. However, the STFT is limited in that its window function is fixed for nonstationary signals, leading to a single fixed resolution. Because it does not use a window function, WVD is VOLUME 10, 2022 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ a high-precision TF analysis method with high TF aggregation and is not constrained by the uncertainty principle. However, this was accompanied by a large number of cross interference terms. In practical applications, owing to the existence of interference terms, the obtained spectrum cannot precisely analyze the characteristic information of the signal. Therefore, the pseudo-Wigner Ville distribution (PWVD) [12] and smooth pseudo-Wigner Ville distribution (SPWVD) [13] have been proposed to eliminate cross terms, but the application effect in seismic signals is still unsatisfactory. In 1996, Stockwell proposed the S-transform (ST) [14]. ST is a phase transformation of a CWT with a special basis wavelet. Simple harmonic waves and Gaussian functions are two important components of ST and CWT. In CWT, harmonic and Gaussian windows perform the same scaling and translation. However, in ST, the simple harmonic part can be used as a telescopic transform in the time domain. The Gauss function can be transformed by translation in addition to the stretching transformation. Under the same conditions, ST has a higher TF aggregation than the other TF analysis methods. The window function of ST can adaptively reduce the width with an increase in frequency. Therefore, it has a higher frequency resolution at low frequencies and a higher time resolution at high frequencies. However, the general ST has no parameters that can be used to change the Gaussian window function, which cannot be changed once it is selected. Therefore, many scholars [15], [16], [17] have proposed a generalized S-transform (GST) with different window functions, which is widely used in non-stationary signal processing [18], [19], such as feature extraction of seismic signals [20], [21], reservoir identification [22], [23], and seismic signal compensation [24], [25].
To further optimize the TF resolution and energy aggregation of the ST, a new variable four-parameter Gaussian window function is proposed. The TF aggregation was optimized, and the best window parameters were selected to achieve better TF characterization.
For this study, the aim is to develop a method with higher time-frequency accuracy, we proposed an adaptive generalized S-Transform (AGST) algorithm. Moreover, AGST has better noise resistance and higher time-frequency resolution when analyzing seismic signals.

II. STANDARD S-TRANSFORM
In 1996, Stockwell solved the limitations of CWT and STFT with ST. The basis function of ST is the product of a simple harmonic wave and a Gaussian function, where the simple harmonic wave only performs scaling transforms in time domain and the Gaussian function performs both scaling and translation transforms simultaneously. For signal x(t) ∈ L 2 (R), where L 2 (R) is the finite energy function space, the ST expression is where x is a continuous signal to be analyzed, f is the frequency, and τ is the location of the window on the time axis. The Gaussian window function is expressed as follows: where the standard deviation function is defined as: The time window function of ST contains a frequency factor that is inversely proportional to the time window width. For high-frequency signals, the window function has a small width and good time resolution. For low-frequency signals, the window has a large width and good frequency resolution. Moreover, ST has the stability of non-destructive reversibility and is unaffected by linear cross terms. Standard ST solves the problems of fixed window shape and single resolution that exist in the traditional TF analysis method. However, the window shape of the standard ST has a fixed relationship with the frequency, and the window function only varies linearly with the frequency. Therefore, ST does not have enough satisfactory time-frequency resolution when analyzing more complex signals.

III. MODIFIED S-TRANSFORM A. GENERLIZED S-TRANSFORM
In view of the shortage of ST, Reference [21] proposed adding two window factors to the standard Gaussian window to enhance the flexibility of the window. The standard deviation function is expressed as: where p is the attenuation factor, and λ is the regulation factor. p and λ are generally determined manually, and many attempts to obtain the values of the parameters are required to achieve the best effect. To simplify the tedious manual attempts as well as make the window shape more diverse based on nonlinear frequency changes and more appropriate to the needs of actual signals, based on previous results, this paper proposes to modify the Gaussian window, increase the window control parameters, and optimize the parameters according to the TF aggregation. The modified window function is where the standard deviation function is Therefore, the modified GST expression is × exp(− The normalized parameter makes the window a unit area, so that the amplitude of the new generalized TF transform has the same meaning as the Fourier transform. Because the new window is a function of frequency with multiple parameters, the shape change of the window is related to the window parameters and frequency f , which means that the new GST is a multiresolution signal analysis method. When [A, B, C, D] = [0, 1, 1, 0], the new window is equivalent to the standard Gaussian window. Fig. 1 reveals the influence of the window parameters. In the time domain, the window length becomes narrower as parameters A, B, and D increase, and longer as parameter C increases. The window amplitude increased as parameters A, B, and D increased and decreased as parameter C increased. Parameters C and D had the largest influence on the window and were used to roughly adjust the window shape, whereas parameters A and B had a relatively small influence on the window and were used to fine-tune the window shape. To ensure the energy aggregation of the TF analysis, the window shape cannot be too wide or narrow. Therefore, the σ new (f ) is restricted as follows: where S is nTs, M is lTs, n and l depend on the signal to be analyzed. To select the appropriate window width, after a large number of experiments, n is taken as 10 and l is 1000 respectively in this paper. The four parameters were restricted as follows To further improve the resolution of TF analysis, we propose a new algorithm named AGST. The TF Concentration Measurement (CM) [26], [27] was used for window parameter optimization as the objective function. For the GST of the signal, CM is expressed as where |GST (t, f )| is expressed as This method is analogous to the definition of a 'peak' in statistics. Different window parameters result in different TF distributions and resolutions. When CM is larger, the TF distribution becomes more concentrated. During window parameter optimization, first, the standard deviation function of GST is limited; then, the GST parameters are selected within the limited conditions, the CM of the TF energy aggregation is calculated, and finally, all parameters within the limited conditions are traversed to obtain the CM set. The parameters are then optimized to screen out the window parameters A opt , B opt , C opt , and D opt , where the CM is the largest. Finally, the optimized parameters were used to obtain the generalized TF transform under the condition of an optimal window. This new TF transform with optimized parameters is known as the AGST.

IV. SIMULATION A. MULTICOMPONENT SIGNAL
The Ricker wavelet is the most commonly used form of a seismic wavelet in the construction of seismic signal models.
According to Fig. 2, the STFT has high energy aggregation and time resolution at low frequencies, but the overall TF energy modes are aliased, making it difficult to accurately distinguish the signal characteristics. Compared with the STFT, the CWT improves the overall TF aggregation, but the TF characteristics of the low-frequency and medium-frequency signals remain difficult to distinguish. At the same time, we also compared the GST in Reference [21]. It is worth noting that both λ A and p in Formula (4) need a lot of manual attempts according to the signal to obtain the best analysis effect. When λ = 1.2, p = 1.1, we get the best analysis results in Fig. 2 (e). TF analysis performed by applying the AGST proposed in this study is shown in Fig. 2 Fig. 2 d, Fig. 2 (e), and Fig. 2 (f), we can see that AGST has the best TF analysis effect. Table 1 shows the quantitative results of the TF aggregation comparison between the AGST and traditional TF analysis methods. The results of the CM also show that the AGST in this study has a satisfactory TF aggregation effect. For an SPWVD, although the TF aggregation is significantly improved, it is difficult to distinguish the signal characteristics owing to number of cross-interference terms. This algorithm has a significant effect on synthetic signals, reducing the signal interference phenomenon that exists in traditional time-frequency analysis methods in the lowfrequency region, and has concentrated TF energy, high time resolution, and clear signal distribution. To prove the satisfactory anti-noise performance of AGST, we set the multicomponent Ricker wavelet signal with Gaussian white noise, and AGST was performed on the signal at 10, 5, 3 and 1 dB. The TF distribution of AGST is displayed in Fig. 3. Although the TF energy of the original signal is still concentrated at a low signal-to-noise ratio (SNR), the distribution of the main components in the synthetic signal can be clearly distinguished. The signal did not become slightly blurred until 1 dB, and the time resolution remains relatively accurate. Table 2 shows a comparison of CM between the AGST and TF analysis methods under different signal-to-noise ratios. Combined with the TF distribution in Fig. 3, as is shown that, under the condition of low SNR, the cross term of the SPWVD is too serious, and the frequency is difficult to  recognize. The signals at high frequency and low frequency of the CWT overlap, and the TF characteristics of signals at high frequency are difficult to express using standard ST and STFT. Therefore, the AGST proposed in this study has satisfactory anti-noise performance.

B. WEDGE MODEL
A wedge seismic signal model was designed with an acquisition time of 500 ms and 60 seismic data points were collected. A Ricker wavelet with a dominant frequency of 30 Hz was used to synthesize the seismic records. At the top and bottom of the model, the reflection coefficients were 0.4 and -0.4. Fig. 4 a displays a synthetic seismic record based on the wedge model. Fig. 4 b shows the TF map obtained by applying the standard ST to the synthetic seismic record, with a TF aggregation CM = 0.0 130. The interference between the synthetic signals is severe in and near the 10th, 30th and 50th seismic data of the wedge model, and the recognition effect of thin-layers is fuzzy. Fig. 4 (c) shows the TF aggregation CM = 0.0 162 obtained by applying the AGST proposed in this paper to the synthetic seismic record. The parameters of the 60 channels data are optimized, and a TF analysis is performed. Subsequently, the dominant frequency was extracted, and an amplitude slice was created. From the amplitude spectrum used to distinguish the thinlayer interface, it clearly shows that the resolution is markedly enhanced and the thin-layer identification becomes more accurate.

C. HORIZONTAL THIN-LAYER
A horizontal thin-layer model was designed. The model has 60 channels and an acquisition time of 600 ms. There is a single reflection surface from channel 1 to channel 60 at 300 ms, with reflection coefficients of 0.4280 ms and 350 ms. Horizontal thin-layers with a reflection coefficient of -0.4 can be found at channels 25 to 60 and channels 7 to 30. The seismic profile was obtained by the convolution of the horizontal thin-layer model with a Ricker wavelet. As shown in Fig. 5 (a), the dominant frequency of the Ricker wavelet was 60 Hz. One of them is extracted from the profile to calculate the dominant frequency, and the TF distribution in Fig. 5 (b) is obtained by the standard ST, with the corresponding CM = 0.0 175. Fig. 5 (c) shows the TF distribution map obtained by the AGST using the method in this study, and its TF aggregation was increased to CM = 0.0 217. Comparing the two figures, Fig. 5 (c) clearly shows the position information and shape characteristics of the thin-layer. Even if the depth difference of the thin-layer between tracks 25 to 60 was only 20 ms, the algorithm could still accurately identify the thin-layer.

D. NATURAL SEISMIC SIGNAL
A natural seismic wave with a time sampling interval of 0.001 s and 601 data points was selected for analysis, see Finally, we extracted a section of natural seismic wave data from 80 seismic channels and analyzed the 10th seismic signal. The sampling interval was set to one, and there were 256 sampling points. There are two obvious wavelets with dominant frequencies of 15 Hz at the 98th and 110th sampling points, and some small wavelets are scattered at other points. Fig. 7 (a) shows the amplitude spectrum of the extracted signal. Fig. 7 (b) shows the TF distribution obtained by performing the standard ST. By applying the AGST, TF distribution is obtained as is shown in Fig. 7 (c). When the two images  are compared, it is shown that AGST significantly improves the time resolution and TF aggregation, as well as the ability to distinguish the distribution of seismic wavelets.

V. CONCLUSION
An AGST for seismic signal analysis is proposed in this paper in response to the shortcomings of traditional time window functions, such as single window shape transformations and difficulty in meeting the needs of TF analysis for complex non-stationary signals. To ensure that the AGST has a higher TF aggregation and TF resolution, the algorithm incorporates a new Gaussian window function with four parameters that makes the window more flexible and uses TF aggregation as the objective function for parameter optimization. The simulation results show that AGST outperforms the traditional method in the thin-layer discrimination of seismic signals, and the recognition is more accurate. In addition, the algorithm exhibits an improved anti-noise performance.