A new approach for extracting the amplitude spectrum of the seismic wavelet from the seismic traces

In reflection seismology, knowing the seismic wavelet is important both for processing seismic data and for modeling the seismic response. There are two approaches to obtain the seismic wavelet. One approach is deterministic and the other is statistic. This work belongs to the second category. A seismic wavelet is determined by the product of amplitude spectrum and phase spectrum. A conventional method uses a two-step procedure to estimate the seismic wavelet. The first step is for the amplitude spectrum, and the second step is for the phase spectrum. So extracting the amplitude spectrum of a seismic wavelet (ASSW) from the amplitude spectrum of a seismic trace is a key step. The commonly used methods are correlation-based method, the log-spectrum-averaging method and spectrum-shaping method. All these methods assume the reflection coefficient sequence is white, which may not be valid under some conditions. In this paper, we propose a new approach to obtain ASSW without the whiteness assumption about reflectivity. We define an operator on a properly chosen function space and prove that the operator is contractive in this space. Then, we convert the problem of estimating the ASSW into one of finding the fixed point of the operator. We give the algorithm in detail based on the contraction operator mapping (COM method). We compare our method with the widely used methods by synthetic signals in which the reflectivity is not white. The results show that our method performs satisfactorily for nonwhite reflection series, on which other methods do not work well. Moreover, our method is robust for the noise and the frequency interval on which COM works.


