**Part 2**

**Applications of the Wavelet Transforms in Geoscience** 

26 Will-be-set-by-IN-TECH

250 Wavelet Transforms and Their Recent Applications in Biology and Geoscience

M C Robini, I E Magnin, H Cattin, and A Baskurt. Two-dimensional ultrasonic flaw

R Haralick, K Shunmugam, and I Dinstein. Textural features for image classification, IEEE Transactions on Systems, Man Cybernetics, SMC, Vol 3, No 1, pp 610-621, 1973. F Wong, R Nagarajan, S Yaacob, A Chekima, and N. E. Belkhamza. An image segmentation

K M Rajpoot andNMRajpoot. Wavelets and support vector machines for texture

A Al-Ataby, W Al-Nuaimy, C R Brett and O Zahran. Automatic detection and classification

British Standards Institution, 'Guide to calibration and setting-up of the ultrasonic time of

on Image Processing, Vol 17, No 9, pp 1709-1720, September 2008.

Ferroelectrics, and Frequency Control, Vol 44, 1997.

pp 328-333, December 2004.

7706, December 1993.

and Its Applications (ISSPA), Vol 1, pp 144-147, August 2001. V N Vapnik. The Nature of Statistical Learning Theory, Springer, 2nd edition, 2000.

detection based on the Wavelet packet transform, IEEE transactions on Ultrasonics,

method using fuzzy-based threshold, International Symposium on Signal Processing

classification, Proceedings of INMIC 2004, 8th International Multitopic Conference,

of weld flaws in TOFD data using wavelet transform and support vector machines, Insight, Journal of the British Institute of NDT, Vol 52, No 11, November 2010. K Huang and S Aviyente. Wavelet feature selection for image classification, IEEE Transaction

flight diffraction (TOFD) technique for the detection, location and sizing of flaws', BS

**13** 

*Algeria* 

**Multiscale Analysis of Geophysical Signals Using the 2D Continuous Wavelet Transform** 

The continuous wavelet transform has becoming a very useful tool in geophysics (Kumar and Foufoula-Georgiou, 1997, Ouadfeul, 2006, Ouadfeul, 2007, Ouadfeul, 2008, Ouadfeul and Aliouane, 2010, Ouadfeul et al, 2010). In Potential field analysis it was used to locate

Martelet et al (2001) have published a paper on the characterization of geological boundaries using the 1D wavelet transform of gravity data, the proposed technique has been applied on the

The Complex continuous wavelet transform has been used for identification of sources of potential fields using the aeromagnetic profiles of the French Guiana (Sailhac et al, 2000).

Sailhac and Gibert (2003) have proposed a new technique of sources identification form potential field with the continuous wavelet transform, it is based on the two-dimensional

A new technique has been proposed by Vallée et al (2004) to estimate depth and model type

Boukerbout et al (2006) have applied the continuous wavelet transform formalism to the special case of anomalies produced by elongated sources like faults and dikes. They show that, for this particular type of anomalies, the two-dimensional wavelet transform corresponds to the ridgelet analysis and reduces to the 1D wavelet transform applied in the Radon domain. Recently Ouadfeul et al (2010), have published a paper that use the continuous wavelet transform on the mapping of geological contacts from geomagnetic data, the analyzing wavelet is the Mexican Hat, the proposed method has shown a good robustness. In this paper we propose a new technique more precise; this last is based on the 2D directional continuous wavelet transform. The analyzing wavelet is the Poisson's Kernel (Martelet et al, 2001). We start firstly by giving the theoretical basis of the 2D directional Wavelet Transform (DCWT) and the analogy between the wavelet transform and the upward continuation, after that the technique has been applied on a synthetic data of prism and a cylinder. The proposed idea has been

and characterize homogeneous causative sources point in 1D (Moreau et al, 1997).

**1. Introduction** 

Himalaya.

wavelets and multipolar approximations.

using the continuous wavelet transform of magnetic data.

Mohamed Hamoudi2, Amar Boudella2 and Said Eladj3 *1Geosciences and Mines, Algerian Petroleum Institute, IAP* 

Sid-Ali Ouadfeul1,2, Leila Aliouane2,3,

*2Geophysics Department, FSTGAT, USTHB* 

*3Geophysics Department, LABOPHYT, FHC, UMBB* 

### **Multiscale Analysis of Geophysical Signals Using the 2D Continuous Wavelet Transform**

Sid-Ali Ouadfeul1,2, Leila Aliouane2,3, Mohamed Hamoudi2, Amar Boudella2 and Said Eladj3 *1Geosciences and Mines, Algerian Petroleum Institute, IAP 2Geophysics Department, FSTGAT, USTHB 3Geophysics Department, LABOPHYT, FHC, UMBB Algeria* 

#### **1. Introduction**

The continuous wavelet transform has becoming a very useful tool in geophysics (Kumar and Foufoula-Georgiou, 1997, Ouadfeul, 2006, Ouadfeul, 2007, Ouadfeul, 2008, Ouadfeul and Aliouane, 2010, Ouadfeul et al, 2010). In Potential field analysis it was used to locate and characterize homogeneous causative sources point in 1D (Moreau et al, 1997).

Martelet et al (2001) have published a paper on the characterization of geological boundaries using the 1D wavelet transform of gravity data, the proposed technique has been applied on the Himalaya.

The Complex continuous wavelet transform has been used for identification of sources of potential fields using the aeromagnetic profiles of the French Guiana (Sailhac et al, 2000).

Sailhac and Gibert (2003) have proposed a new technique of sources identification form potential field with the continuous wavelet transform, it is based on the two-dimensional wavelets and multipolar approximations.

A new technique has been proposed by Vallée et al (2004) to estimate depth and model type using the continuous wavelet transform of magnetic data.

Boukerbout et al (2006) have applied the continuous wavelet transform formalism to the special case of anomalies produced by elongated sources like faults and dikes. They show that, for this particular type of anomalies, the two-dimensional wavelet transform corresponds to the ridgelet analysis and reduces to the 1D wavelet transform applied in the Radon domain.

Recently Ouadfeul et al (2010), have published a paper that use the continuous wavelet transform on the mapping of geological contacts from geomagnetic data, the analyzing wavelet is the Mexican Hat, the proposed method has shown a good robustness. In this paper we propose a new technique more precise; this last is based on the 2D directional continuous wavelet transform. The analyzing wavelet is the Poisson's Kernel (Martelet et al, 2001). We start firstly by giving the theoretical basis of the 2D directional Wavelet Transform (DCWT) and the analogy between the wavelet transform and the upward continuation, after that the technique has been applied on a synthetic data of prism and a cylinder. The proposed idea has been

Multiscale Analysis of Geophysical Signals Using the 2D Continuous Wavelet Transform 255

potential field F at a scale *a* is equivalent to the Upward continue of this field at the level Z=a. Positioning of maxima of the modulus of the continuous wavelet transform at this scale is equivalent to the maxima of the horizontal gradient of the potential field upward at the

> 3/2 2 2 1 1 , <sup>2</sup> <sup>1</sup>

The sharp contrasts that show the potential data are assumed to result from discontinuities or interfaces such as faults, flexures, contrasts intrusive rocks ... For contacts analysis between geological structures, we use usually the classical methods based on the location of local maxima of the modulus of the total (Nabighian, 1984) or the horizontal gradient (Blakely et al., 1986), or the Euler's deconvolution (Reid et al., 1990). This technique allows, in addition to localization in the horizontal plane of contact, an estimate of their depth. The potential field, over a vertical contact, involving the presence of rocks of different susceptibilities is indicated by a low in side rocks of low susceptibility and a high in side rocks of high susceptibility. The inflection point is found directly below the vertical contact. We can use this characteristic of geomagnetic anomalous for localization of abrupt susceptibility change. If the contact has a dip, the maxima of horizontal gradients move in the direction of dip. To determine the dip direction of contacts, we upward the map of the potential field at different altitudes. At each level, the maxima of horizontal gradient are located. If the structures are vertical, all maxima are superposed. However, moving of maxima with the upward indicates the direction of the dip. The potential theory lends

By choosing an appropriate wavelet, measurement of geomagnetic field or its spatial derivatives can be processed as a wavelet transform. Indeed, this analysis unifies various classical techniques: it process gradients that have been upward to a range of altitudes. The expressions of various conventional operations on the potential field are well-designed in the wavelet domain. The most important is the equivalence between the concept of scaling and the upward. Indeed, the wavelet transform of a potential field F0 (x, y) at a certain scale

For a multiscale analysis of contacts, it is sufficient to look for local maxima of the modulus of the continuous wavelet transform (CWT) for different scales to get an exact information

The proposed idea has been applied on synthetic geomagnetic response of a cylinder and a prism, after that the technique is applied on real geomagnetic data of In Ouzzal area, located

*x y*  (4)

*P xy*

**3. Wavelet transform and potential data** 

perfectly to a multiscale analysis by wavelet transform.

1. Upward continue the measured filed at level Z=a\*Z0 2. Calculation of the horizontal gradient in the plane (x,y).

about geological boundaries (Ouadfeul et al, 2010).

on Algerian Sahara. Lets us starting by the synthetic example.

**4. Multiscale analysis geomagnetic data** 

3. Multiplication by a.

a = Z/Z0 can be obtained from measurements made on the level Z0 by:

level Z.

applied at the anomaly magnetic field of In Ouzzal area. The obtained results after applying the proposed filter have been compared with analytic signal solutions (Nabighian, 1984; Roset et al, 1992). We finalize the part by propping a model of contacts for this area and a conclusion.

The second part of this chapter consists to apply a multiscale analysis at the 3D analytic signal of the real aeromagnetic data using the wavelet transform, the proposed technique shows a robustness, where classical methods are fails. The third part is an application of the wavelet transform at the gravity data of an area located in the Algerian Sahara.

The last part consists to propose a new method to reduce the noise effect when analyzing the 3D GPR data by the wavelet transform. We finalize this chapter by a conclusion that resumes the different applications of the conventional and directional wavelet transform and its benefits of the different geophysical data.

#### **2. The 2D directional continuous wavelet transform**

The 2D directional continuous wavelet transform has been introduced firstly by Murenzi (1989). The wavelet decomposition of a given function 2 2 *f L R*( ) with an analyzing wavelet 2 2 *g L R*( ) defined for all a>0 , <sup>2</sup> *b R* , 2 2 *f LR* ( ), 0,2 by:

$$\mathcal{W}\_{\mathcal{G}}f(a,b,\alpha) = \iint\_{\mathbb{R}^2} f(\mathbf{x}) \frac{1}{a} \mathbf{g}(r\_{-\alpha}(\frac{\mathbf{x}-b}{a})) d\mathbf{x} \tag{1}$$

Where *r* is a rotation with an angle in the <sup>2</sup> *R* .

An example of a directional wavelet transform is the gradient of the Gaussian Wavelet *G* . Because the convolution of an image with *G* is equivalent to analysis of the gradient of the modulus of the continuous wavelet transform. Canny has introduced another tool for edge detection, after the convolution of the image with a Gaussian, we calculate its gradient to seek the set of points that corresponding to the high variation of the intensity of the 2D continuous wavelet transform. The use of the gradient of the Gaussian as an analyzing wavelet has been introduced by Mallat and Hwang (1992).

