# Fractal Assessment of the Disturbances of Phosphate Series Using Lacunarity and Succolarity Analysis on Geoelectrical Images (Sidi Chennane, Morocco)

Authors: Abderrahim Ayad; Saad Bakkali
Journal: Complexity
 (2019)

---

## Abstract

In Sidi Chennane phosphate deposit, the disturbances cause serious challenges for the OCP mining engineers during both exploration and the exploitation stages. They are qualified as worst rocky hard to be quantified since they interfere with phosphates series. In this paper, we propose the fractal analysis as a simple efficient tool to quantify the rate of the disturbances by two different methods, namely, the Lacunarity and the Succolarity. The analysis was carried out on eight geoelectrical images of a disturbed area of 50 hectares located in the northern part of Sidi Chennane. The results proved that there is a strong correlation between the disturbances rate and the corresponding fractal value indexes. It appears that the distinction between two disturbed areas is the difference between their corresponding fractal values. This has an important implication to discriminate between the phosphate deposit at high risk of disturbances and the deposit at low risk. The fractal analysis can be thus used as a crucial concept in ranking the prospective zones of phosphate as well as improving the phosphate reserve estimation.


---

## Body

## 1. Introduction

Morocco owns more than three-quarters of phosphate-rock reserves in the world (50 000 billion metric tons). These appreciable phosphate reserves known since 1908 were deposited in several basins, from northeast to southwest of Morocco and have been exploited by l’Office Chérifien des Phosphates (OCP). The large extraction intensity since the 1920s allows Morocco to become the world’s leading producer of phosphate with mining capacities in excess of 32 million tons per year [1, 2].The main phosphate mines exist in the Ouled Abdoun Basin located west of the Atlas Mountains in the central part of Morocco, at about 200 km in the south-east of the capital Rabat [3] (Figure 1). It is the largest basin encompassing an area of 4500 square km, containing 44% of the whole reserves and at least 26.8 billion tons of phosphate.Figure 1
Geographical location of Sidi Chennane mine in the sedimentary basin of Ouled Abdoun.In this study, we focus on Sidi Chennane mine, one of the most important phosphate deposits in the Ouled Abdoun Basin, located near the city of Khouribga. The phosphate series of Sidi Chennane is composed of several distinct phosphate layers alternating with layers of calcareous and argillaceous hardpan. In this deposit, the rocky material is extracted at open-pit and the phosphate mineral is then separated from associated gangue (Marls, limestones, silex, clays, etc.) by a combination of various processing stages involving crushing and screening, washing, and flotation.The tabular structure of the phosphate series of Sidi Chennane suffers frequently from the inclusion of many sterile bodies composed of limestone blocks within an argillaceous matrix considered as waste rocky material that must be removed (Figure2). These bodies are locally qualified as “derangements” or “disturbances” since they disturb the phosphate series regularity which causes many problems during both exploration and the exploitation stages. These bodies have a conical or subcircular form; their diameter varies from 10m to more than 150m. They are hidden under the quaternary cover (topsoil) and cannot be directly mapped from the surface. The geological studies carried out in Sidi Chennane have revealed that the origin of these sterile bodies is referred to natural collapse caused by dissolution of Gypsiferous rocks under the phosphate layers [4, 5]. In fact, the appearance of the sterile bodies evokes a partial dissolution, a weathering phenomenon, even a leaching of carbonates, clays, etc.Figure 2
Geological sections showing disturbed phosphate series fronts in Sidi Chennane quarry.In their interesting paper [6], Boujo notes that the origin of the disturbances is often supergene, even pedological, resulting in a favorable way such as deposits of residual concentration, for example, in Morocco (Sidi Daoui); in Brazil (Olinda), (Congaçari); and in Senegal (Taïba deposit).In general, there are two main types of sterile bodies, based on the type of material, the size, the hardness, or type of contact with the phosphate series. The first type is found throughout the mineral deposit: it appears to be a mixture of clays, marls, limestones, and phosphate with large amounts of cherty limestone. The second type is the most abundant, and it is highly disturbing and appears as an accumulation of phosphate limestone blocks with large nodules of marl, and some fragments of phosphate rock.Moreover, the disturbances extended over the phosphate deposit of Sidi Chennane with densities more or less important. The statistical analysis of the investigated boreholes performed at Sidi Chennane showed that, from 157 boreholes, 53 revealed a “disturbed” phosphate series, i.e., 33% of the investigation boreholes.The presence of the disturbances causes serious problems in the exploitation of open cast in Sidi Chennane: firstly, as they are compact and hard, they disturb seriously the kinematic chain of the exploitation in some yards. Indeed, during the firing stage, the boring-grids have to be tightened. In addition, operators must fill the necessary boreholes with dynamites; thereby explosives and time-consuming increase drastically. Secondly, as we do not know the disturbed volume proportion, the phosphate reserves estimations can be wrong. For these reasons, the quantification of the disturbances rate within the phosphate deposit is considered as a crucial step and a challenging question for the OCP mining engineers [7, 8].The normal phosphate-bearing rock and the disturbances differ mainly by their chemical and physical properties, and as a result, there are strong resistivity contrasts between them. For that reason, the geoelectrical prospecting methods were used as appropriate tools able to image the disturbances hidden under the quaternary cover and differentiate them from the phosphate rocks [9, 10]. The geoelectrical anomalies detected are used as a potential feature to obtain an image of the disturbances in the subsurface.In this work, we intend to numerically quantify the disorder of the disturbed areas displayed on the geoelectrical images. So how to design a useful tool which could effectively distinguish a variety of disturbed areas is the key problem that should be considered.Many analysis methods were used to quantify the texture of varied objects in different scientific research fields. One of them is the fractal geometry coined by Benoît Mandelbrot since 1974 [11]. It is introduced to solve the complex character of objects in which their shapes show an irregularity and roughness [12–15]. Here we present a widely applicable fractal approach to quantify disorder of the disturbances, based on two analysis methods, namely, the Lacunarity (λ) and Succolarity (σ) [16, 17].Lacunarity analysis has been developed to describe and analyze the spatial distribution of gaps in images texture, in the sense that images of different texture have different Lacunarity and viscera [18]. However, Mandelbrot defines the Succolarity (connectivity of the gaps) as the degree of percolation of a virtual fluid inside an image according to the direction of the penetration [19].These fractal analysis methods were applied extensively in many research studies like image processing, geophysics, medicine, and other fields [20–23]. Using this approach of analysis, we tried to point out that the fractal methods may be more appropriate to quantify the disturbed areas since they commonly show an irregular geometry displayed on the geoelectrical images.


