# Gaussian Process-Based Response Surface Method for Slope Reliability Analysis

Authors: Bin Hu; Guo-shao Su; Jianqing Jiang; Yilong Xiao
Journal: Advances in Civil Engineering
 (2019)

---

## Abstract

A new response surface method (RSM) for slope reliability analysis was proposed based on Gaussian process (GP) machine learning technology. The method involves the approximation of limit state function by the trained GP model and estimation of failure probability using the first-order reliability method (FORM). A small amount of training samples were firstly built by the limited equilibrium method for training the GP model. Then, the implicit limit state function of slope was approximated by the trained GP model. Thus, the implicit limit state function and its derivatives for slope stability analysis were approximated by the GP model with the explicit formulation. Furthermore, an iterative algorithm was presented to improve the precision of approximation of the limit state function at the region near the design point which contributes significantly to the failure probability. Results of four case studies including one nonslope and three slope problems indicate that the proposed method is more efficient to achieve reasonable accuracy for slope reliability analysis than the traditional RSM.


---

## Body

## 1. Introduction

Slope stability analysis is an important issue in geotechnical engineering. The application of probabilistic concepts to the slope stability analysis has been pursued over the past decades. On real slope engineering, slope reliability can be obtained via an implicit limit state function with a numerical procedure employing the limit equilibrium methods [1, 2] or finite element method [3, 4]. The classical reliability analysis methods including first-order reliability method (FORM) and second-order reliability method (SORM) [5–8] that require limit state function gradients with respect to the basic variables may encounter great difficulties when direct or analytical differential is not possible. The Monte Carlo simulation (MCS) can handle the implicit performance function. However, the MCS was notorious for its unendurable computational cost.The response surface method (RSM), which usually employs a polynomial function to approximate the unknown implicit performance function, has been used to the slope reliability analysis with the implicit performance function [4, 9–11]. However, the RSM is a rigidly nonadaptive regression technique in the statistical learning perspective [12]. Consequently, the RSM will become computationally impractical for problems involving many random variables, particularly when the involved random variables are mixed or statistically dependent. In addition, there is no guarantee that the fitted surface is a sufficiently close fit in all regions of interest.In recent decades, machine learning methods such as artificial neural networks (ANN) [13–17], radial basis function network (RBFN) [18], and support vector machine (SVM) [19, 20] have been widely used for the reliability analysis. These methods can describe nonlinear and complex interactions among variables in a system and are therefore more suitable for solving reliability problems with highly nonlinear implicit performance. However, when these machine learning methods are used, some difficulties will be encountered. For ANN, it is always difficult to obtain an appropriate network topology and the optimum hyperparameters. In addition, the ANN has some inherent drawbacks such as limitation in solving the small sample problem [21]. The SVM cannot avoid the blindness which is the common phenomenon in the artificial choice of the kernel function and its hyperparameter [21]. Hence, it is desirable to develop an efficient framework for slope reliability analysis.The Gaussian process (GP), based on statistical learning theory and Bayesian theory, is a newly developed machine learning technology [22]. In recent years, the GP has gained considerable attention in the machine learning community and has been successfully applied in solving nonlinear, small samples, and high-dimensional problems [23–28]. The GP has many advantages, such as it does not require a predefined structure, can approximate arbitrary function landscapes including discontinuities and multimodality, has meaningful hyperparameters, and includes a theoretical framework for obtaining the optimum hyperparameter self-adaptively. Furthermore, GP provides an uncertainty measure in the form of a standard deviation for predicted function values.In the present study, a new GP-based RSM for slope reliability analysis was put forward by combining the GP with the FORM. This paper is constructed as follows. The used machine learning technique, namely, Gaussian process regression, is briefly described in Section2. The GP-based RSM method is presented in detail in Section 3. Finally, four case studies are used to verify the flexibility and efficiency of the proposed method.


## 2. Gaussian Process

A Gaussian process is a collection of random variables, any finite set of which has a joint Gaussian distribution. A Gaussian process can be defined by its mean functionmx and the covariance function kx,x′ as(1)fx∼GPmx,kx,x′.There is a training setD of m observations, D=xi,yi∣i=1,…,m, where x is an input vector with n dimensions and y is a scalar output or target. By using Bayesian forecasting, the distribution py∗∣x∗,D of output y∗ given a test input x∗ and a set of training points D can be calculated. Using the Bayesian rule, the posterior distribution for the Gaussian process outputs y∗ can be obtained. By conditioning on the observed targets in the training set, the predictive distribution is Gaussian:(2)y∗∣x∗,X,y∼Ny¯x∗,σ¯x∗,where the mean and variance are defined by(3)y¯x∗=k∗TK+σn2I−1y,(4)σ¯2x∗=kx∗,x∗−k∗TK+σn2I−1k∗T,where the compact forms of the notation setting for the matrix of the covariance functions are k∗=KX,x∗, K=KX,X, σn2 is the unknown variance of the Gaussian noise.Assumingα=K+σn2I−1y, equation (3) can be seen as a linear combination of m covariance functions, each one centered on a training point, by writing(5)y¯x∗=∑j=1mαjkxj,x∗.The Gaussian process procedure can handle interesting models by simply using a covariance function with an exponential term:(6)kxp,xq=σf2exp−xp−xq22l2+σn2δpq,where l denotes the characteristic length scale, σf2 denotes the signal variance, and δpq denotes a Kronecker delta. This function expresses the idea the cases with nearby inputs will have highly correlated outputs. The GP employs a set of hyperparameters θ including the length scale l, the signal variance σf2, and the noise variance σn2. The hyperparameters θ can be optimized based on the log-likelihood framework:(7)L=logpy∣X,θ=−12yTC−1y−12logdetC−n2log2π.The log-likelihood and its derivative with respect toθ can be described as(8)∂L∂θ=−12trC−1∂C∂θ+12yTC−1∂C∂θC−1y,where C=K+σn2I.Hyperparametersθ are initialized to random values in a reasonable range and then optimized using an iterative, such as the conjugate gradient.More detailed theory of the GP can be found in the literature [24].


## 3. GP-Based RSM

### 3.1. Explicit Formation of Reliability Index Using GP

