Since the differential element displacement *μ*1 and *μ*2 have the following relationship:

$$\mathbf{C}\_0^2 \frac{\partial^2 \mu\_1}{\partial \mathbf{x}^2} = \frac{\partial^2 \mu\_1}{\partial t^2} \tag{5}$$

Based on the basic theory of wave dynamics, the wave equation in one-dimensional elastic compression bar is derived.

When the bullet hits the incident bar at a certain impact velocity, a compression strain pulse *ε*1(*t*) will be produced in the incident bar. Under the condition of one-dimensional stress propagation, the stress pulse, also known as the elastic stress wave, propagates forward with velocity C0 in the incident bar. When the incident wave propagates to the interface between the incident bar and the rock specimen, because the wave impedance of the rock is less than that of the incident bar, part of the pulse is reflected into the incident bar to form a reflection unloading strain pulse *εR*(*t*), and the remaining part generates a transmitted compression strain pulse *εT*(*t*) in the transmitted bar [28]. Strain gauges adhered to the incident bar and transmitted bar can measure three types of strain pulses *ε*1(*t*), *εR*(*t*) and *εT*(*t*).

The displacement of the left and right ends of the specimen is *μ*1 and *μ*2, respectively. The velocities of the incident wave, reflection wave, and transmitted wave on the end surface are *vI*, *vR*, and *vT*, respectively. The velocities of the left and right ends are *v*1 and *v*2, respectively. According to the theoretical derivation formula:

$$\begin{cases} \boldsymbol{\upsilon}\_{1} = \boldsymbol{\upsilon}\_{I} + \boldsymbol{\upsilon}\_{R} = -\mathsf{C}\_{0}\boldsymbol{\varepsilon}\_{I} + \mathsf{C}\_{0}\boldsymbol{\varepsilon}\_{R} = \mathsf{C}\_{0}(\boldsymbol{\varepsilon}\_{R} - \boldsymbol{\varepsilon}\_{I})\\ \boldsymbol{\upsilon}\_{2} = \boldsymbol{\upsilon}\_{T} = -\mathsf{C}\_{0}\boldsymbol{\varepsilon}\_{T} \end{cases} \tag{6}$$

Displacements are:

$$\begin{cases} \begin{aligned} \mu\_1 &= \int\_0^t \upsilon\_1 dt = \mathbb{C}\_0 \int\_0^t (\varepsilon\_R - \varepsilon\_I) dt \\ \mu\_2 &= \int\_0^t \upsilon\_2 dt = \mathbb{C}\_0 \int\_0^t \varepsilon\_T dt \end{aligned} \end{cases} \tag{7}$$

Assuming that the original length of the specimen is LS, the average strain of the specimen is:

$$
\varepsilon\_S = \frac{\mu\_1 - \mu\_2}{\mathcal{L}\_S} = \frac{\mathcal{C}\_0}{\mathcal{L}\_S} \int\_0^t (\varepsilon\_I - \varepsilon\_R - \varepsilon\_T) dt \tag{8}
$$

Then, the force *F*1 and *F*2 of sample face 1 and face 2 are:

$$\begin{cases} \quad F\_1 = \mathcal{A}\_0 \mathcal{E}\_0 (\varepsilon\_I + \varepsilon\_R) \\ \quad F\_2 = \mathcal{A}\_0 \mathcal{E}\_0 \varepsilon\_T \end{cases} \tag{9}$$

When the cross-sectional area of the specimen is the same as the area of the elastic compression bar, the stress at both ends of the specimen is:

$$\begin{cases} \ \sigma\_1 = \mathcal{E}\_0(\varepsilon\_I + \varepsilon\_R) \\ \ \sigma\_2 = \mathcal{E}\_0 \varepsilon\_T \end{cases} \tag{10}$$

The average stress *σS* in the sample is:

$$
\sigma\_S = \frac{(\sigma\_1 + \sigma\_2)}{2} = \frac{\mathrm{E\_0(\varepsilon\_I + \varepsilon\_R + \varepsilon\_T)}}{2} \tag{11}
$$

In summary, the stress, strain, and strain rate of the specimen are obtained [29]:

$$\begin{cases} \begin{aligned} \sigma\_{\mathcal{S}}(t) &= \frac{1}{2} \mathbf{E}\_{0} (\boldsymbol{\varepsilon}\_{I}(t) + \boldsymbol{\varepsilon}\_{\mathcal{R}}(t) + \boldsymbol{\varepsilon}\_{T}(t)) \\ \boldsymbol{\varepsilon}\_{\mathcal{S}}(t) &= \int\_{0}^{t} (\boldsymbol{\varepsilon}\_{I}(t) - \boldsymbol{\varepsilon}\_{\mathcal{R}}(t) - \boldsymbol{\varepsilon}\_{T}(t)) dt \\ \dot{\boldsymbol{\varepsilon}}\_{\mathcal{S}}(t) &= \frac{\mathbf{C}\_{0}}{\mathbf{L}\_{\mathcal{S}}} (\boldsymbol{\varepsilon}\_{I}(t) - \boldsymbol{\varepsilon}\_{\mathcal{R}}(t) - \boldsymbol{\varepsilon}\_{T}(t)) \end{aligned} \end{cases} \tag{12}$$

The energy carried by incident wave, reflected wave, and transmitted wave can be calculated by integrating the measured strain. The integral formula is as follows:

$$\mathcal{W} = \text{A}0\rho\text{p}\text{C}0 \int\_0^t \varepsilon^2(t)dt\tag{13}$$

The energy carried by incident wave, reflected wave, and transmitted wave is as follows:

$$\begin{cases} \begin{aligned} \boldsymbol{W}\_{I} &= \mathbf{A}\_{0} \rho\_{0} \mathbf{C}\_{0} \int\_{0}^{t} \boldsymbol{\varepsilon}\_{I}^{2} dt \\ \boldsymbol{W}\_{R} &= \mathbf{A}\_{0} \rho\_{0} \mathbf{C}\_{0} \int\_{0}^{t} \boldsymbol{\varepsilon}\_{R}^{2} dt \\ \boldsymbol{W}\_{T} &= \mathbf{A}\_{0} \rho\_{0} \mathbf{C}\_{0} \int\_{0}^{t} \boldsymbol{\varepsilon}\_{T}^{2} dt \end{aligned} \end{cases} \tag{14}$$

According to the law of conservation of energy, the energy loss in the impact process can be obtained, and the energy absorbed by the specimen can be expressed as:

$$\mathcal{W}\_A = \mathcal{W}\_I - \mathcal{W}\_R - \mathcal{W}\_T \tag{15}$$

where *WI* represents the energy carried by the incident wave; *WR* is the energy carried by the reflected wave; *WT* is the energy carried by the transmitted wave; and *WA* is the energy dissipated in the test, that is, the energy absorbed by the failure of the specimen.

#### *2.2. Model and Principle of One-Dimensional Static and Dynamic Combination Loading*

Before the one-dimensional static and dynamic combined loading test, it is necessary to analyze whether the bar and specimen under axial compression meet the stress wave transmitted theory on which the device depends.

As shown in Figure 3, deep rock mass under one-dimensional stress is often subjected to both static stress and dynamic load. Sample micro-force of static and dynamic combination is shown in Figure 4. According to the assumption of one-dimensional stress wave in the bar [30], the force–deformation relationship of the micro-element under combined loading can be obtained:

**Figure 3.** Schematic diagram of static and dynamic combination.

**Figure 4.** Sample micro-force of static and dynamic combination.

$$-\frac{\partial (P\_S + P\_d)}{\partial \mathbf{x}} \Delta \mathbf{x} = \rho\_0 \mathbf{A}\_0 \Delta \mathbf{x} \frac{\partial^2 \mu}{\partial t^2} \tag{16}$$

A0 and *ρ* are the cross-sectional area and density of the elastic bar, respectively; *u* is the displacement of the micro-element after being stressed; and *PS* and *Pd* are the static and dynamic loading of the sample, respectively.

According to stress, strain, and Hooke's Law:

$$\begin{cases} \ \sigma = \frac{P\_S + P\_d}{A\_0} \\ \ \sigma = \mathcal{E}\_0 \varepsilon \\ \ \varepsilon = -\frac{\partial u}{\partial x} \end{cases} \tag{17}$$

Combining the above equations yields:

$$
\rho\_0 \frac{\partial^2 \mu}{\partial t^2} = \mathcal{E}\_0 \frac{\partial^2 \mu}{\partial x^2} \tag{18}
$$

The velocity of incident stress wave in compression bar can be expressed by Formula (4), then Formula (18) can be expressed as:

$$\frac{\partial^2 u}{\partial t^2} - \mathbb{C}^2 \frac{\partial^2 u}{\partial x^2} = 0 \tag{19}$$

The same fluctuation equations derived from the combined dynamic and static loading experimental system and the conventional SHPB test system illustrate the applicability of the one-dimensional stress wave theory to the one-dimensional combined dynamic and static loading experimental system.

#### *2.3. Experimental Procedure*

In this experiment, the loading axial pressure of marble was determined to be 6 MPa, or 10% of its uniaxial compressive strength, based on its uniaxial compressive strength. The impact of velocity was determined by conducting pre-tests on specimens at room temperature. The pre-test results show that when the impact pressure is 0.4 MPa, the specimen is damaged after three times of impact. Therefore, the impact pressure was set as 0.5 MPa, 0.55 MPa, 0.6 MPa, 0.65 MPa in this experiment. The corresponding average impact velocities were 15.32 m/s, 18.17 m/s, 21.83 m/s and 23.49 m/s, respectively. The corresponding average impact velocities were 15.32 m/s, 18.17 m/s, 21.83 m/s, and 23.49 m/s, while the temperature gradients were 25 ◦C, 200 ◦C, 400 ◦C, 600 ◦C, and 800 ◦C, respectively. For this reason, this experiment was separated into four groups, with five pieces in each group undergoing an impact test at a different temperature. Three parallel tests were designed for each temperature, totaling 60 pieces.

#### **3. Physical Properties of Marble before and after High Temperature**

#### *3.1. Apparent Morphological Characteristics of Specimens before and after High Temperature*

The apparent diagram of marble specimens treated at different temperatures (25~800 ◦C) is shown in Figure 5. The diagram shows that the apparent color of marble specimen heated at 200 ◦C is deepened. When the temperature exceeds 400 ◦C, the color of marble specimen surface gradually becomes lighter, that is, from light gray to milky white, and a large number of black spots appear on the specimen surface [31]. At 600 ◦C, the surface of marble becomes very rough, and many microcracks appear on the surface, indicating that the mineral composition of marble has undergone phase transformation, which destroys the original microstructure of the rock. When the temperature reaches 81,000 ◦C, the color of the specimen becomes white, the volume expansion decreases obviously, and the internal structure of marble has been seriously damaged.

**Figure 5.** Apparent morphology of marble specimens after high-temperature treatment.

#### *3.2. Variation in Specimen Mass and Longitudinal Wave Velocity*

Figure 6a–d shows the changes in geometric size and physical properties of marble before and after high-temperature treatment. Figure 6a depicts the relationship between marble sample height and temperature. The specimen height increases with the increase in temperature on the whole. When heated to 200 ◦C~400 ◦C, the specimen height changed little. When heated to 600 ◦C~800 ◦C, the height of the specimen changes noticeably. When the temperature was 800 ◦C, the height of the specimen increased from 50.23 mm before heating to 51.34 mm, with an increase of 2.2%. The fluctuation curve of marble sample diameter with temperature is shown in Figure 6b The variation in the diameter with temperature is similar to that of the diameter with temperature. The diameter of a specimen increases as the temperature rises. At a room temperature of 200 ◦C, the height of the specimen changed little. At 200 ◦C, the diameter of the specimen increased by only 0.23 mm. When heated to 400 ◦C, 600 ◦C, and 800 ◦C, the diameter of the specimen changed drastically. After a high temperature of 800 ◦C, the diameter increased from 52.03 mm before heating to 54.96 mm, with an increase of 5.6%. This is primarily due to the irreversible thermal expansion of the internal part of the specimen at high temperature. Therefore, the specimen size will not return to its original state after cooling. Due to the decomposition of carbonates when the temperature reaches 800 ◦C, the rock strength rapidly declines, and CO2 gas is emitted during the decomposition process, causing the expansion of the height and diameter of high-temperature marble and loose particles [32]. Figure 6c illustrates the variation curve of sample mass and temperature. Between the room temperature and high temperature of 600 ◦C, the mass change in the sample was not obvious, the loss rate was less than 1%. At 800 ◦C, the mass decreased from 243.59 g to 234.82 g, reduced by 3.6%. Due to the thermal decomposition of some minerals in marble caused by the high temperature, the specimen mass changes significantly [33]. The fluctuation curve of marble longitudinal wave velocity with temperature is depicted in Figure 6d The longitudinal wave velocity of the sample decreases with the increase in the overall temperature of the sample. The longitudinal wave velocity of the sample heated to 400 ◦C changed little, only reducing by 8.6%. When heated to 600 ◦C and 800 ◦C, the longitudinal wave velocity decreased significantly by 41.27% and 70.62%, respectively. There are two main reasons for the decrease in rock longitudinal wave velocity after high temperature. On the one hand, when the rock is subjected to high temperature, the free water in the pore evaporates to water vapor and the volume of the pore increases. The pore has a blocking effect on the propagation of longitudinal wave velocity, resulting in a decrease in wave velocity. On the other hand, there are a large number of cracks in the rock, which will lead to the further expansion of cracks and the generation of new cracks under the action of temperature. The higher the temperature, the higher the number of cracks, so the wave velocity of rock samples decreases significantly after high temperature [34]. The above result shows that temperature has a significant impact on the physical properties of marble, and the greater the temperature, the more obvious the impact.

