Abstract

The estimation of the cross-correlation of shear strength parameters (i.e., cohesion and internal friction angle) and the subsequent determination of the probability of failure have long been challenges in slope reliability analysis. Here, a copula-based approach is proposed to calculate the probability of failure by integrating the copula-based joint probability density function (PDF) on the slope failure domain delimited with the -line. Here, copulas are used to construct the joint PDF of shear strength parameters with specific marginal distributions and correlation structure. In the paper a failure (limit state) function approach is applied to investigate a system characterized by a homogeneous slope. The results show that the values obtained by using the failure function approach are similar to those calculated by means of conventional methods, such as the first-order reliability method (FORM) and Monte Carlo simulations (MC). In addition, an entropy weight (EW) copula is proposed to address the discrepancies of the results calculated by different copulas to avoid over- or underestimating the slope reliability.

1. Introduction

Numerous problems regarding uncertainty in geotechnical engineering exist currently, such as soil and rock properties, environmental conditions, and theoretical models [13]. The application of reliability analysis for addressing uncertainty in geotechnical design is rather common [46]. To determine the reliability of a slope with high accuracy, it is imperative to investigate the relevance of parameters and establish their joint cumulative distribution function (CDF) or joint probability density function (PDF). However, establishing an accurate joint CDF or PDF remains a practical challenge in geotechnical engineering due to the limited available datasets [79]. Copulas may be used to address this problem [10, 11]. The use of copulas has the distinct advantage of enabling the construction of a joint CDF of random variables with any form of marginal distribution, overcoming the limitation of conventional approaches [12, 13] in which nonnormal variables must be converted into normal variables [9, 14].

In recent years, copulas have been applied to different disciplines, including finance [10, 15], biomedicine [16, 17], hydrology [1820], and coastal engineering [21, 22]. In geotechnical engineering, Li et al. [23] used copulas to simulate the bivariate probability distribution of two curve-fitting parameters underlying the load-displacement models of piles. Tang et al. [9, 14] calculated the failure probabilities for an infinite slope and a retaining wall using the direct integration method and investigated the impact of different types of copulas on reliability. Marchant et al. [24] introduced a more general multivariate function for the spatial prediction of soil properties based on copulas. Wu [25, 26] determined the failure probability of a slope using a copula-based sampling method. Huang et al. [27] used a copula-based method to estimate the shear strength parameters of rock mass. Motamedi and Liang [28] proposed a probabilistic methodology for landslide hazard assessment at regional scale based on the copula modelling technique considering the possible dependence between the hazard elements.

Although copulas have been well developed and widely applied to determine the relevance of material properties and to establish the CDF or PDF, it is difficult to determine the performance function of a slope system in the critical state. Analysing the reliability of a slope remains a significant challenge in practice.

With the help of the -line [29], we propose a reliability analysis method based on copulas, which can be applied to any given homogeneous slopes. The remainder of this paper is organized as follows. Section 2 briefly introduces the copula functions and presents the procedure used in the copula-based slope reliability analysis. In Section 3, after finishing the detailed introductions of -line, we present an approach to identify the failure domain of a slope of any shape. The range of shear strength parameters, which can represent the slope failure domain, is also deduced. Section 4 explains how to calculate the failure probability by directly integrating the joint PDF of the considered variables on the failure domain identified by the -line. Section 5 presents a demonstrative application with evaluation of the failure probability for an example slope. Finally, the results calculated from different copulas are compared and discussed for an entropy weight (EW) copula in Section 6.

2. Copulas for the Shear Strength Parameters

2.1. Copula Function

A copula function is a multivariate distribution function with uniform marginal distributions which are all over . For an -dimensional uniform random vector , a copula is expressed as follows [30]:where and are the th uniform random variables on the interval and the corresponding values and represents the probability value.

Using Sklar’s theorem [31], the joint CDF of an -dimensional random variable can be established via the copula function aswhere is a copula function, is the marginal distribution function of the th random variable and are the copula parameters. Copulas offer the possibility of separately modelling the marginal distributions and the dependence structure among the random variables.

