Comprehensive Prediction and Discriminant Model for Rockburst Intensity Based on Improved Variable Fuzzy Sets Approach

: Rockburst intensity prediction is one of the basic works of underground engineering disaster prevention and mitigation. Considering the dynamic variability and fuzziness in rockburst intensity prediction, variable fuzzy sets (VFS) are selected for evaluation and prediction. Here, there are two problems in the application of traditional VFS: (i) the relative membership degree (RMD) calculation process is complex and time-consuming, and the RMD matrix of all indexes can be only obtained by using the RMD function repeatedly; (ii) unreasonable weights of indicators have great impact on the synthetic relative membership degree (SRMD), so it is di ﬃ cult to guarantee the correctness of the ﬁnal prediction result. In view of the above problem, this paper established three simpliﬁed feature relationship expressions of RMD based on VFS principle and used the SRMD function to establish a BP neural network model to optimize SRMD. The improved VFS method is more e ﬃ cient and the prediction results are more stable and reliable than the traditional VFS method. The main advantages are as follows: (1) the improved VFS method has higher computational e ﬃ ciency; (2) the improved VFS method can verify the correctness of RMD at all times; (3) the improved VFS method has higher prediction accuracy; and (4) the improved VFS method has higher fault tolerance and practicability.


Introduction
The frequent occurrence of rockburst disasters has been acknowledged as one of the most serious problems in underground projects all over the world because it directly threatened the safety of underground constructors, equipment, and buildings and even induced mine earthquakes [1]. In recent decades, many reliable measures have been studied to prevent rockburst. For example, Zheng et al. [2] applied borehole pressure relief technology to prevent rockburst in deep-buried mine roadways and achieved good results. Skrzypkowski et al. [3,4] studied a new type of bolt support system, which had good energy absorption and rockburst prevention ability under the condition of tremor. However, in addition to finding reliable prevention and control measures, it is also important to establish an accurate rockburst prediction model in order to resist the occurrence of rockburst disasters [5][6][7][8]. Presently, many scholars have studied the mechanism of rockburst from different perspectives and put forward corresponding prediction methods of rockburst intensity, for example, based on single factor strength theory, rigidity theory, energy theory, catastrophe theory, bifurcation theory, instability theory, Russeenes criterion [7], Wang Yuanhan criterion [8] and Lu Jiayou criterion [9], and random forest classification [10], cloud model theory [11], attribute comprehensive evaluation method [12], artificial neural network [13], matter-element extension theory [14], etc. All of the above methods have predicted rockburst from different perspectives and achieved certain prediction results. However, where D A (x) is the relative difference function of x to A. Hence, the RMD can be denoted as follows:

RMD Function
In VFS, attraction set and extended set are two important concepts to establish RMD function. Let X 0 = [a, b] be the attracting set of x on the real axis (0 < D A (x) ≤ 1) and X = [c, d] be the extended set containing X 0 (X 0 ∈ X). The M represents equilibrium point D A (x) = 1 in the attracting set [a, b]. From the principle of VFS, it is clear that the intervals [c, a] and [b, d] are repelling sets of x (−1 < D A (x) ≤ 0) ( Figure 1).

Principle of VFS
Definition: Let U be the universe of discourse, A be a fuzzy concept in U, and be a random element, ∈U. The RMDs of to the concept of attracting (A) and repelling property (A ) are ( ) and ( ), respectively. There are the following equations among them:

RMD Function
In VFS, attraction set and extended set are two important concepts to establish RMD function. Let 0 = [ , ] be the attracting set of on the real axis (0 < ( ) ≤ 1 ) and = [ , ] be the extended set containing 0 ( 0 ∈ ). The M represents equilibrium point ( ) = 1 in the attracting set [ , ]. From the principle of VFS, it is clear that the intervals [ , ] and [ , ] are repelling sets of (−1 < ( ) ≤ 0) ( Figure 1).
When is on the right side of M, the relative difference function is established as follows: According to Equations (2) and (3), the RMD function is expressed as follows: According to Equations (2) and (4), the RMD function is expressed as follows: ; when x is on the left side of M, the relative difference function is established as follows: When x is on the right side of M, the relative difference function is established as follows: According to Equations (2) and (3), the RMD function is expressed as follows: Appl. Sci. 2019, 9, 3173 4 of 19 According to Equations (2) and (4), the RMD function is expressed as follows: u A (x) = 0.5 In Equations (5)-(8), β is a nonnegative index and, generally, β = 1. Equations (5)-(8) meet the following conditions:

