• The disturbances are one of the main issues encountered in Sidi Chennane phosphate deposit.

• Geophysical mapping is used to detect the disturbed zones hidden under the quaternary cover.

• The multifractal analysis is used to estimate the disturbances rate through the geophysical maps.

Morocco is the world’s leading exporter of sedimentary phosphate with more than 30 percent of world exports and holds more than 50 billion MT of phosphate in reserve, which represents more than 70% of world supply. Phosphate rock has been known in Morocco since 1908, but its economic potential was not discovered until 1917 (Ayad and Bakkali, 2017; Brives, 1908). These phosphate reserves are sufficient to meet several centuries of demand. In this paper, the work was carried out in the Sidi Chennane deposit consisting of a succession of several high-quality phosphate layers. This deposit is part of the Ouled Abdoun basin and it is restricted to the west by the RP11 regional road, to the east by the RP309 national road, to the south by the 226000 meridian, and to the north by the outcrops of the lower phosphate layer (Fig. 1).

Figure 1. Geographical map of Sidi Chennane phosphate deposit located in the sedimentary basin of Ouled Abdoun. The study area is delimited by a red rectangle.

The phosphate series of Sidi Chennane frequently suffers from the inclusion of many sterile bodies called disturbances or derangements consisting of a waste rocky material and characterized with irregular geometrical shapes (Fig. 2)(Ayad et al., 2019; Ayad and Bakkali, 2019). Geological studies carried out in Sidi Chennane have revealed that the origin of these sterile bodies is linked to the dissolution of Senonian gypsiferous layers located under the phosphate series (Ait Taleb et al., 2009).

Figure 2. Section of phosphate sequences fronts demonstrating disturbances.

The presence of disturbances in the phosphate mines causes two major problems: (1) from an economic point of view, it is not easy to determine exactly the part of these bodies in the overall volume of the deposit and one cannot therefore correctly calculate phosphate reserves; (2) from a technical point of view, when encountering a disturbance, it is necessary to increase the number of holes in order to draw up the boundaries of the disturbances, and to increase the explosive power during the blasting operation, thereby explosives and time-consuming increase drastically.

Geophysical methods constitutes an appropriate tool in the Sidi Chennane deposit to cover the areas that may be disturbed (Bakkali, 2005). The anomalies detected can be used as a potential characteristic to create a model of the disturbances and would allow defining these structures before the mining front reaches them. In this work, we aim to numerically describe the visual aspect of the geoelectric images to diagnose them and quantify disturbed areas so that a control strategy can be developed for further estimation of phosphate reserve.

During these decades, many pattern recognition tools have been widely used to quantify complex structures. One of them is fractal and multifractal analysis used to describe the irregular geometry of many natural phenomena (Mandelbrot, 1989; Mandelbrot, 1990). In this paper, we intend to numerically characterize the shape of the disturbances using the following multifractal parameters: The generalized fractal dimension D(q), the singularities of strength α(q), the local dimension f(α), and their conjugate parameter the mass exponent τ(q) as well as f(α)-α spectrum. The analysis was performed on 2D geoelectric images of an area of 50 hectares located in the Sidi Chennane phosphate deposit.

2.1. Fractal and Multifractal

Fractal theory was invented by French mathematician Benoit Mandelbrot in 1975 (Mandelbrot, 1967; Mandelbrot, 1975). Mandelbrot introduced the fractal as an extremely irregular shape for which any suitably chosen part is similar to the whole at increasingly smaller scales, for which the Hausdorff-Besicovitch dimension exceeds its topological dimension as described in Mandelbrot et al. (1983).

The fractal approach has been widely used as an effective spatial analysis in different fields to better understand the fundamental scaling of many structures (Thomas et al., 2012; Tunas et al., 2016; Abidin et al., 2017, Ayad, 2022). Fractal analysis uses various techniques to assess and characterize the complexity of many objects using different fractal parameters such as fractal dimension (Ayad and Bakkali, 2022). However, the patterns of many complex objects cannot be fully described using simple fractal (monofractal) analysis. This analysis suffers from a lack of self-comparability and cannot be used to reveal much detail.

The weakness of fractal analysis can be corrected using multifractal analysis. It is a systematic approach that can help discover the characteristics of all kinds of objects and structures using more than one scaling process (Stanley et al., 1988). It provides more information about the space-filling characteristics than a fractal analysis. Compared to monofractal, a multifractal structure can be characterized as a superposition of homogeneous monofractal structures and cannot be described using the fractal dimension. Multifractal analysis has been widely employed successfully to characterize the complexity of many natural phenomena in different fields of scientific research (Lindinger et al., 2017; Mandal et al., 2017; Wang et al., 2017).