The wavelet transform of a function f, with an analyzing wavelet *g G* is a vector defined for all a>0, 2 2 *f L R*( ) <sup>2</sup> *b R* by:

$$\mathcal{W}\_{\nabla G} f(a, b, \alpha) = \iint\_{\mathbb{R}^2} f(\mathbf{x}) \frac{1}{a} \nabla G(\frac{\mathbf{x} - b}{a}) d\mathbf{x} \tag{2}$$

This transformation relates the directional wavelet transform if we choose *<sup>G</sup> g x* . In fact, we have the following relation:

$$\underbrace{\mathcal{W}\_{\text{\tiny\text{\tiny\text{\tiny\text{\tiny\text{\tiny\text{\tiny\text{\tiny\text{\tiny\text{\tiny\text{\tiny\text{\tiny\text{\tiny\text{\tiny\text{\text}}}}}}}}}}}}\_{\text{\text{\text{\tiny\text{\textquotedblleft}}}}f}f(a,b,b,a) = \mathfrak{u}\_{\alpha}\mathcal{W}\_{\text{\textquotedblleft}\text{\textquotedblright}}f(a,b)}f(a,b)\tag{3}$$

Where *u* is the unite vector in the direction : (cos( ),sin( )) *u* , and "." Is the Euclidian scalar product in <sup>2</sup> *R* .

In this paper we use a new type of directional wavelet based on the Poisson's Kerenel defined in equation (4) as a filer. The modulus of the directional wavelet transform of a potential field F at a scale *a* is equivalent to the Upward continue of this field at the level Z=a. Positioning of maxima of the modulus of the continuous wavelet transform at this scale is equivalent to the maxima of the horizontal gradient of the potential field upward at the level Z.

$$P\left(\mathbf{x}, y\right) = \frac{1}{2\Pi} \frac{1}{\left[\mathbf{x}^2 + y^2 + 1\right]^{3/2}}\tag{4}$$

#### **3. Wavelet transform and potential data**

254 Wavelet Transforms and Their Recent Applications in Biology and Geoscience

applied at the anomaly magnetic field of In Ouzzal area. The obtained results after applying the proposed filter have been compared with analytic signal solutions (Nabighian, 1984; Roset et al, 1992). We finalize the part by propping a model of contacts for this area and a conclusion.

The second part of this chapter consists to apply a multiscale analysis at the 3D analytic signal of the real aeromagnetic data using the wavelet transform, the proposed technique shows a robustness, where classical methods are fails. The third part is an application of the

The last part consists to propose a new method to reduce the noise effect when analyzing the 3D GPR data by the wavelet transform. We finalize this chapter by a conclusion that resumes the different applications of the conventional and directional wavelet transform

The 2D directional continuous wavelet transform has been introduced firstly by Murenzi (1989). The wavelet decomposition of a given function 2 2 *f L R*( ) with an analyzing

*x b <sup>W</sup> <sup>f</sup> a b <sup>f</sup> <sup>x</sup> <sup>g</sup> r dx*

 in the <sup>2</sup> *R* . An example of a directional wavelet transform is the gradient of the Gaussian Wavelet *G* . Because the convolution of an image with *G* is equivalent to analysis of the gradient of the modulus of the continuous wavelet transform. Canny has introduced another tool for edge detection, after the convolution of the image with a Gaussian, we calculate its gradient to seek the set of points that corresponding to the high variation of the intensity of the 2D continuous wavelet transform. The use of the gradient of the Gaussian as an analyzing

The wavelet transform of a function f, with an analyzing wavelet *g G* is a vector defined

2 <sup>1</sup> ( , , ) ( ) ( )) *<sup>G</sup> R <sup>x</sup> <sup>b</sup> W f a b f x G dx*

> *W f ab u W f a b*

In this paper we use a new type of directional wavelet based on the Poisson's Kerenel defined in equation (4) as a filer. The modulus of the directional wavelet transform of a

This transformation relates the directional wavelet transform if we choose *<sup>G</sup>*

2 <sup>1</sup> ( , , ) ( ) ( ( )) *<sup>g</sup>*

*R*

*a a*

*a a*

: (cos( ),sin( )) *u*

  by:

(1)

(2)

(3)

*g*

, and "." Is the Euclidian

*x* 

. In fact,

wavelet transform at the gravity data of an area located in the Algerian Sahara.

and its benefits of the different geophysical data.

is a rotation with an angle

wavelet has been introduced by Mallat and Hwang (1992).

(,, ) (,) *<sup>G</sup> <sup>G</sup>*

is the unite vector in the direction

*x*

Where *r*

Where *u*

scalar product in <sup>2</sup> *R* .

for all a>0, 2 2 *f L R*( ) <sup>2</sup> *b R* by:

we have the following relation:

**2. The 2D directional continuous wavelet transform** 

wavelet 2 2 *g L R*( ) defined for all a>0 , <sup>2</sup> *b R* , 2 2 *f LR* ( ), 0,2

The sharp contrasts that show the potential data are assumed to result from discontinuities or interfaces such as faults, flexures, contrasts intrusive rocks ... For contacts analysis between geological structures, we use usually the classical methods based on the location of local maxima of the modulus of the total (Nabighian, 1984) or the horizontal gradient (Blakely et al., 1986), or the Euler's deconvolution (Reid et al., 1990). This technique allows, in addition to localization in the horizontal plane of contact, an estimate of their depth. The potential field, over a vertical contact, involving the presence of rocks of different susceptibilities is indicated by a low in side rocks of low susceptibility and a high in side rocks of high susceptibility. The inflection point is found directly below the vertical contact. We can use this characteristic of geomagnetic anomalous for localization of abrupt susceptibility change. If the contact has a dip, the maxima of horizontal gradients move in the direction of dip. To determine the dip direction of contacts, we upward the map of the potential field at different altitudes. At each level, the maxima of horizontal gradient are located. If the structures are vertical, all maxima are superposed. However, moving of maxima with the upward indicates the direction of the dip. The potential theory lends perfectly to a multiscale analysis by wavelet transform.

By choosing an appropriate wavelet, measurement of geomagnetic field or its spatial derivatives can be processed as a wavelet transform. Indeed, this analysis unifies various classical techniques: it process gradients that have been upward to a range of altitudes. The expressions of various conventional operations on the potential field are well-designed in the wavelet domain. The most important is the equivalence between the concept of scaling and the upward. Indeed, the wavelet transform of a potential field F0 (x, y) at a certain scale a = Z/Z0 can be obtained from measurements made on the level Z0 by:


For a multiscale analysis of contacts, it is sufficient to look for local maxima of the modulus of the continuous wavelet transform (CWT) for different scales to get an exact information about geological boundaries (Ouadfeul et al, 2010).

#### **4. Multiscale analysis geomagnetic data**

The proposed idea has been applied on synthetic geomagnetic response of a cylinder and a prism, after that the technique is applied on real geomagnetic data of In Ouzzal area, located on Algerian Sahara. Lets us starting by the synthetic example.

Multiscale Analysis of Geophysical Signals Using the 2D Continuous Wavelet Transform 257

Y(m) X(m)

Fig. 1c. Modulus of the continuous wavelet transform of the Prism and the Cylinder plotted

X(m)

Fig. 1b. Anomaly field of the synthetic model of figue1

F(nT)

at the low Scale a=282m.

Y(m)

C(X,Y)

#### **4.1 Application on synthetic geomagnetic data**

The proposed idea has been applied at a synthetic model of a cylinder and prism, parameters of these last are resumed in table1 (See also figure1a). Figure 1b is the magnetic response of this model generated with a grid of dimensions 100mx100m. The first operation is to calculate the modulus of the continuous wavelet transform. The analyzing wavelet is the Poisson's Kernel defined by equation 4.

Figure 1c shows this modulus plotted at the scale a=282m. The second step consists to calculate its maxima, figure 2 is a map of these maxima for all scales (Scales are varied from 282m to 9094m). Solid curves are the exact boundaries of the prism and the cylinder. One can remark that the maxima of the continuous wavelet transform are positioned around the two exact boundaries.


Table 1a. Physical parameters of the Cylinder


Table 1b. Physical parameters of the Prism

Fig. 1a. Physical parameters of a synthetic model composed of Prism and Cylinder.

The proposed idea has been applied at a synthetic model of a cylinder and prism, parameters of these last are resumed in table1 (See also figure1a). Figure 1b is the magnetic response of this model generated with a grid of dimensions 100mx100m. The first operation is to calculate the modulus of the continuous wavelet transform. The analyzing wavelet is

Figure 1c shows this modulus plotted at the scale a=282m. The second step consists to calculate its maxima, figure 2 is a map of these maxima for all scales (Scales are varied from 282m to 9094m). Solid curves are the exact boundaries of the prism and the cylinder. One can remark that the maxima of the continuous wavelet transform are positioned around the

**4.1 Application on synthetic geomagnetic data** 

coordinates of the center (5000, 2500, -250).

coordinates of the center (5000, 7000, -300).

Fig. 1a. Physical parameters of a synthetic model composed of Prism and Cylinder.

2000m

2500m

3000m 3000m 2000m

Ray 1500m. High 2500m. Magnetic Susceptibility K=0.015 SI. F 37000 nT Declination D=0° Inclination I=90°

Width 3000m Length 3000m. High 2000m. Magnetic Susceptibility K=0.01 SI. F 37000 nT Declination D=0° Inclination I=90°

Table 1a. Physical parameters of the Cylinder

Table 1b. Physical parameters of the Prism

the Poisson's Kernel defined by equation 4.

two exact boundaries.

Fig. 1b. Anomaly field of the synthetic model of figue1

Fig. 1c. Modulus of the continuous wavelet transform of the Prism and the Cylinder plotted at the low Scale a=282m.

Multiscale Analysis of Geophysical Signals Using the 2D Continuous Wavelet Transform 259

Fig. 3. Map of the principal subdivisions and principal structural domains of Hoggar, afer

The In Ouzzal terrane (Western Hoggar) is an example of Archaean crust remobilized by a very-high-temperature metamorphism (Ouzeggane and Boumaza, 1996) during the Paleoproterozoic (2 Ga). Structural geometry of the In Ouzzal terrane is characterized by closed structures trending NE-SW to ENE-WSW(figure4) that correspond to domes of charnockitic orthogneiss. The supracrustal series are made up of metasediments and basicultrabasic rocks that occupy the basins located between these domes. In In-Ouzzal area, the supracrustal synforms and orthogneiss domes exhibit linear corridors near their contacts corresponding to shear zones. The structural features in In Ouzzal area (Djemaï et al, 2009), observed at the level of the base of the crust, argue in favour of a deformation taking place entirely under granulite-facies conditions during the Paleoproterozoic. These features are compatible with D1 homogeneous horizontal shortening of overall NW-SE trend that accentuates the vertical stretching and flattening of old structures in the form of basins and domes. This shortening was accommodated by horizontal displacements along transpressive shear corridors. During the Pan-African event, the brittle deformation affected the granulites which were retrogressed amphibolite and greenschists facies (with the development of tremolite and chlorite (Caby and Monie, 2003), in the presence of fluids along shear zones corridors. Brittle deformations were concentrated in the southern boundary of In Ouzzal. An important NW-SE-trending dextral strike-slip pattern has been mapped along which we can see the Eburnean foliation F1 overprinted. This period was also

Caby et al(1981)

**4.2.2 Geological context of In Ouzzal** 

Fig. 2. Mapped contacts obtained by positioning of maxima of the modulus of the 2D DCWT. The solid curves are the exact boundaries of the Prism and the Cylinder.

#### **4.2 Application on real data**

The proposed idea has been applied to the aeromagnetic data of In Ouzzal, it is located in Hoggar. We start by describing the geology of the massif of Hoggar.

#### **4.2.1 Geological setting of Hoggar**

The lens of the massif of Hoggar that occupies the southern part of Algeria (Figure 3) is defined as a wide lapel Precambrian. It integrates mobile areas caught between the West African Craton and East Africa. It extends at 1000 km from East to West and 700km from North to South, It is extending through the branches 'append' of Air-NIGER at South-East and the Adrar des Iforas-MALI at the South West. Its territory is covered to the North East and South in part by training Paleozoic sedimentary of Tassili. The entire range (Hoggar, Air and Adrar des Iforas) forms the Tuareg shield (also known as Tuargi shield). Stable since the Cambrian, this massif belongs to the chain called "Pan" which belts the West African Craton to the east. This old chain is probably the result of the evolution of both intercontinental basins and accretion of micro plates involving the creation of oceans and island arcs. It is surrounded by sediment platform of Paleozoic age.

Fig. 3. Map of the principal subdivisions and principal structural domains of Hoggar, afer Caby et al(1981)

#### **4.2.2 Geological context of In Ouzzal**

258 Wavelet Transforms and Their Recent Applications in Biology and Geoscience

Fig. 2. Mapped contacts obtained by positioning of maxima of the modulus of the 2D DCWT. The solid curves are the exact boundaries of the Prism and the Cylinder.

Hoggar. We start by describing the geology of the massif of Hoggar.

island arcs. It is surrounded by sediment platform of Paleozoic age.

The proposed idea has been applied to the aeromagnetic data of In Ouzzal, it is located in

The lens of the massif of Hoggar that occupies the southern part of Algeria (Figure 3) is defined as a wide lapel Precambrian. It integrates mobile areas caught between the West African Craton and East Africa. It extends at 1000 km from East to West and 700km from North to South, It is extending through the branches 'append' of Air-NIGER at South-East and the Adrar des Iforas-MALI at the South West. Its territory is covered to the North East and South in part by training Paleozoic sedimentary of Tassili. The entire range (Hoggar, Air and Adrar des Iforas) forms the Tuareg shield (also known as Tuargi shield). Stable since the Cambrian, this massif belongs to the chain called "Pan" which belts the West African Craton to the east. This old chain is probably the result of the evolution of both intercontinental basins and accretion of micro plates involving the creation of oceans and

**4.2 Application on real data** 

**4.2.1 Geological setting of Hoggar** 

The In Ouzzal terrane (Western Hoggar) is an example of Archaean crust remobilized by a very-high-temperature metamorphism (Ouzeggane and Boumaza, 1996) during the Paleoproterozoic (2 Ga). Structural geometry of the In Ouzzal terrane is characterized by closed structures trending NE-SW to ENE-WSW(figure4) that correspond to domes of charnockitic orthogneiss. The supracrustal series are made up of metasediments and basicultrabasic rocks that occupy the basins located between these domes. In In-Ouzzal area, the supracrustal synforms and orthogneiss domes exhibit linear corridors near their contacts corresponding to shear zones. The structural features in In Ouzzal area (Djemaï et al, 2009), observed at the level of the base of the crust, argue in favour of a deformation taking place entirely under granulite-facies conditions during the Paleoproterozoic. These features are compatible with D1 homogeneous horizontal shortening of overall NW-SE trend that accentuates the vertical stretching and flattening of old structures in the form of basins and domes. This shortening was accommodated by horizontal displacements along transpressive shear corridors. During the Pan-African event, the brittle deformation affected the granulites which were retrogressed amphibolite and greenschists facies (with the development of tremolite and chlorite (Caby and Monie, 2003), in the presence of fluids along shear zones corridors. Brittle deformations were concentrated in the southern boundary of In Ouzzal. An important NW-SE-trending dextral strike-slip pattern has been mapped along which we can see the Eburnean foliation F1 overprinted. This period was also

Multiscale Analysis of Geophysical Signals Using the 2D Continuous Wavelet Transform 261

400000 440000 480000 520000 560000 600000

400000 440000 480000 520000 560000 600000

X(m)





X(m)

Fig. 5a. Total magnetic field of In Ouzzal

Fig. 5b. Anomaly magnetic field of In Ouzzal.

Y(m)

Y(m)

marked by ductile to brittle deformation along the eastern shear zone bordering the In Ouzzal terrane with steep fracture cleavage (NNW-SSE) and conjugate joint pattern. All these structural features are compatible with an ENE-WSW shortening in relation with the collision between the West African Craton and the Hoggar during the Pan-African orogeny.

1-Archaean granulites ; 2- Gneiss and metasediments ; 3- Gneiss with facies amphibole; 4- Indif gneiss; 5- Paleozoic curvature; 6- Panafrican granite; 7-Volcano-sediments of Tafassasset ; 8- Major faults

Fig. 4. Map Geological map of the Mole In Ouzzal extracted from the map of Hoggar (After Caby et al, 1981).

#### **4.2.3 Data processing**

In this section we have analyzed the aeromagnetic data of In Ouzzal to demonstrate the power of the 2D DCWT method to identify geological contacts. Source codes in C language are developed to calculate the 2D directional continuous wavelet transform and the spatial distribution of its maxima at different scales. The geomagnetic field data are processed with a regular grid of 750m. Figure 5a is a map of total magnetic field and figure 5b represents the anomaly magnetic field *T* . The IGRF75 model is used for calculation of *T* .

Fig. 5a. Total magnetic field of In Ouzzal

marked by ductile to brittle deformation along the eastern shear zone bordering the In Ouzzal terrane with steep fracture cleavage (NNW-SSE) and conjugate joint pattern. All these structural features are compatible with an ENE-WSW shortening in relation with the collision between the West African Craton and the Hoggar during the Pan-African

1-Archaean granulites ; 2- Gneiss and metasediments ; 3- Gneiss with facies amphibole; 4- Indif gneiss; 5- Paleozoic curvature; 6- Panafrican granite; 7-Volcano-sediments of Tafassasset ; 8- Major faults Fig. 4. Map Geological map of the Mole In Ouzzal extracted from the map of Hoggar

In this section we have analyzed the aeromagnetic data of In Ouzzal to demonstrate the power of the 2D DCWT method to identify geological contacts. Source codes in C language are developed to calculate the 2D directional continuous wavelet transform and the spatial distribution of its maxima at different scales. The geomagnetic field data are processed with a regular grid of 750m. Figure 5a is a map of total magnetic field and figure 5b represents the anomaly magnetic field *T* . The IGRF75 model is used for calculation of *T* .

orogeny.

(After Caby et al, 1981).

**4.2.3 Data processing** 

Fig. 5b. Anomaly magnetic field of In Ouzzal.

Multiscale Analysis of Geophysical Signals Using the 2D Continuous Wavelet Transform 263

400000 440000 480000 520000 560000 600000

400000 440000 480000 520000 560000 600000

X(m)

Fig. 7. Mapped maxima of the 2D DCWT for all range of scales

X(m)

Fig. 6. Modulus of the 2D DCWT of the Anomaly magnetic field reduced to the pole

a=2121m a=2611m a=3215m a=3958m a=4873m a=6000m a=7386m

a=9094m

2420000

Y(m)

2440000

2460000

2480000

2500000

Y(m)

2520000

2540000

2560000

2580000

Fig. 5c. Anomaly magnetic field of In Ouzzal after reduction to the pole (RTP)

The anomaly magnetic field is reduced to the pole; parameters of reduction to the pole are illustrated in table2. Figure 5c is magnetic anomaly field after reduction to the pole (RTP).

After RTP the maximum of the anomaly magnetic field will be found at the vertical of the physical structures. Figure 6 represents the modulus of DCWT at the scale a = 2.1 km. The analyzing wavelet is the Poisson's Kernel (See equation4). The next operation consists to calculate maxima of the modulus of the continuous wavelet transform at each scale (scales varied between 2.1 and 9.09 km). At each scale we map points of maxima in the plan. The obtained set of maxima for all ranged scales will give the geometry of geologic contacts (Figure 7).


Table 2. Parameters of the reduction to the pole of the anomaly field of In Ouzzal

400000 440000 480000 520000 560000 600000

The anomaly magnetic field is reduced to the pole; parameters of reduction to the pole are illustrated in table2. Figure 5c is magnetic anomaly field after reduction to the pole

After RTP the maximum of the anomaly magnetic field will be found at the vertical of the physical structures. Figure 6 represents the modulus of DCWT at the scale a = 2.1 km. The analyzing wavelet is the Poisson's Kernel (See equation4). The next operation consists to calculate maxima of the modulus of the continuous wavelet transform at each scale (scales varied between 2.1 and 9.09 km). At each scale we map points of maxima in the plan. The obtained set of maxima for all ranged scales will give the geometry of geologic contacts

Fig. 5c. Anomaly magnetic field of In Ouzzal after reduction to the pole (RTP)

Longitude 3° Latitude 22.5° Elevation 1000m Inclinaton 22.6° Declination -4.38°

Table 2. Parameters of the reduction to the pole of the anomaly field of In Ouzzal

X(m)

2420000

2440000

2460000

2480000

2500000

Y(m)

(RTP).

(Figure 7).

2520000

2540000

2560000

2580000

Fig. 6. Modulus of the 2D DCWT of the Anomaly magnetic field reduced to the pole

Fig. 7. Mapped maxima of the 2D DCWT for all range of scales

Multiscale Analysis of Geophysical Signals Using the 2D Continuous Wavelet Transform 265

obtained results shows robustness of the 2D DCWT. We have applied this technique at the aeromagnetic data of In Ouzzal. Obtained results are compared with the geological map and the analytic signal solutions. One can remark that the 2D directional wavelet transform is able to detect boundaries defined by geologists. Comparison with analytic signal shows that DCWT is able to identify contacts that not exist in the map of contacts defined by AS. The results of this study show that the proposed technique of edge detection based on the wavelet transform is very efficient for geological contacts analysis from maps of

420000 460000 500000 540000 580000

Fig. 9. Localization of magnetic sources obtained by analytic signal. Circles indicate position of peaks of analytic signal amplitude determinate following Blakley and Simpson Algorithm

(1986). Minimal threshold for detection of maxima is fixed to 0.5nT/km.

X(m)

geomagnetic anomalies.

2420000

 Obtained contacts by the 2D DCWT Obtained contacts by Analytic signal

2440000

2460000

2480000

2500000

Y(m)

2520000

2540000

2560000

2580000

#### **4.2.4 Results interpretation**

The obtained contacts by DCWT are compared with the geological map (Figure 8). One can remark that the proposed technique is able to identify contacts that exist in the structural geology map.

Fig. 8. Obtained contacts compared with the structural map of In Ouzzal.

To better ensure our results, we have compared these contacts with the analytic signal (AS) solutions. The analytic signal is very sensitive to noise. This is due to the derivative operator which amplifies noise intensity (Cooper, 2006). For that we use a threshold to eliminate fictitious solutions, in this case the threshold is equal to 0.5nT/km. Figure 9 shows the comparison between contacts identified by the two methods. We can note that the AS is note able to identify a contact that exist firstly in the geological map, where 2D DCWT has identified this last (blue dashed line in figure 9). It is clear that DCWT technique has detected some contacts that not exist in the model proposed by AS (See figure 9).This is due firstly at the difference between the principles of the two techniques. In the AS technique the value of the anomaly field at the level Z=Z0 is obtained by the Hilbert transform (Nabighian, 1984). However in the proposed technique the intensity of anomaly field at Z=Z0 is defined by the modulus of 2D DCWT at the scale a=Z0. The DCWT using the Poisson Kernel is a low pass filter which decrease the noise compared to the AS technique, this last amplify the noise intensity. Figure 10 is the proposed structural model based on the mapped contacts using the 2D DCWT.

#### **4.3 Conclusion**

We have proposed a technique of boundaries identification based on the 2D directional continuous wavelet transform. Firstly we have applied this idea at a synthetic model,

The obtained contacts by DCWT are compared with the geological map (Figure 8). One can remark that the proposed technique is able to identify contacts that exist in the structural

Fig. 8. Obtained contacts compared with the structural map of In Ouzzal.

To better ensure our results, we have compared these contacts with the analytic signal (AS) solutions. The analytic signal is very sensitive to noise. This is due to the derivative operator which amplifies noise intensity (Cooper, 2006). For that we use a threshold to eliminate fictitious solutions, in this case the threshold is equal to 0.5nT/km. Figure 9 shows the comparison between contacts identified by the two methods. We can note that the AS is note able to identify a contact that exist firstly in the geological map, where 2D DCWT has identified this last (blue dashed line in figure 9). It is clear that DCWT technique has detected some contacts that not exist in the model proposed by AS (See figure 9).This is due firstly at the difference between the principles of the two techniques. In the AS technique the value of the anomaly field at the level Z=Z0 is obtained by the Hilbert transform (Nabighian, 1984). However in the proposed technique the intensity of anomaly field at Z=Z0 is defined by the modulus of 2D DCWT at the scale a=Z0. The DCWT using the Poisson Kernel is a low pass filter which decrease the noise compared to the AS technique, this last amplify the noise intensity. Figure 10 is the proposed structural model based on the mapped contacts

We have proposed a technique of boundaries identification based on the 2D directional continuous wavelet transform. Firstly we have applied this idea at a synthetic model,

**4.2.4 Results interpretation** 

geology map.

using the 2D DCWT.

**4.3 Conclusion** 

obtained results shows robustness of the 2D DCWT. We have applied this technique at the aeromagnetic data of In Ouzzal. Obtained results are compared with the geological map and the analytic signal solutions. One can remark that the 2D directional wavelet transform is able to detect boundaries defined by geologists. Comparison with analytic signal shows that DCWT is able to identify contacts that not exist in the map of contacts defined by AS. The results of this study show that the proposed technique of edge detection based on the wavelet transform is very efficient for geological contacts analysis from maps of geomagnetic anomalies.

Obtained contacts by the 2D DCWT

Obtained contacts by Analytic signal

Fig. 9. Localization of magnetic sources obtained by analytic signal. Circles indicate position of peaks of analytic signal amplitude determinate following Blakley and Simpson Algorithm (1986). Minimal threshold for detection of maxima is fixed to 0.5nT/km.

Multiscale Analysis of Geophysical Signals Using the 2D Continuous Wavelet Transform 267

Nabighian (1972, 1984) developed the notion of 2-D analytic signal, or energy envelope, of magnetic anomalies. Roest, et al (1992), showed that the amplitude (absolute value) of the 3-D analytic signal at location (x,y) can be easily derived from the three orthogonal

( ) *TTT A x*

An important comment at this point is the analytic signal can be easily calculated. The x and y derivatives can be calculated directly form a total magnetic field grid using a simple 3x3 filter, and the vertical gradient is routinely calculated using FFT techniques. The analytic signal anomaly over a 2-D magnetic contact located at (*x*=0) and at depth *h* is described by

<sup>1</sup> ( ) ( )

*h x* 

The proposed idea is based on the calculation of maxima of the modulus of the 2D continuous wavelet transform of the amplitude of the analytic signal. Figure 11 shows a

> Calculation of amplitude of the 3D analytic signal (AS) of the potential field F

Calculation of the modulus of the continuous wavelet transform (CWT) of the amplitude of the analytic signal

Calculation of maxima of the modulus of the CWT of the AS for all range of scales

Mapping maxima of the modulus of the CWT

*A x*

 is the amplitude factor 2 2 

Fig. 11. Flow chart of the potential field analysis using the 3D AS and 2D CWT

2 2 2

(5)

(6)

*xyz* 

2 2 0.5

2 sin( )(1 cos( ) sin( ) )) *Md I A*

The proposed idea has been applied on the In Ouzzal area

gradients of the total magnetic field using the expression:

Where : *A x*( ) is the amplitude of the analytic signal at (*x,y*). T is the is the intensity of the observed magnetic field at (*x,y*).

the expression (after Nabighian, 1972):

*h* is the depth to the top of the contact. *M* is the strength of magnetization

*I* is the inclination of the magnetization vector *A* is the direction of the magnetization vector

**5.1 Analytic signal** 

Where

*d* is the dip of the contact

**5.2 The processing algorithm** 

detailed flow chart of this technique.

Geological contacts certain

Geological contacts supposed from multisacle analysis of aeromagnetic data.

#### **5. Mutiscale analysis of the 3D analytic signal using the 2D directional continuous wavelet transform**

Many model of potential field data interpretation are proposed, the horizontal gradient (Blakely and Simpson, 1986) and the analytic signal or full gradient (Nabighian, 1984, Roest et al, 1992) are the classical methods of the last decades. The big weakness of these techniques is sensitivity to noise (Cooper, 2006). Since the presence of noise in the potential field can give fictitious contacts. This is due to the de derivative operator that amplifies noise effect. Usually we use a threshold to eliminate these fictitious contacts. However by this way, many high frequency causative sources will be missed. For this reason we propose in this paper a technique of boundaries delimitation from geomagnetic data. This last is based on the analysis of the amplitude of 3D analytic signal defined by Roest et al (1992) by the directional continuous wavelet transform. After that maxima of the modulus of the CWT for all range of scales are mapped. The set of maxima will give geological boundaries.

Fig. 10. Proposed contacts obtained by positioning of maxima of the 2D DCWT

The proposed idea has been applied on the In Ouzzal area

#### **5.1 Analytic signal**

266 Wavelet Transforms and Their Recent Applications in Biology and Geoscience

2.2 2.4 2.6 2.8 3 3.2 3.4 3.6 3.8

Geological contacts supposed from multisacle analysis of aeromagnetic data.

**5. Mutiscale analysis of the 3D analytic signal using the 2D directional** 

Many model of potential field data interpretation are proposed, the horizontal gradient (Blakely and Simpson, 1986) and the analytic signal or full gradient (Nabighian, 1984, Roest et al, 1992) are the classical methods of the last decades. The big weakness of these techniques is sensitivity to noise (Cooper, 2006). Since the presence of noise in the potential field can give fictitious contacts. This is due to the de derivative operator that amplifies noise effect. Usually we use a threshold to eliminate these fictitious contacts. However by this way, many high frequency causative sources will be missed. For this reason we propose in this paper a technique of boundaries delimitation from geomagnetic data. This last is based on the analysis of the amplitude of 3D analytic signal defined by Roest et al (1992) by the directional continuous wavelet transform. After that maxima of the modulus of the CWT for all range of scales are mapped. The set of maxima will give

Fig. 10. Proposed contacts obtained by positioning of maxima of the 2D DCWT

21.8

Geological contacts certain

**continuous wavelet transform** 

geological boundaries.

22

22.2

22.4

22.6

22.8

23

23.2

Nabighian (1972, 1984) developed the notion of 2-D analytic signal, or energy envelope, of magnetic anomalies. Roest, et al (1992), showed that the amplitude (absolute value) of the 3-D analytic signal at location (x,y) can be easily derived from the three orthogonal gradients of the total magnetic field using the expression:

$$\left| A(\mathbf{x}) \right| = \sqrt{\left(\frac{\partial T}{\partial \mathbf{x}}\right)^2 + \left(\frac{\partial T}{\partial y}\right)^2 + \left(\frac{\partial T}{\partial z}\right)^2} \tag{5}$$

Where : *A x*( ) is the amplitude of the analytic signal at (*x,y*). T is the is the intensity of the observed magnetic field at (*x,y*).

An important comment at this point is the analytic signal can be easily calculated. The x and y derivatives can be calculated directly form a total magnetic field grid using a simple 3x3 filter, and the vertical gradient is routinely calculated using FFT techniques. The analytic signal anomaly over a 2-D magnetic contact located at (*x*=0) and at depth *h* is described by the expression (after Nabighian, 1972):

$$\left| A(\mathbf{x}) \right| = \alpha \frac{1}{\left(\hbar^2 + \mathbf{x}^2\right)^{0.5}} \tag{6}$$

Where is the amplitude factor 2 2 2 sin( )(1 cos( ) sin( ) )) *Md I A h* is the depth to the top of the contact. *M* is the strength of magnetization *d* is the dip of the contact *I* is the inclination of the magnetization vector

*A* is the direction of the magnetization vector

#### **5.2 The processing algorithm**

The proposed idea is based on the calculation of maxima of the modulus of the 2D continuous wavelet transform of the amplitude of the analytic signal. Figure 11 shows a detailed flow chart of this technique.

Fig. 11. Flow chart of the potential field analysis using the 3D AS and 2D CWT

Multiscale Analysis of Geophysical Signals Using the 2D Continuous Wavelet Transform 269

400000 440000 480000 520000 560000 600000

Fig. 14. Obtained contacts by the CWT analysis combined of AS compared with geology.

Fig. 13. Modulus of the CWT of the AS plotted at the scale a=2121m.

X(m)

Y(m)

#### **5.3 Data processing**

The proposed technique has been applied at the anomaly magnetic field without reduction to the pole( figure5b). The first operation consists to calculate the amplitude of the analytic signal of this field. This last is represented in figure 12. The modulus of the continuous wavelet transform of the amplitude of the AS is represented in figure 13; the analyzing wavelet is the Poisson's Kernel. Maxima of the modulus of the CWT are mapped for all range of scales (scales varied from 2121m to 9094m). The set of maxima will give structural boundaries.

#### **5.4 Results interpretation**

Geological contacts obtained by mapping maxima of the CWT of the amplitude the the AS are compared with geological map (See figure 14). One can remark that the proposed method is able to identify contacts defined by geology. Obtained boundaries by the proposed method are compared with contacts of analytic signal(see figure 15); for this last we have eliminated fictitious contacts dues to noise using a threshold of 0.5nT/m. We observe that by using this threshold we have eliminated a lot of contacts, for example the contact defined by the dashed line in figure 15 is identified by the CWT combined with AS, however it is not detected by analytic signal, this is due to the threshold effect, this last has been applied to reduce the noise effect on the AS.

Fig. 12. Amplitude of the analytic signal of the anomaly magnetic field of In-Ouzzal

The proposed technique has been applied at the anomaly magnetic field without reduction to the pole( figure5b). The first operation consists to calculate the amplitude of the analytic signal of this field. This last is represented in figure 12. The modulus of the continuous wavelet transform of the amplitude of the AS is represented in figure 13; the analyzing wavelet is the Poisson's Kernel. Maxima of the modulus of the CWT are mapped for all range of scales (scales varied from 2121m to 9094m). The set of maxima will give structural

Geological contacts obtained by mapping maxima of the CWT of the amplitude the the AS are compared with geological map (See figure 14). One can remark that the proposed method is able to identify contacts defined by geology. Obtained boundaries by the proposed method are compared with contacts of analytic signal(see figure 15); for this last we have eliminated fictitious contacts dues to noise using a threshold of 0.5nT/m. We observe that by using this threshold we have eliminated a lot of contacts, for example the contact defined by the dashed line in figure 15 is identified by the CWT combined with AS, however it is not detected by analytic signal, this is due to the threshold effect, this last has

400000 440000 480000 520000 560000 600000

Fig. 12. Amplitude of the analytic signal of the anomaly magnetic field of In-Ouzzal

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 1.6

X(m)

**5.3 Data processing** 

**5.4 Results interpretation** 

been applied to reduce the noise effect on the AS.

2420000

2440000

2460000

2480000

2500000

2520000

2540000

2560000

2580000

Y(m)

boundaries.

Fig. 13. Modulus of the CWT of the AS plotted at the scale a=2121m.

Fig. 14. Obtained contacts by the CWT analysis combined of AS compared with geology.

Multiscale Analysis of Geophysical Signals Using the 2D Continuous Wavelet Transform 271


290000 310000 330000

3500000

(a) (b)

292045.8 312045.8 332045.8

X(m)

(c)

3520000

3540000

3560000

3580000

Y(m)

Fig. 16. Multiscale analyse of gravity data using the 2D DCWT

(b) Modulus of the 2D DCWT plotted at the low scale a=705m (c) Obtained contacts mapped by sweeping all range of scales

3600000

3620000

3520000

3540000

3560000

3580000

3600000

3620000

290000 310000 330000

3500000

(a) Bouger Anomaly

3520000

3540000

3560000

3580000

3600000

3620000

Fig. 15. Obtained contacts by the proposed method compared with analytic signal solutions.

#### **5.5 Conclusion**

We have proposed a technique of contacts identification based on the maxima of the modulus of the 2D continuous wavelet transform of the amplitude of the 3D analytic signal. Application on the real aeromagnetic data of In Ouzzal shows that the proposed idea is able to identify contacts that are not detected by the analytic signal solutions. By implanting our method we have resolved the ambiguity of application of a threshold at the maxima of analytic signal, which can eliminate a lot of high frequency causative sources. Secondly we have resolved the ambiguity of reduction to the pole that exist in the classical methods of application of the CWT for contacts delimitation, since this last needs a data reduced to the pole, which is not very efficient in the case of high geomagnetic remanance.

#### **6. Multiscale analysis of the gravity data**

We have applied the proposed idea on the gravity data of an area located in the Algerian Sahara. Figure 16a shows the Bouguer anomaly processed with a regular grid of 500mX500m. Firstly a 2D continuous wavelet transforms has been applied, the analyzing wavelet is the Poisson's Kernel. Modulus of the 2D DCWT is presented for the low scale a=705m in figure16b.

420000 450000 480000 510000 540000 570000 600000

X(m)

Fig. 15. Obtained contacts by the proposed method compared with analytic signal solutions.

We have proposed a technique of contacts identification based on the maxima of the modulus of the 2D continuous wavelet transform of the amplitude of the 3D analytic signal. Application on the real aeromagnetic data of In Ouzzal shows that the proposed idea is able to identify contacts that are not detected by the analytic signal solutions. By implanting our method we have resolved the ambiguity of application of a threshold at the maxima of analytic signal, which can eliminate a lot of high frequency causative sources. Secondly we have resolved the ambiguity of reduction to the pole that exist in the classical methods of application of the CWT for contacts delimitation, since this last needs a data reduced to the pole, which is not very efficient in the case of high

We have applied the proposed idea on the gravity data of an area located in the Algerian Sahara. Figure 16a shows the Bouguer anomaly processed with a regular grid of 500mX500m. Firstly a 2D continuous wavelet transforms has been applied, the analyzing wavelet is the Poisson's Kernel. Modulus of the 2D DCWT is presented for the low scale a=705m in figure16b.

2420000

**5.5 Conclusion** 

geomagnetic remanance.

**6. Multiscale analysis of the gravity data** 

2440000

2460000

2480000

Y(m)

2500000

2520000

2540000

2560000

2580000

Fig. 16. Multiscale analyse of gravity data using the 2D DCWT (a) Bouger Anomaly

(b) Modulus of the 2D DCWT plotted at the low scale a=705m

(c) Obtained contacts mapped by sweeping all range of scales

Multiscale Analysis of Geophysical Signals Using the 2D Continuous Wavelet Transform 273

408000 412000 416000 420000 424000

408000 412000 416000 420000 424000

Fig. 19. Modulus of the 2D continuous wavelet transform of the noisy GPR image, plotted at

X(m)

X(m)

3468000

Fig. 18. Noisy 3D GPR image.

3468000

3470000

3472000

Y(m)

the scale a=751m

3474000

3476000

3478000

3470000

3472000

Y(m)

3474000

3476000

3478000

The maxima of the modulus of the continuous wavelet transform have been mapped for all range of scales. Figure 16c shows the map of the obtained boundaries, obtained results are compared with a seismic profile that passes by the prospected region, it exhibits a big correlation.

#### **7. Multiscale analysis of noisy 3D GPR data**

Here we present a new technique of noise effect attenuation in the 3D GPR data analysis using the 2D continuous wavelet transform. Ouadfeul and Aliouane(2010), have presented a technique of multiscale analysis of the 3D GPR data suing the CWT, it has been applied on land topographical GPR data analysis, this last play a high important role in the seismic design (See Ouadfeul and Aliouane, 2010). The proposed technique is very sensitive to noise, to demonstrate this; we have added a 05% as a white noise in the 3D GPR data analyzed in the cited paper. Figure 17 is the map of this GPR data and figure 18 presents the noisy GPR data. Modulus of the 2D CWT plotted at the scale (a=751m) is presented in figure 19. Maxima of the modulus of the continuous wavelet transform are presented in figure20. It is clear that we are not able to identify topographic orientations; this is due to the noise effect on the CWT analysis. For this reason we propose an algorithm to reduce this phenomenon, this last is based on the application of an exponential low pass filter, at the modulus of the 2D CWT for the low range of scales, the filter coefficients are presented in table 03. After application of the low pass filter, maxima of the CWT are mapped for all range of scales. Compraison of of maxima of the initial model (without noise) and the filtered model is presented in figure 21. One can remark that the filtred model is not very far from the original one, and we are now able to identify the dominant topographic orientation.

Fig. 17. Initial 3D GPR image without noise

Fig. 18. Noisy 3D GPR image.

The maxima of the modulus of the continuous wavelet transform have been mapped for all range of scales. Figure 16c shows the map of the obtained boundaries, obtained results are compared with a seismic profile that passes by the prospected region, it exhibits a big

Here we present a new technique of noise effect attenuation in the 3D GPR data analysis using the 2D continuous wavelet transform. Ouadfeul and Aliouane(2010), have presented a technique of multiscale analysis of the 3D GPR data suing the CWT, it has been applied on land topographical GPR data analysis, this last play a high important role in the seismic design (See Ouadfeul and Aliouane, 2010). The proposed technique is very sensitive to noise, to demonstrate this; we have added a 05% as a white noise in the 3D GPR data analyzed in the cited paper. Figure 17 is the map of this GPR data and figure 18 presents the noisy GPR data. Modulus of the 2D CWT plotted at the scale (a=751m) is presented in figure 19. Maxima of the modulus of the continuous wavelet transform are presented in figure20. It is clear that we are not able to identify topographic orientations; this is due to the noise effect on the CWT analysis. For this reason we propose an algorithm to reduce this phenomenon, this last is based on the application of an exponential low pass filter, at the modulus of the 2D CWT for the low range of scales, the filter coefficients are presented in table 03. After application of the low pass filter, maxima of the CWT are mapped for all range of scales. Compraison of of maxima of the initial model (without noise) and the filtered model is presented in figure 21. One can remark that the filtred model is not very far from the original one, and we are now

408000 412000 416000 420000 424000

X(m)

correlation.

**7. Multiscale analysis of noisy 3D GPR data** 

able to identify the dominant topographic orientation.

3468000

Fig. 17. Initial 3D GPR image without noise

3470000

3472000

3474000

Y(m)

3476000

3478000

Fig. 19. Modulus of the 2D continuous wavelet transform of the noisy GPR image, plotted at the scale a=751m

Multiscale Analysis of Geophysical Signals Using the 2D Continuous Wavelet Transform 275

0.0003354626279 0.006737946999 0.01831563889 0.006737946999 0.0003354626279 0.006737946999 0.1353352832 0.3678794412 0.1353352832 0.006737946999 0.01831563889 0.3678794412 1 0.3678794412 0.01831563889 0.006737946999 0.1353352832 0.3678794412 0.1353352832 0.006737946999 0.0003354626279 0.006737946999 0.01831563889 0.006737946999 0.0003354626279

In this chapter we have realize a multiscale analysis of many 2D geophysical signals by the continuous wavelet transform, the goal is to resolve many ambiguities in geophysics.

In geomagnetism the CWT proves to be a very useful tool for structural boundaries delimitation from geomagnetic data reduced to the pole. Multiscale analysis of the 3D analytic signal proves to be very powerful tool to for contact identification, this tool is less sensitive to noise and doesn't require a reduction to the pole, which is not very easy in areas

Wavelet analysis of gravity data shows that this last can be used for mapping of geological

We have proposed a new technique to reduce the noise effect in wavelet analysis of geophysical signals, application on a 3D ground penetrating radar data shows that the proposed method can be used to reduce the noise effect when analyzing geophysical data by

Arneodo,A., Decoster,N., Kestener, P., and Roux,S.G.(2003). A wavelet-based method for

Blakely , R.J., and Simpson, R.W. (1986). Approximating edges of source bodies from

Caby. R, Bertrand. J. M. L, and Black. R. (1981). Pan-African closure and ontinental collision

Caby. R, Monie. P. (2003). Neoproterozoic subduction and differential exhumation of

Cooper, G.R.J. (2006). Interpreting potential field data using continuous wavelet transforms of their horizontal derivatives, Computers & Geosciences 32 (2006) 984–992 Djemaï, S., , Haddoum, H., Ouzegane. K, And Kienast, J-R. (2009). Archaean Series Reworked

P-T Path, Bulletin du Service Géologique National, Vol. 20, n°1, pp. 3 – 29. Kumar, P., Foufoula-Georgiou, E. (1997). wavelet analysis for geophysical applications,

Mallat, S., and Hwang. (1992). Singularity detection and processing with wavelets, IEEE

multi- fractal image analysis: From theoretical concepts to experimental

in the Hoggar-Iforas segment, central Sahara. in Kroner A (ed) Precambrian Plate

western Hoggar (southwest Algeria): new sructural, petrological and

At Proterozoic In Amesmessa (West Hoggar): Cartography, Tectonic Evolution And

Table 3. An exponential low pass filter coefficients

accidents, this information can be used for seismic design.

applications. Adv. Imaging Electr. Phys. 126, 1-92.

Tectonics, Elsevier, Amst. 407-434.

Reviews of Geophysics,35,4, 385-412.

Trans. Info. Theory 38(2), 617-643, 1992.

magnetic or gravity anomalies, Geophysics : 51 ,1494–1498.

geochronological evidence, Journal of African Earth Sciences. 37.

**8. Conclusion** 

with high remanance.

the wavelet transform.

**9. References** 

Fig. 20. Mapped maxima of the modulus of the CWT of the noisy GPR model

Initial 3D GPR data

Noisy 3D GPR data after filtering of the 2D CWT coefficients

Fig. 21. Obtained maxima of the modulus of the CWT of:


Table 3. An exponential low pass filter coefficients

#### **8. Conclusion**

274 Wavelet Transforms and Their Recent Applications in Biology and Geoscience

410000 412000 414000 416000 418000 420000 422000

X(m)

Fig. 20. Mapped maxima of the modulus of the CWT of the noisy GPR model

410000 412000 414000 416000 418000 420000 422000

Noisy 3D GPR data after filtering of the 2D CWT coefficients

Fig. 21. Obtained maxima of the modulus of the CWT of:

X(m)

3470000

3472000

3474000

Y(m)

3470000

Initial 3D GPR data

3472000

3474000

3476000

3478000

Y(m)

3476000

3478000

In this chapter we have realize a multiscale analysis of many 2D geophysical signals by the continuous wavelet transform, the goal is to resolve many ambiguities in geophysics.

In geomagnetism the CWT proves to be a very useful tool for structural boundaries delimitation from geomagnetic data reduced to the pole. Multiscale analysis of the 3D analytic signal proves to be very powerful tool to for contact identification, this tool is less sensitive to noise and doesn't require a reduction to the pole, which is not very easy in areas with high remanance.

Wavelet analysis of gravity data shows that this last can be used for mapping of geological accidents, this information can be used for seismic design.

We have proposed a new technique to reduce the noise effect in wavelet analysis of geophysical signals, application on a 3D ground penetrating radar data shows that the proposed method can be used to reduce the noise effect when analyzing geophysical data by the wavelet transform.

#### **9. References**


**14** 

*Algeria* 

**1D Wavelet Transform and Geosciences** 

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

Ouadfeul and Aliouane (2011) have proposed a technique of lithofacies segmentation based

Cooper et al (2010) have published a paper focused on the fault determination using one

In seismic data processing the 1D WT has been used by many researchers to denoise the

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

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

Chamoli(2009) has analyzed the geophysical time series using the wavelet transform.

has been applied on real well-logs data of a borehole located in Algerian Sahara.

**2. Random seismic noise attenuation using the wavelet transform** 

characterization (Martelet et al, 2001, Ouadfeul et al, 2010).

on processing of well-logs data by the 1D continuous wavelet transform.

seismic data (Xiaogui Miao and Scott Cheadle, 1998, Ouadfeul, 2007).

**1. Introduction** 

dimensional wavelet analysis.

to attenuate this type of noise.

Mohamed Hamoudi2, Amar Boudella2 and Said Eladj3 *1Geosciences and Mines, Algerian Petroleum Institute, IAP* 

Sid-Ali Ouadfeul1,2, Leila Aliouane2,3,

*2Geophysics Department, FSTGAT, USTHB* 

*3Geophysics Department, LABOPHYT, FHC, UMBB* 


## **1D Wavelet Transform and Geosciences**

Sid-Ali Ouadfeul1,2, Leila Aliouane2,3, Mohamed Hamoudi2, Amar Boudella2 and Said Eladj3 *1Geosciences and Mines, Algerian Petroleum Institute, IAP 2Geophysics Department, FSTGAT, USTHB 3Geophysics Department, LABOPHYT, FHC, UMBB Algeria* 

#### **1. Introduction**

276 Wavelet Transforms and Their Recent Applications in Biology and Geoscience

Martelet, G., Sailhac, P., Moreau, F., Diament, M.(2001). Characterization of geological

Moreau, F., D. Gibert, M. Holschneider, and G. Saracco. (1999). Identification of sources of

Moreau , F., Gibert, D., Holschneider M., Saracco G. (1997). Wavelet analysis of potential

Murenzi, R., 1989, Transformée en ondelettes multidimensionnelle et application à l'analyse

Nabighian, M. N. 1(972). The analytic signal of two-dimensional magnetic bodies with

Nabighian, M.N. (1984). Toward a three-dimensional automatic interpretation of potential field data via generalized Hilbert transforms: Fundamental relations: Geophysics, 49 ,957–966. Ouzegane K, Boumaza. S. (1996). An example of ultra-high temperature metamorphism :

Ouadfeul , S. (2006). Automatic lithofacies segmentation using the wavelet transform

Ouadfeul , S. ( 2007). Very fines layers delimitation using the wavelet transform modulus maxima lines WTMM combined with the DWT , SEG SRW , Expanded abstract. Ouadfeul, S.(2008). Reservoir characterization using the continuous wavelet transform

Ouadfeul, S., and Aliouane, L. (2010). Multiscale of 3D GPR data using the continuous

Ouadfeul, S, Aliouane, L., Eladj, S. (2010). Multiscale analysis of geomagnetic data using the

Reid, A.B., Allsop, J.M., Granser, H., Millett, A.J., Somerton, I.W.(1990). Magnetic

Roest, W.R, Verhoef, J, and Pilkington, M. (1992). Magnetic interpretation using the 3-D

Sailhac. P, Galdeano. A, Gibert. D, Moreau. F, and Delor. C. (2000). Identification of sources

Sailhac, P., Galdeano, A., Gibert, D., Moreau, F., Delor, C. (2000). Identification of sources of

Valee, M.A., Keating, P., Smith, R.S., St-Hilaire, C. (2004). Estimating depth and model type using the continuous wavelet transform of magnetic data. Geophysics 69, 191–199.

approximations, Journal of geophysical research, vol. 108, noB5.

Abstracts 29, 1222 (2010); doi:10.1190/1.3513065.

signal analytic, Geophysics, 57:116-125.

Research 105, 19455–19475.

interpretation: Geophysics, 37, 507–517, doi: 10.1190/1.1440276.

the Himalayas. Geophysics 66, 1116–1129.

fields. Inverse Problem, 13, 165-178.

d'images, Thèse Louvain-La-Neuve.

Res., 104, 5003–5013.

Geology, 14: 693-708.

Expanded abstract .

Expanded abstract.

boundaries using 1-D wavelet transform on gravity data: theory and application to

potential fields with the continuous wavelet transform: Basic theory, J. Geophys.

polygonal cross-section: Its properties and use of automated anomaly

orthpyroxene-sillimnite, -garnet, sapphirine-quartz and spinel-quartz parageneses in Al-Mg granulites froum In Hihaou, In Ouzzal, Hoggar, Journal of Metamorphic

modulus maxima lines (WTMM) combined with the detrended fluctuation analysis (DFA), 17the International geophysical congress and exhibition of turkey,

combined with the Self Organizing map (SOM) neural network: SPE ECMOR XI,

wavelet transform, presented in GPR2010,IEEE, doi: 10.1109/ ICGPR.2010.5550177.

continuous wavelet transform. Application to Hoggar (Algeria), SEG Expanded

interpretation in three dimensions using Euler deconvolution: Geophysics, 55, 80–91.

of potential fields with the continuous wavelet transform: Complex wavelets and applications to magnetic profiles in French Guiana, J, Geophy, Res., 105: 19455-75. Sailhac , P., and Gibert , D. (2003). Identification of sources of potential fields with the

continuous wavelet transform: Two-dimensional wavelets and multipolar

potential fields with the continuous wavelet transform: complex wavelets and application to aeromagnetic profiles in French Guiana. Journal of Geophysical 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 and Scott 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.

#### **2. 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.

1D Wavelet Transform and Geosciences 279

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 *z0*. More precisely, we have the following power-law relation (Hermann, 1997**;** Audit *et al.,*

( )

understood as a global indicator of the local differentiability of a function s.

correlated (*h >* 0.5), uncorrelated (*h* = 0.5) or anticorrelated (*h <* 0.5).

*a*

*j*

Hence its contrary transform is 2

2 2

**2.2 The discrete wavelet transform (DWT)** 

The wavelet transform of a function 2 2

is the dilation of

2( ) *<sup>j</sup> a jZ* , then the wavelet is 2 2

( ) ( )\* (*t*) *a t f t*

2 2

There <sup>1</sup> () ( ) *a a <sup>t</sup> <sup>t</sup> a a*

 

ˆ(2 ) (2 ) 1 *j j*

 *wx w* 

() 0 *t dt*

Where *h* is the Hölder exponent (or singularity strength). The Hölder exponent can be

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

*L2(R)* denotes the Hibert space of measurable, square-integrable functions. The function

() ( ) *t LR* is said to be a wavelet if and only when the following condition is satisfied.

() ( ) *t LR* is defined by :

In order to be used expediently in practice a, is scattered as discrete binary system ,i.e. Let

*t*

<sup>1</sup> ( ) ( )\* ( ) ( )\* ( ) 2 2 *<sup>j</sup>*

<sup>1</sup> () ( ) 2 2 *j j j j*

 

*<sup>t</sup> W ft ft t ft* 

( )*t* by the scale factor *s*.

*t* ,its wavelet transform is :

*j j*

(6)

( ) ( )\* ( ) *<sup>j</sup> f t W f t xt* . Where *x(t)* satisfies

  (5)

<sup>0</sup> (, ) *h zo C az a <sup>s</sup>* , when *<sup>a</sup>* <sup>0</sup> (4)

2002):

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.

#### **2.1 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:

$$C\_{\mathcal{S}}(a,b) = \frac{1}{a} \int\_{-\infty}^{+\infty} s(z) \nu^\* \left(\frac{z-b}{a}\right) dz$$

Each family test function is derived from a single function ( ) *z* defined to as the analyzing wavelet according to (Torresiani, 1995):

$$
\psi\_{a,b}(z) = \psi(\frac{z-b}{a})\,'\,\tag{2}
$$

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 :

$$\int\_{-\infty}^{+\infty} z^n \varphi(z) \, dz = 0 \quad \text{for} \quad 0 \le n \le p - 1 \tag{3}$$

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

Siyuan Cao and Xiangpeng Chen (2005) have used the second-generation wavelet transform

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

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

> *<sup>z</sup> <sup>b</sup> <sup>s</sup> <sup>z</sup> <sup>a</sup> <sup>a</sup> <sup>b</sup> Cs* ( ) ( ) <sup>1</sup> ( , )

> > , () ( ) *a b z b*

Where *a R* is a scale parameter, *b R* is the translation and \* is the complex conjugate

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

() 0 0 1 *nz z dz for n p*

 

*a*

*z*

 

*dz <sup>a</sup>*

( ) *z* is generally chosen to be well localized in space (or time)

(3)

, (1)

, (2)

( ) *z* defined to as the analyzing

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

polynomials, up to the degree *n*−1, i.e., to have *n* vanishing moments :

wavelet transform.

for random noise attenuation.

results interpretation and conclusion.

frequency simultaneously is not possible.

wavelet according to (Torresiani, 1995):

of . The analyzing function

**2.1 The continuous wavelet transform (CWT)** 

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 *z0*. More precisely, we have the following power-law relation (Hermann, 1997**;** Audit *et al.,* 2002):

$$\left|\mathbf{C}\_{s}(a, z\_{0})\right| \approx a^{l\text{(zo)}} \quad , \text{ when} \quad a \to 0^{+} \tag{4}$$

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

#### **2.2 The discrete wavelet transform (DWT)**

*L2(R)* denotes the Hibert space of measurable, square-integrable functions. The function 2 2 () ( ) *t LR* is said to be a wavelet if and only when the following condition is satisfied.

$$\int\_{-\infty}^{\infty} \varphi(t)dt = 0$$

The wavelet transform of a function 2 2 () ( ) *t LR* is defined by :

$$
\psi\_{\mathcal{Q}}(t) = f(t)^\* \psi\_{\mathcal{Q}}(t) \tag{5}
$$

There <sup>1</sup> () ( ) *a a <sup>t</sup> <sup>t</sup> a a* 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 2( ) *<sup>j</sup> a jZ* , then the wavelet is 2 2 <sup>1</sup> () ( ) 2 2 *j j j j t t* ,its wavelet transform is :

$$\mathcal{W}\_2^{\ j} f(t) = f(t) \,^\* \boldsymbol{\nu}\_{\underline{\boldsymbol{\nu}}^{j}}(t) = f(t) \,^\* \frac{1}{2^{\dot{\boldsymbol{\nu}}}} \boldsymbol{\nu}(\frac{t}{2^{\dot{\boldsymbol{\nu}}}}) \tag{6}$$

Hence its contrary transform is 2 ( ) ( )\* ( ) *<sup>j</sup> f t W f t xt* . Where *x(t)* satisfies ˆ(2 ) (2 ) 1 *j j wx w* 

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)

$$\begin{aligned} S\_{2^j}f &= S\_{2^{j-1}}f \ast H\_{j-1} \\ \mathcal{W}\_{2^j}f &= S\_{2^{j-1}}f \ast G\_{j-1} \end{aligned} \tag{7}$$

1D Wavelet Transform and Geosciences 281


Discrete wavelet transform of the seismic seismogram T(t)

Densoing of T(t) using the threshold method Td(t) is the denoised seismic seismogram

Fig. 2. Flow chart of the proposed technique of the random noise attenuation form seismic

Plotting of the CWT(t,a) at the scale a=2\*∆T ∆T : is the sampling interval

Calculation of the continuous wavelet transform CWT(t,a) of Td(t) The analyzing wavelet is the modified Mexican Hat.

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

Fig. 1. Graph of the modified Mexican wavelet.

data.

**2.5 Application on synthetic data** 

and scale are presented in figure6.


There 2 *<sup>j</sup> W f* is the wavelet transform coefficients of *f(t).* It approximates *f(t)* on the scale 2j .

H*j* and G*j* are the discrete filters gained by inserting (2*<sup>j</sup>* - 1) zeros into every two samples of H, G. And the relation between G and H is:

$$\mathcal{g}\_k = (-1)^{k-1} \overline{h}\_{1-k} \tag{8}$$

#### **2.3 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.

#### **2.4 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.

$$\nu(t) = \begin{cases} 1 & \text{1...if (|t| \le 1)} \\ & -1 \ne 2 \dots \text{if (1 < |t| \le 3)} \\ & 0 \dots \dots \text{if (|t| > 3)} \end{cases} \tag{9}$$

Being dispersed in time domain farther, a discrete wavelet transform can be obtained. It

*Sf S f H Wf S fG* 

*j j j j*

1 1 1 2 2 1 2 2

*<sup>j</sup> W f* is the wavelet transform coefficients of *f(t).* It approximates *f(t)* on the scale 2j .

1 <sup>1</sup> ( 1)*<sup>k</sup>*

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

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

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

The analyzing wavelet is the modified Mexican Hat, it is defined by equation 9 (See figure1).

1...... ( 1) ( ) 1 / 2..... (1 3) 0....... ( 3)

*if t*

*if t t if t*

The flow chart of the proposed processing algorithm is detailed in figure2.

\* \*

*j j*

 

(7)

*k k g h* (8)


(9)

exists an effective and fast algorithm, it is based on equation (7)

H*j* and G*j* are the discrete filters gained by inserting (2*<sup>j</sup>*

H, G. And the relation between G and H is:

There 2

**2.3 Signal denoising** 

reconstruct the data set.

scale a=2.\*∆T. Where:

∆T is the sampling interval.

**2.4 The denoising algorithm** 

apply a threshold to denoise the seismic trace.


Fig. 1. Graph of the modified Mexican wavelet.

Fig. 2. Flow chart of the proposed technique of the random noise attenuation form seismic data.

#### **2.5 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.

1D Wavelet Transform and Geosciences 283

Fig. 4. Discrete wavelet decomposition of the synthetic seismic seismogram

Fig. 5. Denoised synthetic seismogram using the thresholding method of the DWT

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.


Table 1. Acoustic parameters of the synthetic model.

Fig. 3. Noisy synthetic seismic seismogram.

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

First Layer

Second Layer

Third Layer

Fourth Layer

]. The next operation consists to apply the

 , + 

Thickness 800m Velocity of the P wave 2500m/s

Thickness 700m Velocity of the P wave 3625m/s

Thickness 800m Velocity of the P wave 4000m/s

Thickness +∞ Velocity of the P wave 4500m/s

Table 1. Acoustic parameters of the synthetic model.

Fig. 3. Noisy synthetic seismic seismogram.

this noise contains all the angles [-

sketched method at real data.

Fig. 4. Discrete wavelet decomposition of the synthetic seismic seismogram

Fig. 5. Denoised synthetic seismogram using the thresholding method of the DWT

1D Wavelet Transform and Geosciences 285

0 100 200 300 400 500 Frequency (H z)

0 100 200 300 400 500

Frequency (H z)

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

0,00 0,25 0,50 0,75 1,00 1,25 1,50 1,75 2,00

t(s)

Fig. 8. Seismic seismogram of a raw vertical seismic profile realized in Algeria.

Fig. 7b. Spectral analysis of the residual noise using the Fourier's transform

a=0.004s. This last is presented in figure12 versus the time.

0

**2.6 Application on real data** 





0,00

0,05

0,10

0,15

2

4

Amplitude

6


Angle(deg)

8

Fig. 6. CWT coefficients of the denoised synthetic seismogram suing the DWT. The analyzing wavelet is the modified Mexican.

Fig. 7a. CWT coefficients at the scale a= 0.004s

0 0.5 1 1.5 2 2.5

Fig. 6. CWT coefficients of the denoised synthetic seismogram suing the DWT. The

t(s)




0

50

100

150

200


analyzing wavelet is the modified Mexican.

Fig. 7a. CWT coefficients at the scale a= 0.004s



a(s)



0

Log(a(m))

Fig. 7b. Spectral analysis of the residual noise using the Fourier's transform

#### **2.6 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.

Fig. 8. Seismic seismogram of a raw vertical seismic profile realized in Algeria.

1D Wavelet Transform and Geosciences 287

012

0,0 0,2 0,4 0,6 0,8 1,0 1,2 1,4 1,6 1,8 2,0

t(s)

Fig. 12. Modulus of the CWT at the scale a=0.004s, the analyzing wavelet is the modified

t(s)

Fig. 11. Modulus of the CWT of the denoised VSP seismogram using the DWT.





0

5

10

15

20

25

30

35


Mexican Hat.











a(s)

Log(a(m))

0

Fig. 9. Discrete wavelet decomposition using the Haar of level 5 of the raw VSP seismogram.

Fig. 10. Denoised VSP seismogram using the DWT with the Haar wavelet of level 5 of the VSP seismogram.

Fig. 9. Discrete wavelet decomposition using the Haar of level 5 of the raw VSP

0,00 0,25 0,50 0,75 1,00 1,25 1,50 1,75 2,00

t(s)

Fig. 10. Denoised VSP seismogram using the DWT with the Haar wavelet of level 5 of the

seismogram.


VSP seismogram.




0,00

0,05

0,10

0,15

Fig. 11. Modulus of the CWT of the denoised VSP seismogram using the DWT.

Fig. 12. Modulus of the CWT at the scale a=0.004s, the analyzing wavelet is the modified Mexican Hat.

1D Wavelet Transform and Geosciences 289

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

Reading of the X and Y components from the InterMagnet network

Calculation of the CWT

Calculation of maxima of the CWT at the low scale a0

Fig. 14. Flow chart of the proposed technique of geomagnetic disturbances analysis

We have applied the proposed technique at the horizontal component of the magnetic field of Wingst observatory. Figure 15 presents the fluctuation of this component versus the time.

Estimation of the Holer exponents at maxima of the modulus of the CWT

0 10000 20000 30000 40000

Fig. 15. Horizontal component of the magnetic field of Mai 2002, recorded by the Wingst

t (m in u te )

**3.1 Fractal analysis of geomagnetic disturbances using the CWT** 

14.

H(nT)

observatory

Fig. 13. Spectral analysis of the residual noise using the Fourier's transform.

#### **2.7 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.

#### **3. 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 *Inte*rnational *R*eal-time *Mag*netic Observatory *Net*work (InterMagnet) observatories are used.

Frequency (Hz)

0 100 200

0 50 100 150 200 250

Frequency (Hz)

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

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

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 *Inte*rnational *R*eal-time *Mag*netic Observatory *Net*work (InterMagnet) observatories are

Fig. 13. Spectral analysis of the residual noise using the Fourier's transform.

**3. Solar geomagnetic disturbances analysis using the CWT** 

0,0

sweep the full interval [-π,+π].

seismic data processing software's.

used.

**2.7 Results Interpretation and conclusion** 

2,5

Amplitude

5,0


Angle(deg)

#### **3.1 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.

Fig. 14. Flow chart of the proposed technique of geomagnetic disturbances analysis

We have applied the proposed technique at the horizontal component of the magnetic field of Wingst observatory. Figure 15 presents the fluctuation of this component versus the time.

Fig. 15. Horizontal component of the magnetic field of Mai 2002, recorded by the Wingst observatory

1D Wavelet Transform and Geosciences 291

DST Index

Average Local Holder exponent

2 4 6 8 10 12 14 16 18 20 22 24 26 28 30

t(Jour)

t(day)

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

> *<sup>q</sup> D q q*

( ) *q* is the spectrum of exponent estimated using the function of partition (See Arneodo et

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 D0, D1 and D2 for every 60 minutes (one hour of the month). Obtained results are presented in

(9)

Fig. 18. Average Local Hölder exponent compared with the DST index

**3.2 Generalized fractal dimension and geomagnetic disturbances** 

( ) ( ) ( 1)

**3.2.1 Application on the Hermanus observatory data** 

al, 1955, Ouadfeul et al, 2011).

0,0

figure 19.

0,2

0,4

0,6

0,8

1,0

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 ( *h* 0) .

Fig. 16. Modulus of the 1D CWT of the horizontal component of the geomagnetic field of Hermanus observatory

Fig. 17. Local Hölder exponent estimated by the CWT

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

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





Fig. 16. Modulus of the 1D CWT of the horizontal component of the geomagnetic field of

0 500 1000 1500 2000 2500 3000 3500 4000

0 10000 20000 30000 40000

t(m in u te )

for solar geomagnetic disturbances schedule.

Hölder exponent has a very low value ( *h* 0) .

Hermanus observatory

0,0

Fig. 17. Local Hölder exponent estimated by the CWT

0,5

1,0

Holder exponent

1,5

2,0

2,5

Fig. 18. Average Local Hölder exponent compared with the DST index

#### **3.2 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).

$$D(q) = \frac{\tau(q)}{(q-1)}\tag{9}$$

 ( ) *q* is the spectrum of exponent estimated using the function of partition (See Arneodo et al, 1955, Ouadfeul et al, 2011).

#### **3.2.1 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 D0, D1 and D2 for every 60 minutes (one hour of the month). Obtained results are presented in figure 19.

1D Wavelet Transform and Geosciences 293

**4. Heterogeneities analysis using the 1D wavelet transform modulus maxima** 

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 and Aliouane, 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**;**

<sup>0</sup> (, ) *h zo C az a <sup>s</sup>*

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

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

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

are : the Clay , The sandstone, the Sandy clay and the clayey sandstone.

Sandstone deposits constitute very important hydrocarbon reservoirs.

( )

**lines** 

Audit et al., 2002) :

C a,z s 0 : Is the modulus of the CWT at the depth z0.

h(z0) : is the Hölder exponent at the depth z0

**4.1.1 Geological setting of the Berkine Basin** 

**4.1 Application on real data** 

Where :

a : is the scale

topography.

Fig. 19. Fractal dimensions calculated for each hour of the total field recorded in the period of May 2002. (a): q = 0, (b): q = 1, (c): q = 2

Analysis of obtained results shows that the fractal dimension D0 is not sensitive to magnetic disturbances. However D1 and D2 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).