Using a one parameter copula function, the joint CDF of the shear strength parameters and , corresponding to cohesion and internal friction angle , can be expressed as follows:in which the parameter characterizes the dependence structure between the variables and .

The corresponding joint PDF iswhere and are the PDFs of and and is the copula density function given by

2.2. Slope Reliability Assessment Using Copulas

The establishment of the joint CDF or PDF in terms of copulas for the implementation of the slope reliability analysis involves the following 4 steps.

(1) Estimation of the Marginal Distributions of Each Random Variable. As for fitting copulas for shear strength, the first-line procedure is to determine the appropriate marginal distribution of each random variable. In various previous studies the normal, lognormal, and gamma distributions were selected as the best-fitting marginal distributions for the shear strength parameters ( and ) [25, 32, 33]. The lognormal distribution is defined only for nonnegative values [32]. Kolmogorov-Smirnov (K-S) test commonly used to check whether a set of data follows a certain distribution is adopted to examine marginal distribution [27]. The key procedure of K-S test is computing the test statistic , which can be expressed aswhere is the established marginal distribution functions and is the cumulative frequency distribution; it is defined asin which is the indicator function, equal to 1 if and equal to 0 otherwise. will be accepted if ,where is the critical value of a given significance level . can be obtained by using the established table [34] when and have been determined, and is the number of observed datasets.

(2) Copula Parameter Estimation. Two-dimensional copulas commonly used in geotechnical engineering belong to two different classes. The elliptical class includes the Gaussian and copulas. The Archimedean class includes the Frank, Gumbel, and Clayton copulas.

Most analyses performed on the shear strength parameters report a negative correlation coefficient between and , with value between −0.24 and −0.70. However, a positive correlation was also found in several cases [9, 25]. The shear strength parameters of rocks and soils typically have symmetric structures [9] and asymmetric structures [35, 36]. Considering the correlation characteristics and symmetry between shear strength parameters, the capability of the Gaussian, Frank, Clayton, Gumbel, and Plackett copulas is tested to model the dependence structure of the data [9, 25]. The considered copulas are symmetric copulas, except for the Clayton and Gumbel copulas. In addition, negative correlations may not be allowed in some of those copulas; nonetheless, a simple conversion which negates the values of one variable can achieve a positive value for the correlation.

Copula parameters can be obtained by using maximum likelihood estimation (MLE) [30], the Pearson linear correlation coefficient , and the Kendall or Spearman rank correlation coefficient [11]. In this paper, the Kendall rank correlation coefficient is chosen to estimate the copula parameter, . The parameter is independent of the marginal distribution. For the elliptical class of copulas, the relationship between and is given byFor the Archimedean class of copulas, Kendall’s is given bywhere is the generator function. Kendall’s may be computed by where is the total number of samples, , and is the sign function, which will be 1 for or in all other cases.

Then, the corresponding value is determined by using (8) or (9) depending on the copula class. The parameter of the Plackett copula was determined with the maximum likelihood method [30].

(3) Identification of the Best-Fitting Copula. There are various methods to identify the best-fitting copula model, among these candidate copulas, such as the Akaike information criterion (AIC) [37], root mean square error (RMSE), bias [38], and Cramér-von Mises statistical method [39, 40]. The AIC is the method most commonly used in engineering practice, and it can be defined as follows:In this study, AICc [41] takes a correction into consideration for finite sample sizes, which can be determined by the AIC:where is the number of estimated parameters and denotes the number of observations (for a single parameter, ). The copula with the smallest AICc value is preferred. The Cramér-von Mises statistics, , can be expressed as the sum of the squared differences between the true and empirical copulas for any considered point [39, 40]:where is the copula with the estimated value of the copula parameter, , and is the empirical copula given byin which and are empirical CDFs, is the indicator function, if , , and otherwise. A parametric bootstrapping approach is used to determine values [40]. The copula with the minimum and maximum can be considered as optimal. The RMSE and bias can be determined by the following equations [38]:where is the th observed value and is the th estimated value. The best-fitting model is the one that has the smallest RMSE or bias.

