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

Authors: Peng He; Shang-qu Sun; Gang Wang; Wei-teng Li
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

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. The 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.


---

## Body

## 1. Introduction

With the development of rock mechanics theory, rock quality evaluation has played a key role in rock engineering [1–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]. The 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. The 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. There 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. Thus, 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–19]. Artificial neural networks, support vector machine, and Gradient Boosting Machine [20–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. The key evaluation factors of rough set method may be deleted during the process of attribute reduction. The probability distribution models of evaluation indices could not be considered in the interval theory. The 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–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–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 Figure1). 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. The 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.Figure 1
Location of Laohushan tunnel area and its satellite image [14].


## 2. 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. This 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. These edges can define regions. Other methods divide the image into regions based on color values or texture. There is no specific criteria here. The 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. The next sections describe each of these.


### 2.1. Linear Bunching Extraction Method

#### 2.1.1. Line Segment Detector Theory (LSD)

LSD is a linear-time Line Segment Detector giving subpixel accurate results. It is designed to work on any digital image without parameter tuning. LSD is aimed at detecting locally straight contours on images. This is what we call line segments. Contours are zones of the image where the gray level is changing fast enough from dark to light or the opposite. Thus, the gradient and level-lines of the image are key concepts.The algorithm starts from the computation of the level-line angle at each pixel to produce a level-line field, that is, a unit vector field such that all vectors are tangent to the level-line going through their base point. Then, this field is segmented into connected regions of pixels that share the same level-line angle up to a certain toleranceτ. These connected regions are called line support regions as shown in Figure 2.Figure 2
Diagrammatic sketch of line support regions. (a) Image. (b) Level-line field. (c) Line support regions.
(a)(b)(c)Each line support region (a set of pixels) is a candidate for a line segment. The corresponding geometrical object (a rectangle in this case) must be associated with it. The 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.


#### 2.1.2. 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. Then, the Hough transform is carried out. The line segments after being detected and transformed still are not discontiguous and have different dip angles. The 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. Then, 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 Figure3). The 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.Figure 3
Result of linear group extraction.
(a)(b)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.


### 2.2. Magnetic Tracking Extraction Method

#### 2.2.1. Intelligent Scissors

For the original intelligent scissors [53–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. The contour linking two anchor points, f, e (selected by the user) is computed by searching minimal distance path [53]:(1)C=argminC∑k=1Ncostck,ck+1,where contour C = {c1, c2, … ,cN}, c1 = f, …, cN = e, and the local distance measure is defined as(2)costp,q=wZfZq+wDfDp,q+wGfGq+wpfpq+wIfIq+wOfOq,where fZ (q), fD (p, q), fG (q), fP (q), fI (q), and fO (q) are image 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.Figure 4
Converting an image into a graph for boundary determination. (a) Image, (b) graph, and (c) boundary.
(a)(b)(c)


#### 2.2.2. 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. The 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 points.(2)
Move the seed points along the boundary, and take the current mouse point t as the target point.(3)
Search the optimal path between seed points and target point t by using path searching algorithm [56].(4)
Freeze the searched optimal path betweens 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. These primitives kept are the extracted boundaries of the special geological body or large-wide cracks (see Figure5).Figure 5
Result of magnetic tracing extraction.
(a)(b)For fault zone and weak intercalated layer, the geometrical parameters such as width and length are required. Taking Figure6 as an example, the coordinate values of seed points and new points (x1, x2, …, xn and y1, y2, … , yn) of structural plane boundaries can be obtained by using the above magnetic tracing extraction algorithm. Taking one point yi 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 yi 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.Figure 6
Calculation sketch of the crack gap.Based on the above research, the mean length of side boundariesL can be taken as the extension length of fault zone or weak intercalated layer, which can be expressed as follows: L = (length (e2) + length (e1))/2. Among them, length (e1) and length (e2) 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.


### 2.3. Morphologic Edge Detection Method

#### 2.3.1. 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. Then, an image f could be recognized as the mapping f: E ⟶ T, f ∈ F (E, T).The two fundamental operators, which are the basic mappings of mathematical morphology, are defined as follows [57, 58]:(3)f⊕B=maxu,vfx−u,y−v+Bu,v,fΘB=minu,vfx+u,y+v−Bu,v.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]:(4)f∘B=fΘB⊕B,f•B=f⊕BΘB,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]:(5)AF1f=f∘B•B,AF2f=f•B∘B,AF3f=f∘B•B∘B,AF4f=f•B∘B•B,where AF1, AF2, AF3, and AF4 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.


#### 2.3.2. Edge Detection

Edge detection aims to locate edge pixels in an image, which is used for subsequent image processing [59]. The principle of morphology-based edge detection is to choose structure elements with specific shapes to detect edges in an image. The 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].


