A Novel Model Using Virtual State Variables and Bayesian Discriminant Analysis to Classify Surrounding Rock Stability

To accurately classify the stability of surrounding rock masses, a novel method (VSV-BDA) based on virtual state variables (VSVs) and Bayesian discriminant analysis (BDA) is proposed. The factors inﬂuencing stability are mapped by an artiﬁcial neural network (ANN) capable of recognizing the model of rock mass classiﬁcation, and the obtained output vector is treated as VSVs, which are veriﬁed as obeying a multinormal distribution with equal covariance matrixes by normal distribution testing and constructed statistics. The prediction variance ratio test method is introduced to determine the optimal dimension of the VSVs. The VSV-BDA model is constructed through the use of VSVs and the optimal dimension on the basis of the training samples, which are divided from the collected samples into three situations having diﬀerent numbers. ANN and BDA models are also constructed based on the same training samples. The predictions by the three models for the testing samples are compared; the results show that the proposed VSV-BDA model has high prediction accuracy and can be applied in practical engineering.


Introduction
Rock mass classification is generally considered a usable and practical approach to evaluating the stability of a rock mass in underground engineering [1,2]. Reasonable classification can reflect the mechanical characteristics of the rock mass and provide a reliable basis for the design of underground engineering excavations and supporting systems [3][4][5]. e classification approach has gradually evolved from the beginning with a single index to the current generation with multiple quantitative and/or qualitative indexes.
Rock quality design (RQD) was an efficient method proposed for rock mass quality assessment [6,7]. e rock structure rating (RSR) system was proposed for tunnel support design [8]. RSR was further developed into rock mass rating (RMR), a portion-rating system [9,10]. A rock tunneling quality index called the Q-system has been widely used in underground engineering and is closely related to RMR [11]. In recent times, with the development of numerical analysis and computer science, researchers have extensively studied rock mass classification and stability evaluation approaches: fuzzy analytic hierarchy process [12], artificial neural networks (ANNs) [13,14], distance discriminant analysis [15,16], set pair analysis [17], and so on. e stability of surrounding rocks is determined by many factors including geological conditions, exploitation factors, rock properties, and hydrogeological conditions [18][19][20][21][22][23][24]. ese factors have mutual effects and complicated nonlinear relationships. ANNs can represent nonlinear relationships, but a sufficient number of available training samples is needed [25,26]. Bayesian discriminant analysis (BDA) is a statistical method to classify samples that can use prior probability to increase the accuracy, but the relationships are too simple [27,28].
A novel method combining virtual state variables (VSVs) and Bayesian discriminant analysis (BDA), denoted as the VSV-BDA model, is proposed. e characteristics include the following: (1) the VSVs and prior probability is based to increase the accuracy of prediction; (2) quantile-quantile (QQ) plots and constructed statistics are used to test the normal distribution and equality of covariance instead of assuming a normal distribution and an equal covariance. e remainder of the paper includes three parts. In Section 2, the VSV-BDA approach is proposed with virtual state variables determination, Bayesian discriminant theory, and evaluation of the VSV-BDA. In Section 3, the proposed approach is used in rock mass classification and testing the accuracy. In Section 4, existing problems are analyzed and further study is analyzed.

VSV-BDA Approach
e factors influencing rock mass stability have nonlinear relationships [11], and ANNs can recognize these relationships [29][30][31]. e output of training samples predicted by an ANN model, denoted as Y − , is compared with the actual output by the residual variance ratio (RVR) method to determine the construction [32,33]. en, the first k − 1 layers (except for the output layer) and the regulating unit are used to construct a recognition network [34]. e output vector of the recognition network through nonlinear transformation of the input vector is called the VSVs by ANN, with the i-th variable denoted as v i . e vector with l components is called the virtual state vector, denoted as V, which does not have physical meaning but contains particular information characteristics to classify the stability, as shown in Figure 1. In addition, these variables instead of the influencing factors are used to construct the BDA model to classify the stability. e processes of VSV-BDA are summarized as follows: (1) collect data and group them into training samples and testing samples; (2) construct a multilayer ANN to determine the VSVs (influencing factors); (3) calculate the Mahalanobis distance of VSVs and use the QQ plot to normality test; the data should be transformed to satisfy the multinormal distribution by a Box-Cox transformation; (4) create discriminant functions based on the VSV of the training samples and cross-validation to estimate the accuracy; (5) use the constructed discriminant function to calculate the scores of the testing samples belonging to the collectivities; (6) classify the testing samples associated with the highest posterior probability and calculate the ratio of misclassification.
e process of BDA is illustrated in Figure 2.

