Jurnal 2

International Journal of Rock Mechanics & Mining Sciences 128 (2020) 104279
Contents lists available at ScienceDirect
International Journal of Rock Mechanics and Mining Sciences

journal homepage: http://www.elsevier.com/locate/ijrmms
Rockburst assessment in deep geotechnical conditions using true-triaxial

tests and data-driven approaches
Roohollah Shirani Faradonbeh a, *, Abbas Taheri a, Luis Ribeiro e Sousa b, c, Murat Karakus a
a
School of Civil, Environmental and Mining Engineering, The University of Adelaide, Adelaide, SA, 5005, Australia
b
State Key Laboratory for Geomechanics and Deep Underground Engineering, Beijing, China
c
University of Porto, Porto, Portugal
A R T I C L E I N F O A B S T R A C T
Keywords: Deep underground excavations in mining and civil engineering are subjected to high in-situ stresses which can
Rockburst maximum stress cause rockburst. Rockburst is an instantaneous release of a large amount of strain energy stored in rockmass that
Rockburst risk index can lead to injuries, deaths, and damage to infrastructures. Many studies have been done regarding rockburst,
True-triaxial test
however, there is no practical model to predict the stress level that rockburst occurs (i.e. maximum rockburst
Gene expression programming
Classification and regression tree
stress) and its related risk (i.e. rockburst risk index) based on real rockburst tests, and the main rock mechanical
properties. In this study, a comprehensive database of true-triaxial unloading tests on rocks having a wide range
of properties was compiled. The agglomerative hierarchical clustering (AHC) analysis was carried out on the
original database to evaluate the presence of natural groups and outliers. Then, the stepwise selection and
elimination (SSE) procedure were employed for dimension reduction of the problem and identifying the most
influential attributes on rockburst parameters. Afterward, two robust non-linear algorithms, including gene
expression programming (GEP) and classification and regression tree (CART) were used to develop the predictive
models for rockburst maximum stress and its risk index. The validation verification of the proposed models using
several indices proved the high prediction performance of the developed non-linear models. Finally, a parametric
analysis was carried out to evaluate the influence of each input parameter on the corresponding output. The
proposed models in this study are practical and do not require any presupposition about rockburst mechanism,
which makes them be used easily in practice by engineers at the design and progress stages of the underground
projects.
1. Introduction defined as a dynamic instability phenomenon of the surrounding rock

mass of an underground opening in the highly-stressed zone that is
Deep underground conditions can be characterised by the high level accompanied by a violent release of strain energy stored in the rock
of in-situ stresses, high groundwater pressure, high temperature, and mass.5 Rockbursts can occur either during the excavation (are known as
high brittleness of rocks.1,2 Consequently, the deep underground activ strainbursts) or after excavation (are known as impact-induced or
ities such as mining and tunnelling are usually subjected to different delayed rockbursts) based on the triggering factors in the form of a strip
nature-induced hazards during and after their closure. One of them is of rock slices, rock fall, ejection of rock fragments, with roaring sound.3,6
related to seismic events.3 When the deviatoric stress either in the Due to the violent and unexpected nature of this phenomenon, it may
confined rock mass or near an excavation reaches to the confined or lead to worker injury, damage to mine infrastructure and equipment,
unconfined rock mass strength, the seismic events are occurred and may and possibly economic loss of underground excavation (see Fig. 1).7–9
lead to unstable failure by a shear slip or shear rupture. The radiated Many factors including the rock mechanical properties, excavation ge
energy may then damage the excavation and cause rockburst damage.4 ometry, discontinuities, in-situ and mining-induced stresses, and con
Rockburst is the most severe disaster that threatens the safety of mining struction method, may affect the rockburst triggering.10 Hence,
operation and the stability of the excavation surface. Rockburst can be rockburst prediction and mitigation are one of the most challenging
* Corresponding author.
E-mail addresses: roohollah.shiranifaradonbeh@adelaide.edu.au (R. Shirani Faradonbeh), abbas.taheri@adelaide.edu.au (A. Taheri), sousa-scu@hotmail.com
(L. Ribeiro e Sousa), murat.karakus@adelaide.edu.au (M. Karakus).
https://doi.org/10.1016/j.ijrmms.2020.104279
Received 17 June 2019; Received in revised form 21 February 2020; Accepted 21 February 2020
Available online 27 February 2020
1365-1609/© 2020 Elsevier Ltd. All rights reserved.
R. Shirani Faradonbeh et al. International Journal of Rock Mechanics and Mining Sciences 128 (2020) 104279
problems in rock engineering. In last decades, a considerable investi the platen is removed abruptly from one surface of the sample that is
gation has been carried out on rockburst to understand more about its subjected to radial stress ðσ 3 Þ while the axial stress ðσ2 Þ is held constant
mechanism, prediction, and controlling methods theoretically and and σ 1 is increased gradually until the specimen fails. Due to the supe
experimentally.1,6 In terms of rockburst mechanism, many theories such riority of the modified true-triaxial testing system over the conventional
as the energy theory, the strength theory, the burst liability theory, the testing methods, many researchers recently have used this machine to
bifurcation theory, and the chaos theory have been proposed by re study the influence of different parameters such as moisture con
searchers to study the strain localisation and the stability of the rock tent,31–33 temperature27,34 tunnel axis stress,35 aspect ratio of the rock
mass.1,11,12 sample,36 and unloading rate37,38 on rockburst hazard. During the
From the viewpoint of rockburst prediction, more than one hundred rockburst tests, the rocks experience a severe failure at a specific stress
empirical criteria have been proposed by scholars mostly based on level which is known as rockburst maximum stress ðσ RB Þ (see Fig. 3). The
strength, strain, and strain energy parameters. Moreover, novel and correct prediction of this stress level for different rock types by consid
robust data-mining techniques including the supervised and unsuper ering the most important rock mechanical characteristics can help en
vised algorithms have been used recently by researchers to predict both gineers to recognise rockburst hazards in different in-situ stress
the rockburst occurrence (as a binary problem) and its intensity (as a conditions, to increase the stability of the underground structures as well
four-class problem).5,13 In terms of rockburst control, several techniques as for numerical studies. Furthermore, the predictor model can be used
such as the application of energy-absorbing bolts,14 ground pre to estimate bursting stress when rockburst testing facilities are not
conditioning (e.g. destress drilling, destress blasting, water injection), available. In this regard, He et al.16 proposed three models using mul
and alternative mining methods (e.g. pillarless mining and mining with tiple linear regression (MLR), artificial neural network (ANN), and
sacrifice galleries) have been suggested as the potential solutions for support vector machine (SVM) techniques to predict σRB and rockburst
rockburst mitigation.15 Rockburst also has been investigated by labo risk index (IRB ) based on the database compiled from true-triaxial tests.
ratory tests since these tests can provide useful information about In their study, ANN and SVM as two subsets of soft computing al
rockburst mechanism, the influence of different parameters on bursting gorithms outperformed the common multiple regression model and
behaviour, calibrating the numerical models, and identifying the stress were identified as possible solutions for studying such a complex prob
state prone to severe failure.6,16 Many attempts have been made by lem. Although these techniques can provide models with suitable ac
different researchers to simplify the process of failure using conven curacy, they are known as “black-box” methods i.e. their internal
tional testing methods including uniaxial and biaxial compression structure and calculations are not clear and easy to understand by a
tests,17–21 dynamic uniaxial compression tests,17,22 triaxial compression human. In addition, such algorithms may stick in local minimum during
unloading tests23–25 as well as the true-compression triaxial tests.23,26–28 the training process, and their outputs may not be very reliable. Above
Despite the many excellent efforts made so far, the common testing all, these techniques are not very practical since they cannot offer any
methods cannot properly simulate the stress state of the rock mass where mathematical or visual output to let the users apply them without using
rockburst occurs. As a matter of fact, rockburst usually occurs near the a code.39 It is also worth mentioning that many factors affect rockburst
excavated boundary while the stress state after the excavation is trans hazard and its related parameters, and selecting the most influential
formed from a triaxial equilibrium state ðσ1 > σ 2 > σ 3 Þ to a newly ones during the modelling, can affect the complexity, accuracy, and
redistributed state ðσ3 � 0; σ1 6¼ 0; σ2 6¼ 0Þ. In this new state, the more importantly the reliability of the developed models. To triumph
tangential stress ðσθ Þ increases progressively and may reach the ultimate over such complexities and open the black-box methods, the extraction
strength of the rock element and rockburst may occur (see Fig. 2). More of rules from the model, and the use of visualization techniques are
recently, a modified true-triaxial testing system was introduced by He recommended.39 Contrary to ANN and SVM techniques, there are other
et al.29 which can unload the pressure on one surface of the rock sample powerful algorithms such as genetic programming (GP), gene expression
and replicate the stress concentration in the surrounding rock mass of programming (GEP), and classification and regression tree (CART) that
the free face and measure the released kinetic energy during bursting. can provide more practical outputs. The successful application of these
This apparatus consists of several main components including the main algorithms has been reported by other researchers in mining and
hydraulic loading and unloading system, a data acquisition system that geotechnical engineering fields.40–43 Hence, it is necessary to use
continuously records the changes in forces and displacements, the state-of-the-art modelling techniques to address the mentioned diffi
high-speed cameras for measuring the kinetic energy of the ejected rock culties and develop new models for predicting rockburst maximum
fragments during bursting, an acoustic emission (AE) system to track the stress and its risk index based on field conditions. As it has been sum
process of damage accumulation during the test, and an infrared marised in Fig. 4, this study focuses on the following steps: 1) compiling
monitoring system to measure surface temperature of the sample and a database based on the true-triaxial unloading tests on different rock
inspect the internal damage.16,30 In true-triaxial testing apparatus as can types and performing a broad statistical analysis on it to create a ho
be seen from Fig. 3, the loads are applied simultaneously to the prismatic mogeneous database and to select the most influential parameters based
specimen in all three principal stress directions to reach a pre-defined on an appropriate strategy; 2) Developing genetic-based and decision
in-situ stress state, and the loads are kept constant on the specimen for tree-based models for the prediction of maximum rockburst stress ðσ RB Þ
a while to achieve a uniform stress distribution in the specimen. Finally, and rockburst risk index ðIRB Þ based on the selected input parameters; 3)
Fig. 1. Rock ejection and deformation of the supporting system due to strain bursting.3
2
Fig. 2. The schematic representation of the rock element stress state before and after tunnelling (modified from Su et al.27).
density, rock specific weight, mineral contents of rocks, loading and

