Gaussian Process Model of Surrounding Rock Classification Based on Digital Characterization of Rock Mass Structure and Its Application

As the important data basis of surrounding rock classification, rock mass structural information obtained by traditional image processing and feature extraction algorithms could not be quantitatively analyzed because of the uncertainty and geometric randomness of structural planes. In this paper, based on straight line detection, intelligent scissors, and morphological edge detection algorithms, the multiple interpretation system of rock mass image including linear bunching extraction, magnetic tracking extraction, and multiparameter characterization was researched and developed, and the actual distribution information and the related probability distribution model of structural planes could be obtained directly. On the basis of this, plenty of corresponding random rating-values meeting the probability distribution models of these evaluation indices were gained by Monte Carlo Simulation. ,e distribution probability affiliated with different rock mass grade was attained by inductive statistics, and the robust evaluation of surrounding rock classification could be carried out. Taking the robust results as learning samples, the response model of surrounding rock grade based on Gaussian process classification was established, making the evaluation of surrounding rock subclassification more rapid and robust.


Introduction
With the development of rock mechanics theory, rock quality evaluation has played a key role in rock engineering [1][2][3][4]. Reasonable evaluation result of surrounding rock grade cannot only reflect the physical and mechanical characteristics of rock mass objectively, but also provide reliable data basis for construction method transformation and supporting parameters optimization.
Among various surrounding rock classification methods developed so far, the Rock Mass Rating (RMR) system and Norwegian Geotechnical Institute (NGI) Q system are the most commonly used [5,6]. Moreover, there are also RMI system method by Palmstrom, QTS method by Otakar, Basic Quality (BQ) method [7]; and Hydropower Classification (HC) method [8]. Among them, the Q-system classification method emphasizes on the numbers, roughness, and alterations of joints; therefore, if considering the directions of joints, this method is not suitable [9]. e RMR method can reflect the quality of hard rock mass, but it is not suitable for weak rock mass; furthermore, the evaluation results are correspondingly conservative. e HC method is only used in the classification of surrounding rocks in low ground stress areas instead of high ground stress areas, that is, the rock burst areas [1]. In short, during the complexity of engineering and geological conditions, different rock quality classification methods have different emphasis points.
Classification of surrounding rock is influenced by many factors, such as lithology, rock mass structure, and groundwater. It shows a highly nonlinear relationship between surrounding rock grades and the factors that affect them. In addition, during the impact of uncertainty and complexity of rock mass structural characteristics on surrounding rock classification, these methods often lead to inconsistent and unreliable results. Each indicator has its standard. As we all know, the hardness degree of rock (i.e., uniaxial compression strength of rock) as well as the volumetric joint count of rock mass both can be divided into several grades, but we cannot evaluate the surrounding rock grade based on one single indicator. Moreover, the grade numbers for each indicator are not the same.
ere is a problem with ambiguity and uncertainty here. Taking the RMR method as an example, the grade of each indicator can be calculated by the given criteria or standard, but the surrounding rock grades are classified based on the points approach in the end. us, many new methods and theories are applied to the classification of surrounding rock, for example, the statistical clustering model, gray theoretical model, fuzzy mathematic theory, and extenics theory [10][11][12][13][14][15][16][17][18][19]. Artificial neural networks, support vector machine, and Gradient Boosting Machine [20][21][22][23][24][25][26][27][28][29] have also been introduced in surrounding rock classification. One of the most important advantages of those methods is that they do not need a certain standard.
However, there are still many defects of these methods in practical application. For example, the membership function and weights of fuzzy mathematical model cannot be easily determined. e key evaluation factors of rough set method may be deleted during the process of attribute reduction. e probability distribution models of evaluation indices could not be considered in the interval theory. e diversity factor of set pair analysis method is difficult to determine in practice. ANN model has some limitations such as black box approach, low generalization capability, and overtraining problem [30]. SVM cannot escape from the blindness which is the common phenomenon in man-made choice of kernel function and its parameters.
What is more, as one of the most important indexes to evaluate surrounding rock grade, rock mass fragmentation degree is always considerably difficult to be described or defined precisely and quantitatively. Improvements in digital photogrammetry and image processing technology have prompted gradual maturation in noncontact methods designed for acquiring rock mass structural information [31,32]. Many scholars have obtained favorable results by combining these technologies with key block theory, the stereographic projection method, 3D geological model reconstruction, and fine numerical simulation [33][34][35][36][37][38][39][40][41][42][43][44][45], especially the fast shape-from-focus method proposed by Martišek [46] and the hierarchical-fractal annealing algorithm proposed by Zhou and Xiao [47][48][49][50][51]. During the uncertainty and geometric randomness of structural planes (joints, fissures, and fault zone), the information obtained by traditional image processing and feature extraction algorithms could not be still analyzed quantitatively and used directly.
In this paper, the Laohushan tunnel on the Jinan belt expressway in China is taken as an engineering background (see Figure 1). Based on straight line detection, intelligent scissors, and morphological edge detection algorithm, the multiple interpretation system of rock mass image including linear bunching extraction, magnetic tracking extraction, and multiparameter characterization was researched and developed, and the actual distribution information and the related probability distribution model of structural planes could be obtained directly. On this basis, plenty of corresponding random rating-values meeting the probability distribution models of these evaluation indices were gained by Monte Carlo Simulation. e distribution probability affiliated with different rock mass grade was attained by inductive statistics, and the robust evaluation of surrounding rock classification could be carried out. Taking the robust results as learning samples, the response model of surrounding rock grade based on Gaussian process classification was established, in order to evaluate surrounding rock subgrade more rapidly.

