The Application of MsPSO in the Rockfill Parameter Inversion of CFRD

An intelligent algorithm that simultaneously analyzes multiple rockfill parameters is proposed. First, the paper introduces the operation and monitoring condition of the Shuibuya concrete-faced rockfill dam (CFRD). Then the constitutive rockfill models and the FEM analysis procedure are introduced in this paper. Third, the MsPSO intelligent algorithm was adopted to inverse the rockfill parameters. The recalculated displacement of Shuibuya CFRD using the inversed rockfill parameters is presented, and a satisfactory result was obtained, indicating that the inversion method is correct and effective. The method developed in this paper can be adopted in any geotechnical engineering parameter inversion.


Introduction
Concrete-faced rockfill dams (CFRDs) have rapidly developed in recent decades due to their strong adaptability to topography, geology, and climate and easily available construction materials.CFRDs were constructed with increasingly tall heights in recent years [1,2].Rockfill is used as the support structure and maintains the stability of CFRDs.Creep deformation accounts for the majority of rockfill dam postdeformations.Large rockfill creep deformation causes separation between the slabs and the cushion layer or even cracks in concrete slab [3], thus impairing the integrity of seepage control systems and weakening the structure or even threatening the safety of the dam [4].To avoid this potential risk in tall CFRDs, rockfill creep problem must be taken seriously; it is necessary to understand the deformation characteristics of rockfill to make the rockfill compatible with the concrete face in order to improve the design of tall CFRDs.
Currently, the common method for studying the character of rockfill creep references laboratory creep tests [5][6][7][8].At first, the laboratory creep tests were conducted under low confining pressure.Early in 1985, Parkin [5] presented rate methods and performed laboratory rockfill creep tests on the odometer.In China, the earliest rockfill creep test was completed with stress-controlled triaxial equipment by Shen [6] in 1991, and the exponent creep curve model was the first creep model type.As rockfill creep under low confining pressure laboratory tests cannot reflect the behavior of 200 m high CFRDs, Cheng and Ding [8] proposed a power function rockfill creep curve model under high confining pressure.In their study, many groups of laboratory triaxial creep tests with different confining pressures were performed; the creep deformation over time can be expressed by a power function.Until recently, the power function creep curve model has been applied widely for the creep calculation of tall CFRDs.
It is known that the in situ geomechanical properties of rockfill may differ from those determined by laboratory tests; one of the most obvious differences is called the "scale effect" [9,10].Therefore, to understand the in situ geomechanical properties of rockfill, it is better to perform back-analysis based on constructed CFRDs.Displacement inversion is an effective method to check and modify the parameters of rockfill models.Moreover, the back-analysis results are beneficial to the safety assessment of the project itself and to improving the design of subsequent structures of the same type.
Intelligent algorithm such as the artificial neural network (ANN) and the particle swarm optimization (PSO) algorithm were proposed to solve the geotechnical engineering inversion problem in recent years and achieved good results [11][12][13][14][15][16].Generally, the optimal values of the parameters to be determined are progressively approximated through iteration by minimizing the error function.Because of the complexity of problems in geotechnical engineering, numerical calculation is often used; the time-consuming finite element method (FEM) calculation is performed frequently in the inversion, so the rate of convergence is slow and sometimes the inversion fails due to nonlinear problems on a large scale.Furthermore, the result is affected by the initial values, and a local minimum or premature convergence is likely to be obtained; as a result, the solution is sometimes unstable.Zhou et al. reported the analyses of actual measured deformations resulting from continuous monitoring of the Shuibuya CFRD, which is currently the tallest CFRD in the world [14].The actual settlement records and laboratory test deformation parameters of the main and secondary rockfill were introduced.In this paper, a displacement inversion for rockfill parameters is performed using the multiswarm particle swarm optimization (MsPSO) algorithm [17].This work focused on vertical displacements, and the available monitoring data span nearly 5 years, covering the construction phase, the first filling stage of the reservoir, and 2 years of its operational time [14].

