Introduction
Seismic signal analysis plays a key role in petroleum exploration by helping to enhance information difficult to visualize naked eye. Today several techniques and algorithms are used to interpret seismic 3-D data. Spectral analysis comprises several methods that enhance specific seismic information enabling to solve stratigraphic and structural details (i.e., Rivera-Recillas et al., 2005 and Coconi-Morales et al., 2010), to estimate reservoir dimensions, etc. Fourier transform is commonly used to analyze frequency content of a seismic signal. However, when frequency content varies with time, this tool cannot show time position of the frequency content. Spectral content of seismograms varies significantly with time, i.e., they are non-stationary and require non-standard decomposition methods. The discrete wavelet transform enables decomposing a non-stationary time series in its different frequency components and transforms the time domain information into a time-scale domain where scale is inversely proportional to frequency.
Discrete wavelet transform (DWT) is based in filter bank theory. Convolution of a filter bank with a signal provides frequency rank windows. The filter banks for a particular wavelet must satisfy two conditions; must be of compact support and of zero average. This tool enables separating high frequencies from low frequencies and locating its position in time. A growing number of geophysical studies using DWT have provided satisfactory results; however sometime, separation between scales is not fully achieved due to an intra-scale coupling effect. So that, an appropriate wavelet is needed to conduct a successful data processing based in this technique. There are many wavelets but many present the aliasing problem, here we present a performance analysis of several wavelets with respect to the aliasing effect in multiresolution analysis of seismic signals. We tested the performance of Vaidyanathan wavelet with a real 3-D seismic cube data.
We first briefly introduce the wavelet theory. Then we describe the Vaidyanathan wavelet. We reconstructed seismic sections based in the respective multiresolution analysis of the Boonsville field seismic cube, located in north-central Texas (Bureau of Economic Geology, 1996). The seismic cube is open access and well documented.
Theory
A signal can be expressed in terms of a set of functions with different resolution. Multiresolution based on the discrete wavelet transform generates this function base, to each resolution a certain information content of the signal is associated. The theory of discrete wavelet transform has been exposed in many books (i.e., Hubbard, 1998; Chui, 1992; Dwight and Olejniczak, 2003). This technique can be developed on Daubechies's (1992) pyramidal algorithm where the discrete wavelet transform is obtained by convolving the signal with a quadrature mirror filter (QMF) bank, built from a compact support and zero mean wavelet. The wavelet is dilated to different scales by a factor of two. Translation is done in binary form.
The dilatation and shifting of a wavelet ψ (x), can be expressed as:
where j denote scale and к translation. At small scales, when the wavelet is contracted, high frequencies are displayed; at great scales, when the wavelet is expanded, low frequency contents are obtained. From expression (1) two functions are generated which are employed in the decomposition: a wavelet function as expression (1) and a scalar function as:
The functions ψ(x) j,k and φ(x) j,k generates a sub-set of a vector space that spans signals orthonormal to the analyzed signal.
The respective filter bank is constituted by one high and one low pass filter. The low pass filter is obtained from the wavelet function (ψ), while the high pass filter is obtained from the corresponding scalar function (φ).
So, when the high pass filter is applied to a non-stationary signal detailed coefficients are obtained, and when the low pass filter is applied to the same signal, we obtain smooth or approximate coefficients. The detail coefficients capture the top half frequency content of the data while the smooth coefficients contain the bottom half frequency content. This first step corresponds to the first decomposition level and is named first scale. To generate the next scale the smooth coefficients are used as input signal, and the above described process is repeated (Figure 1). To decompose the signal, a sub-sampling is done, because the translation of a wavelet along the signal is made in a binary form.
The wavelet coefficients can be inverse transformed to exactly reproduce the original time series. This is achieved by using the filter bank in its synthesis form and reversing the sequence of the forward transform algorithm (Figure 2). Because the sub-sampling introduced in the forward decomposition, in reconstructing the signal it is needed an up-sampling by 2, and this can be achieved by adding one zero between two coefficients.
Vaidyanathan wavelet
As mentioned, discrete wavelet transform uses a wavelet to build an analysis based in frequency content, and such a wavelet can be derived from a pair of filters which satisfies the following frequency domain conditions for a perfect reconstruction (Foster et al., 1997):
where L(ω) and
these are known in the literature as a Quadrature Mirror Filter (QMF). In wavelet applications Finite Impulse Response (FIR) filters satisfying conditions (3) and (4), are compact support in the time domain, which is important for space-time operations.
There exists a group of wavelets that satisfies these conditions. The choice of a wavelet is very important for any wavelet domain processing application. In seismic processing it is desirable a wavelet that produces an optimal separation of information between scales and gives rise to a minimum overlap. Thus we need a wavelet which enables a perfect reconstruction and will minimize any artifact that may be introduced in the processing of a signal and appearing in its reconstruction.
There exists a two-channel QMF bank which satisfies the condition for a perfect reconstruction, and ensures a good stop-band of frequencies. Vaidyanathan and Hoang. (1988) introduced this filter. This wavelet is known as Vaidyanathan wavelet (Figure 3). We implemented the multiresolution analysis and conducted a performance test of several wavelets including the Vaidyanathan wavelet.
Programming
The programmed multiresolution analysis was based on the pyramidal algorithm, and on the 1D discrete wavelet transform. The programmed structure comprises five sub-programs to allow an optimal execution (Figure 4).
The objective of the first program is to communicate with the user (i.e., a friendly interface). Several windows enable the user to input the data file, to select a wavelet, as well as the scale or resolution level to be displayed.
The second sub-program distributes this information to other three sub-programs. The third sub program allocates enough space for all of the needed variables. The fourth subprogram contains all information defining each of the wavelets contained in the catalog shown in the first subprogram. Finally, the fifth subprogram performs the wavelet transform by means of the pyramidal algorithm.
This last subprogram analyses the seismic information trace by trace. It handles in this way a 3D data volume. When the user wants to elaborate a time-slice, this program shows the values at the user selected time. The case of one horizon is managed similarly.
Performance assessment
To assess how well isolated is the frequency content associated with a given scale was the goal of this study (i.e., which wavelet does preserve the power spectrum in an optimum way). Figure 5 shows a signal created by summing a series of sines with frequencies between 30 and 211 Hz. This signal encompass a wide frequency range: it has a 511 msec length and a sampling rate of 1 msec. It can be considered a non-stationary signal. This is not a seismic signal but is useful to assess the performance, in multiresolution analysis, of the following wavelets: Haar, Symplet, Coiflet, Daubechies (wavelets employed in Matlab, Misiti et al., 1996) and Vaidyanathan (Vaidyanathan and Hoang, 1988).
Figure 6 shows the power spectrums of the original signal and those corresponding to the second level (or scale) of the multiresolution analysis based on the above mentioned wavelets.
To the second decomposition scale approximately corresponds the frequency content between 125 and 250 Hz of the original signal. The original signal (blue line) power spectrum has a stronger gradient towards approximately 170 Hz. In this portion, the reconstructed spectrums from the multiresolution analysis based on Haar, Symplet, Coiflet, Daubechies wavelets show spurious effects. The exception corresponds to the Vaidyanathan wavelet. At 350 Hz approximately, Haar, Symplet, Coiflet, Daubechies wavelets produce a spurious pike, possible an armonic from the information contained in the frequency range of the second scale. We can see that for the Vaidyanathan wavelet this effect is minimum.
This performance analysis indicates that the Vaidyanathan wavelet best preserves the original signal power spectrum (i.e., it distorts in a minimum degree the spectrum of the original signal), so that, in particular, this wavelet is very well suited to conduct the multiresolution analysis (trace by trace in this study) of 3D seismic data, where it is very important to preserve the original seismic amplitude and no to introduce artifacts.
Example of multiresolution of 3D seismic data
Boonsville 3D seismic data
Boonsville 3D seismic dataset were obtained in the Jack and Wise counties, Fort Worth Basin, north-central Texas (Figure 7). The study area comprises approximately 67 km2. The data are well documented and can be acquired from the Bureau of Economic Geology. Vaidyanathan and Hoang, (1988) wavelet based multiresolution analysis of the Boonsville 3-D seismic cube is now presented.
The data length is two seconds, sampling of the seismic data were done at 1 ms (Bureau of Economic Geology, 1996), with a 500 Hz corresponding Nyquist frequency. The dominat frequency is 57 Hz, with mean velocity of 3,600 m/s. According to Rayleigh criterium, the corresponding vertical resolution only allows to study layers with a thickness > 15 m.
Figure 8 shows the Boonsville seismic 3-D cube layout and the location of the oblique section used to illustrate the performance of the multiresolution analysis. This section was created to pass through two wells (BY11 and BY13 wells), because we wanted a zoom along of the well length to see information at detail. Figure 9 shows the seismic events along this oblique section, in the right lower part it is showed the respective power spectrum. We can shown that approximatly the dominant frequencies range from 30 to 115 Hz. We use one scale within this dominant frequency range to conduct an multiresolution analysis using Vaidyanathan wavelet and the Daubechies wavelet of order 10 (see Figure 6).
Figure 10 shows the reconstructed section using only the fourth scale (frequency content between 31.25 and 62.5 Hz). The respective multiresolution analysis was based on Vaidyanathan wavelet. Figure 11 presents the corresponding fourth scale obtained from the Daubechies wavelet of order 10 based multiresolution analysis.
The power spectra of the reconstructed sections based respectively on the Vaidyadanath and Daubechies wavelets (Figures 10 and 11) indicate that the aliasing effect due to the Vaidyanathan wavelet is not visible (i.e., negligible). However, the Daubechies wavelet of order 10 introduces high frequencies armonics enclosed by black elipsoids (i.e., a noticeable effect).
We can note that the seismic horizons in Figure 10 correlate quite well with the original seismic information. The seismic horizons in Figure 11 correlate fair well with then original seismic horizons. However, seismic horizonts in Figure 10 change its position along the seismic section in a smooth manner, but the seismic horizons in Figure 11 presents jumps. We believe that it is due to high frequency associated with the high frequency armonics introduced by the Daubechies wavelet of order 10 (see the power spectrum of Figure 11).
Finally, Figure 12 shows the FFT based bandpass filtered section of the original data. The band pass filter encompasses the fourth scale frequency range (31.25 - 62.5 Hz). Figures 10 and 12 correlate very well, which illustrates how the Vaidyadanath wavelet minimizes the aliasing effect.
Conclusions
This study reports a performance analysis of the Vaidyanathan wavelet that minimizes the aliasing effect (Figures 5 and 6) (Vaidyanathan and Hoang, 1988). Several wavelets have been tested with a synthetic signal and with a real 3-D seismic dataset. The best results were obtained with Vaidyanathan wavelet.
This study has illustrated how a discrete wavelet transform based multiresolution analysis makes possible separation of the information content of a non-stationary signal in different frequency ranges. This descomposition provides the seismic interpreter frequency information of interest that might not be visible in band-pass filters. This numerical study shows that several wavelets can be used with this technique, but it is important to select the appropriate wavelet, because a bad selection can give rise to spurious effects (i.e., artifacts) due to the overlay between scales, causing that the amplitude of the frequency content of some frequencies be enhanced, and the seismic interpreter can be mislead with these artifacts, and consider them subsoil information.