SRMD Function
In practical applications, the prediction and evaluation objects often have multiple attributes. In order to make the results more realistic, Chen [15] proposed a comprehensive evaluation model to calculate SRMD. The formula is as follows: In Equation (10), n and m are the number of indicators and classification grade, respectively; u ij (x) is the RMD of eigenvalue x of the ith indicator to level j; α and p are the model optimization criterion parameter and the distance parameter, respectively (α = 1 for least absolute criteria and α = 2 for least squares criterion, p = 1 denotes Heming distance and p = 2 denotes Euclidean distance); and w i is the weight of indicator i, n i=1 w i = 1.

Calculating Steps of VFS Method
Step 1: Establishing parameters (c, a, M, b, d) of VFS model Assuming an indicator i is divided into m levels and the boundary eigenvalues of each level are [L i1 , L i2 ], [L i2 , L i3 ], · · · , L im , L i(m+1) (Figure 2). According to previous literatures [13,22,25], the attracting interval a ij , b ij , extended interval c ij , d ij , and point M ij for each level are calculated as follows: (i) For the minimum level, the intervals [a i1 , b i1 ] and [c i1 , d i1 ] and point M i1 were calculated as follows: (ii) For level j (j = 2, 3, · · · , m − 1), the intervals a ij , b ij and c ij , d ij and point M ij were expressed as follows: Appl. Sci. 2019, 9, x FOR PEER REVIEW 5 of 20 Step 2: Calculating RMD Substituting the above parameters (c, a, M, b, d) into Equations (5)- (8) repeatedly can obtain RMD functions corresponding to each level. Then, the measured value of index i is brought into the RMD function to calculate the RMD value of each level, which is expressed as follows: Now, assuming that the prediction and assessment object has n indicators, the RMD matrix of n indicators to m levels can be obtained by repeating the above process n times, which can be expressed as follows: Step 3: Calculating SRMD Based on the obtained RMD matrix, the SRMD of the prediction and assessment object to each level can be obtained by introducing ( ) into Equation (10): After calculating the SRMD, the final prediction and assessment level eigenvalue H can be obtained: The relationship between classification level and eigenvalue H is expressed as follows: In the above calculation process, it can be clearly found that two issues need to be improved when employing the VFS method for multi-index prediction and evaluation. Firstly, how to simplify RMD calculation process? The traditional RMD calculation process is very complicated and timeconsuming; a prediction object needs to run n × m times Equations (5)- (8) repeatedly to get the RMD matrix. Although Equation (9) can reduce some, it still has a large computational burden. Secondly, how to consider the indicator weights of SRMD function so as to obtain high-precision prediction results? In Equation (10), weight assignment is the core issue, which is related to the prediction accuracy of SRMD. However, the weighting is generally carried out in a relatively specific situation in the actual forecasting process. As shown in Table 1, weights are assigned differently based on different occasions, objects, and methods. Weights applicable to one occasion cannot necessarily Step 2: Calculating RMD Substituting the above parameters (c, a, M, b, d) into Equations (5)- (8) repeatedly can obtain RMD functions corresponding to each level. Then, the measured value x of index i is brought into the RMD function to calculate the RMD value of each level, which is expressed as follows: Now, assuming that the prediction and assessment object has n indicators, the RMD matrix of n indicators to m levels can be obtained by repeating the above process n times, which can be expressed as follows: Step 3: Calculating SRMD Based on the obtained RMD matrix, the SRMD of the prediction and assessment object to each level can be obtained by introducing u ij (x) into Equation (10): After calculating the SRMD, the final prediction and assessment level eigenvalue H can be obtained: The relationship between classification level and eigenvalue H is expressed as follows: In the above calculation process, it can be clearly found that two issues need to be improved when employing the VFS method for multi-index prediction and evaluation. Firstly, how to simplify RMD calculation process? The traditional RMD calculation process is very complicated and time-consuming; a prediction object needs to run n × m times Equations (5)- (8) repeatedly to get the RMD matrix. Although Equation (9) can reduce some, it still has a large computational burden. Secondly, how to consider the indicator weights of SRMD function so as to obtain high-precision prediction results? In Equation (10), weight assignment is the core issue, which is related to the prediction accuracy of SRMD. However, the weighting is generally carried out in a relatively specific situation in the actual forecasting process. As shown in Table 1, weights are assigned differently based on different occasions, objects, and methods. Weights applicable to one occasion cannot necessarily apply to other occasions, which hinders the promotion and application of prediction methods. Therefore, it is necessary to find a suitable method to optimize SRMD. Only in this way can we effectively guarantee the prediction accuracy of SRMD, ignore the core role of weight, and improve the fault tolerance and application ability of the method.

