# Prediction of the Height of Water Flowing Fractured Zone Based on the MPSO-BP Neural Network Model

Authors: Xing Xu; Xingzhi Wang; Yuanzhi Li
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

The height of the water flowing fractured zone (HWFFZ) is the key parameter for safe coal mining under water bodies. In order to improve the accuracy of prediction model of the HWFFZ, the grey relational analysis (GRA) is used to quantitatively analyze the influencing factors of the HWFFZ, the main influencing factors of the HWFFZ under the condition of comprehensive mechanized mining and comprehensive top coal caving mining are selected scientifically, and these five factors that include mining depth, mining thickness, inclined length of working face, hard rock lithology proportion coefficient, and coal seam dip angle are selected as the discriminant indexes for predicting the HWFFZ. Nonlinear inertia weight adjustment, linear adjustments of local acceleration factor and global acceleration factor, and Gaussian mutation operation are introduced into particle swarm optimization (PSO) algorithm to improve the performance of PSO algorithm. The modified particle swarm optimization (MPSO) algorithm is used to optimize the initial weights and thresholds of back propagation (BP) neural network, and the prediction model of the HWFFZ based on MPSO-BP neural network algorithm is established according to the measured data. The prediction results are compared with those of PSO-BP neural network model and empirical formulas, the results show that the prediction accuracy of MPSO-BP neural network prediction model is higher than that of PSO-BP neural network prediction model and empirical formulas, and its prediction mean relative error (MRE) is smallest. The model also has good application effect and certain guiding significance for actual production.


---

## Body

## 1. Introduction

The movement and damage of overlying strata caused by coal seam mining are divided into caving zone, fractured zone, and bending subsidence zone from bottom to top according to the damage degree of overlying strata, namely the “upper three zones” model in China [1, 2]. There are many gaps and strong connectivity between rock blocks in the caving zone, which is the channel for the overlying water and water sand to break into the underground. The fractured zone that locates above the caving zone has a water flowing fracture connected with the goaf. The fractured zone spreads to the overlying water, and the water can be imported into the underground. Because both caving zone and fractured zone can conduct water, they are collectively called the water flowing fractured zone. Once the water flowing fractured zone connects with the overlying water bodies, it is likely to lead to the occurrence of water inrush accident in the working face, which will not only cause casualties but also cause the loss of resources [3–5]. The HWFFZ is a key hydrogeological parameter for safe mining under water bodies [6, 7]; therefore, accurately predicting the HWFFZ has important practical significance for safe mining under water bodies.The HWFFZ of working face which will be mined can only be predicted, and the prediction methods mainly include theoretical calculation method, similar material simulation method, numerical simulation method, and empirical formula method [8–10]. Among them, the model of theoretical calculation method is too idealized and deviates greatly from the complex actual geological occurrence conditions [11, 12]; Similar material simulation method requires high accuracy of material proportioning, and it is difficult to simulate some complex geological conditions [13, 14]. The accuracy of numerical simulation method is closely related to the geological condition parameters of the model, and it is difficult to accurately obtain these parameters [15, 16]. Empirical formula method in the Specification for Coal Pillar Retention and Coal Mining under Buildings, Water Bodies, Railways and Main Roadways [17] is mainly used to estimate the HWFFZ. Empirical formulas are only related to mining thickness, in fact, the HWFFZ of coal seam roof is influenced by many factors, such as mining depth, coal seam dip angle, uniaxial compressive strength of roof rock, mining thickness, mining method and inclined length of working face, and so on. In addition, according to the uniaxial compressive strength of roof rock, the roof is divided into four types including hard type, medium hard type, weak type, and extremely weak type. However, due to the multilayer structure of rock stratum, it is difficult to standardize the division of the four types in specific application. In addition, empirical formulas are mainly obtained by fitting the measured data, although they are simple and practical, they also have the disadvantage of large error. Therefore, the accuracy of empirical formulas cannot meet the requirements of efficient and safe production. How to improve the prediction accuracy of the HWFFZ and establish a fast and effective prediction model has always been an important subject for scientific and technological workers.Artificial intelligence is a new computing method, which can approach various inputs and outputs infinitely and establish various nonlinear relationships, so it has broad application prospects [18–23]. At present, many different forms of prediction model of the HWFFZ have been established based on artificial intelligence; these models have high accuracy and have been put into the link of guiding actual production [24–27]. Chen et al. [28], Ma et al. [29], Ma et al. [30], Li et al. [31], and Zhao et al. [32] used multiple factors as input indexes and established BP neural network prediction model of the HWFFZ with different parameters; they successfully applied it in engineering. BP neural network has the advantages of self-learning and self-adapting, strong nonlinear mapping ability, and generalization ability; however, BP neural network has the problems of slow convergence and easy to fall into local minimum, which will have a certain impact on the prediction accuracy of the HWFFZ. At present, in terms of prediction model of the HWFFZ, only Lou and Tan [33] used PSO algorithm to optimize BP neural network and established a prediction model of the HWFFZ based on PSO-BP neural network algorithm in China. The application results show that the model can effectively predict the HWFFZ. As the particles in PSO algorithm all converge to the optimal position, a rapid convergence effect of particle population will be formed, and it is easy to fall into local optimization and premature convergence. At the same time, since there is no relevant reference on improving PSO-BP neural network algorithm in the field of coal mine water disaster prevention, we hope to further enhance the performance of PSO-BP neural network algorithm through the improved methods; the establishment of prediction model of the HWFFZ with faster convergence speed and higher accuracy has certain reference significance for the accurate prediction of the HWFFZ and the promotion of the application range of the model.The remainder of this paper is organized as follows: in Section2, we first give a review of the related works about BP neural network and PSO algorithm, then we propose the MPSO-BP neural network algorithm and present its implementation steps. In Section 3, we first determine the main influencing factors of the HWFFZ by GRA, then we establish the MPSO-BP neural network prediction model based on the measured data of the HWFFZ. In Section 4, we first test the MPSO-BP neural network prediction model and compare with PSO-BP neural network prediction model and empirical formulas, then a case application is provided to examine the effectiveness of the MPSO-BP neural network prediction model. In Section 5, we make some conclusions.


## 2. MPSO-BP Neural Network Algorithm

### 2.1. BP Neural Network

BP neural network was proposed by Rumelhart and McClelland in 1986, and it is also the most widely used multilayer feedforward neural network [34]. BP neural network has strong nonlinear mapping ability, fault tolerance, and robustness. It is mainly composed of input layer, one or more hidden layers and output layer, each layer is fully connected, and neurons in the same layer are not connected. The forward transmission of the working signal and the back propagation of the error signal are its main characteristics, in the forward transmission stage, the working signal is processed by the input layer through the hidden layer, and the final result is output through the output layer through the conversion of the transfer function. If there is a certain error between the actual output and the expected output, it enters the back propagation stage of the error signal and uses the gradient descent method to continuously adjust the weights and thresholds of the network according to the prediction error until the error reaches the expected value, so as to make BP neural network achieve the best performance. The three-layer BP neural network that can approximate any given continuous function with arbitrary accuracy is the most commonly used and mature; its structure is shown in Figure 1.Figure 1
The three-layer BP neural network structure.The three-layer BP neural network can briefly express relationships between input and output with the following equations; here, the threshold is written to the connection weight, then the outputzk of the kth k=1,2,⋯,q hidden layer node and the output yj of the jth j=1,2,⋯,m output layer node are expressed as follows:(1)zk=f∑i=1nvki∗xi,yj=f∑k=1qwjk∗zk,where f⋅ is the transfer function; n is the number of input layer nodes, i=1,2,⋯,n; q is the number of hidden layer nodes; m is the number of output layer nodes; xi is the ith vector of the network input; vki is the connection weight from the ith input layer node to the kth hidden layer node; wjk is the connection weight from the kth hidden layer node to the jth output layer node.Lettj is the expected value of the jth output layer node, error function E is expressed as follows:(2)E=12∑j=1mtj−f∑k=1qwjk∗f∑i=1nvki∗xi2.According to the gradient descent method, the correction amountΔwjk of the connection weight wjk and the correction amount Δvki of the connection weight vki are modified, equations are expressed as follows:(3)Δwjk=−η∗∂E∂wjk,Δvki=−η∗∂E∂vki,where η is the learning step.


### 2.2. MPSO Algorithm

#### 2.2.1. PSO Algorithm

PSO algorithm is a meta heuristic algorithm proposed by Kennedy and Eberhart in 1995 [35]. It belongs to evolutionary algorithm and has the advantages of simplicity, speed, and strong searchability. Each particle in the particle population is a search individual, and the current corresponding position of each particle is a candidate solution. The optimal solution searched by each particle is the current local optimal solution; the global optimal solution is to obtain the best value by comparing the local optimal solutions in the whole particle population. Each particle has two parameters: speed and position; the speed of flight represents the speed of the particle, and the position is the position of the particle, by constantly updating these two parameters, new particles can be obtained, and the global optimal solution can be obtained by continuous iteration.LetZi=zi1,zi2,⋯,ziN be the N-dimensional position vector of the ith (i = 1, 2, …, n) particle, whether it is the optimal solution can be judged by the value of fitness function, Vi=vi1,vi2,⋯,viN is the flight speed of the ith particle, Pi=pi1,pi2,⋯,piN is the best position searched by the ith particle so far, Pg=pg1,pg2,⋯,pgN is the best position searched by the whole particle population so far, particles update their speed and position by tracking Pi and Pg, update equations of the speed and position of the particle are expressed as follows:(4)vijk+1=ω∗vijk+c1∗r1∗pijk−zijk+c2∗r2∗pgjk−zijk,(5)zijk+1=zijk+vijk+1.where k is the current iteration k; i is the number of particles, i=1,2,⋯,n; j is the dimension of the particle, j=1,2,⋯,N; ω is inertia weight; r1,r2∈0,1, both of them are random numbers; c1 is the local acceleration factor, c2 is the global acceleration factor, both of them are nonnegative numbers; vmax is the maximum limiting speed, vijk≤vmax.


