Probabilistic evaluation of ground-support interaction for deep rock excavation using artificial neural network and uniform design

Tunnelling and Underground Space Technology 32 (2012) 1–18

Contents lists available at SciVerse ScienceDirect

Tunnelling and Underground Space Technology journal homepage: www.elsevier.com/locate/tust

Probabilistic evaluation of ground-support interaction for deep rock excavation using artificial neural network and uniform design Qing Lü a,b,⇑, Chin Loong Chan b, Bak Kong Low b a b

College of Civil Engineering and Architecture, Zhejiang University, Hangzhou 310058, PR China School of Civil and Environmental Engineering, Nanyang Technological University, Singapore 639798, Singapore

a r t i c l e

i n f o

Article history: Received 16 October 2011 Received in revised form 3 April 2012 Accepted 11 April 2012 Available online 12 June 2012 Keywords: Rock tunnel Convergence–confinement method Artificial neural network Uniform design FORM SORM

a b s t r a c t An efficient approach is proposed in this paper for probabilistic ground-support interaction analysis of deep rock excavation using the artificial neural network (ANN) and uniform design. The deterministic model is based on the convergence–confinement method. The ANN model is employed as the response surface to fit the real limit state surface. The uniform design table is used to prepare the sampling points for training the ANN and for determining the parameters of the network via an iterative procedure. The probability of failure is estimated from the first-order and second-order reliability method (FORM/SORM) based on the generated ANN response surface and compared with Monte Carlo simulations and polynomial response surface method. The efficiency and the accuracy of the proposed approach are first illustrated with the case of a circular tunnel involving analytical solutions with respect to three performance functions. The results show that the support installation position and the parametric correlations have great influence on the probability of the three failure modes. Reliability analyses involving four-parameter beta distributions are also investigated. Finally, an example of a deep rock cavern excavation is presented to illustrate the feasibility of the proposed approach for practical applications where complex numerical procedures are needed to compute the performance function. Ó 2012 Elsevier Ltd. All rights reserved.

1. Introduction The ground-support interaction analysis of underground rock excavation is usually carried out using a deterministic approach based on the mean or average properties of the rock mass strength and deformation characteristics in tunnel design (e.g., CarranzaTorres and Fairhurst, 2000; Oreste, 2003; Vardakos et al., 2007; Gonzalez-Nicieza et al., 2008). However, due to the inherent uncertainties of the rock mass in nature, a reliability-based analysis and design approach is more rational since it can account for the parametric uncertainties and correlations explicitly, and provides a means of distinguishing the important from unimportant uncertainties of the input parameters. Theories and methods of reliability analysis have been wellestablished in the field of structural reliability (Rackwitz, 2001) but their application to underground excavation is still uncommon. Only a few publications can be found in some recent studies on this topic. For example, Hoek (1998) presented a reliability analysis of a circular tunnel using Monte Carlo simulation (MCS). Li and Low (2010) adopted the algorithm of Low and Tang (2007) for FORM to compute the reliability index of the tunnel in Hoek (1998). ⇑ Corresponding author at: College of Civil Engineering and Architecture, Zhejiang University, Hangzhou 310058, PR China. Tel./fax: +86 57188208732. E-mail address: [email protected] (Q. Lü). 0886-7798/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.tust.2012.04.014

The influences of the support pressure and the in situ stress on the reliability index were investigated. In the above two studies, the performance functions were evaluated from the Duncan Fama (1993) analytical solutions of a circular tunnel based on Mohr–Coulomb yield criterion, thus MCS and FORM can be directly applied on the closed-form performance functions to compute the probability of failure. For many practical applications, the performance function may not be available in explicit form. In these cases, response surface method (RSM) is commonly used to approximate the original implicit performance function with a simple polynomial function. Mollon et al. (2009) presented a reliability analysis of circular tunnels in homogeneous soil based on the RSM. Two numerical models developed in FLAC3D were used to construct the response surface. Recently the present writers published two applications of RSM for reliability analysis of underground rock excavation and ground-support interaction (Lü and Low, 2011; Lü et al., 2011), in which the quadratic polynomial with cross terms combined with an iterative procedure was used to obtain a fitted response surface for the real limit state surface (LSS). FORM/SORM were then performed on the response surface. However, it was found in Mollon et al. (2009) that the quadratic approximation for the LSS cannot converge for the case concerning serviceability limit state; instead, a specific function was suggested as the response surface. The non-convergence issue was also encountered

2

Q. Lü et al. / Tunnelling and Underground Space Technology 32 (2012) 1–18

when a quadratic polynomial without cross terms was employed for certain cases presented in this study, as will be seen in Section 5.1. The inclusion of the cross terms speeds up the iteration process, enables quicker convergence, and improves the accuracy of the approximation, especially for SORM (Lü et al., 2011). However, when all cross terms are included, the total number of the sampling points required to construct the response surface in each iteration will increase from 2n + 1 to (n + 1)(n + 2)/2 (n denotes the number of random variables). The effort required to evaluate the sampling points can be tremendous, when a large number of random variables are involved. In this paper, the artificial neural network (ANN) is used as an alternative form of response surface (referred to as ANN-RS) to approximate the real LSS for probabilistic ground-support interaction analysis of deep rock excavation. The uniform design (UD) is employed to prepare the sampling points for training the ANN and obtaining the unknown parameters of the network. An iterative algorithm that operates on the independent standard normal basic random variables in U-space is presented. Three failure modes with respect to the support capacity, the permissible tunnel convergence and the plastic zone extent are investigated. The sensitivities of the computed failure probability to the support installation position, the parametric correlations and distributions are studied. The purpose of this paper is to illustrate a practical and transparent procedure in a ubiquitous platform to conduct probabilistic analysis of the ground-support system using the ANN and UD approach. If desired, other platforms such as MATLAB can be used. The authors hope that the transparency of operation and clarity of concept of the proposed procedure would reduce the computational and conceptual barriers that limit the application of reliability approaches in rock engineering.

2. Ground-support interaction analyses of a circular tunnel using CCM Several approaches have been proposed for determining the support required to stabilize an underground excavation and predicting the final tunnel wall convergence in tunnel design, among which the convergence–confinement method (CCM) is regarded as one of the most realistic and appropriate way to account for the actual mechanical interaction of ground and support (Carranza-Torres and Diederichs, 2009). The principle of CCM is well-documented in the literature (e.g., Brown et al., 1983; Carranza-Torres and Fairhurst, 2000; Oreste, 2003). A brief description is given below. The CCM for tunnels comprises three main procedures, and three corresponding basic curves, namely the Longitudinal Deformation Profile (LDP), the Ground Reaction Curve (GRC) and the Support Characteristic Curve (SCC), are required, as illustrated in Fig. 1. For a circular tunnel under hydrostatic in situ stress, these curves can be constructed using some established analytical solutions. In this paper the formulas provided by Vlachopoulos and Diederichs (2009) is employed to approximate the LDP, the Carranza-Torres solution (2004) is adopted to construct the GRC and the equations for concrete lining proposed in Oreste (2003) are used to calculate the SCC parameters for the circular tunnel example in Section 3. To illustrate the procedure of ground-support interaction analysis using LDP, GRC and SCC, assume that the support is installed at a distance Ls from the tunnel face and the ground has already converged by an amount given by uin r (point B of LDP in Fig. 1). This uin r is also the starting point (point D) of SCC since the vertical coordinate of LDP and the horizontal coordinate of GRC (SCC) are the same as shown in Fig. 1. As the tunnel face

Fig. 1. Schematic representation of the Longitudinal Deformation Profile (LDP), Ground Reaction Curve (GRC) and Support Characteristic Curve (SCC) and their interaction (after Carranza-Torres and Fairhurst (2000)).

3

Q. Lü et al. / Tunnelling and Underground Space Technology 32 (2012) 1–18

Fig. 2. A circular tunnel under a hydrostatic in situ stress and supported with concrete lining.

advances, the ground and the support deform together at the section under investigation. The pressure on the support increases (from point D of SCC in Fig. 1) and the simulated support pressure provided by the tunnel face decreases (from point F of GRC in Fig. 1). Once the tunnel face has advanced sufficiently far ahead, its supporting effect on the section of concern becomes negligible, and the tunnel support system reaches equilibrium at point K (i.e., the intersection of GRC and SCC). At this point, the pressure pDs defined by point K then represents the final design pressure that the ground transmits to the support and uDr represents the final deformation of the support system. Thus, determining the tunnel wall convergence uin r at the time of support installation and seeking the intersection point of GRC and SCC are the two key operations of the aforementioned procedure, which are automatically executed using VBA program codes in Microsoft Excel created in this study for the circular tunnel. 3. Performance functions and parameters of a circular tunnel According to the above principle of the CCM, the performance functions (limit state functions) of the circular tunnel support system may be given as:

g 1 ðxÞ ¼ pmax  pDs s

ð1Þ

g 2 ðxÞ ¼ emax  u

uDr Rt

ð2Þ

 g 3 ðxÞ ¼ nmax pl

Rpl Rt

ð3Þ

where x is the vector of random variables, pmax is the support capacs ity, pDs is the equilibrium pressure acting on the support (Fig. 1), Rt is the tunnel radius, uDr is the total inward radial displacement of the tunnel (Fig. 1), Rpl is the radius of the plastic zone, emax and nmax are u pl the maximum permissible ratios with respect to uDr =Rt and Rpl/Rt, respectively. In Eqs. (1)–(3) the values of pDs ; uDr , and Rpl are evaluated based on x via the CCM procedure and therefore are random variables, whereas pmax ; emax ; nmax and Rt are design parameters that s u pl are considered as constants. The above performance functions define three different failure modes of unsatisfactory performance of the tunnel support system. The first performance function is the support capacity criterion and becomes negative (i.e., support collapses or yields) when the pmax is s smaller than pDs . The second performance function is the tunnel

