Abstract

In China, gas and oil reserves are very scarce, but coal resources are abundant in the energy architecture, which decides that coal will remain the dominant energy source for a long time in the future. The accurate prediction of the size and extent of surface movement after coal seam mining is of great significance for the safe promotion of production activities in the mine area and the safety of people’s lives and properties in the mine area. The surface movement deformation under thick loose seam conditions indicates the phenomenon of a large subsidence value and influence range. To predict the size and range of surface movement deformation under thick loose layer conditions accurately, a hyperbolic secant model is constructed based on the hyperbolic secant function. For high nonlinearity of the model parameters, the adaptive step fruit fly algorithm (ASFOA) is introduced into the process of solving the model parameters. Simulation experiments are conducted in three aspects: monitoring point antideficiency, antigross error, and parameter stability. The simulation results show that the ASFOA algorithm achieves high accuracy in finding the parameters of the hyperbolic secant model. The hyperbolic secant model was applied to the 11111 working face under the mining conditions of thick loose layer geology in the Huainan mine. The engineering application results indicate that the hyperbolic secant model performs well on the prediction of surface movement deformation under thick loose layer conditions.

1. Introduction

In China, the current energy structure of more coal, less gas, and poor oil have made coal the dominant energy source for a long time [13]. In the past decades, coal has made great contributions to the rapid development of China’s economy, but the large-scale mining of coal resources has caused a series of environmental problems such as surface subsidence, damage to arable land, water pollution, and air pollution [47]. How to accurately predict the surface movement deformation caused by coal mining is important for guiding coal production activities, rebuilding the ecological environment of mining areas, and ensuring the safe use of mining structures.

The loose layer is generally composed of quaternary as well as neotertiary associated sediments or accumulations. When the thickness of the loose layers exceeds 50 m, they are called thick loose layers [8]. Thick loose layers of different thicknesses exist in northeast, north, central, and east China [9, 10]. Compared with the general geological conditions, the mining subsidence characteristics under the geological conditions of thick loose layers show significant differences. They mainly include large subsidence values, slow convergence of surface subsidence basin boundaries, larger subsidence influence range, and more intense movement deformation during the active surface movement period [1113]. The form of movement and compressive capacity of the loose and rocky layers under the influence of mining are quite different. The rock layers are harder and can play a significant role in supporting internal stresses when they are damaged. The hardness of the loose layer is less, which makes the loosened layer act as a huge load on the rock layer without supported capacity when the internal stress of the loose layer changes [14]. The huge difference between the nature of the rock and the loose layer makes the surface movement deformation show special characteristics such as a larger subsidence value and slow convergence of surface subsidence basin boundaries.

The stress state of the surrounding rock is disrupted, and the stress is redistributed after the coal is extracted [15, 16]. During this process, the ground surface moves and deforms. Mining subsidence prediction is one of the core elements of mining subsidence. The current related research on predicting mining subsidence mainly has three categories. The first category is the empirical method based on actual measurement information. It established an empirical formula suitable for the mine by preprocessing and analyzing a large amount of actual measurement data, such as using the typical curve method and the profile function method [17, 18]. Although it achieved high accuracy, the empirical formula established by this method is constrained within this mining area instead of other areas. The second category is the theoretical model method. It generalizes the overlying rock layers into some type of mathematical or mechanical model. Baryakh et al. formulated the rock layers into a viscoelastic model and constructed a dynamic prediction model of surface movement deformation by considering time as the independent variable [19]. Salamon proposed the surface element theory based on the elasticity theory [20]. Hao et al. construct a surface full-section prediction model based on the elastic plate theory [21]. The model established by this method is more general but has lower accuracy with unclear physical meaning of the parameters. The third category is the influence function method. This method is based on the distribution function and the principle of equal influence superposition to construct the surface movement deformation formula [22, 23]. Liu developed the probability integral method based on the stochastic medium, which is most widely used in China because of its good universality, the clear physical meaning of parameters, and high accuracy [24]. However, the probabilistic integral method for predicting surface movement deformation under special geological mining conditions (such as multiple coal seams, steep coal seams, and thick loose layers) is not applicable [2529].