Table 2. Magnetic storms recorded in a Month May 2002

#### **4. 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 and Aliouane, 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) :

$$\left| \mathbb{C}\_s(a, z\_0) \right| \approx a^{h(zo)}$$

Where :

292 Wavelet Transforms and Their Recent Applications in Biology and Geoscience

2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32

(c)

t(day)

D2

Fig. 19. Fractal dimensions calculated for each hour of the total field recorded in the period


Analysis of obtained results shows that the fractal dimension D0 is not sensitive to magnetic disturbances. However D1 and D2 are very sensitive to the solar geomagnetic activity.

One can remark that the major magnetic disturbances are characterized by spikes ( See table

(in 11 , A=37,Kmax=6)

2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32

t(day)

(in 14 , A=29,Kmax=5)

Date Starting hour Importance 11 10.13 A strong storm

14 xx.xx A storm

Table 2. Magnetic storms recorded in a Month May 2002

23 10.48 Violent storm in 23

of May 2002. (a): q = 0, (b): q = 1, (c): q = 2

0,4

2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32

t(day)

0,5

0,6

0,7

D0

0,8

0,9

(a)

2 and Figure 19).

(b)

D1

C a,z s 0 : Is the modulus of the CWT at the depth z0.

a : is the scale h(z0) : is the Hölder exponent at the depth z0

