Non-Stationary Gabor Deconvolution Enhanced Full-Stack Penobscot Seismic Data is in Preparation for Public Download

The Penobscot dataset consisting of pre-stack seismic data and well-logs data is an exercising, open-to-public, oilfield data that is widely used by academics and researchers to test their workflow on reservoir evaluations, stratigraphic inferring, and reserve calculations. To make this data even more useful for such studies, we are enhancing the seismic data to an acceptable degree of resolution. Seismic resolution is primarily limited to wavelet overprints and interferences. Improving the vertical resolution of seismic data is challenged by the attenuation in wavelet propagation. A non-stationary deconvolution technique is needed to compensate for such anelastic effects. Gabor deconvolution is a widely accepted method to increase the resolution of seismic data as it utilizes information from the data seismic itself to minimize subjectivity in incorporating wavelet models. Our preliminary results show that the technique can widen up to twice the initial bandwidth. The resolution is thickened from 45.32 m to 22.15 m based on the tuning thickness analysis. The widening occurs on both sides at low and high frequencies without unnecessary energy increases and noise manifestations. The results show good agreement with the two available wells (L-30 and B-41). Our motivation is to provide higher resolution seismic data thus it can be used in workshops, software demonstrations, and academic discussion among practitioners and researchers. The preliminary enhanced seismic data is currently accessible by request, but it will soon be available to the public. In the future, we shall implement the same analysis to the full pre-stack seismic data so the benefits of Penobscot dataset will gain their optimum.


Introduction
In the current exploitation of oil and gas, exploration targets have increased from thick reservoirs to thin reservoirs which are difficult to resolve by seismic data. Changing reservoir targets requires new techniques to detect thin reservoir layers by increasing the resolution of the seismic data [1].
Seismic resolution is mainly limited to the overprint and wavelet interference associated with the energy from the source which will have an impact on the width or narrowing of the wavelet bandwidth [2]. The resolution of seismic data is influenced by the presence of subsurface anelastic effects which are very complex, making the response and energy inconsistent, so that the signal in the seismic data is non-stationary. The existence of nonstationary characteristics in seismic data, causing interpretation using seismic amplitude will give a problem to the reflector of seismic data which looks blurry so that it will give doubts in determining the stratigraphic boundaries and the results of reservoir mapping. One way to solve this problem is by performing a non-stationary deconvolution technique to compensate for the anelastic effect on the seismic data developed by Gary F. Margrave and Michael P. Lamoureux 2001 [3]. Gabor deconvolution is a widely accepted method for increasing the resolution of seismic data because it uses information from the seismic data itself to minimize subjectivity in aggregating and modeling wavelets that evolve over time assuming a minimum phase. In this paper, the Gabor deconvolution method will be applied to synthetic data and full-stack seismic data to obtain high-resolution seismic data.

Seismic
Seismic traces are obtained from the result of the convolution (*) between the wavelet with the reflection coefficient (RC) plus noise. In the form of an equation it can be written as follows: Where s(t) is a seismic trace, w(t) is Wavelet, r(t) is a reflectivity coefficient, and n(t) is a Noise.

Seismic Resolution
Resolution is the minimum distance of seismic waves to separate two adjacent objects, this is related to the occurrence of seismic wave interference. Seismic resolution is divided into two, namely vertical resolution and horizontal resolution, wherein this study more emphasis is placed on a vertical resolution.
Vertical resolution is the minimum distance of seismic waves to separate two adjacent objects, the ability to separate the two objects or rock reflectors depends on the thickness and wavelength. The minimum thickness that can still be distinguished is called the tuning thickness. The minimum thickness of an object to be able to provide its reflection varies from 1 8 ⁄ to 1 30 ⁄ , it's just that because the noise in the seismic data makes the thickness of the tuning that can be separated is 1 4 ⁄ [2], where the seismic wavelength is can be formulated as follows: Where v is an average seismic velocity in well, and f is the dominant frequency from the area of interest.

Spectral Decomposition
Spectral decomposition is a seismic signal processing method based on time-frequency spectrum analysis, this is because the seismic data is non-stationary or the frequency changes with time due to the absorption of waves by the medium being passed. In carrying out time-frequency analysis on seismic data, it is necessary to transform the seismic data from the time domain to the frequency domain by performing Fast Fourier Transform (FFT) [4]. This spectral decomposition is carried out because a reflection from a thin layer has certain characteristics in the frequency domain that these characteristics cannot usually be seen in the time domain [4].