#### 2.2.2. Improved Methods of PSO Algorithm

Omkar et al. [36] proved that the ergodicity of the basic PSO algorithm in the solution space cannot be guaranteed, so it cannot guarantee the convergence to the global optimization in theory. In PSO algorithm, equation (4) which is about speed evolution equation is mainly composed of three parts, the first part which is ωvijk represents the inertial effect of the particle’s previous speed, and its size represents the trade-off between global search and local search, Part 2 which is c1∗r1∗pijk−zijk represents the behavior of the particle to recognize its own strength, and Part 3 which is c2∗r2∗pgjk−zijk represents the behavior of the particle to recognize the surrounding environment. If equation (4) only retains Parts 2 and 3, the speed of the particle depends on the historical optimal position of the individual and the particle population, and it loses the inheritance of the previous particle speed, resulting in the particle only following the optimal position of the current particle population, thus losing the ability to find a better position and making PSO algorithm fall into a local minimum; therefore, Parts 2 and 3 in equation (4) play a vital role in the performance of PSO algorithm. At the beginning stage of PSO algorithm, particles should search for the optimal solution in a large range, and the movement of particles should be based on their own recognization, in order to maintain the diversity of particles, the inertia weight and the local acceleration factor should be larger, but the global acceleration factor should be smaller. At the later stage of PSO algorithm, the particles should be more controlled by the global optimal solution at this time, and fine local search should be carried out near it, so the inertia weight and the local acceleration factor should be smaller, but the global acceleration factor should be larger. At the same time, it is necessary to ensuring the diversity of the particle population and expand the search space of the particles, which can improve the possibility of the particles to find the better value. According to these ideas, we propose MPSO algorithm by improving the parameters of PSO algorithm.(1)
Nonlinear inertia weight adjustment: PSO algorithm has the problems of premature maturity and oscillation near the optimal position in the later stage. A largerω is conducive to give full play to the global searchability of PSO algorithm, while a smaller ω is conducive to give full play to the local searchability of PSO algorithm [37, 38]. Because linear inertia weigh adjustment which is commonly used has the disadvantage of low search accuracy in the later stage, it cannot balance the global and local searchability [39, 40]. Therefore, the exponential nonlinear adjustment is used to adjust ω [41], equation is expressed as follows:(6)ω=ωmax∗exp−ρ∗kT,whereωmax is the initial inertia weight, which is taken as 0.9; ρ is the control parameter, which is taken as 5.0; k is the current number of iterations; T is the maximum number of iterations.The nonlinear adjustment has largeω at the beginning of the search and strong global searchability, which helps to improve the convergence speed of PSO algorithm. In the middle and later stages of the search, ω decreases slowly and gradually turns to local search; it has the local fine searchability.(2)
Linear adjustments of local acceleration factor and global acceleration factor: local acceleration factorc1 and global acceleration factor c2 determine the step size of particle flight; if the step size is too large, the optimal position may be missed, and if it is too small, the algorithm may fall into local optimization or even nonconvergence [42]. Kennedy and Eberhart [35] suggested that they should be taken as 2.0, some scholars believe that they should be taken as 1.49 [43, 44], and some scholars have proposed other improvement methods [45, 46]. In order to maintain the diversity of particle population and improve the convergence speed, a dynamic adjustment is proposed here, c1 decreases linearly from 2.5 to 0.5, while c2 increases linearly from 0.5 to 2.5, equations are expressed as follows:(7)c1=2.5−2.0∗k−1T−1,(8)c2=0.5+2.0∗k−1T−1,wherek is the current number of iterations; T is the maximum number of iterations.According to the above adjustment,c1 should be greater than c2 at the beginning stage of PSO algorithm, c1 should be smaller than c2 at the later stage of PSO algorithm. That is, particles focus on recognizing their own strength at the beginning stage of PSO algorithm, and particles focus on recognizing the surrounding environment at the later stage of PSO algorithm, so as to ensure the dynamic allocation of particles’ self-recognization ability and global recognization ability.(3)
Gaussian mutation operation: in the process of PSO algorithm, on the one hand, the diversity of the particle population will gradually weaken; on the other hand, PSO algorithm is also in the local fine search stage, the particle speed tends to zero, and it is difficult to jump out of the local optimal position. In order to solve the possible premature convergence problem, based on the idea of mutation in genetic algorithm, a mutation operation is introduced into PSO algorithm to reinitialize particles with a certain probability. Mutation operation expands the particle population search space that is shrinking in the iteration, so that the particle can jump out of the position of the previously searched optimal value and search in a larger space; at the same time, it maintains the particle population diversity and improves the possibility of the PSO algorithm to find a better value [47, 48]. Scholars have designed many mutation operations, such as uniform mutation operation, nonuniform mutation operation, adaptive mutation operation, and normal mutation operation (Gaussian mutation operation), among them. Gauss mutation operation is widely used. Gaussian mutation operation has been proved to enhance the local searchability of meta heuristic algorithm. Qu and He [49] proposed a novel artificial fish-school algorithm based on adaptive Gauss mutation and historical best fish; the tests show that the algorithm is feasible and effective. Mo et al. [50] proposed an artificial glowworm optimization algorithm with Gaussian mutation to overcome the shortcomings of the basic glowworm optimization algorithm in solving the global optimal value of the function; the results show that the improved artificial glowworm algorithm has higher convergence speed and convergence success rate than the basic glowworm optimization algorithm. Gao et al. [51] proposed a multiobjective PSO algorithm based on Gaussian mutation and adaptive reference point fusion; the experimental results show that the improved multiobjective PSO algorithm has better diversity and convergence than the traditional multiobjective PSO algorithm. Zheng et al. [52] proposed the shuffled frog leaping and bat algorithm with Gauss mutation to solve the complex functions; the test results show that the algorithm has better performance and is suitable for various high-dimensional and multiextremum complex function optimization problems in engineering applications. Gaussian mutation operation is derived from Gaussian distribution. Gaussian distribution is an important probability distribution in theoretical research and engineering application; it can generate random variables within a certain range according to the mean and the standard deviation [53–56]. The density function fx of Gaussian distribution is expressed as follows:(9)fx=1σ∗2π∗e−x−μ2/2σ2,wherex is a random variable, x∈−∞,∞; μ is the mean value; σ is the standard deviation.The probability density function image of Gaussian distribution is symmetrical about the expectation, and about 99.8% of the area under the function curve is within the range of about three times the standard deviation on both sides of the mean value, that is, when random variables are generated, the variables can be concentrated near the mean value. Here, a variation thresholdpm is set as the condition of Gaussian mutation operation, when the condition of Gaussian mutation operation is met, namely, the inequality rand≥pm is met, the particle will have a position mutation to make the particle jump out of the current position, pm can be set an appropriate value according to the actual situation, rand is a random number which belong to [0, 1]. The method of adding random disturbance is used to vary the position zijk of the particle, equation of the Gaussian mutation operation is expressed as follows:(10)zijk+1=zijk∗1+0.5δ,where i=1,2,⋯,n; j=1,2,⋯,N; k=1,2,⋯,T; δ is a random variable which obeys the density function of Gaussian distribution with μ=0 and σ = 1.


### 2.3. Implementation of MPSO-BP Neural Network Algorithm

First, MPSO algorithm is used to optimize BP neural network parameters, then BP neural network is used to further optimize its network parameters, so as to establish MPSO-BP neural network algorithm combined with MPSO algorithm and BP neural network. This hybrid algorithm can optimize the initial weights and thresholds of BP neural network and improve the output accuracy of BP neural network through the first global optimization and the later local optimization.The implementation steps of MPSO-BP neural network algorithm are as follows:Step 1. Initialize the particle population. Determine the variableN, the variable n, the variable vmax and the variable pm, initialize the position and speed of particles randomly. Determine ω according to equation (6), determine c1 and c2 according to equations (7) and (8).Step 2. Evaluate the optimization effect of MPSO algorithm. Use the mean absolute error (MAE) as the fitness function, calculate the fitness of the current particle population with the fitness function, then let the individual optimal value of the current particle bePi and the historical optimal value of the current particle population be Pg. MAE is expressed as follows:(11)MAE=1n∗∑l=1nyl−Yl,wheren is the total number of training samples; yl is the prediction value of the lth sample, m; Yl is the measured value of the lth sample, m.Step 3. Update the particle speed according to equation (4).Step 4. Generate the random numberrand, when the inequality r and≥pm is met, generate the new position of the particle according to equation (10), replace the current position of the particle, and calculate the fitness value of particle population, go to Step: 6. Otherwise, go directly to Step: 5.Step 5. Update the position of the particle according to equation (5), calculate the fitness value of particle population.Step: 6. Update the individual historical optimal position and the particle population historical optimal position of each particle.Step: 7. Repeat Step: 3 to Step: 6; if the termination condition is reached, stop the calculation, and output the optimal solutionPg, establish the optimized BP neural network.The flow chart of MPSO-BP neural network algorithm is shown in Figure2.Figure 2
Flow chart of MPSO-BP neural network algorithm.


## 2.1. BP Neural Network