#### 2.3.3. Application Example

This method is applicable to extract these tiny and discontinued structural planes or secondary cracks. As shown in Figure7, 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.Figure 7
Result of morphology edge detection.
(a)(b)


## 2.1. Linear Bunching Extraction Method

### 2.1.1. Line Segment Detector Theory (LSD)

LSD is a linear-time Line Segment Detector giving subpixel accurate results. It is designed to work on any digital image without parameter tuning. LSD is aimed at detecting locally straight contours on images. This is what we call line segments. Contours are zones of the image where the gray level is changing fast enough from dark to light or the opposite. Thus, the gradient and level-lines of the image are key concepts.The algorithm starts from the computation of the level-line angle at each pixel to produce a level-line field, that is, a unit vector field such that all vectors are tangent to the level-line going through their base point. Then, this field is segmented into connected regions of pixels that share the same level-line angle up to a certain toleranceτ. These connected regions are called line support regions as shown in Figure 2.Figure 2
Diagrammatic sketch of line support regions. (a) Image. (b) Level-line field. (c) Line support regions.
(a)(b)(c)Each line support region (a set of pixels) is a candidate for a line segment. The corresponding geometrical object (a rectangle in this case) must be associated with it. The 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.


### 2.1.2. 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. Then, the Hough transform is carried out. The line segments after being detected and transformed still are not discontiguous and have different dip angles. The 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. Then, 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 Figure3). The 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.Figure 3
Result of linear group extraction.
(a)(b)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.


## 2.1.1. Line Segment Detector Theory (LSD)

LSD is a linear-time Line Segment Detector giving subpixel accurate results. It is designed to work on any digital image without parameter tuning. LSD is aimed at detecting locally straight contours on images. This is what we call line segments. Contours are zones of the image where the gray level is changing fast enough from dark to light or the opposite. Thus, the gradient and level-lines of the image are key concepts.The algorithm starts from the computation of the level-line angle at each pixel to produce a level-line field, that is, a unit vector field such that all vectors are tangent to the level-line going through their base point. Then, this field is segmented into connected regions of pixels that share the same level-line angle up to a certain toleranceτ. These connected regions are called line support regions as shown in Figure 2.Figure 2
Diagrammatic sketch of line support regions. (a) Image. (b) Level-line field. (c) Line support regions.
(a)(b)(c)Each line support region (a set of pixels) is a candidate for a line segment. The corresponding geometrical object (a rectangle in this case) must be associated with it. The 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.


## 2.1.2. 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. Then, the Hough transform is carried out. The line segments after being detected and transformed still are not discontiguous and have different dip angles. The 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. Then, 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 Figure3). The 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.Figure 3
Result of linear group extraction.
(a)(b)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.


## 2.2. Magnetic Tracking Extraction Method

### 2.2.1. Intelligent Scissors

