A Grey System Approach for Estimating the Hölderian Regularity with an Application to Algerian Well Log Data

: The Hölderian regularity is an important mathematical feature of a signal, connected with the physical nature of the measured parameter. Many algorithms have been proposed in literature for estimating the local Hölder exponent value, but all of them lead to biased estimates. This paper attempts to apply the grey system theory (GST) on the raw signal for improving the accuracy of Hölderian regularity estimation. First, synthetic logs data are generated by the successive random additions (SRA) method with different types of Hölder functions. The application on these simulated signals shows that the Hölder functions estimated by the GST are more precise than those derived from the raw data. Additionally, noisy signals are considered for the same experiment, and more accurate regularity is obtained using signals processed using GST. Second, the proposed technique is implemented on well log data measured at an Algerian exploration borehole. It is demonstrated that the regularity determined from the well logs analyzed by the GST is more reliable than that inferred from the raw data. In addition, the obtained Hölder functions almost reﬂect the lithological discontinuities encountered by the well. To conclude, the GST is a powerful tool for enhancing the estimation of the Hölderian regularity of signals.


Introduction
The Hölderian regularity is a robust property of physical signals. The understanding of its time/space variations allows to extract meaningful information related to the physical nature of the analyzed signal.
However, the precise estimation of local regularity function of the mBm process is still difficult. This study aims to propose a technique based on grey system theory (GST) to enhance the accuracy of quantifying the local regularity of mBm paths. The suggested technique will be first tested on synthetic mBm paths, and next it will be applied on P-and S-wave velocity logs data measured at the KTB main borehole.
The GST is a concept first introduced by Deng [14]. Since then, it has gained an increasing interest in many research fields thanks to its suitability for analyzing systems whose parameters are partially known. In contrast with the conventional statistical techniques, the grey model needs limited data to predict the behavior of unknown systems. The GST is widely used to establish a forecasting model, to minimize randomness in time series and to increase the regularity of data [15][16][17]. Recently, the grey system theory has been combined with fractal estimation methods for improving results in the prediction or classification issues. In this context, Wu et al. [18] proposed an improved fractal model to forecast mine slope deformation by using the GST. They showed that the proposed method can make a more accurate prediction than the fractal model, additionally with high computational efficiency. Huang et al. [19] proposed a new method, grey detrended fluctuation analysis (GDFA), for calculating the fractal scaling exponent by combining grey system theory with grey detrended fluctuation analysis (DFA). They show that the fractal scaling exponent calculated using the proposed method has good antinoise ability. Chen et al. [20] proposed a different approach based on dual improved generalized fractal box-counting for improving traditional fractal box-counting dimension algorithm in subtle feature extraction of radiation source signals. Zhang et al. [21] considered the problem of signal feature extraction and classifier design under low SNR environment by proposing a method that first evaluates the multifractal dimension of nine modulated communication signals, and then an improved grey relation algorithm is used to recognize the extracted subtle characteristics. However, there are not works in literature combining GST and multifractal models.
In geophysics, the GST has been already applied on seismic datasets in order to minimize data randomness and increase data regularity, and thus to compute a meaningful fractal scaling exponent-based attribute that provides valuable information related to the distribution of sedimentary facies [19]. In this research, the GST technique is implemented for the first time on well log data (P-wave seismic velocity, bulk density, gamma ray, and photoelectric absorption factor). Well log data are complex to interpret since they are normally affected by noise (see, for example, [22]).
The major drawback of this analysis is that the suggested algorithms fail to accurately estimate the regularity functions. In this view, the GST is introduced to achieve best Hölder exponent estimates.

Local Hölderian Regularity
The Hölder exponent is defined for a stochastic process X, whose trajectories are continuous but nowhere differentiable, by Peltier and Lévy-Véhel [1]: This coefficient describes the local regularity strength of a signal. From a geometrical point of view, that means that the increments X(z) − X(z 0 ) in the neighborhood of z 0 are included by a Hölderian envelope defined by |X(z) − X(z 0 )| α X (z 0 ) . The higher α X (z 0 ) value corresponds to a smoother signal at z 0 , and conversely.