BP neural network was proposed by Rumelhart and McClelland in 1986, and it is also the most widely used multilayer feedforward neural network [34]. BP neural network has strong nonlinear mapping ability, fault tolerance, and robustness. It is mainly composed of input layer, one or more hidden layers and output layer, each layer is fully connected, and neurons in the same layer are not connected. The forward transmission of the working signal and the back propagation of the error signal are its main characteristics, in the forward transmission stage, the working signal is processed by the input layer through the hidden layer, and the final result is output through the output layer through the conversion of the transfer function. If there is a certain error between the actual output and the expected output, it enters the back propagation stage of the error signal and uses the gradient descent method to continuously adjust the weights and thresholds of the network according to the prediction error until the error reaches the expected value, so as to make BP neural network achieve the best performance. The three-layer BP neural network that can approximate any given continuous function with arbitrary accuracy is the most commonly used and mature; its structure is shown in Figure 1.Figure 1
The three-layer BP neural network structure.The three-layer BP neural network can briefly express relationships between input and output with the following equations; here, the threshold is written to the connection weight, then the outputzk of the kth k=1,2,⋯,q hidden layer node and the output yj of the jth j=1,2,⋯,m output layer node are expressed as follows:(1)zk=f∑i=1nvki∗xi,yj=f∑k=1qwjk∗zk,where f⋅ is the transfer function; n is the number of input layer nodes, i=1,2,⋯,n; q is the number of hidden layer nodes; m is the number of output layer nodes; xi is the ith vector of the network input; vki is the connection weight from the ith input layer node to the kth hidden layer node; wjk is the connection weight from the kth hidden layer node to the jth output layer node.Lettj is the expected value of the jth output layer node, error function E is expressed as follows:(2)E=12∑j=1mtj−f∑k=1qwjk∗f∑i=1nvki∗xi2.According to the gradient descent method, the correction amountΔwjk of the connection weight wjk and the correction amount Δvki of the connection weight vki are modified, equations are expressed as follows:(3)Δwjk=−η∗∂E∂wjk,Δvki=−η∗∂E∂vki,where η is the learning step.


## 2.2. MPSO Algorithm

### 2.2.1. PSO Algorithm

PSO algorithm is a meta heuristic algorithm proposed by Kennedy and Eberhart in 1995 [35]. It belongs to evolutionary algorithm and has the advantages of simplicity, speed, and strong searchability. Each particle in the particle population is a search individual, and the current corresponding position of each particle is a candidate solution. The optimal solution searched by each particle is the current local optimal solution; the global optimal solution is to obtain the best value by comparing the local optimal solutions in the whole particle population. Each particle has two parameters: speed and position; the speed of flight represents the speed of the particle, and the position is the position of the particle, by constantly updating these two parameters, new particles can be obtained, and the global optimal solution can be obtained by continuous iteration.LetZi=zi1,zi2,⋯,ziN be the N-dimensional position vector of the ith (i = 1, 2, …, n) particle, whether it is the optimal solution can be judged by the value of fitness function, Vi=vi1,vi2,⋯,viN is the flight speed of the ith particle, Pi=pi1,pi2,⋯,piN is the best position searched by the ith particle so far, Pg=pg1,pg2,⋯,pgN is the best position searched by the whole particle population so far, particles update their speed and position by tracking Pi and Pg, update equations of the speed and position of the particle are expressed as follows:(4)vijk+1=ω∗vijk+c1∗r1∗pijk−zijk+c2∗r2∗pgjk−zijk,(5)zijk+1=zijk+vijk+1.where k is the current iteration k; i is the number of particles, i=1,2,⋯,n; j is the dimension of the particle, j=1,2,⋯,N; ω is inertia weight; r1,r2∈0,1, both of them are random numbers; c1 is the local acceleration factor, c2 is the global acceleration factor, both of them are nonnegative numbers; vmax is the maximum limiting speed, vijk≤vmax.


### 2.2.2. Improved Methods of PSO Algorithm

Omkar et al. [36] proved that the ergodicity of the basic PSO algorithm in the solution space cannot be guaranteed, so it cannot guarantee the convergence to the global optimization in theory. In PSO algorithm, equation (4) which is about speed evolution equation is mainly composed of three parts, the first part which is ωvijk represents the inertial effect of the particle’s previous speed, and its size represents the trade-off between global search and local search, Part 2 which is c1∗r1∗pijk−zijk represents the behavior of the particle to recognize its own strength, and Part 3 which is c2∗r2∗pgjk−zijk represents the behavior of the particle to recognize the surrounding environment. If equation (4) only retains Parts 2 and 3, the speed of the particle depends on the historical optimal position of the individual and the particle population, and it loses the inheritance of the previous particle speed, resulting in the particle only following the optimal position of the current particle population, thus losing the ability to find a better position and making PSO algorithm fall into a local minimum; therefore, Parts 2 and 3 in equation (4) play a vital role in the performance of PSO algorithm. At the beginning stage of PSO algorithm, particles should search for the optimal solution in a large range, and the movement of particles should be based on their own recognization, in order to maintain the diversity of particles, the inertia weight and the local acceleration factor should be larger, but the global acceleration factor should be smaller. At the later stage of PSO algorithm, the particles should be more controlled by the global optimal solution at this time, and fine local search should be carried out near it, so the inertia weight and the local acceleration factor should be smaller, but the global acceleration factor should be larger. At the same time, it is necessary to ensuring the diversity of the particle population and expand the search space of the particles, which can improve the possibility of the particles to find the better value. According to these ideas, we propose MPSO algorithm by improving the parameters of PSO algorithm.(1)
Nonlinear inertia weight adjustment: PSO algorithm has the problems of premature maturity and oscillation near the optimal position in the later stage. A largerω is conducive to give full play to the global searchability of PSO algorithm, while a smaller ω is conducive to give full play to the local searchability of PSO algorithm [37, 38]. Because linear inertia weigh adjustment which is commonly used has the disadvantage of low search accuracy in the later stage, it cannot balance the global and local searchability [39, 40]. Therefore, the exponential nonlinear adjustment is used to adjust ω [41], equation is expressed as follows:(6)ω=ωmax∗exp−ρ∗kT,whereωmax is the initial inertia weight, which is taken as 0.9; ρ is the control parameter, which is taken as 5.0; k is the current number of iterations; T is the maximum number of iterations.The nonlinear adjustment has largeω at the beginning of the search and strong global searchability, which helps to improve the convergence speed of PSO algorithm. In the middle and later stages of the search, ω decreases slowly and gradually turns to local search; it has the local fine searchability.(2)
Linear adjustments of local acceleration factor and global acceleration factor: local acceleration factorc1 and global acceleration factor c2 determine the step size of particle flight; if the step size is too large, the optimal position may be missed, and if it is too small, the algorithm may fall into local optimization or even nonconvergence [42]. Kennedy and Eberhart [35] suggested that they should be taken as 2.0, some scholars believe that they should be taken as 1.49 [43, 44], and some scholars have proposed other improvement methods [45, 46]. In order to maintain the diversity of particle population and improve the convergence speed, a dynamic adjustment is proposed here, c1 decreases linearly from 2.5 to 0.5, while c2 increases linearly from 0.5 to 2.5, equations are expressed as follows:(7)c1=2.5−2.0∗k−1T−1,(8)c2=0.5+2.0∗k−1T−1,wherek is the current number of iterations; T is the maximum number of iterations.According to the above adjustment,c1 should be greater than c2 at the beginning stage of PSO algorithm, c1 should be smaller than c2 at the later stage of PSO algorithm. That is, particles focus on recognizing their own strength at the beginning stage of PSO algorithm, and particles focus on recognizing the surrounding environment at the later stage of PSO algorithm, so as to ensure the dynamic allocation of particles’ self-recognization ability and global recognization ability.(3)
Gaussian mutation operation: in the process of PSO algorithm, on the one hand, the diversity of the particle population will gradually weaken; on the other hand, PSO algorithm is also in the local fine search stage, the particle speed tends to zero, and it is difficult to jump out of the local optimal position. In order to solve the possible premature convergence problem, based on the idea of mutation in genetic algorithm, a mutation operation is introduced into PSO algorithm to reinitialize particles with a certain probability. Mutation operation expands the particle population search space that is shrinking in the iteration, so that the particle can jump out of the position of the previously searched optimal value and search in a larger space; at the same time, it maintains the particle population diversity and improves the possibility of the PSO algorithm to find a better value [47, 48]. Scholars have designed many mutation operations, such as uniform mutation operation, nonuniform mutation operation, adaptive mutation operation, and normal mutation operation (Gaussian mutation operation), among them. Gauss mutation operation is widely used. Gaussian mutation operation has been proved to enhance the local searchability of meta heuristic algorithm. Qu and He [49] proposed a novel artificial fish-school algorithm based on adaptive Gauss mutation and historical best fish; the tests show that the algorithm is feasible and effective. Mo et al. [50] proposed an artificial glowworm optimization algorithm with Gaussian mutation to overcome the shortcomings of the basic glowworm optimization algorithm in solving the global optimal value of the function; the results show that the improved artificial glowworm algorithm has higher convergence speed and convergence success rate than the basic glowworm optimization algorithm. Gao et al. [51] proposed a multiobjective PSO algorithm based on Gaussian mutation and adaptive reference point fusion; the experimental results show that the improved multiobjective PSO algorithm has better diversity and convergence than the traditional multiobjective PSO algorithm. Zheng et al. [52] proposed the shuffled frog leaping and bat algorithm with Gauss mutation to solve the complex functions; the test results show that the algorithm has better performance and is suitable for various high-dimensional and multiextremum complex function optimization problems in engineering applications. Gaussian mutation operation is derived from Gaussian distribution. Gaussian distribution is an important probability distribution in theoretical research and engineering application; it can generate random variables within a certain range according to the mean and the standard deviation [53–56]. The density function fx of Gaussian distribution is expressed as follows:(9)fx=1σ∗2π∗e−x−μ2/2σ2,wherex is a random variable, x∈−∞,∞; μ is the mean value; σ is the standard deviation.The probability density function image of Gaussian distribution is symmetrical about the expectation, and about 99.8% of the area under the function curve is within the range of about three times the standard deviation on both sides of the mean value, that is, when random variables are generated, the variables can be concentrated near the mean value. Here, a variation thresholdpm is set as the condition of Gaussian mutation operation, when the condition of Gaussian mutation operation is met, namely, the inequality rand≥pm is met, the particle will have a position mutation to make the particle jump out of the current position, pm can be set an appropriate value according to the actual situation, rand is a random number which belong to [0, 1]. The method of adding random disturbance is used to vary the position zijk of the particle, equation of the Gaussian mutation operation is expressed as follows:(10)zijk+1=zijk∗1+0.5δ,where i=1,2,⋯,n; j=1,2,⋯,N; k=1,2,⋯,T; δ is a random variable which obeys the density function of Gaussian distribution with μ=0 and σ = 1.