Simplifying the RMD Calculation Process
In the actual calculation process, combining Equations (5)-(8) and (11)- (13), it can be found that, when the measured value x of index i is located in different levels, the RMDs exhibit some general characteristics. Detailed analysis is as follows: (1) When x is in the minimum level interval, L i1 ≤ x ≤ L i2 (Figure 3a).
When L i1 ≤ x ≤ L i2 , only RMD functions of level 1 and level 2 are activated; for other levels, RMDs are equal to 0 because x is not in the [c, d] range. According to Equations (11) and (12), the parameters of the VFS model with levels 1 and 2 can be obtained as follows: Thus, a linear relationship between x and parameters c, a, M, b, and d for levels 1 and 2 can be established as shown in Figure 3b. From Figure 3b, it can be seen that x is in the [M i1 , b i1 ] interval for level 1, and for level 2, x is in the [c i2 , a i2 ] interval. Therefore, the RMD of x for levels 1 and 2 can be calculated by Equations (6) and (7) as follows: In Equation (19) In addition, it can be easily found that u i1 (x) + u i2 (x) = 1. The above features can also be presented intuitively in Figure 3c. Based on the above analysis, when the measured value x of index i is within level 1, the characteristics of RMD can be summarized as follows: Appl. Sci. 2019, 9, x FOR PEER REVIEW 7 of 20 (2) When is located in the interval of a mid-level j (j = 2,3, ⋯ , − 1), ≤ ≤ ( +1) ( Figure  4a). When ≤ ≤ ( +1) , only RMD functions of level j and adjacent levels j − 1 and j + 1 are activated; for other levels, RMDs are equal to 0 because is not in the [ , ] range. According to Equation (12), the parameters of the VFS model with levels j, j -1, and j + 1 can be obtained as follows: Level j + 1  (5) and (7), the RMD of x for level j can be calculated as follows: (2) When x is located in the interval of a mid-level j only RMD functions of level j and adjacent levels j − 1 and j + 1 are activated; for other levels, RMDs are equal to 0 because x is not in the [c, d] range. According to Equation (12), the parameters of the VFS model with levels j, j − 1, and j + 1 can be obtained as follows: Level j + 1 Based on Equations (21)-(23), a linear relationship between x and c, a, M, b, and d for levels j, j − 1, and j + 1 can be established as shown in Figure 4b. From Figure 4b, it can be seen that there are two cases for level j when x is in the a ij , b ij interval: (1) a ij ≤ x ≤ M ij and (2) M ij ≤ x ≤ b ij . Combined with Equations (5) and (7), the RMD of x for level j can be calculated as follows: In Equation (24), 0 ≤ 2    (6) and (8): In Equation (25) The above features can also be clearly seen in Figure 4c, which can be summarized as follows: , respectively, and RMD can be calculated according to Equations (6) and (8): In Equation (25) The above features can also be clearly seen in Figure 4c, which can be summarized as follows:  (12) and (13), the parameters of the VFS model for levels m and m − 1 can be obtained as follows: . Thus, a linear relationship between x and c, a, M, b, and d for levels m and m − 1 can be established as shown in Figure 5b. Figure 5b shows that x is in the [a im , M im ] interval for level m and that, for level m − 1, x is in the b i(m−1) , d i(m−1) interval. Therefore, the RMD of x for levels m and m − 1 can be calculated by Equations (5) and (8) as follows: In Equation (27) these features can also be visually found in Figure 5c, summarized as follows: In Equation (27) In the above analysis, only considering the change of measured value between the interval [ 1 , ( +1) ], the situation of ≤ 1 or ≥ ( +1) is not considered because the cognition is exact and not fuzzy in these two cases; its corresponding RMD value is equal to 1, which is expressed as follows: According to the above characteristics (Equations (20), (26), (28), and (29)), combined with Equations (5)-(8), a simplified method of RMD calculation is proposed. The detailed steps are as follows: (1) The VFS model parameters (c, a, M, b, d) corresponding to each level are established according to the criteria of indicator division, and the position of measured value x in the level interval needs to be identified. (2) According to the level interval where x is located, the corresponding activation RMD function is determined. If x is in the minimum or maximum levels, it only needs to run the RMD function once. For level 1, RMD can be calculated by Equation (7), and for the maximum level, RMD can be obtained by Equation (5). If x is in the middle level, only the RMD corresponding to the level itself and one of the adjacent levels need to be calculated for RMD of the level itself, which can In the above analysis, only considering the change of measured value x between the interval L i1 , L i(m+1) , the situation of x ≤ L i1 or x ≥ L i(m+1) is not considered because the cognition is exact and not fuzzy in these two cases; its corresponding RMD value is equal to 1, which is expressed as follows: According to the above characteristics (Equations (20), (26), (28), and (29)), combined with Equations (5)-(8), a simplified method of RMD calculation is proposed. The detailed steps are as follows: (1) The VFS model parameters (c, a, M, b, d) corresponding to each level are established according to the criteria of indicator division, and the position of measured value x in the level interval needs to be identified. (2) According to the level interval where x is located, the corresponding activation RMD function is determined. If x is in the minimum or maximum levels, it only needs to run the RMD function once. For level 1, RMD can be calculated by Equation (7), and for the maximum level, RMD can be obtained by Equation (5). If x is in the middle level, only the RMD corresponding to the level itself and one of the adjacent levels need to be calculated for RMD of the level itself, which can be calculated by selecting an appropriate function from Equations (5) and (7), while for the RMD of the adjacent level, one of Equations (6) and (8) can be selected for calculation. The RMD corresponding to other levels can be obtained by feature Equations (20), (26), and (28). In addition, when x is greater than L i(m+1) or less than L i1 , the RMD function does not need to be run and the RMD is determined directly by Equation (29). (3) By repeating the above steps, the RMD values for other indicators can be calculated and the final RMD matrix can be obtained.

