Multi-Ricker Spectral Modeling in the S-transform Domain for Enhancing Vertical Resolution of Seismic Reflection Data

Abstract A relatively straightforward methodology is presented for extending seismic bandwidth, and hence enhancing the seismic resolution by performing time-variant deconvolution. The generalized S-transform (GST) approach is used in order to properly compute the time-frequency components of the seismic reflection trace. In estimating the time-variant wavelet, a spectral modeling method is proposed, named multi-Ricker spectral approximation (MRA). After obtaining the estimated wavelet spectrum at each time sample, a deconvolution filter can then be built and applied in the S-transform domain. This proposed time-variant seismic enhancement method needs neither information on subsurface attenuation model nor assumption that the subsurface reflectivity is random. It is a data-driven methodology which is based on the seismic data only. This proposed method is validated on a noise-free and noisy synthetics and also applied to a field data. Results show that, after enhancement, overall seismic bandwidth can be extended resulting in a higher vertical resolution. Correlation with VSP corridor stack at well location ensures that the generated reflection details after enhancement are geologically plausible.


Introduction
Frequency bandwidth is one of the important factors in seismic reflection data analysis. A broader bandwidth makes seismic data easier to be interpreted and analyzed which, consequently, gives better result on the inferred subsurface geometry and its derived properties. Many attempts have therefore been made to broaden the seismic frequency bandwidth for resolution enhancement. The common method for enhancing seismic reso-lution is deconvolution, a process used to reverse the convolution effects on recorded seismic data. The model assumes that the seismic trace is a convolution between wavelet(s) and subsurface reflectivity series. The conventional seismic deconvolution methods considered that the propagating wavelet was a time-invariant wavelet which did not evolve with time, thus assuming the seismic trace as stationary signal. However, due to the attenuating nature of the earth subsurface, one should be aware that most seismic reflection data, in reality, consist of time-varying wavelets rather than a single wavelet. This is due to the fact that the longer the wave travels, the greater the attenuation affects the wavelet. Consequently, deeper reflection events will have a shorter bandwidth (broader wavelet), hence lower in vertical resolution compared to the shallower events. The seismic reflection trace should be considered better as a nonstationary signal.
During the last decade, a number of methodologies for enhancing seismic resolution have been proposed. Wang (2006) proposed a stabilized inverse-Q filtering approach in an attempt to broaden the seismic spectrum for resolution enhancement. The method is good when knowledge of the Q-values is available. However, the estimation of an accurate Q-model is difficult in practice. Smith et al. (2008) used Continuous Wavelet Transform (CWT) based methodology in extending the seismic bandwidth. The method commonly used a certain type of wavelets, such as Morlet or Mexican hat wavelets to match with the wavelets of the seismic data. Margrave et al. (2011) used the Gabor transform as the basis for a time-varying deconvolution method. Wang et al. (2013) proposed a spectrum-broadening method from a set of specific Gabor windows. Others used a combination between inverse-Q filtering and time-varying deconvolution approaches as discussed by van der Baan (2008van der Baan ( , 2012. Recently, Djeffal et al. (2016) showed that Margrave deconvolution (Margrave et al., 2011) can robustly be enhanced using the modified S-transform. While the original S-transform, as proposed by Stockwell et al. (1996) used a fixed mother wavelet. Zhou et al. (2016) and Jia et al. ( 2017) proposed a way of enhancing resolution in the generalized S-Transform domain following McFadden et al. (1999), Pinnegar and Mansinha (2003), and Gao et al. (2003).
In this paper a relatively straightforward methodology is presented for extending seismic bandwidth and hence enhancing the seismic resolution using the generalized S-transform. The approach allows us to properly compute the time-frequency components of the seismic trace. A spectral modeling method using Multi-Ricker spectral approximation (MRA) is then proposed in estimating the smooth and slowly varying wavelet spectra in the time-frequency domain. The proposed method needs neither the information on subsurface Q-model nor assumption that a certain type of wavelets exists in the data. It is a data driven methodology which is based on the input seismic data only. This proposed method is validated on a synthetic and applied to a field data to illustrate the performance and effectiveness in enhancing resolution.

Methods and Data
The methodology for seismic resolution enhancement proposed in this paper, shown by a simplified workflow in Figure 1, basically comprises three main steps. The first step is how to decompose nonstationary seismic data. The generalized S-transform is used in this step to accurately map the time-frequency distribution. The second step is how to estimate time-varying (or nonstationary) wavelets in the seismic data. In handling this task, a Multi-Ricker spectral modeling approach is used to fit time-variant wavelet spectra in the S-transform domain. The last step, after obtaining the time-varying wavelet, is how to design and apply the time-varying deconvolution operator to the nonstationary seismic data.

Estimating Time-Frequency Distribution of Nonstationary Seismic Data Using Generalized S-Transform
Even though the Fourier transform (FT) is a fundamental building block in seismic data analysis, the spectral information given by the FT represents the entire seismic trace. In order to show when (in time) a particular frequency component is dominant, one needs information in the time-frequency domain. The simplest method for generating time-frequency history is to use a short-windowed Fourier transform (Allen, 1977) which is often called as Short-Time Fourier Transform (STFT) proposed originally by Gabor (1946). In a continuous form STFT is defined as, Where: X(τ,f) is a two-dimensional time-frequency map with regard to τ f is as the time and frequency variables w(t-τ) is a time shifted rectangular window, x(t) is the input signal is the complex sinusoids in exponential form representing Fourier transformation from time to frequency domains.
In practice, a time signal is divided into shorter segments with equal length and the Fourier transform is computed for each of those segment. Plots of the changing spectra as a function of time, allow us to analyze nonstationary signals.
Another methodology for decomposing a signal into a time-frequency is the Stockwell transform (S-transform) which was originally proposed by Stockwell et al. (1996). The S-Transform is defined as, It looks quite similar to STFT, but instead of using a rectangular window as shown in (1), a Gaussian window with varying width is used. A narrow Gaussian window for handling high frequency events and broad for localizing the low frequency signals. In practical application, the generalized S-transform (GST) is, however, being used more often. The GST is defined as, This shows that the GST for m = 0 is equal to the Gabor transform. Similarly, the S-transform is obtained by setting γ = 1 and m = 1. Using the GST, the shape and width of the Gaussian window can flexibly be customized. To illustrate the effect of changing the parameters in the timefrequency analysis, a modulated sine wave is generated ( Figure 2a). Figures 2b -d show the corresponding time-frequency maps using the Gabor, S, and GST transforms. The temporal resolution in the Gabor transform (γ=0.2, m = 0) is poor in the low frequency range due to the fixed window. The temporal and spectral resolution of the S-transform (γ=1, m=1) is better in the lower frequencies, but loses resolution in the higher frequencies. The GST with γ=0.4 and m=0.7, shows better temporal and spectral resolution because of the use of a smaller temporal window. Parameters in the GST provide the ways to flexibly achieve the optimum time-frequency decomposition result.
To further illustrate the important of the S-transform in estimating the time-frequency distribution of nonstationary wavelets in seismic data, two synthetic signals are generated. These two signals, shown in the first column of Figures 3a -b, are formed by convolving a single reflectivity series with two nonstationary wavelets. The first wavelet has frequency content of 50 Hz to 40 Hz that decays linearly with time. The second wavelet has frequency content of 50 Hz to 15 Hz, but decays exponentially with time. Note on the estimated amplitude spectrum, as shown by the corresponding time-frequency (T-F) maps of those two signals that vary with time differently. It varies with time linearly on the first signal ( Figure  3a), but it varies with time exponentially on the second signal (Figure 3b). Considering the fact that both signals have exactly identical reflectivity series, the observed difference is clearly due Vol. 6 No. 3 December 2019: 223-233 to the amplitude spectrum of the nonstationary wavelets. This example demonstrates that amplitude spectrum of the time-varying wavelet may be viewed and estimated better in the timefrequency domain.

Estimation of Time-Variant Wavelet Using Multi-Ricker Spectral Modeling
On the conventional time-invariant deconvolution, a stationary wavelet is commonly estimated by averaging and smoothing the amplitude spectrum of the seismic traces. However, this is not the case for a nonstationary wavelet. As mentioned in the previous subsections, the nonstationary wavelet amplitude spectrum varies in shape with time. As such, the estimation of the wavelet amplitude spectrum should better be done in the time-frequency domain. Zhou et al. (2016) following Rosa and Ulrych (1991) approximated a smooth wavelet amplitude spectrum with a low order exponentiated polynomial up to seventh order. In this paper, inspired by the original Ricker's paper (Ricker, 1953) and the more recent discussions on Ricker wavelet by Gholamy and Kreinovich (2014), Multi-Ricker wavelet spectral modeling (MRA) is proposed in approximating the smooth wavelet amplitude spectrum at each time sample in the S-transform domain. In this scheme, the signal spectrum is simply described as a linear combination of Ricker wavelet spectra with different peak frequencies as shown in Figures 4 and 5. Each Ricker spectrum is defined as: Where: A(f) is Ricker wavelet spectrum f is frequency f m is peak frequency.
It can be seen that the unknown coefficients c is simply the least square solution shown in equation 5 below: Where: W is the dictionary of Ricker wavelet spectra λ is the damping factor constant S is the seismic signal spectrum at each time sample The damping factor is used to control the smoothness of the estimated wavelet spectrum. The advantage of using MRA is that it can efficiently be applied in the time-frequency domain  Figure 6.

Performing Time-Variant Deconvolution in the S-Transform Domain
After obtaining the estimated wavelet amplitude spectrum at each time sample, a deconvolution filter is then built in the S-transform domain. A deconvolution filter is designed using a simple spectral division scheme defined by the follow-  Apply Least-square Where: F is a deconvolution filter in the S-transform domain which is designed by simply performing division of the input signal I from the desired output spectra, DO.
The spectral division is performed on the amplitude spectra of the signal, not on their complex forms. This is done in order to keep the original phases of the input signals unaltered. However, most seismic data are of band limited natures. They have a limited frequency contents which commonly span from low to high. Outside this frequency bands, the amplitude spectrum goes down to very low values. A small value of prewhitening must be added to prevent division by zero. After the deconvolution filter in the Stransform domain being computed, it is simply multiplied with the amplitude spectrum of the input seismic signal to obtain the deconvolved signal. The enhanced seismic signal in time domain can then be computed from the deconvolved signal using the inverse S-transform.

Application to Synthetic Data
In order to see the effects of the proposed method in enhancing resolution of nonstationary seismic data, the nonstationary convolution model is utilized as described by Margrave et al. (2008) to create a synthetic seismic signal as shown in Figure 7. A synthetic seismogram (Figure 7a) is generated by convolving reflectivity series (Figure 7b) with a set of time-varying wavelets having frequency ranges from 50 Hz to 15 Hz that decays exponentially. The signal after applying the proposed enhancement method is shown in Figure 7c. Comparing the synthetic seismograms before and after enhancement process, it clearly demonstrates the ability of the proposed time-variant deconvolution method in enhancing the nonstationary seismic signal. An experiment with a noisy synthetic is also performed in this study. A noise-free synthetic seismogram (Figure 8a) is added by some additive noise of 10% ( Figure 8b) and of 20% (Figure 8c), to simulate two noisy synthetic seismograms. The signals after applying the enhancement method are shown in Figures 9a -9c. This illustrates how noisy signals can still be handled by the proposed method. This is due to the fact that the MRA tends to estimate the smooth wavelet amplitude spectrum as a red spectrum, i.e. a spectrum that skewed toward the low frequencies (as it is shown in Figures  4 and 5). This behaves as if a low-pass filter is applied to the input data before performing the deconvolution process, it will act like a preconditioning filter. As a result, the proposed method is able to suppress the additive noise and then enhance the resolution. Although the noise is still observed after deconvolution, its presence could still be tolerated.

Application to Field Data
In order to see the ability and stability of the proposed algorithm in enhancing nonstationary field seismic data, the proposed method is applied to a 2D seismic line, shown in Figure 10. In this example, there is a VSP corridor stack data that can be used to validate the resulting seismic data after enhancement. An example of T-F map from one trace shows that frequency content of the seismic data decreases with time ( Figure 11a). The estimated time-varying wavelet obtained using the MRA algorithm is shown in Figure 11b. The resulting seismic trace after enhancement in the time-frequency domain is shown in Figure 11c.
It can be seen that the seismic spectrum ( Figure 11c) is successfully extended as compared to the original T-F map shown in Figure  11a. The comparisons between the 2D seismic data before and after enhancements are shown in Figures 10 and 12. Note that the enhanced seismic ( Figure 12) shows a much better reflectors which are now well resolved compared to the original seismic data before enhancement ( Figure 10). Also note that the enhanced seismic is also properly correlated with the VSP cor-  ridor stack. This ensures that the result of new reflection detail generated after enhancement is geologically plausible.

Conclusions
Time-frequency decomposition using the generalized S-transform has been used for building a time-variant deconvolution algorithm for seismic resolution enhancement. Instead of using polynomial fitting in estimating time-varying seismic wavelets, a simple multi-Ricker spectral approximation is used to fit the wavelet spectrum at each time sample in the S-transform domain. The ability of the method in enhancing the seismic resolution has been shown as applied to synthetic and field data. Furthermore, in the application to field data, better correlation with VSP corridor stack at well location is achieved. This ensures that the result of new reflection detail generated after enhancement is geologically plausible and enables the data to be interpreted and further analyzed for stratigraphic details.