## 2.2.1. PSO Algorithm

PSO algorithm is a meta heuristic algorithm proposed by Kennedy and Eberhart in 1995 [35]. It belongs to evolutionary algorithm and has the advantages of simplicity, speed, and strong searchability. Each particle in the particle population is a search individual, and the current corresponding position of each particle is a candidate solution. The optimal solution searched by each particle is the current local optimal solution; the global optimal solution is to obtain the best value by comparing the local optimal solutions in the whole particle population. Each particle has two parameters: speed and position; the speed of flight represents the speed of the particle, and the position is the position of the particle, by constantly updating these two parameters, new particles can be obtained, and the global optimal solution can be obtained by continuous iteration.LetZi=zi1,zi2,⋯,ziN be the N-dimensional position vector of the ith (i = 1, 2, …, n) particle, whether it is the optimal solution can be judged by the value of fitness function, Vi=vi1,vi2,⋯,viN is the flight speed of the ith particle, Pi=pi1,pi2,⋯,piN is the best position searched by the ith particle so far, Pg=pg1,pg2,⋯,pgN is the best position searched by the whole particle population so far, particles update their speed and position by tracking Pi and Pg, update equations of the speed and position of the particle are expressed as follows:(4)vijk+1=ω∗vijk+c1∗r1∗pijk−zijk+c2∗r2∗pgjk−zijk,(5)zijk+1=zijk+vijk+1.where k is the current iteration k; i is the number of particles, i=1,2,⋯,n; j is the dimension of the particle, j=1,2,⋯,N; ω is inertia weight; r1,r2∈0,1, both of them are random numbers; c1 is the local acceleration factor, c2 is the global acceleration factor, both of them are nonnegative numbers; vmax is the maximum limiting speed, vijk≤vmax.


## 2.2.2. Improved Methods of PSO Algorithm

Omkar et al. [36] proved that the ergodicity of the basic PSO algorithm in the solution space cannot be guaranteed, so it cannot guarantee the convergence to the global optimization in theory. In PSO algorithm, equation (4) which is about speed evolution equation is mainly composed of three parts, the first part which is ωvijk represents the inertial effect of the particle’s previous speed, and its size represents the trade-off between global search and local search, Part 2 which is c1∗r1∗pijk−zijk represents the behavior of the particle to recognize its own strength, and Part 3 which is c2∗r2∗pgjk−zijk represents the behavior of the particle to recognize the surrounding environment. If equation (4) only retains Parts 2 and 3, the speed of the particle depends on the historical optimal position of the individual and the particle population, and it loses the inheritance of the previous particle speed, resulting in the particle only following the optimal position of the current particle population, thus losing the ability to find a better position and making PSO algorithm fall into a local minimum; therefore, Parts 2 and 3 in equation (4) play a vital role in the performance of PSO algorithm. At the beginning stage of PSO algorithm, particles should search for the optimal solution in a large range, and the movement of particles should be based on their own recognization, in order to maintain the diversity of particles, the inertia weight and the local acceleration factor should be larger, but the global acceleration factor should be smaller. At the later stage of PSO algorithm, the particles should be more controlled by the global optimal solution at this time, and fine local search should be carried out near it, so the inertia weight and the local acceleration factor should be smaller, but the global acceleration factor should be larger. At the same time, it is necessary to ensuring the diversity of the particle population and expand the search space of the particles, which can improve the possibility of the particles to find the better value. According to these ideas, we propose MPSO algorithm by improving the parameters of PSO algorithm.(1)
Nonlinear inertia weight adjustment: PSO algorithm has the problems of premature maturity and oscillation near the optimal position in the later stage. A largerω is conducive to give full play to the global searchability of PSO algorithm, while a smaller ω is conducive to give full play to the local searchability of PSO algorithm [37, 38]. Because linear inertia weigh adjustment which is commonly used has the disadvantage of low search accuracy in the later stage, it cannot balance the global and local searchability [39, 40]. Therefore, the exponential nonlinear adjustment is used to adjust ω [41], equation is expressed as follows:(6)ω=ωmax∗exp−ρ∗kT,whereωmax is the initial inertia weight, which is taken as 0.9; ρ is the control parameter, which is taken as 5.0; k is the current number of iterations; T is the maximum number of iterations.The nonlinear adjustment has largeω at the beginning of the search and strong global searchability, which helps to improve the convergence speed of PSO algorithm. In the middle and later stages of the search, ω decreases slowly and gradually turns to local search; it has the local fine searchability.(2)
Linear adjustments of local acceleration factor and global acceleration factor: local acceleration factorc1 and global acceleration factor c2 determine the step size of particle flight; if the step size is too large, the optimal position may be missed, and if it is too small, the algorithm may fall into local optimization or even nonconvergence [42]. Kennedy and Eberhart [35] suggested that they should be taken as 2.0, some scholars believe that they should be taken as 1.49 [43, 44], and some scholars have proposed other improvement methods [45, 46]. In order to maintain the diversity of particle population and improve the convergence speed, a dynamic adjustment is proposed here, c1 decreases linearly from 2.5 to 0.5, while c2 increases linearly from 0.5 to 2.5, equations are expressed as follows:(7)c1=2.5−2.0∗k−1T−1,(8)c2=0.5+2.0∗k−1T−1,wherek is the current number of iterations; T is the maximum number of iterations.According to the above adjustment,c1 should be greater than c2 at the beginning stage of PSO algorithm, c1 should be smaller than c2 at the later stage of PSO algorithm. That is, particles focus on recognizing their own strength at the beginning stage of PSO algorithm, and particles focus on recognizing the surrounding environment at the later stage of PSO algorithm, so as to ensure the dynamic allocation of particles’ self-recognization ability and global recognization ability.(3)
Gaussian mutation operation: in the process of PSO algorithm, on the one hand, the diversity of the particle population will gradually weaken; on the other hand, PSO algorithm is also in the local fine search stage, the particle speed tends to zero, and it is difficult to jump out of the local optimal position. In order to solve the possible premature convergence problem, based on the idea of mutation in genetic algorithm, a mutation operation is introduced into PSO algorithm to reinitialize particles with a certain probability. Mutation operation expands the particle population search space that is shrinking in the iteration, so that the particle can jump out of the position of the previously searched optimal value and search in a larger space; at the same time, it maintains the particle population diversity and improves the possibility of the PSO algorithm to find a better value [47, 48]. Scholars have designed many mutation operations, such as uniform mutation operation, nonuniform mutation operation, adaptive mutation operation, and normal mutation operation (Gaussian mutation operation), among them. Gauss mutation operation is widely used. Gaussian mutation operation has been proved to enhance the local searchability of meta heuristic algorithm. Qu and He [49] proposed a novel artificial fish-school algorithm based on adaptive Gauss mutation and historical best fish; the tests show that the algorithm is feasible and effective. Mo et al. [50] proposed an artificial glowworm optimization algorithm with Gaussian mutation to overcome the shortcomings of the basic glowworm optimization algorithm in solving the global optimal value of the function; the results show that the improved artificial glowworm algorithm has higher convergence speed and convergence success rate than the basic glowworm optimization algorithm. Gao et al. [51] proposed a multiobjective PSO algorithm based on Gaussian mutation and adaptive reference point fusion; the experimental results show that the improved multiobjective PSO algorithm has better diversity and convergence than the traditional multiobjective PSO algorithm. Zheng et al. [52] proposed the shuffled frog leaping and bat algorithm with Gauss mutation to solve the complex functions; the test results show that the algorithm has better performance and is suitable for various high-dimensional and multiextremum complex function optimization problems in engineering applications. Gaussian mutation operation is derived from Gaussian distribution. Gaussian distribution is an important probability distribution in theoretical research and engineering application; it can generate random variables within a certain range according to the mean and the standard deviation [53–56]. The density function fx of Gaussian distribution is expressed as follows:(9)fx=1σ∗2π∗e−x−μ2/2σ2,wherex is a random variable, x∈−∞,∞; μ is the mean value; σ is the standard deviation.The probability density function image of Gaussian distribution is symmetrical about the expectation, and about 99.8% of the area under the function curve is within the range of about three times the standard deviation on both sides of the mean value, that is, when random variables are generated, the variables can be concentrated near the mean value. Here, a variation thresholdpm is set as the condition of Gaussian mutation operation, when the condition of Gaussian mutation operation is met, namely, the inequality rand≥pm is met, the particle will have a position mutation to make the particle jump out of the current position, pm can be set an appropriate value according to the actual situation, rand is a random number which belong to [0, 1]. The method of adding random disturbance is used to vary the position zijk of the particle, equation of the Gaussian mutation operation is expressed as follows:(10)zijk+1=zijk∗1+0.5δ,where i=1,2,⋯,n; j=1,2,⋯,N; k=1,2,⋯,T; δ is a random variable which obeys the density function of Gaussian distribution with μ=0 and σ = 1.


## 2.3. Implementation of MPSO-BP Neural Network Algorithm