#### *3.3. Influence of High Temperature on the Microstructure of Marble*

The change in the internal microstructure of rock can be observed by scanning electron microscope images. Figure 7 illustrates typical SEM images of granite samples subjected to varying temperatures. For comparison purposes, all images have the same magnification (Mag = 500). Even at room temperature (25 ◦C), marble retains its natural cracks. When marble samples were heated to 200 ◦C, the number of microcracks did neither increase nor decrease. Compared to the marble at room temperature, the crack width and number of specimens increased marginally following heat treatment at 400 ◦C, but the increase was minimal. In the specimens of marble subjected to heat treatments at 600 ◦C and 800 ◦C, the number and width of cracks rose dramatically. At 800 ◦C, cracks even penetrated the specimens, and the thermal damage of specimens was serious.

**Figure 6.** Changes of physical properties of marble samples after different temperatures.

**Figure 7.** SEM micrograph of marble after high-temperature treatment.

#### **4. Dynamic Compressive Mechanical Properties of Marble after High Temperature** *4.1. Stress–Strain Properties*

The total stress–strain curves of marble specimens under different loading speeds at room temperature and after different high temperatures are shown in Figure 8.

It can be seen from Figure 8 that the total stress–strain curves of marble specimens treated at different temperatures (*T*) under different loading rates (*v*) can be roughly divided into four stages, namely fissure compaction stage, elastic deformation stage, plastic deformation stage, and failure stage. In the initial compaction stage, the microcracks in marble tend to close under external dynamic load. Therefore, in this stage, the stress-strain curve slightly upward bending, curve slope increases gradually. Under dynamic load, the stress-strain curves of concave stage are not obvious, but it do exist [35,36]. In the elastic deformation stage, the strain grows as the stress increases, and their correlation is approximately linear. At varying temperatures, the slopes of marble specimens vary. In the plastic deformation stage, as stress increases, the slope of the curve gradually reduces, and

the rate of slope reduction for marble at different temperatures varies. When peak stress is attained, the slope of the curve becomes zero. In the failure stage, the stress-strain curve decreases rapidly, and the slope of the curve is negative. At this time, the bearing capacity of marble decreases.

**Figure 8.** Stress–strain curves in marble after high-temperature treatment.

When the impact velocity is constant, the stress-strain curves of marble vary at different temperatures. Specifically, before 400 ◦C, the stress-strain curves of marble had little difference. At this time, the failure stage curve of the specimen has a "drop" phenomenon, indicating that the brittleness of the specimen is obvious within this temperature range. When the temperature exceeds 400 ◦C, the stress-strain curve of the specimen gradually shifts to the right, and the slope of the curve slows down in the failure stage, indicating that the mechanical properties of marble change from brittleness to plasticity. When the temperature is constant, the stress and strain of marble are similar under different impact velocities. With the increase of impact velocity, the specimen shows the strengthening effect of impact velocity, that is, the peak stress increases with the increase of impact velocity.

#### *4.2. The Variation Pattern of Peak Stress*

Figure 9 shows the relationship between the peak stress, impact velocity, and temperature of the sample. At the same temperature, the peak stress of the marble specimen increases as the impact velocity increases, indicating that the impact velocity has a significant strengthening effect. The dynamic peak stress (*σP*) increases linearly with the impact velocity (*v*), and the fitting relationship is:

$$
\sigma\_P = av + b \tag{20}
$$

where *a* and *b* are the fitting parameters. The values of marble specimens at different temperatures are shown in Table 1.

**Figure 9.** Relationship between the dynamic peak stress and the impact velocity.

Figure 9 demonstrates that, at the same temperature, the peak stress of the specimen increases progressively with increasing impact velocity. When the impact velocity remains unchanged, the specimen's peak stress steadily falls as the temperature rises.

There are two primary reasons for the analysis: first, the high temperature causes the expansion stress inside the marble specimen to increase, resulting in the expansion of the original micro-cracks and the generation of new micro-cracks [37]. Due to the fact that the primary components of marble are calcite (CaCO3) and dolomite (MgCa(CO3)2), high temperature will lead to the degradation and decomposition of dolomite structure, thereby reducing the ability of marble to resist external load damage.

The relationships between dynamic peak stress and the impact velocity of marble specimens at various temperatures are listed in Table 1. It can be seen from Table 1 that the minimum fitting correlation coefficient R2 between dynamic peak stress and the impact velocity of marble specimens at different temperatures is 0.87, indicating that the correlation between them is obvious. In the fitting formula, coefficient a represents the rate of peak stress rise with impact velocity. The greater the value of a, the more pronounced the strengthening effect of impact velocity. In general, with the gradual increase in temperature, the strengthening effect of impact velocity on marble peak stress is lower.

#### *4.3. Crack Extension Process and Damage Mode*

Figure 10 illustrates the axial damage process of marble at a speed of 18.17 m/s and a temperature of 400 ◦C using a high-speed camera.

**Figure 10.** Failure process of specimen taken by high-speed camera.

During the early loading stage of the stress wave, no cracks are detected on the specimen surface. At 50 μs, as a result of the reciprocating propagation of the stress wave in the specimen, an obvious crack, referred to as the major crack, appears on the side of the specimen. At 100~150 μs, the crack length and width gradually expand with time due to the ongoing action of the stress wave. At 200~300 μs, cracks begin to gather and penetrate the whole specimen. At 400 μs, the specimen became unstable, and its bearing capacity drastically dropped. The formation of cracks is mainly due to the dynamic loading induced stress concentration at the crack tip, such that the stress value at the crack tip exceeds the tensile strength of the specimen.

Figure 11 depicts the failure modes of marble specimens at 25~800 ◦C under four impact velocities. At the same impact velocity, the crushing degree at 600~800 ◦C is greater than that at 25~400 ◦C. It shows that temperature has a significant effect on the fracture characteristics of the specimen. When the impact velocity is low (15.32 m/s ≤ *v* ≤ 618.17 m/s), the specimen can still maintain a certain bearing capacity at temperatures between 25 ~600 ◦C. At 800 ◦C, the specimen is unstable. At a higher impact velocity (21.83 m/s ≤ *v* ≤ 23.49 m/s), at 25 ◦C ≤ *T* ≤ 400 ◦C, the bulk rate of the specimen is higher and the average particle size is larger. When 600 ◦C ≤ *T* ≤ 800 ◦C, the specimen suffered from crushing failure, with small particle size and relatively uniform distribution.

**Figure 11.** Failure mode of marble after high temperatures.

At the same temperature, with the increase in impact velocity, the fracture surface of the specimen gradually expands, the degree of fragmentation increases, and the fragment size decreases [21].

#### *4.4. Energy Analysis of Rocks under Combined Dynamic and Static Loading*

In rock engineering, the excavation, fragmentation, and disturbance of rock mass inevitably involve the inflow, accumulation, dissipation, and outflow of energy, and the energy changes throughout the entire rock deformation and failure process. Therefore, it is crucial to examine the failure deformation of rock from the perspective of energy [38]. In the SHPB test, The incident energy (*WI*), absorption energy (*WA*), transmitted energy (*WT*), and reflection energy (*WR*) can be calculated using Formulas (14) and (15). Figure 12 depicts the relationship between incident energy, absorbed energy, transmitted energy, reflected energy, and temperature at various impact velocities.

**Figure 12.** Change of energy with temperature under different impact velocity.

Figure 12a–d shows the relationship between the energy and temperature of the specimen under different impact velocities. When the impact velocity is constant, the incident energy increases first and then decreases with the increase of temperature, and reaches the maximum at 400 ◦C. The reflected energy first decreases and then increases with the increase of temperature, and the two are quadratic functions with an upward opening. The transmission energy and absorption energy first increase and then decrease with the increase of temperature, and the two are quadratic functions with an open downward direction. With the increase of impact velocity, the incident energy, reflection energy, transmission energy and absorption energy of the specimen increase. The fitting relationship between temperature and energy at different impact velocities is listed in Table 2.

**Table 2.** Fitting relationship between temperatures and energies at different impact velocities.

#### **5. Conclusions**

At room temperature and high temperatures (25~800 ◦C), the fundamental physical properties of marble were determined. SHPB impact compression test was used to measure the effect of temperature on the dynamic characteristics of marble.

**Author Contributions:** Conceptualization, methodology, software, validation, formal analysis, writing—original draft, writing—review and editing: Y.W.; investigation, data curation: X.Z.; supervision, software: L.H.; supervision, methodology: J.W.; funding acquisition, project administration: X.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was financially supported by the National Natural Science Foundation of China (No. 52164010, No. 52064025), the General Project of Basic Research Program of Yunnan Province (No. 202201AT070178), and Yunnan major scientific and technological special project (202102AG050024). Their support is gratefully appreciated.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data sharing not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

## *Article* **Development of Predictive Models for Determination of the Extent of Damage in Granite Caused by Thermal Treatment and Cooling Conditions Using Artificial Intelligence**

**Naseer Muhammad Khan 1,2,3, Kewang Cao 1,4,\*, Muhammad Zaka Emad 5,\*, Sajjad Hussain 6, Hafeezur Rehman 3,7, Kausar Sultan Shah 7, Faheem Ur Rehman 8 and Aamir Muhammad 9**

**Abstract:** Thermal treatment followed by subsequent cooling conditions (slow and rapid) can induce damage to the rock surface and internal structure, which may lead to the instability and failure of the rock. The extent of the damage is measured by the damage factor (*DT*), which can be quantified in a laboratory by evaluating the changes in porosity, elastic modulus, ultrasonic velocities, acoustic emission signals, etc. However, the execution process for quantifying the damage factor necessitates laborious procedures and sophisticated equipment, which are time-consuming, costly, and may require technical expertise. Therefore, it is essential to quantify the extent of damage to the rock via alternate computer simulations. In this research, a new predictive model is proposed to quantify the damage factor. Three predictive models for quantifying the damage factors were developed based on multilinear regression (MLR), artificial neural networks (ANNs), and the adoptive neural-fuzzy inference system (ANFIS). The temperature (*T*), porosity ( *ρ* ), density (*D*), and P-waves were used as input variables in the development of predictive models for the damage factor. The performance of each predictive model was evaluated by the coefficient of determination (R2), the A20 index, the mean absolute percentage error (MAPE), the root mean square error (RMSE), and the variance accounted for (VAF). The comparative analysis of predictive models revealed that ANN models used for predicting the rock damage factor based on porosity in slow conditions give an R2 of 0.99, A20 index of 0.99, RMSE of 0.01, MAPE of 0.14, and a VAF of 100%, while rapid cooling gives an R2 of 0.99, A20 index of 0.99, RMSE of 0.02, MAPE of 0.36%, and a VAF of 99.99%. It has been proposed that an ANN-based predictive model is the most efficient model for quantifying the rock damage factor based on porosity compared to other models. The findings of this study will facilitate the rapid quantification of damage factors induced by thermal treatment and cooling conditions for effective and successful engineering project execution in high-temperature rock mechanics environments.

**Keywords:** predictive models; damage factor; thermal treatment; computer simulations; ANNs; ANFIS

**MSC:** 68T07

**Citation:** Khan, N.M.; Cao, K.; Emad, M.Z.; Hussain, S.; Rehman, H.; Shah, K.S.; Rehman, F.U.; Muhammad, A. Development of Predictive Models for Determination of the Extent of Damage in Granite Caused by Thermal Treatment and Cooling Conditions Using Artificial Intelligence. *Mathematics* **2022**, *10*, 2883. https://doi.org/10.3390/ math10162883

Academic Editor: Elena Benvenuti

Received: 22 February 2022 Accepted: 6 August 2022 Published: 11 August 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### **1. Introduction**

Temperature is an essential consideration that has a significant impact on a rock's chemical, physical, and mechanical properties. Underground coal gasification, geothermal resources, nuclear waste disposal, coal mine gas explosions, underground engineering, fire reconstruction, and improved oil recovery are examples of rocks being exposed to high-temperature conditions [1–4]. The interactions of rocks in these projects take place for a long time and they are continuously exposed to high-temperature ranges, from 500 to 1500 ◦C [5–16]. The long-term exposure to high temperatures yields voids, pores, and microcracks. It also propagates the lengths of existing microcracks, causing damage to the integrity and stability of rocks. Moreover, it further impacts the physical, chemical, and mechanical characteristics of rocks [4,17–23]. Therefore, it is imperative to evaluate the damages thoroughly induced in rocks from high temperatures or thermal treatments for the safe execution of engineering projects. Comprehensive investigations in research studies have been carried out on different rocks subjected to thermal treatments (for various time exposures and under different cooling conditions) to evaluate rock damage mechanisms, thermal cracking, deformation mechanisms, thermal-induced stresses, strength reductions, and changes in the physical properties under high temperatures [1,3,24–38]. However, long time exposures of rocks and vast contrasts in thermal interactions make it essential to evaluate and quantify the extent of thermal damage for the safe execution of engineering projects.

Numerous researchers have investigated the degree of damage using various approaches. Placido [39] studied concrete samples that were thermally treated up to 500 ◦C by thermoluminescence (*TL*) to measure the degree of damage. Similarly, Chew [40] performed a number of experiments, such as visual observations, a rock sound test, the Schmidt hammer test, and ultrasonic pulse velocity, by thermally treating (up to 500 ◦C) and non-treating concrete samples, and quantifying the extent of damage by *TL*. These findings are instrumental in determining the degree of damage to concrete and rock in civil and mining engineering structures, such as tunnels, highways, buildings, etc. [26,41].