Gabor Transform
This Gabor transform adopts the principle of Short-Time Fourier Transform (STFT) by using a window with a constant width that is shifted with time. The Gabor transform is a Fourier transform of a seismic signal that has been carried out by spectral decomposition using a Gaussian window which makes the seismic signal localized. Mathematically, Gabor transformation can be expressed as follows (Mertins, 1999) Where s (t) is the signal and g (t) is the Gabor analysis window where τ is the window center location. Signals that have been transformed into the frequency domain can be returned to the time domain by inverse Gabor transformation which can be formulated as follows:

Gabor Deconvolution
Starting from the non-stationary nature of seismic data that always changes with time, with the assumption that the earth's subsurface is not a perfectly elastic medium and seismic waves weaken as they propagate below the surface. This attenuation of seismic waves occurs because the absorption of waves by each medium in its path is marked by a changing wavelet shape. The simplest model of attenuation is the constant Q model from Kjartansson (1979) [5], where Q-constant refers to Q which does not depend on frequency but depends on space and time traveled, Margrave (2001) uses Q-constant as an attenuation [3] function which can be formulated as follows: Where Q is the quality factor of the attenuation medium, and H is the Hilbert transform. The Gabor deconvolution is an extension of the Wiener deconvolution process carried out in the frequency domain on the data that has been carried out by the Gabor transformation [3]. This Gabor deconvolution compensates for the attenuation factor in the seismic data caused by the earth filter effect, which can be formulated as follows: Where ( , ) is Gabor transformed signal, ( ) is the source wavelet is in the frequency domain, ( , ) is the attenuation factor of the signal, ( , ) is Gabor deconvolution reflectivity estimation, | ( , )| is the result of smoothing the signal transformed by Gabor.
Smoothing is done for estimation of wavelet propagation, can be done by a convolution between the results of the Gabor transformation with 2D boxcar, Gaussians, triangles, and hyperbolic smoothing [3]. Thus, the amplitude spectrum of the deconvolution operator becomes: Where is the maximum value of | ( , )| and is stability constant with small dimensions for valid logarithmic computations. The deconvolution operator is assumed to have the minimum phase in wavelet propagation calculated from the amplitude spectrum and the Hilbert transform aid so that it can be formulated as The results obtained in the Gabor deconvolution are in the form of the Gabor reflection coefficient which is still in the frequency domain, to return it to the time domain it is necessary to inversion the Gabor transformation using equation 4.

Smoothing Hyperbolic
Hyperbolic smoothing refers to equation 5 by assuming the Q-constant model which represents the attenuation factor in the form of constant τf, by dividing the time-frequency map into hyperbolic lines and calculating the average of each hyperbolic line representing the amplitude lost due to attenuation. Broadly speaking, the attenuation effect is obtained by filling each hyperbolic line with the mean of the interpolated time-frequency map [3,9]. Can be formulated as follows: Where ( , ) is an indicator of the function of ℎ hyperbolic lines for example will be zero except at the point of the hyperbolic line ℎ . N is the total number of hyperbolic lines, so an estimate is obtained from | ( , )| . To estimate | ( )| can be approached by removing the estimated attenuation and running a little smoothing towards the frequency, this is following convolution smoothing ( ) [3,9].

