Fractal and Chaos in Exploration Geophysics

1. Introduction

The fractal analysis has been widely used in exploration geophysics. In gravity and magnetism it is used for causative sources characterization [1,2, 3, 4]. In seismology, the fractal analysis is used for earthquake characterization [5, 6; 7]. In petrophysics the fractal analysis is ued for lithofacies classification and reservoir characterization [8, 9]. We cite for example the paper of Lozada-Zumaeta et al [10], they have analyzed the distribution of petrophysical properties for sandy-clayey reservoirs by fractal interpolation. San José Martínez et al [11] have published a paper for the detection of representative elementary area for multifractal analysis of soil porosity using entropy dimension.

The current chapter is composed of three application of the fractal analysis in geophysics. The first one consists to use the fractal analysis for facies identification from seismic data, the proposed idea is based on the estimation of the generalized fractal dimensions. After that an application to the pilot KTB borehole is realized The second part of this chapter is to apply the same technique but on the well-logs data for another objectives. This last, consist to identify heterogeneities. Application on the Pilot KTB borehole shows a robustness of this last. The third application of the fractal formalism for heterogeneities analysis from synthetic amplitude versus offset (AVO) data.

2. Facies recognition from seismic data using the fractal analysis

The fractal analysis has been widely used for seismic data processing, we cite for example the paper of Rivastava and Sen [11]. They have developed a new fractal-based stochastic inversion of poststack seismic data using very fast simulated annealing. A Simultaneous stochastic inversion of prestack seismic data using hybrid evolutionary algorithm based on the fractal process has been developed recently by Saraswat and Sen [12]

Facies recognition from seismic data using the fractal analysis has becoming a very interesting subject of research; in fact many papers have been published in this topic.

We cite for example the paper of Lopez and Aldana [13], these last have used the wavelet based fractal analysis for facies recognition and waveform classifier from

Oritupano-A field.

Barnes [14] has used the fractal analysis of Fault Attributes Derived from Seismic Discontinuity Data. Nath and Dewangan [15], have established a technique to detect reflections from seismic attributes. Gholamy et al [16], have proposed a technique based on frcatal methods for automatic detection of interfering seismic wavelets.

Here, we present a wavelet-based fractal analysis technique for reflection recognition from seismic data; we used the so-called wavelet transform modulus maxima lines, to calculate the generalized fractal dimensions. We start by describing the principles of the continuous wavelet transform and the WTMM. The next section consists to explain the principle of the seismic data processing algorithm. The proposed idea is applied to the synthetic seismic seismogram of the Germann KTB pilot borehole data. We finalize the paper by a conclusions and a discussion of the obtained results.

2.1. The Wavelet Transform Modulus Maxima lines (WTMM) method

The wavelet transform modulus maxima lines (WTMM) is a multifractal formalism proposed by Arneodo et al [17], the WTMM is based on the continuous wavelet transform. Lets us define an analyzing wavelet$$\psi (z)$$, the continuous wavelet transform of a signal S(z) is given by [18]:

$$C_s(a,b)=\frac{1}{\sqrt{a}}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}s(z){\psi }^{\ast }}(\frac{z-b}{a})dz$$E1

Where a: is a scale parameter.

$$\psi (\frac{z-b}{a})$$Is the dilated version of the analyzing wavelet.$${\psi }^{*}(z)$$ is the conjugate of $$\psi (z)$$

The analyzing wavelet must check the admissibility condition:

$${\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}}\psi (z)dz=0$$ in frequency domain this condition is equivalent to :

$${\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}}\frac{\psi (w)}{\omega }dw=0$$E2

For some data processing requirements $$\psi (z)$$must have N vanishing moments:

$${\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{z}^n}\psi (z)dz=0$$ for all n, 0≤n≤N-1

The first step on the WTMM consists to calculate the modulus of the continuous wavelet transform, after that the maxima of this modulus are calculated. We call maxima of the modulus of the CWT at point $b_0$ if for all $b\to b_0\mathrm{ }\left|C_s(a,b_0)\right|>\left|C_s(a,b)\right|$.

The next step is to calculate the function of partition $Z\left(q,a\right)$. If one call $L\left(b\right)$ the set of maxima of the modulus of the CWT, the function of partition is defined by :