unloading rates of the true-triaxial tests, rockburst maximum stress,
rockburst risk index, test duration and bursting mechanism. Considering
a circular shape for the tunnel crown, the stress concentration factor
equal to 2, and the specific weight of 27 kN/m3 for the overburden rock
mass, the rockburst critical depth (He ) was calculated by the following
equation:
He ¼ 18:52σ RB (1)
The rockburst risk index (IRB ) also was calculated for all the samples
through the following equation:44
H H
IRB ¼ ¼ 0:054 (2)
He σRB
He44 defined a new classification for IRB as shown in Table 1. Based
on this classification, a 56% of the tested samples have low IRB , 13% of
the samples have very high IRB , and the remained 31% of samples have
moderate to high IRB . Since all the foregoing parameters have not been
collected during the rockburst tests, there are some missing values in the
Fig. 3. Loading path for rockburst true-triaxial tests. database. To have a homogeneous database, the missing values (thirty
records) were eliminated from the primary database, and finally, the
validation verification of the developed models; and 4) conducting a results of 109 tests were considered for further analyses. Before devel
parametric analysis to assess the effect of input parameters on the cor oping any model, the presence of natural groups and outliers in the raw
responding outputs. database was evaluated using agglomerative hierarchical clustering
(AHC) analysis. In fact, the presence of outliers and natural groups can
2. Data collection and statistical analysis decrease the generality and liability of the developed models.13,45,46 The
AHC is the most common type of clustering techniques which is used in
In this study, a database containing information about the 139 earth sciences.45 The AHC follows a bottom-up procedure that itera
rockburst laboratory tests conducted on different rock types from 2004 tively creates the single object clusters and then these clusters are
to 2012 at the State Key Laboratory for Geomechanics and Deep Un merged into the larger clusters based on the similarity or dissimilarity
derground Engineering (SKLGDUE), China was compiled. The tested criteria. The common criterion for clustering is “distance”, and this
rock samples were gathered from the depth of 200 m–3375 m. This means that objects in the same cluster have the least distance from each
database consists of many parameters such as rock mechanical proper other, while objects in different clusters are at a great distance from one
ties, in-situ stresses, rock sample depth, rockburst critical depth, rock another. The process of cluster generating and merging is continued
3
Fig. 4. Process of rockburst assessment in this

study (UCS: uniaxial compressive strength; E:
Young’s modulus; ν: Poisson’s ratio; σh1 : hori
zontal in-situ stress; σ h2 : horizontal in-situ stress
in the face to be unloaded; σ v : vertical in-situ
stress; σ RB : rockburst maximum stress; IRB : rock
burst risk index; D: depth, ρ: density, K: hori
zontal pressure coefficient (ratio of average
horizontal stresses to the vertical stress due to
overburden), MLR: multiple linear regression;
VIF: variance inflation factor; R2 : coefficient of
determination).
until all the objects (datasets) are placed in a single cluster or the
Table 1
pre-defined termination condition is satisfied. For measuring the dis
Rockburst risk index classification.16
tance between the objects, the average-linkage function that measures
Rockburst risk index (IRB ) Class the average distance of any object of one cluster from an object of the
IRB < 0:6 Low other cluster was used to form the clusters:47,48
0:6 < IRB � 1:2 Moderate 1 XX
1:2 < IRB � 2:0 High dða; bÞ (3)
jAjjBj a2A b2B
IRB � 2:0 Very high
where A and B are two clusters with the sizes of jAj and jBj, respectively.
a and b are objects from the mentioned clusters and d is the squared
Euclidean distance between two objects.
4
Fig. 5 shows the dendrogram derived from the conducted clustering distinctive behaviour on the modelling process, and the subsequent
analysis by AHC. A dendrogram is a tool that represents the relative size analysis was carried out on the 107 data samples. Fig. 6 shows the his
of the calculated distances at which the objects and clusters are com togram and the descriptive statistics of all the eleven measured param
bined. The objects with the low squared Euclidean distance (high sim eters. This figure illustrates that there are specific ranges of values for
ilarity) are close together and vice versa. The X-axis shows the dataset the parameters in which the predictions are meaningful. Accordingly,
number and the Y-axis shows the rescaled value of the distance. To future predictions with a new database should be conducted only in
prevent Fig. 5 to be crowded and large, the numbers of the datasets have these ranges. Although the current database covers a great range of
been summarised on the X-axis. Clearly can be seen from Fig. 5 that the values for the parameters and includes different rock types, it is possible
whole 109 collected datasets were clustered into one distinct group to extend these ranges in the future by new information obtained from
between the rescaled distances of 0 and 5 except for two cases of 75 and true-triaxial tests to increase the generality of the developed models.
76 which were placed in the second group. By checking the database, it
was found out that the main parameter that caused to grouping is depth, 3. Methods and results
and the members of group 2 belong to the depth of 3375 m which are
known as outliers for the current database. Therefore, these two cases 3.1. Stepwise selection and elimination process
were removed from the database to avoid the influence of their
This section aims to do a systematic stepwise selection and elimi
nation (SSE) analysis to identify the most important parameters on the
outputs and reduce the complexity of the developed models. The process
of parameter reduction also is carried out using the variable pressure
tools of the robust data-mining techniques i.e. GEP and CART. There are
several critical statistical terms which have been used in this study for
the primary assessment of the database and are defined in the following.
Multicollinearity, a high correlation between the independent (predic
tor) variables, can be considered as one of the most prominent chal
lenges for multiple regressions. The existence of this phenomenon may
lead to developing an unstable regression model having high values for
variance and covariance coefficients.49 Variance inflation factor (VIF) is
a statistical index to quantify the extent of the multicollinearity between
the independent (input) parameters. This index is the ratio of model
variance considering several inputs to the variance of the model with a
single input parameter. The VIF lower than 10 shows the non-existence
of multicollinearity.50 Another important index is Sig. (2-tailed) or
p-value of the correlations. The Sig (2-tailed) represents the significance
of the correlation at a prescribed alpha level (5%). The Sig. (2-tailed)
should be less than or equal to 0.05 to reject the influence of chance
factor. The coefficient of determination (denoted by R2 ) is another sta
tistical measure for evaluation of the model performance. This index
interprets the proportion of the output (dependent) variable’s variance
that is predictable from the input (independent) variables. An R2 of 1
indicates that the regression predictions perfectly fit the data.50–52
In the current study, uniaxial compressive strength (UCS), Young’s
modulus (E), Poisson’s ratio (ν), horizontal in-situ stress (σh1 ), horizontal
in-situ stress in the face to be unloaded (σ h2 ), vertical in-situ stress (σv ),
depth (D), density (ρ), and horizontal pressure coefficient (K) are known
as the input parameters for the maximum rockburst stress (σ RB ), while
all the mentioned parameters are considered as inputs for the rockburst
risk index (IRB ). The SPSS software package 25.0 was used for per
forming the statistical evaluations. Initially, the database was fed to the
software, and the Person’s correlation coefficient (r) between the input
parameters as well as between the inputs and the corresponding outputs
was calculated. Table 2 lists the calculated correlation values. As can be
seen from this table, all the inputs significantly correlating with σ RB (i.e.
Sig: ð2 tailedÞ � 0:05), while D (depth) with the Sig: > 0:05 and low
correlation coefficient (r ¼ 0:128) was removed from the input pa
rameters for further modelling of IRB . The elimination of parameters
does not show that they have not any influence on the output, but simply
it means that the effect of those parameters will be minimum in pre
dicting the output. As an initial multicollinearity assessment between
input parameters, no one of the correlations exceeds from the condition
of r > 0:90.53 However, these input parameters may show multi
collinearity when a combination of them are used as regressors in MLR.
Based on the above analysis, all the inputs (except parameter D for IRB )
were retained for multiple linear regression (MLR). The MLR models
with the possible multicollinearity were developed separately using the
Fig. 5. Dendrogram resulting from agglomerative hierarchical clustering selected parameters for both σRB and IRB .
(AHC) analysis.
5
Fig. 6. Histogram of the collected parameters along with the descriptive statistics.
Table 3 shows the model summary, calculated coefficients, and the calculations in this process are like the forward selection and backward
statistical indices for evaluating the developed MLR models. In this procedure.54 Table 4 summarises the results of the stepwise selection
stage, according to Fig. 4, several conditions including VIF < 10, Std: and elimination procedure carried out using the algorithm provided in
error � Coeff: ðBÞ, and Coeff: ðBÞ 6¼ 0 were checked for different inputs the SPSS 25. In this algorithm, the parameters enter the model if the
to retain them for further evaluations. Considering Table 3, for rockburst probability (significance level) of its F value is less than the Entry value
maximum stress (σ RB ), the parameters of K and ρ have VIF > 10 and t- (i.e. 0.05) and are eliminated if the probability is greater than the
significance higher than 0.05, respectively, which shows that the effect Removal value (i.e. 0.100). Entry must be less than Removal, and both
of these parameters on the σRB is insignificant. Therefore, these pa values must be positive. As given in Table 4, during the process of se
rameters were removed for further modelling of σRB . About the rock lection and elimination, the correlation coefficient (R) between the
burst risk index (IRB ), all the VIF values for inputs are less than 10, but measured output and the predicted one was increased from 0.884
the t-significance values of the UCS, σv , σ h1 , and σh2 are higher than 0.05. (model 1) to 0.910 (model 3) for σRB and from 0.754 (model 1) to 0.821
Thus, these parameters also were removed from the input set of IRB . In (model 5) for IRB , respectively. In other words, the parameters of UCS, E,
the next step, two stepwise selection and elimination procedures were and σ v can explain 82.8% (R2 ¼ 0:828) variations in σRB . As such, the
performed using the selected inputs for each dependent parameter. In parameters of σRB , K, E, ν, and ρ can explain 67.4% (R2 ¼ 0:674) vari
this procedure, a parameter which is entered in the model at the initial ations in IRB . Thus, these parameters were known as the most influential
stage of selection may be removed at the later stages. In fact, the ones among the initial inputs to describe the rockburst parameters. The
6
Fig. 6. (continued).
regression coefficients and the collinearity statistics of the best (GEP) as a meta-heuristic algorithm and the classification and regression
SSE-based models are shown in Table 5. In both models, the VIF factor tree (CART) as a subset of decision tree algorithms were selected for
that shows the multicollinearity is lower than 10, the t-significance is discovering the non-linear latent relationships with more accuracy and
lower than 0.05, and the Std. error values are lower than the regression lower estimation error. These algorithms despite the various datamining
coefficients which show the reliability of the SSE process in identifying and soft computing techniques such as ANNs, SVM, etc. can provide
the most influential parameters. practical and easy to use outputs for the engineers and the researchers
Considering the above analyses, the agglomerative hierarchical when the true-triaxial testing machine is not available. A summary of the
clustering (AHC) accompanied by the stepwise selection and elimination modelling procedure by these techniques is presented in the following
(SSE) method could provide a homogeneous rockburst database by sections.
removing the outliers and decreasing the dimensionality of the problem.
This process also can be useful for the complexity reduction of the next 3.2. Non-linear regression analysis
predictive models by applying a few input parameters. Due to the high
non-linear and complex nature of rockburst hazard16,55,56 there is a need Non-linear regression (NLR) attempts to find a function which is a
to use the non-linear data-mining algorithms to provide more accurate non-linear combination of the input parameters using a method of suc
predictive models for rockburst parameters. To do so, two robust cessive approximation.57,58 In geoscience, most of the dependent pa
data-driven approaches including the gene expression programming rameters show a non-linear relationship with the related influential
7
Table 2
Pearson’s correlation (r) coefficient between different parameters.
UCS ρ E ν σv σh1 σh2 D K σRB IRB
UCS r 1
Sig.
ρ r 0.827 1
Sig. 0.000
E r 0.796 0.829 1
Sig. 0.000 0.000
ν r 0.555 0.438 0.594 1
Sig. 0.000 0.000 0.000
σv r 0.660 0.491 0.458 0.367 1
Sig. 0.000 0.000 0.000 0.000
σh1 r 0.737 0.523 0.515 0.420 0.846 1
Sig. 0.000 0.000 0.000 0.000 0.000
σh2 r 0.713 0.549 0.648 0.460 0.685 0.803 1
Sig. 0.000 0.000 0.000 0.000 0.000 0.000
D r 0.450 0.065 0.247 0.563 0.450 0.554 0.430 1
Sig. 0.000 0.509 0.010 0.000 0.000 0.000 0.000
K r 0.667 0.717 0.655 0.362 0.629 0.712 0.791 0.005 1
Sig. 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.961
σRB r 0.820 0.664 0.668 0.484 0.884 0.838 0.762 0.413 0.707 1
Sig. 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
IRB r 0.617 0.539 0.624 0.255 0.619 0.645 0.618 0.128 0.656 0.754 1
Sig. 0.000 0.000 0.000 0.008 0.000 0.000 0.000 0.188 0.000 0.000
* Correlation is significant at the 0.05 level (2-tailed).
Table 3
Summary of the MLR models for σ RB and IRB with the selected parameters (with multicollinearity).
Model summary
Dependent parameter r r-square Adjusted r-square Std. error of the estimate
σRB 0.952 0.906 0.897 19.385