## 2. Overview of the Fractal Geometry

The fractals were invented by a French-American mathematician Benoît Mandelbrot [24]. The term “fractal” came up from the Latin root “fractus” which means “broken” or “fractured” (e.g., irregularly crushed and broken stone). Mandelbrot began his treatise on fractal geometry by explaining that when a fractal object is divided into parts, each part would be a smaller replica of the whole object. This is the central concept of fractal geometry summarized in the notion “self-similarity.” This means that the object is statistically the same, no matter what scale it is viewed at. On part resembles the whole. A classic example is the measurement of the length of the perimeter of the coastline of Britain. If we ask “How long is the coast of Britain,” the answer to this question depends on how the measuring stick length is. For example, if the perimeter is measured with a ruler, as we let the length of the ruler get small, the length of the coast becomes infinite because the surface of the coast shows increasingly details the closer we zoom into it. Many natural objects such as mountains, trees, clouds, river deltas, galaxies, plants, vascular systems, and snowflakes are self-similar after a few orders of magnitude of zooming in. The geometry of these natural objects is called Fractal Geometry.There are fractals called mathematical fractals generated by exact iteration and are infinitely self-similar. This means we can see an exact copy of the fractal at each point of self-similarity. Much effort has been made centuries ago to establish objects that are now described as “strict” and “exact” fractals, for example, the Von Koch curve, the Cantor set, the Sierpiński curve, the Dragon curve, and the Peano curve [25–28].Over the years, fractal analysis becomes a valuable tool introduced complementary to solve the complex character of objects. The main attraction of this approach stems from its ability to offer an appropriate way to quantitatively analyze data in a variety of fields, such as economics, physics, earth sciences, biology, chemistry, and other disciplines [29, 30].


## 3. Materials and Methods