Due to the peculiarities exhibited by the ground movement deformation under the action of thick loose layers, the existing prediction models are difficult to accurately predict the ground movement deformation under thick loose layer conditions. In this paper, we propose a prediction method based on the hyperbolic secant function suitable for surface movement deformation under the condition of the thick loose layer. The prediction formulas for the movement and deformation of the finite mining main section and any point of the main section are derived. For the high nonlinearity of the prediction equation, the ASFOA algorithm is proposed to be introduced into the parameter acquisition. The applicability of the hyperbolic secant model proposed is evaluated through simulation experiments and engineering applications.

2. Hyperbolic Secant Model Construction and Parameters Acquisition

2.1. Construction of Hyperbolic Secant Model

Hyperbolic functions belong to one class of functions similar to circular functions. Hyperbolic functions are widely used in areas such as hanging chain lines and Laplace equations [30]. The hyperbolic secant function expression is shown in equation (1). The hyperbolic secant function curve is like a semi-infinite mined probability integral function curve. Using the hyperbolic secant function as the influence function of unit mining, the expression of the hyperbolic secant function of unit mining can be obtained as shown in equation (2).where R is the major influence radius.

The plot of the subsidence value of the unit mining versus the parameter R is shown in Figure 1. As the value of R increases, the boundary convergence becomes slower, and the maximum subsidence value becomes smaller.

The hyperbolic secant function curve form and the surface movement deformation curve form have high similarity. The hyperbolic secant function is introduced into the field of mining subsidence to construct a new model to predict surface movement deformation. First, we constructed a unit mining subsidence prediction formula based on a hyperbolic secant function. Then, we derived a semi-infinite mining prediction formula followed by a prediction formula for finite mining strike and inclination in the main section. In the end, a prediction formula for any point in the moving basin was applied. The derivation flow is shown in Figure 2.

2.1.1. Prediction of Surface Subsidence during Unit Mining

In practical production activities, when thick loose layers and bedrock have the same medium, the movement deformation values calculated theoretically will have large deviations from the actual ones. In this paper, two hyperbolic secant subsidence basins with different parameters R are superimposed and combined based on a certain ratio, in order to construct a new model for surface movement deformation prediction. The formula for the projected subsidence of the unit mining after the superposition and combination is shown as follows:where P is the ratio factor, R1 and R2 are the major influence radii of main influence radii of two subsidence basins, and R1 and R2 can be obtained by the following equations: and , where H is the mining depth, is the tangent of the main influence angle corresponding to R1, and is the tangent of the main influence angle corresponding to R2.

From Figure 3, the constructed surface movement deformation prediction model is able to describe the subsidence basin formed by unit mining.

2.1.2. Prediction of Surface Movement Basin Strike Main Section Movement Deformation during Semi-Infinite Mining

The establishment of the surface coordinate system, the coal seam coordinate system, the surface coordinate system, and the coal seam coordinate system are shown in Figure 4. Taking the surface point above the working face boundary as the coordinate origin O, the positive direction of the x-axis points to the side of the mining area, the horizontal movement U(x) upward is positive, and the subsidence W(x) downward is positive. The top plate at the coal wall of the working face is the coal seam coordinate origin O1, the s-axis points to the side of the mining area along with the top plate of the working face, and the vertical coordinate axis Z-axis is vertically upward. The subsidence of any point A on the ground caused by unit mining is We(x − s), and the horizontal movement is Ue(x − s).

From Figure 4, the mining unit with the horizontal axis causes the subsidence of any point A on the surface as We(x − s) and the horizontal movement as Ue(x − s). The formula for subsidence is given in equation (4) and the formula for horizontal movement is given in equation (5).where B is a constant.

In semi-infinite mining, the subsidence value of any point A on the ground is the sum of the subsidence values caused by the mining of each unit in the range s = 0⟶+∞. The subsidence of point A is obtained by the integral transformation as shown in the following formula:

The integral transformation of (5) can be derived from the semi-infinite mining strike of the horizontal movement of the main section U(x), the horizontal movement formula is shown in the following formula:

Simplifying equation (7), let

Simplified as

The tilt i(x) is the first-order derivative of the subsidence W(x) and the tilt equation is given in the following formula:

The curvature K(x) is the first-order derivative of the tilt i(x) and the curvature equation is given in the following formula:

The horizontal deformation ε(x) is the first-order derivative of the horizontal movement U(x) and the horizontal deformation equation is given in the following formula:

2.1.3. Prediction of Surface Movement Basin Strike Main Section Movement Deformation during Limited Mining