Virtual State Variable Determination.
Interconnected processing elements, called neurons or cells, are used to construct the ANN, which offers a computational mechanism to acquire, compute, and represent a mapping from one information space layer to another, and to obtain a dataset to represent the relationships. An ANN identifies relationships by focusing on the parallel processing of many simple units, which can describe a complex function when combined. Essentially, the ANN is a gathering of simple processing units that exchange information that can be modified and filtered by the processing units' connections [35].
A multilayer ANN, which consists of an input layer, an output layer, and one or more hidden layers, is used to explain how to determine the VSVs. Factors influencing rock stability are taken as the input vector, represented by x � [x 1 , x 2 , . . . , x l ]. To simplify, a typical representation with one hidden layer whose function is denoted by f(u) is illustrated in Figure 2. e output is the resulting rock mass classification, denoted as Y − (Figure 3). e input dataset of the hidden layer is u j � l i�1 U ij x i , and the output dataset is f j (u j ) with a weight index w j . When the hidden layer function is a wavelet function or a linear function, the output Y − can be denoted as equation (1).
(1) e common transfer function in this paper is the sigmoid function [36], so equation (1) is replaced by equation (2).
where w � 1 and θ is a threshold value. ese parameters can be calculated by the ANN training process [29,37]. As aforementioned, the output of the hidden layer f j (u j ) is used as the VSV (v j ), and the virtual state vector is repre- e optimal dimension is determined by RVR, which is illustrated by equation (3). e initial state of V dimension is Λ, and the calculated RVR is denoted as S e (Λ). When the dimension is decreased or increased, called a state change and denoted as Λ, the RVR is represented by S e (Λ).
where Y − ANN (k) is the prediction by the ANN model, Y(k) is the actual stability level, and N is the number of training samples.
When S e is supposed to follow a normal distribution, to test the availability, a testing variable F is illustrated by For the given significance level α, the processes of RVR are shown as follows: (1) When F > F α (N − 1, N − 1), the RVR has a significant increase, and the prediction performance is worse, so the new state is invalid (2) When F < F 1−α (N − 1, N − 1), the RVR has a significant decrease, and the prediction performance is better, so the new state is valid 2 Shock and Vibration (3) In the third situation, the RVR is changing, but it cannot decide for the better or the worse, so the state remains unchanged

Bayesian Discriminant eory (BDA).
Bayesian discriminant analysis is a probability analysis method with various types of distribution density functions that should be obtained at the beginning. e prior distribution is used to describe the awareness level of training samples before extracting the testing samples; then, the posterior distribution is obtained by modifying the prior distribution from the testing samples.
Suppose there are k collectivities with p member indexes (considering p indexes):G 1 , G 2 , . . . , G k (k ≥ 2), and the covariance matrix Σ i > 0(i � 1, 2, . . . , k). e prior probability of G i is denoted as p i (i � 1, 2, . . . , k) and allocated by the proportion of G i to all collectivities with k i�1 p i � 1, as shown in where n i is the number of the training samples belonging to G i and n is the number of training samples. According to Bayesian theory, the posterior probability where f i (X) is the distribution function of G i . en, the optimal belonging can be obtained as 2.2.1. Normal Distribution Testing. e square difference of the Mahalanobis distance, denoted as d 2 i (V), is used to represent the distance between sample V and collectivities G i .
First k-1 layers ANN Residual variance ratio .

Bayesian discriminant analysis
Recognition networks . .   where μ i is the expect vector of G i . e square differences of the Mahalanobis distances are sorted from the smallest to the largest.

Regulating unit
A QQ plot, which is represented by the points (d 2 (i) , χ 2 (i) ), where χ 2 (i) is the chi-square distribution, is used to test whether the collectivities obey the multinormal distribution. When the points are all near the line passing through the origin with slope equal to one, the collectivities are regarded as obeying the multinormal distribution; otherwise, the data should be transformed to satisfy the multinormal distribution by a Box-Cox transformation. [33].

Discriminant Criterion.
rough the testing by the QQ plot and data transformation, all collectivities would obey the normal distribution; i.e., e discriminant function (W i (V)) of collectivity G i is linear, and the best divisions (R j ) are obtained by equations (10), (11), and (14). where e discriminant function (W i (X)) of collectivity G i is quadratic, and the best divisions are obtained by equations (12)- (14).

Estimation of Parameters.
In fact, μ i and Σ i are unknown, and an unbiased estimation can be obtained from the training samples. Supposing that training samples Σ i is defined as S i .
e statistic used to examine the equality of the covariance matrix is defined as ⎞ ⎠ , n i are not completely equal, 2p 2 + 3p − 1 (k + 1) 6(p + 1)(k − 1) , n i are completely equal.

Shock and Vibration
For the given α, the probability p � P(ξ > χ 2 α (f)) is calculated. If p > α, S i are completely equal; otherwise, they are not completely equal.

