A Refinement of the Generalized Blind Deconvolution Method Based on the Gaussian Mixture and its Application

A two-channel blind deconvolution method is developed and clearly shown to be capable of recovering the unknown reflectivity sequences from two observation data sets. The blind deconvolution technique, based on the Gaussian mixture model for the reflectivity sequence, had previously been developed, but it has had oue shortcoming which is that the order of the inverse filter could not be determined, thus preventing the removal of the source wavelet from the observation data. Here, this imperfection is over­ come by using a technique that identifies the optimum order by using two time sequences along with the cost function. Two synthetic seismograms · are used to examine the technique, and very good results are obtained. We, therefore, apply the method to analyze seismic exploration data and deter­ mine the location of the reflective signal where the reflectivity of each chan­ nel is consistent. This clearly indicates that this enhanced method is, indeed, effective for determining the true layering effect when it is applied to seis­ mic exploration data. (


INTRODUCTION
Based on a convolution model, observed seismic data are considered as the result of a source wavelet convoluting with a reflectivity sequence.Although one objective of seismic processing is to identify a layered structure of the subsurface from a reflected seismic data, the best way to estimate the reflectivity sequence and the source wavelet had never been devel oped until now.Estimating both of these directly from observed seismic data without any prior information about them is known as blind deconvolution.After all, seismic deconvolution is 1 Institute of Seismology, National Chung Cheng University, Chia-Yi, Taiwan, ROG essentially a blind process since it is difficult to quantitatively measure the source wavelet generated by an explosion on the ground.
Built upon the vast experience of geophysicists, the conventional deconvolution method most often assumes that the source wavelet is the minimum phase or that the model of the reflectivity sequence agrees with the Bernoulli-Gaussian model (Robinson 1975;Mendel 1990).Although these assumptions are frequently reasonable in practice, they certainly cannot be made in every situation.The conventional method cannot recover full reflectivity; on the contrary, it only recovers the white all-pass component of the reflectivity (Saggaf and Toksoz 1999).
To resolve the problem of blind deconvolution, a number of algorithms have been proposed.Gurelli and Nikias (1995) formnlated an algorithm based on the eigenvector method to esti mate the orders and the root locations of the channel transform functions for solving deconvolution.Haykin (1994) presented an extensive discussion of some algorithms by using Higher Order Statistics (HOS).Velis and Ulrych (1996) used the fourth-order cumulant func tion to find the mixed phase models and applied it to seismic deconvolution.More recently, Pflug (2000) discussed the relationship between the signal passband, the trispectral domain and the fourth-order cumulant function in the application of the fourth-order deconvolution method.Santamaria et al. (1999) developed an algorithm which takes into consideration the Gaussian mixture model to identify the reflectivity sequence in the context of iterative seismic deconvolution, and this can effectively be used to delete the diffraction effects.Similarly based on an iterative approach, new criteria for the multichannel blind deconvolution algorithm have been proposed by Inouye and Sato (1999).
In this paper, the source wavelet will be removed to estimate the reflectivity functions by applying the adaptive Gaussian mixture algorithm (Santamaria et al. 1999) in conjunction with the two-channel cost function (Furuya and Kaneda 1997) which can obtain the optimal deconvolution filter order.The performance of this algorithm will be tested in detail by syn thetic data.In the final stage, we apply this algorithm to seismic exploration data.

METHODOLOGY
In seismic exploration, the observed seismic data z;(t) are customarily represented by: z;(t)=x1(t)*h(t)+n1(t) , (i=l,2) (1) where x1(t) represents reflectivity function, h(t) is the seismic source wavelet which is possi bly the non-minimum phase, n;(t) is noise, and * denotes convolution.Employing the adap tive Gaussian mixture model, Santamaria et al. (1999) proposed an algorithm to estimate the reflectivity x;(t) of a single channel.This algorithm provides an inverse filter of the source wavelet which makes it possible to remove source effect and eliminate the diffraction effect through a series of iteration calculations.Then, the estimated reflectivity can be determined.This procedure leads to excellent results, especially if the reflectivity follows the Bemoulli Gaussian model.Given that different deconvolution filter orders produce different results, however, the determination of the filter order becomes critical in Santamaria's method.To solve this problem, in this research we use Furuya' s method to determine the most reasonable deconvolution filter order.
Assuming the initial deconvolution filter order is n, we apply the Santamaria' s method to deconvolnte z;(t) and come up with two estimated reflectivity values, x1(t,n) and x,(t,n).
Following the multi-channel inverse filtering algorithm (Miyoshi and Kaneda 1988), we then derive the two-channel inverse filters g1 (t,n) and g2(t,n) for xJt,n) and x2 (t,n), respectively, by solving the equation: where d =[l, 0, 0, 0, ... , 0], and the order of dis 2n .gJt,n) and g,(t,n) are applied to recover the estimated source wavelet h(t,n) which can be written as: (3) In order to evaluate the fitness between the estimated and observed signals, Furuya and Kaneda (1997) proposed the cost function.It is defined as: (5) where a is an arbitrary constant.From the above deductions, we pin-point the optimal deconvolution filter with order n that minimizes the cost function PE(n ) and, accordingly, obtain the best estimations xJt,n) and x2(t,n) of the reflectivities x/t) and x2 (t), respectively.