Geological Regional and Methodology
The location of this research is in Canada, precisely offshore Nova Scotia in the Penobscot fields, with the coordinates 44˚ 07' 46'' N / 60˚ 06' 00'' W. In this area, a survey and acquisition were carried out in 1992. As can be seen in Figure 1 above Penobscot field is at the top of the Scotia Basin, in this field half of the basin is on the continental shelf and part is on the continental slope. The Scotia Basin was formed in the passive continental margin that developed after the rifting period when North America separated from Africa during the Pangea rupture period. Can be seen in Figure 2 above that the stratigraphy of the Scotia Basin is dominated by clastic sedimentary deposits, namely sandstone, and shale and there are carbonate rock inserts.
In this research, non-stationary deconvolution was carried out using the Gabor deconvolution algorithm which was made using python software so that it was necessary to test non-stationary synthetic data. After doing some analysis on the synthetic data starting from the Gaussian window parameter in carrying out the Gabor transformation, smoothing hyperbolic, inverse Gabor transformation and the results of the deconvolution Gabor estimation are similar to synthetic reflectivity and a wider spectrum almost synthetic reflectivity. Then this algorithm is applied to the full-stack data to obtain high-resolution data, but it is necessary to do an analysis first in calculating the wavelet propagation estimate from the trace in the well which will be applied to the entire seismic volume when the Gabor deconvolution is carried out. After that, overlay the reflectivity in the well with the estimated reflectivity of the Gabor deconvolution. Cross plot and tuning thickness tests was carried out and validation of each spectral range of the well. The above process can be described as the workflow in Figure 3 below.  Figure 3. Workflow of the research.
The first step in this study was to test the Gabor deconvolution algorithm, to conduct trials using non-stationary synthetic data with a minimum phase, the to create a random reflection coefficient (1D) with a maximum time of 6 s and a sampling rate of 0.004 s. The wavelet used in this process is the ricker wavelet that has been phase transformed using Hilbert transformation process into a minimum phase wavelet with a dominant frequency of 30 Hz, and a wavelet length of 0.6 s. After that, the signal is carried out for each seismic trace that has been localized. This process is called the Gabor transformation as shown in figure 4a. The next step is carried out hyperbolic smoothing on the data that has been Gabor transformation. This hyperbolic smoothing aims to estimate wavelet propagation that has added an attenuation factor to the synthetic data itself as shown in figure 4b. Workflow application on synthetic trace can be seen in Figure 4 Figure 6. a). The amplitude spectrum results from non-stationary synthetic data b). Amplitude spectrum after Gabor deconvolution c). Crossplot reflectivity true vs Estimation reflectivity Gabor deconvolution.
The analysis of the estimated reflection coefficient in Figure 5.4 (c) is carried out. When compared with the reflection coefficient, it actually has a very strong correlation between the two as shown in the crossplot above. This will also be reflected in the amplitude spectrum of the data that has a bandwidth widening that is almost close to the amplitude spectrum of the synthetic reflection coefficient. By designing suitable parameters, the effect of wavelets on seismic data can be maximally eliminated. It can be said that the Gabor deconvolution process is very effective when selecting parameters that are in accordance with the seismic data, so that this method can answer the limitations of seismic resolution to perform a more detailed subsurface stratigraphic analysis. As well as Gabor deconvolution is very sensitive to the presence of noise, so it is necessary to select suitable parameters so that the presence of noise in the data can be eliminated and lead us to a more precise estimation of the reflection coefficient.

Gabor Deconvolution Analysis of Full-Stack Data
After a series of tests for Gabor deconvolution parameters on synthetic seismics, Gabor deconvolution algorithm is then applied to real seismic data. The data used are 2D arbitrary line seismic data controlled by 2 wells B-41 and L-30. In carrying out the deconvolution process, Gabor assumes that in the seismic data there is wavelet propagation or it can be said that there is still an attenuation effect on the seismic data. Furthermore, the search for the right parameters in real 1 trace data analysis near the well. In order for Gabor deconvolution process on real data to run optimally to obtain a wider seismic data bandwidth, the validation process uses well data to confirm the results of Gabor deconvolution. After analyzing the parameters suitable for real trace data with the maximum Gabor deconvolution results, it is then applied to 2D seismic data using the same parameters for each full-stack seismic trace. In the figure below is a comparison of the before and after real seismic data applied by the Gabor deconvolution algorithm.  Figure 6. (a) 2D arbitrary line original full-stack data and amplitude spectrum with well B-41, L-30 (black trace is synthetic seismic & red trace is log gamma ray), (b) 2D arbitrary line and amplitude spectrum after Gabor deconvolution. Figure 6 above shows the occurrence of seismic data bandwidth widening and the recovery of high frequency lost due to attenuation, because in Gabor deconvolution the effect of wavelet propagation on the data has been removed. The widening of the bandwidth spectrum in seismic data will be reflected in the increase in seismic temporal resolution so that the ability to temporally separate layers will increase.
After checking the correlation between the estimated reflection coefficient after Gabor deconvolution with wells B-41 and L-30, the correlation results were 0.136 and 0.1111. With a very low correlation value, it will not completely result in a wrong estimate because the processes and stages carried out are very careful in selecting the most optimal parameters. It's just that in the real data, there is a limitation in getting the estimation of the subsurface reflection coefficient, the results are not like synthetic data. Based on research, in real seismic conditions, many factors are difficult to eliminate, including the estimated wavelet propagation which is assumed to be the minimum phase when viewed in the field, the minimum assumption is not always constant but varies depending on subsurface complexity, the estimation of the attenuation effect with time-constant frequency along the hyperbolic strip is not well suited for real seismic data, and noise exposure [11].   Figure 7. (a-d) plotting color jet before dan after Gabor deconvolution 2D arbitrary line on full-stack data, (e-h) plotting color grays density before dan after Gabor deconvolution (Black log is the zero-phase synthetic log, and red log is the gamma ray log).
From the image above, it can be seen that the details of the different layers with different plotting colors can be seen in Figure 7 showing the details of the layers more clearly laterally and can be validated with a low gamma ray log representing sand lithology and synthetic seismic with trough events. As in the picture above, the deconvolution results can resolve the non-resolved thin layers undergoing tuning by using color grays density plotting. From the results of plotting using gray density, it can show low frequency in seismic data which is characterized by a more detailed termination and clearer contrast between layer boundaries which can reduce ambiguity in determining layer boundaries, making it possible to perform a more detailed and precise stratigraphy analysis. It can be seen that in TWT 2 -2.1 s the contrast of the layers is very clear for the boundaries of the formation, if it is continued to the well then this area shows a thin sand layer that is characterized by the low log gamma ray and low synthetic seismic.