The schematic diagram for calculating the moving deformation of the main section of the moving basin strike during limited mining is shown in Figure 5. The actual mining length of the main section of the working face is D3, and the inflection offset distance is S3 and S4 due to the cantilevering effect of the overlying rock layer. The formula for calculating the moving deformation of the main section of the working face during finite-element mining is calculated by superposition, and the formula is shown as follows:where l is the effective mining length of the main section of the strike, and the effective mining length is l = D3 − S3 − S4.

2.1.4. Prediction of Surface Movement Basin Inclination Main Section Movement Deformation during Limited Mining

The principal diagram for calculating the moving deformation of the inclined main section of the mobile basin during limited mining is shown in Figure 6. Due to the cantilevering effect of the overlying rock layer, the inflection point offset distance is S1 and S2 for the inflection point offset distance of the lower mountain and the inflection point offset distance of the upper mountain, respectively. The moving deformation of the inclined main section during limited mining is calculated according to the principle of equal influence, as shown in equations (14) and (15), where before the comma represents the variable and after the comma represents the parameter, i.e., y represents the variable, G1 represents the relevant parameter corresponding to the downhill direction, and G2 represents the relevant parameter corresponding to the uphill direction.

Of whichwhere D1 is the inclination length of the working face, α is the coal seam inclination angle, and θ0 is the mining influence propagation angle.

2.1.5. Prediction of Surface Any Point Movement Deformation during Movement Basin

Any point space coordinate system is shown in Figure 7, O1 is the origin of the coal seam coordinate system, O is the origin of the ground coordinate system, and the horizontal projections of the two coordinate systems overlap. The mining of unit B causes the subsidence of any point A on the ground surface as shown in the following equation:

As shown in Figure 7, if the mining area is the area enclosed by D1, D3, s-axis, and t-axis, then the entire mining causes the subsidence of A as shown in the following equation:

In view of the expected derivation process for the movement deformation of the strike and inclination main sections during finite mining, equation (17) can be transformed as follows:

The tilt i(x, y, φ) of any point A(x, y) along the φ-direction is the directional derivative of the sink W(x, y) in the φ-direction, and the formula is given in the following equation:

Similarly, the curvature calculation formula can be obtained as follows:

The horizontal movement calculation formula is as follows:

The horizontal deformation calculation formula is as follows:

2.2. ASFOA Algorithm for Hyperbolic Secant Model Parameters Acquisition
2.2.1. ASFOA Algorithm

The hyperbolic secant model is a highly nonlinear function. Obtaining the globally optimal model parameters is a difficult research problem. The fruit fly algorithm (FOA) is a global optimization search method based on the characteristics of fruit flies searching for food [3133]. The step size of the FOA algorithm is fixed in the process of searching for the global optimal solution. If the step size is set too large, the global optimal solution can be missed easily. If the step size is set too small, the locally optimal solution can be trapped easily as well. For the FOA algorithm, it is difficult to determine the search step length. The ASFOA algorithm is proposed to improve the search step length of the FOA algorithm, and the ASFOA algorithm can realize the adaptive change of the search step length with the number of iterations and the change of the search result. The adaptive search step can effectively avoid the blind search of each fruit fly and overcome the difficulty that the FOA algorithm failed to execute because of the local optimum.

2.2.2. Parameters Acquisition Procedure of Hyperbolic Secant Model Based on the ASFOA Algorithm

The ASFOA algorithm is introduced into the process of obtaining the parameters of the hyperbolic secant model. The process of hyperbolic secant model parameters acquisition based on the ASFOA algorithm is shown in Figure 8. The main process is as follows:(1)Initializing the location (X-axis, Y-axis), the maximum number of iterations (Gmax), population size (N), and search step length (S) of the fruit fly population. The initial iteration values of the parameters are determined based on the empirical values of the mining area.(2)Calculating the flavor concentration value. First, the distance (Di) between the individual fruit fly and the origin is calculated, then the flavor concentration (fi) determination value is obtained. The predicted values are calculated based on the initialization parameters of the hyperbolic secant model and the information on the working surface.(3)Establishing the fitness function (F). We assume that the predicted subsidence and horizontal movements are Wi and Ui, the measured subsidence and horizontal movement are Wm and Um, respectively, and the minimum sum of squares of the differences between the predicted and measured values are used as the criterion for a determination as shown in the following formula:(4)Let the current number of iterations be , the search step of iterations be , and the optimal concentration value of iterations be . Then the adaptive step formula is(5)The optimal model prediction parameters and fruit fly locations at the current number of iterations are updated and recorded.(6)The relationship between the current number of iterations and the maximum number of iterations Gmax is judged, and the optimal solution is the output when the set condition is satisfied, else, we skip from step 6 to step 2 and continue the iterative loop operation.

