# Inverse Parametric Analysis of Seismic Permanent Deformation for Earth-Rockfill Dams Using Artificial Neural Networks

Authors: Xu Wang; Fei Kang; Junjie Li; Xin Wang
Journal: Mathematical Problems in Engineering
 (2012)

---

## 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.


---

## Body

## 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. Section2 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 c, friction φ (or φ0, Δφ), failure ration Rf, modulus number K, modulus exponent n, bulk modulus number Kb, and bulk modulus exponent m. The nonlinear stress-strain relation of rockfill is represented by a hyperbola, whose instantaneous slope is the tangent modulus Et. According to the conventional triaxial tests, Et can be expressed as follows:
(2.1)Et=[1-Rf(1-sinφ)2c×cosφ+2σ3sinφ(σ1-σ3)]2×K×Pa(σ3Pa)n.
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 σ3. 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:
(2.2)φ=φ0-Δφlg(σ3Pa),
where φ0 is the value of φ for σ3=Pa and Δφ is the reduction of φ for a 10-fold increase in σ3. And the bulk modulus can be expressed as
(2.3)B=KbPa(σ3Pa)m.


### 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 G and equivalent damping ratio λ. The key of the model is to determine the relation between maximum dynamic shear modulus Gmax and mean effective principal stress σ0, as well as the variation of G 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 Gmax can be expressed as
(2.4)Gmax=K×Pa×(σ0Pa)n,
where K,n are modulus and exponent usually determined by experiments. Pa is the atmospheric pressure. In the experiments, the three-axis instrument and wave velocity test device were used. In regard to G and λ, a relation curve is achieved firstly through experiments, which describes the variation of the dynamic shear modulus ratio G/Gmax and damping ratio λ with dynamic shear strain γ. Then γ is normalized by reference shear strain γr to make instantaneous G 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 nodei is calculated by simplified Westergaard formula [10]:
(2.5)mwi=ψ9078ρH0yiAi,
where ψ is the angle between the upstream slop and the horizontal plane H0 is depth of water from the surface to the bottom of reservoirs. yi is depth of water from the surface to node i, and Ai is the corresponding control area of node i.