2.2. Geoelectrical Survey

Geoelectric prospecting is a tool for measuring the electrical resistivity of rocks below the surface. The working principle of this geophysical prospecting method is to inject an electric current of intensity I between two electrodes, and the voltage difference value (∆V) is measured between two potential electrodes after the current through the rocks. Then, the apparent resistivity of the rocks is found using the following formula: ρ_a = K (∆V/I), where K is the geometric factor of the instrument which depends on the configuration of the current and potential electrodes. The geoelectrical data used in this study was collected over an area of 50 hectares in the Sidi Chennane phosphate deposit. The measurements were performed by means of a Syscal2 resistivity meter using a 20m × 5m rectangular array. A total of 51 geoelectrical profiles spaced 20m apart were executed to reach a depth of 40m. There were 101 stations at 5m distance for each profile, which makes 5151 measurement stations.

The acquired data was then interpolated and imaged into a 2D contour map (Fig. 3). The map obtained allows defining the anomalies corresponding to the disturbances. It is used as a marker to distinguish between disturbances and phosphate rocks. The anomalies were classified as: (1) normal phosphate rocks of resistivity less than 200Ω.m; (2) disturbances of resistivity that mainly exceeded 200Ω.m.

Figure 3. Geoelectrical contour map of the study area.

2.3. Multifractal Analysis