IRB 0.832 0.692 0.664 0.470
Rockburst maximum stress (σRB )

Parameters Unstandardized coefficients standardized coefficients t Sig. Collinearity statistics
B Std. Error Beta Tolerance VIF
(Constant) 71.247 20.947 3.401 0.001

UCS 0.380 .093 0.354 4.097 0.000 0.130 7.719
ρ 1.844 1.794 0.086 1.028 0.307 0.138 7.271
E .400 0.193 0.144 2.072 0.041 0.200 5.005
ν 125.188 50.304 0.132 2.489 0.015 0.344 2.908
σv 0.594 0.061 0.572 9.688 0.000 0.278 3.592
σh1 0.650 0.246 0.259 2.639 0.010 0.101 9.926
σh2 0.664 0.284 0.186 2.342 0.021 0.154 6.490
D 0.060 0.016 0.333 3.750 0.000 0.123 8.146
K 15.291 5.571 0.303 2.745 0.007 0.080 12.556
Rockburst risk index (IRB )

Parameters Unstandardized coefficients standardized coefficients t Sig. Collinearity statistics
B Std. Error Beta Tolerance VIF
(Constant) 2.766 0.282 9.816 0.000

UCS 0.001 0.002 0.065 0.415 0.679 0.129 7.733
ρ 0.107 0.042 0.375 2.528 0.013 0.144 6.926
E 0.020 0.005 0.550 4.286 0.000 0.193 5.194
ν 3.406 0.932 0.268 3.653 0.000 0.589 1.698
σv 0.003 0.002 0.223 1.535 0.128 0.151 6.629
σh1 0.006 0.005 0.185 1.305 0.195 0.158 6.335
σh2 0.011 0.006 0.231 1.807 0.074 0.194 5.167
K 0.221 0.083 0.327 2.662 0.009 0.210 4.751
σRB 0.011 0.002 0.805 4.689 0.000 0.108 9.284
parameters. So, the non-linear regression analysis has been widely used common NLR technique suffers from several significant drawbacks. In
by researchers in the last decades.41,59,60 The NLR technique is capable NLR, there is no a closed-form and holistic mathematical structure be
of accommodating a broad range of functions including exponential, tween the dependent and the independent parameters as there is in
power, logarithmic, sigmoid, logistic, trigonometric, Gaussian, etc. that multiple linear regression (MLR), while the choice of the model struc
boosts the process of function finding. Another advantage of the NLR is ture is a crucial task to obtain the best solution. In addition, the selection
the efficient use of data, i.e. it can provide reasonable estimates of the and utilizing the suitable mathematical functions from the large library
unknown parameters for a comparatively small data. However, the of functions need an iterative optimization procedure that is not possible
8
Table 4 advance by the user.64 By inspiring from the Darwinian principle of

MLR models developed based on SSE method (without multicollinearity) for σ RB “Survival of the Fittest”65 and the natural evolution, a new subset of soft
and IRB . computing was introduced as the evolutionary algorithm (EA). Gener
Rockburst maximum stress (σRB ) ally speaking, EAs work with a randomly generated population of in
dividuals which are then improved using a group of genetic operators (e.
Stepwise selection and elimination Model summary
g. mutation, crossover and reproduction) and finally, the solutions are
Model Variables Variables Removed r r- Std. encoded into the specific forms such as binary strings in genetic algo
entered square error
rithm. The main differences between EAs are related to the method of
1 σv UCS; E; ν; σh1 ; σh2 ; 0.884 0.781 28.412 presenting the solutions, genetic operators, selection mechanism, and
D the performance measurement method.62,66 Gene expression program
2 UCS E ν, σh1 ; σh2 ; D 0.938 0.880 21.102 ming (GEP)67 is a well-known evolutionary algorithm that inherits two
3 E ν 0.910 0.828 20.399 essential features from its siblings i.e. the use of simple, fixed-length, and
σh1 ; σh2 , D
linear chromosomes with different shapes and sizes from genetic algo
Rockburst risk index (IRB ) rithm (GA) and the expression tree (ET) structure from genetic pro
Stepwise selection and elimination Model summary gramming (GP) that improves the robustness of GEP for solving the
Model Variables Variables Removed r r- Std. non-linear problems.68 The main entities of GEP algorithm are termi
entered square error
nal set (input parameters and constant values), function set (e.g. þ ; ;
1 σRB ρ, E, ν, K 0.754 0.568 0.564 � ; �), fitness function (for evaluating the generated solutions), and
2 K ρ, E, ν 0.773 0.598 0.591 genetic operators (mutation, inversion, transposition, and
3 ν ρ, E 0.784 0.615 0.604 recombination).
4 E ρ 0.805 0.649 0.635 A flowchart detailing the GEP modelling procedure is shown in
5 ρ 0.821 0.674 0.658 Fig. 7. In summary, GEP generates a population of chromosomes (so
lution/individual) by combining the user-defined terminals and func
tions. These chromosomes follow a bilingual and unequivocal
in common NLR modellings. Accordingly, the researchers may have to expression system that is called Karva language.69 The chromosomes
use numerical optimization algorithms to find the best-fitting parame have a specified number of genes (sub -ETs) which are linked together
ters but still, there is a need to define the starting values for the unknown using a linking function (e.g. “=” in Fig. 7 that links two genes of a
parameters in these methods. Inappropriate assigning the starting values chromosome). Each gene contains two parts of head and tail that the
may cause to getting caught in the local minima rather than finding the terminals (inputs) and functions (mathematical functions) are placed in
global minimum that introduces the least squares estimates.57,58,61 For them and the genetic operators are applied to these areas to modify the
these difficulties, the researchers prefer to use a non-linear regression solutions. To have a quick understanding regarding the built-in mathe
form that has been used successfully in similar applications. Hereupon, matical equations of chromosomes, the Karva coded programs are then
the application of intelligent algorithms is needed to cope with these parsed into ETs. Then, the mathematical form of the programs is
issues. In the following sections, the process of rockburst assessment extracted from ETs and their fitness is evaluated by a fitness function. If
using two robust non-linear techniques comprising the gene expression the stopping condition(s) such as reaching to a specific number of iter
programming (GEP) and classification and regression tree (CART) are ations or the desired fitness value is not met, the selected chromosomes
explained. are replicated into a new generation, and the remained ones undergo a
modification process using the genetic operators. The above process is
3.2.1. Rockburst assessment using GEP-based models repeated and finally, the best solution (predictive model) describing the
Soft computing is the relatively new branch of data-mining methods relationship between the input and output parameters is found. More
and can be considered as an alternative to the prevalent hard computing details about the mechanism of genetic operators and GEP algorithm can
methods for solving the real-world problems.62,63 Soft computing tech be found in Ferreira.69 In the current study, the selected inputs from the
niques have been successfully employed in mining, rock mechanics, and SSE analysis were considered as terminal sets to formulate the rockburst
geotechnical problems but despite their good performance, they cannot parameters nonlinearly as follows:
generate practical equations, and their structure needs to be assigned in
Table 5
The statistical parameters of the best SSE-based MLR models.
Rockburst maximum stress (σRB )
Model Unstandardized coefficients Standardized coefficients t Sig. Collinearity statistics
B Std. error Beta Tolerance VIF
(Constant) 6.019 3.080 1.954 0.053

σv 0.650 0.046 0.626 14.196 0.000 0.552 1.811
UCS 0.301 0.069 0.281 4.336 0.000 0.256 3.912
E 0.437 0.152 0.158 2.881 0.005 0.358 2.794
Rockburst risk index (IRB )

Model Unstandardized coefficients Standardized coefficients t Sig. Collinearity statistics
B Std. error Beta Tolerance VIF
(Constant) 2.780 0.278 10.006 0.000