For the original intelligent scissors [53–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. The contour linking two anchor points, f, e (selected by the user) is computed by searching minimal distance path [53]:(1)C=argminC∑k=1Ncostck,ck+1,where contour C = {c1, c2, … ,cN}, c1 = f, …, cN = e, and the local distance measure is defined as(2)costp,q=wZfZq+wDfDp,q+wGfGq+wpfpq+wIfIq+wOfOq,where fZ (q), fD (p, q), fG (q), fP (q), fI (q), and fO (q) are image 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.Figure 4
Converting an image into a graph for boundary determination. (a) Image, (b) graph, and (c) boundary.
(a)(b)(c)


### 2.2.2. 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. The 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 points.(2)
Move the seed points along the boundary, and take the current mouse point t as the target point.(3)
Search the optimal path between seed points and target point t by using path searching algorithm [56].(4)
Freeze the searched optimal path betweens 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. These primitives kept are the extracted boundaries of the special geological body or large-wide cracks (see Figure5).Figure 5
Result of magnetic tracing extraction.
(a)(b)For fault zone and weak intercalated layer, the geometrical parameters such as width and length are required. Taking Figure6 as an example, the coordinate values of seed points and new points (x1, x2, …, xn and y1, y2, … , yn) of structural plane boundaries can be obtained by using the above magnetic tracing extraction algorithm. Taking one point yi 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 yi 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.Figure 6
Calculation sketch of the crack gap.Based on the above research, the mean length of side boundariesL can be taken as the extension length of fault zone or weak intercalated layer, which can be expressed as follows: L = (length (e2) + length (e1))/2. Among them, length (e1) and length (e2) 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.


## 2.2.1. Intelligent Scissors

For the original intelligent scissors [53–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. The contour linking two anchor points, f, e (selected by the user) is computed by searching minimal distance path [53]:(1)C=argminC∑k=1Ncostck,ck+1,where contour C = {c1, c2, … ,cN}, c1 = f, …, cN = e, and the local distance measure is defined as(2)costp,q=wZfZq+wDfDp,q+wGfGq+wpfpq+wIfIq+wOfOq,where fZ (q), fD (p, q), fG (q), fP (q), fI (q), and fO (q) are image 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.Figure 4
Converting an image into a graph for boundary determination. (a) Image, (b) graph, and (c) boundary.
(a)(b)(c)


## 2.2.2. 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. The 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 points.(2)
Move the seed points along the boundary, and take the current mouse point t as the target point.(3)
Search the optimal path between seed points and target point t by using path searching algorithm [56].(4)
Freeze the searched optimal path betweens 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. These primitives kept are the extracted boundaries of the special geological body or large-wide cracks (see Figure5).Figure 5
Result of magnetic tracing extraction.
(a)(b)For fault zone and weak intercalated layer, the geometrical parameters such as width and length are required. Taking Figure6 as an example, the coordinate values of seed points and new points (x1, x2, …, xn and y1, y2, … , yn) of structural plane boundaries can be obtained by using the above magnetic tracing extraction algorithm. Taking one point yi 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 yi 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.Figure 6
Calculation sketch of the crack gap.Based on the above research, the mean length of side boundariesL can be taken as the extension length of fault zone or weak intercalated layer, which can be expressed as follows: L = (length (e2) + length (e1))/2. Among them, length (e1) and length (e2) 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.


## 2.3. Morphologic Edge Detection Method

### 2.3.1. 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. Then, an image f could be recognized as the mapping f: E ⟶ T, f ∈ F (E, T).The two fundamental operators, which are the basic mappings of mathematical morphology, are defined as follows [57, 58]:(3)f⊕B=maxu,vfx−u,y−v+Bu,v,fΘB=minu,vfx+u,y+v−Bu,v.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]:(4)f∘B=fΘB⊕B,f•B=f⊕BΘB,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]:(5)AF1f=f∘B•B,AF2f=f•B∘B,AF3f=f∘B•B∘B,AF4f=f•B∘B•B,where AF1, AF2, AF3, and AF4 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.


### 2.3.2. Edge Detection

Edge detection aims to locate edge pixels in an image, which is used for subsequent image processing [59]. The principle of morphology-based edge detection is to choose structure elements with specific shapes to detect edges in an image. The 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].


### 2.3.3. Application Example

This method is applicable to extract these tiny and discontinued structural planes or secondary cracks. As shown in Figure7, 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.Figure 7
Result of morphology edge detection.
(a)(b)


## 2.3.1. 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. Then, an image f could be recognized as the mapping f: E ⟶ T, f ∈ F (E, T).The two fundamental operators, which are the basic mappings of mathematical morphology, are defined as follows [57, 58]:(3)f⊕B=maxu,vfx−u,y−v+Bu,v,fΘB=minu,vfx+u,y+v−Bu,v.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]:(4)f∘B=fΘB⊕B,f•B=f⊕BΘB,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]:(5)AF1f=f∘B•B,AF2f=f•B∘B,AF3f=f∘B•B∘B,AF4f=f•B∘B•B,where AF1, AF2, AF3, and AF4 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.


## 2.3.2. Edge Detection

Edge detection aims to locate edge pixels in an image, which is used for subsequent image processing [59]. The principle of morphology-based edge detection is to choose structure elements with specific shapes to detect edges in an image. The 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].


