Local regularity analysis of strata heterogeneities from sonic logs

Borehole logs provide geological information about the rocks crossed by the wells. Several properties of rocks can be interpreted in terms of lithology, type and quantity of the fluid filling the pores and fractures. Here, the logs are assumed to be nonhomogeneous Brownian motions (nhBms) which are generalized fractional Brownian motions (fBms) indexed by depth-dependent Hurst parameters H(z). Three techniques, the local wavelet approach (LWA), the average-local wavelet approach (ALWA), and Peltier Algorithm (PA), are suggested to estimate the Hurst functions (or the regularity profiles) from the logs. First, two synthetic sonic logs with different parameters, shaped by the successive random additions (SRA) algorithm, are used to demonstrate the potential of the proposed methods. The obtained Hurst functions are close to the theoretical Hurst functions. Besides, the transitions between the modeled layers are marked by Hurst values discontinuities. It is also shown that PA leads to the best Hurst value estimations. Second, we investigate the multifractional property of sonic logs data recorded at two scientific deep boreholes: the pilot hole VB and the ultra deep main hole HB, drilled for the German Continental Deep Drilling Program (KTB). All the regularity profiles independently obtained for the logs provide a clear correlation with lithology, and from each regularity profile, we derive a similar segmentation in terms of lithological units. The lithological discontinuities (strata’ bounds and faults contacts) are located at the local extrema of the Hurst functions. Moreover, the regularity profiles are compared with the KTB estimated porosity logs, showing a significant relation between the local extrema of the Hurst functions and the fluid-filled fractures. The Hurst function may then constitute a tool to characterize underground heterogeneities. Correspondence to: S. Gaci (saidgaci@yahoo.com)


Introduction
Well logs allow to directly measure physical properties of the earth, and hence provide information exhibiting more reliability and high resolution than surface seismic data.Since the number of wells is limited, it is important to model and analyze the logs in order to extract maximum features of the subsurface.
The fractal behavior of porosity well logs was first noticed by Hewett (1986).Then, Hardy and Beier (1994) found that well logs can be conveniently considered as fractional Brownian motions (fBms).These non-stationary self-affine fractal processes are characterized by a fractal k −β -power spectrum model where k is the wavenumber and the parameter β is related to the Hurst parameter H .Moreover, fBms are monofractal whose complexity is defined by a single global coefficient H , associated to the fractal dimension D, and closely related to the regularity Hölder degree.
Indeed, the Hurst exponent H describes well logs behavior.It is used to know whether the process is persistent ( 1 / 2 <H<1) or anti-persistent (0<H< 1 / 2 ).A persistent process does not show unexpected rapid changes (long memory), i.e. deviations tend to keep the same sign.The limiting case H = 1, reflects X(t) ∝ t, a smooth signal.While an anti-persistent process may exhibit sudden variations, that is, deviations of one sign are generally followed by deviations with the opposite sign.The limiting case H = 0, (short memory), corresponds to a white noise, where fluctuations at all frequencies are equally present.When H is close to 1, the signal is weakly irregular, while for H close to 0, it is highly irregular.The coefficient H can be estimated by various methods, the most popular is the Rescaled-Range (R/S) analysis (Dolan and Bean, 1997;Leonardi andKümpel, 1998, 1999;Malamud and Turcotte, 1999;Li, 2003;Arizabalo et al., 2004Arizabalo et al., , 2006;;Bansal and Dimri, 2005;Chamoli et al., 2007).In addition to the R/S analysis, wavelet-based estimators have been used very successfully for estimating Published by Copernicus Publications on behalf of the European Geosciences Union and the American Geophysical Union.
In this study, a nhBm model is assumed for sonic logs, and the singular behavior is explored using the multiscale wavelet analysis and Peltier Algorithm (PA).The former is based on two approaches: the local wavelet approach (LWA) and the average-local wavelet approach (ALWA).The local spectral exponents derived from LWA provide a singularity pattern at the finest detail of the signal, while using the ALWA, the average-local scaling parameter estimated within a moving window on the scalogram, reflects a global effective regularity linked to the main structural phenomena.On the other hand, PA is based on the local growth of the logs' increment process, and allows for the estimation of the local Hurst parameter of the sonic log from its stochastic component.
This paper is structured as follows.Firstly, a short mathematical description of the methods used for the determination of the local Hurst function is given.Secondly, two synthetic sonic logs with different parameters are generated using the successive random additions (SRA) method.The simulated logs are then used to demonstrate the potential of the suggested analysis techniques and the results are compared.Furthermore, the three techniques are implemented to study multifractionality properties of the P-and S-waves sonic measurements recorded in the KTB pilot and main boreholes.The KTB offers a unique opportunity to investigate the properties of crustal heterogeneities and many studies have been devoted to analyze the well log measurements (Leary, 1991;Maus and Dimri, 1994;Wu et al., 1994;Kneib, 1995;Holliger, 1996;Jones and Holliger, 1997;Li, 1998;Dolan et al., 1998;Zhou and Thybo, 1998;Leonardi andKümpel, 1998, 1999;Marsan andBean, 1999, 2003;Goff and Holliger, 1999;Fedi, 2003;Fedi et al., 2005;Gaci and Zaourar, 2010).

