Integrated parameter inversion analysis method of a CFRD based on multi-output support vector machines and the clonal selection algorithm

Computers and Geotechnics 47 (2013) 68–77

Contents lists available at SciVerse ScienceDirect

Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo

Integrated parameter inversion analysis method of a CFRD based on multi-output support vector machines and the clonal selection algorithm Dongjian Zheng, Lin Cheng ⇑, Tengfei Bao, Beibei Lv State Key Laboratory of Hydrology-Water Resources and Hydraulic Engineering, Hohai University, Nanjing 210098, China National Engineering Research Center of Water Resources Efficient Utilization and Engineering Safety, Hohai University, Nanjing 210098, China College of Water-Conservancy and Hydropower, Hohai University, Nanjing 210098, China

a r t i c l e

i n f o

Article history: Received 18 January 2012 Received in revised form 12 July 2012 Accepted 13 July 2012 Available online 10 August 2012 Keywords: Concrete-faced rockfill dam Multi-output support vector machines Clonal selection algorithm Integrated inversion analysis

a b s t r a c t A geotechnical parameter identification method based on multi-output support vector machines (M-SVMs) and the clonal selection algorithm (CSA) is proposed. Using this method, the rockfill material parameters of a concrete-faced rockfill dam (CFRD) are identified. Based on the Taguchi design, some possible combinations of material parameters are generated within the admissible ranges of material parameters. Then, using these combinations, the displacement of all the observation points of the CFRD is calculated using the finite element method (FEM). Next, different combinations of material parameters are used as input, and the calculated displacement is used as output to train some M-SVM models to simulate the complex relationship between the material parameters and the dam displacement. Integrated inversion analysis, which takes the static and creep properties of rockfill material into account simultaneously, is implemented to achieve a comprehensive understanding of the mechanical properties of the rockfill material. The optimization problem corresponding to the integrated inversion analysis is solved by the CSA, which has global convergence and is robust. During the process of searching for optimal material parameters, the dam displacement is calculated by the M-SVM mapping instead of the FEM, which may greatly reduce the computation time. Based on the observed settlement and finite element model of the CFRD, the inversion analysis method described above is implemented. The results indicate that the parameter identification method for rock-fill material proposed in this article is accurate, has a fast convergence rate and can be applied in practical engineering applications. Ó 2012 Elsevier Ltd. All rights reserved.

1. Introduction Concrete-faced rockfill dams (CFRDs) have become one of the most preferred dam types in hydraulic engineering because they are safe, easily constructed and economic. In many large hydraulic engineering projects [1,2], CFRDs have been adopted as water retaining structures. Therefore, it is of great importance to adopt reasonable and effective methods to predict the mechanical and hydraulic behavior of CFRD’s during normal operation period. In addition to direct analysis of the data collected from CFRDs [3], the inversion analysis method, which integrates monitoring data analysis with structural simulation, has become an important means to study the mechanical properties of CFRDs. The inversion analysis method of Geotechnical parameters includes analytical and numerical methods. The analytical solution of the inversion problem is often difficult to obtain and usually unknown in practice. The numerical geotechnical parameters identification method ⇑ Corresponding author at: Hohai University, Xikang road, Nanjing, Jiangsu Province 210098, China. E-mail address: [email protected] (L. Cheng). 0266-352X/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compgeo.2012.07.006

based on finite element method (FEM) and monitoring data analysis was first proposed by Kavanagh and Clough [4]. Numerical methods include the uncertainty inversion analysis method [5] and the optimization inversion analysis method. Among these methods, the optimization inversion analysis method is currently the most widely used method. Zhou et al. [2] used the genetic algorithm (GA) to identify the static and creep model parameters of rock-fill material, and Zhang et al. [6] adopted the complex method to identify the E-B model parameters of a CFRD. To reduce the inversion analysis calculation time, artificial neural networks (ANNs) can be used instead of FEM calculation to obtain the dam displacement. For example, Zhou et al. [7], Yu et al. [8] and Kang et al. [9] integrated the particle swarm optimization (PSO) algorithm, evolution algorithm (EA) and ant colony optimization (ACO), respectively, with the ANN to identify the material parameters of a dam. However, it should be noted that training an ANN often requires a large number of samples, and ANNs easily get stuck at a local minimum during the training process. Recently, a new nonlinear data modeling method has attracted attention in the machine-learning field: support vector machines (SVMs) [10–12]. SVMs are trained based on the structural risk

69

D. Zheng et al. / Computers and Geotechnics 47 (2013) 68–77

minimization (SRM) principle, which will make it giving a better generalization than ANNs. The training of SVMs is a uniquely solvable quadratic optimization problem, which will overcome the local minimization problems during the training process. SVMs have been successfully applied to many fields of geotechnical engineering [13–19]. For the geotechnical parameter identification problem, Wang et al. [16] and Zhao and Yin [17] integrated the SVM with the PSO, while Zhao and Feng [18] and Feng et al. [19] integrated the SVM with the GA. They all have successfully implemented material parameter inversion analysis with their integrated methods. However, the current SVM model that is most commonly used is a single-output model. In terms of structure inversion analysis, it is necessary to use the monitoring data from multiple observation points to achieve a more accurate understanding of the overall mechanical properties of a material. To process this multi-output problem, the current commonly used method is to establish a single-output SVM for each output independently. However, establishing an independent model for each output could not only greatly increase the computation time but may also increase the error because the correlation between different outputs is not taken into account. Considering the drawbacks of the single-output SVM, a type of multi-output support vector machine (M-SVM) model is used in this paper to simulate the complex relationship between monitoring data and rock-fill material parameters. To achieve a comprehensive understanding of the mechanical properties of the rock-fill material, integrated inversion analysis, which considers the static and creep properties of the material simultaneously, is performed. The clonal selection algorithm (CSA) is adopted to search the optimal material parameters that minimize the deviation between observed and calculated displacement. In the optimization process, the calculated displacement is obtained by the M-SVM mapping instead of the FEM calculation to reduce the computation time. Based on the observed settlement and the finite element model of a CFRD, the inversion analysis method above was implemented. The results demonstrate the good performance of this method. 2. Multi-output support vector machine M-SVM 2.1. Single output SVM model For a series of data (x1, y1), (x2, y2), . . . , (xl, yl) there is a nonlinear relationship between the input vector xi e Rr and the output data yi e R. According to the SRM principle, in the typical single-output SVM model, we try to find a function f(x) that transforms the nonlinear relationship between the input and the output into a linear mapping in feature space.