## 2.3.3. Application Example

This method is applicable to extract these tiny and discontinued structural planes or secondary cracks. As shown in Figure7, 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.Figure 7
Result of morphology edge detection.
(a)(b)


## 3. Reliability Analysis of Surrounding Rock Classification

### 3.1. Project Overview

The 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 Figure1). The 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 m3. 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. The 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.Figure 8
Geological section map of Laohushan tunnel [14].


### 3.2. 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. The 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 Figure9).Figure 9
Probabilistic distribution models and parameters of structural information. (a) Probability histogram of trace length. (b) Probability histogram of crack rate of rock mass. (c) Distribution contour of structural plane sets. (d) Probability histogram of spacing.
(a)(b)(c)(d)


### 3.3. 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. The 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–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 strengthRc 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 mass classification system, plenty of corresponding random rating-values of these evaluation indices are also gained. The 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.Figure 10
Probabilistic distribution models and their relevant parameters of evaluation indexes. (a) Rock uniaxial compressive strength. (b) Distribution contour of RQD. (c) Rock mass integrity indexKv.
(a)(b)(c)Table 1
Membership probability attached to different surrounding rock grades in ZK3 + 420∼380.
Surrounding rock gradeMembership probability obtained by Q-system methodMembership probability obtained by BQ methodMembership gradeI0.01640Max(Pi) = III (i = I, II, …, V)II0.06680.0131III0.66150.6608IV0.21010.3249V0.04520.0012Now, therefore, the reliable probabilities of different surrounding rock grades for the left line of Laohushan tunnel can also be calculated as shown in Table2. 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.Table 2
Comparison of classification results of Laohushan tunnel left line.
Tunnel chainageDesigned gradeActual gradeEvaluated gradeMax membership probabilityDesigned excavation methodThe changed excavation methodZK3 + 810∼790V2V2V20.6377Double-side-drift methodDouble-side-drift methodZK3 +V790∼770V2V1V10.6079Double-side-drift methodCRD methodZK3 + 770∼750V2V1V10.4442Double-side-drift methodCRD methodZK3 + 750∼720IV3IV3IV30.4073CRD methodCD methodZK3 + 720∼700IV3IV3IV2IV30.3564CRD methodCD methodZK3 + 700∼675IV3IV2IV20.4555CRD methodCD methodZK3 + 675∼660IV2IV2IV20.4094CD methodCD methodZK3 + 660∼640IV2IV2IV20.4043CD methodCD methodZK3 + 640∼620IV2IV3IV30.4494CD methodCD methodZK3 + 620∼600IV2IV2IV20.4165CD methodCD method…——————ZK3 + 420∼380IV2III2III20.5107CD methodBenching methodSince the membership probability of level V1 was larger than that of level V2 in ZK3 + 790∼770, the original designed double-side-drift method (DSD) was changed to center cross diaphragm method (CRD). Besides, during the grade-skipping 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 IV3 was the largest, and the membership probability of level IV2 is also larger than that of V1 (see Figure 11). Therefore, the center diaphragm method (CD) and level IV supporting pattern was ultimately suggested in the practical construction. The follow-up monitoring data indicated that the accumulated value of vault settlement and convergence were within the limits permitted by the project.Figure 11
Distribution of reliable probability attached to different surrounding rock levels in Laohushan tunnel.The 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 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.


## 3.1. Project Overview

The 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 Figure1). The 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 m3. 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. The 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.Figure 8
Geological section map of Laohushan tunnel [14].


## 3.2. 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. The 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 Figure9).Figure 9
Probabilistic distribution models and parameters of structural information. (a) Probability histogram of trace length. (b) Probability histogram of crack rate of rock mass. (c) Distribution contour of structural plane sets. (d) Probability histogram of spacing.
(a)(b)(c)(d)