Optimizing SRMD
In Equation (10), a weight is determined according to the actual situation. If the research target already has a reasonable weight assignment, it is taken as the initialization weight. If it does not exist, equal weight is preferable. Based on this, the obtained RMD is brought into Equation (10) to calculate an initial SRMD. Then, a more accurate result is obtained by optimizing the initial SRMD using the BP network model. In this way, the second issue raised in Section 2.3 can be effectively solved. Many studies have shown that the fusion of BP neural network and fuzzy set theory can give better play to their respective advantages [30,31]. In this paper, we use the sigmoid function characteristics of the SRMD function in Equation (10) to establish a BP neural network to optimize the SRMD. The establishment process mainly includes three parts: (1) creating a network structure; (2) setting the activation function and calculating input and output of each layer; and (3) creating a weight adjustment model. The specific establishment process is as follows: (1) Creating a network structure For convenient discussion, a three-layer structure of BP neural network is constructed, as shown in Figure 6. Let the input layer have m nodes, i.e., the initial SRMD corresponds to m levels; the hidden layer have p nodes, i.e., p unit systems; and the output layer have also m nodes corresponding to the optimized SRMD. The error back propagation algorithm of BP network is used to update the connection weights (w ik and w kh ) between network layers.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 10 of 20 be calculated by selecting an appropriate function from Equations (5) and (7), while for the RMD of the adjacent level, one of Equations (6) and (8) can be selected for calculation. The RMD corresponding to other levels can be obtained by feature Equations (20), (26), and (28). In addition, when x is greater than ( +1) or less than 1 , the RMD function does not need to be run and the RMD is determined directly by Equation (29). (3) By repeating the above steps, the RMD values for other indicators can be calculated and the final RMD matrix can be obtained.