### 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 εvr and residual shear strain γr can be written as
(2.6)εvr=cvrlg(1+N),γr=cdrlg(1+N),cvr=c1(γd)c2exp(-c3Sl2),cdr=c4(γd)c5Sl2,
where γd is dynamic shear strain amplitude, Sl is stress level, N is the number of vibrations, and c1, c2, c3, c4, and c5 are experimental parameters. The parameter c3 is assumed to have nothing to do with cνγ, that is, c3=0.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:
(2.7)Δενr=c1(γd)c2exp(-c3Sl2)ΔN1+N,Δγr=c4(γd)c5SlΔN1+N,
where Δεvr, Δγr 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σ0 and stress level Sl.(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εvp, Δγrp according to (2.7). And strain potential of finite elements is obtained according to the following formula:
(2.8){Δξp}={ΔξxΔξyΔξzΔγxyΔγyzΔγzx}=13Δεvp{111000}+Δγrmaxpq{σx-pσy-pσz-p2τxy2τyz2τzx},whereq is generalized shear stress,  p is average principal stress, and {Δξp} is increment of residual strain in Cartesian coordinates.(3)
Calculate the equivalent nodal force with the converted strain potential according to the formula that(2.9){ΔF}=∭v[B]T[D]{Δξp}dV,where[B] is the conversion matrix of strain [D] is the elastic matrix. Then permanent deformations are calculated with the equivalent nodal force applied on finite elements.


## 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 c, friction φ (or φ0, Δφ), failure ration Rf, modulus number K, modulus exponent n, bulk modulus number Kb, and bulk modulus exponent m. The nonlinear stress-strain relation of rockfill is represented by a hyperbola, whose instantaneous slope is the tangent modulus Et. According to the conventional triaxial tests, Et can be expressed as follows:
(2.1)Et=[1-Rf(1-sinφ)2c×cosφ+2σ3sinφ(σ1-σ3)]2×K×Pa(σ3Pa)n.
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 σ3. 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:
(2.2)φ=φ0-Δφlg(σ3Pa),
where φ0 is the value of φ for σ3=Pa and Δφ is the reduction of φ for a 10-fold increase in σ3. And the bulk modulus can be expressed as
(2.3)B=KbPa(σ3Pa)m.


## 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 G and equivalent damping ratio λ. The key of the model is to determine the relation between maximum dynamic shear modulus Gmax and mean effective principal stress σ0, as well as the variation of G 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 Gmax can be expressed as
(2.4)Gmax=K×Pa×(σ0Pa)n,
where K,n are modulus and exponent usually determined by experiments. Pa is the atmospheric pressure. In the experiments, the three-axis instrument and wave velocity test device were used. In regard to G and λ, a relation curve is achieved firstly through experiments, which describes the variation of the dynamic shear modulus ratio G/Gmax and damping ratio λ with dynamic shear strain γ. Then γ is normalized by reference shear strain γr to make instantaneous G 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 nodei is calculated by simplified Westergaard formula [10]:
(2.5)mwi=ψ9078ρH0yiAi,
where ψ is the angle between the upstream slop and the horizontal plane H0 is depth of water from the surface to the bottom of reservoirs. yi is depth of water from the surface to node i, and Ai is the corresponding control area of node i.


## 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 εvr and residual shear strain γr can be written as
(2.6)εvr=cvrlg(1+N),γr=cdrlg(1+N),cvr=c1(γd)c2exp(-c3Sl2),cdr=c4(γd)c5Sl2,
where γd is dynamic shear strain amplitude, Sl is stress level, N is the number of vibrations, and c1, c2, c3, c4, and c5 are experimental parameters. The parameter c3 is assumed to have nothing to do with cνγ, that is, c3=0.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:
(2.7)Δενr=c1(γd)c2exp(-c3Sl2)ΔN1+N,Δγr=c4(γd)c5SlΔN1+N,
where Δεvr, Δγr 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σ0 and stress level Sl.(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εvp, Δγrp according to (2.7). And strain potential of finite elements is obtained according to the following formula:
(2.8){Δξp}={ΔξxΔξyΔξzΔγxyΔγyzΔγzx}=13Δεvp{111000}+Δγrmaxpq{σx-pσy-pσz-p2τxy2τyz2τzx},whereq is generalized shear stress,  p is average principal stress, and {Δξp} is increment of residual strain in Cartesian coordinates.(3)
Calculate the equivalent nodal force with the converted strain potential according to the formula that(2.9){ΔF}=∭v[B]T[D]{Δξp}dV,where[B] is the conversion matrix of strain [D] is the elastic matrix. Then permanent deformations are calculated with the equivalent nodal force applied on finite elements.


## 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 εvr and residual shear strain γr can be written as
(2.6)εvr=cvrlg(1+N),γr=cdrlg(1+N),cvr=c1(γd)c2exp(-c3Sl2),cdr=c4(γd)c5Sl2,
where γd is dynamic shear strain amplitude, Sl is stress level, N is the number of vibrations, and c1, c2, c3, c4, and c5 are experimental parameters. The parameter c3 is assumed to have nothing to do with cνγ, that is, c3=0.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:
(2.7)Δενr=c1(γd)c2exp(-c3Sl2)ΔN1+N,Δγr=c4(γd)c5SlΔN1+N,
where Δεvr, Δγr 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σ0 and stress level Sl.(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εvp, Δγrp according to (2.7). And strain potential of finite elements is obtained according to the following formula:
(2.8){Δξp}={ΔξxΔξyΔξzΔγxyΔγyzΔγzx}=13Δεvp{111000}+Δγrmaxpq{σx-pσy-pσz-p2τxy2τyz2τzx},whereq is generalized shear stress,  p is average principal stress, and {Δξp} is increment of residual strain in Cartesian coordinates.(3)
Calculate the equivalent nodal force with the converted strain potential according to the formula that(2.9){ΔF}=∭v[B]T[D]{Δξp}dV,where[B] is the conversion matrix of strain [D] 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 structurem-k-p, that is, a vector input x(x1,x2,…,xm), k hidden units, and an output vector y(y1,y2,…,yp)  as in Figure 1, the equation that expresses the relationship between the input and output can be written as
(3.1)yj=foutput[∑h=1kwjhfhidden(∑i=1mwhixi+bias1h)]+bias2j,
where m is the number of input units, k is the number of neurons in the hidden layer, xi is the ith input unit, ωhi is the weight parameter between input i and hidden neuron h, ωjh is the weight parameter between hidden neuron h and output neuron j, fhidden is the activation function of the hidden layer, and foutput is the transfer function of output layer.Figure 1
Structure of BP neutral network.The weights were estimated and adjusted in the learning process with an aim of minimizing an error functionEd as follows:
(3.2)ED=∑i=1n(yi-ti)2=∑i=1nei2,
where n is the number of input and output examples of the training dataset and t is the target value. The errors were fed backward through the network to adjust the weights until the error ED 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 ED, 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:
(3.3)wi+1=wi-(JiTJi+μiI)-1JiTei,
where J is the Jacobian matrix of output errors, I is the identity matrix, and μ is an adaptive parameter. When μ=0, 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. Figure2 displays an RBF network with a structure m-k-p that there are m inputs, k hidden units and p outputs. x=(x1,x2,…,xm)T∈Rm is the input vector. W is the weight matrix of inputs and W∈Rk×p,b1,…,bp are offset of output units. y=(y1,…,yp)T∈Rp is the output vector of the network, and fi(∥x-ci∥) is the activation function of the hidden unit i; one common function is the Gaussian function:
(3.4)φ(x)=exp(-x22σ2),
where σ is spread constant, the role of which is to adjust sensitivity of the Gaussian function.Figure 2
Structure of RBF neutral network.The final output of uniti can be expressed as
(3.5)yi=∑j=1kwijexp(∥x-cj∥22σj2)+bi,
where wij is the weight of output layer and σj 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 toN=k2+2k+cp, where k is the factor number and cp is the replicate number of the central point.(2)
α-values depend on the number of variables and can be calculated by α=2(k-p)/4. For two, three, and four variables, it equals 1.41, 1.68, and 2.00, respectively.(3)
All factors are studied in five levels(-α,-1,0,+1,+α).Figure 3
Central composite design schemes for three factors.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:
(3.6)S=∑i=1q-1((Yi+1-Yi)/Y0)/((Pi+1-Pi)/100)(q-1),
where S is the sensitivity factor, Yi is the output of condition i, Y0 is the output derived from the initial parameters, pi is the percentage of parameter variations in condition i compared to the initial parameters, and q is the number of conditions.


## 3.1. BP Neural Networks

For a BPNN with a structurem-k-p, that is, a vector input x(x1,x2,…,xm), k hidden units, and an output vector y(y1,y2,…,yp)  as in Figure 1, the equation that expresses the relationship between the input and output can be written as
(3.1)yj=foutput[∑h=1kwjhfhidden(∑i=1mwhixi+bias1h)]+bias2j,
where m is the number of input units, k is the number of neurons in the hidden layer, xi is the ith input unit, ωhi is the weight parameter between input i and hidden neuron h, ωjh is the weight parameter between hidden neuron h and output neuron j, fhidden is the activation function of the hidden layer, and foutput is the transfer function of output layer.Figure 1
Structure of BP neutral network.The weights were estimated and adjusted in the learning process with an aim of minimizing an error functionEd as follows:
(3.2)ED=∑i=1n(yi-ti)2=∑i=1nei2,
where n is the number of input and output examples of the training dataset and t is the target value. The errors were fed backward through the network to adjust the weights until the error ED 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 ED, 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:
(3.3)wi+1=wi-(JiTJi+μiI)-1JiTei,
where J is the Jacobian matrix of output errors, I is the identity matrix, and μ is an adaptive parameter. When μ=0, 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. Figure2 displays an RBF network with a structure m-k-p that there are m inputs, k hidden units and p outputs. x=(x1,x2,…,xm)T∈Rm is the input vector. W is the weight matrix of inputs and W∈Rk×p,b1,…,bp are offset of output units. y=(y1,…,yp)T∈Rp is the output vector of the network, and fi(∥x-ci∥) is the activation function of the hidden unit i; one common function is the Gaussian function:
(3.4)φ(x)=exp(-x22σ2),
where σ is spread constant, the role of which is to adjust sensitivity of the Gaussian function.Figure 2
Structure of RBF neutral network.The final output of uniti can be expressed as
(3.5)yi=∑j=1kwijexp(∥x-cj∥22σj2)+bi,
where wij is the weight of output layer and σj 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 toN=k2+2k+cp, where k is the factor number and cp is the replicate number of the central point.(2)
α-values depend on the number of variables and can be calculated by α=2(k-p)/4. For two, three, and four variables, it equals 1.41, 1.68, and 2.00, respectively.(3)
All factors are studied in five levels(-α,-1,0,+1,+α).Figure 3
Central composite design schemes for three factors.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:
(3.6)S=∑i=1q-1((Yi+1-Yi)/Y0)/((Pi+1-Pi)/100)(q-1),
where S is the sensitivity factor, Yi is the output of condition i, Y0 is the output derived from the initial parameters, pi is the percentage of parameter variations in condition i compared to the initial parameters, and q 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 Table1. Besides, the linear elastic model was adopted for the calculation of concrete panels, concrete strength grade was C25, and the corresponding material parameters were density ρ=2400kg/m3, E=2.8×104MPa, Poisson’s ratio μ=0.167.Table 1
The principal parameters for static calculation.
Material  zoning
ρ (kg/m3)
φ (°)
Δ
φ (°)
R
f
K
n
K
b
m
Cushion layer
2300
57.51
10.65
0.84
1274
0.44
1260
−0.026
Transition layer
2250
57.63
11.44
0.75
1153
0.38
1085
−0.089
Main rockfill
2160
55.39
10.60
0.75
1120
0.32
490
0.12
Secondary rockfill
2150
55.39
10.60
0.71
1033
0.38
338
0.03
Covering layer
2150
49.00
10.00
0.80
820
0.40
430
0.25To 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 Figure4. 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.Figure 4
3D mesh of the dam for calculation.Figure 5
Continuous contours of vertical displacement of the dam before earthquake (m).Figure 6
Continuous contours of major stress of the dam before earthquake (MPa).


### 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.Figure 7
Epicenter of “5.12” Wenchuan Earthquake and observation stations.Figure 8
Acceleration time histories of dam site (a/g).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 werec1,c2,c4,c5 but c3 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.Table 2
Ranges of parameters to be inversed and initial values.
Parameter variation
c
1
c
2
c
4
c
5
Inversion range of parameters
0.005
0.5
0.1
0.75
~0.01
~1.0
~0.2
~1.50
Initial parameters
0.0064
0.75
0.18
1.33


#### 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 Figures9 and 10.Figure 9
Trend of vertical subsidence change.Figure 10
Result of parameter sensitivity analysis.According to the classification results in [17], parameter sensitivity was graded into four levels: |Si|≥1 highly sensitive, 0.2≤|Si|<1 sensitive, 0.05≤|Si|<0.2 moderately sensitive, and 0≤|Si|<0.05 insensitive.  Figure 10 shows that c1 is a moderately sensitive parameter, c2, c4 are sensitive ones, and c5 is highly sensitive. Besides, Figure 9 shows that the subsidence under earthquake increases with the increase of c1,c4 and decreases with the increase of c2,c5.


#### 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:(4.1)RMSE=∑i=1n∑j=1p(yij-tij)2n.
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:
(4.2)RMSE=∑i=1n∑j=1p[(yij-tij)*Sj]2n,
where n is the sample number of validation dataset, p is the dimension of the output vector, and yij and tij 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 σ=1.6.Figure 11
History of RMSE along with expand constant.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.Figure 12
History of RMSE along with the number of hidden neurons.To test the prediction accuracy of the networks, the outputs of the testing samples were listed in Table3. 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 c5 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.Table 3
Comparison of predicted results and actual values of the test samples’ parameters.
Testing samples
RBF
BP
c
1
c
2
c
4
c
5
c
1
c
2
c
4
c
5
1
Actual value
0.0064
0.7500
0.1850
1.3300
0.0064
0.7500
0.1850
1.3300
Predictive value
0.0064
0.7658
0.1894
1.3263
0.0073
0.8075
0.1851
1.3250
Error/%
0.0000
2.1067
2.3784
−0.2782
14.5467
7.6702
0.0876
0.3763
2
Actual value
0.0058
0.6750
0.1665
1.1970
0.0058
0.6750
0.1665
1.1970
Predictive value
0.0062
0.7226
0.1753
1.1950
0.0070
0.7525
0.1675
1.1895
Error/%
6.8966
7.05185
5.2853
−0.1671
20.7000
11.4833
0.6283
0.6293
3
Actual value
0.0070
0.8250
0.2035
1.4630
0.0070
0.8250
0.2035
1.4630
Predictive value
0.0068
0.8127
0.1954
1.4354
0.0077
0.8645
0.1982
1.4496
Error/%
−2.8571
−1.4909
−3.9803
−1.8865
9.7816
4.7920
2.5910
0.9153To 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 Figure13, 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.Figure 13
Comparison of dam subsidence.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 Table4.Table 4
Backanalysis results of rockfill deformation parameters.
Material
c
1
c
2
c
4
c
5
Rockfill
0.0072
0.8000
0.1342
1.4477The permanent displacement calculated by the inversion parameters is shown in Figure14, 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.Figure 14
Comparison of dam subsidence.Figure 15
The deformation of the dam grid figure.


## 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 Table1. Besides, the linear elastic model was adopted for the calculation of concrete panels, concrete strength grade was C25, and the corresponding material parameters were density ρ=2400kg/m3, E=2.8×104MPa, Poisson’s ratio μ=0.167.Table 1
The principal parameters for static calculation.
Material  zoning
ρ (kg/m3)
φ (°)
Δ
φ (°)
R
f
K
n
K
b
m
Cushion layer
2300
57.51
10.65
0.84
1274
0.44
1260
−0.026
Transition layer
2250
57.63
11.44
0.75
1153
0.38
1085
−0.089
Main rockfill
2160
55.39
10.60
0.75
1120
0.32
490
0.12
Secondary rockfill
2150
55.39
10.60
0.71
1033
0.38
338
0.03
Covering layer
2150
49.00
10.00
0.80
820
0.40
430
0.25To 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 Figure4. 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.Figure 4
3D mesh of the dam for calculation.Figure 5
Continuous contours of vertical displacement of the dam before earthquake (m).Figure 6
Continuous contours of major stress of the dam before earthquake (MPa).


## 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.Figure 7
Epicenter of “5.12” Wenchuan Earthquake and observation stations.Figure 8
Acceleration time histories of dam site (a/g).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 werec1,c2,c4,c5 but c3 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.Table 2
Ranges of parameters to be inversed and initial values.
Parameter variation
c
1
c
2
c
4
c
5
Inversion range of parameters
0.005
0.5
0.1
0.75
~0.01
~1.0
~0.2
~1.50
Initial parameters
0.0064
0.75
0.18
1.33


### 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 Figures9 and 10.Figure 9
Trend of vertical subsidence change.Figure 10
Result of parameter sensitivity analysis.According to the classification results in [17], parameter sensitivity was graded into four levels: |Si|≥1 highly sensitive, 0.2≤|Si|<1 sensitive, 0.05≤|Si|<0.2 moderately sensitive, and 0≤|Si|<0.05 insensitive.  Figure 10 shows that c1 is a moderately sensitive parameter, c2, c4 are sensitive ones, and c5 is highly sensitive. Besides, Figure 9 shows that the subsidence under earthquake increases with the increase of c1,c4 and decreases with the increase of c2,c5.


### 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:(4.1)RMSE=∑i=1n∑j=1p(yij-tij)2n.
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:
(4.2)RMSE=∑i=1n∑j=1p[(yij-tij)*Sj]2n,
where n is the sample number of validation dataset, p is the dimension of the output vector, and yij and tij 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 σ=1.6.Figure 11
History of RMSE along with expand constant.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.Figure 12
History of RMSE along with the number of hidden neurons.To test the prediction accuracy of the networks, the outputs of the testing samples were listed in Table3. 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 c5 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.Table 3
Comparison of predicted results and actual values of the test samples’ parameters.
Testing samples
RBF
BP
c
1
c
2
c
4
c
5
c
1
c
2
c
4
c
5
1
Actual value
0.0064
0.7500
0.1850
1.3300
0.0064
0.7500
0.1850
1.3300
Predictive value
0.0064
0.7658
0.1894
1.3263
0.0073
0.8075
0.1851
1.3250
Error/%
0.0000
2.1067
2.3784
−0.2782
14.5467
7.6702
0.0876
0.3763
2
Actual value
0.0058
0.6750
0.1665
1.1970
0.0058
0.6750
0.1665
1.1970
Predictive value
0.0062
0.7226
0.1753
1.1950
0.0070
0.7525
0.1675
1.1895
Error/%
6.8966
7.05185
5.2853
−0.1671
20.7000
11.4833
0.6283
0.6293
3
Actual value
0.0070
0.8250
0.2035
1.4630
0.0070
0.8250
0.2035
1.4630
Predictive value
0.0068
0.8127
0.1954
1.4354
0.0077
0.8645
0.1982
1.4496
Error/%
−2.8571
−1.4909
−3.9803
−1.8865
9.7816
4.7920
2.5910
0.9153To 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 Figure13, 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.Figure 13
Comparison of dam subsidence.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 Table4.Table 4
Backanalysis results of rockfill deformation parameters.
Material
c
1
c
2
c
4
c
5
Rockfill
0.0072
0.8000
0.1342
1.4477The permanent displacement calculated by the inversion parameters is shown in Figure14, 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.Figure 14
Comparison of dam subsidence.Figure 15
The deformation of the dam grid figure.


## 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 Figures9 and 10.Figure 9
Trend of vertical subsidence change.Figure 10
Result of parameter sensitivity analysis.According to the classification results in [17], parameter sensitivity was graded into four levels: |Si|≥1 highly sensitive, 0.2≤|Si|<1 sensitive, 0.05≤|Si|<0.2 moderately sensitive, and 0≤|Si|<0.05 insensitive.  Figure 10 shows that c1 is a moderately sensitive parameter, c2, c4 are sensitive ones, and c5 is highly sensitive. Besides, Figure 9 shows that the subsidence under earthquake increases with the increase of c1,c4 and decreases with the increase of c2,c5.


## 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:(4.1)RMSE=∑i=1n∑j=1p(yij-tij)2n.
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:
(4.2)RMSE=∑i=1n∑j=1p[(yij-tij)*Sj]2n,
where n is the sample number of validation dataset, p is the dimension of the output vector, and yij and tij 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 σ=1.6.Figure 11
History of RMSE along with expand constant.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.Figure 12
History of RMSE along with the number of hidden neurons.To test the prediction accuracy of the networks, the outputs of the testing samples were listed in Table3. 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 c5 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.Table 3
Comparison of predicted results and actual values of the test samples’ parameters.
Testing samples
RBF
BP
c
1
c
2
c
4
c
5
c
1
c
2
c
4
c
5
1
Actual value
0.0064
0.7500
0.1850
1.3300
0.0064
0.7500
0.1850
1.3300
Predictive value
0.0064
0.7658
0.1894
1.3263
0.0073
0.8075
0.1851
1.3250
Error/%
0.0000
2.1067
2.3784
−0.2782
14.5467
7.6702
0.0876
0.3763
2
Actual value
0.0058
0.6750
0.1665
1.1970
0.0058
0.6750
0.1665
1.1970
Predictive value
0.0062
0.7226
0.1753
1.1950
0.0070
0.7525
0.1675
1.1895
Error/%
6.8966
7.05185
5.2853
−0.1671
20.7000
11.4833
0.6283
0.6293
3
Actual value
0.0070
0.8250
0.2035
1.4630
0.0070
0.8250
0.2035
1.4630
Predictive value
0.0068
0.8127
0.1954
1.4354
0.0077
0.8645
0.1982
1.4496
Error/%
−2.8571
−1.4909
−3.9803
−1.8865
9.7816
4.7920
2.5910
0.9153To 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 Figure13, 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.Figure 13
Comparison of dam subsidence.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 Table4.Table 4
Backanalysis results of rockfill deformation parameters.
Material
c
1
c
2
c
4
c
5
Rockfill
0.0072
0.8000
0.1342
1.4477The permanent displacement calculated by the inversion parameters is shown in Figure14, 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.Figure 14
Comparison of dam subsidence.Figure 15
The deformation of the dam grid figure.


## 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.



---

*Source: 383749-2012-12-27.xml*