Mixed-phase wavelet estimation – A case study

Introduction

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.

Case Study

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.

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 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.

Acknowledgements

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.