#### **4.1 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.

#### **4.1.1 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.

1D Wavelet Transform and Geosciences 295

2850 2900 2950 3000 3050

0 -10

2850 2900 2950 3000 3050 3100 Z(m) Depth(m)

2850 2900 2950 3000 3050 3100

Depth(m)

40 90

Z(m)

<sup>2850</sup> <sup>2900</sup> <sup>2950</sup> <sup>3000</sup> <sup>3050</sup> <sup>3100</sup>

Fig. 22. Modulus of the continuous wavelet transform of the gamma ray log

Z (m )

0

Fig. 23. Skeleton of the modulus of the CWT of the Gr log

1

2

Log(a(m))

3

4

5

Fig. 21. Natural Gamma ray log of Sif-Fatima2 borehole

Log(a(m))

Gr(API)

log(a(m))

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.

#### **4.1.2 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.

Fig. 20. Geographic location of the Berkine Basin

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

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.

Fig. 20. Geographic location of the Berkine Basin

petrophysical characteristics.

**4.1.2 Data processing** 

Fig. 21. Natural Gamma ray log of Sif-Fatima2 borehole

Fig. 22. Modulus of the continuous wavelet transform of the gamma ray log

Fig. 23. Skeleton of the modulus of the CWT of the Gr log

1D Wavelet Transform and Geosciences 297

