# A Copula-Based Method for Estimating Shear Strength Parameters of Rock Mass

Authors: Da Huang; Chao Yang; Bin Zeng; Guoyang Fu
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

The shear strength parameters (i.e., the internal friction coefficientf and cohesion c) are very important in rock engineering, especially for the stability analysis and reinforcement design of slopes and underground caverns. In this paper, a probabilistic method, Copula-based method, is proposed for estimating the shear strength parameters of rock mass. The optimal Copula functions between rock mass quality Q and f, Q and c for the marbles are established based on the correlation analyses of the results of 12 sets of in situ tests in the exploration adits of Jinping I-Stage Hydropower Station. Although the Copula functions are derived from the in situ tests for the marbles, they can be extended to be applied to other types of rock mass with similar geological and mechanical properties. For another 9 sets of in situ tests as an extensional application, by comparison with the results from Hoek-Brown criterion, the estimated values of f and c from the Copula-based method achieve better accuracy. Therefore, the proposed Copula-based method is an effective tool in estimating rock strength parameters.


---

## Body

## 1. Introduction

The shear strength parameters including friction coefficientf (the tangent of friction angle φ) and cohesion c are the most fundamental mechanical parameters for rock mass. They are of very high importance for stability evaluation and support design of rock mass in rock engineering, especially for landslide risk mitigation. Different approaches including laboratory tests, in situ tests, back analysis, and empirical equations can be used to estimate these strength parameters. In the laboratory tests, scale effects are involved since rock specimens are of small sizes [1, 2]. The back analysis method is also recommended by [3, 4] the fact that it might be more reliable than laboratory or in situ test method. Some empirical equations between shear strength parameters and geological indices are widely used due to simplicity and practicality [5–13]. The empirical equations based on famous Hoek-Brown criterion are [13]
(1)f=tansin-1⁡6ambs+mbσ3na-121+a2+a+6ambs+mbσ3na-1,c=σc1+2as+1-ambσ3ns+mbσ3na-11+a2+a1+6ambs+mbσ3na-1/(1+a)(2+a),
where mb is a reduction of the rock material constant mi and obtained by
(2)mb=miexp⁡⁡GSI-10028-14D,
where GSI is the Geological Strength Index and D is a disturbance factor, depending on the degree of disturbance by which the rock mass has been subjected to blast damage and stress relaxation. The D varies from zero for undisturbed in situ rock mass to 1 for very disturbed rock mass. s and a in (1) can be calculated by
(3)a=12+16e-GSI/15-e-20/3,s=exp⁡⁡GSI-1009-3D.In (1), σc is uniaxial compressive strength, σ3n=σ3max⁡⁡/σc, and σ3max⁡ is the upper-limit value of confining stress discussed by Hoek et al. [13].It is important to take full advantage of in situ test data in order to establish the correlativity between geological indices and mechanical parameters of rock mass. With regard to the correlativity of the above variables, an effective way is to establish their joint distribution function. The common statistical methods for obtaining joint distribution function of variables are generally based on two assumptions; that is, the variables have the same marginal distributions and follow the normal distribution [14]. However, rock mass parameters do not always follow the normal distribution especially as the number of samples is small [15]. With respect to small samples, normal information spread estimation method (NISEM) [15] is suitable. Meanwhile, Copula theory [16–21] is a fairly effective method to establish the joint distribution function of nonnormal distribution variables. The advantages of Copula theory are that the marginal distributions and dependence structure of variables can be studied separately and the variables can follow any distributions. The Copula theory has been widely used in many fields such as insurance [16], finance [17], hydrology [18–20], and geotechnical/geological engineering [21].In this paper, the correlations between the rock mass quality indices (i.e., rock quality designation (RQD), rock mass quality (Q), and Rock Block Index (RBI)) and the shear strength parameters (f and c) are investigated based on the results of in situ direct shear tests at Jinping I-Stage Hydropower Station in China. It shows that Q presents the strongest correlation with the shear strength parameters. Then, the joint distribution functions of Q and f, Q and c are proposed based on Copula theory. By means of the functions and the given guarantee rates, the shear strength parameters can be estimated. Compared with those estimated by Hoek-Brown criterion, the estimated values using the proposed Copula-based method are closer to the in situ tests ones.


## 2. Correlations between Rock Mass Quality Indices and Shear Strength

### 2.1. In Situ Test Data

The Jinping Hydropower Station is located on the middle reach of the Yalong River, in Sichuan, China. It has a double-curvature arch dam with a height of 305 m and a total installed capacity of 3300 MW [22, 23]. The Hydropower Station is located in the east edge of the Tibet Plateau tectonic belt. Due to the rapid crust uplift during the Quaternary, the steep slopes exist on both sides of the valley deeply cut by Yalong River. The stratum in the dam zone consists of a series of metamorphic rocks of middle of Triassic period, of which the great majority are marbles with grey or white and fine-grained texture. The main mineral compositions of the marble are calcite and dolomite. The uniaxial compressive strength of the fresh marble is about 60–75 MPa.In order to investigate the shear strength parameters of the representative marbles, 12 sets of in situ direct shear tests with five similar specimens in each set were performed in the exploration adits of Jinping I-Stage Hydropower Station. Figure1 shows the sketch of in situ direct shear test which was performed using flat stacking method [24]. The side length of each cube specimen is 1 m which contains in general two sets of joints (bedding joints and cross joints; see Figure 1). RQD and RBI [25, 26] of each specimen can be obtained using scanline field mapping method. The scanline used for measuring RQD was parallel to shear direction. The shear direction was parallel to the bedding joint along which the shear deformation happened. The five similar specimens in a set were tested at different normal stress σn from 2 to 10 MPa. The shear strength parameters (f and c) of each set are obtained by linear regression according to Mohr-Coulomb criterion. A good fitting straight line between the peak shear stress τ and normal stress σn can be derived from each set test data for which the R2 of fitting degree is greater than 0.9 (see Figure 2). The test results are shown in Table 1, in which the rock mass quality indices are the average of each set samples, respectively.Table 1
Results of in situ direct shear tests and rock mass quality indices of the specimens.
Set number
Rock mass quality indices
Shear strength parameters
RQD/%
RBI
Q
f
c/MPa
1
75.3
17.6
30.5
1.79
1.95
2
89.3
31.5
23.8
1.38
1.98
3
96.8
41.6
38.7
1.17
1.50
4
99.3
31.8
39.7
2.20
2.25
5
96.4
43.4
38.6
1.39
1.80
6
98.1
22.6
26.5
1.60
1.80
7
59.8
8.0
8.1
1.80
2.80
8
56.0
10.3
5.2
1.17
1.20
9
65.9
14.3
8.9
1.23
2.41
10
77.4
13.2
5.2
0.56
3.00
11
67.9
9.0
5.4
1.07
3.25
12
51.8
7.0
4.6
1.15
4.32Figure 1
The sketch of the in situ direct shear test.Figure 2
Linear regression of normal stressσn versus shear stress τ for set number 1 in Table 1.


### 2.2. Correlation Measure