In this study, the fractal analysis was performed on geoelectrical images of an area of 50 hectares located in the northwest part of Sidi Chennane deposit [31, 32]. The geoelectrical anomalies allow to map and constrain the anomalous regions corresponding to disturbances. The successful analysis must be on a priori distinction between the apparent resistivity of the phosphate rocks and the disturbances. In terms of our example, the geoelectrical anomalies are classified into disturbances of apparent resistivity exceeding 200Ω.m against 200Ω.m for the phosphate rocks. As we know that we are working only on the disturbed surfaces, we can most likely leave sedimentary phosphates rocks out of our interpretation. We start by extracting the disturbances boundaries at different apparent resistivity cutoff frequencies (ʋar) using Surfer 12 software. Eight geoelectrical images were produced from 200Ω.m to 340Ω.m using a cutoff frequency of 20Ω.m, considering that using the combination of these images will create a comparative framework for analysis and serve to crosscheck and validate the results (Figure 3).Figure 3
The visual spatial patterns of the disturbances extracted from 200Ω.m to 340Ω.m using a cutoff frequency of 20Ω.m. The cutoff frequenciesʋa.r used are, respectively, (a) ʋa.r = 200Ω.m, (b) ʋa.r = 220Ω.m, (c) ʋa.r = 240Ω.m, (d) ʋa.r = 260Ω.m, (e) ʋa.r = 280Ω.m, (f) ʋa.r = 300Ω.m, (g) ʋa.r = 320Ω.m, and (h) ʋa.r = 340Ω.m.
(a) (b) (c) (d) (e) (f) (g) (h)Before starting the analysis process, the geoelectrical images have to be transformed digitally into black and white images in order to be adapted to the fractal analysis programs requirement. Consider that in the Lacunarity case, the black pixels represent the disturbances, while the white pixels represented the phosphate rocks (Figure4). While, in the Succolarity case, the white pixels represent the disturbances, while the black pixels represented the phosphate rocks (Figure 5).Figure 4
Black and white version of images in Figure3 adapted to the Lacunarity analysis requirement.
(a) (b) (c) (d) (e) (f) (g) (h)Figure 5
Black and white version of images in Figure3 adapted to the Succolarity analysis requirement.
(a) (b) (c) (d) (e) (f) (g) (h)Moreover, the precise analysis requires the choice of adequate resolution for the analyzed images. The wrong resolution in the input binary images may lead the method to yield results not corresponding to the results expected. In light of this issue, we carried out several tests using different resolutions, and we have obtained different fractal values for the same image. The chosen resolution corresponds to 510×255 pixels in the Lacunarity analysis case, and 577×289 pixels in the Succolarity analysis case and this is for all images analyzed.The images show a view of the geometric forms of disturbed areas boundaries. The edges extracted at different cutoff frequenciesʋa.r allow observing to what extent the surface of the disturbance changes from one image to another. Then, for each black and white image, the surface of the disturbances was explored by means of the Lacunarity and Succolarity approaches.


### 3.1. The Lacunarity Analysis

The Lacunarity (λ) was calculated using the most popular method called “gliding box” (or “moving windows”) developed by Allain and Cloitre [33]. The implementation of this method consists of placing a box of size r at the top left corner of a binary image (black and white pixels), and then the black pixels -p- in it were counted. Afterward, the box is glided to the right, and the number of black pixels in it is counted again. This operation is repeated for many times until the box reaches the bottom right corner of the image and this for all possible box size (see Figure 6). This process will be used for obtaining the box occupation frequency distribution B(p,r).Figure 6
The gliding box principal for calculating Lacunarity.Next, Lacunarity for a specific scale is estimated as the ratio of the first momentZ(1) that it is the average sum of values in all possible boxes of size r, and the second moments Z(2) that it is the average squared sum of values in all possible boxes of size r.(1)Z1r=∑ppBp,r(2)Z2r=∑pp2Bp,rUsing (1) and (2), the Lacunarity λ, for a given box size r can be calculated using the following equation:(3)λr=Z2rZ1r2In this study, the Lacunarity was computed using the FracLac 2.1 program, plugin for Image J [34]. FracLac outputs a text file containing the values of the Lacunarity and their corresponding boxes sizes. The data obtained will be represented in a graphic of ln⁡(λ) versus ln⁡(r) generating thus a curve for each image.


### 3.2. The Succolarity Analysis