Optimizing SRMD
In Equation (10), a weight is determined according to the actual situation. If the research target already has a reasonable weight assignment, it is taken as the initialization weight. If it does not exist, equal weight is preferable. Based on this, the obtained RMD is brought into Equation (10) to calculate an initial SRMD. Then, a more accurate result is obtained by optimizing the initial SRMD using the BP network model. In this way, the second issue raised in Section 2.3 can be effectively solved. Many studies have shown that the fusion of BP neural network and fuzzy set theory can give better play to their respective advantages [30,31]. In this paper, we use the sigmoid function characteristics of the SRMD function in Equation (10) to establish a BP neural network to optimize the SRMD. The establishment process mainly includes three parts: (1) creating a network structure; (2) setting the activation function and calculating input and output of each layer; and (3) creating a weight adjustment model. The specific establishment process is as follows: (1) Creating a network structure For convenient discussion, a three-layer structure of BP neural network is constructed, as shown in Figure 6. Let the input layer have m nodes, i.e., the initial SRMD corresponds to m levels; the hidden layer have p nodes, i.e., p unit systems; and the output layer have also m nodes corresponding to the optimized SRMD. The error back propagation algorithm of BP network is used to update the connection weights ( and ℎ ) between network layers. (2) Setting the activation function Sigmoid function is a commonly used activation function of BP networks, which has the characteristics of being s-shaped, asymptotically bounded, completely monotone [32,33]. When α = 2, p = 1, let = ∑ [ (1 − ( ))] =1 and = ∑ =1 ( ) such that = 1 − . Submit into Equation (10) to obtain the following: (2) Setting the activation function Sigmoid function is a commonly used activation function of BP networks, which has the characteristics of being s-shaped, asymptotically bounded, completely monotone [32,33]. When α = 2, p = 1, let d jg = n i=1 w i 1 − u ij (x) and d jb = n i=1 w i u ij (x) such that d jg = 1 − d ib . Submit d jb into Equation (10) to obtain the following: According to Equation (30), v(x) j is a nonlinear function of d jb ; its first and second derivatives can be expressed as follows: dv In Equation (31), 0 ≤ d jb ≤ 1; therefore dv(x) j dd jb > 0 and then v(x) j is a monotonically increasing function of d jb .
According to Equation (32), when d jb = 0.5, According to Equation (30), ( ) is a nonlinear function of ; its first and second derivatives can be expressed as follows: In Equation (31)  (3) Creating a weight adjustment model Before creating a weight updating formula, the activation function (Equation (30)) is used to calculate the input and output of hidden and output layers. In the input layer, the input information (SRMD for each level) of node m is transmitted directly to the hidden layer node k, that is, the input and output of the node are equal ( = ( ) ). For the hidden layer node k, the input and output are as follows: is the connection weight between the input layer and the hidden layer. For the output layer node h, the input and output are as follows: (3) Creating a weight adjustment model Before creating a weight updating formula, the activation function (Equation (30)) is used to calculate the input and output of hidden and output layers. In the input layer, the input information (SRMD for each level) of node m is transmitted directly to the hidden layer node k, that is, the input and output of the node are equal (x j = v(x) j ). For the hidden layer node k, the input and output are as follows: w jk is the connection weight between the input layer and the hidden layer. For the output layer node h, the input and output are as follows: w kh is the connection weight between the hidden and output layers; u h is the output value of the output layer of the BP network model, which is the optimized SRMD for level h.
Assuming that the expected output of the evaluation or prediction object is o h , the square error E between the expected output o h and the actual output u h is as follows: According to Equation (35), the weight-updating formula between input layer node j and hidden layer node k can be deduced as follows [15]: In Equation (36), γ is the learning ratio of the BP network. The weight adjustment value ∆w kh between the hidden layer node k and the output layer node h is as follows: Equations (36) and (37) are the weight adjustment models of the BP network. Applying this model and combining with the iteration algorithm of normal neural network, the connection weight can be determined to minimize the error between the actual output and the expected output.