Pearson linear correlation coefficientγn [27], Spearman’s rank correlation coefficient ρn [28], and Kendall’s rank correlation coefficient τn [29] are widely used to measure how relevant the variables are. They can be calculated by
(4)γn=∑i=1nx1i-x1¯x2i-x2¯ns12·s22,ρn=12n(n-1)(n+1)∑i=1nRiSi-3n+1n-1,τn=n2-1∑i<jsign⁡x1i-x1jx2i-x2j,
where n is the number of observations, x1i(x1j) and x2i(x2j) are the observed values of variables (x1 and x2), x1¯ and x2¯ are sample mean, s12 and s22 are sample variance, Ri and Si are the rank of x1i and x2i, sign⁡=1 if x1i<x1j and x2i<x2j, and sign⁡=0 if x1i=x1j and x2i=x2j; otherwise, sign⁡=-1, i, j=1,2,…,n.Table2 shows the results of the γn, ρn, and τn of test data. There are clearly positive correlations (the same variation trend) between three rock mass quality indices (RQD, RBI, and Q) and f. The correlation coefficients between Q and f are much greater than those between RQD and f, RBI and f. There is clearly negative correlation (changes in the opposite trend) between the three indices and c. The correlation between RBI and c is slightly better than Q and c, while RQD and c have the worst correlation. In summary, Q has the best correlation with shear strength parameters (f and c). In fact, RQD and RBI evaluate the rock mass quality only by the integrality of rock core. While many other factors, such as joint surface roughness and degree of weathering for dominantly oriented joint sets, joint water reduction, stress reduction factor, and volumetric joint count, are fully involved in Q-system, the Q-system, therefore, provides a comprehensive description on rock mass quality [30]. The shear strength parameters of jointed rock mass are generally controlled by complex geological conditions and mechanical behaviors of rock material and so on. The best correlation between Q and the shear strength parameters in this statistical study also indicates that Q-system associated with multiple factors of influence on c and f of jointed rocks more precisely depicts the quality of rock mass.Table 2
Correlation coefficients of the rock mass quality indices and the shear strength.
Variables
Pearsonγn
Spearmanρn
Kendallτn
RQD-f
0.307
0.386
0.303
RBI-f
0.241
0.376
0.303
Q-f
0.550
0.656
0.545
RQD-c
−0.517
−0.397
−0.333
RBI-c
−0.573
−0.635
−0.515
Q-c
−0.563
−0.517
−0.485The correlation measure indicates that, for the in situ tests presented in this paper, the marbles with higher rock mass quality indices show higher friction angle and lower cohesion. Many researchers had also found that, for rocks and soils, there was a negative correlation betweenf and c [31, 32]. If this negative correlation is true between f and c, it is easily understood that a variable has a positive correlation with one (f) of shear strength parameters while a negative correlation must appear with the other (c) based on the pure mathematical logic. As shown in Figure 1, there are some intact rock segments on the sheared bedding joint in the case of relatively low Q value due to cross joint, which might explain why higher Q has higher f and lower c. The shear failure occurs along the weak bedding joint with a relatively low c value on condition that the tested specimen has a high Q value (good rock mass quality). When a low Q value appears, these intact rock block segments on the shear surface must be thus cut through before reaching peak strength, which can increase the cohesion c of the failure surface.


## 2.1. In Situ Test Data

The Jinping Hydropower Station is located on the middle reach of the Yalong River, in Sichuan, China. It has a double-curvature arch dam with a height of 305 m and a total installed capacity of 3300 MW [22, 23]. The Hydropower Station is located in the east edge of the Tibet Plateau tectonic belt. Due to the rapid crust uplift during the Quaternary, the steep slopes exist on both sides of the valley deeply cut by Yalong River. The stratum in the dam zone consists of a series of metamorphic rocks of middle of Triassic period, of which the great majority are marbles with grey or white and fine-grained texture. The main mineral compositions of the marble are calcite and dolomite. The uniaxial compressive strength of the fresh marble is about 60–75 MPa.In order to investigate the shear strength parameters of the representative marbles, 12 sets of in situ direct shear tests with five similar specimens in each set were performed in the exploration adits of Jinping I-Stage Hydropower Station. Figure1 shows the sketch of in situ direct shear test which was performed using flat stacking method [24]. The side length of each cube specimen is 1 m which contains in general two sets of joints (bedding joints and cross joints; see Figure 1). RQD and RBI [25, 26] of each specimen can be obtained using scanline field mapping method. The scanline used for measuring RQD was parallel to shear direction. The shear direction was parallel to the bedding joint along which the shear deformation happened. The five similar specimens in a set were tested at different normal stress σn from 2 to 10 MPa. The shear strength parameters (f and c) of each set are obtained by linear regression according to Mohr-Coulomb criterion. A good fitting straight line between the peak shear stress τ and normal stress σn can be derived from each set test data for which the R2 of fitting degree is greater than 0.9 (see Figure 2). The test results are shown in Table 1, in which the rock mass quality indices are the average of each set samples, respectively.Table 1
Results of in situ direct shear tests and rock mass quality indices of the specimens.
Set number
Rock mass quality indices
Shear strength parameters
RQD/%
RBI
Q
f
c/MPa
1
75.3
17.6
30.5
1.79
1.95
2
89.3
31.5
23.8
1.38
1.98
3
96.8
41.6
38.7
1.17
1.50
4
99.3
31.8
39.7
2.20
2.25
5
96.4
43.4
38.6
1.39
1.80
6
98.1
22.6
26.5
1.60
1.80
7
59.8
8.0
8.1
1.80
2.80
8
56.0
10.3
5.2
1.17
1.20
9
65.9
14.3
8.9
1.23
2.41
10
77.4
13.2
5.2
0.56
3.00
11
67.9
9.0
5.4
1.07
3.25
12
51.8
7.0
4.6
1.15
4.32Figure 1
The sketch of the in situ direct shear test.Figure 2
Linear regression of normal stressσn versus shear stress τ for set number 1 in Table 1.


## 2.2. Correlation Measure