Homogeneous (or fractional) Brownian motion
The fractional Brownian motion (fBm) is the most familiar and classical self-affine model used for well log modeling (Li, 1997;Li and Ulrych, 1999).Developed by Mandelbrot and Van Ness (1968), the fBm B H = {B H (z),z > 0} is a Gaussian process with zero mean, stationary increments and B H (0) = 0. Its covariance is given by: where 0 < H < 1, σ >0, s >0 and E[.] denotes the expected value.
For H = 1/2, the process B H (t) corresponds to the standard Brownian motion.
In this case, the increments For other values of H (H = 1/2), the increments Xz correspond to an integrated fractional white noise characterized by a long-range correlation in the sense that the correlation decay is sub-exponentially (i.e.algebraically): One of the main features of fBm process is that its pointwise Hölder regularity is a function of its Hurst parameter.For every point, its Hölder exponent is almost surely equal to H (Ayache et al., 2005).Recall that, for a stochastic process X, the Hölder exponent at a point z, is defined as: This coefficient reveals the singularity strength of the signal.

Nonhomogeneous Brownian motion
A fBm exhibits a constant Hölder exponent value along its trajectory, and it has stationary-increments property.However, most real phenomena do not satisfy these properties, and need to consider complex processes with non stationary increments of any order.For this reason, nonhomogeneous Brownian motion (nhBm) was introduced (Peltier andLévy-Véhel, 1994, 1995).It is a generalized fBm obtained by replacing Hurst parameter with a function H (z).
The nhBm is a zero mean Gaussian process, whose increments are in general neither independent nor stationary.
Contrarily to fBm, the pointwise Hölder regularity of a nhBm may depend on the location.At each point z 0 , the pointwise Hölder exponent is almost surely equal to its Hurst exponent H (z 0 ) (Peltier and Lévy-Véhel, 1995;Benassi et al., 1997).Thus, a nhBm allows to describe phenomena whose regularity varies in space (or in time).
It is shown that a nhBm is locally asymptotically selfsimilar (lass): at each z, there exists an fBm with exponent H (z) which is "tangent" to the nhBm.In other terms, a path of a nhBm can be seen as the lumping of infinitesimal portions of fBms with well-chosen exponents.