Multiple Interpretation of Structural Information for Rock Mass
Image analysis is a method to extract the data and information from images by a semiautomatic or automatic way, mainly including image processing, image segmentation, edge detection, and pattern recognition. Among them, the image segmentation separates the target part from the original picture based on some certain criteria. Image segmentation is the process of partitioning an image into parts or regions. is division into parts is often based on the characteristics of the pixels in the image. For example, one way to find regions in an image is to look for abrupt discontinuities in pixel values, which typically indicate edges. ese edges can define regions. Other methods divide the image into regions based on color values or texture. ere is no specific criteria here. e point is the purpose of the application of image segmentation. In this work, the certain criteria is to partition the jointed rock mass image into joints and complete rock mass. Since the complexity and randomness of image contour, there is no one image segmentation algorithm that can separate all the outlines and boundaries satisfactorily. Especially for those images of jointed rock mass, some comprehensive and semiautomatic extraction approaches are needed. Moreover, the information obtained by traditional image interpretation method could not be analyzed in quantity.
Specially aimed at the large-scale structural planes, specific geological body (fault zone, weak intercalated layer, etc.), and small random cracks, the multiple interpretation system that embraces linear bunching extraction, magnetic tracking extraction, and multiparameter characterization was researched and developed based on straight line detection, intelligent scissors, and morphologic edge detection algorithm, respectively. e next sections describe each of these. is is what we call line segments.

Mathematical Problems in Engineering
Contours are zones of the image where the gray level is changing fast enough from dark to light or the opposite. us, the gradient and level-lines of the image are key concepts. e algorithm starts from the computation of the levelline angle at each pixel to produce a level-line field, that is, a unit vector field such that all vectors are tangent to the levelline going through their base point. en, this field is segmented into connected regions of pixels that share the same level-line angle up to a certain tolerance τ. ese connected regions are called line support regions as shown in Figure 2.
Each line support region (a set of pixels) is a candidate for a line segment. e corresponding geometrical object (a rectangle in this case) must be associated with it. e principal inertial axis of the line support region is used as the main rectangle direction; the size of the rectangle is chosen to cover the full region. As the length limit, only the basic principle of LSD is introduced. Some authoritative literature [52] can be referenced for more algorithm and implementation procedures.