Pearson linear correlation coefficientγn [27], Spearman’s rank correlation coefficient ρn [28], and Kendall’s rank correlation coefficient τn [29] are widely used to measure how relevant the variables are. They can be calculated by
(4)γn=∑i=1nx1i-x1¯x2i-x2¯ns12·s22,ρn=12n(n-1)(n+1)∑i=1nRiSi-3n+1n-1,τn=n2-1∑i<jsign⁡x1i-x1jx2i-x2j,
where n is the number of observations, x1i(x1j) and x2i(x2j) are the observed values of variables (x1 and x2), x1¯ and x2¯ are sample mean, s12 and s22 are sample variance, Ri and Si are the rank of x1i and x2i, sign⁡=1 if x1i<x1j and x2i<x2j, and sign⁡=0 if x1i=x1j and x2i=x2j; otherwise, sign⁡=-1, i, j=1,2,…,n.Table2 shows the results of the γn, ρn, and τn of test data. There are clearly positive correlations (the same variation trend) between three rock mass quality indices (RQD, RBI, and Q) and f. The correlation coefficients between Q and f are much greater than those between RQD and f, RBI and f. There is clearly negative correlation (changes in the opposite trend) between the three indices and c. The correlation between RBI and c is slightly better than Q and c, while RQD and c have the worst correlation. In summary, Q has the best correlation with shear strength parameters (f and c). In fact, RQD and RBI evaluate the rock mass quality only by the integrality of rock core. While many other factors, such as joint surface roughness and degree of weathering for dominantly oriented joint sets, joint water reduction, stress reduction factor, and volumetric joint count, are fully involved in Q-system, the Q-system, therefore, provides a comprehensive description on rock mass quality [30]. The shear strength parameters of jointed rock mass are generally controlled by complex geological conditions and mechanical behaviors of rock material and so on. The best correlation between Q and the shear strength parameters in this statistical study also indicates that Q-system associated with multiple factors of influence on c and f of jointed rocks more precisely depicts the quality of rock mass.Table 2
Correlation coefficients of the rock mass quality indices and the shear strength.
Variables
Pearsonγn
Spearmanρn
Kendallτn
RQD-f
0.307
0.386
0.303
RBI-f
0.241
0.376
0.303
Q-f
0.550
0.656
0.545
RQD-c
−0.517
−0.397
−0.333
RBI-c
−0.573
−0.635
−0.515
Q-c
−0.563
−0.517
−0.485The correlation measure indicates that, for the in situ tests presented in this paper, the marbles with higher rock mass quality indices show higher friction angle and lower cohesion. Many researchers had also found that, for rocks and soils, there was a negative correlation betweenf and c [31, 32]. If this negative correlation is true between f and c, it is easily understood that a variable has a positive correlation with one (f) of shear strength parameters while a negative correlation must appear with the other (c) based on the pure mathematical logic. As shown in Figure 1, there are some intact rock segments on the sheared bedding joint in the case of relatively low Q value due to cross joint, which might explain why higher Q has higher f and lower c. The shear failure occurs along the weak bedding joint with a relatively low c value on condition that the tested specimen has a high Q value (good rock mass quality). When a low Q value appears, these intact rock block segments on the shear surface must be thus cut through before reaching peak strength, which can increase the cohesion c of the failure surface.


## 3. Copula-Based Method

### 3.1. Copula Theory

Considering multivariate random variablesX1,…,Xi,…,Xk with marginal distribution functions as FXi(xi)=PXi(Xi≤xi), where k is the number of random variables and xi is the value of random variable Xi (X1,…,Xi,…,Xk), FXi(xi) is the marginal distribution function and PXi(Xi≤xi) is the probability when Xi≤xi. The joint distribution function H of the multivariate random variables X1,X2,…,Xk can been expressed as [33]
(5)HX1,…Xi,…Xkx1,…,xi,…,xk=PX1≤x1,…,Xi≤xi,…,Xk≤xk.The Copula function is to connect multivariate probability distributions to their marginal probability distributions. Thus, the multivariate joint distributionH is expressed in terms of its marginal distribution function FXi(xi) and the associated dependence structure C as
(6)HX1,…Xi,…Xkx1,…,xi,…,xk=CFX1x1,…,FXixi,…,FXkxk,
where C, called Copula function, is a uniquely determined mapping whenever FXi(xi) are continuous and captures the essential features of the dependence between random variables. C is essentially the joint distribution of multivariate random variables of X1,…,Xi,…,Xk through their marginal distributions FX1(x1),…,FXi(xi),…,FXk(xk) [33].The joint distribution function established using Copula-based method has two major advantages: the marginal distribution and dependence structure can be studied separately, and the variables can follow any distributions. Therefore, the problem of determiningH reduces to determining marginal distribution of variables and Copula function C.In order to clearly demonstrate the calculation procedure using the present Copula-based method, the calculation flow chart is shown in Figure3. The detailed procedures are as follows.Figure 3
The calculation flow chart of Copula-based method.Step 1.
The correlation between target variables(f,c) and relevant variables (Q, RQD, and RBI) is analyzed according to the correlation coefficients obtained by (4) or else. The best relevant variable (i.e., Q) can be selected, as shown in Section 2.2.Step 2.
Utilize suitable functional forms and methods (such as normal information spread estimation method (NISEM)) to establish marginal distribution functionFXi(xi) of variables which could be accepted according to corresponding test (e.g., Kolmogorov-Smirnov (K-S) test). It is addressed in Section 3.2.Step 3.
Introduce suitable multivariable joint distribution functional forms (i.e., Copula functions) and confirm the optimal Copula functionC indicated by a minimum AIC value. Such is stated in Section 3.3.Step 4.
Estimate the values of target variables(f,c) using the Copula method on condition that a guarantee rate β and the best relevant variable Q have been given, which is presented in Section 4.


### 3.2. Verification of Marginal Distribution Functions

For small samples, the normal information spread estimation method (NISEM) [15] is an effective way to generate probability density functions. The main difference between NISEM and the available classical methods [15, 34] is that the variables are not assumed to follow a distribution (such as normal, log normal, and exponential distribution) in advance. It has been proved that the NISEM used for generating the probability density functions of geotechnical parameters with small samples is more accurate and effective than the other classical methods [15]. The probability density function of random variables Xi by using NISEM is expressed as
(7)fxi=12πnh∑j=1nexp⁡-xi-xij22h2,
where h is the window width of standard normal information spread function and xij is the observed value of variable Xi. By the approaching principle of normal information diffusion, h can be calculated by
(8)h=γximax⁡-ximin⁡n-1,
where ximax⁡ and ximin⁡ are the maximal and minimum observed values of variable Xi. The value of γ is related to n, as shown in Table 3 [35].Table 3
Values ofγ  for different values of n [35].
n
γ
n
γ
3
0.849 321 800
11
1.420 835 443
4
1.273 982 782
12
1.420 269 570
5
1.698 643 675
13
1.420 698 795
6
1.336 252 561
14
1.420 669 671
7
1.445 461 208
15
1.420 693 321
8
1.395 189 816
16
1.420 692 226
9
1.422 962 345
17
1.420 693 101
10
1.416 278 786
>17
1.420 693 101For the variables in Table1, γ=1.420269570 while n=12. Therefore, the estimation of probability density functions can be obtained from (7). The values of 1/[2π0.5nh] and 1/2h2 (parameters of f(xi)) are shown in Table 4.Table 4
Values of1/(2πnh) and 1/2h2.
Variables
1
/
(
2
π
n
h
)
1
/
2
h
2
RQD
0.005 420 722
0.013 293 110
RBI
0.007 073 744
0.022 636 592
Q
0.007 335 735
0.024 344 428
f
0.157 002 612
11.151 315 78
c
0.082 527 014
3.081 091 687By means of the probability density functionsf(xi) (see (7)), the marginal distribution functions FXi(xi) can be established as
(9)FXixi=∫-∞+∞fxidxi.Kolmogorov-Smirnov (K-S) test commonly used to check whether a set of data follows a certain distribution [36] is adopted to test the established marginal distribution functions. The procedure of K-S test involves forming the cumulative frequency distribution F~Xi(xi) and computing the test statistic Dn:
(10)Dn=max⁡FXi(xi)-F~Xi(xi).
The cumulative frequency distribution F~Xi(xi) is defined as
(11)F~Xixi=1n∑i=1nIxij≤xi,
where Ixij≤xi is the indicator function, equal to 1 if xij≤xi and equal to 0 otherwise.The established marginal distribution functionsFXi(xi) will be accepted if Dn<Dn,α, otherwise they will be denied. Where Dn,α is the critical value of a given significance level α (usually α=5% [37]), the Dn,α can be obtained using the table established by [38] when n and α have been determined. Here, the significance level α is 5% and n is 12, so Dn,α=0.375 based on the relevant table in [38]. The K-S test results are shown in Table 5. Obviously, all the established marginal distribution functions can be accepted.Table 5
Results of Kolmogorov-Smirnov test.
Variables
n
α
D
n
,
α
D
n
Accept or reject
RQD
12
0.05
0.375
0.13617
Accept
RBI
12
0.05
0.375
0.09205
Accept
Q
12
0.05
0.375
0.11670
Accept
f
12
0.05
0.375
0.12352
Accept
c
12
0.05
0.375
0.11224
Accept