$Z\left(q,a\right)={{\displaystyle {\sum }_{L\left(b\right)}\left|C_s\left(a,b\right)\right|}}^q$E3

Where q is the moment order.

The spectrum of exponents $\tau \left(q\right)$ is related to the function of partition for law scales by :

$\mathrm{ }Z\left(q,a\right)=a^{\tau \left(q\right)}$ if

$$a\to 0$$E4

Note that the spectrum of exponents is estimated by a simple linear regression of $\mathrm{l}\mathrm{o}\mathrm{g}\left(Z\left(q,a\right)\right)$ versus $\mathrm{l}\mathrm{o}\mathrm{g}\left(a\right)$

The generalized fractal dimension is given by [17]:

$${\mathrm{ }D}_q=\frac{\tau \left(q\right)}{(q-1)}$$E5

One can distinguish three important values of $D_q$:

If q=0: the definition becomes the capacity dimension.

If q=1: the definition becomes the information dimension.

If q=2: the definition becomes the correlation dimension.

2.2. The processing algorithm

The proposed idea consists to apply the so-called the wavelet transform modulus maxima lines (WTMM) method with a moving window of 128 samples, the window center will be moved at each 64 samples of the seismic trace. Note that 128 is the less required number of samples for the convergence of the function of partition. The analyzing wavelet is the complex Morlet [18]. After that, the three generalized fractal dimensions that consist to q=1, 2 and 3 are calculated. The goal is to show that these fractal dimensions can be used for thin bed and facies identification and eliminate the noise effect that can give false stratigraphic formations. The detailed flow chart of the proposed idea is presented in figure 1.

2.3. Application to KTB synthetic seismic data

To check the efficiency of the proposed technique, we have analyzed the Kontinentales Tiefbohrprogramm de Bundesreplik Deutschland (KTB) borehole synthetic seismic seismogram calculated from the well-logs data.

Figure 1.

Flow chart of the multifractal analysis of seismic data using the generalized fractal dimensions

2.3.1. Geological context

The planning for the KTB super-deep borehole project, known as German Continental Deep Drilling Program, began in 1978 and drilling started in 1987 on the test hole KTB-VB (finished 1989) and the main borehole KTB-HB started in 1990 (finished 1994). In only 4 years KTB-HB was drilled down to a depth of 9101 m (29,859 ft). They made a number of very surprising finds and non finds including water deep below the surface.

Geologists had formed a picture of the crust at the Windischeschenbac site by examining rock outcrops and tow dimensional (2D) seismic measurements (see figure 2).

The lithology model shows an alternating of paragneisses, metabasite, amphibolites and hornblende gneiss and alternating layers of gneisses and amphibolites.

The structural model shows fomrtaions dipping 50° and 75° south-southwest over the first 3000m, followed by a rotation of the dip to the east with a much shallower dip of 25° in the fold hinge. These models were built from interpretation of core data and boreholes images such as Fomration MicroScanner images. The formation appears to have twisted and pilled up [19].

Figure 2.

Simplified lithology and structure revealed by KTB borehole. (a) Alternating layers of metamorphic rocks. (b) Structural model shows formations dipping.

2.3.2. Data processing

The proposed idea has been applied on the KTB synthetic seismogram, figures 3a and 3b present the Sonic Velocity of the Primary wave and the formation density, after that the reflectivity function at normal incidence is calculated, for this last equation 06 is used [24].

$${\mathrm{ }Cr}_i=\frac{{\rho }_{i+1}{Vp}_{i+1}-{\rho }_i{Vp}_i}{{\rho }_{i+1}{Vp}_{i+1}+{\rho }_i{Vp}_i}$$E6

Where: ${\rho }_i$ is the formation density at depth $Z_i$

is the Velocity of the P wave depth $Z_i$

Figure 3 shows this reflectivity versus the depth, note that the processed depth interval is [290m, 4000m], with a sampling interval of 0.1524m.

The synthetic seismic seismogram is than calculated using the convolution model, in fact a seismic trace with an emitted source wavelet $W\left(t\right)$ is given by [24]:

$$\mathrm{ }T\left(t\right)=W\left(t\right)\mathrm{*}Cr\left(t\right)+N\left(t\right)$$E7

N(t) is a noise.

Figure 3.

