Two-Dimensional Multifractional Brownian Motion- Based Investigation of Heterogeneities from a Core Image

A core sample is a cylindrical section obtained by driving a hollow tube into the undisturbed medium and withdrawing it with its content. In practice, the sample is pushed more or less unbroken into the tube. Once removed from the tube in the laboratory, it is analyzed by different techniques and equipment depending on the desired type of data. The hole made for the core sample is called the "core hole". A variety of core samplers exist to sample different media under diverse conditions. For instance, sediments or rocks are sampled with a hollow steel tube called a core drill. A scientific coring has been used in the first time for sampling the ocean floor. Then, it is soon exploited to analyze lakes, ice, mud, soil and wood. Cores provide precious information about the evolution of climate, species and sedimentary composition during geologic history. In petroleum engineering, core analysis presents a way of measuring well conditions downhole by studying samples of reservoir rocks. It gives the most accurate estimations of porosity, permeability, fluid saturation and grain density. These measurements help to understand the conditions of the well and its potential productivity. In addition to the basic petrophysical properties estimated from the core, a special core analysis can be undertaken in order to determine permeability, wettability, capillary pressure, and electrical properties. Petrographic studies and sieve analysis can also be carried out in such analysis. In recent years, numerical analysis has been widely used for the investigation of images, since it yields results more objective and reliable than those obtained by conventional methods based on human observations. Fractal analysis has been introduced to examine images texture (Bourissou et al., 1994; Levy-Vehel and Mignot, 1994; Liu and Li, 1997; LevyVehel, 1995, 1997, 1998; Pesquet-Popescu and Levy-Vehel, 2002; Malladi et al., 2003; Tahiri et al., 2005). In this study, we suggest to go beyond the conventional core analysis, and to perform a new approach to extract the maximum features from a core image using a fractal analysis. The conventional fractal model used previously in image processing, the two-dimensional fractional Brownian motion (2DfBm), presents a constant Holder function H, thus does not allow to explore the spatial evolution of the local regularity. To do so, we suggest to


Introduction
A core sample is a cylindrical section obtained by driving a hollow tube into the undisturbed medium and withdrawing it with its content.In practice, the sample is pushed more or less unbroken into the tube.Once removed from the tube in the laboratory, it is analyzed by different techniques and equipment depending on the desired type of data.The hole made for the core sample is called the "core hole".A variety of core samplers exist to sample different media under diverse conditions.For instance, sediments or rocks are sampled with a hollow steel tube called a core drill.A scientific coring has been used in the first time for sampling the ocean floor.Then, it is soon exploited to analyze lakes, ice, mud, soil and wood.Cores provide precious information about the evolution of climate, species and sedimentary composition during geologic history.In petroleum engineering, core analysis presents a way of measuring well conditions downhole by studying samples of reservoir rocks.It gives the most accurate estimations of porosity, permeability, fluid saturation and grain density.These measurements help to understand the conditions of the well and its potential productivity.In addition to the basic petrophysical properties estimated from the core, a special core analysis can be undertaken in order to determine permeability, wettability, capillary pressure, and electrical properties.Petrographic studies and sieve analysis can also be carried out in such analysis.In recent years, numerical analysis has been widely used for the investigation of images, since it yields results more objective and reliable than those obtained by conventional methods based on human observations.Fractal analysis has been introduced to examine images texture (Bourissou et al., 1994;Lévy-Véhel and Mignot, 1994;Liu and Li, 1997;Lévy-Véhél, 1995, 1997, 1998;Pesquet-Popescu and Lévy-Véhel, 2002;Malladi et al., 2003;Tahiri et al., 2005).In this study, we suggest to go beyond the conventional core analysis, and to perform a new approach to extract the maximum features from a core image using a fractal analysis.The conventional fractal model used previously in image processing, the two-dimensional fractional Brownian motion (2D-fBm), presents a constant Hölder function H, thus does not allow to explore the spatial evolution of the local regularity.To do so, we suggest to consider a generalized fractal model, the two-dimensional multifractional Brownian motion (2D-mBm), which presents a regularity varying in space.The mBm model, initially proposed by Peltier andLévy-Véhél (1995), andBenassi et al. (1997), is used in many disciplines: images processing (Bicego and Trudda, 2010), traffic phenomena (Li et al., 2007), geophysics (Wanliss, 2005;Wanliss and Cersosimo, 2006;Cersosimo and Wanliss, 2007;Gaci et al., 2010;Gaci andZaourar, 2010, 2011).For the estimation of the mBm processes' local regularity, we propose three algorithms based on the two-dimensional continuous wavelet transform (2D-CWT).The wavelet coefficients are calculated by Fast Fourier Transform (FFT) using the Morlet wavelet and the Mexican hat for the first and the second algorithms, respectively.However for the third algorithm, the coefficients estimation is carried out using the multiple filter technique (2D MFT) that we generalized to two dimensions (Gaci, 2011), from the one-dimensional case (1D MFT) (Li, 1997;Gaci et al., 2011).This chapter is organized as follows.First, we give a brief theory on 2D-mBm model and the wavelet-based estimators of the local regularity.The potential of the suggested algorithms is then demonstrated on synthetic 2D-mBm paths.The results showed that the 2D MFT algorithm yields the best Hölder exponent estimates.Next, the suggested regularity analysis is extended to digitalized image data of a core extracted from an Algerian borehole.It is shown that the data exhibit a fractal behavior.In addition, the derived regularity maps, obtained with the 2D MFT algorithm, show a strong correlation with the core heterogeneities.