### 3.3. Determination of the Optimal Copula Functions

In order to obtain the optimal Copula functions, the common method of the goodness-of-fit is adopted. Therefore, two steps are taken to determine the optimal Copula function: (1) selecting the suitably qualified Copula functions according to the characteristics of marginal distribution functions and dependent structure of variables and (2) determining the parameters of Copula functions selected.With respect to two-dimensional Copula functions, the parameterθ of Copula functions can be obtained by means of Kendall’s rank correlation coefficient τn, which has been calculated in (4). The relation between τn and θ is [39]
(12)τn=4∫01∫01Cui,vi;θdC(ui,vi;θ)-1,
where ui and vi are the marginal distribution functional values on condition that the values of two variables are equal to x1i and x2i (i.e., ui=FX1(x1i) and vi=FX2(x2i)), while C(ui,vi,θ) is the introduced Copula functional value with the parameter θ.Akaike information criteria (AIC) [40] are commonly used for goodness-of-fit test. AIC are based on entropy principle and are a simple and effective method for goodness-of-fit test [41]. And the optimal Copula is the one which has the minimum AIC value; the value of AIC can be calculated as
(13)AIC=nln⁡1n-t∑i=1nFempx1i,x2i-Cui,vi2+2t,
where t is the number of parameters of Copula function (t=1 for two-dimensional Copula function) and C(ui,vi) is the Copula functional value. Femp(x1i,x2i) is the empirical probability functional value; it can be obtained as
(14)Fempx1i,x2i=PX1≤x1i,X2≤x2i=1n+1∑j=1n∑l=1nNj,l,
where Nj,l is the number of joint observations under the condition of X1≤x1i as well as X2≤x2i.In order to visualize the dependence structure of the measured data (x1,x2) and to obtain an insight on a suitable Copula, the measured data in original space should be transformed into the standard uniform random vector T=(T1,T2). Tang et al. [21] defined the random vector T by adopting the empirical distributions of measured data. T1 and T2 are expressed as
(15)t1i=rank⁡(x1i)n+1,t2i=rank⁡(x2i)n+1,
where rank⁡(x1i) (or rank⁡(x2i), (x3i)) denotes the rank of x1i (or x2i, x3i) among x1 (or x2, x3i) in an ascending order. (t1i,t2i) are the values of variables (T1,T2) or (T1,T3) corresponding to (Q,f) or (Q,c).The scatter plots of empirical cumulative distribution of(Q,f) and (Q,c) are shown in Figure 4. Q and f,Q and c are substantially symmetric about the diagonal lines. Although the data points (Q and c) in right Figure 4(b) seem to be not symmetric as well as those in left Figure 4(a), most of them are located around diagonal line (i.e., symmetrical line) in Figure 4(a). Archimedean Copulas have a wide range of application [19–21, 42] due to several competitive advantages: (1) the ease with which they can be constructed; (2) the great variety of families of copulas which belong to this class; and (3) the many nice properties possessed by members of this class [39]. As results, the symmetrical Archimedean Copulas, Nelsen No. 1, Nelsen No. 2, Nelsen No. 4, Nelsen No. 5, Nelsen No. 12, Nelsen No. 14, Nelsen No. 15, and Nelsen No. 18 [39] are selected to model the positive correlation between Q and f, while Nelsen No. 1, Nelsen No. 2, Nelsen No. 5, Nelsen No. 7, Nelsen No. 8, and Nelsen No. 15 [39] are selected to model the negative correlation between Q and c. By means of AIC calculated by (13), the goodness-of-fit test is used to determine the optimal Copula functions. The chosen symmetrical Archimedean Copula functions and the values of AIC are shown in Table 6. It shows that Nelsen No. 1 has the minimum AIC value for establishing the joint distribution function of Q and f, while Nelsen No. 2 shows the minimum for that of Q and c. Therefore, Nelsen No. 1 and Nelsen No. 2 can be used as the optimal Copula functions to construct the joint distribution functions H of Q-f and Q-c, respectively.Table 6
The parameterθ of 2D Copula functions and the values of AIC.
Copula type
C
(
u
,
v
)
θ (Q-f)
AIC (Q-f)
θ (Q-c)
AIC (Q-c)
Nelsen No. 1
max
⁡
[
(
u
-
θ
+
v
-
θ
-
1
)
-
1
/
θ
,
0
]
2.396
−67.990
−0.653
−61.777
Nelsen No. 2
max
⁡
{
1
-
[
(
1
-
u
)
θ
+
(
1
-
v
)
θ
]
1
/
θ
,
0
}
4.396
−19.240
1.347
−66.184
Nelsen No. 4
e
-
[
(
-
ln
⁡
u
)
θ
+
(
-
ln
⁡
v
)
θ
]
1
/
θ
2.198
−66.341
—
—
Nelsen No. 5
-
1
θ
ln
⁡
[
1
+
(
e
-
θ
u
-
1
)
(
e
-
θ
v
-
1
)
e
-
θ
-
1
]
4.978
−65.175
−5.442
−59.637
Nelsen No. 7
max
⁡
[
θ
u
v
-
(
1
-
θ
)
(
u
+
v
-
1
)
,
0
]
—
—
0.634
−44.658
Nelsen No. 8
max
⁡
[
θ
2
u
v
-
(
1
-
u
)
(
1
-
v
)
θ
2
-
(
θ
-
1
)
2
(
1
-
u
)
(
1
-
v
)
,
0
]
—
—
1.629
−64.809
Nelsen No. 12
{
1
+
[
(
u
-
1
-
1
)
θ
+
(
v
-
1
-
1
)
θ
]
1
/
θ
}
-
1
1.465
−67.648
—
—
Nelsen No. 14
{
1
+
[
(
u
-
1
/
θ
-
1
)
θ
+
(
v
-
1
/
θ
-
1
)
θ
]
1
/
θ
}
-
θ
1.698
−67.287
—
—
Nelsen No. 15
max
⁡
(
{
1
-
[
(
1
-
u
-
1
/
θ
)
θ
+
(
1
-
v
-
1
/
θ
)
θ
]
1
/
θ
}
θ
,
0
)
2.698
−65.360
1.173
−60.240
Nelsen No. 18
max
⁡
{
1
+
θ
ln
⁡
(
e
θ
/
u
-
1
+
e
θ
/
v
-
1
)
,
0
}
2.930
−60.246
—
—The scatter plots of empirical distribution: (a)T1/Q versus T2/f and (b) T1/Q versus T3/c.
(a)
(b)