Application Example.
Although the edges can be well detected by choosing different operators or threshold values, the structure information such as gap, length, and dip cannot be counted directly. For the multiple interpretation system of rock mass image in this paper, the detection operator and the threshold value can be adjusted manually according to different characteristics of jointed rock image, in order to obtain the most reasonable result. en, the Hough transform is carried out. e line segments after being detected and transformed still are not discontiguous and have different dip angles. e structure information cannot be counted too. To simplify and automate the detection of joints, the polar coordinates of detected line segments (ρ i , θ i ) can be counted out firstly. en, those line segments with the same ρ i can be merged and connected, and the contiguous line group with different dip angles can be obtained.
Linear bunching extraction is based on the above LSD algorithm, which is applicable to extract those regular large-scale cracks. Without considering the effect of the crack gap and surface roughness, those groups of cracks with distinct linear features can be quickly extracted (see Figure 3). e basic characteristic parameters of the structural plane such as distance, dip angle, and trace length can be obtained, providing database for statistical analysis of rock mass integrity index and joint number per unit volume.
It is worth noting that the geometrical parameters of the extracted cracks are calculated with the dimension of the unit pixel. To obtain the real values, the process of scale transformation is necessary, which can be realized by the ratio of the image size to the actual value of rock mass. [53][54][55], the basic idea is to formulate the image as a weighted graph where pixels represent nodes with weighted edges connecting each pixel with its 8 adjacent neighbors. As illustrated in Figure 4, the original image is regarded as a graph, in which each pixel is a vertex and each vertex has several edges linking to its neighbor pixels. e contour linking two anchor points, f, e (selected by the user) is computed by searching minimal distance path [53]:

Intelligent Scissors. For the original intelligent scissors
., c N � e, and the local distance measure is defined as features, that is, Laplacian zero-crossing, gradient magnitude, gradient direction, edge pixel value, inside pixel value, and outside pixel value, respectively [53]. Each w(.) is a weighting coefficient for corresponding feature function. Computation of optimal paths corresponding to minimum cost contour segments is accomplished using Dijkstra's optimal path search.

Application Example.
Based on intelligent scissors [53], the magnetic tracing extraction algorithm which is suitable for a specific geological body, such as fault zone and weak intercalated layer, is proposed. e following sections present detailed steps about how to do so.
(1) Open one image of rock mass, click with the mouse the boundary endpoints of the specific geological body or large-wide cracks, and take it as one seed point s. (2) Move the seed point s along the boundary, and take the current mouse point t as the target point. (3) Search the optimal path between seed point s and target point t by using path searching algorithm [56]. (4) Freeze the searched optimal path between s and t, while the length of the optimal path happened to be larger than a given threshold. Take the target point t as a new seed point s, and search the next optimal path according to step 1, until the path searched has been closed. (5) Save all the pixels along the frozen boundary into a graphics layer as primitives. ese primitives kept are the extracted boundaries of the special geological body or large-wide cracks (see Figure 5).
For fault zone and weak intercalated layer, the geometrical parameters such as width and length are required. Taking Figure 6 as an example, the coordinate values of seed points and new points (x 1 , x 2 , . . ., x n and y 1 , y 2 , . . . , y n ) of structural plane boundaries can be obtained by using the above magnetic tracing extraction algorithm. Taking one point y i of any side boundary as the starting point, the distance between the starting point and each point of another side boundary can be calculated respectively, and the minimum value of all the calculated distances is the corresponding gap in point y i of the structural plane. With all the minimum values calculated, the mean gap of the structural plane can be obtained, that is, the mean values of all the minimum distances.
Based on the above research, the mean length of side boundaries L can be taken as the extension length of fault zone or weak intercalated layer, which can be expressed as follows: L � (length (e 2 ) + length (e 1 ))/2. Among them, length (e 1 ) and length (e 2 ) are the numbers of pixels of the two side boundaries. It is worth noting that the geometrical parameters of the specific geological body are calculated with the dimension of the unit pixel. To obtain the real values, the process of scale transformation is necessary.