a)Velocity of the P wave of the pilot KTB borehole versus the depth, (b)Density log of the pilot KTB borehole versus the depth

The emitted wavelet used in this paper is the Ricker, figure 5 shows the graph of this last versus the time. Figure 6 shows the synthetic seismic seismogram calculated using the convolution model. One can remark that is very easy to identify facies limits in this trace.

To check the robustness of the generalized fractal dimensions, the seismic seismogram is now noised with 200% of white noise.

Figure 7 shows the noisy seismic seismogram, on can remark that we are not able to identify the facies limit, for the full depth interval, for example geological formations limited between 2750m and 3000m are hidden by noise. The noisy seismic seismogram is then processed using the proposed algorithm. Figure 8 shows the three fractal dimensions versus the depth compared with the geological map.

Figure 4.

Normal incidence reflectivity versus the depth of the pilot KTB borehole.

Figure 5.

Graph of the Ricker wavelet versus the time.

Figure 6.

Synthetic seismic seismogram of the pilot KTB borehole.

Figure 7.

Noisy seismic seismogram

2.4. Results interpretation and conclusion

One can remark that the principles contacts that exist in geological cross section over the drill are identified by the generalized fractal dimensions (see blue dashed lines in figure 8), however the fractal dimension D₀ (called the capacity dimension) is not always very robust than the information and the correlation dimensions. For example at the depth 460m (first blue dashed line) D₀ is not able to detect this contrast. Otherwise, the generalized fractal dimensions are very robust tools for identification of amplitudes that are due to the lithology variation. For example in the depth interval [2750m, 3000m] the three fractal dimensions have identified many contacts that are totally hidden by noise. The multifractal analysis of seismic data can be used for facies identification, we suggest introducing this last in seismic data processing flow and software.

3. Heterogeneities analysis from well-logs data using the multifractal analysis and the continuous wavelet transform

The fractal analysis has been widely used for heterogeneity analysis from well-llogs data. We cite for example the paper of Kue et al [20]. It shows how texture logs computed from multifractal analysis of dipmeter microresistivity signals can be used for characterizing lithofacies in combination with conventional well logs. Li [21], has analyzed the well-logs data as a Fractional Brownian Motion (fBm) model, he has examined the paradoxical results found in the literature concerning the Hurst exponents. Ouadfeul and Aliouane [08] have published a paper that uses the multifractal analysis for lithofacies segmentation from well-logs data. In this paper, we use the wavelet-based generalized fractal dimensions for heterogeneities analysis. The proposed idea has been applied to a synthetic and real sonic well-logs data of the Kontinentales Tiefbohrprogramm de Bundesreplik Deutschland (KTB) also known as German Continental Deep Drilling Program. We start the paper by describing a wavelet-based multifractal analysis called the Wavelet Tranform Modulus Maxima lines (WTMM) method, after that the processing algorithm of heterogeneities analysis from well-logs data is well detailed. The proposed idea is then applied to the synthetic and real data of the Pilot and the Main KTB boreholes. We finalize the paper by the results discussion and a conclusion.

Figure 8.

.Multifractal Analysis of the noisy seismic seismogram: (a) Seismic seismogram without noise. (b) D0. (c) D1. (d) D2. (e) Geological cross section over the pilot KTB drilling.

3.1. The processing algorithm

The processing algorithm is based on the application of the Wavelet Transform Modulus Maxima lines (WTMM) method at each 128 samples of the signal. A moving window with this last number of samples is used. The window center is moved by 64 samples. At each window the three generalized fractal dimensions D0, D1 and D2 are calculated. Note that 128 samples is the less number of samples where the WTMM can be applied [23]. The analyzing wavelet is the Compex Morlet defined by equation 6:

$V{p}_i$E8

Ouadfeul and Aliouane [8,9], have showed that the optimal value of $$\psi (t)=\exp (i\omega t)\exp (-t²/2)-\sqrt{2}\exp (-\Omega ²/4)\exp (i\omega t)\exp (-t²).$$ for a better estimation of the Hurst exponent is equal to 4.8.

The purpose of this work is to use the generalized fractal dimensions for lithology segmentation and heterogeneities analysis.

3.2. Application to synthetic data

