Blind Deconvolution of Seismic Data Based on the Spearman’s Rho

In this paper, we propose a novel seismic blind deconvolution approach based on the Spearman’s rho in the case of band-limited seismic data with a low dominant frequency and short data records. The Spearman’s rho is a measure of the dependence between two continuous random variables without the influence of the marginal distributions, by which a new criterion for blind deconvolution is constructed. The optimization program for new criterion of blind deconvolution is performed by applying Neidell’s wavelet model to the inverse filter. The noise-free and noisy synthetic data, onshore seismic trace in the Ordos Basin, and offshore stacked section in the Bohai Bay Basin examples show good results of the method.


Introduction
In seismic exploration, a recorded seismic trace is often represented as the result of a superposition of wavelets of constant shape weighted by the reflectivity series, thus it can be modeled as a convolution between the source signature and the reflectivity series.Deconvolution is an important process by which the reflectivity series can be estimated from the recorded trace, vertical resolution of the seismic image will be enhanced.
As the reflectivity series, the wavelet is unknown, only seismic trace is accessible.We have the classical blind deconvolution problem.
Statistical deconvolution appears to be one of powerful tool for dealing with practical situation.Many blind deconvolution techniques, under various assumptions about the seismic wavelet and reflectivity series, have been developed with respect to different criteria in recent decades [1]- [5].
In the case that the wavelet has the zeros on the unit circle, blind deconvolution is an ill posed problem, one of the criterions to apply to inverse filter is to maximize the independence between the samples from the deconvolved output [6].Larue et al. (2006) [7] proposed a new blind single-input single-output (SISO) deconvolution method based on the minimization of the mutual information rate of the deconvolved output.The mutual information is a criterion for measuring independence of random variables without the influence of the marginal distribution, if a stochastic process has a vanished mutual information rate.Its components are mutually independent.
To take the band limited nature of seismic wavelets and the presence of noise into account, van der Baan and Pham (2008) [8] present a modification of the mutual-information rate, whitening the deconvolution output only within the wavelet pass band, to prevent noise amplification, Wiener filtering is introduced.
In the above entropy-based approaches, one key point is to estimate the probability density function of the deconvolution output.However, estimating the PDF by nonparametric method needs certain length of the data to provide a good estimation.A seismic signal is treated to be stationary just within a limited time window due to intrinsic attenuation of the earth medium, thus these entropy-based approaches seem to result in a limited performance in the case of short data records.
In statistics theory, the mutual information is the one of measures to show the dependence between two continuous random variables, unlike linear correlation which is not invariant under almost surely strictly monotone functions, Spearman's rho ( , X Y ρ ), as measures of dependence between two continuous random variables X and Y, satisfy the following properties: • if α and β are almost surely strictly monotone functions on RanX and RanY, respectively, then: We can see that Spearman's rho is a useful tool when dealing with blind deconvolution problem.The overall objective of this paper is to construct a new criterion based on Spearman's rho for seismic deconvolution, and propose a new blind deconvolution algorithm.Furthermore, this new method can handle the data with small samples, which means it has ability to deal with nonstationary data.

Spearman's Rho
Nelsen (1999) proved that if x and y are observations from two random variables X and Y with distribution functions F and G, respectively.Lets u = F(x) and v = G(y), the expression for Spearman's rho can be written in the following form [9]: For the stationary seismic signal { } n x , each sample of the trace has the same distribution functions.The em- pirical distribution function can be applied to calculate the n u corresponding to the sample n x : where { } is the indicator function of event A, N is the length of the data .

The Inverse Filter
Following the work by van der Baan and Pham (2008), we can show that if the wavelet is known, Wiener filter achieves an optimum solution for inverse filter.In the frequency domain, it can be written as: W ω is the wavelet, the superscript of the numerator denotes the complex conjugate, 2 ε is a normal factor.
In the blind deconvolution context, the wavelet is unknown, and need to be estimated.Therefore, we adopt Neidell's wavelet model to constructing an inverse filter, which is derived from practical observation, and its z transform has the form [10]: The relationship between the n and m is: p ω denotes the peak frequency of the amplitude spectrum of this wavelet.
So, if the peak frequency of the amplitude spectrum of real wavelet is known, only one parameter needs to be determined.In practice, we usually do not know this peak frequency.Therefore, we have following value as a substitution for it: where FT denotes Fourier transform.