GRmin < GR <GRThreshold

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

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

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

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

Arneodo , A., Grasseau, G. and Holschneider, M. (1988). Wavelet transform of multifractals,

Arneodo, A. et Bacry E. (1995). Ondelettes, multifractal et turbelance de l'ADN aux

Audit, B., Bacry, E.,Muzy,J-F. and Arneodo , A. (2002). Wavelet-Based Estimators of Scaling

Chamoli, A. (2009). Wavelet Analysis of Geophysical Time Series*, e-Journal Earth Science* 

Charles K. Chui, 1992, An Introduction to Wavelets, (1992), Academic Press, San Diego,

Grossman, A. and Morlet, J. (1985). Decomposition of functions into wavelets of constant

Herrmann, F.J.( 1997). A scaling medium representation, a discussion on well-logs, fractals

shape and related transforms , *in* : Streit , L., ed., mathematics and physics ,lectures

and waves, Phd thesis Delft University of Technology, Delft, The Netherlands,


Hölder exponent can be used as an attribute to seek the fines lithofacies.

geomagnetic disturbances. Obtained results show the robustness of the CWT.

croissances cristalines Diderot editeur arts et sciences, Paris.

on recents results, *World Scientific Publishing , Singapore.*

Are considered as Sandy Clay.

table3 .