## 3.3. 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. The 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–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 strengthRc 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 mass classification system, plenty of corresponding random rating-values of these evaluation indices are also gained. The 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.Figure 10
Probabilistic distribution models and their relevant parameters of evaluation indexes. (a) Rock uniaxial compressive strength. (b) Distribution contour of RQD. (c) Rock mass integrity indexKv.
(a)(b)(c)Table 1
Membership probability attached to different surrounding rock grades in ZK3 + 420∼380.
Surrounding rock gradeMembership probability obtained by Q-system methodMembership probability obtained by BQ methodMembership gradeI0.01640Max(Pi) = III (i = I, II, …, V)II0.06680.0131III0.66150.6608IV0.21010.3249V0.04520.0012Now, therefore, the reliable probabilities of different surrounding rock grades for the left line of Laohushan tunnel can also be calculated as shown in Table2. 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.Table 2
Comparison of classification results of Laohushan tunnel left line.
Tunnel chainageDesigned gradeActual gradeEvaluated gradeMax membership probabilityDesigned excavation methodThe changed excavation methodZK3 + 810∼790V2V2V20.6377Double-side-drift methodDouble-side-drift methodZK3 +V790∼770V2V1V10.6079Double-side-drift methodCRD methodZK3 + 770∼750V2V1V10.4442Double-side-drift methodCRD methodZK3 + 750∼720IV3IV3IV30.4073CRD methodCD methodZK3 + 720∼700IV3IV3IV2IV30.3564CRD methodCD methodZK3 + 700∼675IV3IV2IV20.4555CRD methodCD methodZK3 + 675∼660IV2IV2IV20.4094CD methodCD methodZK3 + 660∼640IV2IV2IV20.4043CD methodCD methodZK3 + 640∼620IV2IV3IV30.4494CD methodCD methodZK3 + 620∼600IV2IV2IV20.4165CD methodCD method…——————ZK3 + 420∼380IV2III2III20.5107CD methodBenching methodSince the membership probability of level V1 was larger than that of level V2 in ZK3 + 790∼770, the original designed double-side-drift method (DSD) was changed to center cross diaphragm method (CRD). Besides, during the grade-skipping 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 IV3 was the largest, and the membership probability of level IV2 is also larger than that of V1 (see Figure 11). Therefore, the center diaphragm method (CD) and level IV supporting pattern was ultimately suggested in the practical construction. The follow-up monitoring data indicated that the accumulated value of vault settlement and convergence were within the limits permitted by the project.Figure 11
Distribution of reliable probability attached to different surrounding rock levels in Laohushan tunnel.The 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 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.


## 4. 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.


### 4.1. Gaussian Process Classification Theory

A Gaussian process (GP) is fully specified by its mean functionm (x) and kernel function k (x, x′), expressed as(6)fx∼GPmx,kx,x′.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.The kernel function characterizes correlations between different data points in the Gaussian process and can be learnt from the data. The inference can be carried out directly under the GP framework by learning a kernel function from the training data. The kernel function studied in this paper is the widely used Radial Basis Function (RBF), which is defined as [73](7)kx,x′=σ2exp−x−x′22l2,where x and x′ are input vector pairs, l is the characteristic length scale, and σ2 is the signal variance. The free parameters, that is, l and σ2, are called hyperparameters of GP model.Assume that we have a datasetD with n observations: D = {(xi, yi), i = 1, 2, …, n}, where x is the input vector of dimension d and y is the class label +1/−1. The 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]:(8)pf′|X,y,x′=∫pf′|X,x′,fpf|X,ydf,where pf|X,y=py|fpf|X/py|X is the latent variable posterior, pf′|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(9)π¯′=py′=+1|X,y,x′=∫sf′pf′|X,y,x′df′,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 Figure12. As the length limit, only the basic principle of GP is introduced. Some authoritative literature [73–75] can be referenced for more algorithm and implementation procedures.Figure 12
Sketch map of classification by GP model.


### 4.2. Selection of Evaluation Indices

The 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]. The 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 JV, crack rate of rock mass Cr, spacing of structural planes d, and groundwater condition Aw 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. The volumetric joint count of rock massJV, which can be easily tested, is selected as the evaluation index of geometric characteristics of structural planes. Meanwhile, the crack rate of rock mass Cr 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. That is, groundwater condition Aw 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.