Framework for Assessment and Prediction
Based on the above analysis, the framework of the assessment and prediction method based on the improved VFS is summarized as follows: • Determine the impact indicators of the forecasting and evaluation object, determine the classification criteria of impact indicators, and divide them into m levels.

Rockburst Prediction Indicators and Cases
The occurrence mechanism of rockburst is complex and there are many influencing factors. The selection of criteria is the key step in the prediction process. Based on the related studies of rockburst [8,34,35], this paper synthetically considered the factors of in situ stress, lithology, energy, and geological structure and selected the ratio of tangential stress to uniaxial compressive strength ⁄ (stress coefficient), the ratio of uniaxial compressive strength to tensile strength ⁄ (brittleness coefficient), and elastic energy index as indicators of rockburst proneness prediction. From the current research situation, most of the literature has selected similar prediction indicators because

Rockburst Prediction Indicators and Cases
The occurrence mechanism of rockburst is complex and there are many influencing factors. The selection of criteria is the key step in the prediction process. Based on the related studies of rockburst [8,34,35], this paper synthetically considered the factors of in situ stress, lithology, energy, and geological structure and selected the ratio of tangential stress to uniaxial compressive strength σ θ /σ c (stress coefficient), the ratio of uniaxial compressive strength to tensile strength σ c /σ t (brittleness coefficient), and elastic energy index w et as indicators of rockburst proneness prediction. From the current research situation, most of the literature has selected similar prediction indicators because they cover both internal and external factors affecting the occurrence of rockburst and are easy to obtain in the laboratory or field [36][37][38]. Also, it is convenient for comparative analysis between different rockburst engineering cases. According to Wang et al. [8], rockburst prediction criteria for these three indicators are established as shown in Table 2. In addition, in order to verify the correctness and effectiveness of the improved VFS method in rockburst prediction, this paper takes 18 large underground engineering rockburst cases collected by Wang et al. [8] as examples to analyze. The basic data is shown in Table 3.