The Shuibuya CFRD
The Shuibuya CFRD was constructed in Badong City in Hubei Province, China, which is 117 km from Enshi City, located upstream, and 92 km from the Geheyan water-power plant downstream.The normal water level of the reservoir is 400 m, the maximum dam height is 233 m, the upstream dam slope is 1 : 1.4, and the downstream integrated dam slope is 1 : 1.4.The dam comprises five rockfill materials, including bedding material (IIA), transition material (IIB), primary rockfill (IIC), secondary rockfill (IID), and downstream rockfill (IIE).The boundary between the primary and secondary rockfill zones begins at an elevation of 380 m in the axes of the dam and ends at the downstream bottom of the dam with a slope of 1 : 0. The vertical displacements inside the dam body were measured by hydraulic overflow settlement gauges.They are distributed in three important cross sections: 0+132, 0+212, and 0+356.More specifically, a total of 11 monitoring lines were placed in the three cross sections, five in section 0+212 (at elevations of 235 m, 265 m, 300 m, 335 m, and 370 m) and three each in section 0+132 and section 0+356 (at elevations of 300 m, 335 m, and 370 m).The layout of the in situ instrument measurement scheme of the 0+212 cross section is shown in Figure 2 [14].

Displacement Calculation
The displacement of Shuibuya CFRD was calculated using the FEM.The instantaneous elastic-plastic deformation of rockfill materials was modeled using the Duncan - model [18], and the time-dependent creep deformation was modeled using the power function creep curve models [8].The tangent modulus  t can be expressed as follows: where  f is failure ratio,  i is initial deformation modulus, which expresses the compression-hardening property of rockfill materials, and  L is stress level.These parameters can be expressed as follows: where Pa is the standard atmospheric pressure,  is modulus number,  is exponent for defining the influence of the confining pressure on the initial modulus,  1 is cohesive strength,  is friction angle, and  1 ,  3 are major and minor principal stresses, respectively.For rockfill materials, the cohesive strength  1 is generally taken as zero.However, the angle of internal friction  changes with compression stress: where  0 and Δ are two constants for rockfill.
The bulk modulus  is given by where  b is bulk modulus number and  1 is bulk modulus exponent.
The - model parameters of the main and secondary rockfill of Shuibuya CFRD were introduced by Zhou et al. [14] and are shown in Table 1.[8] is employed in this paper.This model was proposed through the creep test results which showed that the correlation of the axial creep and time can be fitted by the power function in the double logarithmic:

The Power Function Creep Curve Model. The power function creep curve model in
where  1f is the limit of axial creep strain and  1 is a coefficient related to the axial creep strain rate.The time  is calculated in hours.The limit of axial creep strain  1f that combines both the confining pressure and the stress level viscous responses is given by the following relationship: where  3 is confining pressure,  1f is proportional to  3 , the coefficient  is the hyperbolic function of  L , and  2 and  are the creep parameters.
According to the creep test fitting curve between  1 under various  L and  3 ,  1 is independent of  L , and the relation between  1 and  3 can be expressed by the power function as follows: Thus, the axial creep model of rockfill can be described by ( 5)-( 9) and the creep parameters  2 , , , and  2 .
The relation of the volume creep strain  v () can also be expressed by the power function of time as follows: where  v () is the limit of volume creep strain and  v is a coefficient related to the volume creep strain rate.The limit of axial creep strain  v () is defined as where  v and  v are both functions of the stress level, defined by the following relationship: The limit of volume creep strain can be written as where  v is assumed as a constant:  v = const.Thus, the volume creep model of rockfill can be described by ( 10)-( 13) and the creep parameters   ,   ,   ,   , and  v .
The creep test parameters of the main and the secondary rockfill of Shuibuya CFRD are shown in Table 2.
According to ( 5) and ( 10), the creep strain ratio ε () has the following forms: The time history is divided into several increments Δ, and  1 and  v are obtained through summation: Therefore, the increments of shear creep and volume creep are given by the following, respectively: The relation among shear creep strain  s , axial creep strain  1 , and volume creep strain  v under ordinary conditions of triaxial test is given as follows: The above equation is also suitable for ε s (),  s (), and  sf .If the Prandtl-Reuss flow law is adopted on the  plane, the strain velocity tensor can be expressed as where [ ε ] is a strain velocity tensor, [] is a unit tensor, [] is the deviatoric stress tensor, and  is deviator stress, which is given by Then the creep increment tensor [Δ  ] can be written by So far, the numerical simulation method has been completely derived for CFRDs [15,19].The initial strain technique [20] is used to implement this creep model.The equivalent nodal force caused by rockfill creep can be calculated using the following equation: By loading the equivalent force on the corresponding nodes, the displacement caused by rockfill creep can be calculated.

