Creep parameter inversion for high CFRDs based on improved BP neural network response surface method

The creep parameters of rockfill materials obtained from engineering analogy method or indoor tests often cannot accurately reflect the long-term deformation of high Concrete-Faced Rockfill Dams (CFRDs). This paper introduces an optimized inversion method based on multi-population genetic algorithm-improved BP neural network and response surface method (MPGA-BPNN RSM). The parameters used for inversion are determined by parameter sensitivity analysis based on the statistical orthogonal test method. MPGA-BPNN RSM, validated by root-mean-square error, mean absolute percentage error, squared correlation coefficient (R2), etc., completely reflects the response between the creep parameters and the settlement calculation values obtained by finite element method (FEM). MPGA optimized the objective function to obtain the optimal creep parameters. The results show that the settlement values of Xujixia CFRD calculated by FEM using the inversion parameters has great consistency with the monitored values both in size and in distribution, suggesting that the model parameters obtained by the introduced creep parameter inversion method are feasible and effective. The introduced method can improve the inversion efficiency and the prediction accuracy in FEM applications.


Introduction
High Concrete-Faced Rockfill Dams (CFRDs) have become one of the most popular dams in water conservancy projects due to the advantages of strong terrain capability, low cost, convenient construction, and short construction circle (Sukkarak et al. 2017;Xu et al. 2012). Accurately predicting the deformation caused by the creep effect of the rockfill by FEM can effectively reduce the void and cracking of the concrete panel (Zhang et al. 2015;Zhou et al. 2016). The creep parameters determined by engineering analogy (Huang et al. 2015) and indoor test with size effect (Shao et al. 2020;Zhou et al. 2019) often fail to reflect the actual long-term deformation characteristics of the target dam. Therefore, it is necessary for the inversion of model parameters based on the monitoring data.
The optimized inversion method with RSM aims to transform the parameter inversion problem into an objective function optimization problem, and then establish the mapping relationship between the parameters used for inversion and the FEM calculation values (called RSM, mainly including polynomial RSM and neural network RSM) to replace the iteration process in the FEM, the overall cost of the inversion analysis process can be greatly reduced (Guo et al. 2016;Yao et al. 2019;Yin et al. 2011). BPNN, as a kind of ANN, is an information processing mathematical model with a structure similar to the synaptic connections of human brain. In parameter inversion, if there are a sufficient number of combinations of parameter samples, the BPNN is able to predict the mapping result of any unknown parameter combination. Therefore, BPNN RSM can compensate for the limitations of low-order polynomial response surfaces in nonlinear finite element problems (Li et al. 2019;Ren and Bai 2011). Many researchers (Liu et al. 2020b;Togan 2013;Uzlu et al. 2014) have pointed out that the BPNN is sensitive to the initial weights and deviations, poor global search capability, and easy to fall into the local minimum. When predicting the results of parameter combinations using the trained BPNN, there is usually a significant discrepancy from the expected result. MPGA can solve the premature convergence problem of GA by co-evolving multiple populations with different control parameters, including the both global search and local search capabilities (Choubey and Kharat 2013;Pandey et al. 2014). Wang et al. (2017) predicted the deformation by MPGA-BPNN and various types of measured data (such as water level, temperature, aging, etc.). Yong et al. (2020) used MPGA-BPNN to classify the common faults of planetary gearbox. Zhang et al. (2021) used MPGA-BPNN to optimize the initial localization algorithm (ILA) by error correction.
MPGA-BPNN cited beforehand accomplished the purposes of machine learning and optimization algorithm. The further combination of MPGA-BPNN and FEM for parameter inversion is an unclear but worthy subject in the engineering field. This paper validates the rationality and accuracy of MPGA-BPNN RSM and uses it for parameter inversion based on the optimized inversion method. For the example of the Xujixia CFRD, Sect. 2 briefly introduces the monitoring system and creep model for rockfill materials. Section 3 obtains the sensitive parameters for displacement by Analysis of Variance (ANOVA) based on orthogonal test, which can further reduce the inversion cost. Section 4 introduces the specific implementation process of parameter inversion based on MPGA-BPNN RSM and MPGA optimization theory in detail. First, an objective function for dynamic inversion is determined. Second, three types of test samples are determined to reduce sample defects and train MPGA-BPNN to obtain the response surface. RMSE, MAPE, R 2 , etc., are used to completely validate the rationality of the combination of MPGA-BPNN and FEM. Third, the optimal creep calculation parameters are obtained by MPGA optimization, and the rationality and accuracy of the process are verified in the Xujixia CFRD. Sections 5 and 6 are discussions of some future research and conclusions, improving the efficiency and accuracy of parameter inversion.