First, MPSO algorithm is used to optimize BP neural network parameters, then BP neural network is used to further optimize its network parameters, so as to establish MPSO-BP neural network algorithm combined with MPSO algorithm and BP neural network. This hybrid algorithm can optimize the initial weights and thresholds of BP neural network and improve the output accuracy of BP neural network through the first global optimization and the later local optimization.The implementation steps of MPSO-BP neural network algorithm are as follows:Step 1. Initialize the particle population. Determine the variableN, the variable n, the variable vmax and the variable pm, initialize the position and speed of particles randomly. Determine ω according to equation (6), determine c1 and c2 according to equations (7) and (8).Step 2. Evaluate the optimization effect of MPSO algorithm. Use the mean absolute error (MAE) as the fitness function, calculate the fitness of the current particle population with the fitness function, then let the individual optimal value of the current particle bePi and the historical optimal value of the current particle population be Pg. MAE is expressed as follows:(11)MAE=1n∗∑l=1nyl−Yl,wheren is the total number of training samples; yl is the prediction value of the lth sample, m; Yl is the measured value of the lth sample, m.Step 3. Update the particle speed according to equation (4).Step 4. Generate the random numberrand, when the inequality r and≥pm is met, generate the new position of the particle according to equation (10), replace the current position of the particle, and calculate the fitness value of particle population, go to Step: 6. Otherwise, go directly to Step: 5.Step 5. Update the position of the particle according to equation (5), calculate the fitness value of particle population.Step: 6. Update the individual historical optimal position and the particle population historical optimal position of each particle.Step: 7. Repeat Step: 3 to Step: 6; if the termination condition is reached, stop the calculation, and output the optimal solutionPg, establish the optimized BP neural network.The flow chart of MPSO-BP neural network algorithm is shown in Figure2.Figure 2
Flow chart of MPSO-BP neural network algorithm.


## 3. Establishment of Prediction Model

### 3.1. Analysis of Main Influencing Factors

The formation of water flowing fractured zone of coal seam roof under the influence of mining is a complex mechanical problem. According to the Specification for Coal Pillar Retention and Coal Mining under Buildings, Water Bodies, Railways and Main Roadways and field investigation results, the HWFFZ is mainly influenced by many factors such as mining geological conditions and three-dimensional size design parameters of working face [17, 57–60]. Considering that the total caving method is widely used to manage the roof in China, this paper only studies the influence of different mining technologies on the HWFFZ under the long wall mining method. Therefore, it is determined that the influencing factors of the HWFFZ include four factors of mining geological conditions, such as mining depth, coal seam dip angle, uniaxial compressive strength of roof rock and roof strata combination structure type, and include four three-dimensional size design parameters of working face, such as mining thickness, mining method, number of mining layers, and inclined length of working face. Among them, uniaxial compressive strength of roof rock and roof strata combination structure type are replaced by the hard rock lithology proportion coefficient [6], this coefficient can reflect the overall strength of roof rock stratum and the combination structure of roof rock stratum and avoid the difficulty that the roof type cannot be determined according to the uniaxial compressive strength in the Specification for Coal Pillar Retention and Coal Mining under Buildings, Water Bodies, Railways and Main Roadways [17], and it is easy to obtain. Then influencing factors of the HWFFZ can be set as seven influencing factors including mining depth, mining thickness, hard rock lithology proportion coefficient, inclined length of working face, mining method, number of mining layers and coal seam dip angle. Among them, mining method is qualitative data, and other data are quantitative data; mining methods include blasting mining, comprehensive mechanized mining, comprehensive top coal caving mining and layered mining, with values of 0.2, 0.4, 0.6, and 0.8, respectively. By collecting and sorting 31 field measurement cases of relevant coal mines in China, the measured data of the HWFFZ are obtained. The correlation coefficient between each influencing factor of samples and the HWFFZ is calculated by GRA [61]; the broken line diagram is drawn according to the correlation coefficient, as shown in Figure 3.Figure 3
Broken line diagram of correlation coefficient.The correlation degree and correlation order between each influencing factor and the HWFFZ can be obtained, as shown in Table1.Table 1
The correlation degree and correlation order.
Mining depthMining thicknessHard rock lithology proportion coefficientInclined length of working faceMining methodNumber of mining layersCoal seam dip angle0.78930.80370.76010.79160.76040.77950.73073162547It can be seen from Figure3 and Table 1that all influencing factors have a high degree of correlation and consistency with the HWFFZ, indicating that these seven influencing factors have a certain control effect on the HWFFZ. However, it can also be seen that the correlation between mining thickness and the HWFFZ is the highest, with a correlation value of 0.8037, and the correlation between coal seam dip angle and the HWFFZ is the lowest, with a correlation value of 0.7307. According to the broken line diagram of correlation coefficient between coal seam dip angle and the HWFFZ, it can also be seen that the correlation between coal seam dip angle and the HWFFZ is the worst, and its correlation coefficient fluctuates the most. The index values of the seven influencing factors provided by samples 2 and 19 have less correlation with the HWFFZ, and the correlation has obvious consistency; however, the index values of the seven influencing factors provided by the other 29 groups of samples have great correlation with the HWFFZ, and the correlation has a certain discrete type; therefore, the GRA of measured data of the HWFFZ is effective. For further analysis, at present, the most coal mines have realized comprehensive mechanized mining and comprehensive top coal caving mining, blasting mining, and layered mining are rare, and mining thickness has a strong linear correlation with number of mining layers since the selection of main factors should be representative as much as possible and meet the requirements of linear independence between indicators; therefore, five factors including mining depth, mining thickness, hard rock lithology proportion coefficient, working face slope length, and coal seam dip angle are selected as the main influencing factors of the HWFFZ.


### 3.2. Design of Prediction Model

In total, 38 groups of measured data of the HWFFZ under the condition of comprehensive mechanized mining and comprehensive top coal caving mining of Chinese coal mines are collected as research samples; five factors including mining depth, mining thickness, hard rock lithology proportion coefficient, inclined length of working face, and coal seam dip angle are used as the input index of BP neural network, and the HWFFZ is used as the output of BP neural network. The three-layer BP network structure is adopted, and the number of hidden layer neurons is set to 6, the number of output layer neurons is set to 1. The purelin function is selected as the transfer functions of the hidden layer, the logsig function is selected as the output layer, and the trainlm function is selected as the training function. The training times are set to 1000 times, the training target error is set to 0.00001, and the learning efficiency is set to 0.04.MPSO algorithm and PSO algorithm are used to optimize BP neural network globally; according to the sum of initial weights and thresholds of BP neural network, the dimensionN of the particle is set as 43, in order to expand the search space of particles themselves. Gaussian distribution is used to perform mutation on the particles, and the threshold pm is set as 0.9. The number of particles n, the maximum number of iterations T, and the maximum limiting speed vmax are set empirically. The detailed parameter settings of the PSO algorithm and MPSO algorithm are shown in Table 2.Table 2
Parameter setting of PSO algorithm and MPSO algorithm.
AlgorithmNnTωc1c2PmvmaxPSO435050Linear decreasing22No1.9MPSO435050Exponential decreasingLinear decreasing from 2.5 to 0.5Linear increasing from 0.5 to 2.50.91.9


### 3.3. Training of Prediction Model

We use measured data of the HWFFZ to train BP neural network according to the algorithm flow in Figure2 and then establish the MPSO-BP neural network prediction model by MATLAB software. The establishment of PSO-BP neural network prediction model is similar, which is omitted here.The comparison diagram of fitness evolution curves between MPSO algorithm and PSO algorithm in the process of optimizing BP neural network is shown in Figure4.Figure 4
Comparison of fitness evolution curves.It can be seen from Figure4 that the fitness evolution curve of MPSO-BP neural network algorithm and PSO-BP neural network algorithm shows a decreasing trend. It can also be seen that two optimization algorithms have achieved the purpose of optimizing BP neural network, but MPSO algorithm is better than PSO algorithm. From the first generation to the seventh generation, the fitness value of MPSO-BP neural network algorithm evolution curve is greater than that of PSO-BP neural network algorithm evolution curve. In the seventh generation, the fitness value of MPSO-BP neural network algorithm evolution curve becomes 3.5449, while the fitness value of PSO-BP neural network algorithm evolution curve is 3.5804. Until the 10th generation, the fitness value of MPSO-BP neural network algorithm evolution curve is better than that of PSO-BP neural network algorithm evolution curve. However, from the 11th generation, the fitness evolution curve of MPSO-BP neural network algorithm presents a straight line; the value is 3.0979 and remains until the 21st generation. The fitness evolution curve of PSO-BP neural network algorithm decreases rapidly from 3.5061 to 2.8954 and shows a long straight line, no further evolution occurs until the 37th generation, which shows that PSO-BP neural network algorithm fall into local optimization and do not further evolve for a long time. The fitness evolution curve of MPSO-BP neural network algorithm rapidly decreased from 3.0978 to 2.8614 in the 20th generation, and from the 23th generation, the fitness value of evolution curve is less than that of PSO-BP neural network algorithm, and the fitness evolution curve of MPSO-BP neural network algorithm further evolves in the 43th generation, and the fitness decreases from 2.6779 to 2.3018, while the fitness evolution curve of PSO-BP neural network algorithm enters the second straight line from the 39th generation to the 46th generation, and the fitness value remains at 2.7202 and further decreases to 2.3996 in the 47th generation. It can be seen from the synthesis that fitness evolution curve of PSO-BP neural network algorithm evolves rapidly at the beginning, but it fall into local optimization from the 11th generation to the 37th generation and fails to evolve further, reflecting that the particles gathers at the same position and stops moving and refuses to search the solution space in other fields, although the fitness value is further reduced in the later stage, but PSO algorithm is weaker than MPSO algorithm in the process of global optimization. On the contrary, the MPSO-BP evolution curve keeps evolving, and there is no long-term optimization stagnation, that is, the nonevolution stage, which shows that MPSO algorithm can jump out of the local optimization and further evolve. At the end of the optimization process, its fitness value is less than that of PSO algorithm. This shows that MPSO algorithm is better than PSO algorithm, and the improved strategy makes PSO-BP further improve the performance.The comparison between measured value and prediction value of MPSO-BP neural network model and PSO-BP neural network model is shown in Figure5.Figure 5
Training comparison of the two models. (a) Comparison between measured value and prediction value of PSO-BP model. (b) Comparison between measured value and prediction value of MPSO-BP model.
(a)(b)It can be seen from Figure5(a) that the difference between prediction value of PSO-BP model and measured value is small. In 38 training samples, most of prediction value and measured value basically coincides, but there is a certain gap in individual samples, such as samples 4, 8, 10, 15, 17, and 35; its prediction error ranges from −11.37 to 11.55, but most of the prediction errors are small, and some are close to 0. It can be seen from Figure 5(b) that the difference between prediction value of MPSO-BP model and measured value is also small. In 38 training samples, most of prediction value and measured value basically coincides, and there is also a small gap in individual samples, such as samples 4, 15, 22, and 23; its prediction error ranges from −6.37 to 17.92, but most of the prediction errors are also small, and some are close to 0.In order to further analyze prediction performance of the two models, we contrast prediction stability and prediction accuracy by MRE, MRE is expressed as follows:(12)MRE=1n∗∑l=1nyl−YlYl∗100%,where yl−Yl/Yl∗100% is the relative error (RE) of the lth sample; n is the total number of training samples; yl is the prediction value of the lth sample, m; Yl is the measured value of the lth sample, m.Since two models have certain prediction error value on different samples, we use MRE to analyze stability and prediction accuracy of two models. the MRE of PSO-BP model is 3.58, and the MRE of MPSO-BP model is 3.09, the MRE of MPSO-BP model is smaller than the MRE of PSO-BP model; this shows that MPSO-BP model has stronger prediction ability and higher prediction accuracy than PSO-BP model. Combined with the comparison diagram of fitness evolution curves, it can be seen the training of MPSO-BP model is better than that of PSO-BP model.Although the maximum prediction error of MPSO-BP model is 17.92 on Sample 4, the maximum prediction error of PSO-BP model is 11.55 on Sample 17, the maximum prediction error of MPSO-BP model is larger than the maximum prediction error of PSO-BP model, but the stability of a prediction model cannot be determined by the error of predicting a sample alone, but by the MRE of predicting all samples. The RE of MPSO-BP model is 8.7 on Sample 4; it still has engineering significance. The prediction error of PSO-BP model is 11.32 on Sample 4, and the RE of PSO-BP model is 5.5 on Sample 4. The reason for the large prediction error of two models on this sample may be that the sample is lack of representativeness, such as the error in the collection of main influencing factors and the low measured accuracy of the HWFFZ; these will cause the sample to not reflect the objective law of the HWFFZ development. It also may be the number of samples, especially the samples with large mining thickness (here, there are only two samples with large mining thickness of more than 10 m) need to be further increased to better improve the prediction accuracy of two models. In other words, the validity of samples may have a decisive influence on the prediction accuracy of the model; in order to further ensure the effectiveness of MPSO-BP neural network model, representative samples should be selected. A sample database should be established. At the same time, it is also suggested that the geological and hydrological conditions should be closely related in the selection of main influencing factors, so as to give better play to the potential guiding significance of the prediction model.