3. Simulation of Experiments

3.1. Overview of the Simulated Working Face

Geological mining conditions of the simulated working face: the working face is a rectangular working face with a strike length D3 of 800 m, an inclination length D1 of 200 m, an average mining depth H of 400 m, an average mining thickness m of 3 m coal seam, and an inclination angle α of 4°. According to the experience of the thick loose seam mine area, the parameters of the hyperbolic secant model are selected as follows: ratio factor P is 0.3, subsidence factor q is 1.1, horizontal movement factor b is 0.3, main influence angle tangent tanβ1 and tanβ2 are 4.5 and 1.5, respectively, mining influence propagation angle θ is 86°, and the inflection point offset distance S1 = S2 = S3 = S4 = 30 m. A monitoring line is designed for each strike and the inclination above the mining area. The strike observation line is 1500 m long, and the inclination observation line is 780 m long. There are 51 and 27 monitoring points along the strike and inclination direction, respectively, accumulating 78 monitoring points with a distance of 30 m between the adjacent monitoring points. The maximum number of iterations (Gmax) was set to 100 and the population size (N) of the fruit fly was 10. The relationship between the working face and the monitoring points is shown in Figure 9.

Based on the geological mining conditions and the predicted parameters of the hyperbolic secant model, the measured subsidence and horizontal movement values of the simulated working face can be obtained as shown in Figure 10.

3.2. Multiview Analysis
3.2.1. Monitoring Point Antideficiency Analysis

Monitoring point data collection typically lasts for a long period. Continuous data collection for a long period requires monitoring points to be stored for a long period. However, in the production process, it is common that some of the monitoring points are damaged due to natural or human factors. The damage caused by monitoring points is considered as data missing, which affects the results of parameter acquisition to a great extent.

We used the data as new monitoring point data after missing 10%, 20%, and 30% of the data randomly, respectively. The FOA algorithm and the ASFOA algorithm are used to investigate the effect of missing monitoring points. The absolute error (AE) and relative error (RE) of the parameter acquisition results were calculated. The results of the parameter acquisition are shown in Table 1.

As shown in Table 1, in the case of random missing monitoring point data, the AE error and RE error of the parameters obtained by the FOA algorithm and the ASFOA algorithm increased with the increase in the missing point data. In the cases of 10%, 20%, and 30% missing monitoring point data, respectively, the AE error and RE error of the parameters obtained by the ASFOA algorithm are better than those of the FOA algorithm. In the case of 30% missing monitoring point data, the RE errors of the parameters obtained by the FOA algorithm and the ASFOA algorithm are higher and can hardly meet the needs of engineering practice. Therefore, the ASFOA algorithm should not be used to obtain the hyperbolic secant model parameters when the missing rate of monitoring point data exceeds 20%. In the case of 20% missing monitoring point data, the maximum RE of the parameter obtained by the ASFOA algorithm is 9.2553%, which is mainly due to the small parameter value of parameter tanβ1. Currently, the AE error of parameter tanβ1 is only 0.1388. The maximum RE of the parameter obtained by the ASFOA algorithm is 4.9001% for a 10% missing data rate at monitoring points, which is mainly due to the small parameter value of parameter tanβ1. Currently, the AE error of parameter tanβ1 is only 0.0624. Therefore, in comparison to FOA, the ASFOA algorithm has a better ability to resist the missing error of monitoring points.

3.2.2. Antigross Error Analysis

Since the errors at the inflection point and the maximum subsidence value have the greatest influence on the results, 200 mm of gross error is added at the inflection point and the point of maximum subsidence point, respectively. The parameters of the hyperbolic secant model were obtained using the gross error processed data as shown in Table 2.

