1D Wavelet Transform and Geosciences

In this chapter we present some applications of the wavelet transform in geosciences. The goal is to resolve many crucial problems. A new technique based on the discrete and continuous wavelet transform has been proposed for seismic data denoising. In geomagnetism a wavelet based model for solar geomagnetic disturbances study is established. In petrophysics, we have proposed a new tool of heterogeneities analysis based on the 1D wavelet transform modulus maxima lines (WTMM) method, the proposed tool has been applied on real well-logs data of a borehole located in Algerian Sahara.


Introduction
The one directional Wavelet Transform (WT) is a mathematical tool based on the convolution of a signal with an analyzing wavelet. Despite its simplicity it has been used in various fields of geosciences, in gravity the WT is used for causative sources characterization (Martelet et al, 2001, Ouadfeul et al, 2010. Ouadfeul and Aliouane (2011) have proposed a technique of lithofacies segmentation based on processing of well-logs data by the 1D continuous wavelet transform. Cooper et al (2010) have published a paper focused on the fault determination using one dimensional wavelet analysis. Chamoli(2009) has analyzed the geophysical time series using the wavelet transform.
In seismic data processing the 1D WT has been used by many researchers to denoise the seismic data (Xiaogui Miao andScott Cheadle, 1998, Ouadfeul, 2007).
In this chapter we present some applications of the wavelet transform in geosciences. The goal is to resolve many crucial problems. A new technique based on the discrete and continuous wavelet transform has been proposed for seismic data denoising. In geomagnetism a wavelet based model for solar geomagnetic disturbances study is established. In petrophysics, we have proposed a new tool of heterogeneities analysis based on the 1D wavelet transform modulus maxima lines (WTMM) method, the proposed tool has been applied on real well-logs data of a borehole located in Algerian Sahara.

Random seismic noise attenuation using the wavelet transform
Noise attenuation is a very important task in the seismic data processing field. One can distinguish two types of noises, coherent and incoherent. For the coherent noise we use usually the F-K filter, the deconvolution, The Radon transform…etc. For attenuation of incoherent or random noise, the stack of CDP gathers are one of the seismic data processing steps that improve the S/N (Signal to Noise ratio). We can use a band-pass filter before or after stack to attenuate the random noises. Usually, the Butterworth band-pass filter is used to attenuate this type of noise.
The wavelet transform has becoming a very useful tool in the noise attenuation from seismic data. Prasad (2006) has proposed a technique of ground-roll attenuation from seismic data using the wavelet transform, Xiaogui Miao et al (1998) have published a technique of ground-roll, guided waves, swell noise and random noise attenuation using the discrete wavelet transform.
Siyuan Cao and Xiangpeng Chen (2005) have used the second-generation wavelet transform for random noise attenuation. Ouadfeul (2007) has proposed a technique of random noise attenuation based on the fractal analysis of the seismic data; this technique shows robustness for attenuation of random noises. In this section, we propose a new technique of random noises attenuation from seismic data using the discrete and the continuous wavelet transforms, we start by describing the principles of the continuous and discrete wavelet transforms, after that the processing algorithm of the proposed technique has been detailed. The next step consists to apply this technique at a randomized synthetic seismogram. The proposed technique has been used to denoise a VSP seismic seismogram realized in Algeria. We finalize by the results interpretation and conclusion.

The continuous wavelet transform (CWT)
Here we review some of the important properties of wavelets, without any attempt at being complete. What makes this transform special is that the set of basis functions, known as wavelets, are chosen to be well-localized (have compact support) both in space and frequency (Arnéodo et al., 1988;Arnéodo et al., 1995). Thus, one has some kind of "duallocalization" of the wavelets. This contrasts the situation met for the Fourier's transform where one only has "mono-localization", meaning that localization in both position and frequency simultaneously is not possible.
The CWT of a function s(z) is given by Grossmann and Morlet, (1985) as: Each family test function is derived from a single function () z  defined to as the analyzing wavelet according to (Torresiani, 1995): Where a R   is a scale parameter, b R  is the translation and  * is the complex conjugate of . The analyzing function () z  is generally chosen to be well localized in space (or time) and wavenumber. Usually, ψ(z) is only required to be of zero mean, but for the particular purpose of multiscale analysis ψ(z) is also required to be orthogonal to some low order polynomials, up to the degree n−1, i.e., to have n vanishing moments : According to equation (3), p order moment of the wavelet coefficients at scale a reproduce the scaling properties of the processes. Thus, while filtering out the trends, the wavelet transform reveals the local characteristics of a signal, and more precisely its singularities.
It can be shown that the wavelet transform can reveal the local characteristics of s at a point z 0 . More precisely, we have the following power-law relation (Hermann, 1997;Audit et al., 2002): Where h is the Hölder exponent (or singularity strength). The Hölder exponent can be understood as a global indicator of the local differentiability of a function s.
The scaling parameter (the so-called Hurst exponent) estimated when analysing process by using Fourier's transform (Ouadfeul and Aliouane, 2011) is a global measure of self-affine process, while the singularity strength h can be considered as a local version (i.e. it describes 'local similarities') of the Hurst exponent. In the case of monofractal signals, which are characterized by the same singularity strength everywhere (h(z) = constant), the Hurst exponent equals h. Depending on the value of h, the input signal could be long-range correlated (h > 0.5), uncorrelated (h = 0.5) or anticorrelated (h < 0.5).

The discrete wavelet transform (DWT)
L 2 (R) denotes the Hibert space of measurable, square-integrable functions. The function is said to be a wavelet if and only when the following condition is satisfied.
is the dilation of () t  by the scale factor s.
In order to be used expediently in practice a, is scattered as discrete binary system ,i.e. Let ,its wavelet transform is : Hence its contrary transform is Being dispersed in time domain farther, a discrete wavelet transform can be obtained. It exists an effective and fast algorithm, it is based on equation (7) 1 1 1 22 There 2 j Wf is the wavelet transform coefficients of f(t). It approximates f(t) on the scale 2 j .
H j and G j are the discrete filters gained by inserting (2 j -1) zeros into every two samples of H, G. And the relation between G and H is:

Signal denoising
Thresholding is a technique used for signal and image denoising. The discrete wavelet transform uses two types of filters: (1) averaging filters, and (2) detail filters. When we decompose a signal using the wavelet transform, we are left with a set of wavelet coefficients that correlate to the high frequency sub-bands. These high frequency subbands consist of the details in the data set. If these details are small enough, they might be omitted without substantially affecting the main features of the data set. Additionally, these small details are often those associated with noise; therefore, by setting these coefficients to zero, we are essentially killing the noise. This becomes the basic concept behind thresholding-set all frequency sub-band coefficients that are less than a particular threshold to zero and use these coefficients in an inverse wavelet transformation to reconstruct the data set.

The denoising algorithm
The denoising algorithm is based on the discrete wavelet transform decomposition combined with the continuous wavelet transform. Firstly, discrete wavelet decomposition has been made; the analyzing wavelet is the Haar of level 5 (Charles, 1992). After that we apply a threshold to denoise the seismic trace.
The next step consists to calculate the continuous wavelet transform of the denoised trace obtained by DWT. The final denoised seismic seismogram is the wavelet coefficients at the scale a=2.*∆T. Where: ∆T is the sampling interval.
The analyzing wavelet is the modified Mexican Hat, it is defined by equation 9 (See figure1). The flow chart of the proposed processing algorithm is detailed in figure2.

Application on synthetic data
The proposed technique has been applied at the synthetic seismogram of a geological model with the parameters detailed in table1. The synthetic seismogram is generated with a sampling interval of 2ms.The full recording time is 2.5s. Figure 3 is a presentation of the noisy seismogram versus the time with a 200% of white noise. The discrete wavelet decomposition is presented from level 1 to 5 in figure 4. The denoised seismic seismogram is presented in figure5. One can remark that the major high frequency fluctuations are eliminated by this last operation. The next step consists to calculate the wavelet coefficients; the analyzing wavelet is the modified Mexican Hat . The wavelet coefficients versus the time and scale are presented in figure6.

Discrete wavelet transform of the seismic seismogram T(t)
Densoing of T(t) using the threshold method Td(t) is the denoised seismic seismogram Calculation of the continuous wavelet transform CWT(t,a) of Td(t) The analyzing wavelet is the modified Mexican Hat.
Plotting of the CWT(t,a) at the scale a=2*∆T ∆T : is the sampling interval www.intechopen.com The final denoised seismic seismogram is presented in figure 7a. It is clear that the proposed technique is able to attenuate the random noises from the synthetic seismogram. Figure 7b presents the residual noises in the frequency domain. The spectral analysis of the detrended noise is presented in figure 7b, this last represents the amplitude and the phase spectrums. These spectrums are calculated using the Fourier's transform. Analysis of this figure shows that the residual noise containing both the low and the high frequencies. This is the characteristic of the white noise. The phase spectrum shows that this noise contains all the angles [- , +  ]. The next operation consists to apply the sketched method at real data.

Application on real data
The proposed technique has been tested at a raw seismogram of a vertical seismic profile (VSP) realized in Algeria. Figure 8 is a representation of this seismic seismogram versus the time, the sampling interval is 0.002s. The recording interval is 2.048s. The discrete wavelet decomposition using the Haar wavelet of level 5 is presented in figure 9. The denoised seismic seismogram is showed in figure 10. One can remark that many random types of amplitude have been eliminated by this operation. The last procedure of the processing algorithm consists to calculate the continuous wavelet transform coefficients. This last is presented in figure11. The final denoised seismic trace is the wavelet coefficients at the scale a=0.004s. This last is presented in figure12 versus the time.

Results Interpretation and conclusion
Spectral analysis of the residual noise using the Fourier's transform (See figure 13) shows that the residual noise contains the frequency band [100Hz,150Hz]. Note that the spatial filter during acquisition can attenuate some noises. The phase spectrum shows that this last sweep the full interval [-π,+π].
We have proposed a new technique of random noise attenuation based on the threshold method using the discrete and the continuous wavelet transform. Application on noisy synthetic seismic seismogram shows the robustness of the proposed tool. However the proposed tool doesn't preserve amplitude, to resolve this problem we recommend to apply a gain at the final seismic trace, this last is derived from the raw seismic trace. We suggest integrating this technique of seismic noise attenuation using the wavelet transform in the seismic data processing software's.

Solar geomagnetic disturbances analysis using the CWT
In this part, we use the wavelet transform modulus maxima lines WTMM method as a tool to schedule the geomagnetic disturbances. The proposed idea is based on the estimation of the singularity strength (Hölder exponent) at maxima of the modulus of the 1D continuous wavelet transform of the Horizontal component of the geomagnetic field. Data of International R eal-time Magnetic Observatory Network (InterMagnet) observatories are used.

Fractal analysis of geomagnetic disturbances using the CWT
In this section, we use the continuous wavelet transform as a tool for analyzing the horizontal geomagnetic field component of the InterMagnet observatories. The goal is to to schedule the solar geomagnetic disturbances. The proposed technique is based on the calculation of the modulus of the continuous wavelet transform, after that Hölder exponents are estimated on the maxima of the CWT. The analyzing wavelet is the Complex Morlet (See Ouadfeul and Aliouane, 2011). The flow chart of the proposed technique is detailed in figure  14.  Figure 16 shows the modulus of the CWT, the analyzing wavelet is the Complex Morlet. Estimated Hölder exponents at maxima of the modulus of CWT are presented in figure 17. Figure 18 is the average Hölder exponents at each one hour of time compared with the normalized DST index. One can remark that the Hölder exponent estimated by the CWT is a very robust tool for scheduling solar geomagnetic disturbances. It can be used as an index for solar geomagnetic disturbances schedule.
Same analysis has been applied at the geomagnetic data of Backer Lake, Kakioka, Hermanus and Alibag observatories. A detailed analysis shows that before the magnetic storm we observe a decrease of the Hölder exponent. In the moment of the solar disturbance the Hölder exponent has a very low value ( 0) h  .

Generalized fractal dimension and geomagnetic disturbances
A Generalized Fractal Dimension (GFD) based on the spectrum of exponent calculated using the wavelet transform modulus maxima lines (See Arneodo et al, 1995) method has been used for geomagnetic disturbances schedule. The GFD is calculated using equation (9).
is the spectrum of exponent estimated using the function of partition (See Arneodo et al, 1955.

Application on the Hermanus observatory data
In this section we analyze the total magnetic field of the Hermanus observatory for the month of May 2002. The proposed idea consists to estimate the fractal dimensions D 0 , D 1 and D 2 for every 60 minutes (one hour of the month). Obtained results are presented in figure 19.  Analysis of obtained results shows that the fractal dimension D 0 is not sensitive to magnetic disturbances. However D 1 and D 2 are very sensitive to the solar geomagnetic activity.
One can remark that the major magnetic disturbances are characterized by spikes ( See table 2 and Figure 19).

Heterogeneities analysis using the 1D wavelet transform modulus maxima lines
Here we use a wavelet transform based multifractal analysis, called the wavelet transform modulus maxima lines (WTMM) method. For more information about the WTMM, author can read the book of Arneodo et al(1995) or the paper of Ouadfeul and Aliouane(2011).
The proposed technique is based on the estimation of the Hölder exponents or roughness coefficient at maxima of the modulus of the CWT. The roughness coefficient is related to rock's heterogeneities (Ouadfeul, 2011, Ouadfeul andAliouane, 2011). Estimation of Hölder exponents is based on the continuous wavelet transform, in fact for low scales the Hölder exponent is related to the modulus of the continuous wavelet transform by (Hermann, 1997;Audit et al., 2002) :

Application on real data
The proposed idea has been applied on the natural gamma ray (Gr) log of Sif-Fatima2 borehole located on the Berkine basin. The goal is the segmentation of the intercalation of the sandstone and clay. The main reservoir where the data are recorded is the Trias-Argilo-Gréseux inférieur (TAGI), this last is composed mainly of the four lithofacies units, which are : the Clay , The sandstone, the Sandy clay and the clayey sandstone.

Geological setting of the Berkine Basin
The Berkine basin is a vast circular Palaeozoic depression, where the basement is situated at more than 7000 m in depth. Hercynian erosion slightly affected this depression because only Carboniferous and the Devonian are affected at their borders. The Mesozoic overburden varied from 2000m in Southeast to 3200m in the Northeast. This depression is an intracratonic basin which has preserved a sedimentary fill out of more than 6000 m. It is characterized by a complete section of Palaeozoic formations spanning from the Cambrian to the Upper Carboniferous. The Mesozoic to Cenozoic buried very important volume sedimentary material contained in this basin presents an opportunity for hydrocarbons accumulations. The Triassic province is the geological target of this study. It is mainly composed by the Clay and Sandstone deposits. Its thickness can reach up to 230m. The Sandstone deposits constitute very important hydrocarbon reservoirs.
The SIF FATIMA area where the borehole data are collected is restricted in the the labelled 402b block. It is located in the central part of the Berkine basin (Fig.20). The hydrocarbon field is situated in the eastern erg of the basin characterized also by high amplitude topography.
The studied area contains many drillings. However this paper will be focused on the Sif-Fatima2 borehole. The main reservoir, the Lower Triassic Clay Sandstone labelled TAGI , is represented by fluvial and eolian deposits. The TAGI reservoir is characterized by three main levels: Upper , middle , and lower. Each level is subdivided into a total of nine subunits according to SONATRACH nomenclature (Zeroug et al., 2007).The lower TAGI is often of a very small thickness. It is predominantly marked by clay facies, sometimes by sandstones and alternatively by the clay and sandstone intercalations, with poor petrophysical characteristics.

Data processing
Fluctuation of the gamma ray log are presented in figure 21, the modulus of the CWT of this log is presented in the plan depth-log the scale in figure 22. The next step consists to calculate maxima of the modulus of the CWT at the set of scales. Scales are varied from 0.5m to 246m, the dilatation method used is power of 2. The set of maxima mapped for all scales is called the skeleton of the CWT, this last is presented in figure 23.   Table 3. Lithofacies intervals derived from the GR signal for Sif-Fatima 2 well

Results Interpretation
A preliminary raw lithofacies classification based on the natural gamma radioactivity welllog data was made. First, recall that the maximum value GR max of the data is considered as a full clay concentration while the minimum value GR min represents the full sandstone concentration. The mean value (GR max+ GR min )/2 will then represent the threshold that will be used as a decision factor within the interval studied: -Geological formations bearing a natural GR activity characterized by: GR Threshold < GR <GR max Are considered as Sandy Clay.
-Geological formations with a natural GR activity characterized by:

GR min < GR <GR Threshold
Are considered as a Clayey Sandstone.
The results for Sif Fatima2 borehole are illustrated in figure 24 that shed light on the obtained segmentation. Moreover, the depth distribution of the different facies is given in table3 .
Estimated Hölder exponents at maxima of the modulus of the CWT compared with the classical segmentation based on the gamma-ray log (figure 24). One can remark that the Hölder exponent can be used as an attribute to seek the fines lithofacies.

Conclusion
We have used in this chapter the 1D discrete and continuous wavelet transform to resolve many problems in geosciences. The DWT in combination with the CWT prove that they can be used as tool for seismic data denoising. The continuous wavelet transform can be used for fractal and multifractal analysis of geomagnetic data, the goal is to schedule the solar geomagnetic disturbances. Obtained results show the robustness of the CWT.
We have proposed a new tool based on the wavelet transform modulus maxima lines (WTMM) method for lithofacies segmentation from well-logs data. Comparison with the classical method of segmentation based on the gamma ray log shows that the fractal analysis revisited by the continuous wavelet transform can provide geological details and intercalations.