The Calculation Process of the Improved VFS Method
(1) Determining attracting interval, extended interval, and M point Before calculating the parameters of the VFS model, it should be pointed out that Equations (11)-(13) are available for fixed intervals, while levels I and IV in rockburst prediction criteria are semi-interval forms such as [−∞, L max ] and [L min , +∞]. Here, a pseudo-boundary is obtained using the polynomial regression analysis method, assuming that the trend of increasing or decreasing of boundary eigenvalues is consistent with level L and can be recognized (Figure 9). Based on rockburst classification criteria and the polynomial regression technique, all quantitative boundaries of the levels for all indicators can be obtained and the VFS model parameters a ij , b ij , c ij , d ij , and M ij for all level can be determined as follows: where the first row is the attracting interval of σ θ /σ c to four levels and the remaining rows represent σ c /σ t and w et .
where the first row is the points M of σ θ /σ c to four levels and the remaining rows represent σ c /σ t and w et .
where the first row is the extended interval of σ θ /σ c to four levels and the remaining rows represent σ c /σ t and w et . where the first row is the extended interval of ⁄ to four levels and the remaining rows represent ⁄ and .  As for w et , the eigenvalue (x 3 = 6.6) is greater than the pseudo-boundary value 6.5 of level IV. Equation (29) shows that RMD is equal to 1 for level IV and to 0 for other levels.
Therefore, the RMD matrix of these three indicators in Tianshengqiao II hydropower station can be expressed as follows: 0.400 0.700 0.100 0.000 0.000 0.388 0.724 0.112 0.000 0.000 0.000 1.000 Similarly, the RMD matrix of all rockburst cases can be obtained by repeating the above process.
(3) The Initialized SRMD The comprehensive evaluation model of Equation (10) can transform RMD into SRMD, but the weight of each index and model parameters α and p should be determined first. In this paper, α = 2 and p = 1 fuzzy optimization models are selected and eight different weight assignments are obtained from the references in Table 1. Based on the above analysis, combined with RMDs and eight weight assignment values, the initial SRMD can be calculated through Equation (10). Then, the initial SRMD is introduced into Equation (17) to obtain the level eigenvalue H, as shown in Figure 10a.
The comprehensive evaluation model of Equation (10) can transform RMD into SRMD, but the weight of each index and model parameters and should be determined first. In this paper, = 2 and = 1 fuzzy optimization models are selected and eight different weight assignments are obtained from the references in Table 1. Based on the above analysis, combined with RMDs and eight weight assignment values, the initial SRMD can be calculated through Equation (10). Then, the initial SRMD is introduced into Equation (17) to obtain the level eigenvalue H, as shown in Figure 10a.

(4) The Optimized SRMD
In order to optimize the SRMD, a three-layer BP neural network is established. The nodes of the input layer, the hidden layer, and the output layer are 4, 5, and 4, respectively. The input of the network is the initial SRMD value, and the expected output O corresponds to the actual rockburst strength level of each rockburst case. In order to facilitate classification and calculation, the four strength levels are converted into the following forms: The initial SRMD values of the first 14 cases are selected as training samples of the three-layer BP neural network. Combining the expected output O and the weight adjustment models of Equation (36) and (37), the connection weight between the input layer and the hidden layer and the connection weight ℎ between the hidden layer and the output layer can be calculated (where the learning rate γ is 0.05 and the network error E is 0.03). Then, based on the trained connection weight, the initial SRMD of the last four rockburst cases are input into the network as test samples for training and the network output is the optimized SRMD values of each rockburst case, of which the values are brought into Equation (17) to calculate the level eigenvalue H. The results are shown in Figure  10b.  (4) The Optimized SRMD In order to optimize the SRMD, a three-layer BP neural network is established. The nodes of the input layer, the hidden layer, and the output layer are 4, 5, and 4, respectively. The input of the network is the initial SRMD value, and the expected output O corresponds to the actual rockburst strength level of each rockburst case. In order to facilitate classification and calculation, the four strength levels are converted into the following forms: The initial SRMD values of the first 14 cases are selected as training samples of the three-layer BP neural network. Combining the expected output O and the weight adjustment models of Equation (36) and (37), the connection weight w jk between the input layer and the hidden layer and the connection weight w kh between the hidden layer and the output layer can be calculated (where the learning rate γ is 0.05 and the network error E is 0.03). Then, based on the trained connection weight, the initial SRMD of the last four rockburst cases are input into the network as test samples for training and the network output is the optimized SRMD values of each rockburst case, of which the values are brought into Equation (17)