Mathematical Morphology.
Mathematical morphology has been an important nonlinear tool for image analysis [57,58]. Suppose E represent the Euclidean space R. Let T be the gray levels of the images in R. en, an image f could be recognized as the mapping f: e two fundamental operators, which are the basic mappings of mathematical morphology, are defined as follows [57,58]: B is the structuring element in mathematical morphology. (u, v) and (x, y) represent the pixel coordinates. ⊕ and Θ denote the morphological dilation and erosion operators, respectively.
Another two important morphological operators based on the dilation and erosion are defined as follows [57,58]: where ○ and • demote the morphological opening and closing operators, respectively. Based on these basic morphological operators, the alternating filters are defined as follows [57,58]: where AF 1 , AF 2 , AF 3, and AF 4 denote the four types of morphological alternating filters, respectively. And, other morphological operators, such as the center and toggle operators [57,58], could be also defined based on these basic operators.

Edge Detection.
Edge detection aims to locate edge pixels in an image, which is used for subsequent image processing [59]. e principle of morphology-based edge detection is to choose structure elements with specific shapes to detect edges in an image. e edge detection methods use multidirectional structure elements to detect edges and assemble the edge images with equal coefficients; that is, the contributions of structure elements to the final edge image are equal [60,61].

Application Example.
is method is applicable to extract these tiny and discontinued structural planes or secondary cracks. As shown in Figure 7, the cracks extracted by using the morphologic edge detection algorithm are in good agreement with the actual cracks of the rock mass. It is noteworthy that the extraction information cannot be described or defined quantitatively. In other words, the occurrence information such as dip angle, trace length, and gap cannot be obtained directly. For that reason, the crack rate of the rock mass, which is related to the density and size of structural planes, is introduced. It can reasonably express the integrity of the rock mass and development degree of structural planes, especially for fractured rock mass.
Mathematical Problems in Engineering 5

Project Overview.
e Laohushan tunnel of the Jinan belt expressway has the largest one-way tunnel span among all the highway tunnels under construction in China (as shown in Figure 1). e dimensions of the main tunnel hole (width × height) are 17.608 m × 8.961 m, and the earthwork volume of excavation per meter is greater than 200 m 3 . Due to the influence of the Yanshan movement, some large-scale NNW cracks have developed in the monoclinic structure. From east to west, these cracks are the Wenzu crack, Dongwu crack, and Qianfoshan crack. e geological condition is presented in Figure 8. It is noteworthy that some rock masses exposed after excavation are complete and present thick-layer shapes, but continuous fault zone and weak intercalated layers, which have poor interlayer bonding and are filled with yellow altered contents, have also developed in these rock masses.

Statistics of Rock Mass Structural Information.
Based on the developed system of image acquisition and multiple interpretation of jointed rock mass structural information, the geometric distribution information of structural plane such as joints, cracks and fault zone etc., can be rapidly obtained and quantitatively analyzed. Taking Laohushan tunnel as an example, the images of rock mass are collected for each cycle. e trace length and spacing of structural planes are both obtained by Linear bunching extraction, and the crack rate of random joints is calculated based on morphologic edge detection. By using the multiple interpretation system, the three indices of every image can be obtained. With a large number of surrounding rock images, the probability distribution models of trance length, spacing, and crack rate can be counted out by using the statistical deduction method (as shown in Figure 9).