Tuning Thickness Analysis
From Figure 8 below, the tuning thickness can be calculated from the seismic data before and after Gabor deconvolution, the dominant frequency of the seismic data before and after is 22 Hz, 45 Hz with an average speed of 3988.75 m/s. From the calculation results, it is found that the tuning thickness that can be seismic resolution is 45.32 m or 11.3 ms for data before Gabor deconvolution and 22.15 m or 5.5 ms for data after Gabor deconvolution. It can be concluded that by increasing the width of the seismic data bandwidth, the tuning thickness that can be resolved by the seismic will thinner so that the separation of each layer will increase temporal.

Analysis of The Spectral Range
This analysis is useful to prove whether Gabor deconvolution results can recover all reflection signal information both low and high frequency from seismic data before and after Gabor deconvolution which is validated with synthetic seismic data and gamma ray logs that have been checkshot corrections on data wells then bandpass filter with the same frequency as the seismic data. For data from Gabor deconvolution below, it is necessary to check the respective signal spectrum, low frequency, enhancement result, and high frequency to ensure whether the information seen is a reflector or a result of processing so that it needs to be validated with a checkshot corrected well. In the figure below, it can be seen that the results of the bandpass filter on seismic data with a higher frequency in the frequency domain will be reflected in the time domain with narrower temporal separation so that the reflector is seen in the dominant image is very thin. From the results of the bandpass filter seismic data and synthetic seismic data in the well, it can be seen that in almost every spectral range there is validated signal information with well-logs data, while at high frequencies there is no reflector signal information, except for TWT 0.25 and 1 s. From the analysis of the low and high frequencies above, it can be interpreted that Gabor deconvolution has succeeded in optimally recovering all the information on the reflector signal as well as it is in the frequency range of 5-85 Hz.   Figure 9. (a) Estimation result of 2D arbitrary line full-stack data after Gabor deconvolution at 1-15 Hz spectral range (black trace is a synthetic seismic estimate of well B-41 (CDP 108) and well L-30 CDP (390) whose wavelets are statistically extracted from seismic data with zero phases while red traces are log gamma ray), (b) Original in the spectral range 1-40 Hz, (c) After Gabor deconvolution in the spectral range 40-85 Hz, (d) After Gabor deconvolution in the spectral range 105-124 Hz.

Conclusion
Analysis of non-stationary deconvolution on synthetic and data real using the Gabor technique in windowing and estimating propagation wavelet succeeded in providing increased bandwidth at high and low frequencies up to 5-85 Hz (2x the bandwidth initial). Gabor deconvolution non-stationary showed a significant improvement in resolution and was validated using well log information available across all spectral ranges. Quantitatively there was an increase in the resolution temporal seismic with a tuning thickness of 45.32 m or 11.3 ms to 22.15 or 5.5 ms. With the Gabor technique, limitations of bandwidth data can be overcome, making it easier to interpret stratigraphically in a clear and detailed manner. Gabor deconvolution does not completely eliminate the attenuation effect, but this method considers physical parameters which are calculated directly from the seismic data itself.

Acknowledgments
Thanks to the geophysical engineering lecturers who have shared their advice and experiences. The geophysical engineering laboratory of the Sumatera Institute of Technology has allowed us to use the lab facilities so that research can be realized.