convergence criterion and becomes negative when the inward displacement of the wall at the design point uDr is such that the uDr =Rt ratio is greater than the limiting ratio emax (emax =2% in this study). u u The third performance function is the plastic zone criterion and becomes negative when the extent of plastic zone is too large and exceeds the limiting ratio nmax (nmax = 2 in this study). The first pl pl performance function involves the ultimate limit state of support capacity while the second and the third are akin to serviceability limit states. These three performance criteria are investigated separately below. In the section to follow, a circular tunnel of radius 4 m, that is assumed to be constructed in limestone at a depth of 1000 m subjected to a hydrostatic in situ stress of 27 MPa, is chosen as the first illustrative example (Fig. 2). The rock mass is regarded as an equivalent homogeneous and isotropic material with few discontinuities or with dense joints, thus only stress-controlled instability mechanisms are investigated. The structrurally-controlled instability mechanisms which are formed by intersecting discontinuities are not considered in this study. The Mohr–Coulomb (M–C) yield criterion is one of the most frequently used models for underground excavation analyses, although it cannot reflect the non-linear nature of the rock mass failure envelope. To overcome this shortcoming of the linear M–C criterion, the non-linear Hoek–Brown (H–B) yield criterion, proposed by Hoek and Brown (1980) and refined by Hoek et al. (2002), is widely used to describe the strength of equivalent rock mass in rock engineering. A set of equations for estimating equivalent M–C parameters from H–B inputs are also provided by Hoek (1998). However, this equivalent transformation from H–B paramerters to M–C parameters was found to generate large discrepancy in ultimate bearing capacity (Merifield et al., 2006), slope stability factors of safety (Li et al., 2008) and failure probability of underground rock excavation (Lü and Low, 2011). Hence it is better to use the H–B criterion directly. In this paper the H–B Table 1 Statistical properties of random variables for the illustrative example. Variables

Distributions

l

r

COV

GSI mi rci (MPa) Erm (GPa) p0 (MPa) ks (MPa/m) pmax (MPa) s

Normal Normal LogNormal LogNormal Normal LogNormal LogNormal

35 10 120 11.2 27 1073.7 4.684

3.5 1.5 30 4.3 4.1 245.9 1.037

0.1 0.15 0.25 0.38 0.15 0.23 0.22

4

Q. Lü et al. / Tunnelling and Underground Space Technology 32 (2012) 1–18

Fig. 3. FORM with respect to the support capacity criterion g1(x) (i.e., Eq. (1)) with Ls = 2 m.

criterion is used to compute the stress and deformation distributions of the rock mass surrounding the excavation. For the circular tunnel in hand, a 50 cm thick concrete lining is to be installed continuously 2 m away from the advancing face to stabilize the highly stressed rock mass. Seven parameters are regarded as basic random variables, including the rock mass properties, the in situ stress and the support characteristic parameters. The statistical properties and distributions of these random variables are listed in Table 1. The distributions and the coefficients of variation (COV) of rock mass parameters GSI, mi and rci in Table 1 are based on assumptions and some published data in Kim and Gao (1995), Hoek (1998), and Sari (2009). The statistical properties of Erm are computed from MCS of 105 trials using the following Eq. (4) proposed by Hoek and Diederichs (2006), with the input values of GSI, mi, rci as shown in Table 1, and with a lognormally distributed modulus ratio (MR) of l = 800, COV = 0.15.

 Erm ¼ Ei 0:02 þ

1  D=2 1 þ expðð60 þ 15D  GSIÞ=11Þ

 ð4Þ

where Ei is the Young’s modulus of the intact rock, which can be estimated from MR by Ei ¼ MR  rci , and D is the disturbance factor which is assumed to be 0 in this study. The above MCS result also indicates that Erm is correlated with GSI and rci. Both correlation coefficients are 0.6. The statistical characteristics and distributions of SCC parameters (i.e., pmax and s ks) are also fitted from MCS of 105 trials, using the formulas in Oreste (2003) by assuming that COV of tc equals 0.1, COV of rcc and Ec equals 0.2, and all three variables follow lognormal distributions. Positive correlation between ks and pmax is found in the MCS result s with a correlation coefficient of 0.2. The above generated correlation coefficients are reflected in the correlation matrix R of Fig. 3. The influence of the correlation coefficients on the reliability results will be investigated later in Section 5.4. Note that although this paper focuses on methodology and computation, realistic statistical inputs is of equal importance, because the computed failure probability depends on the inputs. 4. Reliability analysis framework using ANN and uniform design 4.1. Reliability index and FORM In classical FORM, the original correlated basic random vector x in X-space are transformed into independent standard normal vector u in U-space. The reliability index b is defined as:

pffiffiffiffiffiffiffiffiffi b ¼ min uT u

ð5Þ

u2F

where F denotes the failure domain in which the performance function gðuÞ 6 0. The point denoted by u which minimizes Eq. (5) and satisfies u 2 F is the design point (DP), and is also the point on the LSS (g(u) = 0) that is closest to the origin of the U-space. An iterative algorithm may be used to search for u, which has been well-documented in Ang and Tang (1984), Ditlvsen and Madsen (1996), Melchers (1999), Baecher and Christian (2003), Choi et al. (2006), for example. A new algorithm for FORM was proposed by Low and Tang (2007), which interpreted the reliability index based on the perspective of an expanding ellipsoid in the original X-space of the basic random variables, and the reliability index is expressed as:

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi b ¼ min nT R1 n

ð6Þ

x2F

where R is the correlation matrix and n is the vector of unrotated and correlated equivalent standard normal variables which is related to the rotated and uncorrelated equivalent standard normal vector u (Low et al., 2011):

n ¼ Lu

ð7Þ T

where L is a Cholesky lower triangular matrix that satisfies LL = R. In the refined algorithm for FORM (Low and Tang, 2007), the original random variables xi can be computed as a function of ni by the following expression:

xi ¼ F 1 ½Uðni Þ

ð8Þ

where UðÞ is the CDF of the standard normal variable, F 1 ðÞ is the inverse CDF of the original non-normal variable. Based on the reliability index b, the probability of failure Pf can be estimated by:

Pf  1  UðbÞ

ð9Þ

In this paper the algorithm of Low and Tang (2007) is used to compute the DP and the reliability index b for constructing the ANN model, as shown in Fig. 3. For the case in Fig. 3, the failure probability is about 0.104% from FORM. The solution in the x column represents the DP on the LSS. For each parameter, the ratio of the mean value to the DP value is similar in nature to the partial factors in the limit state design (e.g., Eurocode 7) (Low, 2008). However, in a reliabilitybased design one does not specify the partial factors. The DP values x are determined automatically and reflect sensitivities, standard

Q. Lü et al. / Tunnelling and Underground Space Technology 32 (2012) 1–18

5

Fig. 4. Outline of SORM procedure and summary of results with respect to support capacity criterion g1(x) (i.e., Eq. (1)) with Ls = 2 m. Procedural details are available in Lü and Low (2011).

deviations, correlation structure, and probability distributions in a way that prescribed partial factors cannot reflect. Despite the merits and deeper insights provided by a probabilistic approach, one should note that more statistical input information is required, and that the results are only as good as the input information. In addition, the non-zero values of n imply that the reliability index of the support capacity criterion is sensitive to all the seven random variables, with pmax being the highest on the sensitivity scale. s The mean value point, at (35, 10, 120, 11.2, 27, 1074, 4.684) is safe against support collapse; but support failure occurs when the random variable values changed to the values of (34.33, 9.529, 105.2, 8.482, 33.01, 1025, 2.603). The ability of the DP to reflect relative parametric sensitivities automatically is an important desirable feature of a reliability analysis as opposed to a deterministic analysis or a design based on a rigid partial factors. 4.2. SORM In general, the SORM is used to enhance the accuracy of the results obtained from FORM, by fitting a quadratic (second-order) surface to the actual LSS in the vicinity of the DP to capture the non-linearity of LSS. The key concept of SORM is to evaluate the multi-dimensional probability integral over this approximate surface. One approach for fitting such surface is to obtain the principal curvature of LSS at DP. This can be done through a procedure of coordinate rotation and matrix diagonalization, e.g., Choi et al. (2006), Lü and Low (2011), Lü et al. (2011). Another approach for SORM is to adopt a point fitting strategy (Der Kiureghian et al., 1987). A practical procedure that adopts this algorithm in the Microsoft Excel platform is presented in Chan and Low (in press). Once the principal curvatures have been obtained, the probability of failure can be evaluated using different established formulas. In this paper, five SORM formulas (Breitung, 1984; Tvedt, 1983; Koyluoglu and Nielsen, 1994; Cai and Elishakoff, 1994; Zhao and Ono, 1999) are employed to estimate the failure probability for mutual comparison and verification. Second-order reliability analyses are performed based on the explicit ANN function which provides an approximation of the real LSS, as explained in the next section and in the Appendix A. All computations of SORM are executed automatically in Microsoft Excel using the procedures proposed in Lü and Low (2011). The main outputs are summarized in Fig. 4. For the case in hand, the failure probability based on SORM ranges from 0.082% to 0.089%, with an average of 0.087%. For comparison, six MCSs with Latin hypercube sampling (each comprising 200,000 trials) are also performed using the software @RISK (www.palisade.com). The failure probabilities are 0.087%, 0.092%, 0.094%, 0.088%, 0.085% and 0.092%, respectively, averaging 0.089%. Thus the relative percentage error of the FORM, given by the formula (Pf,FORM  Pf,MCS)/Pf,MCS  100%, is about 16.9%. On the other hand, the average failure probability obtained from SORM

has a relative percentage error of only 2.2%. Note that FORM is a necessary step in SORM, which builds on the DP and the reliability index of FORM. The failure domain is convex for the case in Fig. 4 because the predominant principal curvature (i.e., j5) is positive. Thus the failure probability computed by FORM on the linear assumption is overestimated for such cases with convex failure domains, as in Fig. 4. On the contrary, a predominantly negative principal curvature will yield an underestimation of the failure probability by FORM. 4.3. ANN response surface Generally, ANN is a powerful paradigm for mapping the relationship between a set of inputs and one or more outputs by means of training data obtained from either real experiments or numerical simulations. A number of applications of ANN in civil engineering have been reported (e.g., Ghaboussi et al., 1991; Goh, 1994; Juang et al., 1999; Goh and Kulhawy, 2005; Cheng and Li, 2008; Cho, 2009). A typical three-layer network that is adopted as a response surface (referred to as ANN-RS) to approximate the real LSS for reliability analysis in this paper is shown in Fig. 5. The architecture of this ANN is composed of a number of neurons which are connected with weights and logically arranged into the input layer, the hidden layer and the output layer. The presence of the hidden neurons enables the network to represent and compute the complicated associations between patterns (Ghaboussi et al., 1991). It has been theoretically proven that such networks with a single hidden layer can approximate any function, provided sufficient