(Multi)fractional Brownian motion 2.1 Fractional Brownian motion
Fractional Brownian motion (fBm) is one of the most popular stochastic fractal models for studying rough signals.It was introduced by Kolmogorov (1940) and studied by Mandelbrot and Van Ness (1968).A fBm, denoted by B H (t), is the zero-mean Gaussian process with stationary increments.It is parameterized by a constant Hurst parameter H.The fBm is H-self affine, i.e.: Where  means the equality of all its finite-dimensional probability distributions.
The bidimensional isotropic fractional Brownian motion, or Lévy Brownian fractional field, with Hurst parameter H is a centered Gaussian field H B with an autocorrelation function (Kamont, 1996): and .
is the usual Euclidian norm.
For H=1/2, the fractional Brownian motion is reduced to a Wiener process.
The regularity of the 2D-fBm is measured by the pointwise Hölder exponent Indeed, it is shown that almost surely: . Therefore, the higher the H value, the smoother the 2D-fBm paths.

www.intechopen.com
Two-Dimensional Multifractional Brownian Motion-Based Investigation of Heterogeneities from a Core Image 111

Multifractional Brownian motion
Multifractional Brownian motion (mBm) was introduced by Peltier and Lévy-Véhél (1995), and Benassi et al. (1997) by allowing H to vary over time.Even if no longer stationary nor self-similar compared to the fBm, the mBm presents the advantage to be very flexible since the function H(t) can model phenomena whose sample paths display a time changing regularity.
For a continuous function , the isotropic multifractional Brownian field Identically to the 2D-fBm, the local regularity of the 2D-mBm paths is measured by means of the pointwise Hölder exponent.For a differentiable function H, the relation

Two-dimensional continuous wavelet transform
The two-dimensional continuous wavelet transform (2D-CWT) of a signal   , sxy is given by (Chui, 1992;Holschneider, 1995): where   , g x y is the analyzing wavelet, "a" is the scale parameter, "b x " and "b y " are the respective translations according to X-axis and Y-axis.The symbol " -" denotes the complex conjugate.
x y x y Pkxy P xy k k Where is the local spectral exponent which is related to the local Hurst (or Hölder) exponent,

 
, x y is computed as the slope of the scalogram versus the wavenumber in the log-log plan, the   , Hxy value is then derived using the equation (9).

Used analyzing wavelets
The analyzing wavelets used in this application are the Morlet wavelet and the Mexican hat (Fig. 1).This choice is motivated by their adequate properties for the regularity analysis. - The Morlet wavelet:

Wavelet-based estimators of the local regularity
As explained earlier, the computation of the local Hölder exponents H(x,y) requires to calculate the two-dimensional wavelet continuous transform.Here, we suggest three algorithms for the implementation of the 2D-CWT, which are:

FFT-based algorithms
These algorithms are based on the property that wavelet coefficients, expressed by the equation ( 4), can be performed via the Fourier transform using the Morlet wavelet and the Mexican hat: ) ., These algorithms are accurate but slow, and the signal length must be a power of 2.