### 4.3. 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 massJV, crack rate of rock mass Cr, and spacing of structural planes d can be easily measured. Besides, the rock mass quality is divided into G1, G2, G3, G4, and G5 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). Therefore, 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. That 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.TakingX1, X2, X3, X4, and X5 as input vectors and G1, G2, G3, G4, and G5 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. The calculated result can be seen in Table 3.Table 3
The training samples of GPC model for surrounding rock grade.
NumberTunnel chainageCompression strengthσci (MPa)Volumetric joint countJVCrack rateCr (%)Spacingd (m)GroundwaterAwDesigned gradeActual gradeEvaluation grade1ZK3  +  810∼79051.073155.850.1214VVV2ZK3 + 790∼77046.822751.490.1513VVV3ZK3 + 770∼75050.602746.220.2013VVV4ZK3 + 750∼72053.772635.740.2212IV3IV3IV35ZK3 + 720∼70052.642431.370.2310IV3IV32IV36ZK3 + 700∼67550.862227.970.259IV3IV2IV27ZK3 + 675∼66054.152429.310.227IV2IV2IV28ZK3 + 660∼64052.772331.420.235IV2IV2IV29ZK3 + 640∼62046.312747.250.184IV2IV3IV310ZK3 + 620∼60051.431839.830.305IV2IV2IV211ZK3 + 600∼58049.541429.110.345IV2IV1IV112ZK3 + 580∼56056.171826.100.338IV2IV1IV113ZK3 + 560∼54057.221928.280.3011IV2IV1IV114ZK3 + 540∼52056.811622.690.3011IV2IV1IV115ZK3 + 520∼50057.151827.070.299IV2IV1IV116ZK3 + 500∼48057.661621.220.3112IV2IV1IV117ZK3 + 480∼46055.971415.450.3513IV2IV1IV118ZK3 + 460∼44058.061412.960.3712IV2IV1IV119ZK3 + 440∼42063.05139.680.4213IV2III2III220ZK3 + 420∼40064.281113.710.4514IV2III2III221ZK3 + 400∼38063.32912.030.4613IV2III2III2……————————28ZK3 + 260∼24061.91109.240.4213IV2III2III229ZK3 + 240∼22058.3098.750.4714IV2III2III230ZK3 + 220∼20059.66109.270.4513IV2III2III2It 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. The 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 IV3 for ZK3 + 720∼710 and IV2 for ZK3 + 710∼700. On the whole, the GPC model of surrounding rock grade has been trained steadily.The 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. The 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 Table4. 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.Table 4
Surrounding rock grade evaluation results in the right line of Laohushan tunnel.
NumberTunnel chainageCompression strengthσci (MPa)Volumetric joint countJVCrack rateCr (%)Spacingd (m)GroundwaterAwDesigned gradeActual gradeEvaluation grade1YK3 + 818∼79848.333260.710.1113VVV2YK3 + 798∼77846.083258.460.1115VVV3YK3 + 778∼75848.302651.050.1813VIV3IV34YK3 + 758∼73849.072133.630.2110IV2IV3IV35YK3 + 738∼71854.132234.250.199IV2IV2IV26YK3 + 718∼69844.921831.050.287IV2IV2IV27YK3 + 698∼67849.851733.320.264IV3IV2IV28YK3 + 678∼65852.191627.390.275IV3IV1IV19YK3 + 658∼63843.981832.250.248IV3IV2IV210YK3 + 638∼61841.142038.760.2311IV3IV3IV3


## 4.1. Gaussian Process Classification Theory

A Gaussian process (GP) is fully specified by its mean functionm (x) and kernel function k (x, x′), expressed as(6)fx∼GPmx,kx,x′.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.The kernel function characterizes correlations between different data points in the Gaussian process and can be learnt from the data. The inference can be carried out directly under the GP framework by learning a kernel function from the training data. The kernel function studied in this paper is the widely used Radial Basis Function (RBF), which is defined as [73](7)kx,x′=σ2exp−x−x′22l2,where x and x′ are input vector pairs, l is the characteristic length scale, and σ2 is the signal variance. The free parameters, that is, l and σ2, are called hyperparameters of GP model.Assume that we have a datasetD with n observations: D = {(xi, yi), i = 1, 2, …, n}, where x is the input vector of dimension d and y is the class label +1/−1. The 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]:(8)pf′|X,y,x′=∫pf′|X,x′,fpf|X,ydf,where pf|X,y=py|fpf|X/py|X is the latent variable posterior, pf′|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(9)π¯′=py′=+1|X,y,x′=∫sf′pf′|X,y,x′df′,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 Figure12. As the length limit, only the basic principle of GP is introduced. Some authoritative literature [73–75] can be referenced for more algorithm and implementation procedures.Figure 12
Sketch map of classification by GP model.