Displacement monitoring system
In order to perform parameter inversion, the monitoring data must be processed and analyzed first. Five sets of hydraulic overflow settlement gauges (a total of 35 gauges) are used to monitor the internal settlement of Xujixia CFRD at five altitudes of 3392, 3395, 3425, 3430, and 3433 m at three observation sections of D0 ? 083.8, D0 ? 163.8, and D0 ? 223.8 m, as shown in Table 1 and Fig. 1. The parameters of the primary rockfill material mainly affect the settlement deformation of the primary rockfill, while the parameters of the secondary rockfill material mainly affect the settlement deformation of the secondary rockfill (Chi and Zhu 2016). Therefore, the settlement measured at the primary rockfill zone was used to invert the parameters of the primary rockfill material first. After determining the parameters of the primary rockfill material, the settlement measured at the secondary rockfill zone was then used to invert the parameters of the secondary rockfill material. In addition, in order to avoid the influence of different dam materials, the monitored values at the monitoring points in the central area of the primary and secondary rockfill were taken as the actually measured values (as shown in Fig. 1, red circles refer to the selected monitoring points in the primary rockfill zone and blue circles refer to the selected monitoring points in the secondary rockfill zone). Figure 2 shows the 3D FEM model of the Xujixia CFRD. Due to space limitation, we would only elaborate the parameter inversion of the primary rockfill zone.

Creep model for rockfill materials
Based on the dam prototype monitoring data, Shen (Shen and Zhao 1998) proposed a three-parameter model of rockfill creep based on the hysteresis deformation theory. The exponential curve of creep characteristic of rockfill is expressed as: where e t ð Þ is time-dependent creep strain, e f is the final creep deformation when t ? !, a is the ratio of initial creep deformation when t = 0, and exp is the base of natural logarithm.
Assuming the creep deformation of the rockfill is related to the confining pressure and stress level, the total creep deformation of the rockfill can be divided into volume creep e vf depending on the confining pressure and shear creep e sf depending on the stress level. According to the rockfill and clay creep deformation test results, Shen determined that the soil volumetric creep deformation and shear creep deformation could be simulated as: The creep deformation of the rockfill is closely related to the stress state of the rockfill. In the three-parameter creep model, e vf and e sf are, respectively, assumed to relate to the confining pressure and the stress level only, and the final volumetric creep was assumed to have a linear relationship with the confining pressure. However, a large number of experiments have proven that e vf is also related to the shear stress, and the volume creep has a nonlinear relationship with the confining pressure. In order to overcome the shortcomings of the three-parameter creep model, Li et al. (2004) proposed an improved seven-parameter creep model based on the three-parameter creep model, in which the final equations of volume creep and shear creep were revised as follows: The improved model can fully reflect the deformation characteristics of particle crushing and slippage of the rockfill under complex confining pressures, and the final creep value is related to the confining pressure, stress level and shear stress. This model contains seven parameters a, b, c, d, m 1 , m 2 , and m 3 . In this paper, the basis for the creep parameters of the primary and secondary rockfill material was determined by comparing the same project (Wen et al. 2017;Yao et al. 2019;Zhou et al. 2011). Then, the range of the parameters were also determined by changing each parameter value within a reasonable range. Table 2 shows the range of the primary and secondary rockfill creep parameters. Subsequently, orthogonal test was performed on the creep parameters in order to analyze the parameter sensitivity to the dam displacement.
3 Sensitivity analysis of creep model parameters based on the orthogonal test