The Finite Element Calculation.
In this paper the discretized mesh for the FEM analysis is composed of 8318 elements and 9138 nodes.
The interface between the slab and the cushion was simulated using Goodman contact elements, which were also applied to simulate the slab joints and the peripheral joints.The concrete face slab was simulated using a linear elastic model with an elastic modulus of 25.5 GPa and a Poisson ratio of 0.167 [21].The process of construction and water storage is divided into 30 loading steps based on the real construction schedule, as shown in Figure 1.Because the creep deformation of the rockfill occurs while the load is applied, the creep is analysis for the whole period of the CFRD construction and operation process.Basically, each loading step is followed by an equivalent loading step for creep deformation calculation.As a result, the creep process consisted of 29 steps.Thus, there are 59 steps used during FEM analysis, as shown in Figure 3.
This paper used the static finite element calculation program developed by Ganchen Gu et al. with Fortran software and added with independent development of a creep deformation equivalent loading calculation module.The use of the program to calculate dam creep deformation has been previously published [22].

Parameter Sensitivity Analysis
The variability of strength parameters , Δ and damage ratio   are small for many engineering experiences because many methods have been developed to calculate them [23] and the method is essentially perfect; therefore, their experimental values can be used directly.It is known that the main structure of CFRD is composed of the main and the secondary rockfill.It is necessary to inverse the two types of material parameters separately; without the three parameters listed above, there are (7+9−3)×2 = 26 parameters that need to be inversed.The calculation work of inversing all 26 parameters is extremely involved; as a result, it is necessary to analyze the sensitivity of these parameters and eliminate those parameters with low sensitivity.Sensitivity Indicators (SIs) of the parameters can be calculated using the modified Morris method [24] as follows: where   () represents the SI of the rockfill model parameters to the th measurement point's settlement,  is the quantity of calculation times, () is the rockfill model parameters value at the th calculation, (0) is the experimental value of the rockfill model parameters,   () represents the calculated settlement at the th measurement point at the th calculation, and   (0) represents the calculated settlement at the th measurement point using the experimental rockfill model parameters.
The measurement points at the elevation of 300 m were analyzed as an example here.The SIs of the main and secondary rockfill parameters are shown in Figure 4.
Several characteristics are shown by Figure 4.
(1) The SI of - model parameter is much larger than that of the creep model parameter, and there are significant quantitative differences between the two.As a result, the range of - model parameters value should be smaller, while the range of the creep model parameters value should be relatively larger.
(2) The experimental values of the parameters that have relatively lower SIs, such as  1 , ,  2 , and  v , can be used directly.
(3) The measurement points buried on the upstream side are more effected by the main rockfill parameters, as shown by the SI of the secondary rockfill parameters becoming smaller as the measurement points get closer to the upstream side.Thus, the measurement points buried on the far upstream side can be used to inverse the main rockfill model parameters, and the points buried on the far downstream side can be used to inverse the secondary rockfill model parameters; the remaining points can be used to verify the rationality of this method.

The Objective Function
The difference between the calculated and observed displacements at the measuring points could be taken as an objective function.To investigate the vertical displacements inside the dam body, the settlements of the measurement points in cross section 0+212 are analyzed here.As shown by Figure 2, the measurement points in cross section 0+212 were placed in 5 lines (at elevations of 370 m, 335 m, 300 m, 265 m, and 235 m).The weight of the measurement points at the same elevation needs to be averaged to reflect the effect of height.
Because the rockfill creep deformation is time-dependent, several time periods should be analyzed; in this paper, there are 9 time periods, according to Figure 1.In the absence of specific instructions, the weight of the measurement points at different time can be averaged.The mathematical expression of the objective function can be expressed as where  denotes time,  denotes elevation,   is the quantity of measuring points at the th elevation,  denotes a measurement point,  is a group of rockfill model parameters to be inversed,   () is the calculated value of displacement at th measuring point at time , and  *  is the corresponding observed value of displacement at th measuring point at time .
The optimization problem includes the objective function and the scope of variables; the variable here was , and its scope can be described as follows: where  is the quantity of the rockfill model parameters to be inversed and  max  and  min  are the up and down limits of the inversed parameters, respectively.
It should be mentioned that the settlement value of each measurement point at the beginning was not zero, and the meaning of this value was not clear.To overcome this problem during the process of inversion, the paper used the settlement difference between measurement points at different time periods.This is shown in the last part of the paper.