Multifractional Brownian Motion
The mBm is defined as [1]: In contrast with fBm, the pointwise Hölder exponent of W H(t) (t), α W = {α W (t), t ∈ R} is expressed as a function of the location. For each t, this exponent equals with probability one to H(t) [1,[30][31][32]:

Local Estimation of the Hölderian Regularity
Peltier and Lévy-Véhel [33] suggested an algorithm to estimate the local Hölder function H(z) at the point z = i/(n − 1), given by with S k,n (i), the local growth of the increment process, defined by 002C where n is the total number of signal samples, k is a fixed window size, and m is the largest integer not exceeding n/k.

Grey System Theory (GST)
In GST, GM (n,m) designates a grey model where n and m correspond to the order of the difference equation and the number of variables, respectively. Among the existing models, the most popular is GM (1,m) corresponding to first-order differential equation model of the m-type variable. Therefore, the GM (1,1) indicates the first-order differential equation model of a variable that is widely used in predictions because of its computational efficiency.
The calculation procedure of the GM (1,1) model are as follows [14][15][16][17]: (1) Data standardization The first required step of calculation is to standardize the signal values since all the amplitude data should be positive.
where min x (0) and max x (0) are the minimum and maximum of sequence x (0) , respectively.
(2) Modeling of the GM (1,1) is as follows: Selecting a subsequence denoted by Constructing an accumulation generation for the subsequence (2) , . . . , X where X (0) and X . . , N Then, the GM (1,1) is built using Sequence (5): where a and u are fitting parameters. The column a = [a, u] T is determined by the least squares method: (c) Replacing a in Equation (6), the grey model can be expressed by where the function X (1) is fitted to the signal accumulated series X (1) .
or simply computing X (0) (t) by the subtraction (e) Calculating the deviation (grey model error) between the accumulated sequence and the fitting function

Application to Simulated Data
First, the grey system technique is implemented on simulated sonic logs data. These logs are considered as multifractional Brownian motions (mBms) simulated using the procedure proposed by Peltier and Lévy-Véhel [1]. The idea consists of generating N fractional Brownian motions (fBms) B H (i) with the Hurst parameters H(i), i = 1, . . . , N, and then constructing the mBm path W H (i) by setting In our applications, the fBms are simulated using the successive random additions (SRA) algorithm [34].
Five types of Hölder functions, linear H 1 , periodic H 2 , logistic H 3 , synthetic H 4 , and stairs H 5 , are used to create mBm paths. They are defined as follows: 1+e −20(t−0.6) ), In the following, for each type of regularity, a local Hölder function is computed, using the algorithm detailed in Section 2.3, from the raw mBm path and the mBm path processed by the GST theory. The obtained results are presented in Figure 1. As it can be noted, the mean square error (MSE) corresponding to the raw standardized mBm is larger than that related to the mBm path processed by GST.
Then, the convergence of the suggested algorithm is evaluated on the basis of 1000 realizations of the mBm paths created using the Hölder function H 5 . Histograms of Hölder exponent values derived from the standardized raw mBm path (first line) and the mBm path processed by GST (second line) are illustrated in Figure 2. The first, second, third, and fourth columns are related to fixed positions whose indexes are i = N/8, 3N/8, 5N/8, and 7N/8 (N = 2048) which correspond to the theoretical Hölder exponent values: 0.2, 0.4, 0.6, and 0.8, respectively. It can be seen that the resulted histograms approximately follow a Gaussian distribution with small values of standard deviations. Note also that, for a given position, the mean of the Hölder exponent values derived from the mBm paths processed with the GST are closer to the theoretical H values than those which resulted in the standardized raw mBm paths. That confirms, again, the accuracy of the suggested method for the estimation of the local regularity functions. In the following, for each type of regularity, a local Hölder function is computed, using the algorithm detailed in Section 2.3, from the raw mBm path and the mBm path processed by the GST theory. The obtained results are presented in Figure 1. As it can be noted, the mean square error (MSE) corresponding to the raw standardized mBm is larger than that related to the mBm path processed by GST. Then, the convergence of the suggested algorithm is evaluated on the basis of 1000 realizations of the mBm paths created using the Hölder function H5. Histograms of Hölder exponent values derived from the standardized raw mBm path (first line) and the mBm path processed by GST (second line) are illustrated in Figure 2. The first, second, third, and fourth columns are related to fixed positions whose indexes are i = N/8, 3N/8, 5N/8, and 7N/8 (N = 2048) which correspond to the theoretical Hölder exponent values: 0.2, 0.4, 0.6, and 0.8, respectively. It can be seen that the resulted histograms approximately follow a Gaussian distribution with small values of standard deviations. Note also that, for a given position, the mean of the Hölder exponent values derived from the mBm paths processed with the GST are closer to the theoretical H values than those which resulted in the standardized raw mBm paths. That confirms, again, the accuracy of the suggested method for the estimation of the local regularity functions. The sensitivity of the proposed technique to noise level is demonstrated on simulated mBm paths. Here, the theoretical Hölder function considered in simulating mBm paths is a staircase (piecewise) function with nine pieces (stairs) presenting H values from 0.1 to 0.9 with an increment of 0.1. In each piece (stair), a fixed position is chosen on the simulated mBm path.
One thousand realizations of random noise are performed. For each realization, amounts of random noise are added to the noise-free mBm path to give signal-to-noise ratios (SNRs) of different values. In the following, SNR is defined as the ratio of signal power to the noise power (expressed in decibels).
For a given realization, the noisy mBm path is processed using GST. Then, the local H(t) is computed for both noisy and GST-processed mBm paths.
For a fixed SNR value, the H mean values of the noisy and GST-processed mBm paths at the nine fixed positions on mBm paths are obtained by averaging, over the total The sensitivity of the proposed technique to noise level is demonstrated on simulated mBm paths. Here, the theoretical Hölder function considered in simulating mBm paths is a staircase (piecewise) function with nine pieces (stairs) presenting H values from 0.1 to 0.9 with an increment of 0.1. In each piece (stair), a fixed position is chosen on the simulated mBm path.
One thousand realizations of random noise are performed. For each realization, amounts of random noise are added to the noise-free mBm path to give signal-to-noise ratios (SNRs) of different values. In the following, SNR is defined as the ratio of signal power to the noise power (expressed in decibels).
For a given realization, the noisy mBm path is processed using GST. Then, the local H(t) is computed for both noisy and GST-processed mBm paths.
For a fixed SNR value, the H mean values of the noisy and GST-processed mBm paths at the nine fixed positions on mBm paths are obtained by averaging, over the total realizations number, the H values of regularity functions obtained from these mBm paths, respectively, at these positions.
From Figure 3, it can be noted that, for all of the fixed positions associated with different theoretical H values (varying from 0.1 to 0.9), the GST-based method enhances the accuracy of the estimated regularity of the noisy mBm paths for all the considered noise levels. As expected, for all of the considered positions, the Hmean values of the GST-processed mBm paths are increasingly closer to the corresponding theoretical H value as SNR reaches high values. paths, respectively, at these positions.
From Figure 3, it can be noted that, for all of the fixed positions associated with different theoretical H values (varying from 0.1 to 0.9), the GST-based method enhances the accuracy of the estimated regularity of the noisy mBm paths for all the considered noise levels. As expected, for all of the considered positions, the Hmean values of the GST-processed mBm paths are increasingly closer to the corresponding theoretical H value as SNR reaches high values. For illustration, an application of the proposed method to a realization of an mBm path, simulated with a periodic function, is presented in Figure 4. As can be noted, the function H(t) estimated using GST is more accurate than that obtained from the noisy mBm path, and presents the least error estimation value. For illustration, an application of the proposed method to a realization of an mBm path, simulated with a periodic function, is presented in Figure 4. As can be noted, the function H(t) estimated using GST is more accurate than that obtained from the noisy mBm path, and presents the least error estimation value.