Fig. 5. Architecture of the three-layer network used in this paper, with n inputs, h hidden neurons and one output.

6

Q. Lü et al. / Tunnelling and Underground Space Technology 32 (2012) 1–18

hidden neurons are available (Hornik et al., 1989). The network is essentially trained by adapting their weights and biases using optimization methods to minimize the mean square error between the predicted and the target values. In order not to obscure the discussion of the practical application of the method, details of the ANN response surface is presented separately in Appendix A. For the case in Figs. 3 and 4, an ANN of 7-20-1 structure was used which has seven neurons (equal to the number of random variables) in the input layer, 20 neurons in the hidden layer and one neuron in the output layer. The ANN was trained using 55 training data prepared by the uniform design methodology which is discussed next.

4.4. Preparing training data set using uniform design For a reliability analysis based on ANN-RS, selecting the sampling points and determining the number of samples are very important for establishing good ANN-RS. In the published literature, the sampling points are usually selected using either a

Fig. 6. Comparison of sampling point distribution generated by random sampling method and uniform design in a 2D domain.

random sampling method (e.g., Goh and Kulhawy, 2005; Cho, 2009) or an experimental design method (e.g., Wong, 1985). However, the random sampling does not guarantee that the selected samples cover the domain of interest uniformly (see Fig. 6), especially when the number of random variables is large and the number of training data is relatively small (Cheng and Li, 2009). On the other hand, the main drawback of the conventional experimental design technique is the large number of sampling data required to build a sufficiently accurate response surface. For a large number of random variables, the required number of evaluations of the real performance function grows exponentially (Schueremans and Van Gemert, 2005), which may be prohibitively costly for some cases involving complex numerical procedures. Uniform design (UD), proposed by Fang (1980) and Wang and Fang (1981), provides a space filling design method which guarantees that the sampling points are scattered in the domain as uniformly as possible. For instance, in Fig. 6 the 25 points sampled with UD in the range of [1, 1] are uniformly distributed in 25 lattices with each lattice containing one sampling point. In contrast, the same number of points generated by the random sampling method are not that uniformly distributed (Fig. 6). A few lattices contain more than one random point; some others have none. Thus, to ensure that the random samples cover the whole domain adequately, more sampling points are usually required when using the random sampling technique. In general, the merits of UD include: (1) it reduces the number of experimental configuration, especially when the experimental region has many factors and multiple levels (Kuo et al., 2007); (2) it is able to produce samples with high representativeness in the domain (Zhang et al., 1998). Due to these merits, UD has been successfully used in selecting samples for ANN training in Zhang et al. (1998), Kuo et al. (2007), Cheng and Li (2008, 2009). However, in these publications there are still no precise guidelines for determining the amount of the sampling data. For example, Cheng and Li (2008, 2009) reported the application of UD for selecting the ANN training data set in structural reliability analysis and optimization, in which the total training data number is still arbitrarily determined. However, in principle the number of sampling points required to accurately fit the LSS is dependent on the number of random variables and the nonlinearity of the problem considered (Goh, 1995).

Fig. 7. Preparing the training data set for establishing ANN model using uniform design table.

Q. Lü et al. / Tunnelling and Underground Space Technology 32 (2012) 1–18