MsPSO Displacement Inversion Method
6.1.MsPSO Algorithm.The MsPSO algorithm that was proposed in [17] is a heterogeneous search approach based on four subswarms.Their definition and update rules are presented as follows.
(1) Define the basic subswarms  1 and  2 for exploitation search, in which particles implement the classic velocity and position update rules [25]: where th particle is represented by   = ( 1 ,  2 , . . .,   ) in -dimensional space and its velocity is   = ( 1 ,  2 , . . .,   ).Each particle will dynamically adjust its trajectory based on its historical best position (  ) and the best position ( best ) discovered by the whole swarm during the search process.The superscript "(1)/(2)" denotes the basic subswarms and  1 and  2 are both positive constants that are known as cognitive and social parameters, respectively. 1 and  2 are vectors of random numbers generated by the uniform distribution between 0 and 1.
(2) Define the adaptive subswarm  3 , in which the particles can adjust the flight directions adaptively by learning from better particles in the basic subswarms.The velocity of 1 P (2)   i P (2)   n P (3)   1 P (3)   i P (3)   n P (4)   1 P (4)   i P (4)   n P (1)   1 P (1)   i P (1)   n g (1)   best g (2)   best g best g (3)   best g (4)   best th particle, which employs the fitness values from the basic subswarms to update, is updated by the following equation: where  1 and  2 are the fitness values of the particles in the subswarm  1 and  2 , respectively. is the sum of the fitness value: that is,  =  1 +  2 .Similar to  1 and  2 ,  3 and  4 are also vectors of random numbers.The position update rule in  3 is the same as that in (26).
(3) Define the exploration subswarm  4 , in which the velocity and position update rules of th particle are defined as  (4)   ( + 1) =  (1)   ( + 1) +  (2)   ( + 1) −  (3)   ( + 1) , where  1 ,  2 , and  3 are called the impact factors and The cooperation and communication model for the subswarms is shown in Figure 5.In the graph the subswarms  3 and  4 use the information from the basic subswarms  1 and  2 to update.As can be seen, subswarm  4 is the only subswarm that shares all the useful information from the others.The four subswarms cooperate with each other directly or indirectly and all contribute to finding the global optimum  best .
The steps of inversion by MsPSO are as follows.
(1) Initialization: set the maximum number of the iterations (max) = 500 and the swarm size  = 40.The inertia weight  decreased linearly from 0.9 to 0.4, (2) Defining subswarms: specify the size /4 of each subswarm   | =1, It should be mentioned that the fitness of each particle depends on the objective function.Because the FEM calculated values cannot completely be consistent with the actual measured values, the value of  best will never reach 0. It is difficult to determine the value of the ; if the  = 0, the procedure will iterate for (max) = 500 times.However, as early as approximately 30 iterations, the optimal value no longer changes.Therefore, to save procedure time, the termination check of step (5) was changed to set a certain number of iterations.Once optimal value has been kept for the times, the procedure will stop, and the optimal value is the global optimum.In this case the maximum number and the  will no longer need to be set.
The particles in this study are the rockfill parameters, and their numbers and ranges should be determined before inversion.This paper combined the value of experimental parameters shown in Tables 1 and 2 and the sensitivity analysis to perform the work.According to (8), to avoid the denominator becoming negative, the value of parameter  cannot be too large.The number and range of the inversed parameters are shown in Table 3.

Procedures of the Method.
Because the FEM analysis needs half an hour to run once and because FEM analysis is performed frequently during the inversion process, considerable computational resources are needed to perform the inversion, and the overall time cost is unacceptable.Therefore, in this study, an RBF neural network was trained to establish the mapping relationship between the rockfill model parameters and the calculated displacement [22].
Figure 6 shows the flowchart of the displacement inversion method proposed in this study.