As shown in Table 2, the accuracy of hyperbolic secant parameters obtained by both the FOA algorithm and the ASFOA algorithm is reduced under the premise of adding errors at the inflection point and maximum subsidence point. However, the results obtained by the ASFOA algorithm are more accurate than the FOA algorithm. In the case of adding 200 mm roughness at the inflection point, the RE error of the obtained parameters is within 0.3598%∼6.1267%. And in the case of adding 200 mm roughness at the maximum subsidence value point, the RE error of the obtained parameters is within 0.0217%∼5.8163%, which indicates that the roughness affects the difference of the results. But the AE error of the parameters is smaller, which indicates that the ASFOA algorithm has a stronger ability to resist gross error.

3.2.3. Parameter Stability Analysis of the Hyperbolic Secant Model Obtained Using the ASFOA Algorithm

To further illustrate the stability of the ASFOA algorithm to obtain the hyperbolic secant model parameters, the FOA algorithm and the ASFOA algorithm were used to obtain the hyperbolic secant model parameters based on the measured subsidence and horizontal movement values of the simulated working face. To avoid the influence of random errors, 60 parameter-finding experiments were conducted under the same experimental condition. The fluctuation in parameters observed by the FOA algorithm and the ASFOA algorithm is shown in Figures 11 and 12.

From Figure 11(a), it can be seen that the ratio factor P obtained by the ASFOA algorithm basically fluctuates within ±0.2, while the ratio factor P obtained by the FOA algorithm basically fluctuates within ±0.4. From Figure 11(b), it can be seen that the subsidence factor q obtained by the ASFOA algorithm basically fluctuates within ±0.05, while the subsidence factor q obtained by the FOA algorithm basically fluctuates within ±0.1. From Figure 11(c), it can be seen that the horizontal movement factor b obtained by the ASFOA algorithm basically fluctuates within ±0.01, while the horizontal movement factor b obtained by the FOA algorithm basically fluctuates within ±0.05. From Figure 11(d), it can be seen that the main influence angle tangent tanβ1 obtained by the ASFOA algorithm basically fluctuates within ±0.1, while the main influence angle tangent tanβ1 obtained by the FOA algorithm basically fluctuates within ±0.2. From Figure 11(e), it can be seen that the main influence angle tangent tanβ2 obtained by the ASFOA algorithm basically fluctuates within ±0.1, while the main influence angle tangent tanβ2 obtained by the FOA algorithm basically fluctuates within ±0.3. From Figure 11(f), it can be seen that the inflection point offset distance S1 obtained by the ASFOA algorithm basically fluctuates within ±5 m, while the inflection point offset distance S1 obtained by the FOA algorithm basically fluctuates within ±8 m.

From Figure 12(a), it can be seen that the inflection point offset distance S2 obtained by the ASFOA algorithm basically fluctuates within ±4 m, while the inflection point offset distance S2 obtained by the FOA algorithm basically fluctuates within ±9 m. From Figure 12(b), it can be seen that the inflection point offset distance S3 obtained by the ASFOA algorithm basically fluctuates within ±3 m, while the inflection point offset distance S3 obtained by the FOA algorithm basically fluctuates within ±8 m. From Figure 12(c), it can be seen that the inflection point offset distance S4 obtained by the ASFOA algorithm basically fluctuates within ±3 m, while the inflection point offset distance S4 obtained by the FOA algorithm basically fluctuates within ±8 m. From Figure 12(d), it can be seen that the mining influence propagation angle θ obtained by the ASFOA algorithm basically fluctuates within ±0.02°, while the mining influence propagation angle θ obtained by the FOA algorithm basically fluctuates within ±0.1°.

In summary, the hyperbolic secant model parameters obtained by the ASFOA algorithm are significantly better than the result obtained by the FOA algorithm. The hyperbolic secant model parameters obtained by the ASFOA algorithm have a smaller fluctuation range, and the hyperbolic secant model parameters obtained using the ASFOA algorithm have better stability.

4. Engineering Applications

4.1. Overview

Pansidong coal mine well field is in Huainan City, Anhui Province, China. The ground elevation is normally between +21.00∼+23.00 m. A 11111 working face is in the Pansidong mine. The working face has an inclination length of 145 m, a strike length of 411 m, and a retrieval area of 59,595 m2. The average mining depth of the working face is 394 m, the thickness of the loose seam is 336 m, the inclination of the coal seam is 6°, and the mining thickness is 4.8 m. Half of the strike observation line and half of the inclination observation line are set on the surface above the 11111 working face along the direction of strike and inclination. There are 24 monitoring points in the strike observation line and 26 monitoring points in the inclination direction. The diagram of monitoring points on the working face is shown in Figure 13.