The performance function of a slope can be built as follows:(9)gX=Fs−1,where X denotes the random variables with mean μx and standard deviation σx. gX>0 indicates a stable slope, while gX<0 indicates a failed slope. gX=0 indicates that the slope is hovering between stable and unstable. Fs is the safety factor.Assuming that there is a pointx∗=x1∗,x2∗,…,xn∗T located on the limit state surface, which is called as the design point. Here, n is the dimension of x∗. gx is generally a nonlinear function. However, it can be linearized at x∗ by neglecting the second- or upper-order term:(10)ZQ=gx∗+∑i=1n∂gx∗∂xixi−xi∗.Assumingx is statistically uncorrelated, the mean value and standard deviation of ZQ are defined by(11)μZQ=gx∗+∑i=1n∂gx∗∂xiμxi−xi∗,σZQ=∑i=1n∂gx∗∂xi2σxi2.Reliability index can be calculated by(12)β=μZQσZQ=gx∗+∑i=1n∂gx∗/∂xiμxi−xi∗∑i=1n∂gx∗/∂xi2σxi2.According to the FORM, a tangent hyperplane is used to fit the limit state surface at the design point. Therefore, the most important step in the method is to obtain the design point. Reliability indexβ is defined as the distance from the origin to the design point. Hence, the coordinate of the design point x∗ is denoted as(13)xi∗=μxi+βσxicosθxi,i=1,2,…,n,where(14)cosθxi=−∂gx∗/∂xiσxi∑i=1n∂gx∗/∂xi2σxi2.According to equation (5), the performance function can be approximated by the GP model as(15)gx∗≈y¯x∗=∑j=1mαjkxj,x∗.The first-order partial derivatives of the approximate function can be given by(16)∂gx∗∂xi≈∂y¯x∗∂xi=∑j=1mαj∂kxj,x∗∂xi,where ∂kxj,x∗/∂xi=xi∗−xij/li2σf2exp−1/2l2·xj−x∗2.Substituting equations (15) and (16) into (17), the reliability index β is given as(17)β=∑j=1mαjkxj,x∗+μy+∑i=1n∑j=1mαj∂kxj,x∗/∂xiμxi−xi∗∑i=1n∑j=1mαj∂kxj,x∗/∂xi2σxi2.According to equations (15) and (16), (14) can be rewritten as(18)cosθxi=−∑j=1mαj∂kxj,x∗/∂xiσxi∑i=1n∑j=1mαj∂kxj,x∗/∂xi2σxi2.Substituting equations (17) and (18) into (13), one can obtain the coordinate of design point x∗ using GP approximation.The failure probability can be calculated as(19)pf≈Φ−β,where Φ⋅ denotes the cumulative distribution function of the standard normal variable.


### 3.2. Procedure of GP-Based RSM

The main difference between the GP-based RSM and the traditional RSM is that the former employs the GP to approximate the performance function and its first-order partial derivatives simultaneously. Moreover, from the viewpoint of machine learning theory, more knowledge means a better effect on the result of training or prediction of the machine learning model. Thus, the training samples, namely, knowledge set in the machine learning domain, is generated as a dynamic updating one to improve the performance of the GP model in the procedure of GP-based RSM.As shown in Figure1, reliability analysis using the GP-based RSM is described in detail as follows:Step 1. Assume initial values of the design point x∗=x1,x2,…,xn. Usually, mean values of the random variables can be selected as the coordinates of the design point x∗.Step 2. According to the experimental plan of the RSM developed by Bucher and Bourgund [29]; the training samples are generated according to the intersection of the axis xμ1,μ2,…,μn and coordinates of xi=μi and μi±fσi, where μi and σi are the mean and standard deviation of the random variable xi. f, which is an integer, ranges from 1 to 4 and is usually set to 1. Then, values of the performance function yx1,x2,…,xn containing yμ1,μ2,…,μi±fσi,…,σn and yμ1,μ2,…,σn are obtained using the slope analysis code. Thus, the training samples D=xi,yi are established, and the number of training samples is 2n × s + 1. s is the number of selected f.Step 3. If the training target y has a nonzero mean, to improve the efficiency of the training process of the model, y is adjusted by the mean μy of y(20)y^i=yi−μy.Step 4. Train the GP model by learning the training samples D=xi,y^ii=1,…,m. Then, we make predictions of performance function on the design point x∗ and obtain their predictions y¯x∗ according to equation (5). Recall that the training targets were centered, thus we must adjust the predictions by offset:(21)gx∗≈y¯x∗+μy=∑j=1mαkxj,x∗+μy.Step 5. Compute the reliability index of the kth iteration step β(k) using equation (17).Step 6. Compute the values of the new design point according to equation (21).Step 7. Check the convergence criterion for βk−βk−1≤ε (ε is 0.001 in present paper). To improve constantly the reconstructing precision at the important region where the failure probability is contributed significantly, the new design point and its value of the real performance function are taken as a new sample added into the old training sample if the convergence criteria are not satisfied. The number of training samples increases to 2n × s + 1 + k. Go to step 3 and repeat steps 3–6 until convergence is satisfied. If the convergence criterion is satisfied, go to step 8.Step 8. Calculate the probability of failure pf using equation (19).Figure 1
Flow chart of the GP-based RSM.To apply the GP-Based RSM to slope reliability analysis, a general program package is developed usingMATLAB.


## 3.1. Explicit Formation of Reliability Index Using GP