The Succolarity(σ) was originally invented by Mandelbrot since 1977 (Mandelbrot 1983), but he did not define a method for calculating it. To our knowledge, the first attempt for measuring the Succolarity was that of de Melo and Conci [35, 36] using the box-counting approach proposed by Russel [37]. De Melo points out that the Succolarity property is used to indicate the capacity of a virtual liquid to flood a binary image in different orientations. The analyzed image holds thus two kinds for pixels: black pixels (penetrable) considered as voids on the image, which allow the liquid to pass and flood the area, and the white pixels (impenetrable) considered as obstacles to the liquid. In order to apply the algorithm, first, the virtual liquid will percolate the image in a given direction and then in the opposite direction ensuring that all black pixels were flooded: e.g., from the top side to the bottom side (t2b), and then, from the bottom side to the top side (b2t), or from the left side to the right side (l2r) and then, from the right side to the left side (r2l) (see Figure 7).Figure 7
The Succolarity analysis approach. (a) From top to the bottom (t2b). (b) From bottom to the top (b2t). (c) From left to right (l2r). (d) From right to the left (r2l).The flooded image will be then analyzed using the box-counting method. A grid of equal boxes size (BS(n), where n is the number of possible divisions) is placed on the flooded image. Then the occupation percentage (OP) is measured within each box size r. Afterward, the sum of the multiplications of OP(BS(r)) (where r ranges from 1 to n) by the pressure matrix PR(BS(r),pc) are calculated for each box size r (where pc is the position on x or y of the centroid of the box on the scale of pressure). Then, the Succolarity is calculated for each direction (dir) using the following equation:(4)σBSrdir=∑r=1nOPBSr×PRBSr,pc∑r=1nPRBSr,pcThe obtained Succolarity (σ) values will be graphically represented as a function of the dividing factor d on an ln×ln plot. The results can be then analyzed in a given direction or using all flooded directions.In this study, the Succolarity computation was performed using a program (fractals-CxImage) developed by Rafael H. C. de Melo from Federal Fluminense University in Niterói, Rio de Janeiro, Brazil.


## 3.1. The Lacunarity Analysis

The Lacunarity (λ) was calculated using the most popular method called “gliding box” (or “moving windows”) developed by Allain and Cloitre [33]. The implementation of this method consists of placing a box of size r at the top left corner of a binary image (black and white pixels), and then the black pixels -p- in it were counted. Afterward, the box is glided to the right, and the number of black pixels in it is counted again. This operation is repeated for many times until the box reaches the bottom right corner of the image and this for all possible box size (see Figure 6). This process will be used for obtaining the box occupation frequency distribution B(p,r).Figure 6
The gliding box principal for calculating Lacunarity.Next, Lacunarity for a specific scale is estimated as the ratio of the first momentZ(1) that it is the average sum of values in all possible boxes of size r, and the second moments Z(2) that it is the average squared sum of values in all possible boxes of size r.(1)Z1r=∑ppBp,r(2)Z2r=∑pp2Bp,rUsing (1) and (2), the Lacunarity λ, for a given box size r can be calculated using the following equation:(3)λr=Z2rZ1r2In this study, the Lacunarity was computed using the FracLac 2.1 program, plugin for Image J [34]. FracLac outputs a text file containing the values of the Lacunarity and their corresponding boxes sizes. The data obtained will be represented in a graphic of ln⁡(λ) versus ln⁡(r) generating thus a curve for each image.


## 3.2. The Succolarity Analysis