(4) Slope Reliability Analysis. Following three successive steps, the joint CDF or PDF of the random variable can be obtained. Subsequently, the probability of failure can be calculated by integrating the joint PDF established by copulas on the slope failure domain; this approach is called the copula-based direct integration method (CDIM). The determination of the slope failure domain and the integral solution of the probability of failure will be discussed in Sections 3 and 4, respectively.

3. -Line and Slope Failure Domain

3.1. -Line

According to Klar et al. [29], the factor of safety of a slope system can be expressed as follows based on the Mohr-Coulomb criterion:where and are the actual shear strength parameters of the soil slope, and are the mobilized parameters required for the limit equilibrium state, denotes the shear strength and denotes the mobilized shear strength with a factor of safety , and is normal stress. Klar et al. [29] proposed the concept of the -line based on (17) to determine the factor of safety geometrically. The -line (Figure 1) represents the combinations of and at which the slope will exist in limit equilibrium (). The combinations of and above the -line imply a stable condition (), whereas the combinations falling below the -line (namely, in the failure domain symbolized by ) are impossible and represent nonphysical conditions since the system cannot mobilize more than its capacity. The geometrical expression of can be written asin which denotes the distance between point and the origin and the distance between the origin and point is designated by . is an arbitrary point in space -, and is the intersection of the -line and .

In the present paper, the simplified Bishop method is used to calculate the factor of safety. Using the software SLOPE/W, hundreds of runs are made to search the critical state () of a slope with given geometry and unit weight, and then the corresponding and values in the limit equilibrium are obtained. The -line is fitted by the combinations of and in the limit equilibrium state.

For any given shape of homogeneous slope, the - values in the limit equilibrium state can always be used to determine the -line in a unique manner when the slope geometry characteristic and unit weight are given. In other words, the -line can be regarded as an inherent characteristic of the slope used to characterize the slope’s limit equilibrium state. The -line makes it convenient to conduct the reliability analysis of the slope.

3.2. -Line Fitting Shape

The -line determines the failure domain of the slope, and its shape also directly affects the results of the slope reliability analysis. However, there are few studies regarding the shape of the -line. In this paper, several slopes with different slope heights (), slope angles (), and unit weights () are used to calculate the and values in the limit equilibrium, as shown in Table 1. Figure 2(a) and Table 1 clearly indicate that the quadratic polynomial is effective for fitting the -line for general homogeneous slopes, and its coefficient of determination () is generally in the range of 0.92 to 0.99. Considering the value of the quadratic polynomial as an objective reference, the slope angle is more sensitive than the unit weight, for a same value. In contrast, the sensitivity of the slope height is relatively low. Quadratic polynomial fitting has a relatively poor performance when considering high and steep slopes, as shown in Figure 2; in contrast, higher-degree polynomials exhibited a better fitting effect.

3.3. Failure Domain Represented by the Shear Strength Parameters

Figure 1 shows the -line and failure domain of a specific slope, where  m, = 45°, and = 20 kN/m3. The -line can be expressed in terms of quadratic expressions, . Similarly, the failure domain can be expressed as . For , the corresponding value reaches a maximum (i.e., ). Similarly, the failure domain is equivalent to ,   and , if and only if , the value is the maximum, and the maximum value of the cohesive force can be expressed as . The range of shear strength parameters in the failure domain must be determined to compute the probability of failure.

4. Application of the Direct Integration Method

The probability of failure can be accurately determined through the direct integration approach [9], which can mitigate the errors inherent in many reliability methods, such as the first-order reliability method (FORM), second-order reliability method (SORM), and Monte Carlo simulation (MCS). In this study, the direct integration approach is adopted to integrate the copula-based - joint PDF for the -line failure domain (CDIM). Thus, the probability of failure of a slope with any complex shape can be calculated by CDIM.