The performance function of a slope can be built as follows:(9)gX=Fs−1,where X denotes the random variables with mean μx and standard deviation σx. gX>0 indicates a stable slope, while gX<0 indicates a failed slope. gX=0 indicates that the slope is hovering between stable and unstable. Fs is the safety factor.Assuming that there is a pointx∗=x1∗,x2∗,…,xn∗T located on the limit state surface, which is called as the design point. Here, n is the dimension of x∗. gx is generally a nonlinear function. However, it can be linearized at x∗ by neglecting the second- or upper-order term:(10)ZQ=gx∗+∑i=1n∂gx∗∂xixi−xi∗.Assumingx is statistically uncorrelated, the mean value and standard deviation of ZQ are defined by(11)μZQ=gx∗+∑i=1n∂gx∗∂xiμxi−xi∗,σZQ=∑i=1n∂gx∗∂xi2σxi2.Reliability index can be calculated by(12)β=μZQσZQ=gx∗+∑i=1n∂gx∗/∂xiμxi−xi∗∑i=1n∂gx∗/∂xi2σxi2.According to the FORM, a tangent hyperplane is used to fit the limit state surface at the design point. Therefore, the most important step in the method is to obtain the design point. Reliability indexβ is defined as the distance from the origin to the design point. Hence, the coordinate of the design point x∗ is denoted as(13)xi∗=μxi+βσxicosθxi,i=1,2,…,n,where(14)cosθxi=−∂gx∗/∂xiσxi∑i=1n∂gx∗/∂xi2σxi2.According to equation (5), the performance function can be approximated by the GP model as(15)gx∗≈y¯x∗=∑j=1mαjkxj,x∗.The first-order partial derivatives of the approximate function can be given by(16)∂gx∗∂xi≈∂y¯x∗∂xi=∑j=1mαj∂kxj,x∗∂xi,where ∂kxj,x∗/∂xi=xi∗−xij/li2σf2exp−1/2l2·xj−x∗2.Substituting equations (15) and (16) into (17), the reliability index β is given as(17)β=∑j=1mαjkxj,x∗+μy+∑i=1n∑j=1mαj∂kxj,x∗/∂xiμxi−xi∗∑i=1n∑j=1mαj∂kxj,x∗/∂xi2σxi2.According to equations (15) and (16), (14) can be rewritten as(18)cosθxi=−∑j=1mαj∂kxj,x∗/∂xiσxi∑i=1n∑j=1mαj∂kxj,x∗/∂xi2σxi2.Substituting equations (17) and (18) into (13), one can obtain the coordinate of design point x∗ using GP approximation.The failure probability can be calculated as(19)pf≈Φ−β,where Φ⋅ denotes the cumulative distribution function of the standard normal variable.


## 3.2. Procedure of GP-Based RSM

The main difference between the GP-based RSM and the traditional RSM is that the former employs the GP to approximate the performance function and its first-order partial derivatives simultaneously. Moreover, from the viewpoint of machine learning theory, more knowledge means a better effect on the result of training or prediction of the machine learning model. Thus, the training samples, namely, knowledge set in the machine learning domain, is generated as a dynamic updating one to improve the performance of the GP model in the procedure of GP-based RSM.As shown in Figure1, reliability analysis using the GP-based RSM is described in detail as follows:Step 1. Assume initial values of the design point x∗=x1,x2,…,xn. Usually, mean values of the random variables can be selected as the coordinates of the design point x∗.Step 2. According to the experimental plan of the RSM developed by Bucher and Bourgund [29]; the training samples are generated according to the intersection of the axis xμ1,μ2,…,μn and coordinates of xi=μi and μi±fσi, where μi and σi are the mean and standard deviation of the random variable xi. f, which is an integer, ranges from 1 to 4 and is usually set to 1. Then, values of the performance function yx1,x2,…,xn containing yμ1,μ2,…,μi±fσi,…,σn and yμ1,μ2,…,σn are obtained using the slope analysis code. Thus, the training samples D=xi,yi are established, and the number of training samples is 2n × s + 1. s is the number of selected f.Step 3. If the training target y has a nonzero mean, to improve the efficiency of the training process of the model, y is adjusted by the mean μy of y(20)y^i=yi−μy.Step 4. Train the GP model by learning the training samples D=xi,y^ii=1,…,m. Then, we make predictions of performance function on the design point x∗ and obtain their predictions y¯x∗ according to equation (5). Recall that the training targets were centered, thus we must adjust the predictions by offset:(21)gx∗≈y¯x∗+μy=∑j=1mαkxj,x∗+μy.Step 5. Compute the reliability index of the kth iteration step β(k) using equation (17).Step 6. Compute the values of the new design point according to equation (21).Step 7. Check the convergence criterion for βk−βk−1≤ε (ε is 0.001 in present paper). To improve constantly the reconstructing precision at the important region where the failure probability is contributed significantly, the new design point and its value of the real performance function are taken as a new sample added into the old training sample if the convergence criteria are not satisfied. The number of training samples increases to 2n × s + 1 + k. Go to step 3 and repeat steps 3–6 until convergence is satisfied. If the convergence criterion is satisfied, go to step 8.Step 8. Calculate the probability of failure pf using equation (19).Figure 1
Flow chart of the GP-based RSM.To apply the GP-Based RSM to slope reliability analysis, a general program package is developed usingMATLAB.


## 4. Case Studies

### 4.1. Case 1: A Highly Nonlinear Performance Function

The performance function of the first example is defined as [30](22)gx1,x2=2+exp−x1210+x154−x2,where x1 and x2 are independent and obey standard normal distribution with zero mean and unit standard deviation.The training samples and predicted values of the performance function are listed in Table1. There are 5 training samples generated using the proposed method, where f is 1. The optimum hyperparameters of trained GP θ = (3.5197, 5.4426, 4.8151, −9.6465). The small prediction error from the GP indicates that the training is successful. If a vector of the basic random variables is input, the GP can generate an adequate accurate value of the corresponding performance function. Consequently, the GP can represent the true performance function well.Table 1
Training samples and the predicted values using GP for Example 1.
No.
x
1
x
2
g
x
GP approximation
1
−1
0
2.9064
2.9064
2
0
−1
4.0000
4.0000
3
1
0
2.9064
2.9063
4
0
1
2.0000
2.0000
5
0
0
3.0000
2.9999The GP approximation of the performance function based on the initial training samples is(23)gxi≈∑j=1mαikxi,xj+2.9625,where the covariance function kxi,xj denotes the column of matrix K and the values of vector α and matrix K based on the initial training samples are as follows:(24)α=−7.85601.8422−7.8568−1.719114.7597,K=0.95240.71900.75350.71900.89830.95240.71900.39090.76230.95240.71900.8983sym.0.95240.76230.9542.Table2 lists both calculation processes of the design point using the FORM and the proposed method. It clearly illustrates that the real performance and its first derivatives are both fitted very well by the GP model.Table 2
Calculation process of the design point using FORM and the proposed method for Example 1.
Method
Iterations
Design point
First derivative
Reliability index
x
1
∗
x
2
∗
∂
g
/
∂
x
1
∗
∂
g
/
∂
x
2
∗
β
FORM
I
0
0
0
−1
3
II
0
3
0
−1
3
III
0
3
—
—
3
GP-based RSM
I
0
0
0.0000
−1.0875
2.7026
II
0.0000
2.7026
−0.0001
−0.9999
3.0000
III
0.0002
3.0000
−0.0052
−0.9999
3.0000
IV
0.0157
3.0000
—
—
3.0000The failure probability and reliability index obtained by various methods in preexisting studies are listed in Table3. The exact reliability index is 2.9044 obtained by the direct MCS [30]. The reliability index yielded by the FORM, RSM, and GP-based RSM are very close to that by the MCS. The GP-based RSM needs 7 function evaluations, while the RSM and direct MCS require 19 and 106 function evaluations, respectively.Table 3
Results of Example 1 using different methods.
Method
MCS
FORM
RSM
GP-based RSM
Reliability index
2.9040
3.0000
3.0000
3.0000
Probability of failure (×10−3)
1.8420
1.3498
1.3500
1.3498
Number of function evaluation
106
2
19
7Figure2 shows the curves of approximated limit state function from the GP based on final training samples composed of 5 initial samples and 2 new samples and that of real limit state function. As can be seen from Figure 3, the GP can capture the whole shape of the real limit state function, and the design point is found quickly by the use of the GP-based RSM. This shows that the dynamic update of the training samples can significantly improve the efficiency of searching the design point in the proposed method.Figure 2
Compare the curve of limit state function from GP based on the final training samples with that from the real limit state function, where the boxes refer to the design points at iterative step in GP-based RSM: (b) enlarged figure; (a) in the range of −0.001 ≤x1 ≤ 0.004, 2.990 ≤ x2 ≤ 3.020.
(a)
(b)Figure 3
Geometry of slope for Example 2.