The degree of damage was measured by investigating the development and propagation of microcracks using acoustic emission signals, ultrasonic wave velocities, and strain measurements [42–44]. The degree of damage was quantified using a thin section analysis, scanning electron microscope (SEM), and micro-CT scanning [45,46]. However, the aforementioned techniques for quantifying the extent of thermal damage require dedicated equipment for measuring and observing changes in strain, ultrasonic wave velocities, acoustic emission signals, and the internal morphology of a rock specimen. The rock specimens needed for these tests could be prepared according to ISRM Ulusay [47–52] using sophisticated equipment, which could cause a delay in the quantification process.

Based on the required testing time and economics, researchers generally prefer to use empirical, statistical, and other machine learning techniques to predict or estimate the required outputs. These techniques have been developed based on the concept of mathematics. These techniques, due to their versatile nature, have been successfully applied in various engineering fields. Nowadays, these techniques have gained more attention in solving complex rock engineering problems. In this regard, various researchers use statistical techniques, including simple regression (SR), multiple linear regression (MLR), and artificial intelligence (AI) techniques comprising artificial neural networks (ANNs) and the adaptive neuro-fuzzy inference system (ANFIS) to predict the strength and deformational properties of rocks using physical properties. These properties include density, porosity, and ultrasonic wave velocities as input variables [53–65]. These studies suggest that artificial intelligence technique prediction performances are superior to statistical techniques. Moreover, the above-mentioned literature studies show that statistical approaches, such as simple regression analysis (SRA) and multivariate regression analysis (MVRA), have been used to demonstrate a link between physical and mechanical parameters. Additionally, research studies have utilized soft computing methods, such as fuzzy inference system(s) (FIS), artificial neural network (ANNs), their combination, and the adaptive neuro-fuzzy

inference system (ANFIS) to forecast the output parameters (mechanical characteristics) based on the input parameters (physical properties). The findings of these investigations indicate that, when compared to statistical methodologies, soft-computing tools are more accurate for predicting mechanical characteristics. Additionally, Sirdesai, et al. [66] concluded that an ANFIS model was shown to have a higher prediction efficiency than ANNs. Statistical and soft computing techniques, on the other hand, have been utilized to forecast the strength and elastic characteristics of untreated specimens. Sirdesai, et al. [66] predicted the degree of thermal damage from the elasticity modulus using MLR, ANN, and ANFIS techniques on thermally-treated and slow cooling fine-grained Indian sandstone, up to 1000 ◦C, using porosity ( *φ* ), density ( *ρ* ), coefficients of linear and volumetric thermal expansion (EL and EV), and wave velocities (VP and VS). After a comparative analysis of models, they proposed that the ANFIS model is the most suitable for the mentioned purpose because of its better performance over MLR and ANN.

After a detailed study of the literature, it was noticed that there is a lack of research on the prediction of the extent of damage/damage factors from porosity under slow and rapid cooling, the selection of the optimized neuron for the best results, comparative analysis of different algorithm functions in AI, and impoverished performances of ANNs. Thus, it was imperative to address the mentioned gaps by conducting innovative research on rocks other than sandstones using statistical and artificial intelligence to predict the extent of thermal damage.

In this research, the degree of thermal damage *DT* of granitic rocks was predicted under the slow and rapid cooling of thermally-treated granite. The physical properties, such as porosity (Φ), density ( *ρ* ), temperature, and P-wave velocity (*PV*) were used as input variables for multilinear regression (MLR) and artificial intelligence (ANN and ANFIS) techniques. The adequacy of each model was evaluated based on the mean absolute percentage error (MAPE), coefficient of determination (R2), A20 index, the variance accounted for (VAF), and root mean square error (RMSE). The most effective model was proposed to predict the extent of thermal damage for granitic rocks. The novelty of research includes the prediction of damage extent based on porosity under cooling conditions, neuron optimization, comparative analysis of different algorithm functions, and monitoring the impoverished performances of ANNs as compared to ANFIS.

#### **2. Design of the Experimentation Process**

#### *2.1. Sample Preparation*

In this study, granite rock samples were used, which were collected from a quarry located in the Baba G Kandaw district Buner, Khyber Pakhtunkhwa, Pakistan, as shown in Figure 1. The cores were extracted from the bulk rock samples and waxed to preserve their initial mechanical properties and avoid mineralogical and size deterioration. The cylindrical core specimens (with dimensions of 54 × 108 mm) were prepared with high geometric integrity. The ends of each core were polished carefully with a grinding machine until the deviation in the flattening of the core end became less than 0.5–0.05 mm [67]. The cores were then heated to the selected temperatures (25 ◦C, 300 ◦C, 600 ◦C, and 900 ◦C) at a constant rate of 5 ◦C/min [68], followed by cooling to room temperature. Cooling was performed in two different ways: exposed to the air for slow cooling and placed in water for rapid cooling. Based on the cooling conditions, the rock samples were categorized into seven groups, namely, An, Bn, Cn, Dn, En, Fn, and Gn. The samples were kept as reference samples; the temperature was assumed as 25 ◦C; Bn, Cn, and Dn samples were used for slow cooling, and En, Fn, and Gn were used for rapid cooling. According to detailed petrographic investigations of the granite samples, they were mostly composed of perthite feldspar, plagioclase, and quartz, with little biotite, muscovite, and opaque oxides and sulfides. The average rock composition included: K-feldspar (47.29%), quartz (25.18%), plagioclase (24.38%), muscovite (0.15%), biotite (32.03%), and other (0.97%).

**Figure 1.** (**a**) Pakistan map and regional geological map showing the distribution of the Cambrian to Ordovician Swat Granite Gneisses, metasedimentary rocks, ranging in age from Late Proterozoic to Late-to-Middle Mesozoic and Late Carboniferous to Permian Ambela and the Shewa Igneous complex; (**b**) granite boulder; (**c**) micrographs of the thin sections (cross-polarized light).

#### *2.2. Experimental Procedure and Instrument*

A portable ultrasonic nondestructive digital indicating tester (PUNDIT) was used to compute the ultrasonic parameter, such as the ultrasonic P-wave velocity, ultrasonic S-wave velocity, and the transit time for each velocity. It is worth mentioning that the test procedure in this part tested the P-wave velocity of the unheated and post-thermal treatments under different cooling conditions of rock samples. The density and porosity of rock specimens were determined before and after heating. The samples were prepared and tested according to the International Society for Rock Mechanics (ISRM). Ten experimental runs for each physical property (density, porosity, P-waves, and elastic modulus) were conducted before and after the thermal treatment and subsequent cooling conditions. The study steps are shown in Figure 2. The average results obtained from testing the thermally-treated granite rock samples and subsequent cooling conditions (slow and rapid cooling) for density, porosity, and P-waves, are presented in Table 1.

**Table 1.** Thermal damage factors under slow cooling and rapid cooling.

**Figure 2.** Flowchart of the methodology for the prediction of thermal damage.

#### *2.3. Thermal Damage Factor*

The thermal damage factor (*DT*) was used to evaluate the degree of damage induced by the thermal treatment and subsequent cooling conditions. It is a key parameter and can induce rock instability and failure when the severity of the damage is maximum. The thermal damage factor was calculated based on porosity and elasticity using Equation (1) and Equation (2), respectively. The calculated values are given in Table 1.

$$D\_{TP} = \left(1 - \frac{1 - n\_T}{1 - n\_{T\_0}}\right) \times 100\% \tag{1}$$

where *nT* is the porosity after temperature and *nTo* represents the porosity before temperature.

$$D\_{TE} = \left(1 - \frac{E\_T}{E\_0}\right) \tag{2}$$

where *ET* is the elasticity at high temperature and *E*0 is the elasticity at room temperature.

#### **3. Prediction Model**

#### *3.1. MLR Model*

MLR is commonly used to predict the relevant parameters. MLR is an extended version of the simple linear regression used in the multiple predictive variables. It can model the input without variables considering their relationship and form a generalized equation, as shown in Equation (3) [69,70].

$$\mathcal{W} = \mathbb{C} + b\_1 z\_1 + b\_2 z\_2 + b\_3 z\_3 + \dots \dots \dots \dots + b\_n z\_n \tag{3}$$

where the partial regression coefficients are *b*1 to *bn*, *W* is the dependent variable, *C* is constant, and *z*1 to *zn* are the independent variables.

#### *3.2. ANN Model*

Numerous AI techniques are used globally for prediction; ANN is one of them. It mimics the behavior of the human brain. It is generally a useful tool in pattern recognition, clustering data, and fitting a function. Because of its learning capacity, memory simulation, and excellent performance owing to features such as categorizing and filtering noisy data, ANN is a particularly significant sector in geotechnical and mining engineering [64,66,71–79]. Furthermore, it is a promising approach used for solving complicated engineering issues involving enormous amounts of data or several input parameters, making manual solutions more difficult. In general, an ANN is made up of components, such as inputs, outputs, weights, activation, training, and numerous neurons. Experiment-collected test data are multiplied by weights and applied to the existing activation functions. [77,79].

Mathematically, the basic ANN is expressed as

$$\mathbf{N} = f(\mathbf{K}\mathbf{x} + \mathbf{C})\tag{4}$$

where

$$K = K\_{1\prime}K\_{2\prime}K\_{3\prime}K\_{4\prime}\dots\dots\dots\dots\dots\_{\prime N\_n}$$

$$M = m\_{1\prime}m\_{2\prime}m\_{3\prime}m\_{4\prime}\dots\dots\dots\dots\dots\dots\_{\prime m\_n}$$

where *K*, *m*, and *C* refer to weights, input, and bias, respectively. The net (*L*) predicted values are calculated using Equation (5)

$$L = \sum\_{i=1}^{n} (\mathbb{K}\_i \mathcal{M}\_i + \mathbb{C}) \tag{5}$$

In this study, tangent sigmoid was used as the transferred function, which was calculated using Equation (6)

$$y = \tanh\left(L\right)\tag{6}$$

$$output = y = \tanh(L) = \tanh\left(\sum\_{i=1}^{n} (K\_i M\_i + \mathbb{C})\right) \tag{7}$$

Generally, the error of the network expresses the difference between the actual and predicted values. This error, affected by the number of neuron weights in the hidden layers (and its value), either increased or decreased. The error of number points (*En*) can be calculated by using Equation (8)

$$E\_n = Actualul\_{value} - Predicted\_{value} \tag{8}$$

The total error (*ET*) can be calculated by using Equation (9)

$$E\_T = \frac{1}{2} \sum\_n E\_n^2 \tag{9}$$

The efficiencies of the ANN networks were assessed using various learning algorithms. The terms "learning algorithm" and "training algorithm" are used interchangeably. In this research, regarding the five learning functions, each function had its advantages and disadvantages. These functions were Levenberg–Marquardt (LM), BFGS quasi-Newton (BFG), resilient backpropagation (RP), scaled conjugate gradient (SCG), and conjugate gradient with Powell/Beale restarts (CBG); they were used as training functions. Additionally, the performance of ANN was greatly affected by hidden-layer neurons. The performance of the ANN network was evaluated by using feed-forward backpropagation type of neural networks, the tangent sigmoid activation function, and numerous neurons. To achieve

optimum performance of the models, the normalization technique was used by applying Equation (10).

$$X\_{norm} = \frac{(X\_{actual} - X\_{min})}{(X\_{max} - X\_{min})} \tag{10}$$

where *Xnorm* refers to the normalized value, *Xactual* refers to the measured value, *Xmin* and *Xmax* refer to the corresponding minimum and maximum values of the dataset. The different numbers of neurons were used for different training functions followed by comparing the models; the best model was recommended in terms of the training function.

#### *3.3. Adaptive Neural Fuzzy Interface System (ANFIS) Models*

The fuzzy logic system was proposed by Zadeh [80] for the first time, which can be used to predict the solution mechanism for complex engineering problems. Owing to the shortcoming, it could not define a standardized procedure for designing such a system. After the advent of the neural network, Jang [81] introduced a new technique known as an adaptive neuro-fuzzy inference system (ANFIS). This technique mapped the output and input by exploiting the fuzzy systems and neural network learning and reasoning proficiencies. Using a combination of fuzzy logic and neural networks, the ANFIS can effectively solve various complex and non-linear problems in any engineering field. ANFIS uses fuzzy rules to predict the output from inputs; these fuzzy rules are developed during the training process. The ANFIS construct and its FIS membership are derived from training data. Two FISs commonly used Mamdani and Sugeno. The key distinction between the two is that the Sugeno output membership is linear or constant, while the Mamdani output membership is triangular, Gaussian, etc. In the present study, Sugeno FIS was used because it is computationally more efficient than Mamdani. The procedure of ANFIS can be described in the form of FIS, with two inputs (S) and (T) and one output (Z). Subsequently, if–then rules of the two-fuzzy were developed, which are given below:

$$\text{Rule 1: If } \mathbf{S} \text{ is } \mathbf{J}\_1 \text{ and } y \text{ is } \mathbf{P}\_1, \text{ then } \mathbf{Z}\_1 = \mathbf{D}\_1 \mathbf{x} + \mathbf{F}\_1 \mathbf{y} + \mathbf{Q}\_1.$$

Rule 2: If *x* is J2 and *y* is P2, then z2 = D2*x* + F2*y* + Q2

where J1, P1, J2, and J2 refer to the input membership functions for 'S' and 'T'; whereas D1, F1, Q2, D2, F1, and Q2 refer to the output function parameters.

Numerous scholars have detailed the ANFIS model, consisting of five layers [66,82–84]. Moreover, the ANFIS model uses two learning algorithms, backpropagation and hybrid, to optimize the result with the minimal value of error between the predicted and estimated values [79–84]. A hybrid optimization method was used in the present study due to its high prediction results [61,77]. This model shows excellent performance, but the convergence of the model was slow due to the high if–and–or relationship.