f ðxÞ ¼ wT uðxÞ þ b

ð1Þ

where u() is the nonlinear mapping that projects the input to the feature space; the vector w and the scalar b are parameters of the linear mapping. The function f(x) should not only ensure that there is only a small deviation e between the actual target yi and the predicted value from f(xi), but it should also be as flat as possible. This goal can be realized by the following optimization problem [10–12]:

Minimize : kwk2 =2 þ C

l X ðni þ ni Þ

rors exceeding the precision e. This problem can be rewritten as the following unconstrained optimization problem:

Minimize : LP ðw; bÞ ¼ kwk2 =2 þ C

l X Le ðui Þ

ð4Þ

i¼1

in which ui = yi  wTu(xi)  b and Le() is the Vapnik e-insensitive loss function

 Le ðuÞ ¼

0;

kuk1 < e

kuk1  e; kuk1 P e

ð5Þ

It is evident from Eq. (5) that for the single-output SVM, an insensitive area near the predicted value is defined by a L1 norm kk1. When the error between the predicted value f(xi) and the actual target yi exceeds e, the penalty factor C is imposed in the optimization objective function (4). 2.2. The M-SVM model Currently, the multi-output problem y e Rk is commonly solved by establishing k single-output SVM models independently. In the case of two outputs, the e-insensitive area is a square zone [20,21], as shown in Fig. 1, when the two output variables are modeled using the two single-output SVM models. In Fig. 1, the center point O (yi1, yi2) represents the actual target; point A and point B represent two predicted values. It is obvious that the predicted value of point A is better than that of point B. However, if two single-output SVM models are used and the loss function is defined as Eq. (5), the point A will be punished by the penalty factor C in the optimization objective function, while point B will not be, because the deviation between the predicted value A and the actual value O is over e in one direction and the predicted value of B is not. With the increase of output dimension k, this unreasonable situation would be even more serious. In addition, because each output variable is modeled independently, the possible correlation between outputs is not considered, namely, the information in the data is not fully used [22–24]. Meanwhile, for k output variables establishing k SVM models independently would greatly increase the computation time. A unified model is used to solve the multi-output problem. The loss function is a L2 norm kk as follows