σRB 0.008 0.001 0.601 6.749 0.012 0.406 2.463
K 0.158 0.062 0.233 2.556 0.000 0.388 2.578
ν 3.501 0.922 0.276 3.797 0.000 0.612 1.634
E 0.018 0.004 0.497 4.280 0.000 0.239 4.189
ρ 0.091 0.032 0.319 2.830 0.006 0.254 3.933
9
Fig. 7. Process of function finding using GEP algorithm.
σ RB ¼ f ðUCS; E; σ v Þ (4) The RMSEi varies between 0 and infinity, with 0 corresponding to the
ideal. Since the process of selection in GEP algorithm is based on the
IRB ¼ f ðσ RB ; K; υ; E; ρÞ (5) increase of fitness, Eq. (6) cannot be used directly. Thus, the following
expression was used for fitness function which obviously ranges between
The rockburst database was divided randomly into training and
0 and 1000, with 1000 corresponding to the ideal:
testing subsets. The training set (80% of the database) was used to train
the model and discover the relationship between inputs and outputs, and 1
RMSE’i ¼ 1000 � (7)
the remaining datasets were used to validate the performance of the 1 þ RMSEi
proposed models. It should be noted that the influence of using different
On the other hand, to apply the parsimony pressure on future
groups of training and testing datasets were also evaluated on the ac
models, overall fitness was defined as:
curacy of the models. However, no noticeable change in the results was
� �
observed. For evaluating the generated solutions during the GEP 1 Smax Si
RMSE’’i ¼ RMSE’i � 1 þ � (8)
modelling, it is necessary to use a fitness function. As mentioned in 5000 Smax Smin
section 3.1, to propose models with lower complexity, it is possible to
apply variable pressure tools to compress the developed models as much where Si is the size of the GEP program, Smax and Smin are the maximum
as possible by eliminating the parameters which have lower importance and minimum program sizes which are calculated by the following
in a non-linear structure. To this end, the root mean squared error equations:
(RMSE) with parsimony pressure was applied to the GEP models of σRB
Smax ¼ Gðh þ tÞ (9)
and IRB :70 The RMSEi of a chromosome (solution) i is calculated by the
following equation: Smin ¼ G (10)
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
u X
u1 n �2 where G is the number of genes, and h and t are the head size and tail
RMSEi ¼ t Pij Tj (6)
n j¼1 size, respectively.
A group of trigonometric and straightforward mathematical func
where Pij is the predicted value by the chromosome i for the dataset j, tions i.e. f þ; ; *; =; √; Ln; 2;
^ 3;
^ 1=3;
^ sin; cos; tang were selected as
and Tj is the measured value for dataset j. the function set based on the previous non-linear studies using GEP
10
algorithm.46,71,72 The other GEP parameters including the number of

chromosomes, head size, the number of genes, and the values of genetic
operators were changed for different runs to obtain the best solution in
such a way that provides not only high accuracy but also less
complexity. Table 6 presents the architecture of the obtained GEP
models for both rockburst maximum stress (σ RB ) and rockburst risk
index (IRB ). By applying the parsimony pressure to the models, the
density parameter (ρ) was identified intelligently as the low-impact
parameter in the non-linear form of IRB . Therefore, this parameter was
removed by GEP automatically during modelling and the number of
inputs for IRB decreased from 5 to 4. About σRB , the GEP algorithm
identified the three inputs of UCS, E, and σv as the influential parameters
for modelling as formerly proved by SSE analysis. The ability of GEP in
identifying the low-influence parameters and excluding them during
modelling can be considered as an internal sensitivity analysis that
distinguishes GEP from other soft computing techniques. Fig. 8 displays
the variations of the coefficient of determination (R2 ) during 5000
generations (iterations) in both training and testing stages of GEP
modelling for rockburst parameters. According to this figure, after a few
numbers of generations (less than 1000), a rapid increase of R2 for the
generated solutions can be seen which shows the high speed and high
capability of GEP algorithm in function finding. From the generation
1000 to 3500, a gentle enhancement in the quality of solutions are
visible, and finally, the algorithm converges into an optimum value and
its value almost remains constant to reach the stopping condition (i.e.
the pre-defined number of generations: 5000). The obtained R2 values
for training and testing stages of σRB are 0.9266 and 0.9398, respec
tively, while the foregoing values are 0.8824 and 0.9459, respectively
for IRB . Fig. 9 shows the correlation of the experimentally measured
values of σRB and IRB versus the predicted ones by the constructed GEP
models for training and testing data groups. As seen, the data points
Table 6
The architecture of the GEP and CART models.
GEP parameter Setting
σRB IRB
Terminal set UCS; E; σv σRB ; K; υ; E; ρ

Excluded parameter – ρ
Function set þ; ; *; =; √; Ln; 2; ^ p
^ 3; ffi
3 ; sin;
Fig. 8. Improvement of R2 during GEP modelling for (a) σRB and (b) IRB .
cos; tan
Population size 90 100 have almost a uniform distribution around the fitted lines in both
Generation number 5000 5000 GEP-based models which show the goodness-of-fit of the models. The
Head size 9 9
developed models and their performance are discussed in more details in
Number of genes 3 2
Linking function Multiplication ( � Multiplication ( � sections 4 and 5. Eventually, the mathematical forms of the proposed
) ) GEP models for σRB and IRB were extracted from their K-expression and
Fitness function RMSE’’i RMSE’’i ETs as Equations (11) and (12). To avoid the prolongation of the paper,
Parsimony pressure Yes Yes the ETs and their K-expressions have not presented here.
Mutation rate 0.01 0.04 � pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi �� pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi � �
Inversion rate 0.1 0.1 σ RB ¼ 3 σv þ EsinðE σv Þ þ E 3
E þ σv þ EsinðEÞ LnðLnðσv ÞþUCS
Transposition 0.1 0.1
One-point recombination 0.3 0.3 (11)
Two-point recombination 0.3 0.3
6 ffiffi
p p ffiffi �
Gene recombination 0.1 0.1 e Eþ ν K 4
IRB ¼ (12)
CART parameter Setting ðν þ KÞðE νÞ Lnðσ RB Þ
σRB IRB
Initial inputs UCS; E; σv σRB ; K; υ; E; ρ 3.2.2. Rockburst assessment using classification and regression tree (CART)
Excluded parameter – ρ Decision tree as a powerful subset of data-mining techniques has
Minimum number of cases for parent 3 3 been used in different real-world applications for different aims such as
node decision making, classification, prediction, pattern recognition, etc.73,74
Minimum number of cases for child 1 1 A decision tree is a tree comprising a root node (i.e. a parameter that can
node
provide maximum degree of discrimination), some internal nodes rep
Minimum change of impurity level 0.0005 0.0003
Maximum tree depth 6 5 resenting input parameters, branches which link the nodes together and
Number of intervals 10 10 contain the binary questions regarding the internal nodes, and some leaf
Impurity measure LSD LSD nodes representing the solutions (predicted value or a specific class of
Total number of nodes 27 33
the dependent parameter). Each path from the root node to a leaf node
11
can be summarised as a rule that this feature makes the decision tree to
be known as a rule-based algorithm.75 Based on the type of dependent
parameter, i.e. being continuous or categorical, the established tree
structure is nominated as regression tree (RT) or classification tree (CT),
respectively. The decision tree has several subgroups such as ID3,76
C4.5, C5.0, CART, CHAID, Exhaustive CHAID, and QUEST75 which have
been used for different aims by scholars.42,77 Among these techniques,
the CART algorithm introduced by Breiman et al.78 has several advan
tages that distinguish it among other decision tree algorithms. This al
gorithm, despite the parametric statistical techniques (e.g. regression
analyses), is inherently non-parametric (rule-based), i.e. no assumption
is made with the distribution of values of the independent parameters.
On the other hand, CART can handle the highly skewed (multimodal)
quantitative data as well as the qualitative parameters with ordinal or
non-ordinal structures.43,78 In this algorithm, it is not necessary to
eliminate the multicollinearity between the independent parameters.
Moreover, CART algorithm can be applied on a database with no ho
mogeneity. CART also can handle the existence of outliers in the raw
database by isolating them into a separate node. Because of the
mentioned advantages, flexibility, and practical output (tree structure)
of this algorithm, it was used in this study for the prediction of rockburst
parameters obtained from true-triaxial tests. As a matter of fact, since
the output parameters in this study (i.e. σRB and IRB ) are continuous, the
aim is to develop two regression trees (RTs) for each parameter. The
process of RT building in CART algorithm focuses mainly on the three
following components: (1) a group of questions in the form of X � a?
where X is an input parameter and a is a constant value in a range that
the parameter X varies. In CART, the response to this type of question is
“yes” or “no”; (2) the best split on a parameter is determined using a split
criterion; (3) calculation of summary statistics for internal nodes. The
goal in CART modelling is to create sub-nodes (children) which are more
homogeneous and purer than parent nodes based upon the reduction in
impurity or improvement score. The term “pure” is related to the values
of given parameter i.e. in the complete pure node, all cases have a
similar value of the splitting parameter and consequently, the node’s
variance equal to zero. This issue is compared for all the input param
eters and the best improvement is chosen for splitting. This procedure
continues until one of the stopping conditions is triggered.78 In CART,
the least squared deviation (LSD) is used as an impurity measure. The
LSD function for splitting a parent node t into two newly generated
sub-nodes tLðeftÞ and tRðightÞ can be calculated using the following equa Fig. 9. Measured vs. predicted values using the GEP models in training and
tion:78,79 testing stages for: (a) σRB and (b) IRB .
1 X 1 X 1 X
ΦðtÞ ¼ R2 ðtÞ pL R2 ðtL Þ pR R2 ðtR Þ ¼ ½yi yðtÞ�2 pL ½yi yðtL Þ�2 pR ½yi yðtR Þ�2 (13)
NðtÞ iεt NðtL Þ iεtL NðtR Þ iεtR
the number of cases in the sub-nodes resulting from the best splits is less
where R2 ðtx Þ is the weighted variance related to the sub-node (child) tx , than pre-defined minimum child size. The stopping criteria used in this
pL is the proportion of cases in parent node t which are classified in the study for CART models are tabulated in Table 6. These criteria and their
left node (tL ), pR is the proportion of cases in parent node t which are corresponding values were obtained in such a way that the results pro
classified in the right node (tR ), Nðtx Þ is the number of cases classified in vide a good trade-off between the prediction accuracy of regression trees
sub-node tx (xεfR; Lg), yi is the value of the objective parameter for the and their complexity (dimension). All these settings also prevent the
case i, yðtÞ is the mean value of parent node, and yðtx Þ is the mean value models from getting stuck in over-fitting problems. The use of a high
of the sub-node tx . number of maximum tree depth can lead to producing a large tree
The best split is obtained by maximizing the ΦðtÞ showing the structure with high complexity that makes it complicated to use in
reduction of impurity of an RT model. This splitting process leads to practice. Additionally, a maximum number of intervals equal to 10 was
creating a tree structure based on several “if-then” rules that make it considered for both models to let the model break down the initial min-
easy to represent. The splitting process proceeds until each leaf node max range of each input parameter to different ranges during the
meets at least one of the stopping criteria. The stopping criteria include: splitting process. In this study, the CART models for rockburst param
(1) reaching the maximum tree depth; (2) the number of cases (datasets) eters (σRB and IRB ) were developed using a code written in MatLab
in the terminal node is less than the predefined minimum parent size; (3) R2019a software environment. To have the same modelling conditions
for further assessments, the training and testing datasets used for GEP
12
Fig. 10. Regression tree model generated by the CART algorithm for σRB .
were fed again to the CART algorithm. According to Table 6, for rock using this technique were depicted in Fig. 12. As it is clear from this
burst risk index (IRB ) model, like the GEP-based one, the density (ρ) figure, the CART algorithm could predict σRB and IRB with high accuracy
parameter has been excluded from the model since the CART benefits like GEP models both in training and testing stages, and the data samples
from an internal principal component analysis (PCA) that enables it to have an appropriate scatter around the fitted line. The developed models
consider most influential parameters. Figs. 10 and 11 demonstrate the and their performance are discussed in more details in sections 4 and 5.
constructed RTs for σ RB and IRB using CART algorithm, respectively. The
tree model of σRB contains 27 nodes and starts with UCS as a root node, 4. Validation verification
while the IRB model has 33 nodes that starts with σRB as the root node.
The extraction of final predicted values of σRB and IRB from these Based on a logical hypothesis,64,80 there is a good correlation be
regression trees is an easy task. For instance, in Fig. 10, consider the tween the measured and predicted values of a dependent parameter
experimentally measured values of 82.7 MPa, 24.3 GPa, 114.6 MPa, and when the absolute correlation coefficient (jRj) is greater than 0.8, and
108.6 MPa for UCS, E, σv , and σRB , respectively; by tracking the asso the error indices such as root mean squared error (RMSE) and mean
ciated tree structure from the root node (i.e. node 1: UCS) and the path absolute error (MAE) are at low values. There are no definite values for
UCS � 40:45, σv < 169:455, σ v � 63:3, UCS < 151:3, σ v < 143:5, and R2 , RMSE, and MAE indices for all applications, and their appropriate
σ v < 118:4, the tree reaches to the leaf node 22 that predicts the value of values usually in different sciences depend on the range of variations of
116.1 for σ RB . The same process can be done for IRB as well. As stated at the response parameter as well as the sensitivity of the problem.
the beginning of the current section, CART is a rule-based technique i.e. Furthermore, to assess the performance of the developed models in
its internal calculations can be expressed clearly for user/reader by depth, new validation indices have been proposed by other researchers.
means of “if-then” rules, and this prominent feature makes the CART as a Golbraikh and Tropsha,81 defined two indices of k and k’ to validate the
white-box technique unlike the soft computing techniques such as models on testing datasets. In addition, Roy and Roy82 proposed an in
ANNs, SVMs that suffer from the lack of this characteristic. Tables A1 dicator called Rm along with another related parameter of R2O to check
and A2 in Appendix A show the generated rules for each node in the the predictability of the proposed models. The corresponding values of
developed CART models for σ RB and IRB . To have a primary insight
the k, k’ , Rm , and R2O can be calculated based on the measured (hi ) and
regarding the prediction power of the CART models, the scatter plots of
predicted (ti ) values of the output parameters (here are σ RB and IRB ). The
measured values of rockburst parameters versus the predicted ones
mathematical expressions of the above indices and their threshold
13
Fig. 11. Regression tree model generated by the CART algorithm for IRB .
values are listed in Table 7. Taking into account the recommendations, error, and subsequently have this capability to be used in practical ap
at least a slope of the regression lines (i.e. k or k’ ) through the origin plications. However, the GEP models of σ RB and IRB outperformed the
should be close to 1, while k is the slope of the regression line when hi is CART models and have slightly better performance. The next section
plotted versus ti , and k’ is the regression line when ti is plotted versus aims to do a parametric analysis on the selected models (i.e. GEP
hi .81 The squared correlation coefficient between the predicted and models) to appraise the effect of the variation of input parameters on the
measured values (R2O ) should be close to 1. The Rm , then, can be calcu predicted values.
lated by R and R2O values, and a threshold of > 0:5 is recommended for
this index to introduce a model as valid. The foregoing indices were 5. Parametric analysis
calculated for the developed GEP-based and CART-based models, and
To investigate the influence of each input parameter on the predicted
their values are listed in Table 7. Indeed, the indices of R, k, k’ , Rm , and
values of the corresponding output, a parametric analysis was carried
R2O were used to verify the validity of the models in testing stage as
out based on the selected GEP models for σRB and IRB . This analysis also
recommended by other researchers.83,84 Then, the statistical indices of
can be another validation for the GEP models by evaluating how well the
R2 , RMSE, and MAE were calculated to compare the prediction perfor results (predicted values) agree with the physical behaviour of the
mance of the GEP and CART models for σ RB and IRB based on testing
rockburst parameters. To do so, the desired independent parameter
datasets to select the best models. It can be observed from Table 7 that should be varied within its range of values, while other independent
the both proposed models in this study satisfy all the required condi
parameters are constant in their averages. Figs. 13 and 14 plot the
tions, and this guarantees that the derived models are strongly credible i. variation of input parameters against the predicted values for rockburst
e. the results are not based on chance factor. In addition, comparing the
parameters. As it is seen in Fig. 13, the σ RB increases monotonically in a
R2 , RMSE, and MAE values of GEP and CART models show that both non-linear fashion with UCS and σv . This result is expected since with the
proposed models have a high degree of accuracy and low estimation increase of UCS, the capacity of the rock to accumulate the strain energy
14
increases, and finally, bursting occurs at a higher stress level violently.18