#### **4. Results and Discussion**

#### *4.1. MLR Models*

The estimated thermal damage factor based on elasticity and porosity under both cooling conditions are given in Equations (11)–(14).

$$D\_{TER} = 50.6946 + 0.0005T - 0.0032\ \rho \ - 0.0185\ \phi \ - 0.0002P\_V \tag{11}$$

$$D\_{TPR} = 211.0464 + 0.00048T - 0.0265 \text{ } \rho \\ - 0.0789 \text{ } \phi \\ + 0.0001 P\_V \tag{12}$$

$$D\_{TES} = -61.34 + 0.00094T + 0.06195\ \rho \,\, + 0.0214\ \phi\,\, + 0.00093P\_V\tag{13}$$

$$D\_{TPS} = 162.9 + 0.00094T - 0.04117\,\text{\AA}\,\text{ }\rho - 0.06099\,\text{\phi}\,\text{ } + 0.00016P\_V\tag{14}$$

Graphically, the experimental and predicted thermal damage factors under each cooling condition are given in Figure 3. Figure 3a,b illustrate that the experimental and predicted thermal damage factors based on porosity and elasticity under slow cooling conditions yielded R2 values of 0.97 and 0.94, respectively. Similarly, the RMSEs for porosity and elasticity were 0.061 and 0.083, respectively. On the other hand, Figure 3c,d show the experimental and predicted thermal damage factors based on porosity under rapid cooling resulting in R2 and RMSE as 0.97 and 0.056, respectively. In the elasticity case, the R2 and RMSE were 0.94 and 0.076, respectively.

**Figure 3.** Experimental and MLR models; (**a**) porosity-based *DT* under slow cooling; (**b**) elasticitybased *DT* under slow cooling; (**c**) porosity-based *DT* under rapid cooling; and (**d**) elasticity-based *DT* under rapid cooling.

#### *4.2. ANN Models*

#### 4.2.1. Model Design

A total of 2000 networks were generated for various training algorithms, each case comprising 500 networks. In the presence of a tangent sigmoid function as an activation function, the feed-forward backpropagation network with 5 training functions (and up to 100 neurons) was tested in the design for each learning algorithm. The model performance and comparative analysis of the different models were evaluated based on R2, the a20 index, RMSE, MAPE, and VAF.

#### 4.2.2. ANN Code Compilation in MATLAB

This study compiled self-generated code for ANN for n numbers of networks, keeping the same training and activation functions for a single loop as shown in Figure 4. A loop function was introduced in this code, which can run for the desired number of networks. The activation function overall in this code was fixed, which can be changed according to the data nature. In the present case, the code was executed for a hundred networks in a single execution. The number of neurons increased in each successor for each network in a loop; i.e., for network1, there was one neuron, for network2, there were two neurons, and so on. In both the hidden and output layers, the same activation function was used.

**Figure 4.** Flowchart of ANN for thermal damage model prediction.

#### 4.2.3. Network Phases and Regression Models

In this research, the basic structure consisted of four inputs (temperature, porosity, density, and P-waves) and one output (thermal damage factor) in both cooling conditions, as was based on porosity and elasticity, as shown in Figure 5. A total of forty data points were taken as a dataset. The dataset was divided into three parts: training (75%), testing (15%), and validation (15%). Figure 6a,b show the training, validation, and testing for the slow cooling thermal damage factor based on porosity and elasticity. Similarly, Figure 6c,d show a rapid cooling thermal damage factor based on porosity and elasticity. Furthermore, prior to comparing various ANN models of thermal damage in both cooling conditions, it was desired to choose the best optimum training function at the first stage for the studied data. Therefore, five different training functions, namely BFG, RP, SCG, CGB, and LM, were evaluated to select the best performance training function for the optimum ANN prediction model.

**Figure 6.** *Cont*.

**Figure 6.** Regression model *DT*, (**a**) porosity-based *DT* under slow cooling, (**b**) elasticity-based *DT* under slow cooling, (**c**) porosity-based *DT* under rapid cooling, and (**d**) elasticity-based *DT* under rapid cooling.

#### 4.2.4. Model Performance

Table 2 shows the optimum results of each algorithm. In comparison to other functions, the overall efficiency of the LM function was very high in terms of R2, RMSE, the number of neurons, and execution timing. As compared to other algorithms, LM converged data faster. Furthermore, the efficiency of LM for DT based on porosity under slow and fast cooling was better than the elasticity-based (in terms of RMSE).

**Table 2.** Performance evaluation of different training algorithms for *DT* under both cooling conditions.

Further, the LM performance was better in slow cooling than in rapid cooling. Figure 7 shows the performance of LM with variations in the number of neurons for the thermal damage factors in both slow and rapid cooling based on porosity and elasticity. The optimum neurons for LM with high R2 and low RSME with the least convergent time details are presented in Table 2, revealing that the optimum neuron numbers in slow cooling based on porosity and elasticity were 80 and 52, respectively. Similarly, for rapid cooling based on porosity and elasticity, the optimal neuron numbers were 18 and 72, respectively.

**Figure 7.** LM optimum number of neurons for (**a**) porosity-based *DT* under slow cooling, (**b**) elasticitybased *DT* under slow cooling, (**c**) porosity-based *DT* under rapid cooling, and (**d**) elasticity-based *DT* under rapid cooling.

Figure 6a shows that the thermal damage factor based on porosity under slow cooling conditions revealed correlation coefficients of 0.99, 0.98, 0.99, and 0.99 for training, validation, testing, and overall, respectively. Similarly, Figure 6b shows the thermal damage factor based on elasticity under slow cooling conditions revealed correlation coefficients of 0.99, 0.98, 0.97, and 0.99 for training, validation, testing, and overall, respectively. In contrast, for rapid cooling conditions, the thermal damage factors based on porosity and elasticity are shown in Figure 6c,d. These figures also revealed correlation coefficients of 0.99, 0.98, 0.97, 0.98 and 0.99, 0.98, 0.92, 0.97, and 0.98 for training, validation, testing, and overall, respectively.

#### 4.2.5. ANN Predicted Models

The effectiveness of the developed ANN models was assessed by comparing the predicted and actual values, as shown in Figure 8. The porosity and elasticity-based damage factors in slow cooling conditions are shown in Figure 8a,b, with the corresponding correlation coefficient value as 0.99; likewise, the RMSE value was recorded as 0.01 and 0.07 for porosity and elasticity-based damage factor, respectively. It indicates that the porosity-based prediction model outperformed the elasticity-based thermal damage factor in slow cooling conditions. In addition, the thermal damage factor based on porosity in rapid cooling (Figure 8c) showed the correlation coefficient and RMSE value as 0.99 and 0.02, respectively.

**Figure 8.** Experimental ANN model; (**a**) porosity-based *DT* under slow cooling, (**b**) elasticity-based *DT* under slow cooling, (**c**) porosity-based *DT* under rapid cooling, and (**d**) elasticity-based *DT* under rapid cooling.

Similarly, Figure 8d shows the elasticity-based damage factor in rapid cooling conditions that resulted in the correlation coefficient and RMSE values of 0.99 and 0.09, respectively. It is worth mentioning that the damage factor dependent on porosity and elasticity in both cooling conditions had almost the same high coefficient of determination and less RMSE than the elasticity-based value, as shown in Figure 8. Furthermore, the porosity-based prediction model is more accurate in terms of R2, RMSE, the number of neurons, and convergent time.

#### *4.3. ANFIS Models*

Similar to ANN, the dataset was divided into three parts, training (75%), testing (15%), and validation (15%). The flowchart of the ANFIS for thermal damage is shown in Figure 9. The data division pattern was used in both cooling conditions to measure the thermal damage factors (porosity, elasticity). The ANFIS models were trained for up to 50 epochs. The fuzzy interface system (FIS) was generated for models using a sub-clustering algorithm and a hybrid training algorithm for FIS optimum methods. The linear Gaussian function was used to predict the thermal damage factor from input under both cooling conditions. Figure 10 shows the structure of the developed ANFIS model. Table 3 describes the details of each parameter used during the model's development.

**Figure 9.** ANFIS flowchart for thermal damage prediction.

**Figure 10.** Proposed ANFIS model for thermal damage factor.

**Table 3.** The properties of the proposed ANFIS model under both cooling conditions.

Note: mean valid.

The efficiencies of the developed ANFIS models were analyzed by comparing the predicted value to the actual value, as shown in Figure 11. Figure 11a shows the porositybased *DT* in slow cooling, which reveals R2 and RMSE values of 0.98 and 0.07, respectively. Similarly, Figure 11b shows the elasticity-based *DT* in slow cooling, revealing R2 and RMSE values of 0.97 and 0.44, respectively. In contrast, Figure 11c shows the porosity-based *DT* in rapid cooling, which reveals R2 and RMSE values of 0.98 and 0.88, respectively, while Figure 11d shows the porosity-based *DT* in rapid cooling, which shows R2 and RMSE values of 0.94 and 0.88, respectively. Furthermore, the prediction damage based on porosity showed higher R2 and lower RMSE values than the elasticity-based. Hence, porosity-based *DT* in both cooling conditions gives more reliable results than elasticity-based *DT* in terms of high (R2, a20 index) and low RMSE values.

**Figure 11.** Experimental and ANFIS models; (**a**) porosity-based *DT* under slow cooling, (**b**) elasticitybased *DT* under slow cooling, (**c**) porosity-based *DT* under rapid cooling, and (**d**) elasticity-based *DT* under rapid cooling.

#### **5. A Comparative Appraisal of Statistics and Intelligent Technique**

The comparison of correlation efficiencies of various developed models was used in this study to improve the performances of the predicted models. The subsequent performance indices, such as R2, MAPE, RMSE, and VAF, were evaluated. An excellent model can be represented by the following performance indices: R2 = 1, a20 index, MAPE = RMSE = 0, and VAF = 100%. The performance indices were calculated using Equations (15)–(19).