Based on the limit state, the performance function of slope for the -line can be expressed as follows:in which , , and are quadratic polynomial coefficients, is the cohesion,and is the internal friction angle. The limit state equation of the -line corresponds to , and the failure domain corresponds to    (); thus, the probability of failure can obtained as

By substituting the copula-based joint probability density function given by (4) into (20), the expression of can be transformed as

Because of the complexity of the copula function, the above double-integral calculation is often difficult to perform, (21) can be converted into a single integral by defining first derivative of the copula [9], and here and denote and , respectively, such that

By substituting the range of shear strength parameters in the -line failure domain into (22), the final expression of solving is

5. Illustrative Example

5.1. Slope Model and Its -Line Failure Domain

A hypothetical homogeneous slope shown in Figure 3 is employed to investigate the probability of failure using the CDIM, with the slope height () of 30.0 m, the slope angle () of 30.0°, and the unit weight () of 18 kN/m3. Figure 3 shows the -line for the slope model and its expression of the quadratic polynomial fitting ( = 0.9788). From the analysis in Section 3.3, the parameter range of cohesion characterizing the slope failure domain is and = 139.8 kPa.

5.2. Datasets

In this study, 63 datasets of measured shear strength parameters collected from the Xiaolangdi hydropower station [42] are used to examine the marginal distributions and to fit the joint behaviour of the soil shear strength. The K-S test results are listed in Table 2. Obviously, and coincide, respectively, with the normal distribution and lognormal distribution well. The gamma distribution is rejected. Figure 4(a) shows the scatter plot of datasets, and the PDFs and CDFs of different distribution forms are given in Figures 4(b), 4(c), 4(d), and 4(e).

From these 63 datasets, we can conclude that the mean and coefficient of variation of are 65.97 kPa and 0.3, respectively. Similarly, the average value of is 0.42, with the coefficient of variation of 0.15. The Pearson correlation coefficient between and is −0.5; correspondingly, Kendall’s is −0.4035. In the analysis, the shear strength parameters, and , are selected as random variables. The slope shape parameters and unit weight are assumed as deterministic.

Considering some copulas that do not exhibit negative parameter correlations, the variable can be transformed into , where is the average value. It is worth keeping in mind that both and have the same mean and standard deviation. Also the correlation coefficient between and will have the same magnitude but opposite in direction to that of and . The scatter diagram of and is shown in Figure 6. To compute the probability of failure, the -line and its failure domain must be correspondingly transformed. The shaded part in Figure 5 shows the corresponding failure domain of the transformed -line.

5.3. Probability of Failure Using the CDIM

Based on those measured datasets, (10) is employed for the aforementioned copulas to obtain Kendall’s rank correlation coefficient , and then the parameters of the copulas can be calculated using (8), (9), and the maximum likelihood method (MLE) [30], as shown in Table 3. The setting for the optimum copula is the minimum and maximum . In this paper, is introduced to evaluate the goodness of fit of the copulas. The smallest value of indicates the optimum copula function. Using similar methods proposed by different researchers [38, 39, 41], the AICc, , RMSE, and bias of the copulas are computed by (11)–(16), respectively, and tabulated in Table 3. The Frank copula exhibited smaller AICc, , and bias values, whereas the Gaussian copula exhibited smaller RMSE values. The optimum copula is not unique for different evaluation approaches. Figure 6 shows the contour plots of the joint PDFs for the transformed shear strength parameters associated with different copulas. The following points can be made following a careful investigation of Table 3 and Figure 6. Frank copulas showed the best fits, and the Gaussian, Gumbel, and Plackett copulas demonstrated better fits than the Clayton copulas.

Using the CDIM, one can determine the probability of slope failure by substituting the correlation parameter (from Table 3) into (23). The corresponding probability of failure is calculated by using FORM and MCS, considering the performance function represented by the -line, as shown in Table 4. In the analyses of MCS, 10000 random numbers of and are generated by using the method proposed by Ilich [43]. The probabilities of failure computed using the CDIM, FORM, and MCS yield similar results. The maximum relative error values for the Clayton copula and Gumbel copula are 7.67% and 5.42%, respectively.