## 3.1. Copula Theory

Considering multivariate random variablesX1,…,Xi,…,Xk with marginal distribution functions as FXi(xi)=PXi(Xi≤xi), where k is the number of random variables and xi is the value of random variable Xi (X1,…,Xi,…,Xk), FXi(xi) is the marginal distribution function and PXi(Xi≤xi) is the probability when Xi≤xi. The joint distribution function H of the multivariate random variables X1,X2,…,Xk can been expressed as [33]
(5)HX1,…Xi,…Xkx1,…,xi,…,xk=PX1≤x1,…,Xi≤xi,…,Xk≤xk.The Copula function is to connect multivariate probability distributions to their marginal probability distributions. Thus, the multivariate joint distributionH is expressed in terms of its marginal distribution function FXi(xi) and the associated dependence structure C as
(6)HX1,…Xi,…Xkx1,…,xi,…,xk=CFX1x1,…,FXixi,…,FXkxk,
where C, called Copula function, is a uniquely determined mapping whenever FXi(xi) are continuous and captures the essential features of the dependence between random variables. C is essentially the joint distribution of multivariate random variables of X1,…,Xi,…,Xk through their marginal distributions FX1(x1),…,FXi(xi),…,FXk(xk) [33].The joint distribution function established using Copula-based method has two major advantages: the marginal distribution and dependence structure can be studied separately, and the variables can follow any distributions. Therefore, the problem of determiningH reduces to determining marginal distribution of variables and Copula function C.In order to clearly demonstrate the calculation procedure using the present Copula-based method, the calculation flow chart is shown in Figure3. The detailed procedures are as follows.Figure 3
The calculation flow chart of Copula-based method.Step 1.
The correlation between target variables(f,c) and relevant variables (Q, RQD, and RBI) is analyzed according to the correlation coefficients obtained by (4) or else. The best relevant variable (i.e., Q) can be selected, as shown in Section 2.2.Step 2.
Utilize suitable functional forms and methods (such as normal information spread estimation method (NISEM)) to establish marginal distribution functionFXi(xi) of variables which could be accepted according to corresponding test (e.g., Kolmogorov-Smirnov (K-S) test). It is addressed in Section 3.2.Step 3.
Introduce suitable multivariable joint distribution functional forms (i.e., Copula functions) and confirm the optimal Copula functionC indicated by a minimum AIC value. Such is stated in Section 3.3.Step 4.
Estimate the values of target variables(f,c) using the Copula method on condition that a guarantee rate β and the best relevant variable Q have been given, which is presented in Section 4.


## 3.2. Verification of Marginal Distribution Functions

For small samples, the normal information spread estimation method (NISEM) [15] is an effective way to generate probability density functions. The main difference between NISEM and the available classical methods [15, 34] is that the variables are not assumed to follow a distribution (such as normal, log normal, and exponential distribution) in advance. It has been proved that the NISEM used for generating the probability density functions of geotechnical parameters with small samples is more accurate and effective than the other classical methods [15]. The probability density function of random variables Xi by using NISEM is expressed as
(7)fxi=12πnh∑j=1nexp⁡-xi-xij22h2,
where h is the window width of standard normal information spread function and xij is the observed value of variable Xi. By the approaching principle of normal information diffusion, h can be calculated by
(8)h=γximax⁡-ximin⁡n-1,
where ximax⁡ and ximin⁡ are the maximal and minimum observed values of variable Xi. The value of γ is related to n, as shown in Table 3 [35].Table 3
Values ofγ  for different values of n [35].
n
γ
n
γ
3
0.849 321 800
11
1.420 835 443
4
1.273 982 782
12
1.420 269 570
5
1.698 643 675
13
1.420 698 795
6
1.336 252 561
14
1.420 669 671
7
1.445 461 208
15
1.420 693 321
8
1.395 189 816
16
1.420 692 226
9
1.422 962 345
17
1.420 693 101
10
1.416 278 786
>17
1.420 693 101For the variables in Table1, γ=1.420269570 while n=12. Therefore, the estimation of probability density functions can be obtained from (7). The values of 1/[2π0.5nh] and 1/2h2 (parameters of f(xi)) are shown in Table 4.Table 4
Values of1/(2πnh) and 1/2h2.
Variables
1
/
(
2
π
n
h
)
1
/
2
h
2
RQD
0.005 420 722
0.013 293 110
RBI
0.007 073 744
0.022 636 592
Q
0.007 335 735
0.024 344 428
f
0.157 002 612
11.151 315 78
c
0.082 527 014
3.081 091 687By means of the probability density functionsf(xi) (see (7)), the marginal distribution functions FXi(xi) can be established as
(9)FXixi=∫-∞+∞fxidxi.Kolmogorov-Smirnov (K-S) test commonly used to check whether a set of data follows a certain distribution [36] is adopted to test the established marginal distribution functions. The procedure of K-S test involves forming the cumulative frequency distribution F~Xi(xi) and computing the test statistic Dn:
(10)Dn=max⁡FXi(xi)-F~Xi(xi).
The cumulative frequency distribution F~Xi(xi) is defined as
(11)F~Xixi=1n∑i=1nIxij≤xi,
where Ixij≤xi is the indicator function, equal to 1 if xij≤xi and equal to 0 otherwise.The established marginal distribution functionsFXi(xi) will be accepted if Dn<Dn,α, otherwise they will be denied. Where Dn,α is the critical value of a given significance level α (usually α=5% [37]), the Dn,α can be obtained using the table established by [38] when n and α have been determined. Here, the significance level α is 5% and n is 12, so Dn,α=0.375 based on the relevant table in [38]. The K-S test results are shown in Table 5. Obviously, all the established marginal distribution functions can be accepted.Table 5
Results of Kolmogorov-Smirnov test.
Variables
n
α
D
n
,
α
D
n
Accept or reject
RQD
12
0.05
0.375
0.13617
Accept
RBI
12
0.05
0.375
0.09205
Accept
Q
12
0.05
0.375
0.11670
Accept
f
12
0.05
0.375
0.12352
Accept
c
12
0.05
0.375
0.11224
Accept


## 3.3. Determination of the Optimal Copula Functions