New Deconvolution Algorithm
Based on the above results, we construct a new frequency domain criterion corresponding to the inverse filter A for blind deconvolution: ( ) ( ) l denotes the time delay.After determining the p ω , ( ) J A is just the function of n or m, so we have: ( ) , Thus our new deconvolution algorithm consists of the following steps: of the inverse filter with the different parameter n, in which the superscript n denotes the shape parameter of the inverse filter; n u is given by (2), and the inverse filter is given by ( 3) and (4).
• For a suitable l , compute the cost function J corresponding to the shape parameter n of the inverse filter.
• Search the maximum of the cost function to get the optimum inverse filter.

Simulation Experiments
We have simulated the seismic trace by convolving a super-Gaussian reflectivity with a 40 Hz Ricker wavelet.The reflectivity is generated as ( ) b n , where ( ) b n are independent normal random variables with zero mean and variance 0.28.The sample interval of the synthetic trace is 1 ms.
In order to assess the performance of our method for the data with small samples, we set the length of the synthetic trace about 200 ms.The estimation results are plotted in Figure 1 and   To test our algorithm for noisy signals, we employed a noisy synthetic trace.The signal-to-noise ratio of the noisy data is 15 dB.The results are shown in Figure 3.The wavelet (blue line) is well estimated.The original spectrum (blue) ranges 0 Hz to 80 Hz, and the deconvolved spectrum (red) range is 0 Hz to 100 Hz.The spectrum of the deconvolved trace is definite broadened.Meanwhile, the result of the deconvolved trace shows the higher signal-to-noise ratio, which demonstrates an anti-noise feature of our method.

Single Trace Data
We applied our method to a real onshore single seismic data, which was obtained from the Ordos Basin, China.The original single trace consisted of 750 samples.
As shown in Figure 4, the resolution of the deconvolved seismic trace is clearly improved.The original spectrum (blue line) ranges 0 Hz to 120 Hz, and the deconvolved spectrum (red line) ranges 0 Hz to 180 Hz.In contrast with the spectrum of the original trace, the spectrum of the deconvolved trace is obviously broadened in Figure 5.

Stacked Section
An offshore stacked section is to assess the potential of applying our method.The stacked section was obtained from China.There are 50 traces in this seismic section.The sampling interval of data is 1 ms. Figure 6 displays the offshore data and result after deconvolution using our method.As shown in Figure 6, the resolution of the deconvolved section is clearly improved.

Conclusion
We have presented a new blind deconvolution method based on the measure of the dependence between two    continuous random variables-Spearman's rho, by which a new criterion for blind deconvolution can be constructed.Synthetic and two real-data examples have shown the ability of our approach in the presence of the band-limited seismic data with a low dominant frequency, short data records, and low signal-to-noise ratio.

Figure 2 .
As shown in the figures, the wavelet (blue line) is very precise estimated.The original spectrum (red) ranges 0 Hz to 100 Hz, and the deconvolved spectrum (blue) range is 0 Hz to 200 Hz.The spectrum of the deconvolved trace is obviously broadened.

Figure 1 .
Figure 1.The results of deconvolution in the noise-free case.

Figure 2 .
Figure 2. The comparison between the original trace and the deconvolved trace.

Figure 3 .
Figure 3.The results of deconvolution for a noisy synthetic trace.

Figure 4 .
Figure 4.The comparison in the waveform of the original seismic trace and the deconvolved trace.The resolution of the deconvolved seismic trace is clearly improved.

Figure 5 .
Figure 5.The contrast in the amplitude spectrum of the original seismic trace and the deconvolved trace.The spectrum of the deconvolved trace is obviously broadened.

Figure 6 .
Figure6.The offshore data and result after deconvolution using our method.