Introduction
Reflection seismology is used widely in oil and gas exploration, environmental and civil engineers, imaging the lower crust and mantle, etc. The reflection seismology method uses the reflection data received on the earth surface to invert the features of underground medium (Syems 2009).
This method is based on the mapping rules between reflection seismic data and the parameters representing the underground medium. Describing the mapping rules in mathematics is called the representation of the reflection seismology.
There are two approaches for the representation of reflection seismology. One is the deterministic approach and the other is the statistical approach (Robinson and Treitel 2000). The deterministic approach consists of utilizing physical theories of wave propagation involving the solution of integral and differential equations satisfying boundary and initial conditions. The statistical approach utilizes statistical theories of time series leading to the expression of the dynamics as a statistical fact (Robinson and Treitel 2000).
In reflection seismology, knowing the seismic wavelet is important both for processing seismic data and for modeling the seismic response. Corresponding to the two approaches to the representation for reflection seismology, there are two approaches to obtain the seismic wavelet. One is the deterministic approach and the other is a statistic approach.
We first review the deterministic approach. Weglein and Secrest (1990) proposed a wavelet estimation method based on wave theory. Their method can obtain the total wavelet, including the phase and source-array pattern. In addition, the theory allows either an acoustic wavelet (marine) or an elastic wavelet (land), so the wavelet is consistent with the earth model to be used in processing the data. Osen et al (1998) start from the general result by Weglein and Secrest (1990) and suggested a new algorithm for seismic wavelet estimation from marine pressure measurements. Their method requires neither the information of the signature nor any knowledge of the earth below the receivers. Tan (1999) proposed a method to estimate the wavelet spectrum based on the acoustic wave equation. His method does not require statistical assumption, of which the most important are whiteness of the earth's reflectivity series and the minimum-phase character of wavelet.
In this paper, we are interested in the statistical approach which is based on the convolution model of a seismic trace. A seismic wavelet is determined by the product of its amplitude spectrum and the phase spectrum. A conventional method uses a two-step procedure to estimate the seismic wavelet. First step is for the amplitude spectrum, and the other is for the phase spectrum. We list some typical methods as follows. Levy and Oldenburg (1987), Longbottom et al (1988) and White (1988) describe such a technique for stationary data. Baan and Pham develop a technique to estimate the seismic wavelet from the seismic trace in which the bandwidth of the seismic wavelet is narrow to very narrow. That is, the bandwidth of the seismic wavelet is similar to the wavelet principal frequency (Van der Baan and Pham 2008). In their method, first they average the amplitude spectrum of the observed traces and use it as the amplitude spectrum of the seismic wavelet. By setting the phase to zero, they obtain a zero-phase wavelet. Taking the mutual information rate minimum as a criterion and using this zero-phase wavelet as the initial guess in the inversion, they obtain the total seismic wavelet and reflectivity series by solving an inverse problem. In their method, an accurate ASSW is very important since it leads to fast convergence.
Some authors propose another representation of a reflection seismogram (Ricker 1940, Dobrin andSavit 1988). In this representation, the reflection response is the convolution of the reflectivity with a time-varying sequence of wavelets. Each time-varying wavelet in the sequence is associated with one and only one interface. Mathematically, a single interface is associated with each discrete instant on the trace. That is, each interface is associated with its own time-varying wavelet. Robinson calls the representation as time-varying convolution model of the seismic trace (Robinson andTreitel 2000, 2008). Other author call it Nonstationary convolution model (Margrave 2011). As an approximate method of time-varying convolution model, the record is divided into time intervals which may be considered approximately stationary. In the past we proposed several methods which can represent the timevarying convolution model as combination of approximately stationary segments. In these approaches, how to find the ASSW in each time interval is a key (Gao et al 2009, Gao and Zhang 2010, Zhang and Gao 2011, Wang et al 2013.
Under the assumptions that reflection coefficient sequence (RCS) is white, some methods have been proposed to extract the amplitude spectrum, like the correlation-function method (CF method), the log-spectrum-averaging method (LSA method) (Otis andSmith 1977, Van der Baan 2008) and spectrum-shaping method (SS method) (Rosa and Ulrych 1991). However, when the underground medium is inhomogeneous transversally, Van der Baan's method fails. The SS method supposes the ASSW satisfying a sort of function form which contains some parameters that can be determined by means of the least-square fit between the amplitude spectrum of the seismic trace and the function (Rosa and Ulrych 1991). In practice, choice of a proper function is difficult, so does the proper frequency range in which the least-square fit is performed. As mentioned above, these methods work when the RCS is whiteness. However, this assumption may not be true in the real data. Some researchers found that the RCSs meet blueness (Wakden 1985, Walden and Hosken 1986, Painter et al 1995, Saggaf and Robinson 2000. In order to extract the ASSW for a RCS which do not meet the whiteness assumption, we need to develop new methods. Neidell (1991) showed that the ASSW is a smooth and unimodal function. For an explosive source, Ricker's work reaches the same conclusion (Ricker 1953a(Ricker , 1953b. The properties of the ASSW are similar to that of a probability density function (PDF). So extracting the ASSW from the amplitude spectrum of a seismic trace is very similar to estimate the PDF from random variable samples like that pioneered by Scheffnerl and Runde (2003). Inspired by the work of Scheffnerl and Runde (2003), in this paper, we develop a new approach to extract the ASSW. We define an operator and prove that it is contractive. Thus, we convert the estimation of the ASSW into finding the fixed point of the operator. We give the theory of the algorithm in detail based on the contraction operator mapping (COM method). By comparing our method with the widely used methods on synthetic signals, we conclude that our method is valid for nonwhite reflection series. Especially, our method does no need to choose the beforehand function form for ASSW, and is robust for the frequency range on which COM works. This paper is organized as follows. In section 2, we construct a contraction operator, and put the existence and uniqueness of its fixed point in appendices A and B. In section 3, we propose a method to estimate the amplitude spectrum of seismic wavelets based on the proposed operator. In section 4, we validate the proposed method using synthetic data and compare to usually used method. Section 5 is devoted to some conclusions and discussions.

The construction of the contract operator and its fixed point
, : on , . (2.1) Obviously, K is closed. Over the space C a b , [ ] and its subset K, we use the canonical norm It is clear that (2.4) Theorem 1. Let the set K and the operator P be defined as in (2.1) and (2.3). Then, we have the inequality with τ τ being understood to be 1 for τ=0. Therefore, for c p sufficiently small, P is a contractive mapping from K into K; consequently, P has a unique fixed point in K. Moreover, the fixedpoint iteration method A proof of theorem 1 is given in appendix A. The next result is proved in appendix B. [ In the iterative sequence defined in theorem 1, for = ], + f n 1 is a unimodal function by theorem 2, and + f n 1 is smoother then

Estimating the amplitude spectrum of seismic wavelet
In this section, we propose a method to estimate the amplitude spectrum of seismic wavelets based on the contract operator defined in (2.3).

The convolution model of seismic traces
Let s t ( )denote the seismic trace, w t ( ) the seismic wavelet, r t ( ) the reflectivity sequence, and n t ( ) a random noise, where w t ( ), r t ( ), n t ( ) all belong to L R 2 ( ). Assume w t ( ) is compactly supported. Then we have (Robinson and Treitel 2000) Equation (3.1) is called the convolution model of a seismic trace in time domain (Robinson 1967, Robinson andTreitel 2000). Applying Fourier transform to both sides of (3.1), we get , ( ) ( ) ( ) are the Fourier transforms of the corresponding time functions.

Determination of parameters in the operator
In equation (3.2), ignoring the noise, we have Suppose the effective bandwidth of the seismic trace is defined in a b , . where a b and are real numbers. Taking the absolute value of both hand sides of (3.3), we get In order to find the ASSW, we let the operator [ ] P f x ; defined by (2.3) act on ω S 0 ( ). In practice, if we select a proper subinterval in a b , In this case, if we let δ = 0, theorem 1 still holds good. So now Now we discuss how to determine the parameters, α β c , , p in the (3.8). We use L R 2 ( ) norm to measure the difference between ω S 0 ( ) and ω P S ; Noticing that the deference depend on the parameters, α β c , , p nonlinearly. In order to obtain a linear inverse problem, we apply logarithmic transformation to ω S 0 ( ) and ω P S ; ]. Then we construct the objective function . m → is vector which represents the model parameters in the inverse problem. They are defined as For the time being, we ignore constrain term in the (3.9), solve the problem. Its solution in the least-square sense (Aster et al 2004 (3.10)

Estimation of the ASSW of a seismic trace
Based on the operator (3.8) with the parameters from (3.10), we construct an iterative algorithm to extract the amplitude spectrum of the seismic wavelet from the amplitude spectrum of a seismic trace. For the seismic trace s t ( ), its Fourier transform is ω S( ). Let ω S k ( ) be thekth mapping. The iteration algorithm can be expressed as follows, which is its fixed point of operator P. Obviously, this fixed point is the scaled amplitude spectrum of the wavelets of the seismic trace (see (3.5)). Since, our method is based on a contract operator mapping, so we call it COM method.

The algorithm flow of COM method
The COM method consists of three steps as follows: Step 1: Determining the parameters in the operator defined by equation (3.8) using equation (3.10)); Step 2: Extracting the ASSW, ω W 0 ( ), by equation (3.11); We give the algorithm as table 1 Taking Syn.2 (see section 4) as an example, we explain the COM method and investigate its performance.
The amplitude spectrum of the synthetic seismic trace is displayed in figure 1 (a). Figures 1(b)-(f) show the results corresponding to different iterations from 1 to 5 in turn. From these results, we can observe that: (1) Each of the curves produced during the iteration process (i.e. (b)-(f)) has only one peak frequency, as expected. This is because the expression of the right hand side of equation (2.3) corresponds to a unimodal function (appendix B).

Synthetic seismic traces
In order to make a comparison for our method with those widely used methods, we generate some synthetic seismic traces using the convolution model defined in (3.1). Here we ignore the noise.

• Seismic wavelet
We take the Ricker wavelet as the seismic wavelet. In time domain, it is formulated as (Ricker 1953a): here f M is its peck frequency. In all models in this paper, we set = f H z 40 M .

• The Reflection Coefficient
Many researchers explore the nature of reflection series statistically. Robinson (1967) proposed a convolution model, and gave a deconvolution method. In his model, the elements of RCS meet the same Gaussian distribution. Following Robinson, geophysicists Input the amplitude spectrum of the seismic trace ( ) ω S 0 Set = k 1; Step 1: Apply nonlinear transform Step 2: Normalization: Step 3: Step 4: Inverse of the nonlinear transform ( ) / = S S k k p 1 Step 5 find out that the RCS meet non-Gaussian distribution. For example, Painter et al (1995) found that the RCSs from 14 wells in Australia have a statistical character consistent with a non-Gaussian scaling noise model. Walden and Hosken (1986) examined the amplitude distribution of a series of primary reflection coefficients. They concluded that the distribution is always essentially symmetric, but has a sharper central peak and larger tails than Gaussian distribution. As representatives, we take the four kinds of random series as RCS with different distributions: RCS 1, Bernoulli-Gaussian distribution (Rosec et al 2003); RCS 2, α − Stable distribution (Samorodnitsky and Taqqu 1994, Painter et al 1995, Weron 1996, 2004, Pavel et al 2005; RCS 3, blue noise (Walden and Hosken 1986, Painter et al 1995, Saggaf and Robinson 1995; and RCS 4 is obtained by using the well data from the sonic and the density logs.

Compare the COM method with CF and SS methods
In this subsection, we apply the COM method, CF method and SS method to all the four synth etic traces to compare their performance. For the convenience of comparison, we define the misfit function as follows: where ∆ ASSW denotes the value of maximum misfit (VMM).
• Application to Syn.1 Figure 2(e) shows the estimated ASSW by CF method and SS method, respectively. Figure 2(g) plots their misfit functions. The estimated ASSW (dot dash line) by CF method is oscillating and has a significant difference with the true one. The misfit function, ∆ ASSW , has a large value almost everywhere; the VMM has a large value too. The estimated ASSW using the SS method (blue line) also has obvious difference with the true one, and the misfit function displays oscillation as well. However, the VMM for the SS method is much less than that for CF method. It can be expected that the CF method leads to a poor result. This is because CF method only works when the RCS is white, which is violated in Syn.1. Lastly, we investigate the result produced by the COM method. Figures 2(f) and (h) display the results by this method. After five iterations, the estimated ASSW is almost the same as the true one; The VMM corresponding to the misfit function becomes less and less as iteration continues. At the 5th iteration, the VMM is already very little. With the iteration continue, the VMM will approach to zero. We noticed that the SS method is sensitive to the length and the location in the frequency axis of the window in which it works. It is tricky to choose them properly. The COM method has an advantage over the SS method in this respect. It is robust to the length and the location of the window in the frequency axis in which it works. In summary, the COM method has advantages over other two methods in accuracy and robustness. • Application to Syn.2 For the 2nd example, figure 3(e) shows the estimated ASSWs by the CF method and SS method. The comparison between figure 3(e) with figure 2(e) shows that the CF and the SS methods work well. Two methods work better than what they do for Syn. 1. But in figures 3(e) and (g), there is still oscillation. In addition, the oscillations appear at different frequencies for the two methods. The former lies from 0 Hz to 150 Hz, and the latter lies from 0 Hz to 120 Hz. The VMM of the former is approximately 0.01; that for the SS method is less than 0.018. Thus, neither method can estimate the ASSW correctly.
The figures 3(f) and (h) displays the results by COM method. After five iterations, the estimated ASSW is almost the same as the true one. Similar to Syn. 1, the COM method can extract the ASSW accurately. • Application to Syn.3 Figure 4(e) shows the estimated ASSWs by CF and SS method, respectively. Figure 4(  In general, from the application to the four different type synthetic traces, the COM method has extracted the ASSW with much higher accuracy than other methods can do.

Conclusions and discussions
In this paper, we construct an operator to extract the amplitude spectrum of the seismic wavelet by analyzing the structures of the amplitude spectrum of a seismic trace. We prove that the operator is contract and has a uniqueness fix point. Then we propose a contract operator method to extract the amplitude spectrum of a seismic wavelet.
We applied our method to different synthetic traces, which are constructed from different wavelet and reflection coefficients having different probability distributions and compared to the results computed from other existing widely used methods. Our method always performs better than other two methods and produces accurate results for all of them.
On the choice of the deterministic and statistical approaches: In a study similar to that of an amplitude-versus-offset (AVO), we need the total wavelet, including the phase and sourcearray pattern. Under this circumstance, the deterministic approach is appropriate. When we process the post-stack seismic data and carry out reservoir parameters inversion based on the convolution model of seismic trace, it is suitable to use the statistical methods to estimate seismic wavelet. If the amplitude spectrum of seismic wavelet is a unimodal function, the proposed approach is a good option.
In the future, we plan to continue our research in the following aspects.
(1) Compare the performance of the proposed method with the commonly used statistical methods for nonstationary seismic trace.
(2) Compare the result of the proposed method with the deterministic one for synthetic data and field data. Then explore the feasibility of applying the proposed method to the inversion of prestack seismic data.