Wavelet-based estimators
The continuous wavelet transform (CWT) of a signal s(z) is given by (Goupillaud et al., 1985): where g(z) is the analyzing wavelet, "a" is a scale parameter and "b" is translation (the overlain symbol "-" denotes the complex conjugate).
The analyzing wavelet must possess a zero mean and a good localization in the depth-wavenumbers domains.Moreover, for the purpose of a multifractal analysis, it is also required to have n vanishing moments (i.e., to be orthogonal to some low order polynomials, up to the degree n-1): Considering this property, the signal does not need a detrending operation in order to study its stochastic component (Malamud and Turcotte, 1999).In the present work, the choice of Morlet wavelet is motivated by its best simultaneous localisation in depth and in wavenumber, and its large number of vanishing moments.The wavelet transform is a well adapted tool for studying scaling processes.Indeed, the scale invariance property can be made evident provided that the analyzing wavelet decreases quickly enough to zero and has enough vanishing moments (Holschneider, 1995): As already mentioned, the Hölder exponent α(b 0 ) at any position b 0 is almost surely equal to the Hurst parameter H .This result establishes the statistical scale invariance of the CWT of fractal signals, and constitutes the main property of the CWT used in this study.Since the scale "a" is inversely proportional to the wavenumber k: the wavelet transform of the signal can be interpreted in wavenumber-depth (k − z) space instead of a − b space.The scalogram, P (k,z), can be defined as the square of the amplitude spectrum |S(k,z)| 2 : for a fBm, it can be shown that for sufficiently large wavenumbers, it behaves like (Flandrin, 1992): where is the local spectral exponent which depends on the local Hurst parameter H (z) (or on the local Hölder exponent).This relationship is used later to estimate the Hurst function (or the regularity profile).
In this work, we will consider two approaches for the multiscale regularity analysis of the signal:

Local wavelet approach
It consists in estimating the spectral exponent β(z 0 ) at each point z 0 , by estimating the log-log slope of the scalogram P (k,z) with respect to the wavenumber (Li, 1998).The local Hurst parameter is then evaluated according to Eq. ( 11).

Average-local wavelet approach
In addition to LWA, we propose another approach called ALWA.At each position z 0 , the average-local scalogram P ave (k,z 0 ) is calculated by averaging the scalogram values within a moving L-length-window centered in z 0 : Then, we estimate the average-spectral exponent represented by the log-log slope of the average-local scalogram with respect to the wavenumber.As previously, the Hurst parameter is evaluated from the average-spectral exponent using Eq. ( 11).The obtained Hurst curve is smoother and represents a more stable evaluation of the local regularity compared to the previous one resulting by LWA.However, sharp variations are smoothed out.

Peltier Algorithm
The estimation of the local Hurst parameter of a nhBm process {X(j ),0 ≤ j ≤ n} of (n+1) samples is given by Peltier Algorithm (PA) (Peltier and Lévy-Véhel, 1995): www.nonlin-processes-geophys.net/17/455/2010/ Nonlin.Processes Geophys., 17, 455-466, 2010 where S k,n (i) is the local growth of the increment process where k is a fixed window size and m is the largest integer not exceeding n/k.
In contrast to the wavelet-based estimators, this algorithm can not be implemented on the real log data without trends removal.
In this study, the global trends removed from the velocity logs V (z) are considered linear: where V 0 and V 1 are calculated by least-squares fitting procedures.
The analyzed part is the stochastic component s(z), also called fractional fluctuation (Shiomi et al., 1997;Zaourar et al., 2006a), given by: where z is the depth, V 0 and V 1 are coefficients determined by linear regression (least-squares fitting procedures).
It is worth noting that the trend removal and normalization (dividing by the trend) operations are required to get more accurate estimates of Hurst values as explained by Li (2003).The application of only one of these operations separately leads to negative estimates of Hurst values.

Application to synthetic sonic data
First, the synthetic logs are used to show the potential of the three suggested methods LWA, ALWA and PA.The logs are modeled as nhBm processes using the approach suggested by  Peltier and Lévy-Véhel (1995).A nhBm path is generated by simulating N fBms with the parameters H (i), i = 1,...,N, and constructed by setting: Here, the fBms needed to generate a nhBm are modeled by the SRA algorithm (Turcotte, 1997).The latter considers the depth interval 0 ≤ z ≤ 1.First, random values are generated based on a Gaussian probability distribution with zero mean and unit variance, and three of these are placed at abscissa 0, 1 / 2 and 1.Then, the values at the midpoints of above points are linearly interpolated to give five values.Five random values are generated with zero mean and a reduced variance V = 1 n 2H (where n is the number of values generated by interpolation) and added to the above five values.Again these resulted five values are interpolated at midpoints and the same process is repeated until the desired number of points is obtained.
The depth sampling interval is selected as 0.1524 m (the interval used by the KTB sonic measurements; see below).For the purpose of assessing the accuracy of the three methods, we have generated 1000 realizations of a four-layer geological model whose parameters are given in Table 1.For each realization, the Hurst curves are estimated using the different techniques.At the position z = 121.77m, which is located in the second layer characterized by a theoretical Hurst value of 0.4, we established histograms of Hurst values obtained by the three methods (Fig. 1).
For all the methods, the histograms exhibit normal distributions centered at the theoretical H-value with the smallest standard deviations values for PA.For the waveletestimators, the standard deviations for ALWA are less than for LWA (Table 2).On the other hand, the biases of the estimators, calculated as the difference between the Hurst values mean and the theoretical Hurst value, are classified Again, the robustness of PA is checked on a synthetic log corresponding to a multi-layer geological model (Fig. 2a) whose parameters are defined in Table 3. From Fig. 2b, it can be seen that the regularity profile obtained by PA allows for the identification of the twelve geological layers composing the modeled log.The transition between two contiguous layers is marked by a H-value jump discontinuity.On the other hand, the estimated H (z) curve seems to globally coincide with the theoretical curve except for a limited number of depth intervals.