Design of the orthogonal test table
The orthogonal test design is a highly efficient, fast and economical test design method based on orthogonality (Hu and Zhang 2019;Yang et al. 2020), which aims to select some uniformly dispersed, neat and comparable representative points from the overall test for the testing purpose. According to this, we obtained the orthogonal test table based on the sensitivity analysis of the creep model parameters in the following steps: 1. Select test indicators For the displacements of the dam in all the three directions, creep has a greater impact on settlement and horizontal displacement (Sukkarak et al. 2017). Therefore, the settlement and horizontal displacement of the monitoring points were chosen as the test indicators for parameter sensitivity analysis. 2. Determine the test factors and levels The seven parameters (a, b, c, d, m 1 , m 2 , and m 3 ) were taken as the test factors. Each factor was regulated up and down within a reasonable range in order to obtain the corresponding three test levels, as shown in Table 3 Table 4. In this table, the combination of factor levels corresponding to each row refers to one type of test plan. Then, FEM calculation was carried out according to the test designed in the orthogonal table, and the literature  shows the specific FEM calculation process.

Analysis of orthogonal test results
According to the FEM calculation results for all tests, variance analysis was performed on the design of the orthogonal test table (Li et al. 2020;Sun et al. 2016). Further, according to the comparison result of the variance square sum and the error square sum of each factor, F test was performed. Based on the quantitative results of the influence of various factors on the test indicators, we could judge whether the effect of each factor is significant. The process of variance analysis consists of the following three steps: 1. Calculate the total dispersion square sum S T , the dispersion square sum of each factor S A , and the dispersion square sum of tests error S E through the following equations, respectively: where f A and f E represent the degree of freedom of factor A and the total errors, respectively. 3. Calculate the statistics F A of the influence degree for the factor A on the test indicators through the following equation: Setting the significance level as 0.05, the sensitivity of a certain factor to the test indicators can be tested by comparing the value of F 0:05 in the F distribution table with the calculated value of F A , or by comparing the significance P A of the target factor with 0.05. F A [ F 0:05 or P A \0:05 indicates that the target factor has high sensitivity to the test indicator with a statistically significant effect.
By performing variance calculation and analysis on the designed orthogonal test table through Eq. (6-11), the sensitivity analysis result of each factor to the horizontal displacement U1 and settlement U2 of the dam was obtained, as shown in Tables 5 and 6, respectively. According to the judgment criterion F A [ F 0:05 2; 3 ð Þ ¼ 9:55 or P A \0:05, it can be concluded that at the monitoring point CS1-3-02, factors b and d have a significant impact on the test indicator U1 and factors a, b, c, d, m 1 , and m 3 have a significant impact on the test indicator U2; at the monitoring point CS1-3-03, factors b, d, and m 2 have a significant impact on the indicator U1, and factors a, b, c, d, and m 1 have a significant impact on the test indicator U2. Although the factors m 2 have a little impact on the test indicator U2, considering that m 2 is sensitive to the horizontal displacement, overall, all the seven parameters (a, b, c, d, m 1 , m 2 , and m 3 ) were chosen for parameter inversion. ANOVA based on orthogonal test determined the sensitive parameters, which can theoretically further reduce the cost of inversion (although not reflected in this case).

Objective function
In this paper, we established the mathematical expression of the objective function through the following equation by considering the disturbance of external factors and influence of internal factors on the monitoring data and the expected minimum value: where F X ð Þ is the objective function; X ¼ g are the parameters used for inversion; D is the number of parameters used for inversion (D = 7); a is the number of inversion periods; u is the number of external environmental factors affecting the target monitoring point (upstream slope disturbance, dam crest disturbance, downstream slope disturbance, etc.); W u ð Þ i is the weight of the external environmental factor u in the i-th ij is the FEM creep increment value of the parameter set X corresponding to the monitoring point j in the i-th time period; U real ij is the monitoring creep increment value of the monitoring point j in the i-th time period; Þare the lower limit and upper limit of each parameter x d . Considering the time periods of the existing monitoring data, i.e., the second pre-settlement period after the completion of the CFRD (the first 5 months) plus the second-phase slab construction period (the last 3 months), the value of a was calculated as a = 8. In this paper, for external factors, we only considered the impact of the upstream slope slab 5-month construction disturbance, i.e., u = 1, W For internal factors, we only considered the impact of the relative position K X ð Þ is the stiffness matrix; U f g is the nodal displacement array; R f g is the equivalent nodal load array.
The calculated creep incremental value U X ð Þ ij can be obtained by the FEM calculation. In this paper, the FEM calculation process described above, i.e., K X ð Þ U f g ¼ R f g, was replaced by the BPNN. We establish the mapping relationship between the creep parameters and the creep settlement by considering the input creep parameters in the FEM calculation as the input neurons of the BPNN and considering the creep settlement of the selected monitoring points in the FEM calculation results as the output neurons of the BPNN. While establishing the mapping relationship between the creep parameters of the primary rockfill material and the creep settlement at the monitoring points, a total of 64 BPNNs were established (the time period a = 8; the number of monitoring points b = 8). Each BPNN is expressed as net ij . Then, the creep increment can be calculated by generalizing BPNN as follows:

The MPGA-BPNN RSM
A typical BPNN is composed of the input layer, the hidden layer, the output layer, the weights and deviations, and the transfer function. Figure 3 illustrates the two-step calculation procedure of each neuron of the BPNN. Firstly, the initial results of the network are obtained by the forward calculation of randomly assigned weights and deviations. The calculation between the input and output layer of the BPNN can be expressed as below (Liu et al. 2020a where X and H are the number of neurons in the input layer and hidden layer, respectively; b h and b y are the deviations of the hidden layer and output layer, respectively; f hidden and f output are the transfer function of the hidden layer and output layer, respectively; w xh represents the weights between the input layer and hidden layer; w hy represents the weights between the hidden layer and output layer. Then, the calculated output value is compared with the actual output to find the difference between the two, which is referred to as the network output error. Based on the back-propagation learning algorithm, the weights and deviations of the BPNN are constantly revised by minimizing the error, until the network output error is less than the allowable error. The equation for calculating the network output error is as follows: whereŶ n and Y n represent the predicted output result obtained through network training and the actual result, respectively. In the BPNN, each of the input layer, hidden layer, and output layer has a different number of neurons. It has been found that the BPNN would have a better approximation effect with a smaller number of output neurons. Since the approximation effect of the BPNN directly affects the inversion result, the network described in this paper has only one output neuron, as shown in Fig. 4. In order to simulate the mapping relationship between the parameters and the settlement of multiple monitoring points and multiple time sequences in the FEM calculation, we adopted a multi-network system. All networks have the same input neurons, which are the parameters. The output neuron of each network is the settlement of one monitoring point in a certain period of time. The number of neurons in the hidden layer was determined by the following equation: where H hidden is the number of neurons in the hidden layer; H input is the number of neurons in the input layer; H output is the number of neurons in the output layer; h ran is a random number in the range of 1-15 (the value is determined by trial and error for deriving the optimal BPNN). On such basis, it was found that the best topology structure of the BPNN was 7-15-1, as shown in Fig. 4. At the same time, in order to allow the training network to converge as fast as possible, the sample data of the input layer and output layer was normalized to the range of [0.2-0.8]: where y is the normalized data; x is the original sample data; x max and x min are the maximum and minimum value of the original sample data, respectively.

Determination of the transfer function and backpropagation learning algorithm
The widely used TanSig function and Purelin function were chosen as the transfer function for the hidden layer and the output layer, respectively, as shown in Eqs. (18) and (19). The output value of the TanSig function is in the range of [-1, 1]. The TanSig function makes the BPNN nonlinear, which imposes a significant impact on the prediction accuracy. The Purelin function allows expanding the output result of BPNN. The back-propagation learning algorithm uses the Trainlm training algorithm, which is believed to be more robust than the Gauss-Newton method and have a higher convergence rate than the ordinary gradient descent method.

Update of weights and deviations
The training algorithm of the BPNN makes it prone to fall into local minimum, and often cannot obtain the global optimal weights and deviations. The SGA is featured with the advantages of strong robustness and global search capability. By combining SGA with the BPNN, it is possible to derive better weights and deviations. Below shows the calculation process of the SGA: 1. Perform individual binary coding on the potential solution of the problem to establish the mapping relationship between the phenotype and genotype, and then initialize the individual population randomly and perform decoding; 2. Evaluate the fitness of each individual, and then calculate the fitness value of each individual with the fitness function; 3. Perform genetic operations, including selection operation, crossover operation, and mutation operation, and then generate new populations.
However, due to the premature convergence problem of the SGA, a certain number of networks that do not meet expectations would be generated during a large amount of network training. In this paper, the parameter inversion of the primary rockfill material alone requires 64 networks to establish the mapping relationship between the creep parameters and the settlement. If there is any network involving the local minimum and premature convergence problem, the parameter inversion would not lead to satisfactory results. In our case, we applied MPGA to optimize the weights and the deviations of the BPNN. On the basis Fig. 4 Topological structure of the BPNN of SGA, the MPGA introduces the following concepts in order to overcome the premature convergence problem: 1. MPGA introduces multiple populations for simultaneous optimization and search. Different populations are assigned with different control parameters in order to achieve different search purposes; 2. The immigration operator can exchange information between various populations to realize the co-evolution of multiple populations. Thus, the acquisition of the optimal solution is the comprehensive result of the coevolution of multiple populations; 3. The artificial selection operator saves the optimal individuals in the evolutionary generations of various populations to ensure that the optimal individuals produced by various populations in the evolution process are not destroyed or lost, and serve as the basis for algorithm convergence. Figure 5 shows the schematic diagram of the algorithmic structure of MPGA-BPNN. The evolution mechanisms of population 1-N are all SGA. To sum up, Table 7 presents the structural parameter configuration of the MPGA-BPNN for building all the RSMs with the maximum efficiency.

Comparison of performance among BPNN RSM and MPGA-BP RSM
As mentioned earlier, while performing parameter inversion of the primary rockfill material, we need to create 64 NN RSM to establish the mapping relationship between the parameters and the settlement. For each NN RSM, we applied some currently popular methods to generate samples according to the value ranges of the parameters for minimizing the defects of insufficient samples. Specifically, there are 400 groups of training samples, including 18 groups of uniform test design U 18 3 7 ð Þ, 68 groups of orthogonal test design L 18 2 Â 3 7 ð Þ and L 50 5 7 ð Þ, and 314 groups of randomized test design. The test samples are composed of 100 groups of randomized test design. The sample parameters were taken as the input group to carry out FEM calculation so as to obtain the settlement as the corresponding output group. Further, the training and testing samples of the normalized data for each input and output group were generated. Then, by training the generated samples, the NN RSM was established.
In order to validate the performance of the NN RSM, we chose the root-mean-square error (RMSE), mean absolute percentage error (MAPE), and squared correlation coefficient (R 2 ) as the quantitative indicators, which can be calculated according to Eqs. (20)-(22). Besides, by comparing the calculated values obtained by the FEM with the predicted values obtained by the trained NN RSM using the same parameter sample. In terms of the four indicators above, we could judge the accuracy and robustness of the established BPNN RSM and MPGA-BPNN RSM.
whereŶ n is the predicted values obtained by the NN RSM for the n-th time; Y n is the calculated values obtained by the FEM. Figures 6, 7 and 8 evaluate the global accuracy, present the training and testing results of RMSE, MAPE, and R 2 based on the BPNN RSM and MPGA-BPNN RSM, respectively. In addition, Table 8 presents the average and range of the training and testing results of RMSE, MAPE, and R 2 , respectively. Figure 9 evaluates the case accuracy, shows the comparison graph and the error graph between the calculated values obtained by the FEM and the predicted values obtained by the trained NN RSM using the same parameter sample.
From the results above, we can draw the following conclusions: 1. As shown in Table 8, compared with the BPNN RSM, the RMSE and MAPE of the MPGA-BPNN RSM are 0.0153 and 0.0213, respectively, which are significantly reduced. The R 2 of the MPGA-BPNN RSM is around 0.95 both during the training and testing processes, which is about 13.84% higher than that of   suggesting that MPGA-BPNN can better find the optimal weights and deviations globally and has stronger robustness. 3. As shown in Fig. 9, based on the same parameters, the predicted values obtained by the MPGA-BPNN RSM show lower anomaly and are closer to the calculated values than that obtained by the BPNN RSM. Meanwhile, the accuracy and robustness of MPGA-BPNN are better than that of BPNN.
In summary, the MPGA-BPNN RSM can derive more accurate results, indicating that MPGA has solved the defects of poor global search ability and easy to fall into local minimum in BPNN. For the same set of parameters, the responses of MPGA-BPNN RSM show little abnormality, which shows that it can ensure both the efficiency and accuracy of the inversion when applied in the optimized inversion method. Therefore, it is a feasible strategy to retain the best performance MPGA-BPNN RSM and use it to parameter inversion based on the optimized inversion method.

Application of the MPGA algorithm for parameter inversion
It can be seen from the objective function that the parameter inversion problem can be transformed into a function optimization problem. The MPGA algorithm is able to overcome the premature convergence problem of SGA. On the basis of the monitored settlement data and the trained MPGA-BPNN RSM, we used MPGA to optimize the objective function. The parameters of MPGA were as follows: number of populations: 15; population size: 50; generation gap: 0.9; crossover probability: 0.7-0.9; mutation probability: 0.001-0.05; the minimum number of retaining generations of the optimal value: 10. Table 9 shows the optimal creep parameters of the primary and secondary rockfill material, respectively, obtained through inversion analysis. Subsequently, we performed the FEM calculation on the Xujixia CFRD using the creep parameters obtained from  Table 9. Figures 10  and 11 show the comparison between the calculated values and the monitored values of the primary and secondary rockfill material zone. Table 10 shows the relative error between the calculated values and the monitored values. The average standard deviations between the calculated values and the monitored values at the primary and secondary rockfill monitoring points are 1.82 and 2.32, respectively, and the average relative errors are 13.82% and 12.00%, respectively. The results above indicate that the settlement values calculated using the parameters obtained by parameter inversion have great consistency with the monitored values both in size and in distribution. Therefore, the inversion results can well-reflect the creep deformation characteristics of the rockfill material and can be used for fem prediction.
The method introduced in this paper has been well-applied in Xujixia CFRD. The combination of MPGA-BPNN and FEM and used to the optimized inversion method can greatly improve the efficiency and accuracy of inversion. According to the process of this method, if a reasonable objective function is constructed, it can be used in other geotechnical engineering FEM calculations.

Discussions
In this paper, we established the mathematical expression of the objective function by considering the disturbance of external factors and influence of internal factors on the monitoring data and the expected minimum value. How to reasonably assign weights to these factors has an important effect on obtaining the correct creep parameters. In our study, we only roughly considered the influence of the upstream slope faced construction disturbance. The weights of different internal and external factors would directly affect the optimization results. Therefore, future research needs to further address how to properly determine the weights of various internal and external factors.

Conclusions
This paper introduced an optimized inversion method for creep parameters with MPGA-BPNN RSM. ANOVA based on orthogonal test determined the sensitive parameters used for inversion, which can theoretically further reduce the inversion cost. MPGA-BPNN RSM established the relationship between creep parameters and FEM settlement increment, which can improve the accuracy of / 0 = Internal friction angle when the confining pressure is one atmosphere, D/ = Internal friction angle that changes with pressure, R f = Failure ratio, K = Tangent modulus coefficient, n = Tangent modulus index, K b = Volume Modulus coefficient, m = Bulk modulus index Fig. 10 Comparison between the calculated value and monitored value at the primary rockfill monitoring point inversion. MPGA was utilized to optimize the objective function to obtain the optimal creep parameters. The main conclusions are: 1. Compared with the BPNN RSM in terms of RMSE, MAPE, R 2 , etc., (1) The RMSE and MAPE of the MPGA-BPNN RSM are 0.0153 and 0.0213, respectively, suggesting a significant reduction; (2) The R 2 of the MPGA-BPNN RSM is around 0.95, which is about 13.84% higher than that of the BPNN RSM; (3) For the same parameter combination, the predicted values of the MPGA-BPNN RSM show a better consistency with the calculated values. The results show that MPGA-BPNN RSM has stronger robustness and accuracy, and is rational for parameter inversion. 2. By comparing the FEM calculation values using the optimal creep parameters with the monitored values, it is found that the average standard deviation of the primary and secondary rockfill of the Xujixia CFRD are 1.82 and 2.32, respectively, and the average relative errors are 13.82% and 12.00%, respectively, indicating that the settlement calculated by the parameter inversion method introduced in this paper is greatly consistent with the monitored data both in size and in distribution.
The parameter inversion method introduced in this paper has been well-applied on Xujixia CFRD, which can improve the efficiency and accuracy of parameter inversion, and can provide a reference for similar geotechnical engineering FEM calculation.