## 3.1. Analysis of Main Influencing Factors

The formation of water flowing fractured zone of coal seam roof under the influence of mining is a complex mechanical problem. According to the Specification for Coal Pillar Retention and Coal Mining under Buildings, Water Bodies, Railways and Main Roadways and field investigation results, the HWFFZ is mainly influenced by many factors such as mining geological conditions and three-dimensional size design parameters of working face [17, 57–60]. Considering that the total caving method is widely used to manage the roof in China, this paper only studies the influence of different mining technologies on the HWFFZ under the long wall mining method. Therefore, it is determined that the influencing factors of the HWFFZ include four factors of mining geological conditions, such as mining depth, coal seam dip angle, uniaxial compressive strength of roof rock and roof strata combination structure type, and include four three-dimensional size design parameters of working face, such as mining thickness, mining method, number of mining layers, and inclined length of working face. Among them, uniaxial compressive strength of roof rock and roof strata combination structure type are replaced by the hard rock lithology proportion coefficient [6], this coefficient can reflect the overall strength of roof rock stratum and the combination structure of roof rock stratum and avoid the difficulty that the roof type cannot be determined according to the uniaxial compressive strength in the Specification for Coal Pillar Retention and Coal Mining under Buildings, Water Bodies, Railways and Main Roadways [17], and it is easy to obtain. Then influencing factors of the HWFFZ can be set as seven influencing factors including mining depth, mining thickness, hard rock lithology proportion coefficient, inclined length of working face, mining method, number of mining layers and coal seam dip angle. Among them, mining method is qualitative data, and other data are quantitative data; mining methods include blasting mining, comprehensive mechanized mining, comprehensive top coal caving mining and layered mining, with values of 0.2, 0.4, 0.6, and 0.8, respectively. By collecting and sorting 31 field measurement cases of relevant coal mines in China, the measured data of the HWFFZ are obtained. The correlation coefficient between each influencing factor of samples and the HWFFZ is calculated by GRA [61]; the broken line diagram is drawn according to the correlation coefficient, as shown in Figure 3.Figure 3
Broken line diagram of correlation coefficient.The correlation degree and correlation order between each influencing factor and the HWFFZ can be obtained, as shown in Table1.Table 1
The correlation degree and correlation order.
Mining depthMining thicknessHard rock lithology proportion coefficientInclined length of working faceMining methodNumber of mining layersCoal seam dip angle0.78930.80370.76010.79160.76040.77950.73073162547It can be seen from Figure3 and Table 1that all influencing factors have a high degree of correlation and consistency with the HWFFZ, indicating that these seven influencing factors have a certain control effect on the HWFFZ. However, it can also be seen that the correlation between mining thickness and the HWFFZ is the highest, with a correlation value of 0.8037, and the correlation between coal seam dip angle and the HWFFZ is the lowest, with a correlation value of 0.7307. According to the broken line diagram of correlation coefficient between coal seam dip angle and the HWFFZ, it can also be seen that the correlation between coal seam dip angle and the HWFFZ is the worst, and its correlation coefficient fluctuates the most. The index values of the seven influencing factors provided by samples 2 and 19 have less correlation with the HWFFZ, and the correlation has obvious consistency; however, the index values of the seven influencing factors provided by the other 29 groups of samples have great correlation with the HWFFZ, and the correlation has a certain discrete type; therefore, the GRA of measured data of the HWFFZ is effective. For further analysis, at present, the most coal mines have realized comprehensive mechanized mining and comprehensive top coal caving mining, blasting mining, and layered mining are rare, and mining thickness has a strong linear correlation with number of mining layers since the selection of main factors should be representative as much as possible and meet the requirements of linear independence between indicators; therefore, five factors including mining depth, mining thickness, hard rock lithology proportion coefficient, working face slope length, and coal seam dip angle are selected as the main influencing factors of the HWFFZ.


## 3.2. Design of Prediction Model

In total, 38 groups of measured data of the HWFFZ under the condition of comprehensive mechanized mining and comprehensive top coal caving mining of Chinese coal mines are collected as research samples; five factors including mining depth, mining thickness, hard rock lithology proportion coefficient, inclined length of working face, and coal seam dip angle are used as the input index of BP neural network, and the HWFFZ is used as the output of BP neural network. The three-layer BP network structure is adopted, and the number of hidden layer neurons is set to 6, the number of output layer neurons is set to 1. The purelin function is selected as the transfer functions of the hidden layer, the logsig function is selected as the output layer, and the trainlm function is selected as the training function. The training times are set to 1000 times, the training target error is set to 0.00001, and the learning efficiency is set to 0.04.MPSO algorithm and PSO algorithm are used to optimize BP neural network globally; according to the sum of initial weights and thresholds of BP neural network, the dimensionN of the particle is set as 43, in order to expand the search space of particles themselves. Gaussian distribution is used to perform mutation on the particles, and the threshold pm is set as 0.9. The number of particles n, the maximum number of iterations T, and the maximum limiting speed vmax are set empirically. The detailed parameter settings of the PSO algorithm and MPSO algorithm are shown in Table 2.Table 2
Parameter setting of PSO algorithm and MPSO algorithm.
AlgorithmNnTωc1c2PmvmaxPSO435050Linear decreasing22No1.9MPSO435050Exponential decreasingLinear decreasing from 2.5 to 0.5Linear increasing from 0.5 to 2.50.91.9


## 3.3. Training of Prediction Model