Generalized two-dimensional multiple filter technique
The 1D MFT was initially suggested by Dziewonski et al. (1969).It consists on carrying out a decomposition using a Gaussian filter: 2 (, ) where n k is a variable center angular frequency (or wavenumber) of the filter (, ) n Gkk , and  is a shaping parameter of the filter.
In order to overcome the poor ''time'' and ''low-frequency'' domains resolution toward the low frequencies, Li (1997) suggests a varying quality factor Q and a varying bandwidth Δk: Where  is a constant, k 1 and k 2 are the -3 dB points of the Gaussian filter.

Application to simulated 2D-mBm paths
In this section, the suggested estimators of the local regularity are tested on synthetic 2D-mBm paths whose lengths are 256 x 256, generated using the kriging method (Barrière, 2007).The regularity functions and the simulated 2D-mBm paths corresponding to the three theoretical H functions are presented in Figure 2. The larger H value, the smoother the modeled surface.
Using the three algorithms, we have estimated H maps.For the first wavelet-based algorithm, we use the Morlet wavelet with  0x = 0y =8.9443, while for 2D MFT, the selected parameters for the two-dimensional Gaussian filters are: - The shaping factor =40, -The minimal central wavenumber of the filter min = min =2.10 -3 rad/m, - The maximal central wavenumber of the filter max = max =2.10 0.5 rad/m, - The number of the central wavenumbers of the filter N=100, - The sampling rate is selected as 0.1524m.The H maps obtained by the three estimators, presented in Figure 3, show that the regularity estimated by the first wavelet-based algorithm using the Morlet wavelet are better than that calculated by the second algorithm with the Mexican hat.In addition, the suggested 2D MFT provides the best estimations of the regularity maps with the least errors.For this reason, we retain only this estimator in the following.It can be also remarked that all the used algorithms yield large absolute values of the estimation error in the limits of the analyzed 2D-mBms.

Application to digitalized core image data
Here, the local regularity analysis is performed on digitalized core image data.The analyzed core is extracted from a well drilled in an Algerian basin.It is chosen since it represents the main geological features of the studied region (Fig. 4).
The core presents medium to fine quartzitic sandstone, clay and quartz cemented cross stratifications underlined by mud films.The formation is affected by a main fracture F1 with a high angle dip (≈75°) filled with quartz.It is also noted the presence of another fracture F2 sub-parallel to F1 but less important.The processing of the core image requires first its digitalization.The core image is digitalized and codified in gray levels with integer values ranged between 0 and 255.The obtained digitalized core image, illustrated by figure 5, corresponds to a matrix of 3642 x 996 with a sampling rate Δx=Δy=0.0121cm.First, the fractal behavior of the digitalized core image data is inspected on five horizontal and five vertical profiles as shown in Figure 5.For each profile, the Fourier amplitude spectrum is computed and represented in a double logarithmic scale.Then, we estimate the local Hölder function H(x) using an algorithm based on the local growth of the increment process S k,n (i) (Peltier andLévy-Véhél, 1994, 1995;Muniandy et al., 2001;Li et al., 2008;Gaci et al., 2010): where n is the signal X length, k is a fixed window size, and m is the largest integer not exceeding n/k.The local Hölder function H(x) at point is given by The obtained results corresponding to the horizontal profiles and the vertical profiles are respectively exposed in figures 6a and 6b.
It can be noted that all the resulted amplitude spectra exhibit an algebraic decay; that illustrates the fractal properties of the digitalized data.Besides, the analyzed profiles present a varying regularity with the position according to X-and Y-axis.They can be then regarded as paths of multifractional Brownian motions (mBms) (Peltier and Lévy-Véhel, 1995).The variation of H exponent value is related to the local lithological changes of the core composition.The next step consists on establishing regularity maps from the digitalized data using the 2D MFT algorithm.The implementation of the latter algorithm requires the ''reconditioning'' of the data so that the matrix dimensions corresponding to the digitalized data are a power of 2. For the purpose of processing the digitalized data, and considering the limitations of the available computer's capacities, we have splited the obtained matrix (3642 x 996) into two overlapping sub-matrixes whose size is 2048 x 1024.The sub-matrixes are padded by zeros so that their dimensions following Y-axis, initially equal to 996, reach 1024.The parameters selected for the 2D MFT are as follow: -The minimal center wavenumbers of the filter : min = min =2.10 -1 rad/m; - The maximal center wavenumbers of the filter: The other parameters ( and N) are similar to those used in the previous section.
The final regularity map is constructed from the H sub-maps related to the two sub-matrixes (Fig. 7).The H values in the overlapping zone are calculated as the average of the H values corresponding to the H values in the sub-maps.From figure 8, it can be seen that the histograms of H values calculated by 2D MFT follow a normal distribution.For each class, the statistical parameters (mean and standard-deviation) are estimated from the histograms of the gray level values, and the corresponding H values (Table 1).It is worth noting that for the six classes, the statistical parameters of H values, estimated from the histograms, present very close values.These results show that the Hölder exponent value can not characterize a geological facies represented by the gray level, whereas its local variation reflects local lithological changes as explained earlier.