**5. Conclusion** 

intercalations.

**6. References** 

*Phys. Rev. Lett*., 61, pp. 2281-2284.

*India,* 2(IV), 258-275

ISBN 0585470901.

pp.315.

Behavior , *IEEE* ,vol.48, pp. 2938-2954.

Are considered as a Clayey Sandstone.

Fig. 24. Estimated Hölder exponents compared with the classical interpretation based on the Gr log.


Table 3. Lithofacies intervals derived from the GR signal for Sif-Fatima 2 well

#### **4.1.3 Results Interpretation**

A preliminary raw lithofacies classification based on the natural gamma radioactivity welllog data was made. First, recall that the maximum value GRmax of the data is considered as a full clay concentration while the minimum value GRmin represents the full sandstone concentration. The mean value (GRmax+GRmin)/2 will then represent the threshold that will be used as a decision factor within the interval studied:


GRThreshold < GR <GRmax

Are considered as Sandy Clay.

296 Wavelet Transforms and Their Recent Applications in Biology and Geoscience

Fig. 24. Estimated Hölder exponents compared with the classical interpretation based on the

2838.50 -2887.20 Clayey Sandstone with increase of the percentage of clay with

2969.00 -2997.50 Clayey Sandstone becoming clean at the bottom with the

A preliminary raw lithofacies classification based on the natural gamma radioactivity welllog data was made. First, recall that the maximum value GRmax of the data is considered as a full clay concentration while the minimum value GRmin represents the full sandstone concentration. The mean value (GRmax+GRmin)/2 will then represent the threshold that will