Sonic measurements at the KTB boreholes
This section is devoted to the analysis of sonic log data recorded at two scientific deep boreholes: the so called pilot hole (VB, 4 km) and the ultra deep main hole (HB, 9.1 km), drilled for the German Contienental Deep Drilling Program (KTB).The drilling site is located in south-eastern Bavaria (Oberpfalz), about 200 km NNE of Munich (Germany).The two KTB wells, separated by 200 m, are largely cased and accessible to nearly bottom hole.The wells go through the crystalline metamorphic rocks of a Hercynian continental collision zone.The drilled section is dominated by paragneisses and metabasites, and shows the presence of many faults, cataclastic shear zones, aplitic and lamprophyric dikes (Weber and Vollbrecht, 1989).
In the present paper, the whole sonic log data of the pilot borehole (VB) recorded in the depth interval (28.194-3990.137m) is used.Whereas, for the main borehole (HB), only the data measured in the depth interval 290.017-4509.97m is considered.Both logs are recorded with the same sampling interval (0.1524 m).For the sake of clarity, the following notations: VB-Vp, VB-Vs, HB-Vp and HB-Vs are used to indicate the P-wave velocity (Vp), the S-wave velocity (Vs) recorded respectively at the pilot borehole (VB) and the main borehole (HB).The least-squares fitting coefficients corresponding to the analyzed sonic logs are summarized in Table 4.For each sonic log data, the depth-Hurst function H (z)is estimated by the three methods: LWA, ALWA and PA, and is noted respectively as: H LWA (z), H LWA (z) and H PA (z).
Then, using a lithologic sketch of the study area, a correlation between the regularity profiles H (z) and the lithology of the strata is established.The geological formations crossed by the wells in the studied depth intervals are: amphibolite units (AU), gneiss units (GU), and variegated units (VU).