On the other hand, the in-situ stresses, especially the vertical in-situ
stress (σv ) are increased in a linear or non-linear relationship with
depth85 and subsequently, due to a high geo-stress state in deep condi
tions, the σ RB is enhanced. However, there are many fluctuations in σ RB
values with the increase of Young’s modulus (E) of rocks, but in general,
an increment trend can be seen. It should be mentioned that a parameter
may do not show a meaningful relationship solely with the output
parameter, while it can be an influential component in a combination of
other parameters in a non-linear form. As mentioned in the GEP
modelling section, during the modelling procedure, by applying the
variable pressure coefficient, excluding any of the selected three pa
rameters (i.e. UCS, E, and σ v ) from the modelling procedure did not
improve the accuracy and complexity of the model.
Regarding IRB , a non-linear decreasing trend can be observed for its
values with all input parameters of Young’s modulus (E), Poisson’s ratio
(ν), horizontal pressure coefficient (K), and rockburst maximum stress
(σRB ). As can be seen from Fig. 14, with the increase of E until 20 MPa,
the rockburst risk index is decreased suddenly but it remains almost
constant with a further increment of E. Moreover, with the increase of
Poisson’s ratio (ν) in its range of values, the risk value decreases from
0.473 to 0.40 that according to Table 1, the risk of rockburst occurrence
is low. Hence, it seems that the variation of ν has no significant influence
on rockburst risk. However, it is still necessary to do more tests on rocks
with a greater range of ν to check its influence on risk parameter.
Generally, the risk of rockburst occurrence for rocks with low strength
(or lower σRB ) which are in low depth (or higher K) is higher than high-
strength rocks in deep conditions. From the results displayed in Figs. 13
and 14, several non-linear equations between rockburst parameters
(σRB , and IRB ) and their related input parameters (except for σRB E and
IRB E) are extracted as follows:
σ RB ¼ 19:911LnðUCSÞ þ 10:636; R2 ¼ 0:9974 (14)
σ RB ¼ 0:0007σ2v þ 0:7162σv þ 41:092; R2 ¼ 0:9877 (15)
IRB ¼ 0:49e 0:535υ

; R2 ¼ 0:9997 (16)
IRB ¼ 0:0186K 3 þ 0:2124K 2 0:7981K þ 1:1997; R2 ¼ 0:9521 (17)
Fig. 12. Measured vs. predicted values using the CART models in training and IRB ¼ 1:2524σRB0:244 ; R2 ¼ 0:9859 (18)
testing stages for: (a) σ RB and (b) IRB .
Table 7
Statistical indices for the external validation of the developed models.
Item Formula Threshold σRB IRB
GEP CART GEP CART

Pn
1 ðhi hi Þðti ti Þ R > 0:8 0.969 0.957 0.972 0.943
R ¼ qffiffiffiffiffiffiffiffiffiffiffii¼1
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
Pn 2 Pn 2
i¼1 i ðh hi Þ i¼1 ðti ti Þ
Pn
2 ðhi ti Þ 0:85 < k < 1:15 0.934 0.931 0.962 0.981
k ¼ Pi¼1 n 2
i¼1 hi
Pn
3 i¼1 ðhi ti Þ 0:85 < k’ < 1:15 1.046 1.040 1.014 0.971
k’ ¼ P n
i¼1 ti
2
qffi�ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi�ffi
4 Rm ¼ R2 ð1 �R2 R2 �Þ Rm > 0:5 0.763 0.698 0.739 0.596
O
Should be close to 1 0.987 0.986 0.997 0.999
Pn O 2
ðt i h Þ
R2O ¼ 1 i¼1
Pn
i
2
;
i¼1 ðti ti Þ
hOi ¼ k ti
5 R2 Should be close to 1 0.939 0.916 0.946 0.889
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
6 1X n Should be minimum (based on output range) 14.249 16.426 0.195 0.273
RMSE ¼ ðhi ti Þ2
n i¼1
7 1 Xn Should be minimum (based on output range) 9.803 8.041 0.136 0.187
MAE ¼ jhi ti j
n i¼1
hi : measured output; ti : predicted output.
15
based on the collected datasets and a specific range of values for different
parameters. So, for future applications, if the input parameters are out of
these ranges, the proposed models should be adjusted again.
6. Summary and conclusions
As a catastrophic geohazard, rockburst threatens the safety of