( Le ðuÞ ¼

0;

kuk < e

ku  ek2 ; kuk P e

ð6Þ

pffiffiffiffiffiffiffiffi where u ¼ eT e; e = y  f(x); f(x) = WTu(x) + b; W = [w1, . . . , wk]; 1 b = [b , . . . , bk]T. When the loss function shown in Eq. (6) is used for the twooutput problem, the corresponding e-insensitive area is a round zone, as shown in Fig. 1. The unreasonable punishment imposed

ð2Þ

i¼1

8 > < yi  f ðxi Þ 6 e þ ni ; Subject to : f ðxi Þ  yi 6 e þ ni ; > : ni ; ni P 0;

i ¼ 1; 2; . . . ; l

ð3Þ

where kwk is the L2 norm of vector w; the constant C > 0 is a penalty factor; slack variables ni and ni are introduced to allow for some er-

Fig. 1. The e-insensitive area for the two-output problem (the area of the hypercubic and the hyper-spherical insensitive zone are of equal area).

70

D. Zheng et al. / Computers and Geotechnics 47 (2013) 68–77

by the optimal objective function on the prediction error would then be prevented successfully. At the same time, different output variables are uniformly limited by the same constraint condition, and it is not necessary to impose a constraint on each output variable. Thus, the calculation complexity is greatly reduced. Similar to the optimization problem (2) and (3) of the singleoutput SVM problem, the optimization problem of the M-SVM is defined as follows: k l X 1X Minimize : kwj k2 þ C ni 2 j¼1 i¼1

( Subject to :

ð7Þ

kyi  fðxi Þk2 6 e þ ni ; i ¼ 1; 2; . . . ; l ni P 0;

ð8Þ

i ¼ 1; 2; . . . ; l

This can also be rewritten as an unconstrained optimization problem as follows: k l X 1X Minimize : LP ðW; bÞ ¼ kwj k2 þ C Le ðui Þ 2 j¼1 i¼1

ð9Þ

To solve the optimization problem above, an iterative re-weighted least-square (IRWLS) algorithm is adopted. In the iterative process, to obtain the next solution (Wt+1, bt+1) from the previous solution (Wt, bt), LP(W, b) is approximated using the first-order Taylor expansion around (Wt, bt). For more details, refers to the literatures [20–22].

LP ðW; bÞ  L00P ðW; bÞ ¼

k l 1X 1X kwj k2 þ ai u2i þ s 2 j¼1 2 i¼1

ð10Þ

When the penalty factor C and the parameter of the kernel function, such as the parameter r2 of the Radial kernel function, are given, the IRWLS algorithm can be adopted to train the M-SVM and obtain the model parameters B = [b1, . . . , bk] and b = [b1, . . . , bk]T. Fig. 2 shows the flow chart of M-SVM learning and its basic steps are:  Step 1: Initialization: set t = 0, Bt = 0, and bt = 0, and then calculate the corresponding u0i and ai.  Step 2: Solve Eq. (14) to obtain Bs and bs; then, the searching descending direction Ps can be given by

" Ps ¼

#

Bs  Bt s

ð15Þ

t

ðb  b ÞT

 Step 3: Calculate the new Bt+1 and bt+1, i.e.,

"

#

Btþ1 tþ1 T

ðb

Þ

"

¼

#

Bt t

ðb ÞT

þ g Pt

g can be obtained using a back-tracking algorithm, i.e., initially, set g = 1.0, and use the following equation to calculate and determine whether LP(Bt+1, bt+1) < LP(Bt, bt) is right:

LP ðB; bÞ ¼

k l X 1X ðbj ÞT Kbj þ C Le ðui Þ 2 j¼1 i¼1



ui 6 e

0;

ð11Þ

2Cðui  eÞ=ui ; ui > e

Expression (10) shows a typical IRWLS problem. When (Wt, bt) has been obtained, the optimization problem corresponding to LP (W, b) can be transformed into searching the optimal solution of L00P ðW; bÞ. Based on the stationary point conditions @L00P j @b

"

@L00P ðW;bÞ @wj

¼ 0 and

¼ 0 we can obtain the following two expressions:

UT Da U IUT a

aT U

#"

wj

#

j

aT 1

b

" ¼

UT Da yj

aT yj

# ðj ¼ 1; 2; . . . ; kÞ

ð12Þ

where Da = diag(a1, . . . , al), U = [u(x1), . . . , u(xl)]T, a = [a1, . . . , al]T and yj = [yj1, . . . , yjl]. According to representer’s theorem [25], the machine learning problem can be expressed as a linear combination of the training samples in the feature space, i.e.,

wj ¼

l X bji uðxi Þ ¼ UT bj

ðj ¼ 1; 2; . . . ; kÞ

ð13Þ

i¼1

in which bj = [bj1, . . . , bjl] T. Placing this expression into Eq. (12), we can obtain the following expression:

"

K þ D1 a

1

aT K

1T

#"

bj b

j

#

" ¼

yj

aT y j

# ðj ¼ 1; 2; . . . ; kÞ

ð14Þ

(K)ij = k(xi, xj) = u(xi)Tu(xj) is the kernel matrix and should satisfy Mercer’s condition such that (K)ij corresponds to an inner product in the feature space. The following are commonly used kernel functions: (1) Polynomial kernel: k(x1, x2) = ((x1x2) + 1)d. (2) Radial kernel: k(x1, x2) = exp (|x1 - x2|/r2). (3) Sigmoid kernel: k(x1, x2) = tan h (w(x1x2) + h).

ð17Þ

If it is right, then go to next step 4; if not, set g = gf (0 < f < 1) and go back to the beginning of this step again.

where s is the sum of constant terms that do not depend on either W or b; ai is given by

ai ¼

ð16Þ

Fig. 2. The flow chart of M-SVM learning.

71

D. Zheng et al. / Computers and Geotechnics 47 (2013) 68–77

 Step 4: Calculate the new utþ1 and ai. If utþ1 < Tolerance is right i i or t has reached the maximum number of iterations tmax, then stop the calculation; otherwise, set t = t + 1, and go to step 2 again. 2.3. Optimize the model parameters of M-SVM During the training process of M-SVM, model parameters, including the parameter of the kernel function and the penalty factor C, should be optimized adaptively. Set a reasonable range for each model parameters initially, and the optimization on the model parameters is implemented with M-fold cross verification (CV) [26] and the grid-searching method. When M-fold CV is used, the l training samples are randomly split into M mutually exclusive subsets D1, D2, . . . , DM of approximately equal size lt = l/M. Select one of the M subsets as testing set and the other subsets as training set to train the M-SVM model. This training process is repeated M times until all the subsets have been selected as testing dataset once. For the jth training process, the root mean square error RMSEj is calculated by

RMSEj ¼

1X ky  fðxi Þk lt i2D i

ð18Þ

j

The average value

ARMSE ¼

M 1X RMSEj M j¼1

Fig. 3. The flow chart of the CSA.

ð19Þ

is used to evaluate the performance of model. Then within the ranges of model parameters, use some optimization method to search the best model parameters, which minimize the ARMSE. 3. The clonal selection algorithm The geotechnical parameter identification problem based on FEM and monitoring data is essentially equivalent to the following optimization problem [27]:

Minimize : gðxÞ Subject to : KðxÞu ¼ F

ð20Þ

x 2 Dx where g(x) is the optimization objective function, x = (x1, x2, . . . , xr) is the vector of parameters to be identified, K(x)u = F is the governing equation or mathematical model of the physical problem and Dx is the feasible region of parameters. Due to the complex behavior of rock-fill material, there are often many parameters to be identified in the constitutive model. Traditional optimization methods require significant computation time. Some biology inspired optimization methods, such as GA [2], PSO [8], and CSA [28–30], have good performance in this type of optimization problem. Among these, the CSA is a type of optimization method inspired by the intelligent behavior of biological immune systems. In this optimization method, the evolution chain of immune response mechanisms is abstracted to a mathematical optimization process. It has been verified that the CSA has global search capability, global convergence and strong robustness. The flow chart of the CSA is shown in Fig. 3. For a group of parameters (x1, x2, . . . , xr) to be optimized, the CSA has the following steps:  Step 1: initialization – Set t = 0, and the parameters, including antibody population size Np, best antibody size Nb, cloning factor c, radius of neighborhood d are initialized by the user empirically. The feasible domain Dx of parameters is given, and a initial population is randomly generated according to

the uniform distribution abi = (abi1, abi2, . . . , abir) = (xi1, xi2, . . . , xir), i = 1, 2, . . . , Np. Among these, the antibody abi is equivalent to a feasible solution of the optimization problem.  Step 2: selection – Select Nb best antibodies according to affinity. The affinity, which is similar to fitness in the GA, is defined as follows in the inversion analysis problem of this article:

aff ðabi Þ ¼

1 gðabi Þ

ð21Þ

From the above equation, it can be seen that for an antibody, as the value of the corresponding optimization objection function decreases, its affinity increases. The Nb antibodies with the highest affinity are selected to perform the cloning operation in the next step.  Step 3: cloning – The Nb selected antibodies are sorted in descending order according to their affinity. For each selected antibody, as the affinity increases, the cloning size increases. The cloning size Ni for a selected antibody i is given by

  cN p Ni ¼ round þ 0:5 i

ð22Þ

where round() is the rounding function, and the cloning factor c generally is equal to 0.2–0.3. Then, the total number of cloned antibodies Nc is

Nc ¼

Nb X

Ni

ð23Þ

i¼1

 Step 4: mutation – The purpose of the mutation is to realize the local optimization function of the CSA, and the Nc cloned antibodies above are mutated as follows:

Tðabij Þ ¼



abij þ ðrand  0:5Þd; rand < Pi abij ;

rand P Pi

ð24Þ

where abij is the jth component of the ith antibody; rand e [0, 1] is a uniformly distributed random number; and d is the radius of neighborhood initialized at the beginning. The mutation probability Pi of

72

D. Zheng et al. / Computers and Geotechnics 47 (2013) 68–77

the cloned antibody i should be inversely proportional to its affinity. The following expression is used in this article to calculate Pi

Pi ¼

affmax  aff ðabi Þ affmax  affmin

ð25Þ

where affmax and affmin are the maximum and minimum values of affinity, respectively, of the Nc cloned antibodies.  Step 5: reselection – The Nc cloned antibodies that have been mutated and the original Nb selected antibodies are re-selected to inhibit antibodies with low affinity and retain ones with high affinity in the new population. Simultaneously, the antibodies with high similarity are removed to increase the diversity of the antibody population and ensure the global search capability of this algorithm. The similarity can be evaluated according to the concentration of antibodies as follows:

 T s ðabi Þ ¼

1; densðabi Þ 6 T 0; densðabi Þ > T

ð26Þ

in which Ts() is the selection operator, and T is the antibody concentration threshold value, and the concentration could be calculated according to the following equation:

densðabi Þ ¼

NX b þN c

1 Nb þ Nc

aff ðabi ; abj Þ

ð27Þ

j¼1

where aff(abi, abj) is the similarity between two antibodies. The similarity can be calculated based on the Euclidean distance, Hamming distance and information entropy. In this article, the Euclidean distance based method is used:

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u r uX aff ðabi ; abj Þ ¼ t ðabih  abjh Þ2

ð28Þ

 

r3 is the internal friction angle; Pa is the atmoPa spheric pressure; and r1 and r3 are the maximum and minimum principle stresses, respectively. Under unloading and reloading conditions, the elastic modulus Eur is given by

u ¼ u0  Du log

Eur ¼ K ur Pa



r3

n ð30Þ

Pa

The bulk modulus Bt is expressed as follows:

 Bt ¼ K b Pa

r3

m ð31Þ

Pa

Parameters c, u, Du, K, Kur, n, Rf, Kb and m in the above equations can be determined by triaxial tests on the rock-fill material. Elastic matrix D can be calculated using the elastic modulus Et and the bulk modulus Bt as follows:

0

3Bt þ Et

B B 3Bt  Et B B B 3B  E t t 3Bt B B D¼ B 9Bt  Et B 0 B B B0 @ 0

3Bt  Et

3Bt  Et

0

0

3Bt þ Et

3Bt  Et

0

0

3Bt  Et

3Bt þ Et

0

0

0

0

Et

0

0

0

0

Et

0

0

0

0

0

1

C 0 C C C 0 C C C C 0 C C C 0 C A

ð32Þ

Et

After the elastic matrix D has been calculated, the overall stiffness matrix of structure is constructed, and the problem can be solved using a common FEM program. 4.2. Creep model

h¼1

in which abih and abjh are, respectively, the hth component of the ith and jth antibody, and r is the number of optimization variables.  Step 6: refreshing- After the series of operations above, the antibody population is refreshed. The total size of the antibody population remains Np, and the shortage is generated randomly. Each antibody in this new population will be tested to determine if it has reached the convergence limit g(abi) < Tolerance. If any antibody has reached the convergence limit or t has reached the maximum number of iterations tmax, the calculation stops; otherwise, set t = t + 1, and the process goes back to step 2. 4. Constitutive model of rock-fill material

Creep is an important mechanical property of the rock-fill material. Currently, some types of creep models, such as the exponential form model, Maxwell model, Burgers model and nine parameters model, are used to express the creep constitutive relationship of rock-fill material. In this paper, the Burgers model [32] is used. As shown in Fig. 4, the Burgers model is composed of a Maxwell model and a Kelvin model. In terms of the rock-fill material, the parameter EM shown in Fig. 4 corresponds to the elastic modulus of the E-B model calculated from Eq. (29) or Eq. (30). When the Burgers model is adopted to represent the creep constitutive of the rock-fill material, the creep strain increment DeV,t during the period from t0 to t = t0 + Dt is calculated as follows:

4.1. Static constitutive model

     EK 1 Dt DeV;t ¼ 1  exp  Dt Cr  eV;K; t 0 þ Cr EK gK gM

Because the mechanical properties of rock-fill are nonlinear, the paper uses the E-B model [31] to express this nonlinear property. Generally, the E-B model is expressed by an elastic modulus Et and a bulk modulus Bt. The elastic modulus Et is given by

Et ¼ KP a ð1  Rf SÞ2



r3

n

Pa

ð29Þ

where parameters K and Rf are the elastic modulus and damage ratio respectively. The shear stress ratio, or stress level, S is expressed as

S¼

The first item on the right side of Eq. (33) is the Kelvin creep strain increment, and the second item corresponds to the Maxwell model. eV,K,t0 is the Kelvin creep strain at time t0. EK, gK and gM are three parameters of the creep model that can be determined by field-testing. The matrix C is Poisson ratio matrix, and its expression can be found in Qian’s book [32]. When the FEM is adopted to solve the creep problem of the rock-fill material, an explicit program also called the initial strain method is used, and the creep strain is treated as a pseudo-load.

r1  r3 ðr1  r3 Þf

(r1–r3)f is the deviator stress when the tested material fails

ðr1  r3 Þf ¼

2c cos u þ 2r3 sin u 1  sin u

ð33Þ

Fig. 4. Burgers creep model.

73

D. Zheng et al. / Computers and Geotechnics 47 (2013) 68–77

( gðx1 ; x2 ; . . . ; xm Þ ¼

l X k 1X ðSVMij ðx1 ; x2 ; . . . ; xm Þ  dij Þ2 kl i¼1 j¼1

)0:5 ð35Þ

where SVMij(x1, x2, . . . , xm) and dij are, respectively, the calculated and the observed settlement at the ith measuring point at time tj. Now, a detailed step-by-step description of the integrated inversion analysis used to identify the parameters of rock-fill material will be given in the following. The flow chart of this inversion analysis method is shown in Fig. 6. Fig. 5. The sensitivity of dam settlement to E-B model parameters.

One advantage of this method is that the stiffness matrix does not need to be reassembled during the solution process. 4.3. Select inversion parameters In terms of the E-B model of rock-fill material, the parameters

u, Du, K, Kur, n, Rf, Kb and m must be identified. In addition, when considering dam material zoning, the number of parameters to be identified is very large. Therefore, it is necessary to analyze the parameter sensitivity and select several parameters with the greatest sensitivity to perform inversion analysis. Due to the large computational time, it is difficult to adopt a theoretical method to calculate the sensitivity. Morris’ method [33,34] is currently the most commonly used method. By perturbing each parameter independently, Morris’ method adopts the sensitive indicators s to reflect the sensitivity of the parameter:

s¼

q1 X ðY iþ1  Y i Þ=Y 0 i¼1

ðP iþ1  Pi Þ=P0

=ðq  1Þ

ð34Þ

where P0 represents the initial material parameters; Pi is the material parameters of the ith calculation; Yi is the calculated displacement of the ith time; Y0 is the calculated displacement using the initial parameters; and q is the total number of calculation. Given Eq. (34), a homogenous dam is used to calculate the sensitivity of the E-B model parameters in this article. Each parameter of the E-B model is perturbed independently and the number of FEM calculation is q = 6. The sensitivity distribution of dam settlement to all of the E-B model parameters u, Du, K, Kur, n, Rf, Kb and m is shown in Fig. 5. From Fig. 5, it can be seen that in terms of the E-B model, settlement has more sensitivity to the parameters K, u and Kb than to other parameters. Therefore, for the E-B model, we only perform inversion analysis on the parameters K, Kb and u. The three parameters EK, gK and gM of the Burgers creep model are all taken into account. Thus, for each type of rock-fill material, six parameters K, Kb, u, EK, gK and gM must be identified. 5. Integrated inversion analysis of rock-fill material Integrated inversion analysis takes the static and creep constitutive relation of rock-fill material into account simultaneously to achieve a comprehensive understanding of the mechanical properties of rock-fill materials. In addition, based on the results of the integrated inversion analysis, the future settlement of a dam can be predicted by the FEM calculation using the identified parameters. Then, the future state of a dam can be monitored by the predicted settlement values. As mentioned in Section 3, the inversion problem is essentially equivalent to the optimization problem defined by Eq. (20). For the integrated inversion analysis of rock-fill material, the optimization objective function is as follows:

(1) At first, for each material parameter to be identified, a reasonable range should be set according to field testing and engineering analogy and then some different values are selected uniformly within this range. Use all the selected values of all the parameters to generate some combinations of material parameters. As the total number of combinations maybe very large, when each possible combination of material parameters is used to perform FEM calculation, the total computation time will be too long to accept. One way to reduce the computational cost is to use the Taguchi design. In the Taguchi design, some Taguchi’s orthogonal arrays are provided to select some typical combinations of material parameters. Then, instead of using each possible combination of material parameters, only these selected combinations of material parameters will be used to perform FEM calculation. The Taguchi’s orthogonal arrays can be obtained from http://www.york.ac.uk/depts/maths/tables/orthogonal.htm. If it is necessary, some other different combinations of material parameters can be generated randomly to increase the robustness of samples. (2) For all the combinations of material parameters generated above, a non-linear FEM program that takes the static and creep properties of material into account is adopted to calculate the displacement of different monitoring points at time t1, t2,. . ., tp. At any time ti (i = 1, 2, . . . , p), different material parameter combinations and the calculated settlement of all the observation points using these parameter combinations form the training data set of M-SVM. (3) For the time ti (i = 1, 2, . . . , p), use M-fold CV to split the training data set into M mutually exclusive subsets and the MSVM model is trained M times using the method shown in Fig. 2 and the optimal model parameters of M-SVM are searched within a reasonable range to minimize ARMSE. Then using the optimal model parameters and all the training samples, the M-SVM is trained to simulate the relationship between material parameters and settlement at the time ti. For all the time t1, t2, . . . , tp, p M-SVM models with different model parameters are built to reflect the time-varying property of the dam settlement. (4) Based on the monitoring data of settlement from all the observation points, the CSA is adopted to search optimal material parameters within their admissible ranges. The optimal parameters are a combination of material parameters that minimize the objective function (35). During the process of searching optimal material parameters, the p trained M-SVM models above are used instead of the FEM to calculate the settlement of all the measurement points at different times t1, t2, . . . , tp. 6. Case study A comprehensive hydraulic construction is located in the Yellow river of China. As shown in Fig. 7, this construction is mainly composed of three parts: a CFRD, a waterpower generation system and discharge structures. The maximum height of the dam is

74

D. Zheng et al. / Computers and Geotechnics 47 (2013) 68–77

Fig. 6. The flow chart of the integrated inversion analysis.

Fig. 7. The plane layout.

132.2 m. The dead storage level, normal flood level and maximum flood level of the reservoir are 1008.0 m, 2002.0 m and 2005.0 m, respectively. The filling construction of the CFRD began in August 8, 2002, and the whole horizontal cross-section of it was filled simultaneously. The dam was filled to crest in only one stage without stage construction. On October 22, 2003, the dam was filled to the crest, and water storage began in August 8, 2004. The water level remained between 2001 m and 2005 m during the dam operation period. At the maximum cross section of the CFRD, a magnetic settlement gauge with 24 observation points was installed. These observation points were numbered EX2-1–EX2-24 from the bottom to the top of the dam. The layout of these observation points and the material zoning of dam is shown Fig. 8. The FE mesh of the maximum dam cross-section is shown in Fig. 9.

Let us now use the inversion analysis method proposed in this paper to determine the optimal material parameters for the CFRD under study. Only the material parameters of the rockfill material 3BIIand 3C are identified, and the material parameter values of 3BIwas set to be the same as that of 3C because the tested parameter values are similar. For other materials, the tested values of the material parameters were used. The tested values and the admissible ranges for all the material parameters to be identified are shown in Table 1. The observed settlement of 15 different days during the period between August 8, 2002 and November 27, 2005 was used. The 15 different days, include 3 typical days during the dam construction period and 12 typical days during the normal operation period. 25 combinations of material parameters were generated based on the orthogonal arrays (L25: six five-level factors) and 15 combinations were generated randomly. For each of the 40 combinations of material parameters, the settlement at 24 observation points of dam in the 15 different days was calculated using FEM. The FEM calculation was implemented using the commercial software MSC.Marc [35]. In order to make a comparison, for each day of the 15 days, an M-SVM model and 24 single-output SVM models were both trained. For each M-SVM and each single output SVM model, we have used a Radial kernel. The initial ranges of model parameters r2 and C were set as r2 e (210, 24) and C e (25, 213). The model parameters were optimized with 8-fold CV and the grid-searching method. The AMRSE calculated by Eq. (19) of the M-SVM and the 24 single-output SVM models is plotted in Fig. 10. As shown in Fig. 10, the M-SVM is more accurate than the 24 single-output SVM models. During the process of searching optimal material parameters, CSA and GA were both adopted to make a comparison. The MATLAB code of the CSA was developed based the flowchart shown in Fig. 3. The GA was implemented using a MATLAB toolbox GAOT

Fig. 8. Material zoning and layout of magnetic settlement gauge.

75

D. Zheng et al. / Computers and Geotechnics 47 (2013) 68–77

Fig. 9. The finite elements mesh of the maximum dam cross section.

Table 1 The tested values an admissible ranges of material parameters. Model

E-B

Parameters

3C Admissible ranges

Tested values

Admissible ranges

Kb

950 53.6 800

[800, 1800] [45, 60] [500, 1500]

850 49.6 285

[600, 1200] [45, 60] [200, 1000]

EK (Pa) gK (Pa s) gM (Pa s)

5e08 5e16 5e18

[1e8, 1e9] [1e16, 1e17] [1e18, 1e19]

5e08 5e16 5e18

[1e8, 1e9] [1e16, 1e17] [1e18, 1e19]

K

u (°) Burgers

3B II Tested values

Table 2 Inversion analysis results. Model

Parameters

3B II

3C

E-B

K Kb

1717.8 57.3 1172.1

881.5 57.8 485.9

EK (Pa) gK (Pa s) gM (Pa s)

7.39e8 6.79e16 2.72e18

5.33e8 2.16e16 2.31e18

u (°) Burgers

Fig. 10. Accuracy comparison between k single-output SVM and a M-SVM in terms of ARMSE.

Fig. 11. Convergence rate comparison between CSA and GA.

[36]. The convergence rate comparison between the CSA and the GA is shown in Fig. 11. It can be seen that the CSA has a faster convergence rate than the GA in this problem. The identified parameters of the E-B model and the Burgers model of two types of rock-fill materials are shown in Table 2. The settlement calculated by the M-SVM and the FEM using the identified parameters above and the observed settlement when dam filling is completed is shown in Fig. 12. The calculated and observed settlement of two observation points, EX2-10 and EX2-11, with the largest settlement, is shown in Fig. 13. It can be seen that the deviation between the settlement calculated by the FEM and by the M-SVM mapping is very small. Therefore, the system can meet the accuracy requirement when using the M-SVM instead of the FEM to calculate dam settlement. The calculated settlement using the optimal rock-fill material parameters and the observed

Fig. 12. Comparison between the calculated and the observed settlement in the day when dam filling was completed.

settlement have similar distributions along different elevations, show similar trends over time, and deviate little from each other. Therefore, the identified parameters of rock-fill material are reliable.

76

D. Zheng et al. / Computers and Geotechnics 47 (2013) 68–77

Fig. 13. Calculated and observed settlement of two observation points: (a) EX2-10 and (b) EX2-11 in 15 different dates since August 8, 2002.

Table 3 Predicted and observed settlement (unit: cm). Number

1 2 3 4 5 6 7 8 9 10

EX2-10 (elevation: 1930.11 m)

EX2-11 (elevation: 1936.06 m)

Observed

Predicted

Error

Observed

Predicted

Error

49.0 49.2 49.6 49.7 50.0 49.9 50.3 50.5 51.0 51.0

49.5 49.7 50.0 50.3 50.4 50.6 50.8 51.0 51.2 51.4

0.5 0.5 0.4 0.6 0.4 0.7 0.5 0.5 0.2 0.4

48.9 49.4 49.4 49.7 49.9 49.4 50.5 50.4 51.1 50.8

48.8 49.0 49.2 49.4 49.5 49.6 49.9 50.0 50.2 50.4

0.2 0.4 0.2 0.3 0.4 0.2 0.6 0.4 0.9 0.4

Fig. 14. Comparison between the predicted and the observed settlement of two measurement points: (a) EX2-10 and (b) EX2-11 since August 8, 2002.

The settlement of 10 different days during the period between January 10, 2006 and November 27, 2007 is predicted by the FEM calculation using the identified material parameters above. The comparison between the predicted and the observed settlement of two monitoring points EX2-10 and EX2-11with the maximum settlement is shown in Table 3. It can be seen that the prediction error is below 0.9 cm, which indicates good prediction accuracy. Assuming that the monitored settlement data are normally distributed, then the predicted settlement is within the interval ^i  3ri ; y ^i þ 3ri Þ with a probability of 99.74%, where ri and y ^i ðy are the standard deviation and the predicted settlement of measuring point i, respectively. As shown in Fig. 14, the observed settlement is within the prediction interval, which means that the settlement of the dam is normal. 7. Conclusions The M-SVM model adopted in this article has better performance than the traditional single-output SVM model; it is more accurate and reduces the computation time. The CSA has the advantage of a fast convergence rate and strong robustness when

it is used to search optimal material parameters. Based on MSVM and CSA, the integrated inversion analysis takes the static and creep properties of rock-fill material into account at the same time. Then, the identified material parameters can be used to calculate the future displacement of the dam using the FEM. In addition to the rock-fill material, the inversion analysis method proposed in this article can also be applied to other geotechnical problems. A further study, which will integrate the M-SVM with extended Bayesian method to identify geotechnical parameter, is being considered. Acknowledgements The authors are grateful to Fernando Pérez Cruz and his coworkers for making the MATLAB implementation of MSVR freely available at http://www.uv.es/gcamps/code/msvr.htm. This work was supported by the National Natural Science Foundation of China (Grant Nos. 51079086, 51139001), the Program for New Century Excellent Talents in University (Grant No. NCET-110628), the Special Fund of State Key Laboratory of China (Grant No. 2010585212), the Fundamental Research Funds for the Central Universities (Grant No. 2012B07214) and the Ministry of Water Re-

D. Zheng et al. / Computers and Geotechnics 47 (2013) 68–77

sources Public Welfare Industry Research Special Fund Project (Grant No. 201101013). References [1] Wilkins JK, Mitchell WR, Fitzpatrick MD, Liggins TB. The design of Cethana concrete face rockfill dam. In: Proceedings of annual engineering conference. Perth; 1973. p. 199–207. [2] Zhou Wei, Hua Junjie, Chang Xiaolin, Zhou Chuangbing. Settlement analysis of the Shuibuya concrete-face rockfill dam. Comput Geotech 2011;38(2):269–80. [3] Kim YS, Kim BT. Prediction of relative crest settlement of concrete-faced rockfill dams analyzed using an artificial neural network model. Comput Geotech 2008;35(3):313–22. [4] Kavanagh KT, Clough RW. Finite element applications in the characterization of elastic solids. Int J Solids Struct 1971;7(1):11–23. [5] Lee In-Mo, Kim Dong-Hyun. Parameter estimation using extended Bayesian method in tunneling. Comput Geotech 1999;24(2):109–24. [6] Zhang Lei, Zhang Ga, Wang Fuqiang, Zhang Jianmin. Rheological deformation back-analysis for Gongboxia concrete-faced rockfill dam. Rock Soil Mech 2011;32(Suppl. 2):521–5 [in Chinese]. [7] Zhou Wei, Ma Gang, Hu Chao. Long-term deformation control theory of high concrete-face rockfill dam and application. In: Proceedings of 2011 Asia-Pacific power and energy engineering conference. Wuhan, China; 2011. p. 1–4. [8] Yu Yuzhen, Zhang Bingyin, Yuan Huina. An intelligent displacement backanalysis method for earth-rockfill dams. Comput Geotech 2007;34(6):423–34. [9] Kang Fei, Li Junjie, Xu Qing. Ant colony clustering radial basis function network model for inverse analysis of rockfill dam. Chin J Geotech Eng 2009;28(Suppl. 2):3639–44 [in Chinese]. [10] Vapnik VN. An overview of statistic learning theory. IEEE Trans Neural Networks 1999;10(5):988–99. [11] Vapnik VN. The nature of statistic learning theory. New York: Springer; 1995. [12] Chapelle O. Training a support vector machine in the primal. Neural Comput 2007;19(5):1155–78. [13] Goh ATC, Goh SH. Support vector machines: their use in geotechnical engineering as illustrated using seismic liquefaction data. Comput Geotech 2007;34(5):410–21. [14] Samui P. Support vector machine applied to settlement of shallow foundations on cohesionless soils. Comput Geotech 2008;35(3):419–27. [15] Hongbo Zhao. Slope reliability analysis using a support vector machine. Comput Geotech 2008;35(3):459–67. [16] Wang Fuming, Li Xiaolong, Miao Li, Xu Ping. Mechanical parameters identification of surrounding rock based on wavelet SVM. J Hydroelectr Eng 2010;29(3):184–90 [in Chinese]. [17] Zhao Hongbo, Yin Shunde. Geotechnical parameters identification by particle swarm optimization and support vector machine. Appl Math Modell 2009;33(10):3997–4012. [18] Zhao Hongbo, Feng Xiating. Study on genetic-support vector machine in displacement back analysis. Chin J Geotech Eng 2003;22(10):1618–22 [in Chinese].

77

[19] Feng Xiating, Zhao Hongbo, Li Shaojun. A new displacement back analysis to identify mechanical geo-material parameters based on hybrid intelligent methodology. Int J Numer Anal Methods Geomech 2004;28(11):1141–65. [20] Pérez-Cruz F, Camps-Valls G, Soria-Olivas E, Pérez-Ruixo J, Figueiras-Vidal A, Artés-Rodríguez A. Multi-dimensional function approximation and regression estimation. In: Proceedings of artificial neural networks-ICNANN 2002. Madrid; 2002. p. 757–62. [21] Sánchez-Fernández M, de-Prado-Cumplido M, Arenas-García J, Pérez-Cruz F. SVM multiregression for nonlinear channel estimation in multiple-input multiple-output systems. IEEE Trans Signal Process 2004;52(8):2298–307. [22] Tuia D, Verrelst J, Alonso L, Pérez-Cruz F, Camps-Valls G. Multioutput support vector regression for remote sensing biophysical parameter estimation. IEEE Geosci Remote Sens Lett 2011;8(4):804–8. [23] Liu Guangcan, Lin Zhouchen, Yu Yong. Multi-output regression on the output manifold. Pattern Recognit 2009;42(11):2737–43. [24] Vazquez E, Walter E. Multi output support vector regression. In: Proceedings of 13th IFAC symposium on system identification. Rotterdam; 2003. p. 1820–25. [25] Schölkopf B, Smola A. Learning with kernels. Cambridge (MA): MIT Press; 2002. [26] Kohavi R. A study of cross-validation and bootstrap for accuracy estimation and model selection. In: Proceedings of the 14th international joint conference on artificial intelligence. Montreal; 1995. p. 1137–43. [27] Zhou Masheng, Li Yiqian, Xiang Zhihai, Swoboda G, Cen Zhangzhi. A modified extended Bayesian method for parameter estimation. Tsinghua Sci Technol 2007;12(5):546–53. [28] De Castro LN, Von Zuben FJ. The clonal selection algorithm with engineering applications. In: GECCO 2000 workshop proceedings of artificial immune systems and their applications. Las Vegas; 2000. p. 36–7. [29] Al-Sheshtawi KA, Abdul-Kader HM, Ismail NA. Artificial immune clonal selection algorithms: a comparative study of CLONALG, opt-IA, and BCA with numerical optimization problems. Int J Comput Sci Network Secur 2010;10(4):24–30. [30] Cutello V, Narzisi G, Nicosia G, Pavone M. Clonal selection algorithms: a comparative case study using elective mutation potentials. In: Proceedings of the 4th international conference on artificial immune systems. Banff, Canada; 2005. p. 13–28. [31] Duncan JM, Chang CY. Nonlinear analysis of stress and strain in soils. J Soil Mech Found Div, ASCE 1970;96(5):1629–53. [32] Qian Jiahuan, Yin Zongze. Theory and calculation of soil engineering. Beijing: China Water Power Process; 1996 [in Chinese]. [33] Morris MD. Factorial sampling plans for preliminary computational experiments. Technometrics 1991;33(2):161–74. [34] Tomlin AS. The use of global uncertainty methods for the evaluation of combustion mechanisms. Reliab Eng Syst Saf 2006;91(10–11):1219–31. [35] MSC. Software Corporation. MSC. Marc vol. A: theory and user information, Version 2005. Santa Ana, CA: MSC. Software Corporation; 2005. [36] Houck CR, Joines JA, Kay MG. A genetic algorithm for function optimization: a Matlab implementation.