To check the efficiency of the generalized fractal dimensions calculated using the WTMM formalism for boundaries layers delimitation, we have tested this last at a synthetic model formed with four roughness coefficients to model the geological variation of lithology. Each lithology will give a synthetic well-log considered as a fractional Brownian motion (fBm) model. Parameters of the synthetic model are detailed in table 1. Figure 9 shows the generated synthetic model. The generalized fractal dimensions obtained by WTMM analysis of this synthetic well-log are presented in figure 10. We can observe that the fractal dimensions D₁ and D₂ are able to detect exactly the boundaries of each layer; however the capacity dimension D0 is not sensitive to lithology variation. By consequence, D0 cannot be used for lithology segmentation.

Figure 9.

Synthetic well-log data realization with four roughnesses

.

Figure 10.

Generalized fractal dimensions obtained by WTMM analysis of a synthetic well-log data. (a) : D0, (b): D1 (c): D2

.

3.3. Application to KTB boreholes data

The proposed idea is applied on the sonic Primary (P)) wave velocity well-log of the pilot KTB borehole. The goal is to check its efficiency on real well-logs data.

3.3.1. Data processing

Figures 3a show the Primary (P) wave velocity in the Pilot KTB borehole. In this paper, only the depth interval [290m, 4000m] is considered.

Data are recorded with a sampling interval of 0.1524m. The wavelet transform modulus maxima lines method is applied with a sliding window of 128 samples. The window center is moved with 64 samples. The analyzing wavelet is the Complex Morlet (see equation 6). At each window the spectrum of exponents is estimated and the three fractal dimensions are calculated using equation 5. Figure11 shows the three fractal dimensions obtained by multifractal analysis of the Primary wave velocity for the Pilot KTB borehole.

Layer number	01	02	03	04
Thickness (m)	204.80	204.80	204.80	204.80
Roughness Coefficient	0.33	0.51	0.76	0.89
Number of Samples	2048	2048	2048	2048

Table 1.

Parameters of the synthetic well-log data.

3.3.2. Results interpretation and conclusion

Obtained results are compared with the simplified lithology and geological cross section of the area (figure 11). One can remark easily that the capacity dimensions D₀ is not sensitive to the lithology variation. Segmentation models based on D₁ and D₂ are made. It is clear that the information and correlation dimensions are able to detect lithological transitions that are defined by geologists (Black dashed lines in figure 11). Geological model based on the fractal dimensions is proposed (see blue dashed lines in the same figures).

Figure 11.

Generalized fractal dimensions compared with the Geological Cross-Section for the pilot KTB borehole. The analyzed log is the velocity of the P wave.

We have implanted a technique of lithology segmentation based on the generalized fractal dimensions. The proposed idea is successfully applied on the pilot KTB borehole. We recommend application of the proposed method on the potential magnetic and gravity data for contact identification and causative sources characterization. The proposed idea can be applied for reservoir characterization and lithofacies segmentation from well-logs data; it can help for oil exploration and increasing oil recovery.

4. Fractal analysis of 2D seismic data for heterogeneities analysis

The 1D wavelet transform modulus maxima lines WTMM is a multifractal analysis technique based on the summation of the modulus of the continuous wavelet transform(CWT) on its maxima. The obtained function is used to estimate two spectrums, one is the spectrum of exponents and another is the spectrum of singularities. The WTMM was used in various domains to resolve many scientific problems [24][25]. A generalized wavelet transform modulus maxima lines WTMM in the 2D domain are used by many researchers to establish physicals problems [25][26]. The wavelet transform has been applied in seismic image processing; Miao and Moon [27] have published a paper on the analysis of seismic data using the wavelet transform. A New sparse representation of seismic data using adaptive easy-path wavelet transform has been developed by Jianwei et al [28].

The continuous wavelet transform, is used by Pitas et al [29] for texture analysis and segmentation of seismic images. In this paper we process the intercept attribute of 3D synthetic seismic AVO data by the 2D WTMM to establish the problem of heterogeneities. It is very complex and need advanced processing tools to get more ideas about morphology of rocks.

4.1. Brownian fractional motion and synthetic model

For H ∈ (0, 1), a Gaussian process {B_H(t)}_t≥0 is a fractional Brownian Motion if for all t, s ∈  it has [30] :

1. A Mean: E[B_H(t)] = 0

