An accurate estimation of wavelet is crucial in the deconvolution of seismic data. As per the convolution model, the recorded seismic trace is the result of convolution of the earth’s unknown reflectivity series with the propagating seismic source wavelet along with the additive noise. The deconvolution of the source wavelet from the recorded seismic traces provides useful estimates of the earth’s unknown reflectivity and comes in handy as an aid to geological interpretation. This deconvolution process usually involves estimation of a wavelet, before it is removed by digital filtering. Since the earth’s reflectivity and seismic noise are both unknown, the wavelet estimation process is not easy. The statistical methods estimate the wavelet using the statistical properties of the seismic data and are based on certain mathematical assumptions. The most commonly used method assumes that the wavelet is minimum phase and that the amplitude spectrum and the autocorrelation of the wavelet is the same as the amplitude spectrum and the autocorrelation of the seismic traces, within a scale factor, in the time zone from where the wavelet is extracted. With the assumption that the wavelet is minimum phase, an estimation of the wavelet is done from the trace autocorrelation. This method always estimates a minimum phase wavelet and so is suitable for wavelet estimates from seismic data acquired using explosive sources, only if their source signature is purely minimum phase and that is retained through during processing. However, it is not applicable for estimating wavelets from sources giving mixed phase signature. For example, deconvolution of non-minimum phase seismic data with a minimum phase wavelet will leave behind a spurious phase in the data. The reason for this is that the autocorrelation function and the power spectrum mentioned above are second-order statistical measures. They contain no phase information and so cannot identify non-minimum phase signals. Also, these measures work well for Gaussian probability distribution of amplitudes, and so will not yield accurate results for non-Gaussian or non-linear distributions. The non-Gaussianity in the seismic data could arise from the non-minimum phase source signature, noise in the data like swell noise, and a non-linear earth response. Consequently, higher-order statistics have been used for dealing with non- Gaussian distributions. These statistics, known as cumulants, and their associated Fourier transforms known as polyspectra, not only reveal amplitude information, but also the phase information (Mendel, 1991). In this paper, we address this issue through the use of higher-order statistics such that the phase component in the data are more accurately estimated and removed.
Higher-order statistics for wavelet estimation
Cumulant, a higher-order statistical property, preserves the phase information of the wavelet under the assumption that the reflectivity series is a non-Gaussian, stationary and statistically independent random process. The second-order cumulant of a zero-mean process is just the autocorrelation, which as stated above, has no phase information. The thirdorder cumulant is a two dimensional correlation function. For a Gaussian process, all cumulants above the second-order are zero, but are non-zero for non-Gaussian processes. Thus these two statistics are not suitable for recovering a non-minimum phase from a convolution process such as the seismic trace. The fourth-order cumulant is a three-dimensional correlation function which contains information about phase. Just as the Fourier transform of the autocorrelation function yields the power spectrum, similarly, the trispectrum relates to the fourth-order cumulant via the 3D Fourier transformation. Lazear (1993) and Velis and Ulrych (1995) estimate the phase of a wavelet by fourth-order cumulant matching wherein an initial guess for the wavelet is iteratively updated until its fourth order statistics match those of the seismic data.
Misra and Sacchi (2006) suggest the parameterization of the embedded mixed phase wavelet as a convolution of the minimum phase wavelet with an all-pass operator. The allpass operator can further be parameterized as the ratio of a maximum phase time sequence and corresponding minimum phase time sequence with the necessary time delay required to enforce causality in the all-pass operator (Porsani and Ursin, 1998). The denominator term in the paramerization of the allpass operator is a minimum phase sequence whose length and coefficients are unknown. As discussed in a later section, we optimize for the unknown coefficients of the minimum phase sequence keeping the length as a constant parameter.
Seismic data is represented as the convolution of the reflectivity sequence with the unknown wavelet. The unknown wavelet showing mixed phase characteristics is further represented as the convolution of minimum phase wavelet and the all-pass operator. Thus the seismic data is further represented in terms of two convolutions, namely convolution of reflectivity with the minimum phase wavelet which in turn is convolved with the all-pass operator. Deconvolving the data with the minimum phase wavelet increases the bandwidth of the output data and we subsequently refer to it as the whitened data. The whitened data is thus represented as the convolution of the underlying reflectivity series and the allpass operator. Hence it is possible to make an estimation of the underlying reflectivity series by estimating the all-pass operator from the whitened data.
Development of the algorithm
The well known Barlett-Brillinger-Rosenblatt formula (Lazear 1993; Mendel 1991) links the fourth-order cumulant of the seismic trace with the fourth-order moment of the embedded wavelet. For non-Gaussian, statistically independent and identically distributed reflectivity series, the fourth-order cumulant of the seismic trace is equal to, within a scale factor, the fourth-order moment of the wavelet provided that the noise distribution is Gaussian. The optimization procedure described in the following paragraph minimizes the cost function given by the L2-norm between the fourth-order normalized trace cumulant and the fourth-order wavelet moment.
The optimization of the cost function thus involves computation of the normalized 4th order trace cumulant of the whitened data (data obtained after deconvolution with a minimum phase wavelet) and the normalized 4th order moment of the all-pass operator.
The shape of the cost function is unknown and may contain several local minima. Local optimization methods based on gradient computation always proceed to the minimum nearest to the chosen initial model. In the present optimization problem where the shape of the cost function is not known, a global optimization algorithm is a preferred choice. Simulated annealing algorithm with a Metropolis acceptance/ rejection criterion (Misra, 2008) is adopted for the optimization of the cost function. The model parameters for the simulated annealing optimization are the coefficients of the minimum phase sequence in the parameterization of the all-pass operator. The all-pass operator for each of the generated model is computed from equation by taking the ratio of the corresponding maximum phase sequence and the minimum phase sequence.
Application to post stack seismic data
In order to test the stability and reliability of the algorithm, the method outlined above was applied to a seismic data volume from North America. For this volume, the data are further subdivided into smaller zones and for each of these zones a mixedphase wavelet is estimated. The data in each of these subdivisions are deconvolved using the estimated mixed-phase wavelet. Further, an average estimated mixed-phase wavelet is obtained by taking the mean of the individual mixed-phase wavelets. The mean mixed-phase wavelet is then used to deconvolve the data. The results are correlated with the Pwave log curve.
Figure 1 shows the workflow for the method outlined above.
A seismic data volume from a province in North America was picked up for testing the algorithm outlined above. The data volume had been processed using a conventional processing sequence. Four different locations on the post-stack volume were selected for the estimation of the mixed-phase wavelets in the same time interval (500ms) for each.
Figure 2a shows a segment of a seismic section around location 1. Figures 2b, c and d show the estimated minimum-phase wavelet, the estimated all-pass operator and the estimated mixed-phase wavelet at the location 1. The minimum-phase wavelets are estimated from the average autocorrelation of the seismic traces by the Wiener- Levinson algorithm. The estimated minimum-phase wavelet is deconvolved from the data which resulted in broadening of the bandwidth. Figures 2e shows the phase-corrected data for the location 1. Figures 3, 4 and 5 show a similar set of images for locations 2, 3 and 4. Notice that for each set of images, the mixedphase wavelet deconvolved sections exhibit the highest level of detail. Again, it would be advisable to correlate the deconvolved sections with the P-wave log curves to gain some confidence in ascertaining if the resolved reflections correlate well and if they are authentic. Notice in Figure 6, the correlation of the section in Figure 6b is better than the one shown in Figure 6c.
Deconvolution of seismic data with a minimum phase wavelet effectively removes the amplitude spectrum of the wavelet from the data. However, in situations where the minimum-phase assumption about the wavelet is not valid, the deconvolution leaves behind a spurious phase component in the data. The method adopted in estimating and hence removing the spurious phase involves estimation of the coefficients of an all-pass operator from the data that have been whitened by the deconvolution of the minimum-phase wavelet. The whitened data are used to optimize the cost function involving the 4th order normalized trace cumulant and the 4th order moment of the all-pass operators. The optimization procedure uses simulated annealing with the Metropolis acceptance/rejection criterion. The estimated all-pass operator is subsequently convolved with the earlier estimated minimum-phase wavelet to estimate the mixed-phase wavelet in the data. The suggested method is tested on a seismic data set belonging to a province in North America. The data set is subsequently deconvolved with the estimated mixed-phase wavelets. Further, an average mixed-phase wavelet is computed from the individual mixedphase wavelets. The data are then deconvolved with the average mixed-phase wavelet are correlated with the P-wave log data and shows better correlation.
We acknowledge Prof. Mauricio Sacchi for his help and suggestions during the development of the algorithm. We also acknowledge an anonymous provider of the data shown in and Arcis Corporation, Calgary, Canada for support of this work.
About the Author(s)
Somanath Misra was born and raised in the coastal province of Orissa, India. After graduating in Physics major, Somanath joined the Indian School of Mines for his Master’s degree in Applied Geophysics, which he successfully completed in 1991. Thereafter, Somanath preferred to join the exploration department in Orissa, India and worked there for about 10 years. In 2003, he joined the Signal Analysis and Imaging Group at the University of Alberta and completed his Ph.D. in 2008. Since then he is working as a Reservoir Geophysicist in the reservoir services group at Arcis Corporation, Calgary, Canada. He is a member of CSEG and SEG.
Satinder Chopra received M.Sc. and M.Phil. degrees in physics from Himachal Pradesh University, Shimla, India. He joined the Oil and Natural Gas Corporation Limited (ONGC) of India in 1984 and served there till 1997. In 1998 he joined CTC Pulsonic at Calgary, which later became Scott Pickford and Core Laboratories Reservoir Technologies. Currently, he is working as Chief Geophysicist (Reservoir), at Arcis Corporation, Calgary. In the last 26 years Satinder has worked in regular seismic processing and interactive interpretation, but has spent more time in special processing of seismic data involving seismic attributes including coherence, curvature and texture attributes, seismic inversion, AVO, VSP processing and frequency enhancement of seismic data. His research interests focus on techniques that are aimed at characterization of reservoirs. He has published 7 books and more than 200 papers and abstracts and likes to make presentations at any beckoning opportunity. He is the Chief Editor of the CSEG RECORDER, the past member of the SEG ‘The Leading Edge’ Editorial Board, and the Ex-Chairman of the SEG Publications Committee.
He received several awards at ONGC, and more recently has received the AAPG George C. Matson Award for his paper entitled ‘Delineating stratigraphic features via cross-plotting of seismic discontinuity attributes and their volume visualization’, being adjudged as the best oral presentation at the 2010 AAPG Annual Convention held at New Orleans, the ‘Top 10 Paper’ Award for his poster entitled ‘Extracting meaningful information from seismic attributes’, presented at the 2009 AAPG Annual Convention held at Denver, the ‘Best Poster’ Award for his paper entitled ‘Seismic attributes for fault/fracture characterization’, presented at the 2008 SEG Convention held at Las Vegas, the ‘Best Paper’ Award for his paper entitled ‘Curvature and iconic Coherence–Attributes adding value to 3D Seismic Data Interpretation’ presented at the CSEG Technical Luncheon, Calgary, in January 2007 and the 2005 CSEG Meritorious Services Award. He and his colleagues have received the CSEG Best Poster Awards in successive years from 2002 to 2005.
He is a member of SEG, CSEG, CSPG, CHOA (Canadian Heavy Oil Association), EAGE, AAPG, APEGGA (Association of Professional Engineers, Geologists and Geophysicists of Alberta) and TBPG (Texas Board of Professional Geoscientists).
Lazear, G. D., 1993, Mixed-phase wavelet estimation using fourth-order cumulants, Geophysics, 58, 1042-1051.
Mendel, J. M., 1991, Tutorial on higher-order statistics (spectra) in signal processing and system theory: theoretical results and some applications, Proceedings of the IEEE, 79, 278-305.
Misra, S. and Sacchi, M. D., 2006, Nonminimum phase wavelet estimation by nonlinear optimization of all-pass operators, Geophysical Prospecting, 55, 223-234.
Misra, S., 2008, Global optimization with application to Geophysics, PhD thesis, University of Alberta, Canada.
Porsani, M. J., and Ursin, B., 1998,Mixed phase deconvolution, Geophysics, 63, 637-647.
Velis, D. R., and Ulrych, T. J., 1996, Simulated annealing wavelet estimation via fourthorder cumulant matching, Geophysics, 61, 1939-1948.