1 ANNALS OF GEOPHYSICS, 65, 4, DM423, 2022; doi:10.4401/ag-8761 O P E N A C C E S S Improved generalized S‑transform deconvolution for non‑stationary seismic data Chao Sun1, Dengke He*,2, Lijun Shen1 and Liang Sun2 (1) Shandong Provincial Research Institute of Coal Geology Planning and Exploration, Jinan 250104,China; (2) College of Geoscience and Surveying Engineering, China University of Mining and Technology (Beijing), Beijing 100083, China Article history: received November 20, 2021; accepted March 23, 2022 Abstract Improving the vertical resolution is one of the significant tasks for seismic data processing. Most traditional resolution-enhancement techniques assume that the seismic wavelet is time-invari- ant. However, the seismic wavelet varies with seismic wave propagation in the subsurface. To solve this issue, a new spectral-modeling method is proposed to extract the time-varying wavelet using improved generalized S‑transform (IGST) and higher‑order Fourier series. The IGST based on modified time‑window function can effectively improve the resolution of the time‑frequency (t‑f ) spectrum. The high‑order Fourier series is used to fit on the logarithm t‑f spectrum and achieve the high‑precision time‑varying wavelet. The proposed method is composed of four steps in the implementation. Firstly, the seismic data is decomposed by the IGST and converted to the logarithm t‑f domain. Secondly, the time‑varying wavelet spectrum is modeled at each time sample using a higher‑order Fourier series. Thirdly, the boxcar smoothing method is used to smooth the time‑vary- ing wavelet spectrum and extract the time‑varying wavelet with Hilbert transform. Finally, using the time‑varying wavelet spectrum to spectrally balance seismic data to flatten the seismic response. Synthetic and field data examples demonstrate the feasibility of the proposed method in improving the signal‑to‑noise ratio and enhancing the vertical resolution. Keywords: Deconvolution; Time-varying wavelet; Improved generalized S-transform; High-order Fourier series; spectral-modeling method 1. Introduction Improving vertical resolution plays a significant role in seismic processing [Vardar and Alpar, 2016; Zhou et al., 2016; Hao et al., 2019, Lari et al., 2021]. Many methods, such as wiener deconvolution [Robinson and Treitel, 1967], spectral modeling [Rosa and Ulrych, 1991; Xue et al., 2021; Li et al., 2022], and sparse‑spike deconvolution [Velis, 2007; Pereg et al., 2017], have been proposed to improve the vertical resolution. These methods tend to be according to that the wavelet is time‑invariant, which meets the condition of the stationary convolutional model. In fact, the seismic wave propagates through a heterogeneous, viscous-elastic media of subsurface which brings both dispersion Chao Sun et al. 2 and attenuation effects, leading to the seismic wavelet’s variation during seismic wave propagation, forming a non‑stationary seismic data [Margrave et al., 2004; Radad et al., 2015; Zhang et al., 2019]. Many researchers [Bickel and Natarajan, 1985; Hargreaves and Calvert, 1991; Hargreaves 1992; Wang, 2002, 2006; Zhang et al., 2007; Zhao et al., 2019] have provided the inverse Q filtering technique to address the non‑stationary seismic data. But evaluating the accurate Q‑factor for the inverse Q filter is not easy [Corrales et al., 2014; Zhang et al., 2015]. Margrave et al. [2011] presented a t‑f deconvolution methodology that calculates the time‑varying wavelet spectra smoothing the Gabor spectrum of the non‑stationary seismic data. The Gabor transform, a t‑f analysis technology, is used to analyze the non‑stationary signal in the t‑f deconvolution method. In fact, numerous t‑f analysis technologies [Gabor, 1946; Morlet et al., 1982; Stockwell et al. 1996], such as the wavelet transform (WT), short-time Fourier transform (STFT), and S-transform (ST), can depict the frequency variation with time of the signal. The ST, which avoids many disadvantages of the Gabor transform, STFT, and WT, is a popular t‑f analysis technique. Moreover, ST is multiresolution, linear, reversible and lossless. Due to the mother wavelet of S‑transform is fixed, it is difficult to adjust it as needed. Mansinha et al. [1982] defined mother wavelets of ST according to the seismic data requirements and expands it to the generalized S‑transform (GST). Inspired by the study work of Rosa and Ulrych [1991], the reflectivity is generally white and the seismic wavelet spectrum is smooth. The time‑varying wavelet spectrum can be estimated by the spectral modeling method after smoothing the t‑f spectrum. Zhou et al. [2016] used a polynomial to smooth the t‑f spectrum at each time sample for estimating the time‑varying wavelet. Wang et al. [2017] used a filter to smooth the secondary t‑f spectra of non‑stationary seismic data. The spectrum can become smooth after being constrained by these spectral modeling methods. In this paper, we proposed an improved generalized S-transform (IGST) to calculate the t-f component of seismic data based on the non‑stationary hypothesis. Then the time‑varying wavelet spectrum is estimated using the higher‑ order Fourier series in the logarithm t‑f domain. The two techniques are employed to deconvolve the non‑stationary seismic data for improving resolution. 2. The improved generalized S‑transform (IGST) Fourier transform, an important tool for seismic signal analysis, establishes the relationship of the signal’s time and frequency domain. However, they are independent of each other, and cannot represent the amplitude spectrum and phase spectrum at a specific time. To deal with the problem, Gabor [1946] proposed that the signal can be divided into different time windows, and each time window can be analyzed by Fourier transformation, which is Gabor transform and can be written as: (1) where is the input signal in the time domain, means the signal in the t-f domain, 𝜏 and 𝑓 are the parameters of time shift and frequency value respectively, and represents the Gaussian function or the window function and can be expressed as (2) However, once 𝜎 was selected, the time resolution and frequency resolution are invariable. Stockwell et al. [1996] proposed the ST by setting and then (3) Non-stationary seismic deconvolution 3 The inverse ST can recover the signal and is expressed as (4) Due to the mother wavelet of ST is fixed, the low‑frequency signal’s time resolution is generally poor. Therefore, if we denote (5) and substitute equation (5) and equation (2) into equation (1), we can obtain the improved generalized S‑transform (IGST) (6) The parameter of 𝜆 can influence the time width and bandwidth of the window. A larger parameter of 𝜆 the larger the time width, the smaller the bandwidth and the higher the frequency resolution. The parameter of 𝑝 can enhance the impact on the frequency parameter 𝑓 for 𝜆. The special parameter 𝑏 can influence the position on the frequency coordinates, which can improve the time frequency resolution low frequencies. When 𝜆 = 1, 𝑝 = 1, 𝑏 = 0 in equation (6), the IGST becomes the S‑transform. Similarly, when 𝑝 = 0, 𝑏 = 0 in equation (6), the S‑transform turns into Gabor transform. The IGST overcomes the restriction of ST where the Gaussian window function cannot be changed, which can improve the time resolution of the low‑frequency signal by adjusting the parameters 𝜆, 𝑝, and 𝑏. Here, we produce a synthetic signal 𝑥(𝑡) to illustrate the performance via varying the parameters of the IGST. The synthetic signal 𝑥(𝑡) contains four kinds of cosine functions, and their frequencies are 20 Hz, 40 Hz, 80 Hz, and 100 Hz, respectively. The synthetic signal 𝑥(𝑡) is expressed as (7) Figure 1 shows the synthetic signal 𝑥(𝑡). Figure 2a represents the actual t‑f spectrum, which precisely describes the t‑f distribution of the synthetic signal. Figure 2b displays the t‑f spectrum using the ST. Figure 2c is the t‑f spectrum using the IGST. The frequency resolution and time resolution of the ST is poorer compared with the actual t‑f spectrum. The low‑frequency signal in figure 2b has lower time resolution than that in figure 2c. The high‑frequency signal in figure 2b has lower frequency resolution than that in figure 2c. Thus, the IGST has better t‑f resolution than ST. Moreover, we can choose IGST parameters flexibly according to our own needs to obtain the relatively optimum t‑f spectrum. Figure 1. The synthetic signal with four main frequencies. Chao Sun et al. 4 Figure 2. The t‑f analysis example. Left: The true t‑f spectrum. Center: The t‑f spectrum by ST (𝜆 = 1, 𝑝 = 0, 𝑏 = 0). Right: The t-f spectrum by the IGST (𝜆 = 1, 𝑝 = 0.85, 𝑏 = 6). 3. Methods 3.1 Non‑stationary seismic data in logarithm t‑f domain Based on the IGST, we can get the suitable t-f spectrum of the non-stationary data by: (8) where is the t-f components of the non-stationary data 𝑥��, IGST[.] denotes IGST by the equation (6). The amplitude spectrum and phase spectrum are respectively expressed as (9) (10) where and indicate the real and the imaginary part of the t-f spectrum . Next, we convert the amplitude spectrum to the log domain (11) where denotes the logarithm t‑f spectrum. ln[·] means the natural logarithm of the amplitude spectrum . Non-stationary seismic deconvolution 5 3.2 Spectral modeling of non‑stationary seismic data in logarithm t‑f domain The spectrum modeling method was proposed by Rosa and Ulrych [1991] for extracting the wavelet in the frequency domain of seismic trace. It is assumed that the wavelet’s amplitude spectrum is a smooth unimodal curve, and the reflectivity’s amplitude spectrum is generally white due to the geological structures commonly meet the statistical model. The amplitude of the wavelet can be acquired by smoothing the amplitude spectrum of stationary data in the logarithm frequency domain. Since the non‑stationary seismic trace can be considered stationary at each time sample in the t‑f domain. We can extract the wavelet’s amplitude spectrum at each time sample of the logarithm t‑f domain by smoothing the amplitude spectrum of the non‑stationary data at each time sample. Finally, the time‑varying wavelet’s amplitude spectrum can be obtained. Rosa and Ulrych [1991] provided a low order polynomial to fit the wavelet amplitude spectrum. Inspired by Fourier analysis, the Fourier series is utilized to fit the wavelet’s amplitude spectrum in logarithm t‑f domain, which can be expressed as (12) where (13) is the L2 norm, 𝜀 indicates the misfit error waited to be optimized, 𝑎0, 𝑎�, and 𝑏� are coefficient of Fourier series varied with time, N is the order of Fourier series, 𝑘� is a positive constant. The 𝑎0, 𝑎� and 𝑏� can be calculated by the nonlinear least square methodology. Equation (12) represents the FSF method, which can estimate the time‑ varying wavelet spectrum. 3.3 3.3. Spectral smoothing of non‑stationary seismic data in t‑f domain The estimated amplitude spectrum of the time-varying wavelet shows that some components are inconsistent with that of the exact amplitude spectrum. The boxcar smoothing method and hyperbolic smoothing method are used to deal with the problem. The boxcar smooth method can be regarded as a 2D boxcar convolution with the estimated amplitude spectrum, which is: (14) where 𝐴 is the estimated amplitude spectrum, 𝐵 represents the convolution operator, 𝑀 represents the length of time series, 𝑁 indicates the length of frequency series. The convolution operator smooths along with the time and frequency directions and forms a transition zone at the discontinuous points of the amplitude spectrum. The hyperbolic smoothing method uses the constant Q‑value model. In the case of the constant Q‑value model, the attenuation interface of the amplitude spectrum is distributed along the hyperbola, if the amplitude spectrum is divided into N parts according to hyperbola, each part corresponds to different attenuation interface, the middle part of a two‑phase hyperbola is called a hyperbolic belt. Thus, the result of the hyperbolic smoothing method can express as [Margrave et al., 2011] (15) where 𝐼� indicates the nth hyperbolic belt, 𝑏(𝑓0) is a convolver and • is an convolution operation in frequency domain. Chao Sun et al. 6 3.4 Deconvolution in t‑f domain After getting the estimated amplitude spectrum of the time-varying wavelet, it is necessary to design the deconvolution operator. We design a deconvolution operator at each time sample in the t‑f domain. Firstly, the maximum 𝐴max at each time sample window is extracted, which is expressed as (16) Secondly, according to the minimum phase hypothesis, the phase of the estimated time-varying wavelet can be represented as: (17) Where 𝐻[·] represents the Hilbert transform, 𝜇 is a very small constant value, which ensures the stability of deconvolution. Thirdly, the deconvolution operator is (18) Finally, the deconvolution result in the t-f domain can be obtained, which is expressed as (19) Then the time‑domain trace can be reconstructed by the inverse transformation of IGST. 4. Synthetic data example For the purpose of demonstrating the advantage and effectiveness of the IGST deconvolution methodology proposed in this paper, we design a synthetic model. The model consists of a series of horizontal strata and a deeper wedge, which velocity range is between 1600 m/s and 2800 m/s. Figure 3 shows the vertical reflection coefficient and Figure 3. Reflection coefficient model (Left) and its non‑stationary seismic section (Right). Non-stationary seismic deconvolution 7 its non‑stationary seismic data that is convoluted by the vertical incidence reflectivity, the minimum phase wavelet which center frequency is 45 Hz, and attenuation function (Q = 50). The trace interval is 1 m, and the time interval is 1 ms. In order to compare the processing effect, the synthetic section is implemented by the Gabor deconvolution and the IGST deconvolution respectively. Seeing in the Figure 4, the result section by the IGST deconvolution has better resolution than that by the Gabor deconvolution. The wavelength about the coda wave from the event after the Gabor deconvolution is longer than that from the event after the IGST deconvolution. The comparison indicates the IGST deconvolution is more effective than the Gabor deconvolution for non‑stationary seismic data. Figure 4. The result section of Gabor deconvolution (Left) and IGST deconvolution (Right). 5. Field data testing To examine the effectiveness of the IGST deconvolution methodology, post‑stack seismic data is provided in this paper. After implemented the mainly process techniques including static correction, surface‑consistent correction, surface wave suppression, normal move‑out and stack, the resulting seismic profile is dealt by the Gabor deconvolution and the IGST deconvolution respectively for comparison. Seeing in the left subplot of the figure 5, the original seismic events within a yellow box cannot be well distinguished, which means that the resolution of the original seismic profile is low. The Gabor deconvolution is used to process the original seismic section, and the outcome is displayed in the center subplot of the figure 5b, which has a higher resolution than the original section. By comparison, the result section after IGST deconvolution, which is shown in the right subplot of the figure 5, has the highest resolution than that after the Gabor deconvolution, especially the events within the yellow box are more detailed and smoother. Figure 5. Original seismic section (Left), Result sections after Gabor deconvolution(Center) and after IGST deconvolution (Right). Chao Sun et al. 8 6. Conclusions The IGST deconvolution method is proposed to evaluate and remove the effect of time-varying wavelets on non-stationary seismic data, which can improve the signal-to-noise ratio and enhancing the vertical resolution of seismic section. According to the new method, the IGST is firstly used to calculate the t‑f distribution of non‑stationary seismic data, and then convert that to logarithm t‑f domain. Secondly, The higher‑order Fourier series fitting is used to smooth the amplitude spectra at every time sample in the logarithm t‑f domain and estimate the time‑varying wavelet. Finally, the deconvolution operator will be executed. Due to employ the high resolution features about the IGST for time frequency analysis and the higher‑order Fourier series for wavelet fitting, the new deconvolution method proposed can enhance the vertical resolution and compensate attenuation effect effectively, which is tested and exhibited by the synthetic seismic data and actual seismic data in the paper. Author Contributions: All authors are equal contributors to this work. Funding: This research was funded by the Natioanl 111 Project (B18052), China University of Mining and Technology (Beijing) Yueqi Distinguished Scholar Program (2019JCA01) and Shandong Province Geological Exploration Fund Project (Lukanzi (2019) No. 13, (2020) No. 19). References Bickel, S.H. and R. R. Natarajan (1985). Plane‑wave Q deconvolution, Geophysics, 50, 1426‑1439. Corrales, Á., F. Cabrera and L. Montes (2014). Enhancing time resolution by stabilized inverse filter and Q estimated on instantaneous spectra, Earth Sci. Res. J., 18, 63‑67. Gabor, D. (1946). Theory of communication. Part 1: The analysis of information, Journal of the Institution of Electrical Engineers‑Part III: Radio Communication Engineering, 93, 429‑441. Hao, Y., H. Huang, J. Gao and S. Zhang (2019). Inversion‑Based Time‑Domain Inverse Q‑Filtering for Seismic Resolution Enhancement, IEEE Geosci. Remote S., 16, 1934‑1938. Hargreaves, N.D. and A. J. Calvert (1991). Inverse Q filtering by Fourier transform, Geophysics, 56, 519‑527. Hargreaves, N.D. (1992). Similarity and the inverse Q filter: Some simple algorithms for inverse Q filtering, Geophysics, 57, 944‑947. Lari, H. H. and A. Gholami (2021). Nonstationary deconvolutive Randon transform, Geophysics, 86, 1JA‑X4. Li, Y.X., J. S. Shen, W. Y. Cai and X. H. Zhou (2022). Fractured formation evaluation by seismic attenuation derived from array acoustic log waves based on modified spectral ratio method and an extended Biot’s poroelastic model, J. Petrol. Sci. Eng., 209, 109838. Mansinha, L., R. G. Stockwell, R. P. Lowe, M. Eramian and R. A. Schincariol (1997). Local S Spectrum Analysis of 1‑D and 2‑D Data, Phys. Earth Planet. In., 103, 329‑336. Margrave, G.F., M. P. Lamoureux and D. C. Henley (2011). Gabor deconvolution: Estimating reflectivity by nonstationary deconvolution of seismic data, Geophysics, 76, W15‑W30. Morlet, J., G. Arens, E. Fourgeau and D. Glard (1982). Wave propagation and sampling theory—Part I: Complex signal and scattering in multilayered media, Geophysics, 47, 203‑221. Pereg, D., I. Cohen and A. A. Vassiliou (2017). Multichannel sparse spike inversion, J. Geophys. Eng., 14, 1290‑1299. Radad, M., A. Gholami and H. R. Siahkoohi (2015). S‑transform with maximum energy concentration: Application to non‑stationary seismic deconvolution, J. Appl.Geophys., 118, 155‑166. Robinson, E. A. and S. Treitel (1967). Principles of digital Wiener filtering, Geophys. Prospect., 15, 311‑332. Rosa, A. and T. Ulrych (1991). Processing via spectral modeling, Geophysics, 56, 1244‑1251. Stockwell, R.G., L. Mansinha and R. Lowe (1996). Localization of the complex spectrum: the S transform, IEEE T. Signal Proces., 44, 998‑1001. Vardar, D. and B. Alpar (2016). High‑Resolution Seismic Characterization of Shallow Gas Accumulations in the Southern Shelf of Marmara Sea, Turkey, Acta Geophys., 64, 589‑609. Velis, D. R. (2007). Stochastic sparse‑spike deconvolution, Geophysics, 73, R1‑R9. Non-stationary seismic deconvolution 9 Wang, D.‑Y., J. P. Huang, X. Kong, Z. C. Li and J. Wang (2017). Improving the resolution of seismic traces based on the secondary time‑frequency spectrum, Appl. Geophys., 14, 236‑246. Wang, Y. (2002). A stable and efficient approach of inverse Q filtering, Geophysics, 67, 657‑663. Wang, Y. (2006). Inverse Q‑filter for seismic resolution enhancement, Geophysics, 71, V51‑V60. Xue, J., C. G. Cai, H. M. Gu and H. M. Luo (2021). Seismic Resolution Enhancement by Spectral Shaping using shaping‑ regularized inversion, IEEE Geosci. Remote S., 1‑1. Zhang, G., X. Wang and Z. He (2015). A stable and self‑adaptive approach for inverse Q‑filter, J. Appl. Geophys., 116, 236‑246. Zhang, P., Y. Dai, Y. Tan, H. Zhang and C. Wang (2019). A time‑varying wavelet extraction method using EMD and the relationship between wavelet amplitude and phase spectra, Chinese J. Geophys‑Ch., 62, 680‑696. Zhang, X., L. Han, F. Zhang and G. Shan (2007). An inverse Q‑filter algorithm based on stable wavefield continuation, Appl. Geophys., 4, 263‑270. Zhao, Y., N. Mao and J. Xu (2019). Generalized stable inverse Q filtering, J. Appl. Geophys., 169, 214‑225. Zhou, H., C. Wang, K. J. Marfurt, Y. Jiang and J. Bi (2016). Enhancing the resolution of non‑stationary seismic data using improved time‑frequency spectral modelling, Geophys. J. Int., 205, 203‑219. *CO R R E S PO N D I N G A U T H O R: Dengke H E, College of Geoscience and Surveying Engineering, China University of Mining and Technology (Beijing), Beijing 100083, China e‑mail: 201120@cumtb.edu.cn mailto:201120@cumtb.edu.cn