### 4.2. Case 2: A Two-Layered Slope

As illustrated in Figure3, the second example considers a two-layered slope which has been studied by Xu and Low [4]; Cho [14]; Hassan and Wolff [31]; Bhattacharya et al. [32]; and Chowdhury and Rao [33]. Property parameters of the soil related to the slope, including the friction angle and cohesion, are taken as random variables, as listed in Table 4. These random variables are assumed to obey normal distribution variables in the numerical analysis. Unit weight of the soil is assumed to be 19 kN/m3. By numerical analysis using SLIDE [24], the factor of safety Fx of circular slip surface was calculated using the Bishop and Spencer methods, respectively.Table 4
Properties of stochastic variables of Example 2.
Random variable
Mean
Standard deviation
c
1 (kN/m2)
38.31
7.662
c
2 (kN/m2)
23.94
4.788
φ
2 (deg)
12
1.2There are 7 training samples (f = 1) established for training the GP model as shown in Table 5. The safety factors of the other 10 samples created randomly were predicted using the trained GP model. The safety factors calculated by the Bishop code and Spencer method and values predicted by the GP are shown in Figure 4. The GP model can reflect the relationship between the random variables and the safety factor of the slope. Then, the performance functions in the explicit form using the GP model were established to approximate the implicit real performance functions.Table 5
The initial learning samples of Example 2.
Sample nos.
c
1 (kN/m2)
c
2 (kN/m2)
φ
2 (deg)
g
x
GP approximation
Bishop
Spencer
Bishop
Spencer
1
30.648
23.94
12
0.4576
0.4596
0.4576
0.4596
2
38.31
19.152
12
0.5412
0.5506
0.5412
0.5506
3
38.31
23.94
10.8
0.6047
0.6133
0.6051
0.6133
4
45.972
23.94
12
0.7795
0.7857
0.7795
0.7857
5
38.31
28.728
12
0.7442
0.7622
0.7442
0.7622
6
38.31
23.94
13.2
0.6792
0.6955
0.6796
0.6954
7
38.31
23.94
12
0.6443
0.6564
0.6433
0.6564Figure 4
Comparison between safety factors calculated by (a) Bishop method and GP prediction and (b) Spencer method and GP prediction.
(a)
(b)Table6 listed the results of the probabilistic analyses obtained from the previous literature for the same problem, while Tables 7 and 8 present the results obtained by the GP-based RSM. It can be observed from Tables 7 and 8, the result of reliability index by the early study shows that the exact solution is about 2.2 [33], and the reliability index estimation by the GP-based RSM with 7 initial samples is 2.249 from the Bishop method and 2.269 from the Spencer method, respectively. The results are quite in agreement with each other.Table 6
Results of probabilistic analysis from different methods for Example 2.
Method
Reliability index
Failure probability
Hassan and Wolff [31]: MVFOSM, Spencer method, noncircular slip surface
2.869
—
Bhattacharya et al. [32]: MVFOSM, Spencer method, noncircular slip surface by direct search
2.861
—
Crum [34]
2.23
1.29 × 10−2
Xu and Low [4]: response surface method, FEM
2.18
1.46 × 10−2
Xu and Low [4]: spreadsheet-based direct implementation of Spencer’s method with the search for the critical slip surface
2.20
1.39 × 10−2
Cho [14]: ANN-based response surface method with FORM
2.189
1.43 × 10−2
Chowdhury and Rao [33]: direct Monte Carlo (1 × 106 samples) using the Bishop method, circular slip surface
2.203
1.38 × 10−2
Chowdhury and Rao [33]: second-order HDMR, Bishop method
2.228
1.29 × 10−2
Chowdhury and Rao [33]: second-order HDMR, Spencer method
2.198
1.39 × 10−2
GP-based RSM (f = 1): Bishop method
2.249
1.23 × 10−2
GP-based RSM (f = 1): Spencer method
2.269
1.16 × 10−2Table 7
Results of probabilistic analysis by the Bishop method for Example 2.
Number of initial training samples
7 (f = 1)
13 (f = 1, 2)
19 (f = 1, 2, 3)
25 (f = 1, 2, 3, 4)
Reliability index
2.249
2.251
2.253
2.251
Design point
c1 (kN/m2)
21.198
21.070
21.055
21.074
c2 (kN/m2)
24.374
23.784
23.908
23.852
φ2 (deg)
11.701
12.080
12.073
12.097
F
s of design point
1.008
1.002
1.001
1.002Table 8
Results of probabilistic analysis by the Spencer method for Example 2.
Number of initial training samples
7 (f = 1)
13 (f = 1, 2)
19 (f = 1, 2, 3)
25 (f = 1, 2, 3, 4)
Reliability index
2.269
2.265
2.253
2.251
Design point
c1 (kN/m2)
20.962
20.958
21.093
21.116
c2 (kN/m2)
24.318
23.790
23.861
23.845
φ2 (deg)
11.852
12.007
12.199
12.208
F
s of design point
1.001
1.001
1.007
1.008Furthermore, Tables7 and 8 indicate that the change of the numbers of initial training samples slightly influences the accuracy of estimated values of probability of failure. It is worth noting that a larger number of initial training samples may mean higher robust and computation cost of the proposed method at the same time. The proposed method needs 12 function evaluations, while the ANN [14], MCS, first-order HDMR-based response surface method, and second-order HDMR-based response surface method [33] requires 25, 106, 13, and 61 number of original stability analysis, respectively. The larger the number of the function evaluations, the more the time needed to finish the computation is, the less efficient the used method is. Consequently, the GP-based RSM method is relatively efficient.The slip surfaces from the Bishop method and Spencer method at the mean point and at the design point are shown in Figures5 and 6, respectively. The slip surfaces at the mean point are very different from that at the design point. The values of the factor of safety Fs at the design point are both very close to 1; in another words, gX≈0, which mean the slope just reaches the limited equilibrium status.Figure 5
Example 2: results of the stability analysis by the Bishop method. Slip surface at (a) the mean point (Fs = 1.644) and (b) the design point (Fs = 1.008).
(a)
(b)Figure 6
Example 2: results of the stability analysis by the Spencer method. Slip surface at (a) the mean point (Fs = 1.656) and (b) the design point (Fs = 1.001).
(a)
(b)