Conclusion
In this study, the 2D-mBm has been successfully used for a local Hölder regularity-based modeling of the core image.We have presented three methods for estimating the local regularity.The first and the second ones are FFT-based algorithms using, respectively, the Morlet wavelet and the Mexican hat, while the third method is obtained by extending the one-dimensional multiple filter technique to 2 dimensions (2D MFT).The application of these methods on synthetic 2D-mBm paths showed that the 2D MFT yields the best estimations of the H functions.The analysis of profiles extracted from the digitalized core image data reveals a fractal behavior.Furthermore, the regularity maps obtained by 2D MFT from the digitalized data can characterize heterogeneities from the analyzed core.Although a Hölder exponent value does not describe a specific geological facies, its local variation reflects the lithological changes (faults, breaks, stratifications, etc.).The presented analysis must be undertaken on a large number of cores in order to establish a relation between a geological facies, the corresponding gray level and H values.
exponent H and a positive factor . If s is a stochastic process, the two sides of the relation follow the same law.
www.intechopen.comBytaking the scale "a" inversely proportional to the wavenumber k:1 ak  , the wavelet coefficients will be expressed in (k, b x , b y ) plane.The scalogram can be defined as the square of the amplitude spectrum: xy

Fig. 3 .
Fig. 3. Regularity functions obtained by the three algorithms from the 2D-mBm paths, represented in Fig.2, simulated with three types of H function: (a) bilinear (b) logistic (c) periodic Line 1: theoretical H function; Lines 2 & 3: regularity functions estimated using FFT-based algorithms with, respectively, the Morlet wavelet and the Mexican hat; Line 4: regularity function estimated using 2D MFT (=40).The columns from left to right represent: (1) the estimated H function , (2) the estimation error , calculated as the difference between the estimated H function and the theoretical H function.

Fig. 5 .
Fig. 5. Positions of the horizontal and the vertical profiles on the digitalized core image.The horizontal profiles numbered from 1 to 5 (in blue) correspond to the respective positions y = 0.0120 m, 0.0361 m, 0.0603 m, 0.0845 m and 0.1087 m, while the vertical profiles numbered from 1 to 5 (in red) correspond to the respective positions x = 0.0241 m, 0.1208 m, 0.2175 m, 0.3142 m and 0.4109 m.

Fig. 6 .
Fig. 6.Investigation of the fractal properties of the five horizontal (a) and vertical (b) profiles extracted from the digitalized core image.The five lines in (a) (resp.(b)), from top to bottom, correspond to the respective horizontal (resp.vertical) profiles numbered from 1 to 5. Left: the profile of the digitalized core image data, middle: the amplitude spectrum module of the data with respect to the wavenumber in the log-log scale, right: the local H exponent.

Fig. 7 .
Fig. 7.A regularity map (b) obtained by 2D MFT from the digitalized core image data (a).The regularity map (b), corresponding to the data of the whole core, is obtained from the regularity maps (c) and (d), related to the two ''sub-zones'' of the core image.From Figure7, it can be seen that the analyzed data present a varying regularity in the XY plan.It is again confirmed that the digitalized core image data can be modeled as a

Fig. 8 .
Fig. 8. Histograms of the digitalized data values extracted from the core image, and the corresponding H values estimated by 2D MFT, for the six classes.The six columns from left to right correspond respectively to the classes 1 to 6.
to 2 dimensions.The idea consists on decomposing the two-dimensional signal using a Gaussian filter(, , ) Three types of Hölder function H are chosen:

Table 1 .
Statistical parameters estimated from the histograms (gray level and corresponding H values) related to the six classes.