Results and discussion
The regularity profiles obtained from the sonic logs recorded at the pilot and the main boreholes are presented, respectively, in Figs. 3 and 4. Considering a geological section crossing the KTB wells, a lithological segmentation is carried out on the estimated Hurst curves.The variation of the resulted Hurst values with depth allowed to put in evidence the following lithological units: AU1-AU2-AU3 (for amphibolite units), GU1-GU2-GU3 (for gneiss units), and VU1-VU2 (for variegated units), which are delimited by a gray line.In addition, five (resp.six) faults: F1-F5 (resp.F1-F6), shown in blue line, are identified, respectively, in Figs. 3 and  4.
From these figures, it can be noticed that the Hurst functions H LWA (z) seem to be rough which makes them difficult to interpret; thus the local lithological changes and the different faults contacts are not clearly detected.This is why our analysis is focused on the Hurst functions obtained by the other methods (H ALWA (z) and H PA (z)).In order to give a statistical significance to the changes of the Hurst coefficient value, the Hurst functions are presented with statistical error bars.
The Hurst functions H ALWA (z) presented in Figs.3d1, d2, and 4d1, d2, are smoother than H LWA (z), but both Hurst functions exhibit a similar behavior.However, the H PA (z) curves are characterized by a slightly different behavior and a high smoothness level compared to the previous ones (Figs. 3e1,e2,and 4e1,e2).The Hurst values obtained from all the analyzed velocity logs are generally less than 0.5.
The established lithological subdivision showed that the layers' boundaries correspond to the local extrema (maxima or minima) of the Hurst function even if they are not noticeable on the sonic logs.All the amphibolite units are well identified on the Hurst curves estimated by ALWA and PA for all the logs.Their bounds correspond to variations of the seismic wave velocity.However, for the gneiss units, only GU1 is identifiable on the regularity profiles and the velocity logs.The GU2 is not put in evidence on the Hurst function of HB-Vp log estimated using PA, while GU3 is not detected on the regularity profiles of HB-Vp and VB-Vs logs obtained, respectively, by ALWA, and by both methods (PA and ALWA).This unit is also not recognized on the regularity profiles of VB-Vp and HB-Vs logs estimated by PA.It is worth noting that the limits of GU2 and GU3 do not correspond to any abrupt changes on the velocity logs.Regarding the variegated units, these are noticeable on the regularity profiles, but not on the velocity logs.The subunit VU1 is identified on all the regularity profiles, while VU2 is recognized on the Hurst functions obtained by ALWA from VB-Vp, VB-Vs and HB-Vs logs.The latter unit is also put in evidence by ALWA and PA on the regularity profiles corresponding to HB-Vp log.
The main problem met during the lithological segmentation is the presence of some peaks and lows of the regularity profile in between the geological subunits.These values originated from the fluid-filled fractures because the seismic wave velocity behavior with depth indicated the dominating effect of the fluid-filled fractures, cataclastic fracture zones and heterogeneity (Harjes et al., 1997;Pechnig et al., 1997;Baisch and Harjes, 2003;Rabbel et al., 2004).The Hurst values are compared with the porosity values estimated by Pechnig et al. (1997) at the pilot and the main boreholes as shown in Figs. 5 and 6.The authors interpreted porosities in the range of the calibration data (<4%) as true matrix porosities whereas higher porosities are considered as indications of fluid-filled fracture zones.It is established that the local extrema values of the regularity profiles correspond to the fractures zones characterized by high values of porosity.
As mentioned above, the amphibolites units are completely identifiable on the analyzed velocity logs, in contrast with Gneiss and Variegated Units.Only the GU1 unit is reflected by the velocity change on all the sonic logs.This point may be explained by the fact that amphibolites exhibit faster velocities (Vp and Vs) than gneiss and variegated units, characterized by close velocity values which make them difficult to differentiate on the sonic logs (Goff and Holliger, 1999).
It is important to outline that the faults (F1, F4, F5 and F6) are identified on the Hurst functions .Their locations correspond generally to the local minima of the regularity profile.It can also be observed that each fault is accompanied by abrupt velocity changes except for the fault F4 on the Hurst logs obtained from VB sonic logs.
In addition, the faults F2 and F3 do not correspond to any change of velocity but they are put in evidence on some regularity profiles.The fault F2 is identified on the Hurst functions estimated from HB-Vp and HB-Vs logs by ALWA, and on the Hurst functions estimated from VB-Vp and VB-Vs logs using ALWA and PA.On the other hand, the fault F3 is outlined by ALWA on the Hurst logs estimated from VB-Vs and HB-Vp logs, and by ALWA and PA from VB-Vp log.This result may be explained by the fact that these faults (F2 and F3) are located completely within the gneiss unit GU2.For recall, our algorithms deal with the regularity of the signal and not with its amplitude variation.Then, each regularity fluctuation of the velocity is captured by the Hurst exponent, and constitutes the efficiency of detecting some lithological heterogeneities not identified on log measurements, such as the case of the faults F2 and F3.Now, we illustrate the relationship between the obtained Hurst values and the lithology by investigating H-value mean versus the three methods for each lithological unit.The amphibolites are characterized by H-values mean varying from 0.185 to 0.379 for LWA, from 0.082 to 0.297 for ALWA and from 0.398 to 0.478 for PA.However Fedi et al. (2005) estimated H-values for this lithological unit between 0.25 and 0.3 for density and magnetic susceptibility logs, and between 0.45 and 0.7 for self potential and resistivity logs.
Finally, the VU are described by H-values mean ranges between 0.164 and 0.348 for LWA, 0.073 and 0.274 for ALWA, 0.367 and 0.485 for PA.While Fedi et al. (2005) estimated Hurst values by a multiscale analysis in an interval from 0.3 to 0.45 for all the used logs.
The accuracy of the results may be assessed by using standard deviation values calculated for Hurst curves according to each lithological unit.For AU, this parameter ranges between 0.073 and 0.129 for LWA, 0.020 and 0.109 for ALWA and 0.008 and 0.025 for PA.While for GU, it ranges between 0.088 and 0.158 for LWA, 0.037 and 0.135 for ALWA, and 0.004 and 0.048 for PA.Regarding VU, the value of this parameter varies from 0.067 to 0.185 for LWA, 0.019 to 0.146 for ALWA, and from 0.006 to 0.028 for PA.We can say that the PA results are the closest ones to those obtained by Fedi et al. (2005) and exhibit the most accuracy.
For comparison, the global spectral exponents obtained by different researches from KTB logs are gathered in Table 5.The spectral exponents related to our study are calculated from average Hurst values which are obtained by the three methods (LWA, ALWA and PA) for each sonic log.They are slightly different from other results, except those estimated by Leonardi and Kümpel (1998) using the R/S analysis.These authors found that the H values estimated by the R/S technique are consistently much greater than those determined by the power spectrum technique.Feder (1988) demonstrated that R/S yields biased results if H is not close to 0.7.In addition, Li (2003) assigned the discrepancy between the results of the two techniques to two main causes.First, well logging data are generally fBm processes, thus applying R/S directly on raw well-log data instead of their incremental series always leads erroneously to a high estimated H value (H >0.85).However, the analysis of incremental series by R/S gives the H estimate close to that obtained by power spectrum technique.The second cause is the sensitivity of R/S to the modifications affecting well log data.Indeed, the incremental series gives the best estimate of H value because these have the shortest transient zone and the highest correlation coefficient from least-squares regression.