Discussion
From Figure 10b, it can be seen that the predicted results of 18 rockburst cases based on the improved VFS method are basically consistent with the actual rockburst strength grade, except for No. 3, 5, and 9. No. 3, 5, and 9 belong to mixed-grade rockburst, and the actual rockburst grades are III~IV, II~III, and III~IV, respectively.
The new method improves the traditional VFS method from two aspects: simplifying the RMD calculation process and optimizing the SRMD. Then, what are the advantages of the new method compared with the traditional VFS method after such improvement? The following are discussed in detail: (1) The improved VFS method has higher computational efficiency: Equation (41) can clearly show that the simplified RMD calculation method can reduce the number of RMD function operations at least two times compared with the traditional RMD calculation method. With the number increase of classifications level, indicators, and samples, this advantage will become more prominent.
Since the distribution characteristics of RMD functions at different levels are known in advance (Figures 3c, 4c and 5c), the simplified RMD calculation method is simpler and more direct for RMD calculation, which greatly improves the operation efficiency of the improved VFS method. (2) The improved VFS method can verify the correctness of RMD calculation results at all times: RMD is the core theme of the improved VFS method, which concerns whether the final prediction results are correct or not. Therefore, the effective guarantee of the correctness of RMD calculation results is the premise of obtaining high-precision prediction results. According to RMD relationship characteristic Equations (20), (26), and (28), the range of traditional RMD value [0, 1] is reduced. For example, the eigenvalue of σ θ /σ c in Tianshengqiao II hydropower station is located in the interval of level II, so 0.5 ≤ u 12 (x) ≤ 1, 0 ≤ u 11 (x) ≤ 0.5, and 0 ≤ u 13 (x) ≤ 0.5. In addition, u 11 (x) + u 13 (x) = 0.5. These features can be used to verify the correctness of RMD calculation results. (3) The improved VFS method has higher prediction accuracy: From Figure 10a, it can be found that the prediction results based on the traditional VFS method are not ideal regardless of the weight. For example, the predicted results of No. 2, 4, and 11 are totally inconsistent with the actual strength grade. However, the prediction accuracy is significantly improved after the optimization of the initial SRMD using the BP neural network. This phenomenon shows that the combination of BP neural network and the VFS method can improve the classification and prediction ability of the model. (4) The improved VFS method has higher fault tolerance and practicability: From Figure 10a, we can also find that the prediction results of traditional VFS method are easily affected by the weight values and that different weight values may lead to different prediction results, resulting in misjudgment. For example, the prediction results of No. 2-4, 8-12, 14, and 18 rockburst cases span different levels in different weights. However, after the initial SRMDs are optimized by using the fuzzy BP optimization neural network model, the prediction results are basically consistent under different weights, which is within the same strength level, except for No. 3. The above results show that the improved VFS method has higher fault tolerance and anti-jamming ability and that its dependence on weight is low, so it has better practicability.

Conclusions
From the above analysis, we can draw the following conclusions.
(1) This research improved the traditional VFS method from two aspects: (i) simplifying the RMD calculation process and (ii) optimizing SRMD, which was applied to the prediction of rockburst strength. Good results have been achieved. (2) Compared to the traditional VFS method, the improved VFS has a clearer, more efficient calculation process and a more credible and stable prediction. The improved VFS method simplified the RMD calculation process by using the characteristic relationship of RMD at different levels and verified the correctness of RMD calculation results through these characteristics at all times. Besides, the improved VFS also uses the BP neural network to optimize the SRMD, which improves the prediction accuracy of SRMD. By this way, the influence of weight change on SRMD is also effectively avoided. Therefore, the improved method has higher fault tolerance rate and anti-jamming ability. (3) The original index data, the rockburst occurrence, and the intensity in underground projects all have certain dynamic variability and fuzziness, so it is difficult to express the rockburst criterion with an accurate relational expression. A more reasonable results can be obtained by using RMD and SRMD to predict. However, the application of improved VFS to rockburst prediction is still in the phase of theory and there are still some problems to be further explored, for example, how to reasonably construct RMD function, how to select and calculate the layers of BP neural network, the number of nodes and the connection weight matrix, etc. to make the rockburst prediction results more in line with the actual situation, especially for the accurate prediction of mix-level and intermediate-level rockburst.

Conflicts of Interest:
The authors declare no conflict interest.