# The Leading Edge

- Copyright © 2010 Society of Exploration Geophysicists

## Abstract

Accurate wavelet estimation 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. Because the Earth's reflectivity and seismic noise are both unknown, the wavelet estimation process is not easy. 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 are 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 (stationary assumption).

Accurate wavelet estimation 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. Because the Earth's reflectivity and seismic noise are both unknown, the wavelet estimation process is not easy. 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 are 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 (stationary assumption).

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 is retained during processing. However, it is not applicable for estimating wavelets from sources giving a mixed-phase signature. For example, deconvolution of nonminimum-phase data with a minimum-phase wavelet will leave 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 nonminimum-phase signals. Also, these measures work well for Gaussian probability distribution of amplitudes, and so will not yield accurate results for non-Gaussian distributions.

Non-Gaussianity in the seismic data could arise from a nonminimum-phase source signature, noise in the data like swell noise, and a nonlinear 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, reveal not only amplitude information, but also phase information (Mendel, 1991).

In this paper, we address this issue through the use of higher-order statistics such that the phase components in the data are more accurately estimated and removed.

## Higher-order statistics for wavelet estimation

The 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 third-order cumulant is a two-dimensional correlation function. For a Gaussian process, all cumulants above the second-order are zero, but are nonzero for non-Gaussian processes. Thus these two statistics are not suitable for recovering a nonminimum 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 transform. Lazear (1993) and Velis and Ulrych (1996) 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 all-pass 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 parameterization of the all-pass 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 and keep the length as a constant parameter.

Seismic data are 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 a minimum-phase wavelet and the all-pass operator. Thus the seismic data are 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 which we subsequently refer to as the whitened data. The whitened data are thus represented as the convolution of the underlying reflectivity series and the all-pass 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 L_{2}-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 fourth-order trace cumulant of the whitened data (data obtained after deconvolution with a minimum-phase wavelet) and the normalized fourth-order moment of the all-pass operator. In the present context, 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. Thus in the present optimization problem where the shape of the cost function is not known, a global optimization algorithm is a preferred choice. A 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 models is computed by taking the ratio of the corresponding maximum-phase sequence and the minimum-phase sequence.

## Application to poststack 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 zone a mixed-phase wavelet is estimated. The data in each subdivision 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 P-wave log curve.

Figure 1 shows the workflow for the method outlined above.

## Case study

The data volume picked for testing the algorithm had been processed using a conventional sequence. Four different locations on the poststack volume were selected for the estimation of the mixed-phase wavelets in the same time interval (500 ms) for each.

Figure 2a shows a segment of a seismic section around location 1. Figures 2b–d show the estimated minimum-phase wavelet, the estimated all-pass operator, and the estimated mixed-phase wavelet at 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 broadened the bandwidth. Figure 2e shows the phase-corrected data for location 1. Figures 3, 4, and 5 show a similar set of images for locations 2–4. Notice that, for each set of images, the mixed-phase-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 confidence in ascertaining if the resolved reflections correlate well and if they are authentic. Notice in Figure 6 that the correlation of the section in Figure 6b is better than the one in Figure 6a.

## Conclusions

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. The method adopted in estimating and hence removing the spurious phase involves estimation of the coefficients of an all-pass operator from 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 fourth-order normalized trace cumulant and the fourth-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 from 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 mixed-phase wavelets. The data are then deconvolved with the average mixed-phase wavelet and correlated with the P-wave log data.

## Acknowledgments

We acknowledge Mauricio Sacchi for his help and suggestions during the development of the algorithm. We also acknowledge an anonymous provider of the data shown in the example and Arcis Corporation, Calgary, Canada for support of this work.

## Footnotes

Coordinated by Alan Jackson