In order to obtain the optimal Copula functions, the common method of the goodness-of-fit is adopted. Therefore, two steps are taken to determine the optimal Copula function: (1) selecting the suitably qualified Copula functions according to the characteristics of marginal distribution functions and dependent structure of variables and (2) determining the parameters of Copula functions selected.With respect to two-dimensional Copula functions, the parameterθ of Copula functions can be obtained by means of Kendall’s rank correlation coefficient τn, which has been calculated in (4). The relation between τn and θ is [39]
(12)τn=4∫01∫01Cui,vi;θdC(ui,vi;θ)-1,
where ui and vi are the marginal distribution functional values on condition that the values of two variables are equal to x1i and x2i (i.e., ui=FX1(x1i) and vi=FX2(x2i)), while C(ui,vi,θ) is the introduced Copula functional value with the parameter θ.Akaike information criteria (AIC) [40] are commonly used for goodness-of-fit test. AIC are based on entropy principle and are a simple and effective method for goodness-of-fit test [41]. And the optimal Copula is the one which has the minimum AIC value; the value of AIC can be calculated as
(13)AIC=nln⁡1n-t∑i=1nFempx1i,x2i-Cui,vi2+2t,
where t is the number of parameters of Copula function (t=1 for two-dimensional Copula function) and C(ui,vi) is the Copula functional value. Femp(x1i,x2i) is the empirical probability functional value; it can be obtained as
(14)Fempx1i,x2i=PX1≤x1i,X2≤x2i=1n+1∑j=1n∑l=1nNj,l,
where Nj,l is the number of joint observations under the condition of X1≤x1i as well as X2≤x2i.In order to visualize the dependence structure of the measured data (x1,x2) and to obtain an insight on a suitable Copula, the measured data in original space should be transformed into the standard uniform random vector T=(T1,T2). Tang et al. [21] defined the random vector T by adopting the empirical distributions of measured data. T1 and T2 are expressed as
(15)t1i=rank⁡(x1i)n+1,t2i=rank⁡(x2i)n+1,
where rank⁡(x1i) (or rank⁡(x2i), (x3i)) denotes the rank of x1i (or x2i, x3i) among x1 (or x2, x3i) in an ascending order. (t1i,t2i) are the values of variables (T1,T2) or (T1,T3) corresponding to (Q,f) or (Q,c).The scatter plots of empirical cumulative distribution of(Q,f) and (Q,c) are shown in Figure 4. Q and f,Q and c are substantially symmetric about the diagonal lines. Although the data points (Q and c) in right Figure 4(b) seem to be not symmetric as well as those in left Figure 4(a), most of them are located around diagonal line (i.e., symmetrical line) in Figure 4(a). Archimedean Copulas have a wide range of application [19–21, 42] due to several competitive advantages: (1) the ease with which they can be constructed; (2) the great variety of families of copulas which belong to this class; and (3) the many nice properties possessed by members of this class [39]. As results, the symmetrical Archimedean Copulas, Nelsen No. 1, Nelsen No. 2, Nelsen No. 4, Nelsen No. 5, Nelsen No. 12, Nelsen No. 14, Nelsen No. 15, and Nelsen No. 18 [39] are selected to model the positive correlation between Q and f, while Nelsen No. 1, Nelsen No. 2, Nelsen No. 5, Nelsen No. 7, Nelsen No. 8, and Nelsen No. 15 [39] are selected to model the negative correlation between Q and c. By means of AIC calculated by (13), the goodness-of-fit test is used to determine the optimal Copula functions. The chosen symmetrical Archimedean Copula functions and the values of AIC are shown in Table 6. It shows that Nelsen No. 1 has the minimum AIC value for establishing the joint distribution function of Q and f, while Nelsen No. 2 shows the minimum for that of Q and c. Therefore, Nelsen No. 1 and Nelsen No. 2 can be used as the optimal Copula functions to construct the joint distribution functions H of Q-f and Q-c, respectively.Table 6
The parameterθ of 2D Copula functions and the values of AIC.
Copula type
C
(
u
,
v
)
θ (Q-f)
AIC (Q-f)
θ (Q-c)
AIC (Q-c)
Nelsen No. 1
max
⁡
[
(
u
-
θ
+
v
-
θ
-
1
)
-
1
/
θ
,
0
]
2.396
−67.990
−0.653
−61.777
Nelsen No. 2
max
⁡
{
1
-
[
(
1
-
u
)
θ
+
(
1
-
v
)
θ
]
1
/
θ
,
0
}
4.396
−19.240
1.347
−66.184
Nelsen No. 4
e
-
[
(
-
ln
⁡
u
)
θ
+
(
-
ln
⁡
v
)
θ
]
1
/
θ
2.198
−66.341
—
—
Nelsen No. 5
-
1
θ
ln
⁡
[
1
+
(
e
-
θ
u
-
1
)
(
e
-
θ
v
-
1
)
e
-
θ
-
1
]
4.978
−65.175
−5.442
−59.637
Nelsen No. 7
max
⁡
[
θ
u
v
-
(
1
-
θ
)
(
u
+
v
-
1
)
,
0
]
—
—
0.634
−44.658
Nelsen No. 8
max
⁡
[
θ
2
u
v
-
(
1
-
u
)
(
1
-
v
)
θ
2
-
(
θ
-
1
)
2
(
1
-
u
)
(
1
-
v
)
,
0
]
—
—
1.629
−64.809
Nelsen No. 12
{
1
+
[
(
u
-
1
-
1
)
θ
+
(
v
-
1
-
1
)
θ
]
1
/
θ
}
-
1
1.465
−67.648
—
—
Nelsen No. 14
{
1
+
[
(
u
-
1
/
θ
-
1
)
θ
+
(
v
-
1
/
θ
-
1
)
θ
]
1
/
θ
}
-
θ
1.698
−67.287
—
—
Nelsen No. 15
max
⁡
(
{
1
-
[
(
1
-
u
-
1
/
θ
)
θ
+
(
1
-
v
-
1
/
θ
)
θ
]
1
/
θ
}
θ
,
0
)
2.698
−65.360
1.173
−60.240
Nelsen No. 18
max
⁡
{
1
+
θ
ln
⁡
(
e
θ
/
u
-
1
+
e
θ
/
v
-
1
)
,
0
}
2.930
−60.246
—
—The scatter plots of empirical distribution: (a)T1/Q versus T2/f and (b) T1/Q versus T3/c.
(a)
(b)


## 4. Application to Estimating Shear Strength Parameters

### 4.1. Guarantee Rate