$$\mathcal{R}^2 = \frac{\sum\_{i=1}^n \left( y\_i \right)^2 - \sum\_{i=1}^n \left( y\_i - k\_i' \right)^2}{\sum\_{i=1}^n \left( y\_i \right)^2} \tag{15}$$

$$\text{MAPE} = \frac{1}{2} \sum\_{i=1}^{n} \left| \frac{y\_i - k\_i'}{y\_i} \right| \times 100 \tag{16}$$

$$\text{RMSE} = \sqrt{\frac{\sum\_{i=1}^{n} (y\_i - k\_i')}{n}} \tag{17}$$

$$\text{VAF} = \left[1 - \frac{\text{var}(y - k')}{\text{var}(y)}\right] \times 100\tag{18}$$

$$\text{a20} - \text{index} = \frac{k^{20}}{M} \tag{19}$$

where *y* is the actual value, *k* is the predicted value, *k*20 is the ratio of the original and predicted values in the range of 0.80–1.20, and *M* is the total datasets.

As a result, the LM-based ANN model was selected compared to LMR and ANFIS. Table 4 illustrates the performance indices. The performance index values demonstrate that the ANN model performed better than the MLR and ANFIS approaches in the current study. The adequacy of ANN is higher than ANFIS and MRL. The most significant limitation of ANFIS over ANN observed in the present study was that the FIS generator training took a long time, particularly as the number of inputs and epochs increased, while ANN executed too fast. Additionally, Figure 12 depicts the predicted and experimental thermal damage factor dependent on porosity and elasticity under both cooling conditions. In both cases, the ANN showed a better prediction comparatively than ANFIS and MLR. Furthermore, all models gave high accuracies at low and high temperatures (below 200 ◦C and greater than 600 ◦C). This fluctuation in prediction and measure value was due to a nonlinear increase in thermal damage. The ANN model is better than other models and overlaps with the experimental curve, as shown in Figure 12.

**Table 4.** Performance indices of the developed models.

**Figure 12.** Experimental ANN, MLR, and ANFIS models; (**a**) porosity-based *DT* under slow cooling, (**b**) elasticity-based *DT* under slow cooling, (**c**) porosity-based *DT* under rapid cooling, and (**d**) elasticitybased *DT* under rapid cooling.

#### **6. Limitations and Future Works**

Khan et al. [77,78] claim that rock behavior varies depending on the region. In this study, the granite rock of a specific area was used for thermal damage. The study can be generalized by considering multiple rocks in different areas. However, in future research, care should be taken regarding thermal damage prediction as the rock behaviors are very sensitive and depend on multiple parameters, including mineralogy, physical, and mechanical properties. The proposed model of this study can also predict DT when rock input parameters, such as temperature, density, porosity, and P-waves, are available in the same range. To increase forecasting accuracy, future models should be trained with more datasets. This research focuses on traditional linear regression model(s) (MLR) and two artificial intelligence approaches (AI) (ANN and ANFIS). Using other approaches, such as decision tree, random forest (RF), K-nearest neighbor (KNN), ANN, and the ANFIS model to anticipate DT values may be investigated in the future. We may also present a greater database of non-destructive rock index tests to provide a more sophisticated intelligence approach since generalization is a crucial feature of predictive models.

#### **7. Conclusions**

The study investigated the damage factors from the porosity and elastic modulus under slow and rapid cooling. It was observed that the extent of damage to granite increased with the increase in temperature. However, the predominant damages to granite rock were observed in rapid cooling compared to slow cooling. Three predictive models were developed to quantify the extent of damages induced by thermal treatment and, subsequently, cooling conditions based on MLR, ANN, and ANFIS. In order to predict thermal damage, physical properties, such as temperature, density, porosity, and P-waves, were used as input variables. The efficiencies of the models were evaluated based on R2, the A20 index, RMSE, MAPE, and VAF. The optimum training function was determined based on five different training functions, namely BFG, RP, SCG, CGB, and LM, for porosity and elasticity under different cooling conditions, to achieve the best ANN prediction model. It was revealed that the optimum result of the LM algorithm performance is better in terms of R2, RMSE, the number of neurons, and execution timing. After a comparative analysis of predictive models, it has been suggested that the ANN-based predictive model is more efficient in the prediction of damage factors as compared to MLR and ANFIS. It has also been concluded that the efficacy of the prediction for the damage factor based on porosity is more predominant than the prediction of the damage factor based on elasticity. Therefore, it is recommended that damage factor prediction based on porosity should be used in the future. The findings of this study will facilitate the rapid quantification of damage factors induced by thermal treatment and cooling conditions for effective and successful engineering project execution in high-temperature rock mechanics environments. In the future, we will evaluate thermal damage in the presence of infrared radiation characteristics, acoustic emissions, and AI applications. These techniques will be adopted to obtain the real-time damage factor, crack intensity, propagation, and direction due to thermal heat.

**Author Contributions:** N.M.K. and S.H. contributed to the research, designed the experiments, and wrote the paper. K.C. and M.Z.E. conceived/were responsible for the research. A.M., H.R., F.U.R. and K.S.S. reviewed and revised the paper. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** All the data and models employed and/or generated during the study appear in the submitted article.

**Acknowledgments:** The authors wish to acknowledge Qiupeng Yuan for his moral support and help.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

## *Article* **Research on** *b* **Value Estimation Based on Apparent Amplitude-Frequency Distribution in Rock Acoustic Emission Tests**

**Daolong Chen 1,2, Changgen Xia 1, Huini Liu 1, Xiling Liu 1,2,\* and Kun Du 1**

	- Xuzhou 221116, China **\*** Correspondence<EMAIL_ADDRESS>
**Abstract:** The rock acoustic emission (AE) technique has often been used to study rock destruction properties and has also been considered an important measure for simulating earthquake foreshock sequences. Among them, the AE *b* value is an essential parameter for the size distribution characteristics and probabilistic hazard analysis of rock fractures. Variations in *b* values obtained in rock AE tests and earthquakes are often compared to establish analogies in the damage process and precursory analysis. Nevertheless, because the amplitudes measured on the sample boundary by an acoustic sensor (apparent amplitude) are often used to estimate the *b* value, which cannot descript the source size distribution, it is necessary to develop a method to obtain the size distribution characteristics of the real source from the apparent amplitude in doubly truncated distribution. In this study, we obtain AE apparent amplitudes by applying an attenuation operator to source amplitudes generated by a computer with an underlying exponential distribution and then use these simulated apparent amplitudes to perform a comparative analysis of various *b* value estimation methods that are used in earthquakes and propose an optimal *b* value estimation procedure for rock AE tests through apparent amplitudes. To further verify the reliability of the newly proposed procedure, a *b* value characteristics analysis was carried out on a non-explosive expansion agent rock AE test and transparent refractive index experiment with red sandstone, marble, granite, and limestone. The results indicate that mineral grains of different sizes and compositions and different types of discontinuities of rock specimens determine the rock fracture characteristics, as well as the *b* value. The dynamic *b* values decreased linearly during the loading process, which confirms that variations in the *b* value also depend on the stress. These results indicate that the newly proposed procedure for estimating the *b* value in rock AE tests based on apparent amplitudes has high reliability.

**Keywords:** rock acoustic emission; apparent amplitude distribution; *b* value; completeness amplitude; bootstrap

**MSC:** 74-05; 86-05

#### **1. Introduction**

The power law size distribution relationship of source energy *E* or seismic moment *M*0 is an intrinsic characteristic of the frequency-size distribution in statistical seismology; it can well record the spatial and temporal distribution of rock fractures from a large number of small-scale ruptures to fewer large-scale ruptures, and it has been widely used in seismic research. Since there is a logarithmic relationship between the local magnitude *M* and the source energy *E* or seismic moment *M*0, we can conclude that the local magnitudefrequency distribution obeys the exponential Gutenberg–Richter (G–R) law [1], which has

**Citation:** Chen, D.; Xia, C.; Liu, H.; Liu, X.; Du, K. Research on *b* Value Estimation Based on Apparent Amplitude-Frequency Distribution in Rock Acoustic Emission Tests. *Mathematics* **2022**, *10*, 3202. https:// doi.org/10.3390/math10173202

Academic Editor: James M. Buick

Received: 29 June 2022 Accepted: 25 July 2022 Published: 5 September 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

also been widely used for probabilistic seismic hazard analysis [2–7]. The G–R law is expressed on a logarithmic scale given by

$$
\log\_{10}(N) = a - bM \tag{1}
$$

where *a* and *b* are constants and *N* is the number of earthquakes that occur in a specific time window with magnitude ≥ *M*. More importantly, parameter *a* reflects the size of the time window of observation; slope *b* is an essential tool in seismotectonic studies and seismicrisk analysis within the same time window in a certain area [5,8], which is often referred to as the *b* value. More recently, the *b* value in the G–R law has also been interpreted as an indicator of the applied shear stress and material heterogeneity [9–14]. Thus, its correct computation represents an important challenge in seismology and rock mechanics [15–18].

The selected discontinuities and missing earthquake events in the magnitude-frequency distribution are the main effects on stable *b* value estimation, which is why the two ends of the magnitude-frequency distribution deviate from the G–R law. For some authors, the right and left end points deviating from the G–R law correspond to the magnitude of completeness *M*c (which is defined as the lowest magnitude at which 100% of the events in a space–time volume are detected [8,19] and the auxiliary magnitude *M*0, respectively. Many different procedures for correctly estimating *M*c and *M*0 have been proposed [8,19–25]. Figure 1 shows that with the value of the assumed *M*c starting from the minimum magnitude in the catalog and increasing gradually, the corresponding *b* value and goodness-of-fit change significantly as *M*c ≤ 2 and tend to be stable as *M*c ≥ 2. In view of this phenomenon, some researchers hope to select a sufficiently large *M*c to estimate the *b* value, but for the statistical value of *M*c and *b*, which will reduce a large number of low-magnitude events and further lead to a decrease in the space–time resolution of variations in *M*c and reliability and robustness of the *b* value estimation [26,27]. Therefore, accurately determining *M*c has become the key to stably estimating the *b* values.

**Figure 1.** The effects of variations of assumed *M*c on *b* value estimation [26]. The FMD and LSR are abbreviations of frequency-magnitude distribution and least-squares regression, respectively.

In rock mechanics, the acoustic emission (AE) technique is often used to study the destruction properties of rocks by recording elastic wave information radiated by crack initiation, propagation, and penetration during rock deformation [28–31]. Additionally, the AE test method is also an important means to simulate earthquake foreshock sequences and study focal mechanisms [14,32–42]. Therefore, the space–time variation characteristics of the *b* value obtained in rock AE deformation tests have been used to simulate earthquake precursor characteristics [43–47]. However, unlike the magnitude used in Equation (1) for *b* value estimation in seismology, AE equipment records the high-frequency elastic wave signal of the small-scale rupture, and the AE amplitude is the apparent amplitude measured by sensors at the sample boundary after attenuation from the seismic source [25]. The corresponding apparent amplitude-frequency distribution does not represent the

source size distribution. Thus, the frequency-size distribution law of the source signals collected by the sensors in the rock AE test will be changed owing to the elastic wave attenuation, and the same deviation will appear at both ends of the amplitude-frequency as the magnitude-frequency distribution [33,37], which would affect the size of the real *b* value. (Here, the analogy with the earthquake is used to define the left and right deviation points of the amplitude-frequency distribution of the rock AE test as the completeness amplitude point *A*c and auxiliary discontinuous amplitude point *A*0, respectively.) Although some researchers have long been concerned about the influence of attenuation on the amplitude-frequency distribution in rock AE tests and have also proposed some corresponding compensation methods to obtain the equivalent AE magnitude with the same significance as the magnitude of the earthquake to analyze rock *b* value characteristics [44,48], the equivalent amplitude distribution still cannot fully represent the real frequency-size distribution of rock cracks. To solve this problem, Liu [25] used a statistical method to prove that the apparent amplitude-frequency distribution retains the source frequency-size distribution characteristics, wherein the key to estimating the real *b* value is to properly truncate the apparent amplitude-frequency distribution. A new *b* value estimation method called the Fisher optimal split and global search algorithm (FGS) was proposed to identify the log-linear segment from the apparent amplitude-frequency distribution of the rock AE test for *b* value estimation. In addition, because AE acquisition equipment is very sensitive to the interference of test conditions, such as environmental noise and current signals, a high threshold value of the signal acquisition is generally set for laboratory rock AE tests. Thus, the completeness amplitude *A*c is usually ignored. Therefore, the completeness of the amplitude data should also be considered when considering how to obtain the source size distribution characteristic parameters using the apparent amplitude-frequency distribution in the estimation of the rock AE *b* value.

Based on the discussion above, firstly, we carried out a synthetic AE simulation test to compare and analyze the applicability of the completeness magnitude estimation methods commonly used in earthquakes for the estimation of rock AE completeness amplitude and proposed an optimal procedure for rock AE *b* value estimation by combining Bootstrap [8] and the FGS method which is used for estimating the characteristic parameters of the source size distribution from the apparent amplitude-frequency distribution. Then, we designed a static dilation rock rupturing AE test to further verify the reliability of the newly proposed optimal procedure of *b* value estimation based on the relationship between the *b* value and rock microscopic composition and stress to provide a reliable and accurate *b* value estimation procedure for laboratory rock AE tests. As a result, this research can provide new insights and methods in the analysis of the precursory characteristics in laboratory rock AE tests and rock mass engineering.

#### **2. Optimal** *b* **Value Estimation Procedure Based on Apparent Amplitude-Frequency Distribution**

In rock AE tests, a high threshold is generally set to remove noise interference, and AE equipment will also define the upper limit amplitude, which will result in a doubly truncated distribution of apparent amplitude frequency, so the completeness amplitude *A*c in the *b* value estimation is usually ignored. In this section, generate AE synthetic data with apparent amplitude and select the estimation method of completeness magnitude *M*c commonly used in seismic research to obtain an optimal algorithm for determining completeness magnitude *A*c, which is a key step for *b* value estimation. Then, combined with the nonparametric statistical Bootstrap method, we compared the obtained optimal algorithm of *A*c with the FGS method and determined an optimal procedure of the *b* value estimated for the apparent amplitude data.

#### *2.1. Synthetic Catalogues of Rock AE Apparent Amplitude*

Because the true underlying completeness amplitude *A*c and *b* values are not known in a laboratory rock AE test, we designed a specific simulation scheme to randomly generate synthetic AE data that can be used to clearly compare five *A*c estimation methods [8,23–25]. According to the previous statistical proof and synthetic data generation method [25,49], we designed to generate data arrays with the same length of source amplitude and amount of attenuation, which are all in decibels with a round-off interval of 1 dB. The apparent amplitude after attenuation can be obtained by the subtraction of two randomly arranged arrays of source amplitude and the amount of attenuation. The specific simulation schemes are as follows: firstly, as the *b* value of most papers is equal to about 1 [50], we generated a source synthetic catalog of *i* = 1, 2, ... , *N* events with amplitudes *Ai* in decibels by randomly sampling an underlying Gutenberg–Richter distribution with *b* = 1.0666 and std = 8 (standard deviation) which varies between limits *Ai min* ≤ *Ai* ≤ *Ai max* and has a probability density function *p*(*Ai*). Then, we generated an attenuation amount catalog *δAi* also in decibels with the same length of *Ai*, which varies between limits *δAi min* ≤ *δAi* ≤ *δAi max* and has a probability density function *p*(*δAi*) that obeys the Poisson distribution (or other forms including normal, exponential, Gamma, and random uniform distributions). Finally, the catalog obtained by the subtraction of randomly arranged *Ai* and *δAi* was used to model the amplitude observed at the sample boundary *Ai obs* = *Ai* − *δAi*. In other words, the interval of the apparent amplitude that still follows the exponential distribution was [*Ai min* − *δAi min*, *Ai max* − *δAi max*]. In this paper, we set *Ai min* = 50 dB, *Ai max* = 109 dB, *δAi min* = 1 dB, and *Ai max* = 10 dB to make the apparent amplitude range of synthetic data close to that in normal rock AE experiment. As a matter of fact, the selection of the ranges for source amplitude and attenuation generation has no effect on the results [49]. In addition, to minimize the effect of data volume on deviation discussion, a data volume of 100,000 was generated. Figure 2A shows the apparent amplitude-frequency distribution of generated data.

#### *2.2. Determination of Optimal Estimation Method for Ac and A0*

In studies on the earthquake sequence, determining the completeness magnitude *M*c is the priority of the seismic sequence analysis. In fact, the core of various seismic *b* value estimation methods is the algorithm for searching for *M*c. As some theories in seismic research are often used in rock AE tests, these *M*c estimation methods can also be applied to estimate the completeness amplitude *A*c of the rock AE apparent amplitude-frequency distribution. Common *M*c estimation methods are as follows:

**Figure 2.** *Cont*.

**Figure 2.** Schematic of five methods for determining the completeness amplitude. The amplitude data are the synthetic AE catalog generated by Section 2.1. (**A**) MAXC. (**B**) GFT, the figure shows the variations of residuals and goodness-of-fit with *A*c, and the horizontal dashed lines indicate 90% and 95% confidence. (**C**) MBS, *b* value, Δ*b*, and δ*b* variations trend with cut-off amplitude, the vertical dashed line is *A*c determined by MBS method. Δ*M* = 0.05, d*M* = 0.25. (**D**) MBASS. (**E**) FGS.

These methods have different algorithms for determining the completeness magnitude *M*c. As shown in Figure 2A, the MAXC approach identifies the magnitude corresponding to the maximum curvature of the cumulative magnitude-frequency distribution as *M*c. In fact, this point is also the magnitude corresponding to the maximum frequency of the incremental magnitude-frequency distribution. The GFT method in Figure 2B determines *M*c by comparing the goodness-of-fit *R* of the fitted frequency-magnitude distribution with that of the actual magnitude-frequency distribution, and the *M*c corresponding to *R* ≥ 95% confidence is taken as the completeness magnitude (*R* ≥ 90% confidence can also be accepted as a completeness magnitude in the actual complex earthquake catalog). The MBS method selects the starting point of magnitude where the change in the *b* value tends to be stable as *M*c (Figure 2C), that is:

$$
\Delta b = |b\_{\text{ave}} - b| \le \delta b \tag{2}
$$

where *bave* is the average estimated *b* value from each magnitude within the magnitude interval [*M*co, *M*co + *d*M] and *δb* is the uncertainty of the *b* value proposed by Shi [51]. *bave* and *δb* can be obtained from Equations (3) and (4):

$$b\_{\rm avc} = \sum\_{M\_{\rm co}}^{M\_{\rm co} + dM} b(M\_{\rm co}) \bullet \frac{\Delta M}{dM} \tag{3}$$

$$\delta b = 2.3b^2 \sqrt{\frac{\sum\_{i=1}^{N} \left(M\_i - \langle M \rangle \right)^2}{N(N-1)}} \tag{4}$$

where *b* is the *b* value of the current magnitude *M*co, *N* is the number of events, *M*i is the magnitude corresponding to event *N*i, Δ*M* is the bin width, and <*M*> is the average magnitude of all events greater than *M*co. The FGS method integrates the Fisher optimal split and global search algorithms to determine the log-linear segment of the incremental amplitude-frequency distribution. As shown in Figure 2E, the FGS fully considers the influence of data volume and goodness-of-fit on the estimation results.

Because the *b* value is a statistical parameter, the *b* value estimation will be more stable and accurate with the increase in data volume, and this conclusion has already been discussed by many researchers [8,50,52]. Here, we compared the accuracy of *b* value estimation using synthetic AE amplitude data generated in Section 2.1 to explore the differences in the above five completeness amplitude estimation methods under different data volumes. We randomly generated a source amplitude that obeyed an exponential distribution with a theoretical *b* value of 1.0666 and a data volume of 100,000 as [50, 109] dB and assumed that the attenuation obeyed the Poisson distribution with the interval of [1, 10] dB. According to the theoretical proof of Liu [25], the amplitude interval that still obeys the exponential distribution after attenuation is [49, 99] dB. Therefore, the theoretical completeness amplitude *A*c of the synthetic apparent amplitude data is 2.45. Then, Bootstrap was used to extract 1000 samples from apparent amplitude data with data volumes of 100, 200, 300, 500, 800, 1000, 2000, 3000, 5000, 8000, 10,000, 30,000, 50,000, and 100,000. Figure 3 shows the variation in average *A*c and average *b* value of 1000 Bootstrap samples with different data volumes; the error bar in the figure is the Bootstrap 95% confidence limit.

As shown in Figure 3, the accuracy of *A*c estimated by the five methods was positively correlated with the accuracy of the *b* value, which indicates that selecting an appropriate *A*c is very important for the correct *b* value estimation. This also proves the availability of generated synthetic data to simulate real AE data to a certain degree. It can be seen from the stability of the estimation results that the goodness-of-fit of all methods exceeded 0.9 when the data volume was 3000. Because GFT-90% uses the cumulative amplitudefrequency distribution to estimate *A*c, the goodness-of-fit can reach 0.9 just for a data volume of 100. An interesting phenomenon is that the accuracy of the *A*c and *b* values decreased with an increase in the data volume, and only when the data volume was greater than 300, the results of 1000 Bootstrap samples could be successfully searched. MBASS is more dependent on the amount of data and requires at least 5000 data volumes to successfully search the results of 1000 Bootstrap samples. Therefore, these two methods are not suitable for estimating the *A*c and *b* values. The confidence limits and uncertainty of MAXC and MBS were consistent when the data volume was less than 3000. However, when the data volume was greater than 3000, the confidence limits of the MAXC and FGS methods began to gradually decrease to 0, while the mean *A*c and *b* value estimation of MBS exceeded the theoretical value, and the uncertainty of *A*c was evident. In addition, because both *A*c and *b* are statistical values, the stability of the statistical results largely depends on the amount of data. To better search the log-linear segment, the doubly truncated amplitude-frequency distribution was accepted by FGS, which fully considers the amplitude distribution characteristics. Therefore, the uncertainty of the Bootstrap confidence limit was higher than that of other methods when the data volume was small. However, it can also be seen that the mean *A*c and *b* values of the Bootstrap samples obtained by FGS were still accurate.

**Figure 3.** The variation of *A*c and *b* values estimated by the five methods with different synthetic data volumes. The theoretical *A*c value is 2.45, theoretical *b* value is 1.0666, and the data volume is 100,000. Estimated *A*c and *b* values are the mean value of 1000 Bootstrap samples. The error bar is the 95% Bootstrap confidence limit. (**a**–**e**) MAXC, GFT-90%, MBS, MBASS, and FGS, respectively.

The column distribution of *A*c obtained by 1000 Bootstrap samples of various completeness amplitude estimation methods is shown in Figure 4. This shows that MAXC was the most stable and reliable, followed by the FGS method. Although the estimation result of the GFT method also seemed very stable, it clearly underestimated the theoretical *A*c, which is inconsistent with the actual situation. Similarly, MBASS not only had unstable results but also overestimated *A*c by more than 15% of the samples. The MBS also overestimated *A*c and reached more than 10% of the samples when the data volume was 10,000 and 50,000. From the above analysis, it can be seen that the estimation methods of completeness magnitude in seismology are also applicable for completeness amplitude *A*c

in small-scale rock AE tests, despite some differences in searching ability, accuracy, and stability of estimation results among different methods.

**Figure 4.** Estimation *A*c of 1000 Bootstrap samples of five methods with three data volume at (**a**) 10,000, (**b**) 50,000, and (**c**) 100,000.

As can be seen from the amplitude-frequency distribution in Figure 2, the rock AE data had an auxiliary discontinuous amplitude point *A*0 at the end of the high amplitude segment, which was similar to *M*0 in the earthquake sequence. Therefore, when estimating the *b* value based on the G–R law, the high amplitude data greater than *A*0 would inevitably affect the accuracy of the *b* value estimation. As shown in Figure 2D,E, both the MBASS and FGS can estimate the auxiliary discontinuous point of the magnitude-frequency distribution; therefore, which one is more suitable for estimating the auxiliary discontinuous point is worth further discussing. Figure 5 clearly shows that *A*0 estimated by MBASS was smaller than *A*c. Therefore, we continued to estimate the third discontinuous point. However, only a few of the third discontinuous points met the requirements at the amplitude of 90 dB. By contrast, the FGS method was able to successfully find both end discontinuities, and the estimation *A*0 was very close to the theoretical *A*0 at the amplitude of 99 dB. Compared with Figure 5a,b, the log-linear segment of the apparent amplitude-frequency distribution between *Ac* and *A0* identified by the FGS method was much more stable and reasonable than that of the MBASS method. Therefore, the FGS was more suitable for estimating the auxiliary discontinuous point *A0* at the right end of the apparent amplitude-frequency distribution in the rock AE.

**Figure 5.** Search for *A*c and *A*0 of apparent amplitude-frequency distribution by FGS and MBASS. In order for MBASS to estimate 1000 Bootstrap samples, we generated 100,000 data volumes. (**a**) FGS. (**b**) MBASS; we also counted the number of third discontinuities of 1000 Bootstrap samples.

#### *2.3. The Optimal b Value Estimation Procedure through Apparent Amplitude-Frequency Distribution for Rock AE Tests*

Because the Bootstrap approach can obtain a reliable estimation and avoid outliers, here, combined with the Bootstrap, we proposed to use MAXC to estimate the completeness amplitude *Ac* and FGS to estimate the auxiliary discontinuous point *A*0, respectively. However, after determining *Ac* and *A*0, another issue worth discussing was which regression method we needed to use to estimate the *b* value of the apparent amplitude-frequency distribution between *Ac* and *A*0.

Most papers still use least-squares regression (LSR), which assumes that the frequency data errors are Gaussian, to estimate the *b* value. However, frequencies based on count data have Poisson sampling uncertainties, which cause bias when using LSR for *b* value estimation. Thus, a generalized linear model (GLM) subject to Poisson error can provide a more accurate fit of count data [45,49,53]. Here, we also used the data generation in Section 2.1 to compare LSR assuming a Gaussian error and GLM assuming a Poisson error with apparent amplitude data between *Ac* and *A*0 at a 95% confidence limit. As shown in Figure 6, the confidence intervals of LSR and GLM indicate significantly different changing trends. The confidence intervals of GLM gradually narrow as the amplitude decreases, and the whole amplitude show a "trumpet" shape. However, the confidence interval of LSR remains parallel from small to large amplitudes, which evidently does not conform to the characteristics of the amplitude-frequency distribution. Maximum likelihood estimation (MLE) is currently one of the most popular methods for estimating the *b* value [54], and it can be seen from Table 1 that the standard deviation and bias of the *b* value obtained by this method are also very small. Therefore, once the apparent amplitude-frequency distribution between *A*c and *A*0 is determined, the MLE is also a good choice for regression [49]. However, it is worth noting that the MLE does not have an accurate confidence interval of the doubly truncated amplitude-frequency distribution and cannot further analyze the uncertainty of the estimation results. Thus, we chose GLM regression to fit log-linear amplitude of apparent amplitude-frequency between the left end point *A*c identified by MAXC and the right end point *A*0 identified by FGS, and here we named this *b* value estimation procedure as MFBG (MAXC-FGS-Bootstrap-GLM).

**Figure 6.** Ninety-five percent confidence limits of log-linear segment identified by MFBG from synthetic data with 100,000 data volumes. (**a**) GLM regression assuming Poisson error. (**b**) LSR regression assuming Gaussian error.

#### **3. Application of MFBG in** *b* **Value Estimation of Rock AE Test**

To evaluate the performance of the new *b* value estimation procedure MFBG, a nonexplosive fracturing agent expansion was designed to conduct an AE test on red sandstone, marble, granite, and limestone. The rock broke and formed a specific failure surface because the non-explosive fracturing agent was injected into the three pre-drilled holes in the middle of the specimen (the mass ratio of expansion agent to water is 5:1.7). This experimental design ensured that the AE signals collected by the sensor were all generated by rock expansion fractures and did not rely on location technology to identify valid rupturing data, and a PCI-2 AE system used in this test with six sensors that were tightly attached to the two sides of rock specimen which parallel to the failure surface to collect the AE signals radiated during rock failure, and the parameter settings of AE equipment are shown in Table 2. The specific experimental process, rock sample size, and sensor distribution are shown in Figure 7.

Six sensors were used in the test: to completely remove the signals of non-rock fractures and ensure sufficient data volume, we set the AE signals arriving at the same time every four channels as the rock fracture signal. In this way, we obtained 15 channel combinations, each of which had 10 *b* values, and used MFBG to estimate the mean value of the *b* values for each channel. Table 3 shows the mean *b* values of 1000 Bootstrap samples from the six sensors for red sandstone, marble, granite, and limestone. We can see that the *b* values of the four types of rocks decrease successively because the four types of rocks have different scales of mineral particles, mineral composition, and discontinuity, which lead to different failure scales under the expansion force.

**Table 2.** Parameter settings of AE device.

**Figure 7.** Schematic of non-explosive fracturing agent expansion rock AE test. Non-explosive fracture agent was injected into three boreholes of the specimen, and the fracture surface was formed in the center of the three boreholes. The label of specimen surface a' is parallel to a, and b' is parallel to b [25].

**Table 3.** *b* values of 6 channels of red sandstone, marble, granite, and limestone.

Generally, red sandstone with smaller mineral particles binding tightly is beneficial to the stress increase and the final large-scale fracture formed, but it also limits the initiation and propagation of rock cracks, resulting in a larger *b* value, as the number of large-scale fractures is far less than that of small-scale fractures. In contrast, marble and granite with larger mineral particles and more defects will have more complex heterogeneity and internal structure characteristics, which is unbeneficial to the stress increase and the final large-scale fracture formation, but this also provides the opportunity for crack propagation and penetration, resulting in the generation of more large-scale fractures with smaller *b* values during the rock failure. In particular, the mineral particles composition of limestone is also small and binds tightly, but unlike red sandstone, there is usually a large range of joints in limestone, which largely control the scale of rock failure, resulting in the smallest *b* value than other rocks. Altogether, the estimated *b* values of various types of rock samples are clearly different, which shows that the *b* value depends on the material heterogeneity [55–57], which is also the basic idea for verifying the effectiveness of MFBG.

To explore in more detail the reasons why the *b* values of different types of rocks showed an interval distribution, we use cross polar light technology to carry out a transparent refractive index experiment on rock slices with a thickness of 0.03 mm and the size of 1.7 mm × 1.3 mm to observe rock microstructure and further analyze the relationships between AE *b* value characteristics and rock microscopic composition. As shown in Figure 8, red sandstone is mainly composed of fine-grained quartz, which is tighter than marble and composed of larger dolomite and calcite. Therefore, the *b* value of red sandstone was slightly larger than that of marble. Granite has various mineral particles and defects or voids, which make it prone to large-scale fractures. Limestone is mainly composed of calcite and is even tighter than red sandstone. However, owing to the numerous joints created during deposition, more large-scale fractures are generated. Therefore, the *b* value of red sandstone was the largest, followed by marble and granite, and the *b* value of limestone was the smallest. The results of the four types of rock specimens with *b* values from microstructural characteristics are the same as the estimated *b* value using MFBG, which indicates that this method is accurate and stable for the estimation of the rock AE *b* value.

**Figure 8.** Microstructure of rocks obtained from transparent refractive index experiment. (**a**–**d**) Red sandstone, marble, granite, and limestone, respectively. Qt—quartz; Dol—dolomite; Cal—calcite; Kfs—potash feldspar; Bt—biotite; Pl—plagioclase [58].

The temporal variation characteristics of the *b* value are often used for seismic hazard analysis, crack scale failure description, and damage accumulation assessment in rock; normally, the *b* value is negatively correlated with the stress. As shown in Figure 9, we used the new *b* value estimation procedure recommended in Section 2.3 to estimate the AE *b* value with temporal variation and combined energy, amplitude, and strain to further verify the effectiveness of the MFBG through the internal relationship between the *b* value and crack scale development during rock failure. Because the *b* value is a statistical value, the accuracy of its estimation results is greatly affected by the data volume. Therefore, in order to improve the representativeness and the readability of the analysis results, here we only conduct a special analysis on the granite with the largest data volume collected. Figure 9c shows that the rock failure under expansion stress was manifested as a continuous increase in deformation macroscopically, which has experienced the entire process of compaction of existing defects and microcracks, initiation, propagation, and interpenetration of new cracks, and finally, the formation of the main fracture [59,60]. In this process, the appearance of energy and amplitude signals was usually triggered by a rupturing scale, which in turn led to a decrease or increase in the *b* value under the constraints of statistical laws. Furthermore, Figure 9 shows that the temporal variation characteristics of *b* values decreased continuously, and especially for the time before 16,550 s and after 16,600 s that *b* value almost linearly decreased. This is consistent with Scholz's

laboratory experiments [14], which also show that the *b* value in the size distribution of AE events decreases linearly with differential stress. Therefore, the MFBG method estimated the *b* value in the AE tests can accurately describe the size distribution characteristic of rock failure.

**Figure 9.** *Cont*.

**Figure 9.** Temporal variation of AE *b* value, energy, amplitude, and strain of (**a**) red sandstone, (**b**) marble, (**c**) granite, and (**d**) limestone. Since the development of non-explosive expansion agent was very slow and few signals were generated in the early stage of the tests, the starting time of the experimental analysis selected here was different.

#### **4. Conclusions**

To research the accurate *b* value estimation procedure based on apparent amplitude distribution of rock AE, this paper performs a comparative analysis of various *b* value estimation methods that are used in earthquakes by simulated apparent amplitudes. A new *b* value estimation procedure was proposed, and its reliability is also further verified on rock AE testing and a transparent refractive index experiment. The following conclusions are drawn from this study:

Attenuation causes the two ends of the apparent amplitude-frequency distribution to deviate from the log-linear relationship. Of course, there also retains a finite amplitude interval that still obeys the G–R law, with the *b* value being the same as the source. Therefore, we must search for the log-linear segment in the apparent amplitude-frequency distribution, which can represent the source size distribution characteristics. Furthermore, we generally set the acquisition threshold and maximum output amplitude of the AE equipment, resulting in a doubly truncated exponential distribution for the rock AE frequency-size distribution. We also define the left and right end amplitude points deviating from the G–R law as the completeness amplitude *A*c and auxiliary discontinuous amplitude *A*0, which correspond to the completeness magnitude *M*c and auxiliary discontinuous magnitude *M*0 in the earthquake catalog. Additionally, more attention must be paid to the completeness amplitude *A*c and auxiliary discontinuous amplitude *A*0 to obtain the source size distribution through the *b* value of the apparent amplitude frequency.

The estimation method of completeness magnitude commonly used in earthquakes was also suitable for identifying the completeness amplitude in rock AE. Especially, the MAXC can better estimate *A*c, and the FGS can better estimate *A*0. Moreover, we combined the Bootstrap approach to propose a new *b* value estimation procedure named MFBG, which can better fit the apparent amplitude distribution without any attenuation compensation, and the effectiveness of the new method was verified by the relationships between rock crack size distribution and mineral grains and internal structure characteristics under laboratory rock AE tests and transparent refractive index experiments. This study can provide new insights and methods for studying the precursory characteristics of laboratory rock tests and rock mass engineering through the variation of rock AE *b* value.

**Author Contributions:** Conceptualization: X.L. and C.X.; methodology: D.C.; investigation: D.C. and H.L.; software: C.X. and K.D.; writing—original draft: D.C. and X.L.; writing—review and editing: X.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the National Natural Science Foundation of China (Grant 42172316) and Natural Science Foundation of Hunan Province (2021JJ30810), and the Research Fund of The State Key Laboratory of Coal Resources and Safe Mining, CUMT (SKLCRSM21KF005). The authors are grateful for the financial support provided by the funds.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data associated with this research are available and can be obtained by contacting the corresponding author.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

## *Article* **Pollutant Migration Pattern during Open-Pit Rock Blasting Based on Digital Image Analysis Technology**

**Jiangjiang Yin 1, Jianyou Lu 1,\*, Fuchao Tian 2 and Shaofeng Wang 1,\***

**\*** Correspondence<EMAIL_ADDRESS>(J.L<EMAIL_ADDRESS>(S.W.)

**Abstract:** Previous studies have revealed that toxic gases and dust (smoke dust) are the most common pollutants generated by the blasting operations in open-pit mines, which might lead to a threat to the environment's condition, health and safety, and properties protection around the blasting site. In order to deal with the problems, a pollution evaluation system was established based on the fractal dimension theory (*Dbox(P)*) and grayscale average algorithm (*Ga*) in digital image-processing technology to recognize and analyze the distributions of the smoke-dust cloud, and subsequently determine the pollution degrees. The computation processes of *Dbox(P)* and *Ga* indicate three fitted correlations between the parameters and diffusion time of smoke dust. Then, a pollution index (*Pi*) is put forward to integrate the global and local features of *Dbox(P)* and *Ga*, and develop a hazard classification mechanism for the blasting pollutants. Results obviously denote three diffusion stages of the pollutants, mainly including generation stage, cloud-formation stage, and diffusion stage. In addition, it has been validated that the proposed system can also be utilized in single-point areas within a whole digital image. Besides, there are variation trends of the thresholds T1 and T2 in binarization with the diffusion of pollutants. With this identification and evaluation system, the pollution condition of smoke dust can be obviously determined and analyzed.

**Keywords:** rock-blasting; smoke-dust pollution; digital image processing; fractal dimension; grayscale average; hazards classification

**MSC:** 54H30

#### **1. Introduction**

With the evolution in industrial fields, the mining industry is confronted with a huge energy demand, and is gradually emerging as the critical support for technological development and progress [1–3]. To adapt to community challenges, ore-deposit extractions are constantly expanding, especially for open-pit mines [4,5]. Blasting is a conventional excavation measure for the ores crushing on the ground, and also an option for constructing roads, tunnels and other constructions in mine projects [6]. However, during blasting operations, there are undesirable consequences engendered, such as seismic waves, explosive flying scatter (fly rocks), air shockwave, noise, toxic gases, blasting dust, etc. [6,7].

Within these hazards, blasting dust, the PM2.5 and PM10, and some toxic gases, including CO, NOx, SO2, H2S, and so on, become prominent components of environmental pollutants in blasting sites [8]. Most of these toxic gases are released in a range of concentrations due to environmental and technical factors [9]. Inhaling a certain amount of these toxic gases during their diffusion period can be fatal [10]. Furthermore, prior studies have denoted that it can reach thousands of milligrams per cubic meter for the local blasting dust concentration in surface mining [11]. The combination of blasting dust with toxic gases, which is collectively referred to as smoke dust, has attracted extensive attention owing to its salient hazardousness on air quality in mining sites and detrimental effects on human

**Citation:** Yin, J.; Lu, J.; Tian, F.; Wang, S. Pollutant Migration Pattern during Open-Pit Rock Blasting Based on Digital Image Analysis Technology. *Mathematics* **2022**, *10*, 3205. https:// doi.org/10.3390/math10173205

Academic Editor: Radu Tudor Ionescu

Received: 21 July 2022 Accepted: 1 September 2022 Published: 5 September 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

health [12,13]. It has been demonstrated that a fairly high percentage of miners suffer from respiratory diseases or reduction in respiratory-system immunity due to long-term exposures to high-concentration dust environments, especially pneumoconiosis [14,15]. In addition, the accumulation of dust reduces the humidity of the atmosphere, which then increases the ambient temperature [16]. Meanwhile, it is noteworthy that smoke dust reduces visibility in the mining site, resulting in high working risks and a significant reduction in production efficiency [17]. Large amounts of dust pollution can depose on the ground, causing vegetation degradation and animal habitat destruction [18]. With the characteristic of high concentration and wide range of diffusion it might spread into the water area for dust clusters, thereby contaminating the downstream areas [8]. Studies over the past decades have provided important information on the generation and diffusion mechanisms of smoke dust in open-pit mines. Numerical simulations, for instance, Fluent software, LS-DYNA etc., were conducted to figure out the smoke dust movement, diffusion rate, and the time consuming to float out of the pit, validating the effects of vortex action and the vertical height in the pit on coal-mine dust, and optimizing blasting parameters at the same time [19,20]. Inspired by the approaches, Tang et al. [21] established a circulating accumulation emission model to analyze dust diffusion, including dust emission, retention, and diffusion, respectively. Besides, pertinent artificial intelligence systems, such as particle swarm optimization (PSO) and the long short-term memory (LSTM), were employed to predict and control dust concentration [22–24]. In terms of the evaluation system of efficiency of dust reduction with multiple techniques, there are also studies [25,26] in virtue of some tools (e.g., computational fluid dynamics model (CFD): ANSYS/Fluent 10.0 software, ANSYS, Inc., Pittsburgh, Pennsylvania, USA) mainly aiming at monitoring the dust diffusion, and then guaranteeing the mitigation of the smoke-dust pollution.

Furthermore, there are increasing concerns that, exposed to high-concentration smoke dust in open-pit mines, measurement equipment has been a crucial factor for immediately monitoring the pollution condition. Attalla et al. [9] adopted non-dispersion infrared (NDIR) and differential optical absorption spectroscopy (DOAS) to examine the diffusion of fugitive NOx and other pollutant gases. In addition, dust sample sensors, which are commonly based on ultrasound, optical, and electrochemical sensing elements, have also emerged as conventional methods for assessing contaminant gases and particulate fractions of dust emission [2]. With the widespread application of unmanned aerial vehicles (UAVs), researches were performed to load a series of sensors in the rotary-winged platforms, such as CO2, NO2, SO2, PM10, etc. [27]. When in the application process, it is essential for the sampling collection location and source scales to be considered owing to the differences in fixed-wing and multi-rotor UAVs, and thus the long response time, limited measurement distance, and the accuracy of sensing elements hindered further use of UAVs.

Recent developments in the field of machine vision have led to a renewed interest in various industries due to their merits in non-contact measurements, full-visualization monitoring, and affordable cost [28–31]. In spite of the fact that the related studies on these vision-based models are currently at an early stage, digital image technology, mainly based on grayscale features and morphological features, has been widely applied in pollution analysis and control [32–35]. Dust-pollution measurement models established on the basis of grayscale features attempted to develop a relationship between the dust concentration and grayscale features extracted from dust images. Yue et al. [30] presented a photobased method for analyzing air pollution using grayscale data. Thereafter, a logarithmic relation was acquired between the grayscale values and the powder volume by Grasa et al. [34], providing a calibration curve without any fitting constants. On the other hand, morphological feature models are inherently associated with the connections between the morphological features (such as the particle shape, size distribution, and surface structure) and dust concentration. R. Davies [36] and M. Taylor [37] were convinced that parameters of concentration and morphology were important in particle detection, and the latter presented mathematical models to characterize the particles. M. Paphael et al. [38] estimated the mean particle size utilizing optical parameters, on-line turbidity, and solid concentrations. It is also notable for the fractal dimension, another aspect of morphological features, to describe the space occupancy and the irregularities of granulated materials [39–42]. To figure out the effects of fractal dimension on the porous texture of coal dust images, Mahamud et al. [41] conducted experiments to measure particles exposed to different temperatures. Zhang et al. [42] apply fractal dimension to record the particle size of crushed specimens and obtain a legible regularity of size distribution in different stress conditions. Nevertheless, traditional grayscale or morphological models only considere the local features of pollutants but ignore the global features of their distributions, resulting in a one-sided evaluation of the pollution environment.

To further present the main advantages and limitations of the technology and equipment of previous studies on the pollutant migration that happened during rock-blasting operations, an overall comparation between them is summarized as the Table 1:

**Table 1.** Comparation of previous studies on analysis of smoke dust.

Overall, all the evidence reviewed here seems to require a comprehensive evaluation system and distribution analysis for smoke dust pollution caused by blasting operations in open-pit mines.

In this paper, a vision-based model was established based on the fractal dimension theory (*Dbox(P)*) and a grayscale average algorithm (*Ga*) for smoke dust evaluation in openpit mines. Framed pictures were first collected from the video recorded by a UAV, and then, through a series of image processing operations, *Dbox(P)* and *Ga* were separately calculated, the results of which were fitted by logarithmic and linear functions, respectively. For the sake of considering both global and local characteristics of pollutants, and fully adopting simplified expressions to characterize the various generation stages of smoke dust in blasting sites, a pollution index (*Pi*) was proposed to integrate the *Dbox(P)* and *Ga*, developing a hazard classification mechanism corresponding to four-level regions. With this evaluation system, the pollution condition of smoke dust can be obviously developed and analyzed, and eventually adequate preparations in response time and emergency preparedness are capable of eliminating detrimental effects on health protection and property safety, as well as improving the environmental quality around the blasting sites.

#### **2. Image Processing Method Based on Grayscale Average and Fractal Dimension**

Digital image processing has been successfully applied in the engineering field depending on the difference in the morphological characteristics or the color features of the research object against the homologous background [44–46]. Based on the technical principles, and the available devices, it is possible to achieve the monitoring and evaluation of the smoke-dust. The specific procedures of the algorithm in this research can be illustrated as a flowchart in Figure 1:

**Figure 1.** Flowchart of the image processing algorithms for extracting features of pollutants, and establishing an evaluation mechanism for hazards classification of smoke dust during blasting operations.

As shown in Figure 1, the proposed method can be divided into two parts, including the modules of image-processing operations and properties analysis, which primarily aim to acquire the feature image of the pollutants, and relevant calculation results about their distribution properties. According to the flowchart, pollution conditions in blasting sites can be analyzed and evaluated through a series of digital images. The specific application procedures are briefly presented below:

Step one: collect pollutants images of smoke dust at the blasting site, and then the software reads the images to form a three-dimensional digital matrix Cm×n×3, which includes pixel values *cij*k and homologous locations (*i*, *j*); m, n represent the image size in height and width, and the third dimension of the matrix represents three color channels: R, G and B, respectively.

Step two: based on the image matrix Cm×n×3, extract the main regions where pollutants distribute via Equations (1) and (2) below in accordance with the color characteristics of the smoke dust [47]. As is usual, the smoke dust pollutants in non-coal mines appear a brighter color, and thus, compared to several dark backgrounds, the pixel values of pollutants tend to be higher. Through the screening method, a series of pixels representing superfluous regions will be removed, whereby the images can be prepared for subsequent image processing operations.

$$C\_{sum} = c\_{i\bar{j}1} + c\_{i\bar{j}2} + c\_{i\bar{j}3} \tag{1}$$

$$
\mathbb{C}\_{sum} > a \tag{2}
$$

Step three: fulfill further image-processing operations with pertinent programs. Image graying is enforced using the empirical formula shown in Figure 1, the matrix result of which is denoted as *Gm*×*n*, and by virtue of image graying, the color dimensions are converted into two dimensions, viz., the image size. Then, binarization was served to recognize the regions of interest, in which a double-threshold method (T1, T2) was put forward to improve the performance of the feature extractions, and in this way, pixels within the gray image that are, respectively, lower or higher than T1 and T2 are assigned as 0 (black in color), while the intermediate areas between the thresholds are assigned as 1 (white in color). Noise reduction is then adopted to eliminate excrescent noise points caused by the forming transversion of pollutant images [48].

Step four: calculate the *Dbox(P)* and *Ga* of the smoke-dust images to quantitively analyze the distribution changes of on-site pollutants. In the process, the two-dimensional box counting method was developed to acquire the fractal property of the pollutants [40]. Within the fractal calculation processes, the binary images are traversed by a slew of windows, the side lengths of which are decided by the original picture size, and are incremental pixel by pixel until the width and height reach half of the original length. Then, for each window, keep a track of the number of square boxes *N(ε)* intersecting the pollutant parts as well as the lengths of the boxes *ε*. Afterwards, instead of the mathematical definition of fractal dimension theory as Equation (3), a regression method is used to analyze two groups of the box data through logarithmic processing, which can be expressed as Equation (4). Obviously, the slope of the regressed equation *Dbox(P)* is the fractal dimension of the smoke dust.

$$D\_{\text{box}}(P) = \lim\_{\in \to 0} \frac{\log N(\in)}{\log 1/\in} \tag{3}$$

$$\log N(\in) = D\_{\text{box}}(P) \* \log 1/\in + \mathcal{C} \text{ ( $\mathcal{C}$  is a regressed constant)} \tag{4}$$

Meanwhile, *Ga* of the binary images is calculated based on Equation (5), which can be defined as the division operation between the pixels denoting pollutants and site background, respectively. When in the calculation process, a grayscale distribution curve is plotted to develop the pollutants distribution along the directions of horizontal strike.

$$\mathbf{G}\_{a} = \frac{\sum\_{i=1}^{h} \sum\_{j=1}^{w} b\_{i,j}}{h \times w} \tag{5}$$

where *h* and *w* represent the height and width of smoke-dust images.

Step five: through the computation of the *Dbox(P)* and *Ga* of the pollutant, a pollution index *Pi* is provided to achieve a comprehensive consideration of *Dbox(P)* and *Ga*, and the integration expression can be depicted as Equation (6):

$$Pi = D\_{\text{box}}(P) + G\_a \tag{6}$$

Then, according to the change tendency of *Pi*, it can be more accurate to evaluate the pollution condition from a quantitative aspect, and with a series of pollutant images collected from a whole blasting project, the index *Pi* is capable of dividing into four regions of pollution levels, Level I~IV, viz. Therefore, the pollutant condition of smoke dust can be dramatically characterized no matter the period of monitoring or throughout the whole working process by this research.

#### **3. Algorithm Application of** *Dbox(P)* **and** *Ga* **with a Fugitive Dust Picture**

For the sake of further demonstrating the proposed evaluation method of pollutant distribution in a blasting site, an image sample with fugitive dust shown as Figure 2 was applied to validate the effectiveness of the image-processing algorithm, and measurement of *Dbox(P)* with *Ga*. Firstly, the original color image of fugitive dust as Figure 2a was read, where in the image dust particles emerge distinctly with homologous pixels when compared to the black background. With the application of the image-processing method, the region of the floating dust was arguably extracted to produce a binary image, which can be visualized as Figure 2b. Based on the binary image, the calculation processes of *Dbox(P)*, *Ga* are also displayed in Figure 2c,d, respectively. As for the grayscale distribution, an ergodic window of 438 × 10 pixels in size was employed to cover the binary image along the direction of the horizontal, and with the traversal calculations, the grayscale distribution, cumulative distribution curve, and grayscale average were obtained, the results of which are displayed in Figure 2e. According to the grayscale parameters, *Ga* of the whole image is 0.2490, and the dust objects primarily mainly distribute at 220~550 pixels of the horizontal strike. Similarly, the calculation procedures of *Dbox(P)* are devised as Figure 2d using the two-dimensional box counting method. With a series of measuring boxes covering the binary image, the computation procedures are revealed as the fitted equation in Figure 2f, and the *Dbox(P)* of the dust image is 1.977, the *R*2 of which is 0.9998, embodying a precise fitting effect when applying it to research objects of fine grains.

**Figure 2.** Description of the proposed method for the *Dbox(P)* and *Ga* of the fugitive dust picture, including (**a**) original image of fugitive dust, (**b**) binary image through a series of image-processing operations, (**c**,**d**) calculation process of grayscale properties and *Dbox(P)* and (**e**,**f**) results of grayscale distribution, *Dbox(P)* and *Ga* of the fugitive dust.

In order to verify the accuracy of calculation results in *Dbox(P)* and *Ga* with the proposed binary dust image handled from the original image, the grayscale average distribution of the gray image, remaining quite complete image information with a 256-level grayscale form, was devised to analyze the change tendency of dust pixels as a comparison to the binary image, and quantitatively appraise the effectiveness on object segmentation within the application of image processing. The analysis results can be shown as Figure 3.

As is mentioned with the grayscale parameters, dust particles are primarily distributed at 220~550 pixels of the horizontal strike in accordance with the binary image. It can be perceived that the developed trend of the grayscale average from gray image tends to be consistent, performing a fairly smooth waveform shape. In addition, it is almost coincident with the high-distribution region of dust particles in any binary image or gray image. While there is a hysteresis effect for the grayscale parameter of gray image at the boundaries between research objects and background, producing two hysteresis intervals in the rising and falling stages of the grayscale average, the area of the hysteresis intervals accounts for about 5.58% of the total enclosed area by the grayscale distribution line of gray image. It is imperative to note the significance turning points, data turning point G1 (90, 25.16) of gray image embodying a hysteresis effect contrasting to B1 (210, 7.37 × 10−5) of binary image for the first hysteresis interval, whereas G2 (560, 51.10) illustrates an almost the same horizontal position with B2 (560, 3.08 × 10−5). A possible explanation for this might be the low pixel values of edge dust particles in the image (i.e., grayscale value). As a whole, in spite of the existence of hysteresis intervals, the application image of fugitive dust presents a validity in the object recognition of dust clouds and other fine-grained substances from a quantitative aspect.

**Figure 3.** Result comparation on grayscale average between the binary image and grayscale image, including (**a**) original image of fugitive dust, (**b**) grayscale image, (**c**) grayscale average of binary and gray image for algorithm validation.

#### **4. Results on the** *Dbox(P)* **with** *Ga* **and the Pollution Evaluation System of Blasting Operations**

A field experiment for analyzing smoke dust distribution, and evaluating the pollution condition based on the proposed method, was conducted in a deep open-pit mine with multi-step terrain production in Xinjiang, China. Shallow hole blasting was adopted for the ores crushing along the hillside wall, and the whole blasting operations with partial diffusion were recorded by a UAV platform flying directly above the blasting region with a fixed height. The platform can be predetermined in an autonomous direction, and illustrates high-temporal sampled pictures in real time.

A total of 95 images were extracted from the video recorded by the UAV with an interval of 1 s, which are 640 pixels in width and 368 pixels in height, integrally revealing the processes of blasting preparation, exploding and diffusion of smoke dust, respectively, which are shown in Figure 4.

**Figure 4.** Frame pictures of pre-blasting, exploding stage and diffusion stages extracted from the UAV video data.

From the images in Figure 4, it can be found that there is an obvious contrast between smoke dust and the field background throughout the periods of blasting operation. At the initial stage, debris and dirt deposits on the ground were lifted up due to the shock wave of blasting, and with the gravity effect, the smoke dust was then dispersed suffering from the wind floating and other meteorological conditions. At this stage, most distribution regions of the pollutants perform as partially bright in color, while regions representing the ground surface and the bench tend to be dim. For a clearer description, a three-dimensional elevation grayscale image was put forward to show the grayscale distributions and the differences in elevation, as depicted in Figure 5. The elevation range in Figure 5e extracted from the region presents an exceedingly grayscale difference between the pollutants and ground surface. By this method, feature extraction of a complicated environment between the pollutants and background can be fulfilled through a series of image-processing algorithms proposed in this research.

**Figure 5.** Grayscale elevation map of smoke dust with homologous calculation processes, including (**a**) original image of smoke dust, (**b**) gray image, (**c**) image fragment of gray image, (**d**) digital matrix of the image fragment, and (**e**) grayscale elevation map based on the image fragment.

To ulteriorly assess the effectiveness of the proposed method in blasting operations, the 95 images from the blasting work were processed by the image-processing algorithm, and then served to determine the *Dbox(P)* and *Ga*. Within the pollutant images, the computation processes begin from the 16th image because the 1st~15th images belong to the prior stage of blasting, and smoke dust is barely distributed in these images. Partial image-processing results of the 15th, 30th, 45th, 60th, 75th and 90th images are shown as Figure 6:

**Figure 6.** Partial visualization results of the image-processing operations, including original images and binary images, corresponding to the 15th, 30th, 45th, 60th, 75th and 90th images.

Following the image-processing results it can be discerned that at the first stage of blasting, a large amount of smoke and dust are engendered to form smoke-dust clouds, mainly due to the explosive at each charge point, such as the 30th picture. At this stage, dust particles start to float due to the shockwave caused by the explosion, which is much larger than the gravity of the dust particles, friction between dust particles and rock particles, and adhesion of dust against other items. Then, with the completion of the blasting operation, the smoke-dust clouds gradually develop, when in this process the smoke dust continues to move upward as a consequence of the inertia effect, and the motion track of a large amount of smoke dust is superimposed, subsequently generating a mushroom or mushroom-like cloud. In addition, from a qualitative perspective for the pollutant distribution, the amount of smoke dust appears to have an upward trend with the explosion process, and then enters into a volatility phase, (60)~(75), for instance. Eventually, as the smoke dust continues to spread, which mainly interfered with meteorological conditions in this period, especially for the disturbed wind flow, the distribution is confronted with a certain degree of reduction.

To further describe the amount of pollutant, *Dbox(P)* and *Ga* for the 15th~95th smoke dust images were calculated, and the results with corresponding calculating time are demonstrated in Table 2.

**Table 2.** Calculation results of fractal dimension and grayscale average for partial pollutant images.

According to the results of *Ga* and *Dbox(P)*, a relatively legible data trend is presented, which can be summarized as increasing to a peaking value, maintaining a stable fluctuation within a certain range, and then decreasing with time, respectively. Moreover, it is found

that there is a discrepancy in calculating time between *Ga* and *Dbox(P)*, the reasons for which may be explained by the varying levels of program complexity.

For the sake of further analyzing the *Dbox(P)* and *Ga* of the collected frame pictures from the video footage of field blasting, computations of all images by the proposed method were carried out and the specific results with related fitting curves can be denoted in Figure 7:

**Figure 7.** Calculating results of *Ga* and *Dbox(P)* with corresponding fitting curves. *G*1, *G*2, *G*3 and *D*1, *D*2, *D*3 represent fitting curves of grayscale average and fractal dimension, respectively.

It can be easily discovered from Figure 7 that there are characteristic distributions for 80 image samples (16th~95th) on the basis of the results, which accord with a pattern consistent with the qualitative description, mainly including three stages: increasing, stable fluctuations, and decreasing. The fitted equations for the *Ga* and *Dbox(P)* are listed in Table 3.