To analyze the spatial pattern of disturbances, two sets of parameters are needed, including global and local parameters. The global parameters include the generalized correlation dimension D(q) and the mass exponent τ(q), while the local parameters include the singularity exponent (Lipschitz–Hölder exponent) α(q), and the local dimension (or Hausdorff dimension) f(α) of the set supporting this exponent (Chhabra and Jensen, 1989; Chhabra et al., 1989). The generalized correlation dimension D(q) is the basic multifractal parameter, based on the Renyi entropy. D(q) basically deals with the variation of mass varies as a function of ε (resolution or box size) in a digitized map. To calculate it, the spatial pattern of a fractal structure will be amplified step by step to some arbitrary exponent, q (-∞Hentschel and Procaccia (1983):

(1)$$D\left(q\right)=\underset{\epsilon \to 0}{\lim }\frac{I_q\left(\epsilon \right)}{1n\left(1/\epsilon \right)}$$

Where q is a continuous variable called the order of moment and Iq(ε) defines the Renyi entropy.

(2)$$I_q\left(\epsilon \right)=1n\sum _{i=1}^{N\left(\epsilon \right)}Pi{\left(\epsilon \right)}^q$$

Where Pi(ε) is the distribution, indicating the ratio of the number of elements appearing in the i^th box (rectangular grid) of linear scale ε to the number of the set, N(ε) represents the number of nonempty boxes.

The most frequently used generalized dimensions are D0 for q= 0; D₁ for q= 1 and D₂ for q= 2, which are respectively called the capacity dimension (or box counting dimension), information dimension (or entropy dimension) and correlation dimension.

The generalized dimension, D(q), is related to the other sets of multifractal exponents. Therefore, D(q), is obtained from the relationship with the second global parameter. τ(q) is known as the mass exponent function of q, which generalizes the relation given by:

(3)$$\tau \left(q\right)=D_q\left(q-1\right)$$

The Hausdorff fractal dimension f(α) is a scale exponent of boxes with a common α, called Holder or singularity exponent. A graph of singularity spectrum can be plotted using f(α) versus α. The functions α and f(α) can be derived by the Legendre transformation as (Halsey et al., 1986):

(4)$$\alpha \left(q\right)=\frac{d\tau \left(q\right)}{dq}=D_q+\left(q-1\right)\frac{d{D}_q}{dq}$$

(5)$$f\left(\alpha \right)=q\alpha \left(q\right)-\tau \left(q\right)=q\alpha \left(q\right)-\left(q-1\right)D_q$$

The Legendre transform is mainly used to relate local parameters with global parameters. This suggests that the two sets of parameters can be converted into each other. If we calculate global parameters, we will obtain local parameters and vice versa.

In practice, from the τ(q) or D(q) curves, one can decide whether multifractality exists or not. If the curve τ(q) is a straight line or D(q) is a constant, then the structure is monofractal. However, if τ(q) or D(q) have convex shapes, then the structure is multifractal.

The multifractal process was performed using open source software called FracLac, a plugin for ImageJ (Karperien, 2004). The program uses the box-counting method to gather information about the distribution of pixel values in 2D binary images. Note that the box counting method consists to cover the image with a grid of decreasing box size, and then the boxes that include part of the figure are counted (Russel and Hanson, 1980). The collected data thus became the basis for generating graphs describing the multifractal scaling in order to show how the disturbances pattern behaves from one image to another.

In order to create a comparative analysis framework to obtain an accurate classification of the disturbed areas by the multifractal approach, we left sedimentary phosphates rocks out of our interpretation and we retain only the disturbed areas. Thus, the edges of the disturbances were traced and delineated by extracting their boundaries at different apparent resistivity cut-off frequencies (ʋ_ar). Eight geoelectric images from ʋ_ar=200Ω.m to ʋ_ar=340Ω.m were created from the original image using a ʋ_a.r cut-off frequency of 20Ω.m. This allows observing how much the surface of the disturbances changes from one image to another. These maps were then converted to white and black binary images using Surfer 12 software. Where black pixels correspond to disturbed areas while white pixels represent sedimentary phosphates rocks (Fig. 4).

Figure 4. Spatial patterns of the disturbances extracted at eight cutoff frequencies through the geoelectrical contour map and their corresponding black and white images.

Before starting the analysis process, the main points of the multifractal approach must be explained. Multifractal spectrums have the potential to find differences in spatial information. It can provide us with multiple viewing angles from which to examine the surface more carefully. The global and local parameters were calculated for the order of moment -10 < q < 10. It can be noted that q can be used as a spatial analysis adjuster, by changing the value of q; we can analyze coarse structures and fine elements at different levels of view. The differences and similarities can be reflected by the multifractal spectra of the local and global parameters D(q), τ(q), α(q), and f(α), as shown in Fig. 5, 6, 7, and 8.

Figure 5. The generalized correlation dimension D(q) spectrums of the different images.

Figure 6. The mass exponent τ(q) spectrums of the different images.

Figure 7. The singularity exponent α(q) spectrums of the different images.

By using the global parameters (D(q) and τ(q)) (Fig. 5 and 6), we can analyze the whole spatial pattern of the disturbances. The D(q) curves are decreasing monotonic functions: D(q) ˃ D(q + 1). Roughly speaking, the functions D(q) are decreasing for negative values of q, while for positive values of q the curves are increasing. If q > 1, the high density pattern will be amplified, and therefore, the value of D(q) will decrease. Otherwise, if q < 1, the low density pattern will be amplified, and thus the value of D(q) will increase. In general, small differences in the disturbed surface may lead to significant changes in the resulting curves.

Moreover, for more precision, the simplest method for judging the results consists of comparing the capacity dimension, the information dimension, and the correlation dimension. Thus, since the spectrum D(q) is a monotonically decreasing function which decreases with increasing of the order of moment q, the capacity dimension D₀ is necessarily greater than the information dimension D₁ which is greater than the correlation dimension D₂. Table 1 shows the results of D₀, D₁ and D₂ of all analyzed images. It is obvious that, for all images, the inequality D₀ > D₁ > D₂ is always established and respects the monotonically decreasing principle of D(q). The range of values (D₀, D₁ and D₂) of the different analyzed images is between 0 and 2, i.e. 0

Table 1 . The calculated results of the capacity dimension, correlation dimension, and information dimension as a function of the corresponding rate of disturbances.

Images	Disturbances rate %	Capacity dimension D0	Information dimension D1	Correlation dimension D2
υ200Ω.m	17.70	1.558	1.457	1.410
υ220Ω.m	12.30	1.336	1.215	1.169
υ240Ω.m	8.96	1.269	1.143	1.094
υ260Ω.m	7.43	1.269	1.115	1.050
υ280Ω.m	6.40	1.341	1.189	1.129
υ300Ω.m	5.60	1.214	1.058	0.988
υ320Ω.m	4.81	1.195	1.038	0.964
υ340Ω.m	4.10	1.161	0.995	0.916

The capacity dimension D₀ is 1.558 for υ_200Ω.m, and it decreases to 1.161 for υ_340Ω.m, the information dimension D₁ is 1.457 for υ_200Ω.m, and it decreases to 0.995 for υ_340Ω.m, and the correlation dimension D₂ is 1.410 for υ_200Ω.m, and it decreases to 0.916 for υ_340Ω.m.

For the second global parameter, the mass exponent τ(q), show nonlinear curves. It is a monotonically increasing function τ(q)< τ(q + 1). Both left tails and the right ones of the curves show different slopes, indicating significant multifractal variability. The difference between these curves reveals the dissimilarity of the disturbed areas within the different images. In contrast, τ(q) tends to be linear for such pattern with lowest variation. In addition, when the statistical moment q = 1 the mass exponent is zero. In this case, there is a significant overlap of the curves. Moreover, it was noted that although the changes were not robust, these spectra could reveal significant differences.

By using the local parameters (f(α) and α(q)), we can analyze the small details of the disturbances (Fig. 7 and 8). This is different from the global parameters (D(q), τ(q)), which analyze the whole pattern. The singularity exponent curves α(q) show a decreasing monotonous behaviors α(q)˃ α(q + 1), while the spectra f(α) have the maximum f[α(q)]max=D₀ when q = 0. Indeed, if q≤ 0, f(α) is a monotonically increasing quantity, and if q≥0, f(α) is a monotonically decreasing quantity. For q≤ 0, we have f [α(q)]≤f [[α(q + 1)] , while for q≥ 0, f [α(q)] ≥f [[α(q+ 1)]. The spectral behaviors of these parameters indicate a varying degree of multifractality for the analyzed images. It seems appropriate for estimating the degree of disturbances, which can allow a better classification of phosphate reserves (undisturbed).

Figure 8. The local dimension f(α) spectrums of the different images.

If the relationship between the local parameters f(α) and α(q) is plotted, an additional unimodal multifractal spectrum f(α)-α can be created (Fig. 9).

Figure 9. The singularity spectrums f(α)-α of the different images.

The f(α)-α spectra show an asymmetric shape, where the left tails are mostly shorter than the right ones. In all cases, these curves are found to be continuous with convex parabolas. When we comparing the spectrums f(α), we notice that their shapes are different, since all of the analyzed images are subject to different disturbances patterns. The dissimilarity of these curves indicates the possibility of associating these multifractal parameters to the rate of disturbances within the phosphate deposits.

Since the f(α)-α spectra does not provide a unique value, but a set of values, it is important to find the index which could be a discriminator in our application. For this purpose, the next step is to extract the indexes that describe these spectra. Several useful indices have been estimated in this work: α_min and α_max are respectively the minimum and maximum values of α, and f(α_min) f(α_max) (Table 2).

Table 2 . Results of the parameters extracted from the singularity spectrum as a function of the corresponding rate of disturbances.

Images	Disturbances rate %	αmin	αmax	f(αmin)	f(αmax)
υ200Ω.m	17.70	1.311	2.351	1.073	0.587
υ220Ω.m	12.30	1.087	2.126	0.840	0.540
υ240Ω.m	8.96	1.006	2.057	0.755	0.484
υ260Ω.m	7.43	0.942	2.077	0.693	0.541
υ280Ω.m	6.40	1.031	2.177	0.787	0.460
υ300Ω.m	5.60	0.885	2.017	0.636	0.278
υ320Ω.m	4.81	0.868	2.004	0.661	0.257
υ340Ω.m	4.10	0.809	1.930	0.545	0.216

According to Table 1, values of the estimated indexes can symbolize a discrimination of the analyzed images. The image ʋ₂₀₀ which corresponds to the highest disturbance rate (17.683%) presents the highest values (α_min, α_min, f(α_min) and f(α_max)). However, the image ʋ₃₄₀ which corresponds to the lowest disturbance rate (4.104%) presents the lowest multifractal values (α_min, α_min, f(α_min) and f(α_max)). We conclude that the increase in the surface occupied by the disturbances implies an increase in the variability of α, and vice versa.

In this research, we systematically analyze the potential utility of the multifractal approach to characterize the disturbed areas of the phosphate series of Sidi Chennane. Two sets of multifractal parameters are used: One is the global parameters, contain the generalized correlation dimension D(q) and the mass exponent τ(q), and the second is the local parameters include the singularity exponent α(q), and the Hausdorff fractal dimension f(α). The work was carried out on the geoelectric images of an area located in the northern part of the Sidi Chennane phosphate deposit, Ouled Abdoun, Morocco.

Eight geoelectric images were analyzed to make a comparative analysis and validate the potentiality of the multifractal approach. Clear results are obtained and suggest the strong sensitivity to characterize the spatial patterns of the disturbances through the geoelectric images. Comparing the similarities and differences can help determine how the surface of the disturbed areas changes, which may differentiate the phosphate deposit at high-risk of disturbances and the deposit at low-risk. The multifractal analysis approach can be used as an important analytical tool to make the best exploration and exploitation planning and obtain an automatic and accurate classification of phosphate reserves.