The Succolarity(σ) was originally invented by Mandelbrot since 1977 (Mandelbrot 1983), but he did not define a method for calculating it. To our knowledge, the first attempt for measuring the Succolarity was that of de Melo and Conci [35, 36] using the box-counting approach proposed by Russel [37]. De Melo points out that the Succolarity property is used to indicate the capacity of a virtual liquid to flood a binary image in different orientations. The analyzed image holds thus two kinds for pixels: black pixels (penetrable) considered as voids on the image, which allow the liquid to pass and flood the area, and the white pixels (impenetrable) considered as obstacles to the liquid. In order to apply the algorithm, first, the virtual liquid will percolate the image in a given direction and then in the opposite direction ensuring that all black pixels were flooded: e.g., from the top side to the bottom side (t2b), and then, from the bottom side to the top side (b2t), or from the left side to the right side (l2r) and then, from the right side to the left side (r2l) (see Figure 7).Figure 7
The Succolarity analysis approach. (a) From top to the bottom (t2b). (b) From bottom to the top (b2t). (c) From left to right (l2r). (d) From right to the left (r2l).The flooded image will be then analyzed using the box-counting method. A grid of equal boxes size (BS(n), where n is the number of possible divisions) is placed on the flooded image. Then the occupation percentage (OP) is measured within each box size r. Afterward, the sum of the multiplications of OP(BS(r)) (where r ranges from 1 to n) by the pressure matrix PR(BS(r),pc) are calculated for each box size r (where pc is the position on x or y of the centroid of the box on the scale of pressure). Then, the Succolarity is calculated for each direction (dir) using the following equation:(4)σBSrdir=∑r=1nOPBSr×PRBSr,pc∑r=1nPRBSr,pcThe obtained Succolarity (σ) values will be graphically represented as a function of the dividing factor d on an ln×ln plot. The results can be then analyzed in a given direction or using all flooded directions.In this study, the Succolarity computation was performed using a program (fractals-CxImage) developed by Rafael H. C. de Melo from Federal Fluminense University in Niterói, Rio de Janeiro, Brazil.


## 4. Results and Discussion