Instability would happen if the strength parameters of rock mass are overvalued, while economic waste could be caused if they are undervalued. Assuming that the estimated values of shear strength parameters aref1 and c1, the actual values are f and c; the guarantee rates βf and βc are then defined as the probability of f>f1 and c>c1, respectively, which can be expressed as
(16)βf=Pf>f1,βc=Pc>c1.The estimated value is overvalued if guarantee ratesβ are low, while being undervalued if high. In consequence, the estimated values can be evaluated by means of the guarantee rates β. On the other hand, the estimated f and c with a certain guarantee rates β can be calculated by solving the conditional probability of the optimal Copula functions. Based on the optimal Copula functions C of Q-f and Q-c established in Section 3.3, βf and βc can be calculated as(17a)βf=Pf>f1∣Q=Q0=1-Pf≤f1∣Q=Q0=1-CU≤u∣V=v=1-lim⁡Δv→0⁡Cu,v+Δv-Cu,vΔv=1-∂∂vCu,vv=v,(17b)βc=Pc>c1∣Q=Q0=1-Pc≤c1∣Q=Q0=1-C(U≤u∣V=v)=1-lim⁡Δv→0⁡Cu,v+Δv-Cu,vΔv=1-∂∂vCu,vv=v,where Q0 is the measured value of Q.For direct shear test, in general, the predominant slope method and the least square method [43] are used to statistically process the test data. Using the least square method, the best linear fitting relations of normal stress σn and shear stress τ for each class of rock mass quality are obtained, and the slope of fitting straight line is f, and the intercept at σn axial is c based on Mohr-Coulomb criterion. The predominant slope method is obviously different from the least square method although they are all based on Mohr-Coulomb criterion and linear fitting method. For the predominant slope method, the upper- and lower-limit values of cohesion c should be conducted (Figure 5); the lower one is often suggested. Considering the obvious differences of Q values in Table 1, 1–6 sets are grouped as Class II while 7–12 sets being grouped as Class III in rock mass quality. Their scatter plots of normal stress σn versus shear stress τ are shown in Figure 5. The estimated results using the two methods are shown in Table 7.Table 7
The estimated values off and c using predominant slope and least square methods.
Class of rock mass quality
Set number
Predominant slope method
Least square method
f
c
f
c
II
1–6
1.35
2.00
1.38
3.16
III
7–12
1.07
1.50
0.89
4.19The scatter plots ofσn versus τ: (a) Class II and (b) Class III.
(a)
(b)As the prescribed standard [44] for the hydroelectric engineering in China, the suggested design parameters must have a guarantee rate β of more than 0.8. Therefore, it could be suggested that the β of more than 0.95 are too high, while those of less than 0.20 are too low. Relative to the 12 sets test results in Table 1, the guarantee rates βf and βc (by (17a) and (17b)) associated with the estimated f and c (in Table 7) are shown in Table 8. Most of them are not acceptable: all βc are too low using least square method, which means that the estimated c are overvalued; the βc are larger than 0.95 for 7–12 sets using predominant slope method, which means that the estimated c are too undervalued. The data points must be located around the fitted straight line as far as possible on account that the least square method focuses on the best linear fitting of normal stress σn and shear stress τ (Figure 5), which can cause the relatively low guarantee rates β defined as (17a) and (17b), while the data pairs of σn and τ are above the lower-limit fitting line (Figure 5) when the lower-limit values are adopted using the predominant slope method, which can lead to the large guarantee rates β. Therefore, the guarantee rate of the value c is low by the least square method but high by the predominant slope method.Table 8
The guarantee rates of estimatedf and c using the predominant slope and least square methods.
Set number
Q
Predominant slope method
Least square method
β
f
β
c
β
f
β
c
1
30.5
0.762
0.397
0.721
0.037
2
23.8
0.650
0.564
0.600
0.070
3
38.7
0.857
0.187
0.828
0.013
4
39.7
0.868
0.206
0.839
0.015
5
38.6
0.857
0.209
0.827
0.015
6
26.5
0.698
0.499
0.652
0.055
7
8.1
0.796
0.963
0.963
0.054
8
5.2
0.532
0.980
0.867
0.098
9
8.9
0.834
0.957
0.971
0.046
10
5.2
0.532
0.980
0.867
0.098
11
5.4
0.557
0.979
0.879
0.094
12
4.6
0.453
0.982
0.819
0.110


### 4.2. Estimated Shear Strength Parameters Using Copula-Based Method

The Copula-based fitting formulas established above are based on the test data of the 12 sets before; therefore, these tested shear strength parameters cannot be used to examine the extensional suitability of the Copula function obtained. The additional in situ direct shear tests are thus performed for the other 9 sets of the marbles with similar properties to check the application of Copula-based method to estimating shear strength parameters of rock mass. TheQ and GSI of the tested samples are also investigated at the same time. The test results are shown in Table 9.Table 9
The results of other 9 sets of in situ direct shear tests.
Set number
Rock mass quality indices
Shear strength parameters
Q
GSI
f
c/MPa
13
4.7
57.8
1.04
1.6
14
5.3
59.6
0.98
2.4
15
6.7
61.5
0.85
2.31
16
25.4
75.8
1.43
2.1
17
6.0
60.3
1.05
3.9
18
7.3
62.1
1.07
2.14
19
14.6
70.2
1.25
2.66
20
23.4
76.3
1.32
2.2
21
33.1
80.3
1.93
1.68The estimatedf and c using Copula-based method have the guarantee rate β of 0.8. The Hoek-Brown criterion (see (1) to (3), the disturbance factor D=0), as a comparative method, is also adopted to calculate shear strength parameters (in Figure 6). The indicator to assess the predicted suitability of Copula-based method and Hoek-Brown criterion is the root mean square errors (RMSE) obtained by
(18)RMSEf=1n∑i=1nfi-fi′2,RMSEc=1n∑i=1nci-ci′2,
where fi and ci are obtained from the in situ shear test and fi′ and ci′ are derived from Copula-based method or Hoek-Brown criterion. Clearly, the smaller the RMSE is, the more suitable the estimated result is.The shear strength parameters derived from Copula-based method and Heok-Brown criterion: (a) friction coefficientf and (b) cohesion c.
(a)
(b)Figure6 indicates that the estimated shear strength parameters by Copula-based method are closer to in situ test values than those by Hoek-Brown criterion for the rock mass specimens in Table 9, which means that the proposed Copula-based method would be a more effective tool in estimating shear strength parameters of the marbles involved in this paper. The Hoek-Brown criterion with a wide range of applications is empirically established on the basis of the test results of a large number of different types of rock mass; the calculated accuracies of f and c are however seriously affected by the uncertain test conditions of σ3n (see (1)) and undetermined disturbance factor D and so forth (see (1)–(3)) [45, 46]. However, the optimal Copula-based fitting functions established in the paper are on the basis of probabilistic analysis on the in situ direct shear test results of the 12 sets of similar marbles under the same test conditions, which should show better estimations of shear strength parameters of the rock mass with similar geological and mechanical properties.


## 4.1. Guarantee Rate