Evaluation of the VSV-BDA.
e resubstitution method and cross-validation method can also be used to estimate the reliability of the constructed discriminant criterion [33,[38][39][40]. For high accuracy, cross-validation is used, the principle of which is to choose one sample as a testing sample and use the rest of the collectivities as training samples to construct the discriminant criteria. e constructed criteria are used to classify the testing sample. e processes are shown as follows: (1) Choose one sample from G 1 as the testing sample and construct the discriminant criterion with the other n − 1 samples. (2) Classify the testing sample with the criterion constructed in process (1). (3) Repeat processes (1) and (2), and define the misclassification number as N 1 after all the samples of G 1 are tested. (4) Repeat the processes (1), (2), and (3) for collectivity G i (1 ≤ i ≤ k) and define the misclassification number as N i after all the samples of G i are tested. (5) e ratio of misclassification (η) can be calculated by

Influencing Factors of Rock Stability and Sample
Collection. Considering the typical factors, the convenience of factors can be obtained and compared in practical engineering. Five factors are selected as influencing factors, including the RQD, uniaxial axial compressive strength of the rock (R w ), rock mass integrity index (K v ), coefficient of structural surface strength (K f ), and groundwater discharge (ω) [41]. e rock stability is divided into five levels presented by collectivities G 1 ∼ G 5 : G 1 is level I of rock stability; G 2 is level II; G 3 is level III; G 4 is level IV; and G 5 is level V. e discriminant criterion belonging to G i is W i (i � 1, 2, . . . , 5). e surrounding rock classifications in the second-stage project of the Guangzhou pumped storage power station are collected to construct the VSV-BDA model and validation, as shown in Table 1 [41].

VSV-BDA Construction.
To test the generalization of the VSV-BDA model and ensure the reliability of the results, the situation with too many training samples and too few testing samples should be avoided. e collected samples in Table 1 are divided into three different situations to construct VSV-BDA models with ratios of training samples to testing samples equal to 20 : 17, 25 : 12, and 30 : 7. Taking the first ratio as an example, no. 1-20 are assigned as the training samples, and no. 21-37 are assigned as the testing samples. e method in Section 2 is applied. An ANN with 4 layers is constructed to determine the VSVs. e significance level is 0.05 (α � 0.05). e input layer and two hidden layers of ANN compose the first 3 layers of the recognition network. e five influencing factors listed in Table 1 are taken as the input vector. e transfer function implemented in the two hidden layers is the sigmoid function. When the number of training samples is smaller than the dimensions of the VSVs, the covariance estimation may be ill-posed [42]. Due to the limitation of collected samples, the number of VSVs can be decreased by increasing the number of hidden layers [32]. e optimal dimension of VSVs is determined by RVR, which has a structure of 5 × 6 × 5 × 1. en, the five VSVs obtained from the second hidden layer are used to construct the VSV-BDA model by the BDA theory.
An ANN [41] and BDA are also used to create the model with the training samples to obtain the classification functions or variables.

Validation and Comparison.
e testing samples are predicted by the three models to gain the output and compare the accuracy. e outputs of these 3 situations are shown in

Input layer
Hidden layer Output layer . . Figure 3: Pattern recognition network.   e average relative error ratios of testing samples are 5.88% for ANN, 3.92% for BDA, and 1.96% for VSV-BDA, as shown in Figure 4.

Shock and Vibration
In situation 2, the inaccurate prediction of samples occurs for sample 26 by ANN and by BDA. e average relative error ratios of testing samples are 3.03% for ANN and 3.03% for BDA, as shown in Figure 5.
In situation 3, the predictions by all three methods are correct. e three models have the same accuracy, which indicates that they can yield the correct result, as shown in Figure 6.
Compared with the ANN, the VSV-BDA model can use prior probability to increase the accuracy. Compared with BDA, the proposed method can represent complex relationships between the factors. From the aforementioned points, the VSV-BDA model has higher accuracy than the ANN and BDA models, especially when the training samples are insufficient. e more training samples are available, the more accurate the prediction. For the three methods, as the number of training samples increases, the accuracy increases.

Application of VSV-BDA Discriminant Criteria.
Due to the limitation of space, only situation 1 is explained in this paper. e samples in Table 1 Table 1. e results illustrate, through equation (23), that the covariances of the three collectivities are equal, so the discriminant functions of rock stability are linear functions.
e training samples are used to examine the accuracy of the discriminant criteria by the cross-validation method, as shown in Figure 4. e probability of misclassification is 0, which is perfect for classifying the rock mass stability. e testing samples are no. 21-37. e discriminant criteria, equations (25)- (27), are used to classify the testing samples. e output of sample no. 26 is level III but should be level IV, so the average relative error ratio is 1.96%.

Conclusions
(1) e accuracy of the VSV-BDA model in rock stability classification depends on the prior probability, probability density, and complex relationships among the VSVs. To satisfy the actual conditions, the QQ plot and statistic are introduced to test the probability distribution and the equality of the covariance instead of the assumed multinormal distribution and equality used in other studies, which could yield a more reasonable distribution and criterion to classify the stability. In the future, further study should be carried out to select the influencing factor and enhance the VSV-BDA model in practical engineering.

Data Availability
All data have been included in the manuscript.

Conflicts of Interest
e authors declare that they have no conflicts of interest.

Authors' Contributions
Jinglai Sun contributed to conceptualization, methodology, and writing-original draft. Zhaofei Chu contributed to validation, writing-review, and editing. Darui Ren and Yu song contributed to investigation. Baoguo Liu contributed to funding acquisition. Shaogang and Xinyang Guo contributed to data curation. 8214049).   Shock and Vibration