We use measured data of the HWFFZ to train BP neural network according to the algorithm flow in Figure2 and then establish the MPSO-BP neural network prediction model by MATLAB software. The establishment of PSO-BP neural network prediction model is similar, which is omitted here.The comparison diagram of fitness evolution curves between MPSO algorithm and PSO algorithm in the process of optimizing BP neural network is shown in Figure4.Figure 4
Comparison of fitness evolution curves.It can be seen from Figure4 that the fitness evolution curve of MPSO-BP neural network algorithm and PSO-BP neural network algorithm shows a decreasing trend. It can also be seen that two optimization algorithms have achieved the purpose of optimizing BP neural network, but MPSO algorithm is better than PSO algorithm. From the first generation to the seventh generation, the fitness value of MPSO-BP neural network algorithm evolution curve is greater than that of PSO-BP neural network algorithm evolution curve. In the seventh generation, the fitness value of MPSO-BP neural network algorithm evolution curve becomes 3.5449, while the fitness value of PSO-BP neural network algorithm evolution curve is 3.5804. Until the 10th generation, the fitness value of MPSO-BP neural network algorithm evolution curve is better than that of PSO-BP neural network algorithm evolution curve. However, from the 11th generation, the fitness evolution curve of MPSO-BP neural network algorithm presents a straight line; the value is 3.0979 and remains until the 21st generation. The fitness evolution curve of PSO-BP neural network algorithm decreases rapidly from 3.5061 to 2.8954 and shows a long straight line, no further evolution occurs until the 37th generation, which shows that PSO-BP neural network algorithm fall into local optimization and do not further evolve for a long time. The fitness evolution curve of MPSO-BP neural network algorithm rapidly decreased from 3.0978 to 2.8614 in the 20th generation, and from the 23th generation, the fitness value of evolution curve is less than that of PSO-BP neural network algorithm, and the fitness evolution curve of MPSO-BP neural network algorithm further evolves in the 43th generation, and the fitness decreases from 2.6779 to 2.3018, while the fitness evolution curve of PSO-BP neural network algorithm enters the second straight line from the 39th generation to the 46th generation, and the fitness value remains at 2.7202 and further decreases to 2.3996 in the 47th generation. It can be seen from the synthesis that fitness evolution curve of PSO-BP neural network algorithm evolves rapidly at the beginning, but it fall into local optimization from the 11th generation to the 37th generation and fails to evolve further, reflecting that the particles gathers at the same position and stops moving and refuses to search the solution space in other fields, although the fitness value is further reduced in the later stage, but PSO algorithm is weaker than MPSO algorithm in the process of global optimization. On the contrary, the MPSO-BP evolution curve keeps evolving, and there is no long-term optimization stagnation, that is, the nonevolution stage, which shows that MPSO algorithm can jump out of the local optimization and further evolve. At the end of the optimization process, its fitness value is less than that of PSO algorithm. This shows that MPSO algorithm is better than PSO algorithm, and the improved strategy makes PSO-BP further improve the performance.The comparison between measured value and prediction value of MPSO-BP neural network model and PSO-BP neural network model is shown in Figure5.Figure 5
Training comparison of the two models. (a) Comparison between measured value and prediction value of PSO-BP model. (b) Comparison between measured value and prediction value of MPSO-BP model.
(a)(b)It can be seen from Figure5(a) that the difference between prediction value of PSO-BP model and measured value is small. In 38 training samples, most of prediction value and measured value basically coincides, but there is a certain gap in individual samples, such as samples 4, 8, 10, 15, 17, and 35; its prediction error ranges from −11.37 to 11.55, but most of the prediction errors are small, and some are close to 0. It can be seen from Figure 5(b) that the difference between prediction value of MPSO-BP model and measured value is also small. In 38 training samples, most of prediction value and measured value basically coincides, and there is also a small gap in individual samples, such as samples 4, 15, 22, and 23; its prediction error ranges from −6.37 to 17.92, but most of the prediction errors are also small, and some are close to 0.In order to further analyze prediction performance of the two models, we contrast prediction stability and prediction accuracy by MRE, MRE is expressed as follows:(12)MRE=1n∗∑l=1nyl−YlYl∗100%,where yl−Yl/Yl∗100% is the relative error (RE) of the lth sample; n is the total number of training samples; yl is the prediction value of the lth sample, m; Yl is the measured value of the lth sample, m.Since two models have certain prediction error value on different samples, we use MRE to analyze stability and prediction accuracy of two models. the MRE of PSO-BP model is 3.58, and the MRE of MPSO-BP model is 3.09, the MRE of MPSO-BP model is smaller than the MRE of PSO-BP model; this shows that MPSO-BP model has stronger prediction ability and higher prediction accuracy than PSO-BP model. Combined with the comparison diagram of fitness evolution curves, it can be seen the training of MPSO-BP model is better than that of PSO-BP model.Although the maximum prediction error of MPSO-BP model is 17.92 on Sample 4, the maximum prediction error of PSO-BP model is 11.55 on Sample 17, the maximum prediction error of MPSO-BP model is larger than the maximum prediction error of PSO-BP model, but the stability of a prediction model cannot be determined by the error of predicting a sample alone, but by the MRE of predicting all samples. The RE of MPSO-BP model is 8.7 on Sample 4; it still has engineering significance. The prediction error of PSO-BP model is 11.32 on Sample 4, and the RE of PSO-BP model is 5.5 on Sample 4. The reason for the large prediction error of two models on this sample may be that the sample is lack of representativeness, such as the error in the collection of main influencing factors and the low measured accuracy of the HWFFZ; these will cause the sample to not reflect the objective law of the HWFFZ development. It also may be the number of samples, especially the samples with large mining thickness (here, there are only two samples with large mining thickness of more than 10 m) need to be further increased to better improve the prediction accuracy of two models. In other words, the validity of samples may have a decisive influence on the prediction accuracy of the model; in order to further ensure the effectiveness of MPSO-BP neural network model, representative samples should be selected. A sample database should be established. At the same time, it is also suggested that the geological and hydrological conditions should be closely related in the selection of main influencing factors, so as to give better play to the potential guiding significance of the prediction model.


## 4. Test and Analysis

### 4.1. Test of Prediction Model

Ten groups of test samples are used to test the two trained neural network models. At the same time, in order to verify whether the established model has the advantages of high prediction accuracy, empirical formulas in the Specification for Coal Pillar Retention and Coal Mining under Buildings, Water Bodies, Railways and Main Roadways are used for prediction [17]. According to the definition of hard rock lithology proportion coefficient [6], the roof of 10 groups of test samples all belongs to medium hard type. Empirical formulas under the condition of medium hard roof type are shown as follows:(13)H=100∑M1.6∑M+3.6±5.6,(14)H=20∑M+10,where H is the prediction value of the HWFFZ, m; M is the mining thickness, m.The comparison between the prediction results of the two neural network models and empirical formulas and the measured values is shown in Figure6.Figure 6
Test comparison of between measured value and prediction value.It can be seen from Figure6 that the most prediction results of two neural network models are close to measured values, and empirical formulas have obvious prediction errors for samples 1, 5, 6, and 10. Especially for Sample 6, mining depth of working face is 590 m and mining thickness is 9.97 m; although the cumulative mining thickness limited by the Specification for Coal Pillar Retention and Coal Mining under Buildings, Water Bodies, Railways and Main Roadways is no more than 15 m, empirical formulas can be used to calculate the HWFFZ. It is found that the REs of empirical formulas are 71.56 and 63.24, respectively, both of them are large; therefore, the prediction results have lost their guiding significance for production. However, two neural network models can predict with the high precision. The RE of MPSO-BP model is 5.47, and the RE of PSO-BP model is 6.22. The prediction result of MPSO-BP model is the closest to measured value. In order to further analyze the stability and generalization of these models, the MRE is used for statistical analysis. The MRE of MPSO-BP model is 3.53, MRE of PSO-BP model is 4.3, the MRE of equation (13) is 20.71, and the MRE of equation (14) is 21.23; it can be seen that MPSO-BP model has the best stability and generalization.The prediction accuracies of empirical formulas are low, which are related to the compilation background of the Specification for Coal Pillar Retention and Coal Mining under Buildings, Water Bodies, Railways and Main Roadways. Measured data based on the Specification for Coal Pillar Retention and Coal Mining under Buildings, Water Bodies, Railways and Main Roadways mainly came from measured values of blasting mining, general mining and layered mining in the 1950s and 1980s, and mining depth is small, mostly about 300 m, generally not more than 500 m. Limited by the production technical conditions at that time, the thickness of coal seam stratified mining in the statistical samples is generally less than 3 m. Since the 1980s, great changes have taken place in coal mining methods, mining depth and inclined length of working face. In particular, comprehensive mechanized mining developed after the 1980s, including the later comprehensive top coal caving mining, have been widely used in China. In addition, now the shallow coal resources are gradually mined out, gradually moving towards deep mining, at present, mining depth has reached more than 1000 m. In addition, with the expansion of mining scope, Chinese coal mining base has gradually shifted from the East to the West, resulting in changes in mining geological conditions. Based on the above analysis, empirical formulas for calculating the HWFFZ are obtained on the basis of a large number of measured data, which meet the requirements of coal mining design under water bodies to a certain extent in China. This can be seen from the REs of predicting the HWFFZ on samples 2, 3, 4, 7, and 9. The RE of equation (13) on Sample 3 is 0, and the RE of equation (14) on Sample 2 is only 1.96. Mining depth of Sample 2 is 264.5 m, and mining thickness is 2.8 m. Mining depth of Sample 3 is 290 m, and mining thickness is 2.6 m. These are consistent with the most measured data background used to fit empirical formulas. Therefore, empirical formulas still have certain guiding significance to a certain extent. However, empirical formulas can only be used when the mining depth is 300 m, and the mining thickness is about 3 m. With the change of mining methods, the increase of mining depth and the change of mining geological conditions, empirical formulas cannot meet the requirements of coal mine safety production; therefore, the MPSO-BP model proposed in this paper has a certain practical significance, and it is also an effective means for predicting the HWFFZ.Here, we put forward further requirements for improving the prediction accuracy of MPSO-BP model. Empirical formulas can occasionally get high-precision prediction values based on only mining thickness, which is mainly due to the collection of representative samples and a large number of samples, so as to fit empirical formulas reflecting the objective law of the HWFFZ development. In order to improve the accuracy of the MPSO-BP model, it is necessary to improve the representativeness of the collected samples and try to establish a rich sample base. The neural network has a strong nonlinear mapping ability and can establish a high-precision prediction model. It can be seen mining thicknesses of Samples 2 and 7 are the same, which leads to the same prediction results of empirical formulas. The mining thicknesses of Samples 3 and 4 are the same, which leads to the same prediction results of empirical formulas. However, the two neural network models can predict the HWFFZ with high accuracy based on more factors, which is not available in empirical formulas. Therefore, it is necessary to further study the influencing factors of the HWFFZ and consider the main factors of the HWFFZ as comprehensively as possible, so as to provide a better basis and foundation for establishing prediction model of the HWFFZ, then a more accurate MPSO-BP prediction model can be established, which is also the direction of further research and efforts.


### 4.2. Case Application