Conclusions
The logs multifractionality may be characterized by the local Hurst functions estimated using LWA, ALWA and PA.The application of these techniques on synthetic models showed that PA led to the best results.The estimated Hurst values are the closest ones to the theoretical values with the least error, and the transition between two adjacent sequences is marked by a H-value jump discontinuity.Furthermore, the application is extended to real sonic logs recorded at KTB wells (VB-Vp, VB-Vs, HB-Vp and HB-Vs).The regularity profiles obtained from the logs allowed to identify the lithological discontinuities (the faults and the bounds of the geological sequences crossed by the wells), even if the latter do not appear on the sonic logs.Indeed, the lithological changes correspond to the local extrema of the Hurst function.Some peaks and lows of the Hurst values located within lithological units are related to the fluid-filled fractures.Overall, the Hurst function derived from sonic logs is strongly correlable with the lithology of the drilled formations.It may then constitute an attribute for characterizing the subsoil heterogeneity.An extended study to other borehole measurements is in progress in order to establish the correspondence between the Hurst value and the geological features of stratigraphic sequences.

Figure 1 .Fig. 2 .
Figure 1.Histograms of the Hurst values obtained by the three methods: (a) LWA (b 2 (c) PA, from 1000 realizations of a synthetic nhBm s(z) corresponding to a four-laye 3 geological model at the position (z=121.77m),located in the second layer (theoretic 4 =0.4). 5

Figure 4 .
Figure 4. Results obtained from sonic measurements recorded at the main borehole (HB). 1

Fig. 6 .
Fig. 6.A comparison between the Hurst functions obtained from sonic logs recorded at the main borehole (HB) and the estimated porosity log.The same abbreviations as in Fig. 5.

Table 1 .
Simulated four-layer geological model parameters.

Table 2 .
Statistical parameters of the obtained Hurst values corresponding to Fig. 1.

Table 3 .
Simulated multi-layer geological model parameters.PA, LWA and ALWA.Definitely PA presents the smallest values of the standard deviations and the bias.These results show that PA is the best algorithm yielding an accurate estimation of Hurst values.

Table 5 .
Spectral exponents for the KTB physical logs obtained from previous works.