### 4.3. Case 3: A Three-Layered Slope

A complex slope with three different soil layers is studied. The cross section of the nonhomogeneous slope with a height of 10 m and an inclination of 2 : 1 is presented in Figure7. The unit weight γ = 19.5 kN·m−3. The corresponding property parameters of the soil related to the slope, including the cohesion c and friction angle φ, are considered as random variables. These random variables obey the normal distribution. The values of the mean and standard deviation of the random variables are listed in Table 9.Figure 7
Geometry of the slope of Example 3.Table 9
Properties of stochastic variables of Example 3.
Soil type
c (kPa)
φ (deg)
Mean value
Standard deviation
Mean value
Standard deviation
1
0
0
38
4
2
5.3
0.7
23
3
3
7.2
0.2
20
3The calculations of the safety factor of the slope are conducted using the Bishop method and the Spencer method by the SLIDE code [35]. 11 initial training samples (f = 1) are listed in Table 10 while the safety factor of the circular slip surface was calculated using the Bishop and Spencer methods, respectively.Table 10
The initial learning samples of Example 3.
Sample nos.
φ
1 (deg)
c
2 (kN/m2)
φ
2 (deg)
c
3 (kN/m2)
φ
3 (deg)
g
x
Bishop
Spencer
1
38
5.3
23
7.2
20
0.4052
0.3747
2
34
5.3
23
7.2
20
0.3493
0.3251
3
38
4.6
23
7.2
20
0.4009
0.3705
4
38
5.3
20
7.2
20
0.3744
0.3491
5
38
5.3
23
7
20
0.3994
0.3691
6
38
5.3
23
7.2
17
0.3009
0.2685
7
42
5.3
23
7.2
20
0.4507
0.4192
8
38
6
23
7.2
20
0.4096
0.3788
9
38
5.3
26
7.2
20
0.4343
0.4008
10
38
5.3
23
7.4
20
0.4111
0.3803
11
38
5.3
23
7.2
23
0.5033
0.4757The results of the slope reliability analysis using different methods are listed in Table11. A sampling size of 106 is considered in the direct MCS by numerical analysis using SLIDE [35] to estimate the failure probability Pf. The failure probability estimation by the GP-based RSM is in good agreement with that of the MCS. The proposed method needs 12 function evaluations, while the number of function evaluations for the SVM-based RSM [20] and MCS are 19 and 106, respectively. This indicates that the proposed GP-based RSM is much more efficient to achieve reasonable accuracy for slope reliability analysis than the SVM-based RSM [20].Table 11
Results comparison of Example 3.
Method
SLIDE [35] MCS (1 × 106 samples)
SVM-based RSM [20]
GP-based RSM
Bishop
Spencer
Bishop
Spencer
Bishop
Spencer
Reliability indexβ
3.542
3.287
3.4293
3.0481
3.4853
3.2118
Number of function evaluation
106
106
19
19
12
12In addition, the slip surfaces from the Bishop method and Spencer method at the design point obtained using the proposed method are both presented in Figure8.Figure 8
Example 3: results of the stability analysis. Slip surface at the design point by (a) the Bishop method (Fs = 0.992) and (b) the Spencer method (Fs = 0.997).
(a)
(b)


### 4.4. Example 4: Congress St. Cut Model