GRThreshold < GR <GRmax

Depth interval (m) Petrographical description

3062.50 -3082 Sandy clay

**4.1.3 Results Interpretation** 

depth

2887.20- 2913.00 Clay, sometimes slightly sandy

be used as a decision factor within the interval studied:

2913.00-2950.50 Metric alternating of Clayey Sandstone and clay

intercalation of sandy clay 2997.50 -3017 Thick layer of clay sometimes slightly sandy 3017 -3062.50 Metric alternating of sandstone and clay

Table 3. Lithofacies intervals derived from the GR signal for Sif-Fatima 2 well


2950.50-2969.00 Clay with the presence of a layer of Clayey Sandstone

Gr log.


$$\mathbf{GR}\_{\text{min}} \rightharpoonup \mathbf{GR} \lessapprox \mathbf{GR}\_{\text{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.

#### **5. 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.

#### **6. References**

	- *http://www.intermagnet.org/Welcom\_f.php*

Martelet, G., Sailhac, P., Moreau, F., Diament, M. (2001). Characterization of geological

Morris Cooper, S., Tianyou, L., Ndoh Mbue,I. (2010). Fault determination using one dimensional wavelet analysis, *Journal of American Science*, 6(7):177-182. Ouadfeul, S. (2007). Very fines layers delimitation using the wavelet transform modulus maxima lines(WTMM) combined with the DWT, *SEG SRW*, Antalya, Turkey. Ouadfeul, S. (2011). Analyse multifractal des signaux géophysiques, Editions universitaires

Ouadfeul, S., Aliouabe, L. (2011). Automatic lithofacies segmentation using the wavelet

Ouadfeul, S., Eladja, S., Aliouane, L. (2010). Structural boundaries form geomagnetic data

Prasad, N. B. R. (2006). Attenuation of Ground Roll Using Wavelet Transform, 6th *International Conference & Exposition on Petroleum Geophysics,* Kolkata. Siyuan Cao and Xiangpeng Chen. (2005). The second-generation wavelet transform and its

Torréasiani, B. (1995). Analyse continue par onde lettes, Inter Editions / CNRS Edition. Xiaogui Miao and Scott Cheadle. (1998). Noise attenuation with Wavelet transforms, *SEG* 

Zeroug, S. , Bounoua , N., Lounissi , R. (2007). Algeria Well Evaluation Conference

*http://www.slb.com/resources/publications/roc/algeria07.aspx* 

analysis *, Arabian Journal of Geosciences*, DOI: 10.1007/s12517-011-0383-7. Ouadfeul, S., Aliouane, L. (2011). Multifractal Analysis Revisited by the Continuous Wavelet

*Journal of Applied Physics and Mathematics*, Vol.1, No.1.

transform modulus maxima lines combined with the detrended fluctuation

Transform Applied in Lithofacies Segmentation from Well-Logs Data, *International* 

using the continuous wavelet transform, Arabian Journal of geosciences, DOI:

application in denoising of seismic data, *Applied geophysics*, Vol.2, No 2, 70-74, DOI:

boundaries using 1-D wavelet transform on gravity data: theory and application to

International Real-time Magnetic Observatory Network, *http://www.intermagnet.org/Welcom\_f.php*

the Himalayas. *Geophysics,* 66, 1116–1129.

europeennes, ISBN 978-613-1-58257-8.

10.1007/s12517-010-0273-

10.1007/s11770-005-0034-4.

*expanded abstract*, 17, 1072.

### *Edited by Dumitru Baleanu*

This book reports on recent applications in biology and geoscience. Among them we mention the application of wavelet transforms in the treatment of EEG signals, the dimensionality reduction of the gait recognition framework, the biometric identification and verification. The book also contains applications of the wavelet transforms in the analysis of data collected from sport and breast cancer. The denoting procedure is analyzed within wavelet transform and applied on data coming from real world applications. The book ends with two important applications of the wavelet transforms in geoscience.

Wavelet Transforms and Their Recent Applications in Biology and Geoscience

Wavelet Transforms and

Their Recent Applications in

Biology and Geoscience

*Edited by Dumitru Baleanu*

Photo by lixai / iStock