Results of the Inversion.
The values of the inversed rockfill model parameters and the corresponding objective values are shown in Table 4.
The recalculated displacement values of Shuibuya CFRD using the inversed rockfill model parameters are shown in Figures 7-9.
It should be mentioned that some of the measurement data are not ideal, such as the data at 3-1#.As 3-1# is closer to the upstream side compared with 3-2#, its settlement should be smaller in theory; however, the data show the opposite.This type of measurement data cannot be used for inversion.Among the measurement points in Figures 7-9, the measurement data at 3-2# and 3-3# were used to inverse the main rockfill model parameters, and the measurement data at 3-7# and 3-8# were used to inverse the secondary rockfill model parameters; the other data were used to verify the rationality of the method in this study.By comparing the results, it can be seen that the plots depict a similar trajectory of calculated settlements and measured settlements, indicating that the method of inversion in this study is correct and effective.
As shown in Table 4, the inversed values of - model parameters are approximately equal to their experimental values, while there is a certain gap between the inversed values of creep model parameters and their experimental values.
The creep model proposed by Cheng and Ding in [8] was based on the laboratory creep test, which lasted for approximately 1000 hours.According to (10) and the corresponding experimental values of rockfill creep model parameters, the following can be obtained.
For the main rockfill,  v (1000) =  vf (1 − 1000 −0.0678 ) = 0.374 vf . (29) Height of the dam Reservoir water level As a result, when the laboratory creep test treats the volume creep strain after 1000 hours as the limit volume creep, a contradiction occurred.In laboratory test, rockfill commonly takes a few hours to complete most of the creep deformation.However, the observed creep deformation from actual engineering lasts for years or even decades.The most important factor that causes the difference is restricting the size of rockfill samples in the laboratory test.Rockfill creep deformation is a process of particle rolling, slipping, rearranging, and the sharp angle crushing.All of these processes end very soon in laboratory tests; furthermore, because the size of the rockfill sample is limited, it is difficult to measure the post creep deformation.Therefore, the observed limit creep deformation is incomplete, which is the reason for the contradiction mentioned above.
This paper concludes that the major problem of laboratory rockfill creep test is the measurement of the limit creep deformation, which causes the relevant parameters to have problems during the calculation of the actual dam creep deformation.As shown in ( 8) and ( 13),  2 , ,   ,   ,   , and   are positively related to the limit creep deformation.Because the measured limit creep deformation by laboratory test is incomplete, the obtained value of these parameters will be smaller than their actual value.Thus, it is reasonable that the values of the inversed parameters are larger than the experimental values.

Conclusion
Some conclusions were obtained in this study.
(1) The MsPSO algorithm used in this paper has 36 evolution times, and the recalculated values of FEM using the inversed parameters are in good agreement with the actual measured values.Therefore, the algorithm has a good effect in terms of convergence accuracy and speed.The main reason for this result is that MsPSO with four related subswarms can search extensively and effectively, which overcomes the problem of falling into a local optimum in the search process.As long as those engineering problems can be simplified down to optimization problems, this method can be used efficiently and accurately.
(2) There is a certain gap between the inversed parameters and the experimental parameters.It was found that the measurement of the limit creep deformation is incomplete in laboratory tests, so the corresponding parameters do not apply to actual engineering.Thus, the method of obtaining rockfill creep model parameters through back-analysis according to the actual measured data appears to be more reasonable.The measured data The calculated data

2 .
The construction of the dam began in 2002 and was completed in 2007 and the filling of the reservoir was conducted in several sequential steps over the period from October 2006 to November 2008.The concrete dam face was constructed over three periods: first stage (January 2005 to March 2005), second stage (January 2006 to April 2006), and third stage (January 2007 to March 2007).The maximum cross section and the construction process of the Shuibuya CFRD are shown in Figure 1, where the solid lines represent the material partition and the dotted lines represent the construction partition.

Figure 2 :
Figure 2: Layout of in situ instrument measurement scheme in the maximum cross section of the Shuibuya CFRD.

Figure 3 :
Figure 3: The load steps in the FEM analysis.

Figure 4 :
Figure 4: SI of the parameters to the measurement points settlement at the elevation of 300 m.

Figure 6 :
Figure 6: Displacement inversion method based on an RBF neural network and MsPSO.

Figure 7 :
Figure 7: Cumulative settlements over time of the main rockfill area at an elevation of 300 m.

Figure 8 :
Figure 8: Cumulative settlement over time of the secondary rockfill area at an elevation of 300 m.

Table 1 :
The - model parameters of the Shuibuya CFRD materials.

Table 2 :
[14]meters of the power function creep curve model of the Shuibuya CFRD materials[14].

Table 3 :
Ranges of the parameters to be inversed.

Table 4 :
The values of the inversed parameters.
Figure 9: Cumulative settlement over time of the center point at elevations of 235 m, 265 m, and 335 m.