Instability would happen if the strength parameters of rock mass are overvalued, while economic waste could be caused if they are undervalued. Assuming that the estimated values of shear strength parameters aref1 and c1, the actual values are f and c; the guarantee rates βf and βc are then defined as the probability of f>f1 and c>c1, respectively, which can be expressed as
(16)βf=Pf>f1,βc=Pc>c1.The estimated value is overvalued if guarantee ratesβ are low, while being undervalued if high. In consequence, the estimated values can be evaluated by means of the guarantee rates β. On the other hand, the estimated f and c with a certain guarantee rates β can be calculated by solving the conditional probability of the optimal Copula functions. Based on the optimal Copula functions C of Q-f and Q-c established in Section 3.3, βf and βc can be calculated as(17a)βf=Pf>f1∣Q=Q0=1-Pf≤f1∣Q=Q0=1-CU≤u∣V=v=1-lim⁡Δv→0⁡Cu,v+Δv-Cu,vΔv=1-∂∂vCu,vv=v,(17b)βc=Pc>c1∣Q=Q0=1-Pc≤c1∣Q=Q0=1-C(U≤u∣V=v)=1-lim⁡Δv→0⁡Cu,v+Δv-Cu,vΔv=1-∂∂vCu,vv=v,where Q0 is the measured value of Q.For direct shear test, in general, the predominant slope method and the least square method [43] are used to statistically process the test data. Using the least square method, the best linear fitting relations of normal stress σn and shear stress τ for each class of rock mass quality are obtained, and the slope of fitting straight line is f, and the intercept at σn axial is c based on Mohr-Coulomb criterion. The predominant slope method is obviously different from the least square method although they are all based on Mohr-Coulomb criterion and linear fitting method. For the predominant slope method, the upper- and lower-limit values of cohesion c should be conducted (Figure 5); the lower one is often suggested. Considering the obvious differences of Q values in Table 1, 1–6 sets are grouped as Class II while 7–12 sets being grouped as Class III in rock mass quality. Their scatter plots of normal stress σn versus shear stress τ are shown in Figure 5. The estimated results using the two methods are shown in Table 7.Table 7
The estimated values off and c using predominant slope and least square methods.
Class of rock mass quality
Set number
Predominant slope method
Least square method
f
c
f
c
II
1–6
1.35
2.00
1.38
3.16
III
7–12
1.07
1.50
0.89
4.19The scatter plots ofσn versus τ: (a) Class II and (b) Class III.
(a)
(b)As the prescribed standard [44] for the hydroelectric engineering in China, the suggested design parameters must have a guarantee rate β of more than 0.8. Therefore, it could be suggested that the β of more than 0.95 are too high, while those of less than 0.20 are too low. Relative to the 12 sets test results in Table 1, the guarantee rates βf and βc (by (17a) and (17b)) associated with the estimated f and c (in Table 7) are shown in Table 8. Most of them are not acceptable: all βc are too low using least square method, which means that the estimated c are overvalued; the βc are larger than 0.95 for 7–12 sets using predominant slope method, which means that the estimated c are too undervalued. The data points must be located around the fitted straight line as far as possible on account that the least square method focuses on the best linear fitting of normal stress σn and shear stress τ (Figure 5), which can cause the relatively low guarantee rates β defined as (17a) and (17b), while the data pairs of σn and τ are above the lower-limit fitting line (Figure 5) when the lower-limit values are adopted using the predominant slope method, which can lead to the large guarantee rates β. Therefore, the guarantee rate of the value c is low by the least square method but high by the predominant slope method.Table 8
The guarantee rates of estimatedf and c using the predominant slope and least square methods.
Set number
Q
Predominant slope method
Least square method
β
f
β
c
β
f
β
c
1
30.5
0.762
0.397
0.721
0.037
2
23.8
0.650
0.564
0.600
0.070
3
38.7
0.857
0.187
0.828
0.013
4
39.7
0.868
0.206
0.839
0.015
5
38.6
0.857
0.209
0.827
0.015
6
26.5
0.698
0.499
0.652
0.055
7
8.1
0.796
0.963
0.963
0.054
8
5.2
0.532
0.980
0.867
0.098
9
8.9
0.834
0.957
0.971
0.046
10
5.2
0.532
0.980
0.867
0.098
11
5.4
0.557
0.979
0.879
0.094
12
4.6
0.453
0.982
0.819
0.110


## 4.2. Estimated Shear Strength Parameters Using Copula-Based Method

The Copula-based fitting formulas established above are based on the test data of the 12 sets before; therefore, these tested shear strength parameters cannot be used to examine the extensional suitability of the Copula function obtained. The additional in situ direct shear tests are thus performed for the other 9 sets of the marbles with similar properties to check the application of Copula-based method to estimating shear strength parameters of rock mass. TheQ and GSI of the tested samples are also investigated at the same time. The test results are shown in Table 9.Table 9
The results of other 9 sets of in situ direct shear tests.
Set number
Rock mass quality indices
Shear strength parameters
Q
GSI
f
c/MPa
13
4.7
57.8
1.04
1.6
14
5.3
59.6
0.98
2.4
15
6.7
61.5
0.85
2.31
16
25.4
75.8
1.43
2.1
17
6.0
60.3
1.05
3.9
18
7.3
62.1
1.07
2.14
19
14.6
70.2
1.25
2.66
20
23.4
76.3
1.32
2.2
21
33.1
80.3
1.93
1.68The estimatedf and c using Copula-based method have the guarantee rate β of 0.8. The Hoek-Brown criterion (see (1) to (3), the disturbance factor D=0), as a comparative method, is also adopted to calculate shear strength parameters (in Figure 6). The indicator to assess the predicted suitability of Copula-based method and Hoek-Brown criterion is the root mean square errors (RMSE) obtained by
(18)RMSEf=1n∑i=1nfi-fi′2,RMSEc=1n∑i=1nci-ci′2,
where fi and ci are obtained from the in situ shear test and fi′ and ci′ are derived from Copula-based method or Hoek-Brown criterion. Clearly, the smaller the RMSE is, the more suitable the estimated result is.The shear strength parameters derived from Copula-based method and Heok-Brown criterion: (a) friction coefficientf and (b) cohesion c.
(a)
(b)Figure6 indicates that the estimated shear strength parameters by Copula-based method are closer to in situ test values than those by Hoek-Brown criterion for the rock mass specimens in Table 9, which means that the proposed Copula-based method would be a more effective tool in estimating shear strength parameters of the marbles involved in this paper. The Hoek-Brown criterion with a wide range of applications is empirically established on the basis of the test results of a large number of different types of rock mass; the calculated accuracies of f and c are however seriously affected by the uncertain test conditions of σ3n (see (1)) and undetermined disturbance factor D and so forth (see (1)–(3)) [45, 46]. However, the optimal Copula-based fitting functions established in the paper are on the basis of probabilistic analysis on the in situ direct shear test results of the 12 sets of similar marbles under the same test conditions, which should show better estimations of shear strength parameters of the rock mass with similar geological and mechanical properties.


## 5. Conclusion

Using the proposed Copula-based method, the shear strength parameters with a givenβ can be estimated provided that Q has been obtained. The calculated guarantee rate β can be also used to analyze the suitability of the estimated values from the statistical methods for direct shear test data (such as predominant slope method and the least square method). For the marbles tested in this paper, compared with Hoek-Brown criterion, the Copula-based method is a more effective tool in estimating shear strength parameters.The paper only considers the relationships between two of the parameters. Actually, the Copula functions could be of multidimensions and can be used to evaluate the parameters related with multiple variables. Furthermore, the Copula functions can be used for reliability analysis of estimated rock mass parameters.



---

*Source: 693062-2014-11-04.xml*