workers and infrastructures in deep geotechnical conditions. In this
study, considering the importance of the stress level that rockburst oc
curs for different rock types in real stress circumstances, two important
rockburst parameters including the maximum rockburst stress (σ RB ) and
rockburst risk index (IRB ) were formulated using the information ob
tained from true-triaxial unloading tests and two robust data-driven
approaches. A comprehensive strategy was applied to the compiled
database using the correlation analysis, the agglomerative hierarchical
clustering (AHC) technique, and the stepwise selection and elimination
(SSE) procedure to provide a homogeneous database free from any
outliers, natural groups, and especially, to identify the most influential
parameters on σRB and IRB . Then, new non-linear models were developed
using robust algorithms of gene expression programming (GEP) and
classification and regression tree (CART). Finally, a parametric analysis
was conducted to study the variation of σ RB and IRB with the change of
input parameters. The conclusions obtained from this study are pre
sented in the following.
The correlation analysis, AHC, SSE, and multiple regression analysis
techniques, as recommended and implemented in the current study,
have presented promising results by dimension reduction (i.e. elimi
nating the redundant input parameters) and choosing the statistically
significant parameters that affect the rockburst parameters (i.e. σRB and
IRB ). This procedure simplifies the rockburst assessment at the field
scale. The obtained dendrogram by AHC analysis (Fig. 5) showed that
there is no natural group in the compiled database except for two data
samples that were identified as outliers and subsequently were elimi
nated from the original database. Therefore, the database was identified
as a homogeneous database for further analysis. The statistically ascer
tained and selected dominant parameters affecting σ RB and IRB in deep
geotechnical conditions have been found as UCS, E, and σv for the
maximum rockburst stress (σ RB ), and ρ, K, E, υ and σ RB for the rockburst
risk index (IRB ) by SSE and multicollinearity analyses.
The proposed non-linear GEP-based and CART-based models by
providing the mathematical functions and visual patterns could unravel
the latent relationships between the rockburst parameters (i.e. σ RB and
IRB ) and their corresponding influential parameters, and the validity of
these models was proved based on several statistical indices (Table 7).
However, the GEP-based models with the values of 0.9398, 14.2493,
and 9.8025 for the performance indices of R2 , RMSE, and MAE for σ RB
and the values of 0.9459, 0.1947, and 0.1365 for the foregoing perfor
mance indices for IRB outperformed the CART models. The results show
that the used robust techniques can be useful tools for solving the high-
complex non-linear problems which are common in mining and
geotechnical projects. The results obtained from the parametric analysis
on the proposed GEP models for σ RB and IRB showed that the stress level
that rockburst occurs increases monotonically with the increase of UCS
and σ v . As such, the risk of rockburst occurrence showed a downward
non-linear trend with K, E, υ and σ RB parameters. On the other hand, the
parametric analysis revealed that there are strong correlations between
the rockburst parameters and their input parameters which show that
Fig. 13. Parametric analysis of σRB on GEP model. the selected inputs are potential indicators for assessing and predicting
rockburst phenomenon in deep underground openings.
According to the above results, clearly can be seen that there is a good
correlation between the rockburst parameters and the inputs. These Declaration of competing interest
equations provide a series of simple equations for calculating σ RB and IRB
based on the single rock mechanical parameter as a primary assessment. The authors declare that they have no known competing financial
These equations may be relevant to investigate rockburst potential. In the interests or personal relationships that could have appeared to influence
end, it is necessary to mention that the developed models in this study are the work reported in this paper.
16
Fig. 14. Parametric analysis of IRB on GEP model.
Table A1 Table A1 (continued )

If-then rules of the constructed regression tree for σRB . Node Rule
18 If UCS in [40.45, 151.3) and σv in [63.3, 169.455) then σRB ¼ 123.561 in