In this example, the confessed Congress St. Cut model studied by Chowdhury and Xu [36] is considered. The soil profile is presented in Figure 9. Two sets of circular slip surfaces are considered. The first set consists of potential failure surface tangential to the lower boundary of the Clay 2 layer (failure mode I), while the second considers slip surfaces tangential to the lower boundary of Clay 3 (failure mode II). The property parameters of soil are listed in Table 12 and are assumed to obey the Gaussian distribution in the numerical analysis. The Bishop method is used to estimate the safety factor from the numerical analysis using SLIDE [35]. For evaluating the failure probability Pf, f = 1 is selected, which results in 7 initial training samples (Table 13). Without any iteration, totally 7 actual stability analyses using the Bishop method are needed for the GP-based RSM. It can be observed from Table 14 that the failure probability estimated by the GP-based RSM agrees well with that reported in the literature [33]. However, the total numbers of slope stability analysis for the proposed method is 12, while that for the MCS, first-order HDMR, and second-order HDMR are 5000, 13, and 61, respectively [33] This indicates that the proposed method is very efficient to achieve reasonable accuracy. The slip surface of failure mode I and failure mode II at the design point obtained using the proposed method are both presented in Figure 10.Figure 9
Geometry of slope of Example 4.Table 12
Material properties of the soil of Example 4.
Material
c (kPa)
φ (deg)
r (kN/m3)
Mean
Standard deviation
Mean
Standard deviation
Clay 1
68.1
6.6
0
—
21
Clay 2
39.3
1.4
0
—
22
Clay 3
50.8
1.5
0
—
22
Sand
0
—
0
—
21
Hard clay
200
—
35
—
—Table 13
The initial learning samples of Example 4.
Sample nos.
c
1 (kN/m2)
c
2 (kN/m2)
c
3 (kN/m2)
g
x
Failure mode I
Failure mode II
1
68.1
39.3
50.8
0.1083
0.0578
2
61.5
39.3
50.8
0.0798
0.0398
3
68.1
37.9
50.8
0.0793
0.0505
4
68.1
39.3
49.3
0.1083
0.0382
5
74.7
39.3
50.8
0.1368
0.0758
6
68.1
40.7
50.8
0.1373
0.0652
7
68.1
39.3
52.3
0.1083
0.0775Table 14
Results of probabilistic analyses using different methods for Example 4.
Method
Failure mode
Probability of failure
Chowdhury and Rao [33]: direct Monte Carlo (5000 samples) using the Bishop method, circular slip surface
I
0.0037
II
0.0175
Chowdhury and Rao [33]: first-order HDMR, Bishop method
I
0.0041
II
0.0171
Chowdhury and Rao [33]: second-order HDMR, Bishop method
I
0.0044
II
0.0166
Chowdhury and Xu [36]
I
0.0048
II
0.0131
GP-based RSM: Bishop method
I
0.0039
II
0.0185Figure 10
Example 4: results of the stability analysis by Bishop method. Slip surface at the design point: (a) failure model I,Fs = 1.000 and (b) failure model II, Fs = 1.000.
(a)
(b)


## 4.1. Case 1: A Highly Nonlinear Performance Function

The performance function of the first example is defined as [30](22)gx1,x2=2+exp−x1210+x154−x2,where x1 and x2 are independent and obey standard normal distribution with zero mean and unit standard deviation.The training samples and predicted values of the performance function are listed in Table1. There are 5 training samples generated using the proposed method, where f is 1. The optimum hyperparameters of trained GP θ = (3.5197, 5.4426, 4.8151, −9.6465). The small prediction error from the GP indicates that the training is successful. If a vector of the basic random variables is input, the GP can generate an adequate accurate value of the corresponding performance function. Consequently, the GP can represent the true performance function well.Table 1
Training samples and the predicted values using GP for Example 1.
No.
x
1
x
2
g
x
GP approximation
1
−1
0
2.9064
2.9064
2
0
−1
4.0000
4.0000
3
1
0
2.9064
2.9063
4
0
1
2.0000
2.0000
5
0
0
3.0000
2.9999The GP approximation of the performance function based on the initial training samples is(23)gxi≈∑j=1mαikxi,xj+2.9625,where the covariance function kxi,xj denotes the column of matrix K and the values of vector α and matrix K based on the initial training samples are as follows:(24)α=−7.85601.8422−7.8568−1.719114.7597,K=0.95240.71900.75350.71900.89830.95240.71900.39090.76230.95240.71900.8983sym.0.95240.76230.9542.Table2 lists both calculation processes of the design point using the FORM and the proposed method. It clearly illustrates that the real performance and its first derivatives are both fitted very well by the GP model.Table 2
Calculation process of the design point using FORM and the proposed method for Example 1.
Method
Iterations
Design point
First derivative
Reliability index
x
1
∗
x
2
∗
∂
g
/
∂
x
1
∗
∂
g
/
∂
x
2
∗
β
FORM
I
0
0
0
−1
3
II
0
3
0
−1
3
III
0
3
—
—
3
GP-based RSM
I
0
0
0.0000
−1.0875
2.7026
II
0.0000
2.7026
−0.0001
−0.9999
3.0000
III
0.0002
3.0000
−0.0052
−0.9999
3.0000
IV
0.0157
3.0000
—
—
3.0000The failure probability and reliability index obtained by various methods in preexisting studies are listed in Table3. The exact reliability index is 2.9044 obtained by the direct MCS [30]. The reliability index yielded by the FORM, RSM, and GP-based RSM are very close to that by the MCS. The GP-based RSM needs 7 function evaluations, while the RSM and direct MCS require 19 and 106 function evaluations, respectively.Table 3
Results of Example 1 using different methods.
Method
MCS
FORM
RSM
GP-based RSM
Reliability index
2.9040
3.0000
3.0000
3.0000
Probability of failure (×10−3)
1.8420
1.3498
1.3500
1.3498
Number of function evaluation
106
2
19
7Figure2 shows the curves of approximated limit state function from the GP based on final training samples composed of 5 initial samples and 2 new samples and that of real limit state function. As can be seen from Figure 3, the GP can capture the whole shape of the real limit state function, and the design point is found quickly by the use of the GP-based RSM. This shows that the dynamic update of the training samples can significantly improve the efficiency of searching the design point in the proposed method.Figure 2
Compare the curve of limit state function from GP based on the final training samples with that from the real limit state function, where the boxes refer to the design points at iterative step in GP-based RSM: (b) enlarged figure; (a) in the range of −0.001 ≤x1 ≤ 0.004, 2.990 ≤ x2 ≤ 3.020.
(a)
(b)Figure 3
Geometry of slope for Example 2.


## 4.2. Case 2: A Two-Layered Slope