## 4.2. Selection of Evaluation Indices

The 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]. The 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 JV, crack rate of rock mass Cr, spacing of structural planes d, and groundwater condition Aw 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. The volumetric joint count of rock massJV, which can be easily tested, is selected as the evaluation index of geometric characteristics of structural planes. Meanwhile, the crack rate of rock mass Cr 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. That is, groundwater condition Aw 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.


## 4.3. 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 massJV, crack rate of rock mass Cr, and spacing of structural planes d can be easily measured. Besides, the rock mass quality is divided into G1, G2, G3, G4, and G5 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). Therefore, 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. That 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.TakingX1, X2, X3, X4, and X5 as input vectors and G1, G2, G3, G4, and G5 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. The calculated result can be seen in Table 3.Table 3
The training samples of GPC model for surrounding rock grade.
NumberTunnel chainageCompression strengthσci (MPa)Volumetric joint countJVCrack rateCr (%)Spacingd (m)GroundwaterAwDesigned gradeActual gradeEvaluation grade1ZK3  +  810∼79051.073155.850.1214VVV2ZK3 + 790∼77046.822751.490.1513VVV3ZK3 + 770∼75050.602746.220.2013VVV4ZK3 + 750∼72053.772635.740.2212IV3IV3IV35ZK3 + 720∼70052.642431.370.2310IV3IV32IV36ZK3 + 700∼67550.862227.970.259IV3IV2IV27ZK3 + 675∼66054.152429.310.227IV2IV2IV28ZK3 + 660∼64052.772331.420.235IV2IV2IV29ZK3 + 640∼62046.312747.250.184IV2IV3IV310ZK3 + 620∼60051.431839.830.305IV2IV2IV211ZK3 + 600∼58049.541429.110.345IV2IV1IV112ZK3 + 580∼56056.171826.100.338IV2IV1IV113ZK3 + 560∼54057.221928.280.3011IV2IV1IV114ZK3 + 540∼52056.811622.690.3011IV2IV1IV115ZK3 + 520∼50057.151827.070.299IV2IV1IV116ZK3 + 500∼48057.661621.220.3112IV2IV1IV117ZK3 + 480∼46055.971415.450.3513IV2IV1IV118ZK3 + 460∼44058.061412.960.3712IV2IV1IV119ZK3 + 440∼42063.05139.680.4213IV2III2III220ZK3 + 420∼40064.281113.710.4514IV2III2III221ZK3 + 400∼38063.32912.030.4613IV2III2III2……————————28ZK3 + 260∼24061.91109.240.4213IV2III2III229ZK3 + 240∼22058.3098.750.4714IV2III2III230ZK3 + 220∼20059.66109.270.4513IV2III2III2It 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. The 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 IV3 for ZK3 + 720∼710 and IV2 for ZK3 + 710∼700. On the whole, the GPC model of surrounding rock grade has been trained steadily.The 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. The 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 Table4. 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.Table 4
Surrounding rock grade evaluation results in the right line of Laohushan tunnel.
NumberTunnel chainageCompression strengthσci (MPa)Volumetric joint countJVCrack rateCr (%)Spacingd (m)GroundwaterAwDesigned gradeActual gradeEvaluation grade1YK3 + 818∼79848.333260.710.1113VVV2YK3 + 798∼77846.083258.460.1115VVV3YK3 + 778∼75848.302651.050.1813VIV3IV34YK3 + 758∼73849.072133.630.2110IV2IV3IV35YK3 + 738∼71854.132234.250.199IV2IV2IV26YK3 + 718∼69844.921831.050.287IV2IV2IV27YK3 + 698∼67849.851733.320.264IV3IV2IV28YK3 + 678∼65852.191627.390.275IV3IV1IV19YK3 + 658∼63843.981832.250.248IV3IV2IV210YK3 + 638∼61841.142038.760.2311IV3IV3IV3


## 5. Conclusion

(1)
The 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)
The 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)
The measured data during the tunnel construction should be selected as much as possible. The 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.



---

*Source: 5264072-2020-10-19.xml*