Node Rule 26.7% of cases
1 – 19 If UCS in [151.3, 234.1) and σv in [63.3, 169.455) then σRB ¼ 150.500 in
2 If UCS in [6.2, 40.45) then σRB ¼ 19.178 in 47.7% of cases 5.8% of cases
3 If UCS in [40.45, 234.1) then σRB ¼ 124.173 in 52.3% of cases 20 If σv in [63.3, 143.5) and UCS in [40.45, 151.3) then σRB ¼ 120.990 in 24.4%
4 If σv in [7.2, 37) and UCS in [6.2, 40.45) then σRB ¼ 18.658 in 46.5% of cases of cases
5 If σv in [37, 40) and UCS in [6.2, 40.45) then σRB ¼ 40 in 1.2% of cases 21 If σv in [143.5, 169.455) and UCS in [40.45, 151.3) then σRB ¼ 150.550 in
6 If σv in [22.3, 169.455) and UCS in [40.45, 234.1) then σRB ¼ 115.993 in 2.3% of cases
47.7% of cases 22 If σv in [63.3, 118.4) and UCS in [40.45, 151.3) then σRB ¼ 116.100 in 14.0%
7 If σv in [169.455, 282) and UCS in [40.45, 234.1) then σRB ¼ 208.025 in 4.7% of cases
of cases 23 If σv in [118.4, 143.5) and UCS in [40.45, 151.3) then σRB ¼ 127.511 in
8 If σv in [22.3, 63.3) and UCS in [40.45, 234.1) then σRB ¼ 89.331 in 15.1% of 10.5% of cases
cases 24 If σv in [63.3, 125.75) and UCS in [151.3, 234.1) then σRB ¼ 150.467 in 3.5%
9 If σv in [63.3, 169.455) and UCS in [40.45, 234.1) then σRB ¼ 128.371 in of cases
32.6% of cases 25 If σv in [125.75, 169.455) and UCS in [151.3, 234.1) then σRB ¼ 150.550 in
10 If UCS in [40.45, 91) and σv in [22.3, 63.3) then σRB ¼ 63.863 in 9.3% of 2.3% of cases
cases 26 If E in [14.1, 20.7) and σv in [169.455, 282) and UCS in [40.45, 234.1) then
11 If UCS in [91, 234.1) and σv in [22.3, 63.3) then σRB ¼ 130.080 in 5.8% of σRB ¼ 177 in 1.2% of cases
cases 27 If E in [20.7, 66.7) and σv in [169.455, 282) and UCS in [40.45, 234.1) then
12 If UCS in [40.45, 77.75) and σv in [22.3, 63.3) then σRB ¼ 70.700 in 5.8% of σRB ¼ 218.367 in 3.5% of cases
cases
13 If UCS in [77.75, 91) and σv in [22.3, 63.3) then σRB ¼ 52.467 in 3.5% of Table A2
cases
If-then rules of the constructed regression tree for IRB .
14 If σv in [22.3, 44.05) and UCS in [40.45, 77.75) then σRB ¼ 61.075 in 4.7% of
cases
15 If σv in [44.05, 63.3) and UCS in [40.45, 77.75) then σRB ¼ 109.200 in 1.2% Node Rule
of cases
1 –
16 If σv in [22.3, 34) and UCS in [91, 234.1) then σRB ¼ 145.533 in 3.5% of cases 2 If σRB in [10.6, 34.45) then IRB ¼ 1.729 in 44.2% of cases
17 If σv in [34, 63.3) and UCS in [91, 234.1) then σRB ¼ 106.900 in 2.3% of cases 3 If σRB in [34.45, 255.5) then IRB ¼ 0.435 in 55.8% of cases
(continued on next column) (continued on next page)
17
Table A2 (continued ) 8 He J, Dou L, Gong S, Li J, Ma Z. Rock burst assessment and prediction by dynamic
and static stress analysis based on micro-seismic monitoring. Int J Rock Mech Min Sci.
Node Rule
2017;93:46–53. https://doi.org/10.1016/J.IJRMMS.2017.01.005.
4 If K in [0.222, 0.62) and σRB in [10.6, 34.45) then IRB ¼ 1.982 in 27.9% of 9 Weng L, Huang L, Taheri A, Li X. Rockburst characteristics and numerical simulation
cases based on a strain energy density index: a case study of a roadway in Linglong gold
5 If K in [0.62, 2.245) and σRB in [10.6, 34.45) then IRB ¼ 1.295 in 16.3% of mine, China. Tunn Undergr Space Technol. 2017;69:223–232. https://doi.org/
10.1016/j.tust.2017.05.011.
cases
10 Sousa LR. Risk Analysis for Tunneling Projects. 2010.
6 If ν in [0.07, 0.28) and K in [0.222, 0.62) and σRB in [10.6, 34.45) then IRB ¼
11 Shi XZ, Zhou J, Dong L, Hu HY, Wang HY, Chen SR. Application of unascertained
2.651 in 8.1% of cases measurement model to prediction of classification of rockburst intensity. Chin J Rock
7 If ν in [0.28, 0.37) and K in [0.222, 0.62) and σRB in [10.6, 34.45) then IRB ¼ Mech Eng. 2010;29(1):2720–2726.
1.706 in 19.8% of cases 12 Tang C, Wang J, Zhang J. Preliminary engineering application of microseismic
8 If ν in [0.07, 0.165) and K in [0.222, 0.62) and σRB in [10.6, 34.45) then IRB monitoring technique to rockburst prediction in tunneling of Jinping II project.
¼ 3.009 in 1.2% of cases J Rock Mech Geotech Eng. 2010;2(3):193–208. https://doi.org/10.3724/SP.
9 If ν in [0.165, 0.28) and K in [0.222, 0.62) and σRB in [10.6, 34.45) then IRB J.1235.2010.00193.
¼ 2.591 in 7.0% of cases 13 Shirani Faradonbeh R, Taheri A. Long-term prediction of rockburst hazard in deep
10 If ν in [0.28, 0.335) and K in [0.222, 0.62) and σRB in [10.6, 34.45) then IRB underground openings using three robust data mining techniques. Eng Comput. June
2018:1–17. https://doi.org/10.1007/s00366-018-0624-4.
¼ 1.937 in 8.1% of cases
14 He M, Gong W, Wang J, et al. Development of a novel energy-absorbing bolt with
11 If ν in [0.335, 0.37) and K in [0.222, 0.62) and σRB in [10.6, 34.45) then IRB
extraordinarily large elongation and constant resistance. Int J Rock Mech Min Sci.
¼ 1.545 in 11.6% of cases 2014;67:29–42. https://doi.org/10.1016/j.ijrmms.2014.01.007.
12 If K in [0.62, 1.56) and σRB in [10.6, 34.45) then IRB ¼ 1.351 in 15.1% of 15 Saharan MR, Mitri H. Destress blasting as a mines safety tool: some fundamental
cases challenges for successful applications. Elsevier Procedia Eng. 2011;26:37–47. https://
13 If K in [1.56, 2.245) and σRB in [10.6, 34.45) then IRB ¼ 0.556 in 1.2% of doi.org/10.1016/j.proeng.2011.11.2137.
cases 16 He M, e Sousa LR, Miranda T, Zhu G. Rockburst laboratory tests database -
14 If ν in [0.27, 0.28) and K in [0.62, 1.56) and σRB in [10.6, 34.45) then IRB ¼ application of data mining techniques. Eng Geol. 2015;185:116–130. https://doi.org/
1.194 in 5.8% of cases 10.1016/j.enggeo.2014.12.008.
15 If ν in [0.28, 0.33) and K in [0.62, 1.56) and σRB in [10.6, 34.45) then IRB ¼ 17 Wang J-A, Park HD. Comprehensive prediction of rockburst based on analysis of
strain energy in rocks. Tunn Undergr Space Technol. 2001;16(1):49–57. https://doi.
1.450 in 9.3% of cases
org/10.1016/S0886-7798(01)00030-X.
16 If σRB in [34.45, 56.8) then IRB ¼ 0.892 in 8.1% of cases
18 Singh SP. The influence of rock properties on the occurrence and control of
17 If σRB in [56.8, 255.5) then IRB ¼ 0.357 in 47.7% of cases rockbursts. Min Sci Technol. 1987;5(1):11–18. https://doi.org/10.1016/S0167-9031
18 If E in [3.5, 13.9) and σRB in [34.45, 56.8) then IRB ¼ 0.445 in 2.3% of cases (87)90854-1.
19 If E in [13.9, 43.1) and σRB in [34.45, 56.8) then IRB ¼ 1.071 in 5.8% of cases 19 Cook NGW. The basic mechanics of rockbursts. J South Afr Inst Min Metall. 1963;64
20 If E in [13.9, 33.7) and σRB in [34.45, 56.8) then IRB ¼ 1.204 in 4.7% of cases (3):71–81.
21 If E in [33.7, 43.1) and σRB in [34.45, 56.8) then IRB ¼ 0.540 in 1.2% of cases 20 Linkov AM. Rockbursts and the instability of rock masses. Int J Rock Mech Min Sci
22 If K in [0, 1.721) and σRB in [56.8, 255.5) then IRB ¼ 0.488 in 19.8% of cases Geomech Abstr. 1996;33(7):727–732. https://doi.org/10.1016/0148-9062(96)
23 If K in [1.721, 5.866) and σRB in [56.8, 255.5) then IRB ¼ 0.265 in 27.9% of 00021-6.
21 Barquins M, Petit J-P. Kinetic instabilities during the propagation of a branch crack:
cases
effects of loading conditions and internal pressure. J Struct Geol. 1992;14(8-9):
24 If σRB in [56.8, 130.6) and K in [0, 1.721) then IRB ¼ 0.533 in 14.0% of cases
893–903. https://doi.org/10.1016/0191-8141(92)90021-N.
25 If σRB in [130.6, 255.5) and K in [0, 1.721) then IRB ¼ 0.380 in 5.8% of cases 22 Bagde MN, Petro�s V. Fatigue properties of intact sandstone samples subjected to
26 If E in [14.1, 58.6) and σRB in [56.8, 130.6) and K in [0, 1.721) then IRB ¼ dynamic uniaxial cyclical loading. Int J Rock Mech Min Sci. 2005;42(2):237–250.
0.512 in 12.8% of cases https://doi.org/10.1016/J.IJRMMS.2004.08.008.
27 If E in [58.6, 74.1) and σRB in [56.8, 130.6) and K in [0, 1.721) then IRB ¼ 23 Chen J, Su G. True triaxial experimental study on hard rock under high geo-stress
0.759 in 1.2% of cases condition based on excavation and support. In: 2010 International Conference on
28 If E in [14.1, 36.05) and K in [1.721, 5.866) and σRB in [56.8, 255.5) then IRB Mechanic Automation and Control Engineering, MACE2010. 2010:1465–1467. https://
¼ 0.331 in 15.1% of cases doi.org/10.1109/MACE.2010.5536121.
29 If E in [36.05, 71) and K in [1.721, 5.866) and σRB in [56.8, 255.5) then IRB ¼ 24 Wang X, Huang R. Analysis of deformation and failure features characteristics of rock
under unloading conditions and their effects on rock burst. Mt Res. 1998;16(4):
0.185 in 12.8% of cases
281–285.
30 If E in [14.1, 29.7) and K in [1.721, 5.866) and σRB in [56.8, 255.5) then IRB ¼
25 Xu LS. Research on the experimental rock mechanics of rockburst under unloading
0.292 in 12.8% of cases condition. J Chongqing Jianzhu Univ. 2003;22(1):1–4.
31 If E in [29.7, 36.05) and K in [1.721, 5.866) and σRB in [56.8, 255.5) then IRB 26 Cheon DS, Keon S, Park C, Ryu C. An experimental study on the brittle failure under
¼ 0.547 in 2.3% of cases true triaxial conditions. Tunn Undergr Space Technol. 2006;21(3-4).
32 If σRB in [56.8, 143.6) and E in [36.05, 71) and K in [1.721, 5.866) then IRB ¼ 27 Su G, Chen Z, Ju JW, Jiang J. Influence of temperature on the strainburst
0.205 in 9.3% of cases characteristics of granite under true triaxial loading conditions. Eng Geol. 2017;222:
33 If σRB in [143.6, 255.5) and E in [36.05, 71) and K in [1.721, 5.866) then IRB 38–52. https://doi.org/10.1016/J.ENGGEO.2017.03.021.
¼ 0.133 in 3.5% of cases 28 Wang Y, He M, Liu D, Gao Y. Rockburst in sandstone containing elliptic holes with
varying axial ratios. Ann Mater Sci Eng. 2019;2019.
29 He MC, Miao JL, Feng JL. Rock burst process of limestone and its acoustic emission
characteristics under true-triaxial unloading conditions. Int J Rock Mech Min Sci.
2010;47(2):286–298. https://doi.org/10.1016/j.ijrmms.2009.09.003.
30 He M, Xia H, Jia X, Gong W, Zhao F, Liang K. Studies on classification, criteria and
control of rockbursts. J Rock Mech Geotech Eng. 2012;4(2):97–114. https://doi.org/
10.3724/SP.J.1235.2012.00097.
References 31 Chen G, Li T, Wang W, Zhu Z, Chen Z, Tang O. Weakening effects of the presence of
water on the brittleness of hard sandstone. Bull Eng Geol Environ. 2019;78(3):
1471–1483. https://doi.org/10.1007/s10064-017-1184-3.
1 Zhao J, Li XB, Shi XZ. Long-term prediction model of rockburst in underground
32 Liu X, Liang Z, Zhang Y, Liang P, Tian B. Experimental study on the monitoring of
openings using heuristic algorithms and support vector machines. Saf Sci. 2012;50
rockburst in tunnels under dry and saturated conditions using AE and infrared
(4):629–644. https://doi.org/10.1016/j.ssci.2011.08.065.
monitoring. Tunn Undergr Space Technol. 2018;82:517–528. https://doi.org/
2 Ranjith PG, Zhao J, Ju M, De Silva RVS, Rathnaweera TD, Bandara AKMS.
10.1016/J.TUST.2018.08.011.
Opportunities and challenges in deep mining: a brief review. Engineering. 2017;3(4):
33 Sun X, Xu H, Zheng L, He M, Gong W. An experimental investigation on acoustic
546–551. https://doi.org/10.1016/J.ENG.2017.04.024.
emission characteristics of sandstone rockburst with different moisture contents. Sci
3 Feng XT, Chen B, Feng G, Zhao Z, Zheng H. Description and engineering
China Technol Sci. 2016;59(10):1549–1558. https://doi.org/10.1007/s11431-016-
phenomenon of rockbursts. In: Rockburst : Mechanisms, Monitoring, Warning, and
0181-8.
Mitigation. first ed. Elsevier - Health Sciences Division; 2017:3–17.
34 Akdag S, Karakus M, Taheri A, Nguyen G, Manchao H. Effects of thermal damage on
4 Cai M, Kaiser P. Rockburst Support Reference Book—Volume I: Rockburst Phenomenon
strain burst mechanism for brittle rocks under true-triaxial loading conditions. Rock
and Support Characteristics. Laurentian University; 2018.
Mech Rock Eng. June 2018:1–26.
5 Zhou J, Li X, Mitri HS. Evaluation method of rockburst: state-of-the-art literature
35 Su G, Jiang J, Zhai S, Zhang G. Influence of tunnel Axis stress on strainburst: an
review. Tunn Undergr Space Technol. 2018;81:632–659. https://doi.org/10.1016/j.
experimental study. Rock Mech Rock Eng. 2017;50(6):1551–1567. https://doi.org/
tust.2018.08.029.
10.1007/s00603-017-1181-7.
6 He M, Ren F, Liu D. Rockburst mechanism research and its control. Int J Min Sci
36 Zhao XG, Cai M. Influence of specimen height-to-width ratio on the strainburst
Technol. 2018;28(5):829–837. https://doi.org/10.1016/j.ijmst.2018.09.002.
characteristics of Tianhu granite under true-triaxial unloading conditions. Can
7 Li T, Cai MF, Cai M. A review of mining-induced seismicity in China. Int J Rock Mech
Geotech J. 2015;52(7):890–902. https://doi.org/10.1139/cgj-2014-0355.
Min Sci. 2007;44(8):1149–1171. https://doi.org/10.1016/j.ijrmms.2007.06.002.
18
37 Zhao XG, Wang J, Cai M, et al. Influence of unloading rate on the strainburst 61 Motulsky HJ, Ransnas LA. Fitting curves to data using nonlinear regression: a
characteristics of beishan granite under true-triaxial unloading conditions. Rock practical and nonmathematical review. Faseb J. 1987;1(5):365–374.
Mech Rock Eng. 2014;47(2):467–483. https://doi.org/10.1007/s00603-013-0443-2. 62 Alavi AH, Hasni H, Lajnef N, Chatti K, Faridazar F. Damage detection using self-
38 Li D, Zhao F, Zheng M. Fractal characteristics of cracks and fragments generated in powered wireless sensor data: an evolutionary approach. Measurement. 2016;82:
unloading rockburst tests. Int J Min Sci Technol. 2014;24(6):819–823. https://doi. 254–283. https://doi.org/10.1016/J.MEASUREMENT.2015.12.020.
org/10.1016/J.IJMST.2014.10.014. 63 Mitchell TM. Does machine learning really work? AI Mag. 1997;18(3):11–20.
39 Cortez P, Embrechts MJ. Opening black box data mining models using sensitivity 64 Alavi AH, Gandomi AH. A robust data mining approach for formulation of
analysis. In: 2011 {IEEE} Symposium on Computational Intelligence and Data Mining geotechnical engineering systems. Eng Comput. 2011;28(3):242–274. https://doi.
({CIDM}). IEEE; 2011. https://doi.org/10.1109/cidm.2011.5949423. org/10.1108/02644401111118132.
40 Hasanipanah M, Faradonbeh RS, Armaghani DJ, Amnieh HB, Khandelwal M. 65 Nazari A, Pacheco Torgal F. Modeling the compressive strength of geopolymeric
Development of a precise model for prediction of blast-induced flyrock using binders by gene expression programming-GEP. Expert Syst Appl. 2013;40(14):
regression tree technique. Environ Earth Sci. 2017;76(1). https://doi.org/10.1007/ 5427–5438. https://doi.org/10.1016/J.ESWA.2013.04.014.
s12665-016-6335-5. 66 Ferreira C. Gene expression programming in problem solving. In: Soft Computing and
41 Armaghani DJ, Faradonbeh RS, Rezaei H, Rashid ASA, Amnieh HB. Settlement Industry. London: Springer London; 2002:635–653. https://doi.org/10.1007/978-1-
prediction of the rock-socketed piles through a new technique based on gene 4471-0123-9_54.
expression programming. Neural Comput Appl. 2016. https://doi.org/10.1007/ 67 Ferreira C. Gene expression programming in problem solving. In: Roy R, K€ oppen M,
s00521-016-2618-8. Ovaska S, Furuhashi T, Hoffmann F, eds. Soft Computing and Industry: Recent
42 Khandelwal M, Armaghani DJ, Faradonbeh RS, Yellishetty M, Majid MZA, Applications. London: Springer London; 2002:635–653. https://doi.org/10.1007/
Monjezi M. Classification and regression tree technique in estimating peak particle 978-1-4471-0123-9_54.
velocity caused by blasting. Eng Comput. 2017;33(1):45–53. https://doi.org/ 68 Power HE, Gharabaghi B, Bonakdari H, Robertson B, Atkinson AL, Baldock TE.
10.1007/s00366-016-0455-0. Prediction of wave runup on beaches using Gene-Expression Programming and
43 Salimi A, Faradonbeh RS, Monjezi M, Moormann C. TBM performance estimation empirical relationships. Coast Eng. 2019;144:47–61. https://doi.org/10.1016/J.
using a classification and regression tree (CART) technique. Bull Eng Geol Environ. COASTALENG.2018.10.006.
2016. https://doi.org/10.1007/s10064-016-0969-0. 69 Ferreira C. Gene Expression Programming: Mathematical Modeling by an Artificial
44 He M. The mechanism of rockburst and its countermeasure of support. Consult Rep Intelligence. Springer; 2006.
Key Technol Safe Rapid Constr Jinping LI Hydropower Stn High Overburden Long Tunnels. 70 Roy R, K€ oppen M, Ovaska S, Furuhashi T, Hoffmann F. Soft Computing and Industry:
2009:23–28. Recent Applications. 2002. https://doi.org/10.1007/978-1-4471-0123-9.
45 Hudaverdi T. Application of multivariate analysis for prediction of blast-induced 71 Kayadelen C. Soil liquefaction modeling by genetic expression programming and
ground vibrations. Soil Dynam Earthq Eng. 2012;43:300–308. https://doi.org/ neuro-fuzzy. Expert Syst Appl. 2011;38(4):4080–4087. https://doi.org/10.1016/j.
10.1016/j.soildyn.2012.08.002. eswa.2010.09.071.
46 Faradonbeh RS, Monjezi M. Prediction and minimization of blast-induced ground 72 Hoseinian FS, Faradonbeh RS, Abdollahzadeh A, Rezai B, Soltani-Mohammadi S.
vibration using two robust meta-heuristic algorithms. Eng Comput. 2017:1–17. Semi-autogenous mill power model development using gene expression
47 Kaufman L, Rousseeuw PJ. Finding Groups in Data: An Introduction to Cluster Analysis. programming. Powder Technol. 2017;308. https://doi.org/10.1016/j.
vol 344. John Wiley & Sons; 2009. powtec.2016.11.045.
48 Saxena A, Prasad M, Gupta A, et al. A review of clustering techniques and 73 Kantardzic M. Data Mining: Concepts, Models, Methods, and Algorithms. John Wiley &
developments. Neurocomputing. 2017;267:664–681. https://doi.org/10.1016/J. Sons; 2003.
NEUCOM.2017.06.053. 74 Hasanipanah M, Faradonbeh RS, Amnieh HB, Armaghani DJ, Monjezi M. Forecasting
49 Sayadi AR, Lashgari A, Paraszczak JJ. Hard-rock LHD cost estimation using single blast-induced ground vibration developing a CART model. Eng Comput. 2017;33(2).
and multiple regressions based on principal component analysis. Tunn Undergr Space https://doi.org/10.1007/s00366-016-0475-9.
Technol. 2012;27(1):133–141. https://doi.org/10.1016/j.tust.2011.08.006. 75 Mahjoobi J, Etemad-Shahidi A. An alternative approach for the prediction of
50 James G, Witten D, Hastie T, Tibshirani R. An Introduction to Statistical Learning. vol significant wave heights based on classification and regression trees. Appl Ocean Res.
112. Springer; 2013. 2008;30(3):172–177. https://doi.org/10.1016/J.APOR.2008.11.001.
51 Montgomery DC, Peck EA, Vining GG. Introduction to Linear Regression Analysis. 76 Quinlan JR. Induction of decision trees. Mach Learn. 1986;1(1):81–106. https://doi.
Wiley; 2012. org/10.1007/BF00116251.
52 Kumar Sharma S, Rai P. Establishment of blasting design parameters influencing 77 Ghasemi E, Kalhori H, Bagherpour R. Stability assessment of hard rock pillars using
mean fragment size using state-of-art statistical tools and techniques. Measurement. two intelligent classification techniques: a comparative study. Tunn Undergr Space
2017;96:34–51. https://doi.org/10.1016/J.MEASUREMENT.2016.10.047. Technol. 2017;68:32–37. https://doi.org/10.1016/j.tust.2017.05.012.
53 Hemmateenejad B, Yazdani M. QSPR models for half-wave reduction potential of 78 Breiman L, Friedman J, Stone CJ, Olshen RA. Classification and Regression Trees. CRC
steroids: a comparative study between feature selection and feature extraction from press; 1984.
subsets of or entire set of descriptors. Anal Chim Acta. 2009;634(1):27–35. https:// 79 Bevilacqua M, Braglia M, Montanari R. The classification and regression tree
doi.org/10.1016/j.aca.2008.11.062. approach to pump failure rate analysis. Reliab Eng Syst Saf. 2003;79(1):59–67.
54 Sarkhosh M, Ghasemi JB, Ayati M. A quantitative structure- property relationship of https://doi.org/10.1016/S0951-8320(02)00180-1.
gas chromatographic/mass spectrometric retention data of 85 volatile organic 80 Smith NG. Probability and Statistics in Civil Engineering. vol 244. Collins Prof Tech
compounds as air pollutant materials by multivariate methods. Chem Cent J. 2012;6 Books; 1986. http://ci.nii.ac.jp/naid/10007808566/en/. Accessed May 29, 2019.
(suppl 2). https://doi.org/10.1186/1752-153X-6-S2-S4 (Suppl 2):S4-S4. 81 Golbraikh A, Tropsha A. Beware of q2!. J Mol Graph Model. 2002;20(4):269–276.
55 Shirani Faradonbeh R, Taheri A. Long-term prediction of rockburst hazard in deep https://doi.org/10.1016/S1093-3263(01)00123-1.
underground openings using three robust data mining techniques. Eng Comput. 2018. 82 Roy PP, Roy K. On some aspects of variable selection for partial least squares
https://doi.org/10.1007/s00366-018-0624-4. regression models. QSAR Comb Sci. 2008;27(3):302–313. https://doi.org/10.1002/
56 Pu Y, Apel DB, Xu H. Rockburst prediction in kimberlite with unsupervised learning qsar.200710043.
method and support vector classifier. Tunn Undergr Space Technol. 2019;90:12–18. 83 Mohammadzadeh SD, Bolouri Bazaz J, Vafaee Jani Yazd SH, Alavi AH. Deriving an
https://doi.org/10.1016/j.tust.2019.04.019. intelligent model for soil compression index utilizing multi-gene genetic
57 Archontoulis SV, Miguez FE. Nonlinear regression models and applications in programming. Environ Earth Sci. 2016;75(3):262. https://doi.org/10.1007/s12665-
agricultural research. Agron J. 2015;107(2):786–798. 015-4889-2.
58 Bethea RM. Statistical Methods for Engineers and Scientists. Routledge; 2018. 84 Soleimani S, Rajaei S, Jiao P, Sabz A, Soheilinia S. New prediction models for
59 Jahed Armaghani D, Faradonbeh RS, Momeni E, Fahimifar A, Tahir MM. unconfined compressive strength of geopolymer stabilized soil using multi-gen
Performance prediction of tunnel boring machine through developing a gene genetic programming. Measurement. 2018;113:99–107. https://doi.org/10.1016/j.
expression programming equation. Eng Comput. 2017. https://doi.org/10.1007/ measurement.2017.08.043.
s00366-017-0526-x. 85 Wagner H. Deep mining: a rock engineering challenge. Rock Mech Rock Eng. April
60 Ghasemi E. Particle swarm optimization approach for forecasting backbreak induced 2019. https://doi.org/10.1007/s00603-019-01799-4.
by bench blasting. Neural Comput Appl. 2017;28(7):1855–1862. https://doi.org/
10.1007/s00521-016-2182-2.
19