Application to Algerian Well Log Data
In this section, the aforementioned technique is implemented on well log data recorded in an Algerian exploration well located in southwestern Algeria. The analyzed logs are P-wave seismic velocity (Vp, m/s), bulk density (rhob, g/cm 3 ), gamma ray (GR, API), and photoelectric absorption factor (PEF). The sampling depth of all the measurements is

Application to Algerian Well Log Data
In this section, the aforementioned technique is implemented on well log data recorded in an Algerian exploration well located in southwestern Algeria. The analyzed logs are P-wave seismic velocity (Vp, m/s), bulk density (rhob, g/cm 3 Figure 5 presents the raw and the GST-processed log data. The local Hölder exponent logs from well logs are calculated using the algorithms detailed in Section 2, from the raw data (in red) and data processed by the GST (in blue) ( Figure 6). As can be observed, the lithological discontinuities, which are either layer boundaries or thin rock beds occurring within the studied depth intervals, are marked on the regularity profiles by sudden variations of H value.
Note, also, that the regularity logs computed using GST are more precise than those obtained from the raw logs, and lead to a more accurate lithological segmentation.
In addition, the correlation coefficients between the computed regularity logs are increased after application of GST on the investigated well logs (Table 1).    The H logs are computed using the algorithm detailed in Section 2.3, from the raw data (in red) and data processed by the GST (in blue). Note that the local extrema of the H logs correspond to the lithological discontinuities.
It can be noted that on the regularity profiles derived from the "porosity-dependent" logs, almost all the lithological discontinuities, which are either layer boundaries or thin rock beds occurring within the studied depth intervals, are marked by jumps in H value.  The H logs are computed using the algorithm detailed in Section 2.3, from the raw data (in red) and data processed by the GST (in blue). Note that the local extrema of the H logs correspond to the lithological discontinuities.
It can be noted that on the regularity profiles derived from the "porosity-dependent" logs, almost all the lithological discontinuities, which are either layer boundaries or thin rock beds occurring within the studied depth intervals, are marked by jumps in H value.
It is worthy to underline that for Vs log, the H logs obtained by both regularity estimation techniques (without and with GST) are matching across the studied depth intervals except for 700-900 m, 1650-1750 m, 2100-2200 m, 2650-2750 m, and 3050-3200 m. However, for Vp log, the regularity log obtained by the velocity log processed by GST presents relatively similar H values to those derived from the raw velocity data across almost all the investigated depth range, except within the depth interval of 2100-2200 m. The obtained results show that the GST allows to highlight local heterogeneities that cannot be seen from the regularity logs estimated from the raw velocity logs. Therefore, the regularity estimation obtained from the velocity well logs processed by the GST is more reliable than that yielded from the raw log data.

Conclusions
This study presents a new approach for accurately estimating Hölderian regularity using grey system theory (GST). The application on synthetic logs with different Hölder functions, simulated by the SRA method, shows that the regularity estimated from logs processed by the GST is more precise than that derived from the raw logs, even in the presence of noise. In addition, the suggested technique was implemented on Algerian exploration well log data. The estimated Hölder exponent logs allowed us to perform a lithological segmentation.