6. Discussion

As previously discussed, the optimum copula is not unique for different evaluation methods. Therefore, AICc, , RMSE, and bias are considered as the evaluation index, and the weights of the selected five copulas are calculated and tabulated in Table 5 using the entropy analysis method [44]. The bivariate EW copula for the shear strength parameters ( and ) can be expressed in terms of the weighted combination and density function as follows:where , , and are the th selected weight, density function, and parameter of the copula function, respectively, . Similarly, the term in (24) is integrated for the -line failure domain. Consequently, the probability of failure of a slope system is determined for the considered EW copula.

In addition, according to the example slope model, five-group factors of safety with different and values under the same coefficient of variation are given in Table 5 to further investigate the variation characteristics of the probability of failure with the factor of safety under different copulas. The probabilities of failure of the example slope under different copulas, including the EW copula, and the ratios () of the probabilities of failure calculated by the five copulas compared to that calculated by the EW copula are also shown in Table 5. For a more intuitive comparison, Figure 7 shows a bar graph of the probability of failure of the slope corresponding to each copula function under different factors of safety.

The probability of failure clearly decreases with increases in the factor of safety, but the probability of failure differs notably between each copula function with different factors of safety. The computed results of the different copulas show insignificant discrepancies for a relatively small factor of safety. Nonetheless, the variation increased rapidly for higher factors of safety. For example, the probability of failure of the Plackett copula is 9-fold higher than that of the Clayton copula (when ); in addition, the value of gradually deviated from “1” and showed progressive discrepancies. For a higher factor of safety, the Gaussian and Clayton copulas demonstrated relatively small probabilities of failure. In contrast, the Frank and Plackett copulas exhibited relatively higher probabilities of failure. Therefore, the Frank and Plackett copulas are relatively conservative, whereas the Gaussian and Clayton copulas overestimate the reliability of the slope.

The different structures of the copulas lead to different probabilities of failure. In particular, for a low probability of failure (i.e., high factor of safety), the result of the reliability analysis shows greater sensitivity to the type of copula function. The EW copula combines the results of different copula functions and effectively controls the discrepancies between the probabilities of failure of the slopes. Thus, the EW copula also can avoid over- or underestimating the slope reliability.

7. Conclusions

In this paper, a CDIM was used to compute the probability of failure of a slope system. The failure domain of a homogeneous slope was determined by using the -line. In addition, the bivariate joint CDFs or PDFs for shear strength parameters were established based on the copula function. Based on the main findings of this paper, the following conclusions are noted:(1)For a certain homogeneous slope, the -line characterizes its limit equilibrium state. This phenomenon can be regarded as the inherent property of the slope. A quadratic polynomial can accurately fit the homogeneous slope -line, and the coefficient of determination ranges from 0.92 to 0.99. Higher-order polynomials exhibit a better fitting effect for high and steep slopes. The slope failure domain determined using the -line can be represented by shear strength parameters ( or ).(2)The computed probabilities of failure using the CDIM, FORM, and MCS yield highly similar results. The maximum relative error values for the Clayton copula and Gumbel copula are 7.67% and 5.42%, respectively. This result validated the rationality of CDIM method. The CDIM provides a further approach to slope reliability analysis. However, further research is required to compute the probability of failure for heterogeneous slopes.(3)The computed results of the copula functions showed minor discrepancies for relatively small factors of safety. Nonetheless, the variation increased rapidly for higher values of factor of safety. The EW copula combines the results of different copula functions and effectively controls the discrepancies between the probability of failure of slopes. The EW copula can also eliminate bias in slope reliability estimates.

Competing Interests

The authors declare that they have no competing interests.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Project no. 51439003) and the Public Welfare Special Scientific Research Project of the China Ministry of Water Resources (Project no. 201401029).