As illustrated in Figure3, the second example considers a two-layered slope which has been studied by Xu and Low [4]; Cho [14]; Hassan and Wolff [31]; Bhattacharya et al. [32]; and Chowdhury and Rao [33]. Property parameters of the soil related to the slope, including the friction angle and cohesion, are taken as random variables, as listed in Table 4. These random variables are assumed to obey normal distribution variables in the numerical analysis. Unit weight of the soil is assumed to be 19 kN/m3. By numerical analysis using SLIDE [24], the factor of safety Fx of circular slip surface was calculated using the Bishop and Spencer methods, respectively.Table 4
Properties of stochastic variables of Example 2.
Random variable
Mean
Standard deviation
c
1 (kN/m2)
38.31
7.662
c
2 (kN/m2)
23.94
4.788
φ
2 (deg)
12
1.2There are 7 training samples (f = 1) established for training the GP model as shown in Table 5. The safety factors of the other 10 samples created randomly were predicted using the trained GP model. The safety factors calculated by the Bishop code and Spencer method and values predicted by the GP are shown in Figure 4. The GP model can reflect the relationship between the random variables and the safety factor of the slope. Then, the performance functions in the explicit form using the GP model were established to approximate the implicit real performance functions.Table 5
The initial learning samples of Example 2.
Sample nos.
c
1 (kN/m2)
c
2 (kN/m2)
φ
2 (deg)
g
x
GP approximation
Bishop
Spencer
Bishop
Spencer
1
30.648
23.94
12
0.4576
0.4596
0.4576
0.4596
2
38.31
19.152
12
0.5412
0.5506
0.5412
0.5506
3
38.31
23.94
10.8
0.6047
0.6133
0.6051
0.6133
4
45.972
23.94
12
0.7795
0.7857
0.7795
0.7857
5
38.31
28.728
12
0.7442
0.7622
0.7442
0.7622
6
38.31
23.94
13.2
0.6792
0.6955
0.6796
0.6954
7
38.31
23.94
12
0.6443
0.6564
0.6433
0.6564Figure 4
Comparison between safety factors calculated by (a) Bishop method and GP prediction and (b) Spencer method and GP prediction.
(a)
(b)Table6 listed the results of the probabilistic analyses obtained from the previous literature for the same problem, while Tables 7 and 8 present the results obtained by the GP-based RSM. It can be observed from Tables 7 and 8, the result of reliability index by the early study shows that the exact solution is about 2.2 [33], and the reliability index estimation by the GP-based RSM with 7 initial samples is 2.249 from the Bishop method and 2.269 from the Spencer method, respectively. The results are quite in agreement with each other.Table 6
Results of probabilistic analysis from different methods for Example 2.
Method
Reliability index
Failure probability
Hassan and Wolff [31]: MVFOSM, Spencer method, noncircular slip surface
2.869
—
Bhattacharya et al. [32]: MVFOSM, Spencer method, noncircular slip surface by direct search
2.861
—
Crum [34]
2.23
1.29 × 10−2
Xu and Low [4]: response surface method, FEM
2.18
1.46 × 10−2
Xu and Low [4]: spreadsheet-based direct implementation of Spencer’s method with the search for the critical slip surface
2.20
1.39 × 10−2
Cho [14]: ANN-based response surface method with FORM
2.189
1.43 × 10−2
Chowdhury and Rao [33]: direct Monte Carlo (1 × 106 samples) using the Bishop method, circular slip surface
2.203
1.38 × 10−2
Chowdhury and Rao [33]: second-order HDMR, Bishop method
2.228
1.29 × 10−2
Chowdhury and Rao [33]: second-order HDMR, Spencer method
2.198
1.39 × 10−2
GP-based RSM (f = 1): Bishop method
2.249
1.23 × 10−2
GP-based RSM (f = 1): Spencer method
2.269
1.16 × 10−2Table 7
Results of probabilistic analysis by the Bishop method for Example 2.
Number of initial training samples
7 (f = 1)
13 (f = 1, 2)
19 (f = 1, 2, 3)
25 (f = 1, 2, 3, 4)
Reliability index
2.249
2.251
2.253
2.251
Design point
c1 (kN/m2)
21.198
21.070
21.055
21.074
c2 (kN/m2)
24.374
23.784
23.908
23.852
φ2 (deg)
11.701
12.080
12.073
12.097
F
s of design point
1.008
1.002
1.001
1.002Table 8
Results of probabilistic analysis by the Spencer method for Example 2.
Number of initial training samples
7 (f = 1)
13 (f = 1, 2)
19 (f = 1, 2, 3)
25 (f = 1, 2, 3, 4)
Reliability index
2.269
2.265
2.253
2.251
Design point
c1 (kN/m2)
20.962
20.958
21.093
21.116
c2 (kN/m2)
24.318
23.790
23.861
23.845
φ2 (deg)
11.852
12.007
12.199
12.208
F
s of design point
1.001
1.001
1.007
1.008Furthermore, Tables7 and 8 indicate that the change of the numbers of initial training samples slightly influences the accuracy of estimated values of probability of failure. It is worth noting that a larger number of initial training samples may mean higher robust and computation cost of the proposed method at the same time. The proposed method needs 12 function evaluations, while the ANN [14], MCS, first-order HDMR-based response surface method, and second-order HDMR-based response surface method [33] requires 25, 106, 13, and 61 number of original stability analysis, respectively. The larger the number of the function evaluations, the more the time needed to finish the computation is, the less efficient the used method is. Consequently, the GP-based RSM method is relatively efficient.The slip surfaces from the Bishop method and Spencer method at the mean point and at the design point are shown in Figures5 and 6, respectively. The slip surfaces at the mean point are very different from that at the design point. The values of the factor of safety Fs at the design point are both very close to 1; in another words, gX≈0, which mean the slope just reaches the limited equilibrium status.Figure 5
Example 2: results of the stability analysis by the Bishop method. Slip surface at (a) the mean point (Fs = 1.644) and (b) the design point (Fs = 1.008).
(a)
(b)Figure 6
Example 2: results of the stability analysis by the Spencer method. Slip surface at (a) the mean point (Fs = 1.656) and (b) the design point (Fs = 1.001).
(a)
(b)


## 4.3. Case 3: A Three-Layered Slope