Membership Probability Calculation.
In this paper, the Basic Quality method and Q-system method are selected to evaluate the surrounding rock grades. Among them, BQ method provides one necessary and fundamental basis for the exploration, design, and quota compilation of rock engineering construction. It combines both the qualitative and quantitative methods to determine the basic quality of rock mass and then takes the characteristics of the specific engineering into account to obtain the rock mass classification. e reliability of BQ method has been verified in water conservancy and hydropower, transportation, railway, and mining projects since it was implemented in China [62].
To be simple, detailed evaluation procedures and criterions could be referenced [5,7,[63][64][65][66][67][68]. Besides, all the correction factors of BQ method and the stress reduction factor of Q-system method could be calculated and analyzed by interval-number mathematics. As the length limit, some related literature [69,70] can be referenced for implementation procedures.
Based on the above rock structural plane probability distribution model, as well as the measured rock uniaxial compressive strength R c and rock mass index RQD shown in Figure 10, lots of random numbers (n � 10000) meeting the probability distribution models and relevant parameters of these evaluation indices can be produced by the Monte Carlo Simulation. According to the evaluation process of each rock y 7 x 1 x 2 x 3 x 4 x 5 x 6 x 7 y n-1 x n-1 y n Figure 6: Calculation sketch of the crack gap.
mass classification system, plenty of corresponding random rating-values of these evaluation indices are also gained. e surrounding rock evaluation results and these probability distribution characteristics by the two rock mass classification systems can be attained easily by the inductive statistics (see Table 1) [71]. With full consideration of the probability distribution model of all the evaluation indices, the membership probability attached to level III is the largest for surrounding rock in ZK3 + 420∼380. Now, therefore, the reliable probabilities of different surrounding rock grades for the left line of Laohushan tunnel can also be calculated as shown in Table 2. It indicates that the evaluated results obtained by membership probabilities are generally in good agreement with the field-observed results. Compared with the original design grade of surrounding rock, the actual rock mass quality is much better. What is more, the reliable probability values attached to different surrounding rock grades can reflect the classification information more intuitively. Based on it, the geological investigators can make a more reasonable assessment of the changing trend of rock mass quality, providing a precise data basis for the transformation of construction method and optimization of support parameters. Since the membership probability of level V 1 was larger than that of level V 2 in ZK3 + 790∼770, the original designed double-side-drift method (DSD) was changed to center cross diaphragm method (CRD). Besides, during the gradeskipping phenomenon in Zk3 + 720, the CRD method was selected as the optimal excavation way by the construction unit. According to the reliability method of surrounding rock classification, the membership probability attached to different surrounding rock grades were calculated. It indicated that the membership probability of level IV 3 was the largest, and the membership probability of level IV 2 is also larger than that of V 1 (see Figure 11). erefore, the center diaphragm method (CD) and level IV supporting pattern was ultimately suggested in the practical construction. e follow-up monitoring data indicated that the accumulated value of vault settlement and convergence were within the limits permitted by the project. e reliability analysis theory was introduced to build the functions of different surrounding rock grades, and the reliable probabilities attached to different evaluation grades were calculated by the Monte Carlo method. Since the uncertainty and discreteness of the evaluation indices were sufficiently taken into account, the calculated reliability values by the reliability analysis method of surrounding    Mathematical Problems in Engineering rock subclassification could make a more steady and robust evaluation. However, this method was based on a great deal of statistical analysis, resulting in a waste of time, manpower, and financial resources. For those construction units and design units, all that is badly needed is a convenient and robust method for rock mass quality evaluation, which can meet the requirements of rapid construction.

Gaussian Process Classification Model of Surrounding Rock
Unlike conventional models, the Gaussian process is a novel machine learning model based on rigorous statistical learning theories and characterized by the self-adaptive determination of optimized hyperparameters [72,73]. It is thus suitable for complicated regression cases involving high   dimension, small sample population, and nonlinearity. Besides, the reliable probabilities calculated by the above method provide robust learning samples and evaluation standards for establishing the Gaussian process classification (GPC) model of surrounding rock grade.

Gaussian Process Classification eory.
A Gaussian process (GP) is fully specified by its mean function m (x) and kernel function k (x, x′), expressed as For notational simplicity, it is common to consider the mean function of the GP to be zero [73]. A GP model is a nonparametric, probabilistic (Bayesian) model in function space. One can think of a GP as defining a distribution over functions and inference taking place directly in the space of functions. e kernel function characterizes correlations between different data points in the Gaussian process and can be learnt from the data. e inference can be carried out directly under the GP framework by learning a kernel function from the training data. e kernel function studied in this paper is the widely used Radial Basis Function (RBF), which is defined as [73] where x and x′ are input vector pairs, l is the characteristic length scale, and σ 2 is the signal variance. e free parameters, that is, l and σ 2 , are called hyperparameters of GP model. Assume that we have a dataset D with n observations: D � {(x i , y i ), i � 1, 2, . . ., n}, where x is the input vector of dimension d and y is the class label +1/−1. e input d × n data matrix is denoted as X. Predictions for new inputs x′ are made out of this given training data using the GP model. GP binary classification is performed by first calculating the distribution over the latent function f corresponding to the test case [73]: where p(f|X, y) � p(y|f)p(f|X)/p(y|X) is the latent variable posterior, p(f ′ |X, x ′ , f) is the predictive posterior with respect to possible latent functions, and values of this could lie anywhere within the range of (−∞, +∞). So, the probabilistic prediction is made by where s can be any sigmoid function that "squashes" the prediction output to guarantee a valid probabilistic value within the range of [0, 1] [73]. For the multiclass classification problem, we can treat each one class as being independent of the others and apply binary classification individually to each (one) class versus the rest classes as shown in Figure 12. As the length limit, only the basic principle of GP is introduced. Some authoritative literature [73][74][75] can be referenced for more algorithm and implementation procedures.

Selection of Evaluation Indices.
e evaluation indices of common surrounding rock classification methods can be categorized into three parts: surrounding rock structure development feature (including strength of intact rock, shape distribution feature, and self-developmental condition of structural planes), geology environment feature (including in situ stress and groundwater developmental status), and engineering factor [66,68,71,76]. e selection of evaluation index should follow the mainstream specification or standard of surrounding rock classification and also be combined with the geology environment feature of tunnelling engineering itself. Based on the full consideration of the technical feasibility and the methodological operability, the uniaxial compression strength of rock σ ci , volumetric joint count of rock mass J V , crack rate of rock mass C r , spacing of structural planes d, and groundwater condition A w are selected as the evaluation indices of GPC model for surrounding rock grade for Laohushan tunnel [71,76,77].
Among them, the uniaxial compression strength of rock represents the strength of intact rock. e volumetric joint count of rock mass J V , which can be easily tested, is selected as the evaluation index of geometric characteristics of structural planes. Meanwhile, the crack rate of rock mass C r and spacing of structural planes d can also reflect the rock fragmentation degree. Considering the specific quantified requirements of the Gaussian process classification model, the evaluation criteria of groundwater condition by RMR method can be referenced. at is, groundwater condition A w can be determined by groundwater quantity, groundwater pressure, and seepage condition. For the sake of simplicity, some literature [6,64,77] can be referenced for more details. Since most of the burial depth of Laohushan tunnel is 20-60 m, the in situ stress would not be considered here. However, it does not mean the in situ stress can be ignored under any circumstances. Actually, the in situ stress should be taken into full consideration, especially for deep depth tunnel and underground engineering. Each specific tunnel engineering should be considered on a case-by-case basis.

Result Analysis.
With the developed system of image process and multiple interpretation of jointed rock mass structural information, the evaluation indices such as volumetric joint count of rock mass J V , crack rate of rock mass C r, and spacing of structural planes d can be easily measured. Besides, the rock mass quality is divided into G 1 , G 2 , G 3 , G 4, and G 5 five levels based on the above evaluation result of exposed rock mass and original designed excavation methods in Laohushan tunnel (the number of surrounding rock levels can be set to more than five levels only with one more binary classification process). erefore, the evaluation of surrounding rock grade in Laohushan tunnel can be treated as a problem of quinary classification, which can be divided into five binary classifications. at is, we can treat each one class as being independent of the others and apply binary classification individually to each (one) class versus the rest classes.
Taking X 1 , X 2 , X 3 , X 4 , and X 5 as input vectors and G 1 , G 2 , G 3 , G 4 , and G 5 as output vectors, the first 30 evaluation results of surrounding rock grade can be taken as learning samples, and the remaining 10 groups are predicting samples (see Table 3). During the search process of optimal superparameters, the convergence criteria of the conjugate gradient method can be controlled by maximum iteration. Taking (0, 0) as the original value, the optimal superparameters can be obtained by five times binary classification. e calculated result can be seen in Table 3.
It indicates that the evaluation results obtained by the GPC model are in good agreement with the actual surrounding rock grade except for ZK3 + 720∼700. e reason is that ZK3 + 720∼700 in Laohushan tunnel happens to be the changing section of rock mass quality. Actually, the surrounding rock grade evaluated is IV 3 for ZK3 + 720∼710 and IV 2 for ZK3 + 710∼700. On the whole, the GPC model of surrounding rock grade has been trained steadily. e learning result also indicates that the measured data during the tunnel construction should be selected as much as possible, in order to build a more robust training set. e amount of data is a critical factor that affects the evaluation results for any machine learning algorithms. With limited data, it is hard to obtain an accurate and robust model, especially for unexpected geological disasters. Hence, the trained model shall be used with caution. Besides, with the change of the lithology and geologic structure along the tunnel rapidly, the training samples will be influenced, and the predicted result accuracy will decrease. Generally, we will update the training sample information in time and delete the noise samples, that is, the place where geological conditions change suddenly.
Based on the established GPC model, the rock mass quality evaluation results of the right line in Laohushan tunnel can be seen in Table 4. By comparing the evaluation grades and in situ results, the GPC model has encouraging effects. It is demonstrated that the Gaussian process classification model based on multiparameters representation of rock mass structural information is reliable and efficient for rock mass quality evaluation, providing precise guidance for rock mass quality rapid evaluation during the whole process of tunnel construction.

Conclusion
(1) e multiple interpretation algorithms of rock mass structural information including linear bunching extraction, magnetic tracking extraction, and multiparameter characterization were researched, and, based on this, the geometric parameters of structural planes could be obtained and analyzed rapidly and accurately.
(2) e reliability analysis theory and Monte Carlo method were introduced to build the function of different surrounding rock grades. Since this method takes into account the uncertainty and discreteness of the evaluation indices sufficiently, the calculated reliability indices can make a more steady and robust evaluation of surrounding rock grade. (3) For meeting the requirements of rapid construction, the response model of surrounding rock grade based on Gaussian process classification was established, making the evaluation of surrounding rock subclassification more rapid and robust.  (4) e measured data during the tunnel construction should be selected as much as possible. e amount of data is a critical factor that affects the evaluation results for any machine learning algorithms. With limited data, it is hard to obtain an accurate and robust model, especially for unexpected geological disasters. Hence, the trained model shall be used with caution. (5) Although the feasibility and the methodological operability of selection of evaluation indices are considered fully, the aspect of theoretical basis is not quite there yet in this work. Besides, some index about self-development situation of joints such as the roughness is not considered in the proposed method. Each specific tunnel engineering should be considered on a case-by-case basis.

Data Availability
All data included in this study are available upon request by contact with the corresponding author.

Conflicts of Interest
e authors declare no conflicts of interest.

Authors' Contributions
Peng He wrote the original draft. Shangqu Sun was responsible for methodology. Gang Wang performed investigation. Wei-teng Li reviewed and edited.