The validity of this study has been verified by applying Lacunarity and Succolarity fractal methods. The data computed were then plotted in bilogarithmic graphs ln-ln. In the Lacunarity case, graphs display theln⁡(λ) values on the vertical axis while the ln⁡(r) values are depicted on the horizontal axis. In all cases, Lacunarity was calculated for box sizes ranging from 5 to 289 and produce a curve which is concave downwards with a high correlation coefficient R2> 0.9. Figure 8 summarizes the results.Figure 8
Logarithmic scatter plots of Lacunarity values as a function of gliding box size r using as input the images in Figure4. (a) ʋ200Ω.m. (b) ʋ220Ω.m. (c) ʋ240Ω.m. (d) ʋ260Ω.m. (e) ʋ280Ω.m. (f) ʋ280Ω.m. (g) ʋ320Ω.m. (h) ʋ340Ω.m.
(a) (b) (c) (d) (e) (f) (g) (h)The shape of the illustrated spectra (Figure8) showing monotone decreasing behavior that can be adjusted as a declined function toward zero as the box size increased. The comparison of the generated curves indicates differences in an “integral Lacunarity,” i.e., Lacunarity overall range of scales. For small boxes size, the Lacunarity is high due to reduced occupied space. For larger boxes size, the Lacunarity is lower due to the largest occupied space (disturbed areas).As one can see, the curve corresponding to imageʋ340Ω.m has the largest spectrum, while that of image ʋ200Ω.m has the shortest spectrum. This can be explained by the fact that the image ʋ200Ω.m has the lowest degree of void spaces and corresponds thus to the highest disturbed surface. While image ʋ340Ω.m has the highest degree of void spaces and corresponds to the lowest disturbed surface. This remarkable difference can be confirmed if one analyzes the results and at the same time observing the images depicted in Figure 4.In the Succolarity case, the analysis was done along the four privileged directions indicated by de Melo: from the top side to bottom side (t2b); from the bottom side to top side (b2t), from the right side to left side (r2l); and from the left side to right side (l2r). Then, the Succolarity values obtained were plotted as a function of their corresponding dividing factor d, in bilogarithmic graphln⁡(σ)×ln⁡(d). The resulting curves were elaborated in detail in Figure 9.Figure 9
Logarithmic scatter plots of Succolarity values a function of the corresponding dividing factor d using as input the images in Figure5. (a) ʋ200Ω.m. (b) ʋ220Ω.m. (c) ʋ240Ω.m. (d) ʋ260Ω.m. (e) ʋ280Ω.m. (f) ʋ300Ω.m. (g) ʋ320Ω.m. (h) ʋ340Ω.m.
(a) (b) (c) (d) (e) (f) (g) (h)Characteristics that could be understood by these plots (Figure9) are that the values from the vertical analysis have similar behavior and that the values for the horizontal analysis are still the same. In other words, the values of Succolarity calculated both from the top-to-bottom and from the bottom-to-top directions were identical, as were the left-to-right and right-to-left directions almost match. For both the two directions, the resulting curves vary considerably from the image ʋ200Ω.m to image ʋ340Ω.m. The difference depends on the space occupied by the disturbances in each image. This conclusion can be easily explained and reinforced when it is linked to the visual inspection of the digitized images depicted in Figure 5.The Succolarity analysis has given essential information as a powerful method to texture classification purpose since the studied images exhibit different disturbed surfaces. In addition, the main advantage of using this method is because it gives not one but different ways to interpret the results (different directions).To validate our study and get an accurate understanding of the analysis, the disturbed surface was measured. Then the disturbances rate (%) estimated for each image is correlated with the corresponding Lacunarity and/or Succolarity value (Figure10).Figure 10
Summary of the Lacunarity and Succolarity values as a function of the disturbances rate.According to Figure10, the correlations were statistically significant. The rate of the disturbances varies as a function of the Lacunarity and Succolarity values, suggesting that these fractal variables were good discriminators. Moreover, it has noted that while changes in the surface of the disturbances were not robust, these fractal parameters could reveal significant changes.Considering the result of this study, the fractal analysis may become a suitable and a potential tool to get an accurate diagnosis for differentiating the phosphate deposit at high risk of disturbances and the deposit at low risk.The corresponding analytical procedure can be used as a potential tool that may help the definition of optimum mining strategy management and making decisions. Furthermore, since the challenging question of deposit exploitation feasibility is the quantification of phosphate reserves and resources, this analysis may help the geometricians to improve the estimation of the phosphate reserves in the deposit. This is compared with the traditional geometric methods that have many sources of errors and require much time and effort.Although the above fractal analysis approach was successfully applied and gave important results, we suggest a new perspective of analysis based on the entropy method [38]. The most important reason for discussing the entropy in the context of this study is to consider it as an alternative to the fractal methods for further studies. The idea comes from the inherent relationships that occur between entropy and fractal methods. For example, the Lacunarity and Succolarity analysis, as well as the spatial entropy, can be measured using the same method (the box-counting method). This is interesting as an opportunity to establish a framework for comparative analysis between the different methods.The entropy is an important measure which describes the degree of disorder of systems. It has been widely used as one of the central concepts in several research fields [39, 40]. In the context of image processing, the entropy could be used as a measure of the amount of disorder and complexity. For an image which is completely disrupted the entropy is usually large, whereas for an image which is very highly ordered the entropy is low. The minimum value is zero and it is achieved when image pixels value is constant in any location, that is, all of the pixels have the same gray level while the maximum value depends on the number of grayscales (how many bits used to represent the grayscale).In order to make a spatial analysis of the disturbances, the concept of entropy may successfully be applied. To achieve this task, the images displayed in Figure3 can be used in order to provide a comparative framework of analysis, showing how the different images have different entropy values, the idea being that a more irregular surface must have higher entropy. Consequently, the more disturbed image typically will be associated with higher entropy and vice versa.Furthermore, according to [41], the entropy values may depend on the space filled by the studied object. This suggests that entropy can be used to describe the space filled by the disturbances from one image to another. Hence, the difference in the entropy between different images can be used as an indicator for differentiating the phosphate deposit at high risk of disturbances and the deposit at low risk. The goal is to provide useful means to quantify the disturbances and understand how the entropy index may help in ranking the phosphate deposits.


## 5. Conclusions

Throughout this work, we studied the potential utility of Lacunarity and Succolarity methods as new fractal parameters used to quantify the disturbed areas displayed on geoelectrical images. A detailed examination of eight images was achieved in order to create an interesting comparative framework for analysis. The article also brings the principle and the methodological frame of the Lacunarity and Succolarity methods including formulas by which they were calculated.This study highlights the sensitivity of the fractal analysis methods to differentiate between various disturbed areas. It was shown how the rate of the disturbances is closely associated with the corresponding fractal value. This provides to the best of our knowledge a quantitative means for automatic discrimination between the phosphate deposit at high risk of disturbances and the deposit at low risk.Since the OCP Group managers are eager for developing tools to quantify the disturbed areas, the Lacunarity and the Succolarity have met managers’ needs since they produce results of practical significance. This may lead to an interesting assumption about using fractal analysis as a valuable addition to classical measures in order to get accurate phosphate reserve estimation and make the best exploration planning. Furthermore, this study has opened a window of opportunity for the development of other automatic analysis tools to estimate the rate of disturbances.



---

*Source: 9404567-2019-07-02.xml*