A complex slope with three different soil layers is studied. The cross section of the nonhomogeneous slope with a height of 10 m and an inclination of 2 : 1 is presented in Figure7. The unit weight γ = 19.5 kN·m−3. The corresponding property parameters of the soil related to the slope, including the cohesion c and friction angle φ, are considered as random variables. These random variables obey the normal distribution. The values of the mean and standard deviation of the random variables are listed in Table 9.Figure 7
Geometry of the slope of Example 3.Table 9
Properties of stochastic variables of Example 3.
Soil type
c (kPa)
φ (deg)
Mean value
Standard deviation
Mean value
Standard deviation
1
0
0
38
4
2
5.3
0.7
23
3
3
7.2
0.2
20
3The calculations of the safety factor of the slope are conducted using the Bishop method and the Spencer method by the SLIDE code [35]. 11 initial training samples (f = 1) are listed in Table 10 while the safety factor of the circular slip surface was calculated using the Bishop and Spencer methods, respectively.Table 10
The initial learning samples of Example 3.
Sample nos.
φ
1 (deg)
c
2 (kN/m2)
φ
2 (deg)
c
3 (kN/m2)
φ
3 (deg)
g
x
Bishop
Spencer
1
38
5.3
23
7.2
20
0.4052
0.3747
2
34
5.3
23
7.2
20
0.3493
0.3251
3
38
4.6
23
7.2
20
0.4009
0.3705
4
38
5.3
20
7.2
20
0.3744
0.3491
5
38
5.3
23
7
20
0.3994
0.3691
6
38
5.3
23
7.2
17
0.3009
0.2685
7
42
5.3
23
7.2
20
0.4507
0.4192
8
38
6
23
7.2
20
0.4096
0.3788
9
38
5.3
26
7.2
20
0.4343
0.4008
10
38
5.3
23
7.4
20
0.4111
0.3803
11
38
5.3
23
7.2
23
0.5033
0.4757The results of the slope reliability analysis using different methods are listed in Table11. A sampling size of 106 is considered in the direct MCS by numerical analysis using SLIDE [35] to estimate the failure probability Pf. The failure probability estimation by the GP-based RSM is in good agreement with that of the MCS. The proposed method needs 12 function evaluations, while the number of function evaluations for the SVM-based RSM [20] and MCS are 19 and 106, respectively. This indicates that the proposed GP-based RSM is much more efficient to achieve reasonable accuracy for slope reliability analysis than the SVM-based RSM [20].Table 11
Results comparison of Example 3.
Method
SLIDE [35] MCS (1 × 106 samples)
SVM-based RSM [20]
GP-based RSM
Bishop
Spencer
Bishop
Spencer
Bishop
Spencer
Reliability indexβ
3.542
3.287
3.4293
3.0481
3.4853
3.2118
Number of function evaluation
106
106
19
19
12
12In addition, the slip surfaces from the Bishop method and Spencer method at the design point obtained using the proposed method are both presented in Figure8.Figure 8
Example 3: results of the stability analysis. Slip surface at the design point by (a) the Bishop method (Fs = 0.992) and (b) the Spencer method (Fs = 0.997).
(a)
(b)


## 4.4. Example 4: Congress St. Cut Model

In this example, the confessed Congress St. Cut model studied by Chowdhury and Xu [36] is considered. The soil profile is presented in Figure 9. Two sets of circular slip surfaces are considered. The first set consists of potential failure surface tangential to the lower boundary of the Clay 2 layer (failure mode I), while the second considers slip surfaces tangential to the lower boundary of Clay 3 (failure mode II). The property parameters of soil are listed in Table 12 and are assumed to obey the Gaussian distribution in the numerical analysis. The Bishop method is used to estimate the safety factor from the numerical analysis using SLIDE [35]. For evaluating the failure probability Pf, f = 1 is selected, which results in 7 initial training samples (Table 13). Without any iteration, totally 7 actual stability analyses using the Bishop method are needed for the GP-based RSM. It can be observed from Table 14 that the failure probability estimated by the GP-based RSM agrees well with that reported in the literature [33]. However, the total numbers of slope stability analysis for the proposed method is 12, while that for the MCS, first-order HDMR, and second-order HDMR are 5000, 13, and 61, respectively [33] This indicates that the proposed method is very efficient to achieve reasonable accuracy. The slip surface of failure mode I and failure mode II at the design point obtained using the proposed method are both presented in Figure 10.Figure 9
Geometry of slope of Example 4.Table 12
Material properties of the soil of Example 4.
Material
c (kPa)
φ (deg)
r (kN/m3)
Mean
Standard deviation
Mean
Standard deviation
Clay 1
68.1
6.6
0
—
21
Clay 2
39.3
1.4
0
—
22
Clay 3
50.8
1.5
0
—
22
Sand
0
—
0
—
21
Hard clay
200
—
35
—
—Table 13
The initial learning samples of Example 4.
Sample nos.
c
1 (kN/m2)
c
2 (kN/m2)
c
3 (kN/m2)
g
x
Failure mode I
Failure mode II
1
68.1
39.3
50.8
0.1083
0.0578
2
61.5
39.3
50.8
0.0798
0.0398
3
68.1
37.9
50.8
0.0793
0.0505
4
68.1
39.3
49.3
0.1083
0.0382
5
74.7
39.3
50.8
0.1368
0.0758
6
68.1
40.7
50.8
0.1373
0.0652
7
68.1
39.3
52.3
0.1083
0.0775Table 14
Results of probabilistic analyses using different methods for Example 4.
Method
Failure mode
Probability of failure
Chowdhury and Rao [33]: direct Monte Carlo (5000 samples) using the Bishop method, circular slip surface
I
0.0037
II
0.0175
Chowdhury and Rao [33]: first-order HDMR, Bishop method
I
0.0041
II
0.0171
Chowdhury and Rao [33]: second-order HDMR, Bishop method
I
0.0044
II
0.0166
Chowdhury and Xu [36]
I
0.0048
II
0.0131
GP-based RSM: Bishop method
I
0.0039
II
0.0185Figure 10
Example 4: results of the stability analysis by Bishop method. Slip surface at the design point: (a) failure model I,Fs = 1.000 and (b) failure model II, Fs = 1.000.
(a)
(b)


## 5. Conclusions

A GP-based RSM was developed to predict failure probability of slope. In this method, the training samples for establishing the GP model are first generated by the design method of the classical RSM method. The limit state function and its first-order partial derivative are then approximated by the trained GP model. Finally, failure probability is estimated using the FORM. In addition, the viewpoint of more knowledge means the better effect on the result of training or prediction of the machine learning model, and the use of the renewed continuously training samples may improve the performance of the GP model for approximating the limit state function around the design point. Thus, the GP significantly reduces the number of required training samples. Otherwise, GP also shows good performance to approximate the limit state function and then provides accurate estimation of the failure probability when connected with the FORM. Four numerical examples including both slope and nonslope problems illustrate the flexibility and efficiency of the proposed method. Compared with the traditional RSM, the GP-based RSM is much more efficient to achieve reasonable accuracy for slope reliability analysis. Moreover, the present method can directly take advantage of existing slope software without modification and thus are convenient to be used for practitioner engineers. However, it should be noted that the proposed algorithm is intended as a possible complement rather than a replacement for existing classical methods.



---

*Source: 9185756-2019-03-27.xml*