For convenient use of UD, some uniform design tables are available on the internet (http://www.math.hkbu.edu.hk/UniformDesign) and in Fang (1994), which are also used in this paper for preparing the sampling points uj and computing g(uj) (Fig. 7). In Fig. 7, a UD table denoted as U15(157) is presented as an illustrative example, in which the superscript ‘7’ represents the number of random variables (factors) and ‘15’ denotes the number of experiments (levels) for each random variate. The sampling range (-3.0 to 3.0 in Fig. 7) is uniformly divided by the number of levels and numbered as li (i = 1, 2, . . ., 15). The input training data (rows denoted by u1 to u15 in Fig. 7) is obtained as follows. To acquire the value for a specific entry in the input data table, lookup the corresponding in-

7

dex value in the UD table located at the same row and column. For example, for the third entry in u2, the index value in row 2 column 3 of the UD table is 11 (see Fig. 7). Under column ‘‘li’’, the 11th level is l11 = 1.29. Add this li value to the sampling central value in the row denoted by u along the same column to get the input data value. This is repeated for all other entries. The sampling central values are initially assigned to be zeros (i.e., at the mean-value point in U-space) for the first iteration, and then shifted to the tentative DP for subsequent iterations. In addition, the sampling range need not stay constant at [3, 3] for all iterations. There is some advantage in progressively reducing the width of the sampling range in subsequent iterations. This refinement is explained in Section 5.

Fig. 8. The procedure of the probabilistic analysis of tunnel support system using ANN-RS and uniform design.

8

Q. Lü et al. / Tunnelling and Underground Space Technology 32 (2012) 1–18

When the training inputs uj are obtained, the xj in original X-space are computed from Eqs. (7) and (8), with which the performance function g(xj) (i.e., g(uj)) at each sampling point xj is evaluated using Eqs. (1)–(3), depending on the failure mode investigated. The above procedures for training data preparation are coded as a VBA macro in Excel for automatic execution. In view of the above-mentioned drawbacks in the existing research, an iterative algorithm is proposed herein in which the total number of sampling points is gradually increased and hence automatically determined by the convergence requirement of the reliability analysis. The proposed framework of the probabilistic analysis of ground-support system using ANN and UD is shown in Fig. 8. The basic idea of this iterative procedure is similar to that based on the polynomial RSM presented by Lü and Low (2011) and Lü et al. (2011), in which a quadratic polynomial is employed as the response surface (referred to as polynomial-RS) to fit the real LSS. The key difference is that in this paper, the polynomial function is replaced by the ANN-RS. The merits and advantages of the proposed approach will be illustrated in Section 5. 5. Probabilistic analyses for circular tunnel using the proposed approach In this section, the probabilistic analyses for the supported circular tunnel of Fig. 2 are carried out using the proposed procedures of ANN-RS, with respect to the three failure modes (i.e., performance function Eqs. (1)–(3)). The distance Ls from the support installation position to the face is assumed to be 2 m. The influence of the support installation position on the failure probability of different modes will be studied later. To establish the ANN-RS, 15 input training data uj (j = 1, 2, . . ., 15) are initially sampled around the mean-value point within a range of [3, 3] by means of UD. The performance function g(xj) (i.e., g(uj)) is then evaluated at each sampling point based on the xj values that are transformed from uj via Eqs. (7) and (8), to produce the output of the training data. Based on the input and output training data sets of uj and g(uj), the ANN is trained to establish a ~ðuÞ. The reliability index b and the correspondtentative ANN-RS g ing DP ut are then computed by using Microsoft Excel’s built-in constrained optimization routine Solver based on this tentative response surface. In subsequent iterations, eight additional sampling points are selected using the tentative DP ut as the center. Based on a total of 23 training data, a new ANN-RS is constructed for computing a new b and a new DP. Generally, repeating the iterative procedure for about four to six times, one would arrive at a converged b and a final ANN-RS, as shown in Table 2. The final ANNRS is obtained based on a total of 55 sampling points: 15 points from the first iteration, and 8 points from each of the subsequent five iterations.

The sampling number N of the first iteration is 15 for the case in hand (Table 2) but it can be as many as 30 for the practical application of a rock cavern in Section 6. This number should not be too small (compared with the number of random variables) because too few training data will hinder the ANN training process. In subsequent iterations, the sampling size is decreased from 15 to 8 to reduce the amount of real performance function evaluations. In Table 2, the sampling range k is progressively decreased from 3 to 0.02 to enhance the accuracy of the probabilistic results, especially for SORM. As the reliability index progressively converges, shrinking the sampling range will concentrate the new sampling points in the vicinity of the final DP, which will improve the fitting precision of the final ANN-RS in the region around the DP. This region is of most importance to the curvature fitting algorithm of SORM assessment. The k value is allowed to go as low as 0.02 (see Table 2) because for the circular tunnel the performance function is evaluated from an analytical solution. If numerical procedures are involved in performance function evaluation, a very small sampling range is not desirable because different sampling points may end up with identical performance function values and result in repetitions of input data. In such cases, sampling range should be sufficiently large, as shown in the cavern example in the next section in which the minimum k value is 0.1. The validity of the trained neural network is verified by comparing the predictions from ANN-RS with the real performance function evaluations at 30 testing points sampled by UD centering at the final DP within a range of [0.5, 0.5] and the mean-value point within a range of [3, 3], respectively, as shown in Fig. 9. A good correlation is observed between the predictions from ANN-RS and the actual values computed from the analytical solutions for both scenarios in Fig. 9, which implies that the ANN model obtained has good generalization quality and is able to represent the actual LSS in the region of interest. 5.1. Compared with MCS and polynomial RSM Table 3 shows the results obtained from different methods with respect to the three failure modes. The results from MCS are the average values of six separate simulations with Latin hypercube sampling (each comprising 200,000 trials), which provides the benchmark for comparison. The major drawback of MCS is the large number of trials required for problems with very small failure probability. For example, even though an analytical solution is used for the case in Table 3, a MCS of 200,000 trials will take more than 2 h because numerical procedures including Newton–Raphson method and fourth-order Runge–Kutta method are involved in the deterministic analyses of the ground-support interaction (Carranza-Torres, 2004). Even with 200,000 trials, the results still have a COV of 7.49% for the performance function g1, based on pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi the formula ð1  P f Þ=NP f (Shooman, 1990). In contrast, SORM

Table 2 Reliability results at each iteration for the example of a circular tunnel with respect to three failure modes in which Ls = 2 m. Iteration 1

Iteration 2

Iteration 3

Iteration 4

Iteration 5

Iteration 6

g1

b Pf,FORM (%) Pf,SORM (%)

2.793 0.261 0.213

3.056 0.112 0.092

3.079 0.104 0.086

3.079 0.104 0.083

3.079 0.104 0.087

– – –

g2

b Pf,FORM (%) Pf,SORM (%)

2.494 0.631 0.630

2.611 0.451 0.478

2.619 0.441 0.479

2.619 0.441 0.470

2.619 0.441 0.459

2.619 0.441 0.455

b Pf,FORM (%) Pf,SORM (%) Sampling number, N Sampling range, k

2.351 0.935 0.919 15 3

2.437 0.740 0.738 8 0.5

2.472 0.671 0.669 8 0.2

2.473 0.670 0.740 8 0.1

2.473 0.670 0.730 8 0.05

2.473 0.670 0.731 8 0.02

g3

9

Q. Lü et al. / Tunnelling and Underground Space Technology 32 (2012) 1–18

Prediction for g2 using ANN-RS, g2,ANN

(a) 0.010

0.005

presented in Table 3. The polynomial RSM includes first-order and second-order polynomial-RS with and without cross terms. Generally, the results from SORM based on ANN-RS and secondorder polynomial-RS with cross terms (RSM-3) are in very good agreement with those obtained from MCS. However, the total number of the sampling points needed for ANN-RS is only about 1/3  1/2 of that for RSM-3 (Table 3). The first-order polynomialRS (RSM-1) can also achieve a converged solution of the same reliability index and the same DP values, and requires as few sampling points as the ANN-RS for the case in hand. However, more samplings are needed to provide quadratic approximations at the converged b and the DP to allow SORM analysis. Note that when the second-order polynomial-RS without cross terms (RSM-2) is used, inaccurate and inconsistent outcomes are obtained for different performance functions. As shown in Table 3, for the support capacity criterion (g1), the Pf from SORM on RSM-2 are still acceptable. However, for the plastic zone criterion (g3) the results from SORM on RSM-2 are worse than that from FORM; it has a relative percentage error of 43.7% when compared with MCS. For the tunnel convergence criterion (g2), the quadratic polynomial-RS without cross terms is not able to obtain a converged result. In summary, for the specific example illustrated herein, the proposed method of using FORM/SORM on the ANN-RS shows its efficiency and accuracy in the reliability results when compared with the MCS and the polynomial RSM.

30 testing points sampled by UD centered at the DP in a range of [-0.5,0.5], with respect to the tunnel convergence criterion, g2

0.000

g2,ANN=0.993g2+0.029

-0.005

R2=0.999

-0.010 -0.010

-0.005

0.000

0.005

0.010

Actual g2 from real limit state function

Prediction for g2 using ANN-RS, g2,ANN

(b)

0.0225

0.0150

g2,ANN=1.003g2+0.215 R2=0.982

5.2. Performing MCS on obtained ANN-RS

0.0075

0.0000

30 testing points sampled by UD centered at the mean value point in a range of [-3,3], with respect to the tunnel convergence criterion,g2

-0.0075

-0.0150 -0.0150

-0.0075

0.0000

0.0075

0.0150

0.0225

Actual g2 from real limit state function Fig. 9. Comparison of the predictions for g2 using ANN-RS with those from real limit state function in the region centered (a) at the design point within a range of [0.5, 0.5], and (b) at the mean value point in a range of [3, 3].

based on ANN-RS needs only few minutes to obtain almost the same results. The relative percentage errors of SORM on ANN-RS to MCS are only 2.2% for g1, 0.9% for g2 and 0.4% for g3, computed from (Pf,SORM  Pf,MCS)/Pf,MCS  100%. For comparison, Pf results from direct FORM (based on bFORM of the actual performance functions and Eq. (9)) and the polynomial RSM using the procedures proposed in Lü et al. (2011) are also

As mentioned above, MCS can be very time-consuming when a large number of trials are required and the performance function assessment involves iterative numerical procedures. However, if MCS is implemented on the simple and explicit response surface instead of the real performance function, the computational cost will be significantly reduced. Table 4 compares the failure probability obtained from MCS based on the real performance function computations, the ANNRS predictions and the second-order polynomial-RS with cross terms (RSM-3) predictions. As shown, the MCS results based on the ANN-RS model agree well with the benchmark results of MCS based on the real performance function evaluations. Furthermore, a simulation of 200,000 trials on ANN-RS takes only half a minute. However, if the second-order polynomial response surface is used to perform MCS, one needs to be cautious with the results. As shown in Table 4, the failure probability corresponding to performance functions g2 and g3 are significantly different from those benchmark results on the real performance function evaluations. The reason can be intuitively appreciated from Fig. 10 in which the approximate LSS of g2 generated from the polynomial-RS has two curves in the 2D U-space of random variables GSI and rci. The lower curve is practically coincidental with the true LSS; the upper curve does not fit the true LSS at all and is redundant. Had the two curves been plotted over a sufficiently wide domain, one

Table 3 Comparison of reliability results from different approaches for the circular tunnel with respect to three failure modes in which Ls = 2 m. Direct FORM

g1 g2 g3

Pf (%) N Pf (%) N Pf (%) N

0.104 – 0.441 – 0.670 –

MCS

0.089 2  105 0.451 2  105 0.734 2  105

RSM-1

RSM-2

FORM

ANN-RS SORM

FORM

FORM

SORM

FORM

RSM-3 SORM

0.104 47 0.441 55 0.670 55

0.087 47 0.455 55 0.731 55

0.104 40 0.441 56 0.670 48

0.104 60 Non-convergence

0.085 60

0.671 75

1.055 75

0.104 108 0.441 144 0.670 144

0.087 108 0.447 144 0.728 144

In the table, Direct FORM represents Pf computed from bFORM of the actual performance functions and Eq. (9); RSM-1 denotes first-order polynomial-RS, RSM-2 denotes second-order polynomial-RS without cross terms and RSM-3 denotes second-order polynomial-RS with cross terms; N is the total amount of the sampling points.

10

Q. Lü et al. / Tunnelling and Underground Space Technology 32 (2012) 1–18

Table 4 Comparison of the failure probability from MCS on analytical solution, ANN-RS and second-order polynomial-RS with cross terms (RSM-3).

a

Method

g1 (%)

g2 (%)

g3 (%)

MCS on real PerFunca MCS on ANN-RS MCS on RSM-3

0.089 0.091 0.089

0.451 0.446 39.25

0.734 0.731 1.027

PerFunc denotes performance function.

second-order derivatives of g~ðuÞ at the DP, where the response surface fits the true LSS almost perfectly; the redundant section of the response surface is not taken into account in FORM/SORM computations. 5.3. Influence of the support installation position Fig. 11 shows the influence of the support installation position on the failure probability of the three performance functions in

6.0

Actual LSS 4.5

u3 (σci)

3.0

Polynomial-RS Pseudo response surface from polynomial-RS

safe domain

1.5

True response surface 0.0

Mean-value point -1.5

failure domain

True LSS

Design point

-3.0

Tunnel convergence criterion,g2 -4.5 -4.5

-3.0

-1.5

0.0

1.5

u1 (GSI) Fig. 10. Actual limit state surface and the polynomial response surface in U-space of GSI and rci, with respect to the tunnel convergence criterion, g2.

would see that both are part of the same tilted and elongated ellipse. All points (i.e., combinations of GSI and rci) on the ellipse sat~ðuÞ ¼ 0, where g ~ðuÞ is the second-order polynomial response isfy g surface function. However, since the response surface is constructed based on sampling points within a small region around the DP, only the lower curve as shown in Fig. 10 can be expected to satisfy g(u) = 0, where g(u) is the true performance function. Those points on the upper curve should give g(u) > 0 because they obviously lie in the safe domain. For the case in hand, the lower and upper curves are near each other and they encompass the mean-value point. This creates an incorrect safe domain within the boundary of the ellipse and an incorrect failure domain outside of it. When MCS is carried out on g~ðuÞ, part of the safe domain (the region above the upper curve in Fig. 10) is erroneously regarded as unsafe. This explains why the failure probability based on MCS on g~ðuÞ is much higher than that based on MCS on g(u). Mathematically the quadratic polynomial function in a 2D space could be a parabola, hyperbola or ellipse of which the form and the shape depend on the coefficients. Therefore second-order polynomial response surface may have different properties for different cases. For g3, although the failure probability of 1.027% from MCS using polynomial-RS has about 40% relative percentage error to the true value, the overestimation is less serious than it is for g2. This is because for this case the redundant portion of the response surface is very far away from the actual LSS and therefore imposes relatively less influence on the results of MCS. As for the first performance function g1, the results from MCS on the polynomial-RS is the same as the benchmark result, because for this case the redundant response surface is in the failure domain and also very far away from the actual LSS. As such, there is negligible effect on the MCS result. It is worth mentioning that FORM/SORM analyses on g~ðuÞ are not affected by the above pitfall. This is because in FORM/SORM, the probability of failure is computed based on the first-order/

Fig. 11. The influence of the support installation position on the failure probability Pf with respect to different failure modes of the ground-support system.

11

Q. Lü et al. / Tunnelling and Underground Space Technology 32 (2012) 1–18 Table 5 Influence of correlation coefficients on reliability index and failure probability with respect to different performance function, in which Ls = 2 m. Correlation coefficient

1 2 3 4 5

Support capacity criterion, g1

Convergence criterion, g2

Plastic zone criterion, g3

qGE

qrE

qkp

bFORM

Pf,FORM (%)

DPf (%)

bFORM

Pf,FORM (%)

DPf (%)

bFORM

Pf,FORM (%)

DPf (%)

0 0.6 0 0 0.6

0 0 0.6 0 0.6

0 0 0 0.2 0.2

2.939 2.989 2.967 3.019 3.079

0.165⁄ 0.140 0.151 0.127 0.104

– 15.0 8.5 23.0 36.9

3.267 2.935 2.840 3.267 2.619

0.054⁄ 0.167 0.226 0.054 0.441

– 207.2 315.8 0 712.3

2.399 2.434 2.447 2.399 2.473

0.822⁄ 0.748 0.721 0.822 0.670

– 9.0 12.3 0 18.5

In the table, DPf is computed by ðP f  P f Þ=P f  100%; in which P f denotes the values in cells with⁄.

Eqs. (1)–(3). For the support capacity criterion, the failure probability decreases with increasing distance Ls (Fig. 11a). In contrast, the influences of the distance Ls on the tunnel convergence criterion and plastic zone criterion have opposite tendencies (Figs. 11b and c). The support installation position has a very significant influence on the failure probability of all three failure modes because the convergence that occurred prior to the support installaD tion (uin r ) and the final equilibrium pressure on the support (ps ) strongly depends on the support installation position (see Fig. 1). Take the support capacity criterion as an example. When the support is installed at the position of the face itself (i.e., Ls = 0), as the face advances, the probability of exceeding the support capacity is

(a) 25

54.3%. When the support is installed at a distance of 2 m behind the face, the probability is dramatically reduced to 0.089%. Note that in Figs. 3 and 11a the reliability index becomes negative when the distance Ls is less than a certain critical value Ls,cr (0.134 m in this case). This implies that the mean value point is already in the unsafe domain before Solver search for the minimum b when Ls < Ls,cr, and the reliability index so obtained should be treated as negative (i.e., Pf > 50%), as explained in Low (2008). Further, when Ls < 0.51 m in Fig. 11b and when Ls < 0.49 m in Fig. 11c, the probability of excessive tunnel convergence and exorbitant plastic zone become constant. This is because the SCC is assumed to conform to an ideal elastic-perfectly-plastic behavior as shown in Fig. 1. Thus, the final equilibrium pressure on the support should not exceed its capacity when installed very close to the face.

Plastic zone criterion, g3, ρGΕ =0, ρkp=0, Ls=2 m 5.4. Influence of correlation coefficients

Erm/ GPa

20

Limit state surface ρσΕ = 0

15

ρσΕ = 0.6

10 Design point

5

Mean value point 0 40

60

80

100

120

140

160

180

σci/ MPa

(b) 16

Tunnel convergence criterion, g2, ρGΕ=0, ρkp=0, Ls=2 m

Erm/ GPa

13

ρσΕ = 0.6

Mean value point

10 ρσΕ = 0

7

4

Limit state surface

Design point 1

0

40

80

120

160

200

240

σci/ MPa Fig. 12. Interpretation of influence of correlation coefficients on the reliability index with respect to (a) positively inclined LSS of plastic zone criterion g3; (b) negatively inclined LSS of tunnel convergence criterion g2.

In the preceding analyses, the rock mass elastic modulus Erm is positively correlated with GSI and rci respectively and the support stiffness ks is positively correlated with the support capacity pmax s (Fig. 3). The correlation coefficients in the correlation matrix R of Fig. 3 are approximated from MCS as described in Section 3. The analyses below investigate the influence of such correlations on the computed reliability index and failure probability. As shown in Table 5, assuming that the distance from the face to the support installation position Ls = 2 m, for the support capacity criterion g1 and the plastic zone criterion g3, the reliability index b for the case where all random variables are uncorrelated (scenario 1) is smaller than those with positive correlations (scenarios 2–5). On the other hand, for the convergence criterion g2, the reliability index based on uncorrelated random variables is larger than those based on correlated random variables. This complicated influence of the correlation coefficients on the reliability results can be intuitively appreciated using the perspective of the expanding dispersion ellipse. Fig. 12 plots the equivalent dispersion ellipses and the limit state surfaces in the space of rci and Erm, by setting the other random variables at their DP values (since the original seven-dimensional equivalent hyper-ellipsoid can only be perceived in the mind’s eye). As shown in Fig. 12, the positively correlated rci  Erm dispersion ellipse has a positively inclined major axis in the coordinate system of rci and Erm, whereas the uncorrelated dispersion ellipse has its major and minor axis parallel to the coordinate axes. Thus, for a positively inclined LSS (i.e., g1 and g3 in this example; only g3 is presented in Fig. 12a), the positively correlated dispersion ellipse will have to expand more (i.e., larger reliability index) than the uncorrelated dispersion ellipse before the point of tangency (DP) on the LSS is reached. As for the negatively inclined LSS, g2, the positively correlated dispersion ellipse will have a smaller reliability index than the uncorrelated one. In addition, in Table 5 the degree of influence of the three correlation coefficients (i.e., qGE, qrE and qkp) on the reliability results

12

Q. Lü et al. / Tunnelling and Underground Space Technology 32 (2012) 1–18

Table 6 Reliability index and failure probability using beta distribution for GSI, mi and p0 in which Ls = 2 m.

g1 g2 g3

b

Pf,FORM

3.075 2.613 2.461

0.105 0.449 0.692

(%)

Pf,SORM 0.086 0.448 0.719

(%)

Pf,MCS

(%)

0.088 0.448 0.724

are different with respect to different failure modes, which mainly depends on the sensitivity of the reliability index to the random variables. For example, for the support capacity criterion g1, the correlation coefficient of most significant influence is qkp between ks and pmax , although its value 0.2 is less than qGE and qrE of 0.6. s This is because for the performance function g1 the reliability index is most sensitive to pmax , based on the n value in Fig. 3, as diss cussed in Section 4.1. However, the qkp of same value has no effect on the reliability index for the convergence criterion g2 and plastic zone criterion g3, as shown in Table 5 scenario 4. This is because the value of the DP n with respect to these two performance functions are (0.848, 0.765, 1.179, 2.437, 1.429, 0.048, 0.010) and (1.093, 1.157, 1.573, 0.137, 0.850, 0.083, 0.017), respectively, which indicates that the resulting reliability index is not sensitive to pmax for its DP value is close to its mean value s (i.e., ni value is close to zero), thus the correlation between ks and pmax has no influence on the reliability index and the inferred s failure probability.

5.5. Reliability analysis using four-parameter beta distribution In the above analyses, the two-parameter normal distribution, which is symmetrical and theoretically has a possible range from negative infinity to positive infinity, is used for GSI, mi and p0 (Table 1). However, for a parameter that obeys normal distribution and admits only positive values, the problem of the parameter taking negative values (which is irrational) is often encountered in reliability analysis and MCS especially when the COV of the parameter is 0.25 or bigger (e.g., Mollon et al., 2009; Lü et al., 2011). The lognormal distribution is usually suggested in lieu of the normal distribution to exclude the negative values. Alternatively, the four-parameter (a1, a2, min, max) beta distribution is more versatile and can also be used instead of the normal distribution. The parameters of the beta distribution

are determined from the mean value l and standard deviation r of a normal distribution (e.g., Evans et al., 2000) with min = 0 to exclude the negative values of random variables, and a1 = a2 to make the distribution symmetrical. For example, if the beta distribution is used with parameters of (49.5, 49.5, 0, 70) for GSI, (21.7, 21.7, 0, 20) for mi and (21.7, 21.7, 0, 54) for p0. The reliability index and the failure probability from FORM, SORM and MCS (Table 6) are close to those in Table 3 in which the normal distribution is used for GSI, mi and p0. Thus the beta distribution is a good substitute for the normal distribution and has advantage over the latter in being bounded in the positive realm. 6. Application for an underground rock cavern The application considered in this section is an example of an underground cavern that is constructed in a good quality rock mass at a depth of 880 m, adapted from a technical note by Hoek (2011) on cavern reinforcement and lining design. The deterministic numerical model was developed using the finite difference code FLAC3D (www.itascacg.com), with which the performance functions are evaluated. The probabilistic assessment is conducted with respect to two performance functions. 6.1. Description of the deterministic numerical model The cavern has a span of 16 m and a height of 25.6 m and is excavated in a homogeneous rock layer of jointed granodiorite which has deterministic properties listed in Fig. 13. The in situ stress field has been measured and it is found that the vertical stress is equal to 23.8 MPa (i.e., equal to overburden stress). The maximum horizontal stress perpendicular to the cavern axis is 1.2 times of the vertical stress while the intermediate horizontal stress is equal to the vertical stress. The material model used in this paper for jointed rock mass is an elastic-perfectly-plastic model of a fixed yield surface defined by the generalized Hoek–Brown failure criterion and non-associated plastic flow rule with a zero dilation angle. Hardening or softening is not considered. The cavern is excavated from an adit in the center of the heading. This adit will then be slashed out on either side to form the complete heading and this will be followed by benching down in six stages to excavate the body of the cavern, as shown in Fig. 14.

Fig. 13. Cross-section of the underground cavern and properties of rock mass and support.

Q. Lü et al. / Tunnelling and Underground Space Technology 32 (2012) 1–18

13

Fig. 14. Sequential excavation and installation of reinforcement and lining, after Hoek (2011).

In order to stabilize the highly stressed rock mass surrounding the cavern, 6 m long 782 kN capacity and 9 m long 1043 kN capacity cable tiebacks are installed as rock mass reinforcement in the cavern roof and sidewalls, respectively, 2 m away from the face in each stage (Fig. 14). These tiebacks are installed on a grid spacing of 1.5  1.5 m with the final 30% end of the cable grouted to form an anchor. The remaining length between the anchor and the faceplate is de-bonded from the grout by enclosing it in a plastic sleeve and therefore is free to deform independently of the rock. Note that the reinforcing mechanism of the cable tieback is different from that of the fully grouted rockbolt. For the latter, the deformation of the bolt cannot be separated from the deformation in the adjacent rock since they are bonded together by grout (see Stille et al., 1989; Indraratna and Kaiser, 1990; Peila and Oreste, 1995, for example). Furthermore, a 30-cm thick shotcrete layer with a 28-day uniaxial compressive strength of 40 MPa and a deformation modulus of 7000 MPa is applied as a final lining on the cavern roof and sidewalls (Fig. 14), sufficiently far away (e.g., 50 m) from the face so that it is not damaged by ongoing deformation as the cavern is excavated downward (Hoek, 2011). To compute the final cavern wall displacement and the mobilized axial forces in the tiebacks for reliability analysis, the following three major deterministic numerical procedures are performed according to the principle of CCM. (1) Determine the amount of cavern wall displacement prior to support installation using a 3D numerical model. In Hoek (2011), the empirical LDPs proposed by Vlachopoulos and Diederichs (2009) were used to establish the relationship between cavern displacement and distance from the face. However, the formulation for the LDP calculation in Vlachopoulos and Diederichs (2009) paper were originally developed for a circular tunnel and the resulting LDPs significantly depend on the ratios of plastic zone radius to tunnel radius. Applying these formulas for an almost square cavern makes the results merely approximate. This will inevitably introduce errors in the results of ground-support interaction analyses and eventually influence the successive performance functions assessment in reliability analysis. For this reason, a 3D numerical model was developed using FLAC3D (Fig. 15a) to compute the LDPs for the case in hand. Symmetry was utilized to reduce model size and save computing time. As shown in Fig. 15a, the whole model was

discretized into 28,080 elements using a radial meshing scheme and was constrained on the boundaries in the normal direction. (2) Using a 2D numerical model and the core material softening technique, determine the modulus reduction that yields the computed amount of cavern displacement prior to support installation in step (1). The core material softening technique is chosen herein because it is more versatile than the field stress relaxation method for multiple adjacent tunnelling issues such as the example under discussion (Hoek, 2011). The 2D numerical model has identical dimension and mesh as the above 3D model in x–z plane. The discretization of this 2D model resulted in a total of 1170 elements. (3) Build a refined 2D model and reduce the core modulus based on the results in step (2) to simulate the advancement of the face. Install tiebacks and shotcrete according to the excavation sequence and determine the final cavern displacement on the wall and the mobilized axial forces in tiebacks. As shown in Fig. 15b, the dimension of this refined 2D model is identical to the 3D model in Fig. 15a in x and z direction and is 1.5 m in y direction which is in accordance with the longitudinal spacing of the tiebacks. All boundaries are constrained in the normal direction which is the same as the preceding two procedures. The whole model was refined with higher mesh density for calculation accuracy and was discretized into 9018 elements. To compute the performance function values at the sampling points, all numerical simulations are automatically driven by a coded procedure using FISH language in FLAC3D. 6.2. Performance functions and random variables Two performance functions are investigated for this example:

g 4 ðxÞ ¼ F max  Ft t

ð10Þ

 ux g 5 ðxÞ ¼ umax x

ð11Þ

The first performance function is tieback capacity criterion and becomes negative (i.e., reinforcement provided by tieback ruptures or yields) when the maximum mobilized axial forces Ft is larger than the tieback capacity F max (782 kN in the roof and 1043 kN in t the wall). In practice, failure of an end-anchored tieback may occur

14

Q. Lü et al. / Tunnelling and Underground Space Technology 32 (2012) 1–18

(a)

(b)

70 m

Tieback 160 m 160 m

26 m

z y

64 m z

x

120 m 80 m

x 80 m

Fig. 15. (a) 3D numerical model for determining the amount of cavern displacement prior to support installation; (b) refined 2D numerical model for evaluating the cavern final displacement and the mobilized tieback axial forces based on the results obtained in procedures (1) and (2) described in Section 6.1.

in two modes. The first one is due to the yielding of the cable material and the second one is due to an insufficient capacity of the cable-rock anchorage (Wang and Garga, 1992). For simplicity, only the first failure mode is considered in this study. The capacity of the anchorage is supposed to be always larger than the cable tensile strength. The second performance function is cavern displacement criterion and becomes negative as the computed maximum horizontal displacement ux, which generally occurs at the center of the vertical sidewall, exceeds a threshold of umax ¼ 150 mm. x The tieback bearing capacity criterion represents the ultimate limit state whist the cavern displacement criterion is akin to serviceability limit state. The rock mass properties and the in situ stress are regarded as basic random variables with statistical properties as shown in Table 7. In addition, the parameter of Erm is assumed to be positively correlated with GSI and rci respectively, with a correlation coefficient of 0.6 in each case.

6.3. Reliability analyses for the rock cavern using the proposed approach The reliability analyses of this rock cavern are performed with respect to Eq. (10) and (11) using the approach proposed in this paper. The performance functions at the sampling points are evaluated using the numerical procedures described in Section 6.1. The training data for establishing the ANN-RS are sampled at the mean-value point and the tentative DPs by means of uniform design. The ANN model has a 7-20-1 structure, same as that of the circular tunnel example. Table 7 Statistical properties of rock mass and in situ stresses. Random variables

Distribution

l

r

COV

GSI mi rci (MPa) Erm (MPa) rzz (MPa) krx kry

Normal Normal Lognormal Lognormal Normal Normal Normal

50 23 125 17 23.8 1.2 1.0

5 2.3 18.75 4.25 2.38 0.12 0.1

0.10 0.10 0.15 0.25 0.1 0.1 0.1

In the table, krx and kry represent the ratios of in situ stress in x direction and y direction to the vertical overburden stress rzz.

For this application, generally five to six iterations would produce a converged solution. Table 8 shows the reliability results for the cavern with respect to the two failure modes. The sampling range and sampling size for each iteration are also presented in this table. The DP values in n of Table 8 imply that the probability of exceeding the capacity of the tiebacks in the wall is most sensitive to the GSI, rci, Erm and rzz while less sensitive to the mi, krx and kry. For the tiebacks in the roof and for the cavern wall displacement, the reliability index and the inferred failure probabilities are sensitive to all random variables except for mi and kry. The sensitivities reflected by the DP n show the relative importance of the uncertainty in each random variable and therefore can be used to identify random variables whose uncertainty has significant influence on the reliability index (e.g., Erm has great influence on all three failure modes in Table 8) and those which can be replaced by deterministic values (e.g., mi and kry). The maximum mobilized axial force Ft of the tiebacks installed in the wall as computed by FLAC3D is equal to 469 kN at the mean value point, which is much less than its capacity F max of t 1043 kN and implies a factor of safety (FOS) of 2.22 calculated by the formula FOS = F max =F t . However, even with a FOS of 2.22 t there is still a probability of exceeding the tieback capacity of about 1% (from SORM in Table 8) due to the uncertainties in the parameters. The tieback with the maximum mobilized axial force is located at the foot of the wall in bench 1, because this tieback is activated in the early stage of the excavation and is therefore significantly influenced by the ongoing deformation as the cavern is excavated downward. The maximum mobilized axial force of the roof tiebacks (located near the roof center of the heading) is 269 kN with a corresponding FOS of 2.91 at the mean-value point. The probability of exceeding roof tieback capacity is about 0.156% from SORM. As for the cavern displacement criterion (Eq. (11)), the maximum horizontal displacement of the wall is 67.5 mm at the mean-value point and is close to 150 mm at the DP, both are located at the center of the wall (see Fig. 16). The computed probability of excessive horizontal displacement is about 0.531% from FORM and 0.546% from SORM (Table 8). Note that the CPU time for the performance function computation at each sampling point for this example ranges from 60 to 90 min (on a PC of Intel Core 2 Quad CPU [email protected] GHz and 3.0 GB of RAM). Thus the average total CPU time is about

15

Q. Lü et al. / Tunnelling and Underground Space Technology 32 (2012) 1–18 Table 8 Reliability results for the example of cavern with respect to different performance functions. Failure modes

b

Pf,FORM

Pf,SORM

Tieback capacity criterion in wall

2.355

0.926%

1.001%

Tieback capacity criterion in roof

3.026

0.124%

0.156%

Cavern displacement criterion

2.555

0.531%

0.546%

Final design point in U-space (u) and X-space (x)

u n x u n x u n x

Sampling range k and sampling numbers N for each of iteration Iteration 1 2 k 3 1.0 N 30 10

GSI

mi

rci (MPa)

Erm (GPa)

rzz (MPa)

krx

kry

1.481 1.481 42.60 1.608 1.608 41.95 1.458 1.458 42.71

0.115 0.115 22.74 0.208 0.208 23.47 0.071 0.071 22.84

1.263 1.263 102.4 1.202 1.202 103.3 1.364 1.364 100.9

0.870 2.107 9.817 1.322 2.386 9.161 1.059 2.254 9.469

0.978 0.978 26.08 1.194 1.194 26.60 0.925 0.925 25.96

0.176 0.176 1.221 1.382 1.382 1.365 0.746 0.746 1.290

0.016 0.016 0.998 0.003 0.003 1.001 0.048 0.048 0.995

3 0.5 10

4 0.3 10

5 0.2 10

6 0.1 10

Fig. 16. Horizontal displacement distribution at (a) the mean value point and (b) the design point (unit: m).

90–100 h depending on the performance function investigated and the iterations required. Monte Carlo simulations are not performed for this numerical example, because even 10,000 trials will take at least 10,000 h (417 days), and the result would still have more than 20% error in the failure probability with 95% confidence interval estimated by the equation in Shooman (1990). Thus for such time-consuming practical application that needs complex numerical procedures to evaluate the performance function, the proposed ANN-RS based approach is more efficient and feasible than MCS. 7. Summary and conclusions This paper proposes an efficient approach for probabilistic ground-support interaction analysis of deep rock excavation using the artificial neural network (ANN) and uniform design (UD). The three-layer feed-forward back-propagation network is employed as the response surface (referred to as ANN-RS) to approximate the actual limit state surface. The uniform design table is adopted to prepare the training data for determining the parameters of the ANN-RS via an iterative procedure. The failure probabilities with respect to three different failure modes are estimated from FORM

and SORM based on the fitted ANN-RS and compared with Monte Carlo simulations (MCS) and polynomial RSM. The efficiency and feasibility of the proposed approach are illustrated with two examples analyzed by an analytical solution and numerical procedures. Some major conclusions are summarized below. Using the random variables in U-space as basic variables and transforming them from U-space into X-space to establish ANN-RS makes the proposed procedure concise and convenient to perform first-order and second-order reliability analysis on the generated ANN-RS. Using uniform design to prepare the training data for establishing ANN-RS reduces the sampling size and ensure that the samples are uniformly distributed in the domain with high representativeness. The distance from the tunnel face when support is installed has a significant influence on the probability of the three failure modes under investigation. It was found that increasing the support installation distance Ls from the face will reduce the probability of exceeding the support capacity, but increase the probability of exceeding the permissible tunnel convergence and the probability of excessive plastic zone. Comparison between analyses using correlated and uncorrelated parameters indicates that the influence of the correlation

16

Q. Lü et al. / Tunnelling and Underground Space Technology 32 (2012) 1–18

on the reliability analysis results depends on the inclination of the LSS and the sensitivities of the reliability index to the correlated parameters. The four-parameter beta distribution is a good substitute of the normal distribution for excluding unrealistic negative values that may be encountered in reliability analysis or MCS. Compared with MCS, the proposed approach needs only a small number of deterministic analyses but is as accurate as MCS, and thus has advantage in terms of efficiency in probability calculation over MCS. This is illustrated in the two examples presented in the paper. Note that running MCS to analyze the rock cavern example is not practical since complex and time-consuming numerical procedures are often involved in the analyses. Furthermore, the proposed procedure is also more accurate and efficient than the second-order polynomial RSM for the circular tunnel example, and is therefore attractive to the designer who is conducting reliability-based analysis and design of ground support systems for underground excavations. The ANN and UD approach outlined in this paper was used to assess the reliability of underground excavations considering parametric uncertainties in the rock and support properties. The intrinsic variation in the design parameters are accounted for in a rational manner, and the risks associated with geotechnicalrelated hazards, such as support collapse or yielding, tunnel instability or squeezing, are quantitatively assessed. The information obtained can help designers and decision-makers to better understand the impacts of the various sources of uncertainty on the tunnel support design and to take effective countermeasures to reduce the underlying risks to an acceptable level. Therefore, the proposed method can play a significant role in tunnelling risk management (e.g., Eskesen et al., 2004). This paper confines itself to probabilistic uncertainty in which uncertain parameters are described by random variables. Other non-probabilistic uncertainties may also have considerable influences on the reliability of the tunnel support system, such as the irregular tunnel boundary due to overbreaks which can affect the stresses in the shotcrete lining (Borio and Peila, 2009). Fuzzy logic approaches that suggested by Sadeghi et al. (2010), Barpi and Peila (2012), for example, may be used to take into account such subjective and linguistically expressed uncertainty. Acknowledgments The authors would like to thank the Jurong Town Corporation of Singapore for providing research funding, which enable the first author (Lü Qing) to work as a post-doctoral research fellow and the second author (Chan Chin Loong) as a project officer at Nanyang Technological University. This research was also supported by Zhejiang Province Key Innovation Team Support Program (No. 2009R50050) and Zhejiang Province Special Scientific and Technological Project on ‘‘Methodology and modeling of tunnel excavation’’ (2010C13029). Appendix A. Details of ANN response surface, architecture and training The most widely used ANN is the multi-layer feed-forward back-propagation network (see Fig. 5). The neurons are the basic processing units and their outputs can be mathematically expressed as:

" aj ¼ f

n X wj;i ui þ bj

# ðA:1Þ

i¼1

where ui are the inputs, aj are the outputs, wj,i and bj are the weights and bias of the neuron j, f() is the transfer function (or activation function). The outputs of the neurons are then transmitted as inputs

to the neurons in the subsequent layer based on the same computation in Eq. (A.1). Finally, the final output from the neuron in the output layer is the network output g~ðuÞ. The transfer functions used in this paper are the logistic sigmoid function for the hidden layer and the pure linear function for the output layer, respectively. The logistic sigmoid function is expressed as:

f ðuÞ ¼

1 1 þ expðuÞ

ðA:2Þ

Eq. (A.2) has an interesting and useful feature, for its derivative can be expressed in terms of itself (Deng et al., 2005):

df ðuÞ ¼ f ðuÞð1  f ðuÞÞ du

ðA:3Þ

The number of neurons in the input and output layers are directly determined by the number of system inputs and outputs. However, the selection of the number of hidden neurons is more difficult because using too few neurons impedes the training process and prevents the correct mapping of inputs to outputs, while using too many neurons tends to overfit the data and requires more computation (Hagan et al., 1996). There is no general rule for determining the appropriate number of hidden neurons. A trial-and-error process is adopted in this study to determine the optimal number, starting with a small initial number of hidden neurons, then increasing the number and repeating the training until there is no significant improvement in the network performance (Goh, 1995). An ANN is essentially trained through adaption of their weights and biases using optimization methods. Most of these methods use either the gradient of the network performance or the Jacobian of the network errors with respect to the weights (Hagan et al., 1996). The gradient and the Jacobian are computed using the back-propagation algorithm. During the training phase, the ANN reads the input and output values in the training data set and optimizes the values of the weights and biases to reduce the mean square error between the predicted and the target values. The error is minimized through an iterative process until the network reaches a converged result. In the application of the back-propagation algorithm, two distinct passes take place in each training cycle: the input information propagated forward through the network and the errors at the output layer propagated backward through the network. Optimization of the weights and biases is then carried out using the back-propagated error. Details of the ANN training algorithm can be found in Hagan et al. (1996), for example. Some commercial codes such as MATLAB (http://www.mathworks.com) and many others have integrated these training algorithms for convenient use of the ANN. This study makes use of the MATLAB neural-network toolbox to train the ANN and to determine the parameters of the network. One issue that may occur during neural network training is overfitting which means that the trained network fit the training data very closely but yield large errors to new data. To overcome this problem and improve network generalization, the Bayesian regularization is introduced into the back-propagation algorithm. In the algorithm of Bayesian regularization, the weights and biases of the network are assumed to be random variables with specific distributions. The regularization parameters are related to the unknown variances associated with these distributions. Statistical techniques are then used to estimate these parameters. With Bayesian regularization, the resulting trained network shows good generalization qualities (Mackay, 1992; Cheng and Li, 2009). However, the Bayesian regularization algorithm generally works best when the network training data are confined within the range of [1, 1]. Thus the neural network training would be more efficient if the system inputs and outputs are normalized. In this paper,

Q. Lü et al. / Tunnelling and Underground Space Technology 32 (2012) 1–18

the sampling points uj of the independent standard normal basic random variables in U-space and the corresponding performance function g(uj) that is computed from either analytical solutions or numerical simulations are considered as inputs and outputs in the training data sets, respectively. All input training data are first normalized to the range of [1, 1] before presented to the network by the following equations:

u0i ¼ u0i;min þ

u0i;max  u0i;min ðui  ui;min Þ ui;max  ui;min

ðA:4Þ

where ui are the original inputs, u0i are the normalized inputs, ui,max and ui,min are the original maximum and minimum of ui, u0i;max and u0i;min are the upper and lower range values of u0i , with u0i;max =1 and u’i;min = 1. Thus Eq. (A.4) can be recast as:

ui ¼ ki u0i þ C i

ðA:5Þ

where the coefficient ki ¼ ðui;max  ui;min Þ=2 and the constant Ci = (ui,max + ui,min)/2. Following the same approach, the normalized output training data g’(u) can be expressed as:

gðuj Þ ¼ kg g 0 ðuj Þ þ C g

ðA:6Þ

where the coefficient kg ¼ ðg max  g min Þ=2 and the constant Cg = (gmax + gmin)/2 in which gmax and gmin are the original maximum and minimum of the g(uj). For a three-layer network that is adopted to approximate the real LSS for reliability analysis in this paper (Fig. 5), the ANN response surface can be explicitly expressed as:

g~ðuÞ ¼ kg

h n X X 1 w21;j f 1 w1j;i u0i þ bj j¼1

!!

! þ

2 b1

þ Cg

ðA:7Þ

i¼1

where the transfer function f1() is the logistic sigmoid function defined by Eq. (A.2). Thus the first-order and second-order derivatives of the ANNRS, g~ðuÞ, can be computed from Eq. (A.3) to Eq. (A.4) using the chain partial differentiation rule (Deng et al., 2005). The first-order derivative gives:

 h  h @ g~ X @ g~ @ g~0 @f 1 @u0i kg X ¼ ¼ w1j;i w21;j fj1 ð1  fj1 Þ 0 0 1 ~ @f @ui @ui @ui @ g k i j¼1 j¼1

ðA:8Þ

P  1 n 1 0 where fj1 ¼ f 1 i¼1 wj;i ui þ bj . And the second-order derivative gives:

" # h   @ 2 g~ @ g~ kg X ¼ w1j;i w21;j fj1 1  fj1 @ui @uk @uk ki j¼1 ¼

h h   i kg X w1 w1 w2 f 1 1  fj1 1  2fj1 ki kk j¼1 j;i j;k 1;j j

ðA:9Þ

Once the ANN-RS is obtained via the iterative procedure presented in Section 4, Eq. (A.8) and Eq. (A.9) may be used to compute the vector a and the matrix B in Fig. 4 for SORM assessment in Section 4.2. References Ang, A.H.S., Tang, W.H., 1984. Probability Concepts in Engineering Planning and Design. Decision, Risk and Reliability, vol. II. Wiley, New York. Baecher, G.B., Christian, J.T., 2003. Reliability and Statistics in Geotechnical Engineering. John Wiley, San Francisco. Barpi, F., Peila, D., 2012. Influence of the tunnel shape on shotcrete lining stresses. Comput.-Aided Civ. Inf. 27 (4), 260–275. Borio, L., Peila, D., 2009. How tunnel boundary irregularities can influence the stresses in a shotcrete lining. GEAM. Geoingegneria Ambientale e Mineraria 47 (2), 45–52. Breitung, K., 1984. Asymptotic approximations for multinormal integrals. J. Eng. Mech. 110 (3), 357–366.

17

Brown, E.T., Bray, J.W., Ladanyi, B., Hoek, E., 1983. Ground response curves for rock tunnels. J. Geotech. Eng. Div. 109 (1), 15–39. Cai, G.Q., Elishakoff, I., 1994. Refined 2nd-order reliability analysis. Struct. Saf. 14 (4), 267–276. Carranza-Torres, C., 2004. Elasto-plastic solution of tunnel problems using the generalized form of the Hoek–Brown failure criterion. Int. J. Rock Mech. Min. Sci. 41, 629–639. Carranza-Torres, C., Diederichs, M., 2009. Mechanical analysis of circular liners with particular reference to composite supports. For example, liners consisting of shotcrete and steel sets. Tunn. Undergr. Sp. Tech. 24 (5), 506–532. Carranza-Torres, C., Fairhurst, C., 2000. Application of the convergence– confinement method of tunnel design to rock masses that satisfy the Hoek– Brown failure criterion. Tunn. Undergr. Sp. Tech. 15 (2), 187–213. Chan, C.L., Low, B.K., in press. Practical second-order reliability analysis applied to foundation engineering. Int. J. Numer. Anal. Met. http://dx.doi.org/10.1002/ nag.1057. Cheng, J., Li, Q.S., 2008. Reliability analysis of structures using artificial neural network based genetic algorithms. Comput. Method. Appl. Mech. Eng. 197 (4548), 3742–3750. Cheng, J., Li, Q.S., 2009. A hybrid artificial neural network method with uniform design for structural optimization. Comput. Mech. 44 (1), 61–71. Cho, S.E., 2009. Probabilistic stability analyses of slopes using the ANN-based response surface. Comput. Geotech. 36 (5), 787–797. Choi, S.K., Grandhi, R.V., Canfield, R.A., 2006. Reliability-Based Structural Design. Springer, London. Deng, J., Gu, D.S., Li, X.B., Yue, Z.Q., 2005. Structural reliability analysis for implicit performance functions using artificial neural network. Struct. Saf. 27 (1), 25–48. Der Kiureghian, A., Lin, H.-Z., Hwang, S.-J., 1987. Second-order reliability approximations. J. Eng. Mech. 113 (8), 1208–1225. Ditlvsen, O., Madsen, H.O., 1996. Structural Reliability Methods. John Wiley & Sons, Chichester. Duncan Fama, M.E., 1993. Numerical modelling of yield zones in weak rocks. In: Hudson, J.A. (Ed.), Comprehensive Rock Engineering. Pergamon, Oxford, 49– 75pp. Eskesen, S.D., Tengborg, P., Kampmann, J., Veicherts, T.H., 2004. Guidelines for tunnelling risk management: international tunnelling association, working group no. 2.. Tunn. Undergr. Sp. Tech. 19 (3), 217–237. Evans, M., Hastings, N., Peacock, B., 2000. Statistical Distributions, third ed. J. Wiley, New York, 221. Fang, K.T., 1980. The uniform design: application of number-theoretic methods in experimental design. Acta Math. Appl. Sin. 3, 363–372. Fang, K.T., 1994. Uniform Design and Uniform Design Tables. Science Press, Beijing (in Chinese). Ghaboussi, J., Garrett, J.H., Wu, X., 1991. Knowledge-based modeling of material behavior with neural networks. J. Eng. Mech.—ASCE 117 (1), 132–153. Goh, A.T.C., 1994. Seismic liquefaction potential assessed by neural networks. J. Geotech. Eng. Div. 120 (9), 1467–1480. Goh, A.T.C., 1995. Back-propagation neural networks for modeling complex systems. Artif. Intell. Eng. 9 (3), 143–151. Goh, A.T.C., Kulhawy, F.H., 2005. Reliability assessment of serviceability performance of braced retaining walls using a neural network approach. Int. J. Numer. Anal. Met. 29 (6), 627–642. Gonzalez-Nicieza, C., Alvarez-Vigil, A.E., Menendez-Diaz, A., Gonzalez-Palacio, C., 2008. Influence of the depth and shape of a tunnel in the application of the convergence–confinement method. Tunn. Undergr. Sp. Tech. 23 (1), 25–37. Hagan, M.T., Demuth, H.B., Beale, M.H., 1996. Neural Network Design. PWS, Boston. Hoek, E., 1998. Reliability of Hoek–Brown estimates of rock mass properties and their impact on design. Int. J. Rock Mech. Min. Sci. 35 (1), 63–68. Hoek, E., 2011. Cavern Reinforcement and Lining Design (A Note Prepared for RocNews). . Hoek, E., Brown, E.T., 1980. Empirical strength criterion for rock masses. J. Geotech. Eng.—ASCE 106 (9), 1013–1035. Hoek, E., Diederichs, M.S., 2006. Empirical estimation of rock mass modulus. Int. J. Rock Mech. Min. Sci. 43 (2), 203–215. Hoek, E., Carranza-Torres, C., Corkum, B., 2002. Hoek–Brown failure criterion–2002 edition. In: Hammah, R., et al. (Eds.), 5th North American Rock Mechanics Symposium and 17th Tunnelling Association of Canada Conference: NARMSTAC 2002, Toronto, pp. 267–271. Hornik, K., Stinchcombe, M., White, H., 1989. Multilayer feedforward networks are universal approximators. Neural Networks 2 (5), 359–366. Indraratna, B., Kaiser, P.K., 1990. Analytical model for the design of grouted rock bolts. Int. J. Numer. Anal. Met. 14 (4), 227–251. Juang, C.H., Chen, C.J.X., Tien, Y.M., 1999. Appraising cone penetration test based liquefaction resistance evaluation methods: artificial neural network approach. Can. Geotech. J. 36 (3), 443–454. Kim, K., Gao, H., 1995. Probabilistic approaches to estimating variation in the mechanical properties of rock masses. Int. J. Rock. Mech. Min. 32 (2), 111–120. Koyluoglu, H.U., Nielsen, S.R.K., 1994. New approximations for SORM integrals. Struct. Saf. 13 (4), 235–246. Kuo, Y., Yang, T., Peters, B.A., Chang, I., 2007. Simulation metamodel development using uniform design and neural networks for automated material handling systems in semiconductor wafer fabrication. Simul. Model. Pract. Theor. 15 (8), 1002–1015. Li, H.Z., Low, B.K., 2010. Reliability analysis of circular tunnel under hydrostatic stress field. Comput. Geotech. 37 (1), 50–58.

18

Q. Lü et al. / Tunnelling and Underground Space Technology 32 (2012) 1–18

Li, A.J., Merifield, R.S., Lyamin, A.V., 2008. Stability charts for rock slopes based on the Hoek–Brown failure criterion. Int. J. Rock Mech. Min. Sci. 45 (5), 689–700. Low, B.K., 2008. Practical reliability approach using spreadsheet. In: Phoon, K.K. (Ed.), Reliability-Based Design in Geotechnical Engineering: Computations and Applications. Taylor & Francis, pp. 134–168 (Chapter 3). Low, B.K., Tang, W.H., 2007. Efficient spreadsheet algorithm for first-order reliability method. J. Eng. Mech. 133 (12), 1378–1387. Low, B.K., Zhang, J., Tang, W.H., 2011. Efficient system reliability analysis illustrated for a retaining wall and a soil slope. Comput. Geotech. 38 (2), 196–204. Lü, Q., Low, B.K., 2011. Probabilistic analysis of underground rock excavation using response surface method and SORM. Comput. Geotech. 38 (8), 1008–1021. Lü, Q., Sun, H.Y., Low, B.K., 2011. Reliability analysis of ground-support interaction in circular tunnels using response surface method. Int. J. Rock. Mech. Min. 48 (8), 1329–1343. Mackay, D.J.C., 1992. Bayesian interpolation. Neural Comput. 4 (3), 415–447. Melchers, R.E., 1999. Structure Reliability Analysis and Prediction, second ed. Wiley, New York. Merifield, R.S., Lyamin, A.V., Sloan, S.W., 2006. Limit analysis solutions for the bearing capacity of rock masses using the generalised Hoek–Brown criterion. Int. J. Rock Mech. Min. Sci. 43 (6), 920–937. Mollon, G., Dias, D., Soubra, A.H., 2009. Probabilistic analysis of circular tunnels in homogeneous soil using response surface methodology. J. Geotech. Geoenviron. 135 (9), 1314–1325. Oreste, P.P., 2003. Analysis of structural interaction in tunnels using the covergenceconfinement approach. Tunn. Undergr. Sp. Tech. 18 (4), 347–363. Peila, I.D., Oreste, I.P.P., 1995. Axisymmetrical analysis of ground reinforcing in tunnelling design. Comput. Geotech. 17 (2), 253–274. Rackwitz, R., 2001. Reliability analysis – a review and some perspectives. Struct. Saf. 23 (4), 365–395.

Sadeghi, N., Fayek, A.R., Pedrycz, W., 2010. Fuzzy monte carlo simulation and risk assessment in construction. Comput.-Aided Civ. Inf. 25 (4), 238–252. Sari, M., 2009. The stochastic assessment of strength and deformability characteristics for a pyroclastic rock mass. Int. J. Rock Mech. Min. Sci. 46 (3), 613–626. Schueremans, L., Van Gemert, D., 2005. Benefit of splines and neural networks in simulation based structural reliability analysis. Struct. Saf. 27 (3), 246–261. Shooman, M.L., 1990. Probability Reliability: An Engineering Approach, second ed. R.E. Krieger Publishing Company, Malabar. Stille, H., Holmberg, M., Nord, G., 1989. Support of weak rock with grouted bolts and shotcrete. Int. J. Rock Mech. Min. Sci. 26 (1), 99–113. Tvedt, L., 1983. Two Second-Order Approximations to the Failure Probability. Det Norske Veritas, Oslo. Vardakos, S.S., Gutierrez, M.S., Barton, N.R., 2007. Back-analysis of Shimizu tunnel no. 3 by distinct element modeling. Tunn. Undergr. Sp. Tech. 22 (4), 401–413. Vlachopoulos, N., Diederichs, M.S., 2009. Improved longitudinal displacement profiles for convergence confinement analysis of deep tunnels. Rock Mech. Rock Eng. 42 (2), 131–146. Wang, Y., Fang, K.T., 1981. A note on uniform distribution and experimental design. Chin. Sci. Bull. 26, 485–489. Wang, B.L., Garga, V.K., 1992. A numerical model for rock bolts. In: Kaiser, P.K., McCreath, D.R. (Eds.), Rock Support in Mining and Underground Construction. Taylor & Francis, Rotterdam, pp. 57–66. Wong, F.S., 1985. Slope reliability and response-surface method. J. Geotech. Eng. Div. 111 (1), 32–53. Zhang, L., Liang, Y.Z., Jiang, J.H., Yu, R.Q., Fang, K.T., 1998. Uniform design applied to nonlinear multivariate calibration by ANN. Anal. Chim. Acta 370 (1), 65–77. Zhao, Y.G., Ono, T., 1999. New approximations for SORM: Part 1. J. Eng. Mech. 125 (1), 79–85.