Theoretical Test
In order to verify the validity of the method, two synthetic data sets are used, and these are shown in Fig. 1.A simple Ricker wavelet is regarded as a source function that is individually convoluted with two different reflectivities with the Gaussian distribution.These two results after convolution are displayed in Figs.le, e.Using this method, we analyze two synthetic seismograms to obtain the optimum deconvolution filter order.The searching range of the deconvolution filter order, an odd number, is from 7 to the terminal point.According to Santamaria et al. (1999), the terminal point can be determined by monitoring the iteration times of each deconvolution filter order.As the order increases, the iteration times stabilize at a certain value.Here we select the order as the termi nal point of searching range at which the iteration times of the two channels have a stable value.Figure 2 shows the relationship between the deconvolution filter orders and the itera tion times of both channels.The circles indicate the terminal points.Figure 3 shows the cost function PE(n) vs. the optimum order, and the minimum value of the cost function appears at the order 39.
Figures 4, 5 show the comparisons between the synthetic reflectivity ((a)) and the esti mated reflectivities ((b) -(i)) obtained from different deconvolution filter orders for channels I and 2. It is apparent that different orders lead to different results.The best fit between the synthetic and estimated reflectivities is obtained when the deconvolution filter order, n, is 39.
It is also shown that it works well for getting the estimated reflectivity when the deconvolution filter order n is equal to 69 (Fig. 5).However, good results cannot be obtained in Fig. 4 when n is 69.Some differences are still obvious between the synthetic and the estimated reflectivities.
This indicates that the method can ascertain the optimum order for obtaining the most reason able estimated reflectivity of the synthetic seismogram.

Application
At this stage, we employ our method to analyze a seismic exploration data set which was collected at a field test site in Oklahoma and is shown in Fig. 6a.The sampling time of each trace is 0.004 sec, and the distance between any two receivers is 50 m.An automatic gain correction (AGC) with a 500 ms window is applied to each trace (Fig. 6a).According to Shieh and Herrmann (1990), the data should have had a reflection signal at about 0.55 seconds, but it was mixed with a strong ground roll.As mentioned earlier, the seismic deconvolution method (Santamaria et al. 1999) can eliminate the diffraction effect which creates scattering noise.Thus, we use this method to identify the reflectivity of each trace and locate the position where the reflectivity of each trace is consistent.The reflection signals appear at .that position, since each reflector represents a change in the acoustic impedance at the interfaces.Deconvolution filter order Fig. 2. Relationship between the deconvolution filter orders and the iteration times of the deconvolution algorithm for channels 1 and 2. The small circles represent the position of the terminal points. 120 In order to calculate the deconvolution filter order for each channel, the neighboring two traces were chosen and applied in the proposed method.The optimum deconvolution filter order is estimated while the cost function is minimum.Figure 7 shows the calculated optimal deconvolution filter order for each channel of seismic data in Fig. 6a.The channel number 1-19 is relative to the seismic trace from bottom to top.When those orders are employed, each trace is deconvoluted to get the reflectivity function shown in Fig. 6b.It is clear that the major reflectivity position for each trace is consistent and appears at about 0.55 second.Of particular interest is that in analyzing the same data set, Shieh and Herrmann (1990) demonstrated that a polarization filter could reject ground roll, and the desired reflection signal at about 0.55 seconds.In other words, an excellent agreement found between our results and those of Shieh and Herrmann (1990) indicates that the method employed here provides an invaluable tool to esti mate the true position of a reflection signal.
Deconvolution filter order le, e.The optimal deconvolution filter order, corresponding to the mini mum value of the cost function, n, is equal to 39.

CONCLUSIONS
In this paper, we employ Furuya's (1997) method to solve the problem resulting from the fact that the order of the deconvolution filter cannot be ascertained by using another trace of the same source.Once the optimum order is found, we utilize it to obtain the estimated reflectivity.We need to consider information about the subsurface when the correct reflectivity is obtained.The theoretical test shows that the method can achieve the blind deconvolution well to identify the estimated reflectivity.
The method developed here is based on the reconstruction of reflectivity with a zero mean Gaussian mixture distribution while the source wavelet is non-minimum phase.The Gaussian mixture distribution is composed of both a broad and a narrow Gaussian model.The broad one displays the major layering effect, whereas the narrow one models the scattering noise.In other words, the blind deconvolution algorithm proposed by Santamaria et al. (1999) has the ability to eliminate the backscattering as well as the noises and models the mayor layering reflectivity effects.Thus, we employ this method to isolate the reflectivity from a set of reflected seismic signals.There is very strong ground roll to cover the desired reflection   signal.These data are deconvoluted by using the technique presented here.The results unam biguously show that the major layering effect appears at about 0.55 sec from the distribution of reflectivity.According to Shieh and Herrmann (1990), the data set has a reflection signal at 0.55 sec, and ground roll is regarded as the scattering effect.Therefore, we can confirm that the technique is a valuable tool to determine the true layering effect when it is applied to seismic expl9ration data.Channel number of seismic data Fig. 7. Optimal deconvolution filter order for each channel of the seismic in Fig. 6a.The channel number 1-19 is relative to the seismic trace from bottom to top.
Figure le is obtained by convoluting the source function (Fig. la) with the reflectivity function shown in Fig. lb.Similarly, Fig. le shows the results after convoluting the same source function (Fig. la) with the other reflectivity function shown in Fig. Id.Besides this, 20-dB white noise is added to Figs. le, e.

Fig . 1 .
Fig . 1. Example showing the deconvolution of two signals.The Ricker wavelet, shown in (a), is selected as the source time function.(b) and ( d) are two different reflectivities with the Gaussian distribution.(c) is the result of (a) convoluting with (b).Similarly, (e) is the result of (a) convoluting with (d).In addition, 20dB noise is added to (c) and (e).

Fig. 3 .
Fig. 3. Cost function PE(n) of the two synthetic seismograms shown in Figs.

Fig. 6 .
Fig. 6.Seismic exploration data and their deconvolution results.The seismic data set shown in (a) was collected at a field test site in Oklahoma.All traces had an automatic gain correction (AGC) with a 500 ms window.The estimated reflectivity for each trace using the deconvolution method developed here is shown in (b ).The major reflection signals are consis tent and appear at about 0.55 sec.