The overall structural form of 121101 working face in Huainan coalfield is a monoclinal structure, the mining method of working face is comprehensive mechanized mining, the coal seam trend is near east-west and inclines to the south, the average thickness of coal seam is 3.63 m, and the average coal seam dip angle is 14°, the strike length of the transportation roadway of working face is 1481 m, the floor elevation is −479.4 m∼−514.5 m, the strike length of the air return roadway of working face is 1522.1 m, the floor elevation is −529.5 m∼−565.4 m, the inclined length of working face is 192.8 m, and the corresponding ground elevation of working face is +24.9 m∼+25.6 m. The roof aquifer of 11-2 Coal Seam in the mine is a sandstone aquifer; the thickness, fracture development degree, and water bearing capacity of the sandstone layer vary with different regions. The total thickness of sandstone within the estimated HWFFZ range of working face roof is 29.8 m∼45.5 m. The water in sandstone aquifer may suddenly gush out during mining process and cause mine flooding. In order to grasp the HWFFZ of working face, effectively guide the prediction and prevention of mine roof water and ensure the mining safety of working face, it is needed to carry out the research on the prediction of the HWFFZ.According to the overview of 121101 working face, five main parameters of prediction model of the HWFFZ are extracted. The location of the measured drilling site arranged underground is taken as the basis for mining depth selection, which is more targeted; according to the floor elevation and ground elevation of the drilling site, mining depth near the drilling site is determined to be 535 m. According to the development of coal seams near the drilling site, mining thickness is taken as 3.63 m. According to the comprehensive histogram of No. 1 borehole in the drilling site of working face, hard rock lithology proportion coefficient is calculated to be 0.33. Inclined length of working face is taken as 192.8 m. According to the coal seam development near the drilling site, coal seam dip angle is taken as 14°. The above parameter values are input into the MPSO-BP model for output prediction of the HWFFZ, and the prediction result is 50.53 m. The prediction value of equation (13) is 44.18 m, and the prediction value of equation (14) is 48.11 m. The HWFFZ of working face was measured with drilling double end plugging leak detection device [62], and the measured value is 53.6 m. The RE of MPSO-BP model is 5.73, the RE of equation (13) is 17.57, and the RE of equation (14) is 10.24; therefore, it can be seen that the MPSO-BP model has the highest accuracy, and it has a certain reliability for guiding production.


## 4.1. Test of Prediction Model

Ten groups of test samples are used to test the two trained neural network models. At the same time, in order to verify whether the established model has the advantages of high prediction accuracy, empirical formulas in the Specification for Coal Pillar Retention and Coal Mining under Buildings, Water Bodies, Railways and Main Roadways are used for prediction [17]. According to the definition of hard rock lithology proportion coefficient [6], the roof of 10 groups of test samples all belongs to medium hard type. Empirical formulas under the condition of medium hard roof type are shown as follows:(13)H=100∑M1.6∑M+3.6±5.6,(14)H=20∑M+10,where H is the prediction value of the HWFFZ, m; M is the mining thickness, m.The comparison between the prediction results of the two neural network models and empirical formulas and the measured values is shown in Figure6.Figure 6
Test comparison of between measured value and prediction value.It can be seen from Figure6 that the most prediction results of two neural network models are close to measured values, and empirical formulas have obvious prediction errors for samples 1, 5, 6, and 10. Especially for Sample 6, mining depth of working face is 590 m and mining thickness is 9.97 m; although the cumulative mining thickness limited by the Specification for Coal Pillar Retention and Coal Mining under Buildings, Water Bodies, Railways and Main Roadways is no more than 15 m, empirical formulas can be used to calculate the HWFFZ. It is found that the REs of empirical formulas are 71.56 and 63.24, respectively, both of them are large; therefore, the prediction results have lost their guiding significance for production. However, two neural network models can predict with the high precision. The RE of MPSO-BP model is 5.47, and the RE of PSO-BP model is 6.22. The prediction result of MPSO-BP model is the closest to measured value. In order to further analyze the stability and generalization of these models, the MRE is used for statistical analysis. The MRE of MPSO-BP model is 3.53, MRE of PSO-BP model is 4.3, the MRE of equation (13) is 20.71, and the MRE of equation (14) is 21.23; it can be seen that MPSO-BP model has the best stability and generalization.The prediction accuracies of empirical formulas are low, which are related to the compilation background of the Specification for Coal Pillar Retention and Coal Mining under Buildings, Water Bodies, Railways and Main Roadways. Measured data based on the Specification for Coal Pillar Retention and Coal Mining under Buildings, Water Bodies, Railways and Main Roadways mainly came from measured values of blasting mining, general mining and layered mining in the 1950s and 1980s, and mining depth is small, mostly about 300 m, generally not more than 500 m. Limited by the production technical conditions at that time, the thickness of coal seam stratified mining in the statistical samples is generally less than 3 m. Since the 1980s, great changes have taken place in coal mining methods, mining depth and inclined length of working face. In particular, comprehensive mechanized mining developed after the 1980s, including the later comprehensive top coal caving mining, have been widely used in China. In addition, now the shallow coal resources are gradually mined out, gradually moving towards deep mining, at present, mining depth has reached more than 1000 m. In addition, with the expansion of mining scope, Chinese coal mining base has gradually shifted from the East to the West, resulting in changes in mining geological conditions. Based on the above analysis, empirical formulas for calculating the HWFFZ are obtained on the basis of a large number of measured data, which meet the requirements of coal mining design under water bodies to a certain extent in China. This can be seen from the REs of predicting the HWFFZ on samples 2, 3, 4, 7, and 9. The RE of equation (13) on Sample 3 is 0, and the RE of equation (14) on Sample 2 is only 1.96. Mining depth of Sample 2 is 264.5 m, and mining thickness is 2.8 m. Mining depth of Sample 3 is 290 m, and mining thickness is 2.6 m. These are consistent with the most measured data background used to fit empirical formulas. Therefore, empirical formulas still have certain guiding significance to a certain extent. However, empirical formulas can only be used when the mining depth is 300 m, and the mining thickness is about 3 m. With the change of mining methods, the increase of mining depth and the change of mining geological conditions, empirical formulas cannot meet the requirements of coal mine safety production; therefore, the MPSO-BP model proposed in this paper has a certain practical significance, and it is also an effective means for predicting the HWFFZ.Here, we put forward further requirements for improving the prediction accuracy of MPSO-BP model. Empirical formulas can occasionally get high-precision prediction values based on only mining thickness, which is mainly due to the collection of representative samples and a large number of samples, so as to fit empirical formulas reflecting the objective law of the HWFFZ development. In order to improve the accuracy of the MPSO-BP model, it is necessary to improve the representativeness of the collected samples and try to establish a rich sample base. The neural network has a strong nonlinear mapping ability and can establish a high-precision prediction model. It can be seen mining thicknesses of Samples 2 and 7 are the same, which leads to the same prediction results of empirical formulas. The mining thicknesses of Samples 3 and 4 are the same, which leads to the same prediction results of empirical formulas. However, the two neural network models can predict the HWFFZ with high accuracy based on more factors, which is not available in empirical formulas. Therefore, it is necessary to further study the influencing factors of the HWFFZ and consider the main factors of the HWFFZ as comprehensively as possible, so as to provide a better basis and foundation for establishing prediction model of the HWFFZ, then a more accurate MPSO-BP prediction model can be established, which is also the direction of further research and efforts.


## 4.2. Case Application

The overall structural form of 121101 working face in Huainan coalfield is a monoclinal structure, the mining method of working face is comprehensive mechanized mining, the coal seam trend is near east-west and inclines to the south, the average thickness of coal seam is 3.63 m, and the average coal seam dip angle is 14°, the strike length of the transportation roadway of working face is 1481 m, the floor elevation is −479.4 m∼−514.5 m, the strike length of the air return roadway of working face is 1522.1 m, the floor elevation is −529.5 m∼−565.4 m, the inclined length of working face is 192.8 m, and the corresponding ground elevation of working face is +24.9 m∼+25.6 m. The roof aquifer of 11-2 Coal Seam in the mine is a sandstone aquifer; the thickness, fracture development degree, and water bearing capacity of the sandstone layer vary with different regions. The total thickness of sandstone within the estimated HWFFZ range of working face roof is 29.8 m∼45.5 m. The water in sandstone aquifer may suddenly gush out during mining process and cause mine flooding. In order to grasp the HWFFZ of working face, effectively guide the prediction and prevention of mine roof water and ensure the mining safety of working face, it is needed to carry out the research on the prediction of the HWFFZ.According to the overview of 121101 working face, five main parameters of prediction model of the HWFFZ are extracted. The location of the measured drilling site arranged underground is taken as the basis for mining depth selection, which is more targeted; according to the floor elevation and ground elevation of the drilling site, mining depth near the drilling site is determined to be 535 m. According to the development of coal seams near the drilling site, mining thickness is taken as 3.63 m. According to the comprehensive histogram of No. 1 borehole in the drilling site of working face, hard rock lithology proportion coefficient is calculated to be 0.33. Inclined length of working face is taken as 192.8 m. According to the coal seam development near the drilling site, coal seam dip angle is taken as 14°. The above parameter values are input into the MPSO-BP model for output prediction of the HWFFZ, and the prediction result is 50.53 m. The prediction value of equation (13) is 44.18 m, and the prediction value of equation (14) is 48.11 m. The HWFFZ of working face was measured with drilling double end plugging leak detection device [62], and the measured value is 53.6 m. The RE of MPSO-BP model is 5.73, the RE of equation (13) is 17.57, and the RE of equation (14) is 10.24; therefore, it can be seen that the MPSO-BP model has the highest accuracy, and it has a certain reliability for guiding production.


## 5. Conclusions

In order to improve the accuracy of PSO-BP neural network prediction model of the HWFFZ and effectively prevent water inrush accident of mine roof, an improved PSO algorithm is proposed, and the MPSO-BP neural network algorithm is established. According to measured data of the HWFFZ, the GRA is used to analyze the main influencing factors of the HWFFZ. The five main influencing factors screened can be used as the input of BP neural network. This analysis not only effectively reduces the number of input indicators but also provides a basis for improving the output accuracy of the MPSO-BP neural network prediction model. Based on the MPSO-BP neural network algorithm, the nonlinear relations between five main influencing factors and the HWFFZ are established. Through training and test, it can be seen that the MPSO-BP neural network prediction model has the highest prediction stability and prediction accuracy than the PSO-BP neural network prediction model and empirical formulas in the Specification for Coal Pillar Retention and Coal Mining under Buildings, Water Bodies, Railways and Main Roadways. Application results also show that MPSO-BP neural network prediction model can meet the needs of engineering and has certain guiding significance for production. This paper provides a new method for predicting the HWFFZ.



---

*Source: 2133695-2022-08-22.xml*