Jurnal 2

Uploaded by

Copyright:

Available Formats

Jurnal 2

Uploaded by

Document Information

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Jurnal 2

Uploaded by

Copyright:

Available Formats

International Journal of Rock Mechanics & Mining Sciences 128 (2020) 104279

Contents lists available at ScienceDirect

International Journal of Rock Mechanics and Mining Sciences

Rockburst assessment in deep geotechnical conditions using true-triaxial

1. Introduction defined as a dynamic instability phenomenon of the surrounding rock

density, rock specific weight, mineral contents of rocks, loading and

Fig. 4. Process of rockburst assessment in this

* Correlation is significant at the 0.05 level (2-tailed).

Dependent parameter r r-square Adjusted r-square Std. error of the estimate

σRB 0.952 0.906 0.897 19.385

Rockburst maximum stress (σRB )

(Constant) 71.247 20.947 3.401 0.001

Rockburst risk index (IRB )

(Constant) 2.766 0.282 9.816 0.000

Table 4 advance by the user.64 By inspiring from the Darwinian principle of

Model Unstandardized coefficients Standardized coefficients t Sig. Collinearity statistics

B Std. error Beta Tolerance VIF

(Constant) 6.019 3.080 1.954 0.053

Rockburst risk index (IRB )

(Constant) 2.780 0.278 10.006 0.000

Fig. 7. Process of function finding using GEP algorithm.

algorithm.46,71,72 The other GEP parameters including the number of

Terminal set UCS; E; σv σRB ; ​ K; ​ υ; ​ E; ​ ρ

increases, and finally, bursting occurs at a higher stress level violently.18

σ RB ¼ 19:911LnðUCSÞ þ 10:636; R2 ¼ 0:9974 (14)

σ RB ¼ 0:0007σ2v þ 0:7162σv þ 41:092; R2 ¼ 0:9877 (15)

IRB ¼ 0:49e 0:535υ

IRB ¼ 0:0186K 3 þ 0:2124K 2 0:7981K þ 1:1997; R2 ¼ 0:9521 (17)

GEP CART GEP CART

hi : measured output; ti : predicted output.

6. Summary and conclusions

As a catastrophic geohazard, rockburst threatens the safety of

Fig. 14. Parametric analysis of IRB on GEP model.

Table A1 Table A1 (continued )

18 If UCS in [40.45, 151.3) and σv in [63.3, 169.455) then σRB ¼ 123.561 in

You might also like

Terminal set UCS; E; σv σRB ; K; υ; E; ρ