Abstract
This paper investigates the potential application of artificial neural networks in permanent deformation parameter identification for rockfill dams. Two kinds of neural network models, multilayer feedforward network (BP) and radial basis function (RBF) networks, are adopted to identify the parameters of seismic permanent deformation for Zipingpu Dam in China. The dynamic analysis is carried out by three-dimensional finite element method, and earthquake-induced permanent deformation is calculated by an equivalent nodal force method. Based on the sensitivity analysis of permanent deformation parameters, an objective function for network training is established by considering parameter sensitivity, which can improve the accuracy of parameter identification. By comparison, it is found that RBF outperforms the BP network in this problem. The proposed inverse analysis model for earth-rockfill dams can identify the seismic deformation parameters with just a small amount of sample designs, and much calculation time can be saved by this method.
1. Introduction
The dynamic response of rockfill dams under earthquake actions, mainly including absolution acceleration and permanent deformation, attracts more and more attention from engineers. The former is used to assess the dynamical load and seismic resistance of the dam. The latter is adopted to provide a basis for the dam height reserved during the design phase. So the prediction of permanent deformation is an essential problem in the seismic design for rockfill dams. As a result, it is important to select model parameters of rockfill dams reasonably, which makes for improving the accuracy of the numerical calculation.
However, it is not easy to carry out. Because of the difference of construction technology and construction quality and so on, the spatial distribution of material properties is considerably random in each project. Along with the development of construction technology, the maximum particle size of materials of rockfill dams may be bigger and bigger, but the model parameters are usually measured in the laboratory by using samples with much smaller size. The prepared experimental samples in the laboratory are different from the construction conditions. Therefore, the mechanical properties of samples determined in the laboratory may be more or less differing from those in situ. And then the stress and deformation acquired with laboratory parameters deviate far from the actual values inevitably. As a consequence, it is necessary to take measures to get model parameters in accord with the results of dam observation and make an accurate evaluation of dam safety and stability after that. Displacement backanalysis is an effective method to check and modify the parameters of soils.
In recent years, reverse analysis is mainly based on two methodologies: optimization algorithms and neural networks [1–5]. There are three types of optimization algorithms that have been frequently used in inverse analysis. The first type is gradient-based direct search algorithms, such as Levenberg-Marquardt method. The second type is the relatively simple direct search methods making no use of gradient, such as the simplex search method. The last type is kinds of intelligent global search algorithms, such as genetic algorithms, differential evolution, particle swarm optimization, and ant colony optimization. The first and the second type algorithms both have an advantage of estimating the solutions in relatively short computational time, but the results are affected by the initial values, and premature convergence is likely to happen. As an alternative to the direct search algorithms, intelligent global search algorithms are being widely adopted in reverse analysis, but they have a disadvantage of being time consuming.
In the geotechnical engineering field, intelligent backanalysis methods based on artificial neural networks (ANNs) and genetic algorithms [6, 7] are often adopted. As for generic algorithms, the range and trial values of the undetermined parameters should be given before the analysis, and then the time-consuming finite element method (FEM) calculation is performed repeatedly, so it is hard to ideally solve complicated nonlinear problems with a lot of finite elements. That is why it has been primarily used for seeking answers to static problems and two-dimensional problems so far. They need relatively few iterations and finite elements. Comparatively speaking, the strong nonlinear relationship between the known and unknown quantity in geotechnical engineering can be mapped well by using ANNs. And neural network approach can obtain inversion parameters quickly with just a small amount of sample designs.
The aim of this paper is to present an inverse analysis model for seismic permanent deformation parameters of earth-rockfill dams based on artificial neural networks. Section 2 presents the theories of forward computational models for static analysis, dynamic response analysis, permanent deformation analysis, and sensitivity analysis of design parameters. Section 3 introduces backprorogation neural networks (BPNNs) and radial basis function neural networks (RBFNNs). Section 4 shows the performance of both ANNs. Finally, conclusions are made in Section 5.
2. The Mathematical Model
2.1. Static Analysis
Duncan and Chang’s E-B model [8] is used to simulate the mechanical behavior of the rockfill materials. It is a nonlinear elastic model with wide application and is characterized by seven parameters: cohesion , friction (or , ), failure ration , modulus number , modulus exponent , bulk modulus number , and bulk modulus exponent . The nonlinear stress-strain relation of rockfill is represented by a hyperbola, whose instantaneous slope is the tangent modulus . According to the conventional triaxial tests, can be expressed as follows: The E-B model follows the Mohr-Coulomb criterion, and the wider the range of pressure involved the greater the curvature of the Mohr-Coulomb envelopes, since friction becomes smaller with the increase of minor principal stress . So as to coarse-grained soil, friction is not constant everywhere in dams. This variation in property may be represented by an equation of the form: where is the value of for and is the reduction of for a 10-fold increase in . And the bulk modulus can be expressed as
2.2. Dynamic Analysis
Equivalent linear elastic model [9] is used to simulate the dynamic properties of the earth-rock mixtures with two basic characteristics: nonlinearity and hysteresis. Soils have been deemed to be viscoelasticity in the model, and the dynamic stress-strain relationship is reflected with equivalent shear modulus and equivalent damping ratio . The key of the model is to determine the relation between maximum dynamic shear modulus and mean effective principal stress , as well as the variation of and along with the amplitude of dynamic shear strain. Based on a large number of experiments conducted by China Institute of Water Resources and Hydropower Research (IWHR), maximum dynamic shear modulus can be expressed as where are modulus and exponent usually determined by experiments. is the atmospheric pressure. In the experiments, the three-axis instrument and wave velocity test device were used. In regard to and , a relation curve is achieved firstly through experiments, which describes the variation of the dynamic shear modulus ratio / and damping ratio with dynamic shear strain . Then is normalized by reference shear strain to make instantaneous and easy to get through interpolation.
The effects of hydrodynamic pressure have to be considered when analyzing the dynamic interaction between dam and water in reservoir, since reservoirs may not always operate at a low water level during an earthquake. An ideal way to consider the effects is to divide finite element grids regarding the water and dam body as a whole, and interface element is utilized between contacting surfaces. However, it requires the computer with sufficient memory for the process and is very time consuming. Besides, stiffness coefficient of the interface element is hard to determine. So the additional mass method has been widely used so far. The effects of dynamic water pressure on seismic response of the dam are taken into account with the equivalent additional mass, and the dynamic analysis is done by adding the equivalent mass and the mass of the dam itself.
In this paper, equivalent additional mass is calculated by the lumped-mass method. The equivalent additional mass focusing on node is calculated by simplified Westergaard formula [10]: where is the angle between the upstream slop and the horizontal plane is depth of water from the surface to the bottom of reservoirs. is depth of water from the surface to node , and is the corresponding control area of node .
2.3. Permanent Deformation Analysis
2.3.1. The Model of Permanent Deformation
There have been several models for permanent deformation calculation at present, such as IWHR model [11], Debouchure model, Shen Zhujiang model [12], and improved models about them [13, 14], among which Shen Zhujiang model owns the broadest application so that the model is expressed in an incremental form, and permanent deformations in various conditions including different vibrations, dynamic shear strain, and stress levels can be calculated with only a set of parameters. Compared with the other models, Shen Zhujiang model not only has a consideration on both shearing deformation and volume deformation, but also is easier to use. Residual volumetric strain and residual shear strain can be written as where is dynamic shear strain amplitude, is stress level, is the number of vibrations, and , , , , and are experimental parameters. The parameter is assumed to have nothing to do with , that is, .
However, studies in recent years have suggested that the deformations calculated by Shen Zhujiang model are larger than actual performance, and it is adverse to an accurate evaluation of the seismic behavior of faced rockfill dams. So it is necessary to appropriately improve the model. Zou Degao focused on the influence of stress level on the residual shear deformation and presented an improved model based on a large number of cycle triaxial experiments [13]. When the earthquake-induced permanent deformations are calculated with FEM, the improved model can be expressed as an incremental form as follows: where , are increments of residual volumetric strain and residual shear strain, respectively.
2.3.2. The Method of Permanent Deformation Analysis
The major ways of overall seismic deformation analysis are based on the method of equivalent nodal force and modulus soften model [15] now. The calculation process of modulus soften model is relatively complicated, so the method of equivalent nodal force is a better choice for the permanent deformation analysis; the ideal about it is that the residual strain during an earthquake is determined firstly by an empirical formula, then the residual strain is converted to equivalent node force of finite elements, and the contributions of residual strain to the dam are replaced by the displacement calculated with the equivalent nodal force. The procedure comprises the following three steps.(1) Perform static calculation for the dam with midpoint incremental method, to get the initial static stress and stress level .(2) Perform dynamical calculation through the approach on the basis of equivalent linear viscoelastic model, to get the dynamic stress of soil, then convert the dynamic stress to stress state in laboratory, and get residual strain potential , according to (2.7). And strain potential of finite elements is obtained according to the following formula: where is generalized shear stress, is average principal stress, and is increment of residual strain in Cartesian coordinates.(3) Calculate the equivalent nodal force with the converted strain potential according to the formula that where is the conversion matrix of strain is the elastic matrix. Then permanent deformations are calculated with the equivalent nodal force applied on finite elements.
3. Permanent Deformation Parameters Inversion Using Artificial Neural Networks
3.1. BP Neural Networks
For a BPNN with a structure , that is, a vector input , hidden units, and an output vector as in Figure 1, the equation that expresses the relationship between the input and output can be written as where is the number of input units, is the number of neurons in the hidden layer, is the th input unit, is the weight parameter between input and hidden neuron , is the weight parameter between hidden neuron and output neuron , is the activation function of the hidden layer, and is the transfer function of output layer.
The weights were estimated and adjusted in the learning process with an aim of minimizing an error function as follows: where is the number of input and output examples of the training dataset and is the target value. The errors were fed backward through the network to adjust the weights until the error was acceptable for the network model. Once the ANN is satisfied in the training process, the synaptic weights will be saved and then used to predict the outcome for the new data. To minimize , optimal parameters of weights and biases have to be determined. One of the algorithms for solving this problem is the Levenberg-Marquardt (LM) algorithm. This algorithm is a modification of the Newton algorithm for finding optimal solutions to a minimization problem. The weights of an LMNN are calculated using the following equation: where is the Jacobian matrix of output errors, is the identity matrix, and is an adaptive parameter. When , it becomes the Gauss-Newton method using the approximate Hessian matrix. If is large enough, the LM algorithm becomes a gradient descent with a small step size (the same as in the standard backpropagation algorithm).
3.2. RBF Neural Networks
RBFNN is a kind of feedforward neural network and generally consists of three components: input layer, hidden layer, and output layer. Figure 2 displays an RBF network with a structure that there are inputs, hidden units and outputs. is the input vector. is the weight matrix of inputs and are offset of output units. is the output vector of the network, and is the activation function of the hidden unit ; one common function is the Gaussian function: where is spread constant, the role of which is to adjust sensitivity of the Gaussian function.
The final output of unit can be expressed as where is the weight of output layer and is spread constant of the base function.
3.3. Sample Set Designs
The prediction accuracy of neural networks is determined by sample quality to some extent, the samples used for inversion have to be accurate and balanced, so as to be representative enough, and it is better to reflect inner characters of the model system. In this paper, the methods of central composite design [16] and orthogonal design are utilized to design samples. The central composite design was presented by Box and Wilson. It consists of the following three parts: a full factorial or fractional factorial design, a central point, and an additional design, often a star design in which experimental points are at a distance from its center. Figure 3 illustrates the full central composite design for optimization of three variables. Full uniformly rotatable central composite designs present the following characteristics:(1)require an experiment number according to , where is the factor number and is the replicate number of the central point.(2)-values depend on the number of variables and can be calculated by . For two, three, and four variables, it equals 1.41, 1.68, and 2.00, respectively.(3)All factors are studied in five levels ).
As a result, samples designed by both methods embody not only inner but also outer conditions of the cube within certain realms, and they have good density distribution and representativeness.
3.4. Parameter Sensitivity Analysis
Sensitivity analysis can estimate the influence of parameter variation on outputs and make us have a more intuitive understanding about the parameters to be considered. Morris method [17] was more applied in sensitivity analysis. It figured out the influence of arguments on the dependent variable through disturbing a parameter and keeping the others the same. The corrected Morris method changed arguments at a fixed step size, calculated the value of influence in each condition, and then took the average like this: where is the sensitivity factor, is the output of condition , is the output derived from the initial parameters, is the percentage of parameter variations in condition compared to the initial parameters, and is the number of conditions.
4. Example Analysis
4.1. Brief Introduction to the Project
The Zipingpu dam is located in a valley, 9 km away in northwest from the Chengdu City, Sichuan province. It is one of the high CFRDs more than 150 m in China, with a maximum height of 156 m and the crest length 663.77 m. It encountered the Sichuan 8.0-magnitude earthquake which was higher than its actual design level. The dam body emerged obvious damage, and it provided rich and precious materials [18, 19] for earthquake engineering research on CFRDs.
4.2. Results of Static Calculation
Static parameters were directly selected from the experimental results coming from the IWHR, considering as well the results from Professor Zhu Cheng who partly backanalyzed the parameters of the project by immune genetic algorithm. Integrating both results, the parameters of rockfill were determined in Table 1. Besides, the linear elastic model was adopted for the calculation of concrete panels, concrete strength grade was 25, and the corresponding material parameters were density , MPa, Poisson’s ratio .
To simulate the process of construction and impoundment of the CFRD, midpoint incremental method and multistage loading process were used in the calculations. The dam was meshed into 6772 finite elements with total 6846 nodes, as shown in Figure 4. And the main results of a typical cross-section in rockfill zone were shown as in Figures 5 and 6, which were at operational water level before the earthquake.
4.3. Results of Dynamic Calculation
Due to influence of all kinds of factors, Zipingpu dam had no acceleration recordings of the actual principle shock of base rock. According to analysis [20, 21] of many scholars, the actual input ground motion was likely to be more than 0.5 g. Referring to the attenuation relationship [22] given by Yu et al., and considering the hanging wall and footwall effects, Professor Zhu Cheng deduced that the peak accelerations of dam site in three directions, respectively, were 0.52 g in east-west, 0.46 g in north-south, and 0.43 g in vertical. And based on the materials concerned in [22], the relative position of seismometer stations and dam site was shown in Figure 7. Finally, the seismograms of Wolong Station (051WCW) closer to the dam site were selected as the ground motion input; meanwhile the acceleration time histories were scaled in accordance with the peak accelerations as aforementioned, as shown in Figure 8.
According to the calculation results, the basic frequency of dam vibration was about 1.65 Hz and the maximum acceleration responses at dam crest were 0.86 g along the river, 0.74 g in vertical, and 1.36 g along the dam axial, which were consistent with the analysis results obtained by Kong et al. [21]. Under the three-dimensional earthquake, the maximum acceleration response lay in the downstream dam crest, and rockfill slid when its acceleration response exceeded yield acceleration, which qualitatively explained the phenomenon on Zipingpu dam that some grains in the downstream dam crest loosened and tumbled.
4.4. Backanalysis of the Project
Though the dam has different material zoning, material sources are same the and the construction control indexes are very near; therefore dynamic properties of different zones will not vary much. To reduce the workload, all rockfills were thought to have the same set of permanent deformation parameters, they were but was out of the inversion range for that it was a constant as mentioned in (2.6). Reference [23] listed several project cases at home and abroad and gave the test values of corresponding permanent deformation parameters [24–27]. However, parameters of different projects varied entirely; as a result, a set of initial parameters were determined firstly referring to the similar engineering and the inversion range of parameter was preliminarily determined as Table 2.
4.4.1. Parameter Sensitivity Analysis
Disturb a parameter with a fixed step size 10% and keep the other parameters the same. And then sensitivity analysis of the dam on maximum subsidence was performed, and the result was shown in Figures 9 and 10.
According to the classification results in [17], parameter sensitivity was graded into four levels: highly sensitive, sensitive, moderately sensitive, and insensitive. Figure 10 shows that is a moderately sensitive parameter, , are sensitive ones, and is highly sensitive. Besides, Figure 9 shows that the subsidence under earthquake increases with the increase of and decreases with the increase of .
4.4.2. Comparison of RBF and BP Networks in Inversion
In order to improve the generalization ability of neural network, a training method [28] was taken that samples were partitioned into several subsets, and then the network was trained and appraisal was done at the same time with the change of spread constant. Samples were partitioned into three subsets here for training, validation, and testing. The training and validation datasets were used to determine synoptic weights of the network model whereas the testing dataset is used to evaluate the prediction results. If the performance index (generally mean square error, (MSE)) of errors of training dataset was satisfied, then determine the optimal network according to that of verification dataset.
There were 90 training samples in all generated by the methods of central composite design and orthogonal design, and 9 samples generated at random among which 6 ones were used for verification, and the others were for testing. With observation displacement as input vector and permanent deformation parameter as an output vector, the neural network could be trained. One of the common ways to evaluate performance of the network was by error of mean square root (RMSE). However, backanalysis of permanent deformation parameters by neural network was a problem of multi-input-output, if the network performance was still evaluated like that: Then it could not reflect the error influence of different parameters; that is to say, the results calculated by inversion parameters might deviate considerably from the actual value since the forecasting error of highly sensitive parameter was relatively big, though the total error was small. So in this paper, an objective function for evaluating network performance was put forward to improve fitting accuracy of high sensitive parameter, where sensitivity factors were taken as weights of error, as follows: where is the sample number of validation dataset, is the dimension of the output vector, and and are final parameter output and the actual parameter value, respectively.
The key of RBF network training is to determine the neural number of hidden layer. It has been a common method at present to gradually increase the neural number automatically by checking the output error, starting from zero. So the process of establishing a network is also the one of training. Spread constant was adjusted until the prediction accuracy of testing samples met the engineering requirements, and then the trained network can be practically used for parameter inversion. It is beneficial for the generalization of neural network to adjust spread constant, the higher the , the more smooth the fitting, but not the higher the better, if is too high, all outputs are likely to become common and the base functions tend to overlap completely. Its value generally depends on the distance between the input vectors. To analyze the effectiveness of the training, the transition of the RMSE was traced with the change of as shown in Figure 11. It can be seen that the RMSE decreases rapidly in the initial stages and almost remains the same after .
The process of BP network training was similar to the RBF network; after MSE of training dataset was satisfied, then determine the optimal network according to RMSE of verification dataset. The maximum number of training epochs was set to 3000. The MSE goal was set to 0.001. The initial value of adaptive parameter was set to 0.001. increased by 10 and decreased by 0.1 until performance value reduced, and its maximum value is set to 1e10. During the training phase, the data were processed several times to see whether any changes occurred due to the assignment of random initial weights. Figure 12 describes the results that the value of RMSE decreases from the largest value of 0.98 with one hidden neuron to the value of 0.18 with 11 hidden neurons. RMSE was then stable at this value with an increasing number of hidden neurons.
To test the prediction accuracy of the networks, the outputs of the testing samples were listed in Table 3. For the three samples, the prediction accuracy of all parameters by RBFNN is around 5% whereas not all the results by BPNN are satisfactory. However, for both networks, the forecasting error of parameter is generally smaller, which indicates that it has a pronounced effect on improvement of highly sensitive parameters with the function of error evaluation, in the problem of multiple parameter inversion.
To further check the error precision of influence on subsidence, the predictive parameters by RBFNN were used for finite element calculations. Owing to space reasons, for the second sample only, the results were compared with the actual value as in Figure 13, and the observation points were on central axis of the typical section. It can be seen that the computation values are very close, which indicates that the prediction effect by the RBF network can basically meet engineering accuracy.
With actual observation displacement of Zipingpu dam as input vector, permanent deformation parameters of rockfill were backanalyzed by the trained RBF network, and the results are shown in Table 4.
The permanent displacement calculated by the inversion parameters is shown in Figure 14, and the marked values are actual displacements of the dam. It can be seen that the vertical displacements increase with the increment of dam height, which accords with the law of the measured settlement, and both values are also very close. Figure 15 shows the deformed mesh for finite element calculations, where deformation has already been enlarged to see clearly. Besides, the dam section becomes smaller and the slopes of both upstream and downstream shrink inward, which embodies the macroscopic character of rockfill under earthquake action. Due to the upstream water load, the maximum of the horizontal permanent deformation points to the downstream of the dam but the earthquake-induced permanent deformation is predominantly vertical settlement. In summary, the results are qualitatively rational, which indicates that it is feasible to calculate the earthquake-induced permanent deformation by the method of equivalent nodal force on the basis of improved Shen Zhujiang model.
5. Conclusions
In the process of conventional back/inverse analysis, it is necessary to perform FEM analysis frequently. For a large-scale multiple parameter and nonlinear problem such as backanalysis of displacement of an earth-rockfill dam, the workload is so daunting that sometimes the backanalysis cannot be carried out. In this paper, ANNs were introduced in the field of dynamic parameter inversion of earth-rockfill dam, BP and RBF networks were compared, then on the basis of RBFNN, a backanalysis model for earthquake-induced permanent deformation parameters was proposed, and it was used for backanalysis of parameters for Zipingpu CFRD. The results indicate the following:(1)The RBFNN model appears more robust and efficient than BPNN model for backanalysis of earthquake-induced permanent deformation parameters of the earth-rockfill dam. Due to the assignment of random initial weights, the structure of BPNN is difficult to determine, network training results are unstable, and it can easily be trap in a local optimum; all of the problems are still to be resolved, and so RBFNN is a better choice. (2)It is an easy and effective way to backanalysis of earthquake-induced permanent deformation parameters of the earth-rockfill dam with RBF network model. In this way, the inversion range of parameters is determined according to the similar engineering, and several samples are generated through the experimental design methods; then after the neural network is trained, the residual deformation parameters can be acquired quickly by putting in actual permanent deformations of the dam.(3)The existing theory and method of dynamic calculation can basically reflect the earthquake-resistant behavior of earth-rockfill dam, but it is still an extremely complex subject to research; many factors such as seismic input and calculation model may lead to great gap between calculated value and actual value. Thereby, further study on dynamic analysis method will make the inversion of permanent deformation parameters more meaningful.
Acknowledgments
This work was supported by National Natural Science Foundation of China (51109028 and 90815024) and the Fundamental Research Funds for the Central Universities (DUT11RC(3)38).