2. A Covariance: E[B_H(t)B_H(s)] = (1/2) {|t|^2H + |s|^2H - | t – s|^2H}

H is the Hurst exponent [31].

We suppose now that we have a geological model of two layers the first is homogonous with the following physical parameters:

1. The density is calculated using the Gardner model [32]: $$Vp\approx 3500\left(\raisebox{1ex}{$m$}\!\left/ \!\raisebox{-1ex}{$s$}\right.\right)$$

2. The velocity of the shear wave is estimated using Castagna Mud-rock line [33]: $$\rho \approx 1.741\times V{P}^{0.25}=2.38\left(\raisebox{1ex}{$g$}\!\left/ \!\raisebox{-1ex}{$cc$}\right.\right)$$

The second is a heterogeneous model with the parameters detailed below.

We suppose that the velocity of the P wave, the velocity of the shear wave and the density of the synthetic model are a Brownian fractional motion model versus the azimuth$$Vs=0.8621XVp-1172.4=1744.95(m/s)$$. $$\theta$$

A generation of the three geological parameters is represented in figures 12.a, 12.b, and 12.c.

We suppose that the three mechanical parameters are in the following range of intervals:

$$0\le \theta \prec 180$$E10

Where Vp(m/s) is the velocity of the P wave, Vs(m/s) is the velocity of the shear wave and $$\begin{array}{l}2000\prec Vp\le 6000\\ 1000\le Vs\le 4000\\ 2.1\le \rho \le 3.\end{array}$$ is the density.

We calculated the reflection coefficients at the null offset R0, which are depending to the azimuth or to X and Y coordinates, obtained results are represented in figure 13.

Figure 12.

Physical parameters of a synthetic heterogonous layer generated randomly versus azimuth (a) Velocity of the P wave (b) Velocity of the S wave (c) Density of the model.

Figure 13.

Response of a model of two layers the first is homogeneous and the second is heterogeneous

4.2. The 2D wavelets transform modulus maxima lines WTMM

The 2D wavelet transform modulus maxima lines WTMM is a signal processing technique introduced by Arneodo and his collaborators in image processing [26][27].It was applied in medical domain as a tool to detect cancer of mammograms and in image processing [25]. The flow chart of this method is detailed in figure 14.

Figure 14.

Flow chart of the 2D wavelet transform modulus maxima lines WTMM

4.3. Application to synthetic model

We applied the proposed technique at the synthetic model proposed above. Figure 15 presents the modulus of the wavelet coefficients at the lower scale a=2.82m. Figure 16 presents the phase proposed by Arneodo et al [31]; the Skelton of the modulus of the continuous wavelet transform is presented in figure 17 and the local Hölder exponents estimated at each point of maxima is presented in figure 18 [31]. One can remark that the Hölder exponents map can provide information about reflection coefficient behavior. So it can be used as a new seismic attribute for lithology analysis and heterogeneities interpretation.

Figure 15.

Modulus of the wavelet transform plotted at the low dilatation

Figure 16.

Phase of the 2D continuous wavelet transform

Figure 17.

Skeleton of the module of the 2D wavelet transform at the scale (a=2.82m)

Figure 18.

Map of local Hölder exponents

4.4. Results interpretation and conclusion

Analysis of the obtained results shows that the 2D WTMM analysis can enhance seismic data interpretation. Hölder exponents map (figure 18) is a good candidate for reservoir heterogeneities analysis. This map is an indicator of geological media roughness.

Conclusion

Application of the proposed technique to an AVO synthetic heterogeneous model shows that the 2D wavelet transform modulus maxima lines (WTMM) is able to give more information about reservoirs rock heterogeneities. The local Hölder exponent can be used as a supplementary seismic attribute to analyze reservoir heterogeneities. The proposed analysis can help the hydrocarbon trapping research and analysis, fracture detection and fractured reservoirs analysis. We suggest application of the proposed philosophy at real seismic AVO data and its attribute, for example the Intercept, The gradient, the Fluid Factor proposed by Smith and Gildow [34] and the attribute product Intercept*Gradient.

5. Conclusion

As a conclusion of the current chapter, we can say that the fractal analysis has showed it’s powerful to help exploration geophysics. Application to synthetic and real seismic and well-logs data shows clearly that this analysis of the earth response considered as a chaotic system can help geosciences.