4.2. Hyperbolic Secant Model Parameter Acquisition

The inversion of the parameters of the probability integral method and hyperbolic secant model was carried out based on the surface movement observation data of the monitoring stations. The parameters obtained are shown in Table 3, and the predicted value of the probability integral model, the predicted value of the hyperbolic secant model, and the measured value are shown in Figure 14.

As shown in Figure 14, the subsidence and horizontal movement values predicted by the hyperbolic secant model are closer to the measured values compared to those predicted by the probability integral model.

5. Discussion

As shown in Figure 14, there is still some variability and discrepancies between the measured values and predicted values of the hyperbolic secant model. The possible reasons for these discrepancies are shown in and mainly in.(1)The superimposed influence of the surrounding working face. Because in the process of mining practice, the surface monitoring station is often affected by the superposition of multiple working faces or the superposition of overlying working faces. This superposition on working faces makes it very difficult to acquire the surface movement observation value of the working face as it increases the observation error of the monitoring point.According to the geological mining information provided by the mine, the mining situation around the 11111 working face is shown in Figure 15. A 11113 working face is above the 11111 working face. The 11113 working face was completed on May 30, 2013. The 11213 working face is to the north of the 11111 working face. The 11213 working face was completed on December 30, 2015. The presence of mining activities on the 11113 and 11213 workings makes it impossible to obtain accurate values of surface movement deformation caused by mining on the 11111 workings.(2)The measurement error of the monitoring point. Due to the complexity of the site observation environment (high-diving level collapse area, cultivation area, mountainous area, etc.), there are some inevitable errors in the process of measurement of monitoring points.Pansidong coal mine is located in the eastern part of China and is a typical coal-grain complex area. The monitoring points placed above the working face are often disturbed by farmers’ farming. Some monitoring points are even damaged or directly dug out by farmers. The 11111 working face part of the damaged monitoring points is shown in Figure 16. The 11111 working face surface movement deformation measurement process of monitoring point disturbance damage makes the measurement results inevitably biased.(3)The error of the hyperbolic secant model itself. In the process of formula derivation, the default engineering geology is homogeneous, but the actual mining environment is complex and variable (faults and water-richness, etc.), which is a complex and variable environment that is not considered. The complexity of the mining environment makes it the hyperbolic secant model hard to predict the accurate value of surface movement deformation caused by mining activities.

The hyperbolic secant model in this paper does not consider the water loss consolidation settlement of the aquifer. The 11111 working face of the lower part of the loose layer is composed of aquifer and water barrier interaction, which is typical of a work face of the thick loose aquifer. The following research work can be taken into account the subsidence generated by water loss consolidation based on the research in this paper.

6. Conclusions

(1)Constructing unit mining impact function based on hyperbolic secant function. The hyperbolic secant model is proposed based on the unit mining impact function and the idea of the combined model. Based on the principle of superposition of equal influence of mining subsidence, the prediction formulas of the main section are derived at semi-infinite mining, the main section at finite mining, and any point.(2)The simulation experimental results show that the parameters of the hyperbolic secant model using the ASFOA algorithm have good stability and accuracy. The engineering application examples show that the hyperbolic secant model has high accuracy for the prediction of the thick loose layer mining area.(3)The superimposed influence of the surrounding working surface, the measurement error of the monitoring point, and the error of the model itself all affect the accuracy of the prediction results. The following research work can take into account the subsidence generated by water loss consolidation based on the research in this paper.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare no conflicts of interest.

Authors’ Contributions

Jinman Zhang wrote the manuscript. Jiewei Li and Liangji Xu verified the formula derivation and data-processing parts. Ruirui Xu entered the equations in the manuscript. Caiya Yue suggested changes to the graphs of the manuscript.

Acknowledgments

The support of the National Natural Science Foundation of China (41472323), the Foreign Science and Technology Cooperation Program of Anhui Province (201904311020015), and Startup Fund of Liaocheng University for Doctoral Research (318052115) are gratefully acknowledged. The authors are grateful to Pan Sidong Mine and Anhui University of Science and Technology for their data support. In addition, The authors would like to give special thanks to Dr. Zhou Bang for his great help during the debugging process of the code of this paper.