# Back Analysis of the Permeability Coefficient of a High Core Rockfill Dam Based on a RBF Neural Network Optimized Using the PSO Algorithm

Authors: Shichun Chi; Shasha Ni; Zhenping Liu
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

It is important to determine the seepage field parameters of a high core rockfill dam using the seepage data obtained during operation. For the Nuozhadu high core rockfill dam, a back analysis model is proposed using the radial basis function neural network optimized using a particle swarm optimization algorithm (PSO-RBFNN) and the technology of finite element analysis for solving the saturated-unsaturated seepage field. The recorded osmotic pressure curves of osmometers, which are distributed in the maximum cross section, are applied to this back analysis. The permeability coefficients of the dam materials are retrieved using the measured seepage pressure values while the steady state seepage condition exists; that is, the water lever remains unchanged. Meanwhile, the parameters are tested using the unstable saturated-unsaturated seepage field while the water level rises. The results show that the permeability coefficients are reasonable and can be used to study the real behavior of a seepage field of a high core rockfill dam during its operation period.


---

## Body

## 1. Introduction

The earth-rock dam design is widely adopted locally and abroad because of low investment, locally produced raw materials, and simple construction. As a result of complex situations, many influencing factors, and poor management, a number of accidents have occurred over the past several years. Many accidents have made it clear that the disaster caused by a flood once a dam breaks is unimaginable. According to the International Commission on Large Dams [1, 2], the number of earth-rock dam breaks accounted for 70 percent of the total number of dam breaks; 25 percent of the earth-rock dam breaks were caused by infiltration or deformation and the percentage was up to 40 percent in China [3]. Seepage failure is one of the main causes of collapse of an earth-rock dam; therefore, it is very important to use mature theories and methods to analyze and evaluate leak and deformation conditions.Generally speaking, there are three types of methods used in permeability coefficient identification: analysis, testing, and back analysis methods. An analysis method, which is based on numerous hypotheses, is difficult to adapt to a dam with a complex structure. Some permeability coefficients can be measured in the field; however, such measurements are usually scarce and are expensive to obtain. Therefore, a back analysis method to fill in the missing parameters is typically conducted before an analysis of the seepage flow. During the past three decades, extensive research has been conducted to study the problem of permeability coefficient estimation for seepage flow models. This research resulted in the development of numerous techniques that employ search methods to find parameter values that minimize the difference between the observed and calculated hydraulic head values. A brief review of the approaches employed in these techniques is provided herein.In 1986, research about numerical methods for permeability coefficient identification was conducted [4]. Tsurumi et al. [5] used a finite element method to identify parameters in groundwater hydrology. Zijlstra and Dane [6] used a quasi-Newton method to identify hydraulic parameters in layered soils. El Harrouni et al. [7] presented a method of estimating the groundwater parameter using optimization. Prasad and Rastogi [8] studied an inverse method of estimating the hydraulic conductivity values based on a genetic algorithm. Zio [9] studied the inverse problem of parameter estimation in groundwater models using artificial neural networks. At home, He and Zheng [10] deduced a series of formulas about numerical inversion of Laplace transform solutions for the dynamics of porous flow. Liu and Wang [11] studied an improved genetic algorithm for back analysis of seepage parameters. Liu et al. [12, 13] studied the back analysis of seepage with an ANN based on alternative and iterative algorithm ANN based on the simulated annealing Gauss-Newton algorithm. Li et al. [14] proposed ant colony optimization to estimate aquifer parameters.In recent years, neural networks (NNs) have been successfully used in the field of civil engineering to directly map nonlinear complex relations. As is stated in [15], the mapping capabilities of a NN are strictly related to the nonlinear component found in the activation function of the neurons. A feed forward neural network (FFNN) is a NN where the inner architecture is organized in subsequent layers of neurons, and the connections are made according to the following rules: every neuron of a layer is connected to all (and only) the neurons of the subsequent layer. The commonly used backpropagation algorithm for FFNN training suffers from slow learning speed and being liable to be trapped in a local minimum. The radial basis function (RBF) neural networks are normally considered as universal approximators [16–18], which can approximate any function to an arbitrary accuracy, provided that the size of network is not constrained. Applications of RBF neural network are widespread and can be found in prediction of time series [19], function approximation problems [20], hydrological modeling [21, 22], and system modeling [23]. General description of the RBF neural network can be found in many of the standard ANN text books [24]. The present work proposes the use of RBF neural network (choosing Gaussian function as its activation function) for mapping the complex nonlinear relations between water heads and permeability coefficients. The RBF neural network has a simple structure, succinct training, fast convergence speed, and the ability to represent any complicated nonlinear function relations. It creates a radial basis network one neuron at a time. Neurons are added to the network until the mean square error falls beneath an error goal or a maximum number of neurons are reached. As you know, the more the input vectors are, the more the neurons are needed to achieve a better fitting effect and forecast precision. It easily leads to dimension disaster [25]. To avoid this phenomenon, this paper selects an optimizing algorithm to determine the RBF neural network structural parameters, which can achieve a better fitting effect and forecast precision with fewer hidden layer neurons.Optimum parameters of RBF neural network can be determined by many methods of varying complexity, some of which are noniterative clustering [26], k-means clustering [27], and orthogonal least squares (OLS) [28, 29]. Evolutionary algorithms also have been very popular in determining the parameters of RBF networks. A genetic algorithm was used in [30] to evolve the network for function approximation. The genetic algorithm population consisted of individual kernels instead of entire networks. A similar strategy has been used in [31] which makes use of an OLS procedure and singular value decomposition to evaluate the contribution of individual kernels. In [32] the number of kernels is also determined dynamically using a genetic algorithm. Other methods (e.g., [33]) evolve an entire population of RBF neural networks. PSO is similar to the genetic algorithm in the sense that they are both population-based search approaches and that they both depend on information sharing among their population members to enhance their search processes using a combination of deterministic and probabilistic rules. Although PSO and the GA on average yield the same effectiveness (solution quality), PSO is more computationally efficient (uses less number of function evaluations) than the GA [34]. Its small computational requirement makes it a good candidate for solving optimization problems. Solutions computed by PSO are derived from both local and global searches. This allows PSO to consider solutions near the vicinity (local) of the starting point, in addition to a possible solution located within a global region. The search process for new point(s) includes consideration of the previous best solutions. Thus, PSO is used to train the neural network to obtain the optimal values of the structure parameters for the RBF neural network in this paper.In this paper, a back analysis model is proposed with the radial basis function neural network optimized using the particle swarm optimization algorithm (PSO-RBFNN) and the technology of finite elements for solving the saturated-unsaturated seepage field. The recorded osmotic pressure curves of osmometers, which are distributed in the maximum cross section, are applied to this back analysis. The permeability coefficients of the dam materials are retrieved from the prototype data measured while the dam was in the steady state seepage condition; that is, the water lever remained unchanged. Meanwhile, the parameters are tested using the unstable saturated-unsaturated seepage field while the water level rises. The results show that the permeability coefficients are reasonable and can be used to study the real behavior of a seepage field of a high core rockfill dam while operating.


## 2. Mathematical Model

The objective of the inverse model or operator is to obtain optimal estimates of unknown permeability coefficients that minimize the difference between observed and calculated hydraulic head values. In this case, the optimization mathematical model of the back analysis of permeability coefficients is as follows [35]:(1)Min⁡S=1n∑i=1nwiHi-Hi∗2i=1,2,…,n(2)Subject  to:kj_≤kj≤kj¯j=1,2,…,m.In the above model, (1) is the objective function which achieves minimum of the mean squared error of water heads; Hi∗ and Hi are the ith measured and computed water head, respectively; wi is the ith weight (equal weights were used in this study); n is the number of water head measuring points. Equation (2) gives the range of permeability coefficients, kj; kj¯ and kj_ are the maximum and minimum of kj, respectively.


## 3. Seepage Flow Model

### 3.1. Partial Differential Water Flow Equations [36]

The flow of water through both saturated and unsaturated soil follows Darcy’s law which states that(3)q=ki,where q is the specific discharge; k is the hydraulic conductivity; and i is the gradient of total hydraulic head.According to the water continuity and Darcy’s law, a two-dimensional unsteady flow through an earthfill dam can be described by the Richards equation as(4)∂∂xkxθw∂H∂x+∂∂ykyθw∂H∂y+Q=∂θw∂t,where H is the total head; kx and ky are the hydraulic conductivity in the x- and y-directions, respectively; θw is the volumetric water content; Q is the applied boundary flux; and t is the time.Equation (4) can be employed to simulate two-dimensional unsteady state water flow through nonhomogeneous, anisotropic, saturated-unsaturated porous media receiving lateral flow. It states that the difference between the flow (flux) entering and leaving an elemental volume at a point in time is equal to the change in storage of the soil systems. More fundamentally, it states that the sum of the rates of change of flows in the x- and y-directions plus the external applied flux is equal to the rate of change of the volumetric water content with respect to time.Under steady state condition, the flux entering and leaving an elemental volume is the same at all times. The right side of (4) consequently vanishes and the equation reduces to(5)∂∂xkxθw∂H∂x+∂∂ykyθw∂H∂y+Q=0.Changes in volumetric water content are dependent on changes in the stress state and the properties of the soil. The stress state for both saturated and unsaturated conditions can be described by two state variables. These stress state variables areσ-ua and ua-uw where σ is the total stress, ua is the pore-air pressure, and uw is the pore-water pressure. It is usually assumed that the total stress and the pore-air pressure are both constants during transient processes. This means that σ-ua remains constant and has no effect on the change in volumetric water content. Changes in volumetric water content are consequently dependent only on changes in the ua-uw stress state variable, and, with ua remaining constant, the change in volumetric water content is a function only of pore-water pressure changes. As a result, the change in volumetric water content can be related to a change in pore-water pressure by the following equation:(6)∂θw=mw∂uw,where mw is the slope of the storage curve.The total hydraulic head,H, is defined as(7)H=uwγw+y,where γw is the unit weight of water; uw is the pore-water pressure; and y is the elevation.Equation (7) can be rearranged as uw=γw(H-y).Substituting (7) into (6) gives the following equation:(8)∂θ=mwγw∂H-y,which now can be substituted into (4), leading to the following expression:(9)∂∂xkx∂H∂x+∂∂yky∂H∂y+Q=mwγw∂H∂t.


### 3.2. Finite Element Water Flow Equations [36]

Applying the Galerkin method of weighted residual to the governing differential equation, the finite element method for two-dimensional seepage equation can be derived as(10)τ∫ABTCBdAH+τ∫AλNTNdAH,t=q∫LNTRdL,where [B] is the gradient matrix; [C] is the element hydraulic conductivity matrix; {H} is the vector of nodal heads; N is the vector of interpolating function; q is the unit flux across the edge of an element; τ is the thickness of an element; t is the time; λ is the storage term for a transient seepage equal to mwγw; A is a designation for summation over the area of an element; and L=a is a designation for summation over the edge of an element.In an abbreviated form, the finite element seepage equation can be expressed as(11)KH+MH,t=Q,where [K] is the element characteristic matrix; [M] is the element mass matrix; Q is the element applied flux vector.


## 3.1. Partial Differential Water Flow Equations [36]

The flow of water through both saturated and unsaturated soil follows Darcy’s law which states that(3)q=ki,where q is the specific discharge; k is the hydraulic conductivity; and i is the gradient of total hydraulic head.According to the water continuity and Darcy’s law, a two-dimensional unsteady flow through an earthfill dam can be described by the Richards equation as(4)∂∂xkxθw∂H∂x+∂∂ykyθw∂H∂y+Q=∂θw∂t,where H is the total head; kx and ky are the hydraulic conductivity in the x- and y-directions, respectively; θw is the volumetric water content; Q is the applied boundary flux; and t is the time.Equation (4) can be employed to simulate two-dimensional unsteady state water flow through nonhomogeneous, anisotropic, saturated-unsaturated porous media receiving lateral flow. It states that the difference between the flow (flux) entering and leaving an elemental volume at a point in time is equal to the change in storage of the soil systems. More fundamentally, it states that the sum of the rates of change of flows in the x- and y-directions plus the external applied flux is equal to the rate of change of the volumetric water content with respect to time.Under steady state condition, the flux entering and leaving an elemental volume is the same at all times. The right side of (4) consequently vanishes and the equation reduces to(5)∂∂xkxθw∂H∂x+∂∂ykyθw∂H∂y+Q=0.Changes in volumetric water content are dependent on changes in the stress state and the properties of the soil. The stress state for both saturated and unsaturated conditions can be described by two state variables. These stress state variables areσ-ua and ua-uw where σ is the total stress, ua is the pore-air pressure, and uw is the pore-water pressure. It is usually assumed that the total stress and the pore-air pressure are both constants during transient processes. This means that σ-ua remains constant and has no effect on the change in volumetric water content. Changes in volumetric water content are consequently dependent only on changes in the ua-uw stress state variable, and, with ua remaining constant, the change in volumetric water content is a function only of pore-water pressure changes. As a result, the change in volumetric water content can be related to a change in pore-water pressure by the following equation:(6)∂θw=mw∂uw,where mw is the slope of the storage curve.The total hydraulic head,H, is defined as(7)H=uwγw+y,where γw is the unit weight of water; uw is the pore-water pressure; and y is the elevation.Equation (7) can be rearranged as uw=γw(H-y).Substituting (7) into (6) gives the following equation:(8)∂θ=mwγw∂H-y,which now can be substituted into (4), leading to the following expression:(9)∂∂xkx∂H∂x+∂∂yky∂H∂y+Q=mwγw∂H∂t.


## 3.2. Finite Element Water Flow Equations [36]

Applying the Galerkin method of weighted residual to the governing differential equation, the finite element method for two-dimensional seepage equation can be derived as(10)τ∫ABTCBdAH+τ∫AλNTNdAH,t=q∫LNTRdL,where [B] is the gradient matrix; [C] is the element hydraulic conductivity matrix; {H} is the vector of nodal heads; N is the vector of interpolating function; q is the unit flux across the edge of an element; τ is the thickness of an element; t is the time; λ is the storage term for a transient seepage equal to mwγw; A is a designation for summation over the area of an element; and L=a is a designation for summation over the edge of an element.In an abbreviated form, the finite element seepage equation can be expressed as(11)KH+MH,t=Q,where [K] is the element characteristic matrix; [M] is the element mass matrix; Q is the element applied flux vector.


## 4. Prediction Model of RBFNN Optimized by PSO

### 4.1. Radial Basis Function Neural Network (RBFNN)

The radial basis function neural network (RBFNN) used in this study is a feed forward neural network with one radial basis layer. It can uniformly approximate any continuous function with a specified accuracy [16–18]. General description of the RBFNN can be found in many of the standard ANN text books [23]. Figure 1 is a radial basis network with R inputs [37].Figure 1
A radial basis network withR inputs.Here the net input to the radbas transfer function is the vector distance between its weight vectorw and the input vector p, multiplied by the bias b (the dist box in this figure accepts the input vector p and the single row input weight matrix and produces the dot product of the two).The transfer function for a radial basis neuron is(12)radbasdist⁡=e-dist2.The most often used RBF is Gaussian function that is characterized by a center and a width, which performs the nonlinear transformation in the hidden layer by measuring the Euclidean distance between the input vectorpj and the radial basis functional center as follows:(13)radbaspj-ci=exp⁡-12σi2pj-ci2i=1,2,…,m.In this function, pj-ci is Euclidean norm; pj is the jth input vector; ci is the center of ith basis function with the same dimension of pj; σi is the ith neuron radial basis function width parameter; i is the number of hidden layer nodes; j=1,2,…,J, j is the number of input vectors, and J is the total number of input vectors.The entire radial basis function neural networks consist of three layers: an input layer, a hidden radial basis layer, and an output linear layer. The network architecture is shown in Figure2 where R is the number of elements in input vector; S1 is the number of neurons in radial basis layer; S2 is the number of elements in output vector.Figure 2
Radial basis function neural network architecture.This paper adopts Gauss function as the activation function; thus, the linear function of output layer is(14)yl=∑i=1kwilexp⁡-12σi2pj-ci2l=1,2,…,S2.In this function, wil is the connection weight between hidden layer and output; yl is the output value of the lth node. Thus it can be seen that the center ci, width parameter σi, and the connection weight wil are very important parameters that affect the accuracy of prediction.Look in detail at how the radial basis layer operates. Each neuron’s weighted input is the distance between the input vector and its weight vector, calculated with dist. Each neuron’s net input is the element-by-element product of its weighted input with its bias, calculated with netprod. Each neuron’s output is its net input passed through radbas. Figure2 shows the network architecture adopted by this study. The hidden layer in the network applies a nonlinear transformation of the input space to the hidden space using hidden neurons (processing units) with a RBF as a nonlinear transfer function to operate on the input data. All of the details of designing this network are built into design functions newrbe and newrb in neural network toolbox [37], and their outputs can be obtained with sim. This paper chooses newrb as it is more efficient. The call for functions is(15)net=newrbPtrain,Ktrain,Kpredict=sim⁡net,Preal.In the above two functions, Ptrain is normalized training samples of the water heads input vectors and Ktrain is normalized training samples of the seepage parameters output vectors. The function newrb takes matrices of Ptrain and Ktrain and returns the desired network. Thus, using the RBFNN mapping, the nonlinear relation between the water heads and seepage parameters is established to approximate the implicit expressions.The RBFNN configuration is formulated as a minimization problem with respect to the number of hidden layer nodes, the center locationsci, width σi, and the connection weights wij. In this paper, the basic idea of searching optimum parameters ci, σi, and wij is as follows: the number of neurons in hidden layer begins at 1; this paper adopts particle swarm optimization (PSO) algorithm as training method to search parameters ci, σi, and wij that minimize the difference between real and predicted matrices of the seepage parameters; the error of the network is checked, and if the error goal is not met, the next neuron is added (the function newrb creates neurons one at a time, and, at each iteration, the input vector that results in lowering the network error the most is used to create a radbas neuron); this procedure is repeated until the error goal is met, or the maximum number of neurons is reached.


### 4.2. Particle Swarm Optimization Algorithm (PSO)

The particle swarm optimization (PSO) was proposed by Kennedy and Eberhart [38]. PSO was inspired from birds flocking and their social behavior. Like other evolutionary algorithms, the PSO generates a random population of solutions which are called particles in this algorithm. Each particle demonstrates a point in the decision space. The particles direction of movement through the search space is governed by the best-known positions of each particle, pbest, and the global best position found during the movement of the swarm, gbest. The position of each particle i in the search space encodes values of the variables of the optimization problem. The objective function, F, of the optimization problem is evaluated for each particles position in the search space. The position that gives the best value of the objective function is pbesti position. The global best position, gbest, is the particle position that has the best evaluation of the objective function in the population. In essence, PSO employs a swarm of particles or possible solutions that fly through the feasible solution space to explore optimal solutions. It can be easily implemented and it is computationally inexpensive. In this paper, the RBFNN combined with the particle swarm optimization is applied to the back analysis of the permeability coefficients. The flow chart of the radial basis function neural network optimized using the particle swarm optimization algorithm is shown in Figure 3.Figure 3
Flow chart of RBF neural network optimized by particle swarm optimization algorithm.The velocity of a particle, in each dimension of the search space, is based on the positions of the best-known solutions,pbesti and gbest, and the current position. For each iteration, n, the velocity in each dimension is given as [39](16)Vidn+1=ωVidn+c1r1Pidn-Xidn+c2r2Pgdn-Xidn,where d=1,2,…,D; i=1,2,…,N, where N is the size of the swarm; r1 and r2 are a random number uniformly distributed in [0,1]; c1 and c2 are positive constant parameters called acceleration coefficients. It has been found that when c1=2 and c2=2, the velocities provide particle movement that improve overall procedure effectiveness; Xidn is the position of particle “i” in the past iteration; Vidn is the velocity in the past iteration; ω is inertia weight that controls the impact of previous velocities and is usually considered to be equal to 1; and n=1,2,… determines the iteration number.Using the new velocity calculated from (16), Vidn+1, and the current position of the particle Xidn, a new position is calculated as(17)Xidn+1=Xidn+Vidn+1,where Xidn+1 and Vidn+1 are position and velocity in the current iteration, respectively.After the particle has moved to the next positionXidn+1, the process starts over with the evaluation of the objective function. This continues for 1 to nmax⁡, number of iterations, or until some convergence criteria are met. As the particles move through the search space, they begin to “swarm” around a position that optimizes the value of the objective function of the problem. In this paper, the objective function which determines the RBFNN performance is represented by mean squared error (MSE), expressed as(18)MSE=1m∑i=1nki-k^i2i=1,2,…,m,where ki and k^i are actual and predicted normalized permeability coefficient vectors, respectively, and m is the total number of data sets.The PSO-RBFNN optimization framework is implemented in MATLAB. The swarm population size is set to 30, and one particle represents a three-dimensional vector (the elements areci, σi, and wij). The input and output vectors of RBFNN are all normalized. Given that the input data space can only span from 0 to 1, the RBF center ci can logically be constrained by the range [0,1]. According to the practice of the writers, the width σi and the connection weights wij are constrained by the ranges [1,30] and [0,100], respectively. Particle swarm is firstly initialized randomly within the given scope of the parameters. Then the particles move according to the flow chart that is shown in Figure 3. The end conditions are that the number of iterations reaches 100 or until the mean squared error is less than or equal to 0.0001.


## 4.1. Radial Basis Function Neural Network (RBFNN)

The radial basis function neural network (RBFNN) used in this study is a feed forward neural network with one radial basis layer. It can uniformly approximate any continuous function with a specified accuracy [16–18]. General description of the RBFNN can be found in many of the standard ANN text books [23]. Figure 1 is a radial basis network with R inputs [37].Figure 1
A radial basis network withR inputs.Here the net input to the radbas transfer function is the vector distance between its weight vectorw and the input vector p, multiplied by the bias b (the dist box in this figure accepts the input vector p and the single row input weight matrix and produces the dot product of the two).The transfer function for a radial basis neuron is(12)radbasdist⁡=e-dist2.The most often used RBF is Gaussian function that is characterized by a center and a width, which performs the nonlinear transformation in the hidden layer by measuring the Euclidean distance between the input vectorpj and the radial basis functional center as follows:(13)radbaspj-ci=exp⁡-12σi2pj-ci2i=1,2,…,m.In this function, pj-ci is Euclidean norm; pj is the jth input vector; ci is the center of ith basis function with the same dimension of pj; σi is the ith neuron radial basis function width parameter; i is the number of hidden layer nodes; j=1,2,…,J, j is the number of input vectors, and J is the total number of input vectors.The entire radial basis function neural networks consist of three layers: an input layer, a hidden radial basis layer, and an output linear layer. The network architecture is shown in Figure2 where R is the number of elements in input vector; S1 is the number of neurons in radial basis layer; S2 is the number of elements in output vector.Figure 2
Radial basis function neural network architecture.This paper adopts Gauss function as the activation function; thus, the linear function of output layer is(14)yl=∑i=1kwilexp⁡-12σi2pj-ci2l=1,2,…,S2.In this function, wil is the connection weight between hidden layer and output; yl is the output value of the lth node. Thus it can be seen that the center ci, width parameter σi, and the connection weight wil are very important parameters that affect the accuracy of prediction.Look in detail at how the radial basis layer operates. Each neuron’s weighted input is the distance between the input vector and its weight vector, calculated with dist. Each neuron’s net input is the element-by-element product of its weighted input with its bias, calculated with netprod. Each neuron’s output is its net input passed through radbas. Figure2 shows the network architecture adopted by this study. The hidden layer in the network applies a nonlinear transformation of the input space to the hidden space using hidden neurons (processing units) with a RBF as a nonlinear transfer function to operate on the input data. All of the details of designing this network are built into design functions newrbe and newrb in neural network toolbox [37], and their outputs can be obtained with sim. This paper chooses newrb as it is more efficient. The call for functions is(15)net=newrbPtrain,Ktrain,Kpredict=sim⁡net,Preal.In the above two functions, Ptrain is normalized training samples of the water heads input vectors and Ktrain is normalized training samples of the seepage parameters output vectors. The function newrb takes matrices of Ptrain and Ktrain and returns the desired network. Thus, using the RBFNN mapping, the nonlinear relation between the water heads and seepage parameters is established to approximate the implicit expressions.The RBFNN configuration is formulated as a minimization problem with respect to the number of hidden layer nodes, the center locationsci, width σi, and the connection weights wij. In this paper, the basic idea of searching optimum parameters ci, σi, and wij is as follows: the number of neurons in hidden layer begins at 1; this paper adopts particle swarm optimization (PSO) algorithm as training method to search parameters ci, σi, and wij that minimize the difference between real and predicted matrices of the seepage parameters; the error of the network is checked, and if the error goal is not met, the next neuron is added (the function newrb creates neurons one at a time, and, at each iteration, the input vector that results in lowering the network error the most is used to create a radbas neuron); this procedure is repeated until the error goal is met, or the maximum number of neurons is reached.


## 4.2. Particle Swarm Optimization Algorithm (PSO)

The particle swarm optimization (PSO) was proposed by Kennedy and Eberhart [38]. PSO was inspired from birds flocking and their social behavior. Like other evolutionary algorithms, the PSO generates a random population of solutions which are called particles in this algorithm. Each particle demonstrates a point in the decision space. The particles direction of movement through the search space is governed by the best-known positions of each particle, pbest, and the global best position found during the movement of the swarm, gbest. The position of each particle i in the search space encodes values of the variables of the optimization problem. The objective function, F, of the optimization problem is evaluated for each particles position in the search space. The position that gives the best value of the objective function is pbesti position. The global best position, gbest, is the particle position that has the best evaluation of the objective function in the population. In essence, PSO employs a swarm of particles or possible solutions that fly through the feasible solution space to explore optimal solutions. It can be easily implemented and it is computationally inexpensive. In this paper, the RBFNN combined with the particle swarm optimization is applied to the back analysis of the permeability coefficients. The flow chart of the radial basis function neural network optimized using the particle swarm optimization algorithm is shown in Figure 3.Figure 3
Flow chart of RBF neural network optimized by particle swarm optimization algorithm.The velocity of a particle, in each dimension of the search space, is based on the positions of the best-known solutions,pbesti and gbest, and the current position. For each iteration, n, the velocity in each dimension is given as [39](16)Vidn+1=ωVidn+c1r1Pidn-Xidn+c2r2Pgdn-Xidn,where d=1,2,…,D; i=1,2,…,N, where N is the size of the swarm; r1 and r2 are a random number uniformly distributed in [0,1]; c1 and c2 are positive constant parameters called acceleration coefficients. It has been found that when c1=2 and c2=2, the velocities provide particle movement that improve overall procedure effectiveness; Xidn is the position of particle “i” in the past iteration; Vidn is the velocity in the past iteration; ω is inertia weight that controls the impact of previous velocities and is usually considered to be equal to 1; and n=1,2,… determines the iteration number.Using the new velocity calculated from (16), Vidn+1, and the current position of the particle Xidn, a new position is calculated as(17)Xidn+1=Xidn+Vidn+1,where Xidn+1 and Vidn+1 are position and velocity in the current iteration, respectively.After the particle has moved to the next positionXidn+1, the process starts over with the evaluation of the objective function. This continues for 1 to nmax⁡, number of iterations, or until some convergence criteria are met. As the particles move through the search space, they begin to “swarm” around a position that optimizes the value of the objective function of the problem. In this paper, the objective function which determines the RBFNN performance is represented by mean squared error (MSE), expressed as(18)MSE=1m∑i=1nki-k^i2i=1,2,…,m,where ki and k^i are actual and predicted normalized permeability coefficient vectors, respectively, and m is the total number of data sets.The PSO-RBFNN optimization framework is implemented in MATLAB. The swarm population size is set to 30, and one particle represents a three-dimensional vector (the elements areci, σi, and wij). The input and output vectors of RBFNN are all normalized. Given that the input data space can only span from 0 to 1, the RBF center ci can logically be constrained by the range [0,1]. According to the practice of the writers, the width σi and the connection weights wij are constrained by the ranges [1,30] and [0,100], respectively. Particle swarm is firstly initialized randomly within the given scope of the parameters. Then the particles move according to the flow chart that is shown in Figure 3. The end conditions are that the number of iterations reaches 100 or until the mean squared error is less than or equal to 0.0001.


## 5. Application: A Case Study of the Nuozhadu Dam

### 5.1. Dam Technical Properties

The Nuozhadu dam site is located on the Lancang River approximately 350 km southwest of Kunming, China. The dam site is located in a V-shaped valley surrounded by magnificent hills. The site geology is uniform and mainly consists of granite and T2M1 sand shale. The maximum cross section is shown in Figure4.Figure 4
The maximum cross section of the Nuozhadu core rockfill dam.As is shown in Figure4, the Nuozhadu dam is a rockfill dam with a mixed earth material and a gravel core. The elevation from the lowest foundation surface is 261.5 meters, and the reservoir volume is 237.03 × 108 m3. As shown in Figure 4, the dam body consists of 6 material partitions: core-wall, filter material I, filter material II, fine rockfill material, coarse rockfill material I, and coarse rockfill material II. The core-wall is built with a 10 m width at the top, upstream slope and downstream slope = 1 : 0.2, the top elevation = 820.5 m, and the lowest foundation surface elevation = 560.0 m. To prevent concentrated leakage and contact scouring, a reinforced concrete cushion was used at the contact areas among the core-wall and bedrocks, filter material I and bedrocks, and filter material II and bedrocks. Moreover, a 2 m thick layer of high plasticity clay was set between the cushion and the core-wall. The bottom of the cushion was treated with consolidated grouting [40, 41].


### 5.2. Permeability Coefficients of the Nuozhadu Dam

Considering the compacted construction process of the core-wall and the influence of the material gravity, this paper suggests that the vertical permeability coefficient of the core-wall is smaller than the horizontal permeability coefficient [42]. Thus, the core-wall material was regarded as anisotropic. The permeability values of the body and foundation were based on laboratory and in situ experimental results. Based on these results, the permeability values of the filter material II and the fine rockfill material are both variables ranging from 5 × 10−2 to 1 cm/s, which is 5 orders of magnitude larger than the core-wall; the permeability of the coarse rockfill I material is a variable ranging from 0.1 to 2 cm/s, which is 6 orders of magnitude larger than the core-wall; permeability values of the filter material I and the coarse rockfill material II are variables ranging from 5 × 10−4 to 1 × 10−2 cm/s, which is 3 orders of magnitude larger than the core-wall. This proves that the permeability values of the filter material II, the fine rockfill material, and the coarse rockfill material I are 5 to 6 orders of magnitude larger than the core-wall; thus, the permeability values of these three materials will only have a marginal effect on the seepage flow of the dam. Therefore, the core-wall (anisotropic), the filter material I, and the coarse rockfill material II, which are considered to be isotropic materials, are selected for the back analysis of the permeability coefficients.Permeability values of the materials in different parts of the body and foundation of the Nuozhadu dam are presented in Table1. The permeability values shown in Table 1 were determined based on laboratory and in situ experimental results on the dam and foundation materials.Table 1
Permeability coefficient of the Nuozhadu dam (cm/s).
Material
Range of permeability coefficient
Material
Value of permeability coefficient
Core (x-direction)
1 × 10−7~1 × 10−5
Filter material
0.2
Core (y-direction)
1 × 10−8~1 × 10−5
Fine rockfill material
0.2
Filter material I
5 × 10−3~5 × 10−2
Coarse rockfill material I
1
Coarse rockfill material II
5 × 10−3~5 × 10−2
Cofferdam
1
Weakly weathered rock
5 × 10−5
Cutoff grouting
5 × 10−6
Slightly weathered rock
1 × 10−5


### 5.3. Measured Hydraulic Head

In the body of the dam, 7 osmometers are distributed in the core-wall (DB-C-P-27, DB-C-P-35, and DB-C-P-43), filter material I (DB-C-P-24, DB-C-P-33), and coarse rockfill material II (DB-C-P-21, DB-C-P-40) and have continuous and relatively regular measurements that are used for the back analysis of the permeability coefficients. The upstream and downstream time-varying reservoir water lever curves, filling process curves of some partitions, and measured time-varying hydraulic head curves of the selected osmometers are shown in Figure5.Figure 5
Filling process curves of some partitions and measured time-varying hydraulic head curves.As is shown in Figure5, the measurement points of the core-wall (DB-C-P-27, DB-C-P-35, and DB-C-P-43) have high excess pore-water pressures before the storage, while the measured hydraulic heads of the measurement points distributed in the filter material I and coarse rockfill material II have synchronized changes with the upstream reservoir level. This suggests that the excess pore-water pressure arose because of the compaction procedure during construction of the core-wall and not because of a head difference between the upstream and downstream. Consequently, the excess pore-water pressure measured when core-wall filling was completed, which is shown in Table 2, should be deducted from the measured value in the back analysis of the core-wall permeability coefficients. The correction time-varying hydraulic head curves are shown in Figure 6.Table 2
Excess pore-water pressure head of measure points distributed in the core-wall.
Osmometer number
DB-C-P-43
DB-C-P-35
DB-C-P-27
Elevation (m)
738.11
701.78
659.40
Measured hydraulic head on December 21, 2012 (m)
787.50
793.69
774.72
Excess pore water pressure head (m)
45.39
91.91
115.32Figure 6
Measured and correction time-varying hydraulic head curves.In this paper, the permeability coefficients of the dam materials are retrieved using the average hydraulic head (as shown in Table3) measured under the condition of steady seepage (2012.12.21~2013.6.20). Meanwhile, the parameters are tested by measuring the hydraulic heads (as shown in Table 4) for the condition of an unstable saturated-unsaturated seepage while the water lever rose (2013.6.21~2013.10.21).Table 3
Average hydraulic head under steady seepage condition.
Osmometer number
DB-C-P-27
DB-C-P-35
DB-C-P-43
DB-C-P-24
DB-C-P-33
DB-C-P-21
DB-C-P-40
Excess pore water pressure (m)
115.32
91.91
45.39
—
—
—
—
Measured hydraulic head (m)
763.35
783.97
777.49
772.91
773.12
773.80
773.39
Correction hydraulic head (m)
648.03
692.06
732.10
772.91
773.12
773.80
773.39Table 4
Measured hydraulic head under unsteady seepage condition.
Osmometer number
DB-C-P-27
DB-C-P-35
DB-C-P-43
DB-C-P-24
DB-C-P-33
DB-C-P-21
DB-C-P-40
Excess pore water pressure (m)
115.32
91.91
43.39
—
—
—
—
August 19, 2013
Measured hydraulic head (m)
763.15
790.68
787.84
788.35
789.63
789.06
792.35
Correction hydraulic head (m)
647.83
698.77
742.45
788.35
789.63
789.06
792.35
September 16, 2013
Measured hydraulic head (m)
770.66
798.14
796.32
801.79
803.53
801.20
804.85
Correction hydraulic head (m)
655.34
706.23
750.93
801.79
803.53
801.20
804.85


### 5.4. Back Analysis of the Permeability Coefficients

The flow of water through soil is one of the fundamental issues in geotechnical and geoenvironmental engineering. In engineering, the more important issue is the pore-water pressure. The pore-water pressure, whether positive or negative, has a direct bearing on the shear strength and volume change characteristics of the soil. This paper needs to simulate the saturated-unsaturated seepage field. It is no longer acceptable to take a simplified approach and ignore unsaturated flow above the phreatic surface. It is mandatory to deal with unsaturated flow in typical situations such as transient flow problems. Modeling the flow of water through soil with a numerical solution can be very complex. Natural soil deposits are generally highly heterogeneous and nonisotropic. In addition, boundary conditions often change with time and cannot always be defined with certainty at the beginning of an analysis. Furthermore, when a soil becomes unsaturated, the coefficient of permeability becomes a function of the negative pore-water pressure in the soil. The pore-water pressure is the primary unknown and needs to be determined, so iterative numerical techniques are required to match the computed pore-water pressure and the material property, which makes the solution highly nonlinear. These complexities make it necessary to use some form of numerical analysis to analyze seepage problems for all, but the simplest cases [36]. This paper uses commercial software SEEP/W which can consider saturated-unsaturated unsteady seepage in numerical modeling.Broadly speaking, there are three main parts to a finite element analysis. The first is discretization: dividing the domain into small areas called elements. The second part is specifying and assigning material properties. The third is specifying and applying boundary conditions [36]. The two-dimensional finite element (FEM) method is used to calculate and analyze the seepage field. The boundary of the calculation domain is as follows. In the x-direction, considering the upstream cofferdam bottom is very wide, a quarter of the height of the dam is chosen in the upstream dam foundation, while a distance equal to the height of the dam is chosen in the downstream dam foundation. In the z-direction, the depth of the dam foundation is the same as the elevation. According to the real condition of the material partitions of the dam body and its foundation, the finite element mesh is set up using the automatic meshing method for the dam. After discretization, the finite element mesh of the dam and its foundation is shown in Figure 7, and the solid black panes indicate the measuring points. The steps of permeability coefficients inversion are as follows.Figure 7
2D finite element meshes and measure points of Nuozhadu dam.Step 1 (generate 50 sets of samples of permeability coefficient).
In this paper, based on the good-point set theory, the number-theoretic method (NTM) is proposed to establish normalized permeability coefficients with good diversity by a uniform design [43]. The program of the NTM is written in a MATLAB environment, and 50 sets of normalized permeability coefficients are generated (as shown in Table 5).Table 5
Partial normalized permeability coefficient.
Samples
Partition
Core-wall (x-direction)
Core-wall (y-direction)
Filter material I
Coarse rockfill material II
1
0.979
0.654
0.915
0.811
2
0.957
0.322
0.830
0.622
3
0.936
0.102
0.746
0.433
4
0.915
0.615
0.661
0.244
5
0.893
0.304
0.576
0.054
6
0.872
0.169
0.491
0.865
7
0.851
0.575
0.406
0.676
8
0.829
0.285
0.322
0.487
9
0.808
0.091
0.237
0.298
10
0.787
0.614
0.152
0.109
11
0.766
0.266
0.067
0.920
12
0.744
0.118
0.982
0.731
13
0.723
0.494
0.898
0.542
14
0.702
0.247
0.813
0.353
15
0.680
0.082
0.728
0.163
⋮
⋮
⋮
⋮
⋮
47
0.147
0.106
0.608
0.436
48
0.126
0.049
0.523
0.247
49
0.104
0.016
0.438
0.058
50
0.084
0.060
0.353
0.869Combining the ranges of the inversed permeability coefficients (as shown in Table1) with the 50 sets of normalized data, 50 sets of permeability coefficients evenly distributed using the reverse-normalization equation are obtained, which are given as(19)k=kmin⁡+kmax⁡-kmin⁡x,where kmax⁡ and kmin⁡ are the maximum and minimum values of k, respectively, and x is a normalized value (as is shown in Table 2).Step 2 (generate samples).
As shown in Figure5, the upstream and downstream water level had almost no change from December 2012 to June 2013 and the average water level (upstream water level = 774.02 m; downstream water level = 601.42 m) is taken as the hydraulic boundary condition. Then, 50 sets of reverse-normalized permeability coefficients are entered into the finite element model and the seepage field is simulated. After that, the hydraulic heads at the measuring points are extracted and normalized. Then they are taken as input vectors of the PSO-RBFNN. Partial samples are shown in Table 6.Table 6
Partial samples from the numerical simulation (m).
Partition
Osmometer number
Hydraulic head (m)
1
2
3
4
⋯
48
49
50
Core-wall
DB-C-P-27
645.44
645.54
645.66
645.59
⋯
645.32
645.58
645.69
DB-C-P-35
687.86
687.84
687.84
687.82
⋯
687.81
687.82
687.85
DB-C-P-43
731.13
730.96
730.81
730.83
⋯
731.13
730.16
730.80
Filter material I
DB-C-P-24
773.98
773.99
774.00
774.02
⋯
774.01
774.02
773.99
DB-C-P-33
773.96
773.97
773.98
774.01
⋯
774.00
774.01
773.97
Coarse rockfill material II
DB-C-P-21
774.00
774.00
774.01
774.02
⋯
774.01
774.02
774.00
DB-C-P-40
774.01
774.01
774.02
774.02
⋯
774.02
774.02
774.01Step 3 (establish the RBF neural network optimized by the PSO algorithm).
(a) Partial parameters initialization: the number of particlesN=30; acceleration coefficients c1=c2=2; inertia weight ω=1; the mean squared error goal MSE=0.0001; and the maximum number of iterations nmax⁡=100. Particle swarm is firstly initialized randomly within the given scope of the parameters; that is, ci∈[0,1], σi∈[1,30], and wij∈[0,100].
(b) Generally, in order to test the generalization capacity of feed forward neural networks, the number of testing samples should not be less than 20% of the number of training samples. Thus, taking 40 sets of normalized hydraulic heads as input vectors and the associated normalized permeability coefficients for the hydraulic heads as the output vectors, the prediction model of the PSO-RBFNN, introduced in part 4, is established. The seven-dimensional input normalized hydraulic heads vectors are being mapped to four-dimensional output normalized permeability coefficients vectors. The number of neurons in the hidden layer is problem dependent. In this study, when the number of neurons is added to twenty eight, the mean squared error between actual and predicted normalized permeability coefficients falls beneath an error goal (1 × 10−4). The PSO algorithm fitness function is (18) and while the number of neurons is added to twenty eight, the best particle fitness value converges following the iteration shown in Figure 8. Using the trained PSO-RBFNN, take the last 10 sets of normalized hydraulic heads as input vectors to predict the normalized permeability coefficients, and calculate the mean squared error between the actual and predicted normalized permeability coefficients. The testing result shows the trained PSO-RBFNN hold a satisfactory generalization.Figure 8
The best particle fitness value converging following the iteration.Step 4 (inverse permeability coefficient).
Choosing the correction average hydraulic head that is shown in Table3 as input vector, the PSO-RBFNN is used to inverse the permeability coefficient. The predicted results are listed in Table 7.Table 7
Back analysis results and indoor testing results of the permeability coefficients.
Partition
Core-wall (x-direction)
Core-wall (y-direction)
Filter material I
Coarse rockfill material II
k
Prediction
1.85 × 10−6
3.92 × 10−7
1.71 × 10−2
3.79 × 10−2
k
experiment
5 × 10−6
1.51 × 10−2
2.97 × 10−2Step 5 (evaluate the accuracy of the inversed permeability coefficient).
Choose the steady seepage field (during the period from December 2012 to June 2013) to be the initial seepage field with the hydraulic boundary condition: upstream water level = 774.02 m and downstream water level = 601.42 m; the model simulated the unsteady saturated-unsaturated seepage field during the period from June 2013 to September 2013. The hydraulic boundary condition is shown in Figure5. The dam materials’ nonsaturated property should be considered, except the cofferdam, antiseepage curtain, and dam foundation, which are all underwater. Generally, the soil-water characteristic curve and permeability coefficient curve are representative of the materials’ nonsaturated property. The soil-water characteristic curve is estimated using the saturated volumetric water content and a spline function. The permeability coefficient curve is estimated using the saturated permeability coefficient, soil-water characteristic curve, and Fredlund-Xing model [44, 45]. Based on the required saturated volumetric water content saturated permeability coefficient listed in Table 8, the fitting soil-water characteristic curve and permeability coefficient curve are shown in Figures 9 and 10. The simulated and measured hydraulic head of the steady and unsteady seepage field (the selected moment: August 19, 2013, and September 16, 2013) at measuring points are listed in Table 9, which also contains absolute error and relative error of the simulated and measured hydraulic head. Figure 11 shows the comparison diagram of the measured and simulated hydraulic head for the condition of steady and unsteady seepage. Figures 12 and 13 are histograms of the absolute error and relative error of the simulated and measured hydraulic head.Table 8
Unsteady saturated-unsaturated seepage parameters of the Nuozhadu dam.
Partition
θ
saturated
K
saturated (cm/s)
Materials
θ
saturated
K
saturated (cm/s)
Core-wall (x-direction)
0.22
1.97 × 10−6
Fine rockfill material
0.32
0.2
Core-wall (y-direction)
0.22
3.81 × 10−7
Coarse rockfill material I
0.28
1
Filter material I
0.37
1.68 × 10−2
Coarse rockfill material II
0.25
3.85 × 10−2
Filter material II
0.41
0.2
Coffer dam
0.28
1Table 9
Measured and simulated hydraulic head under the condition of steady and unsteady seepage.
Partition
Coarse rockfill material II
Filter material I
Core-wall
Osmometer number
DB-C-P-21
DB-C-P-40
DB-C-P-24
DB-C-P-33
DB-C-P-43
DB-C-P-35
DB-C-P-27
Abscissa (m)
529.9
629.0
629.9
638.1
669.5
680.0
698.5
Steady seepage
Measured hydraulic head (m)
773.80
773.39
772.91
773.12
732.10
692.06
648.03
Simulated hydraulic head (m)
774.00
773.99
773.97
773.98
730.01
687.80
644.67
Absolute error (m)
0.20
0.60
1.06
0.86
−2.09
−4.26
−3.36
Relative error (%)
0.09
0.28
0.50
0.40
−1.22
−3.23
−3.82
Unsteady seepage (on August 19, 2013)
Measured hydraulic head (m)
789.06
792.35
788.35
789.63
742.45
698.77
647.83
Simulated hydraulic head (m)
790.26
790.25
790.21
790.22
746.17
695.72
648.73
Absolute error (m)
1.20
−2.10
1.86
0.59
3.72
−3.05
0.90
Relative error (%)
0.52
−0.90
0.82
0.26
2.03
−2.20
1.03
Unsteady Seepage (on September 16, 2013)
Measured hydraulic head (m)
801.20
804.85
801.79
803.53
750.93
706.23
655.34
Simulated hydraulic head (m)
802.86
802.85
802.80
802.81
752.97
701.27
650.86
Absolute error (m)
1.66
−2.00
1.10
−0.72
2.04
−4.96
−4.48
Relative error (%)
0.69
−0.82
0.42
−0.30
1.07
−3.39
−4.70
Note: absolute error = simulated hydraulic head − measured hydraulic head; relative error = absolute error/(measured hydraulic head − the lowest elevation of foundation surface, 560 m).Figure 9
Soil-water characteristic curve.Figure 10
Permeability coefficient curve.Figure 11
Measured and simulated hydraulic head under the condition of steady and unsteady seepage.Figure 12
Bar graph of absolute error.Figure 13
Bar graph of relative error.As shown in Figure11, the measured and simulated hydraulic head agree well and reflect basically the distribution and variation of the steady and unsteady seepage fields in the dam, which validate the proposed model and inversion method. Table 9 shows that, for the three different boundary conditions of the seepage field, the maximum absolute error between the simulated and measured hydraulic head is 4.96 m and more than 70% of the absolute error is less than 2.1 m; the maximum relative error is 4.7% and more than 70% of relative error is less than 1.1%; moreover, the average absolute relative error is 1.37%, which strongly indicates that the inversed permeability coefficients are reasonable.


## 5.1. Dam Technical Properties

The Nuozhadu dam site is located on the Lancang River approximately 350 km southwest of Kunming, China. The dam site is located in a V-shaped valley surrounded by magnificent hills. The site geology is uniform and mainly consists of granite and T2M1 sand shale. The maximum cross section is shown in Figure4.Figure 4
The maximum cross section of the Nuozhadu core rockfill dam.As is shown in Figure4, the Nuozhadu dam is a rockfill dam with a mixed earth material and a gravel core. The elevation from the lowest foundation surface is 261.5 meters, and the reservoir volume is 237.03 × 108 m3. As shown in Figure 4, the dam body consists of 6 material partitions: core-wall, filter material I, filter material II, fine rockfill material, coarse rockfill material I, and coarse rockfill material II. The core-wall is built with a 10 m width at the top, upstream slope and downstream slope = 1 : 0.2, the top elevation = 820.5 m, and the lowest foundation surface elevation = 560.0 m. To prevent concentrated leakage and contact scouring, a reinforced concrete cushion was used at the contact areas among the core-wall and bedrocks, filter material I and bedrocks, and filter material II and bedrocks. Moreover, a 2 m thick layer of high plasticity clay was set between the cushion and the core-wall. The bottom of the cushion was treated with consolidated grouting [40, 41].


## 5.2. Permeability Coefficients of the Nuozhadu Dam

Considering the compacted construction process of the core-wall and the influence of the material gravity, this paper suggests that the vertical permeability coefficient of the core-wall is smaller than the horizontal permeability coefficient [42]. Thus, the core-wall material was regarded as anisotropic. The permeability values of the body and foundation were based on laboratory and in situ experimental results. Based on these results, the permeability values of the filter material II and the fine rockfill material are both variables ranging from 5 × 10−2 to 1 cm/s, which is 5 orders of magnitude larger than the core-wall; the permeability of the coarse rockfill I material is a variable ranging from 0.1 to 2 cm/s, which is 6 orders of magnitude larger than the core-wall; permeability values of the filter material I and the coarse rockfill material II are variables ranging from 5 × 10−4 to 1 × 10−2 cm/s, which is 3 orders of magnitude larger than the core-wall. This proves that the permeability values of the filter material II, the fine rockfill material, and the coarse rockfill material I are 5 to 6 orders of magnitude larger than the core-wall; thus, the permeability values of these three materials will only have a marginal effect on the seepage flow of the dam. Therefore, the core-wall (anisotropic), the filter material I, and the coarse rockfill material II, which are considered to be isotropic materials, are selected for the back analysis of the permeability coefficients.Permeability values of the materials in different parts of the body and foundation of the Nuozhadu dam are presented in Table1. The permeability values shown in Table 1 were determined based on laboratory and in situ experimental results on the dam and foundation materials.Table 1
Permeability coefficient of the Nuozhadu dam (cm/s).
Material
Range of permeability coefficient
Material
Value of permeability coefficient
Core (x-direction)
1 × 10−7~1 × 10−5
Filter material
0.2
Core (y-direction)
1 × 10−8~1 × 10−5
Fine rockfill material
0.2
Filter material I
5 × 10−3~5 × 10−2
Coarse rockfill material I
1
Coarse rockfill material II
5 × 10−3~5 × 10−2
Cofferdam
1
Weakly weathered rock
5 × 10−5
Cutoff grouting
5 × 10−6
Slightly weathered rock
1 × 10−5


## 5.3. Measured Hydraulic Head

In the body of the dam, 7 osmometers are distributed in the core-wall (DB-C-P-27, DB-C-P-35, and DB-C-P-43), filter material I (DB-C-P-24, DB-C-P-33), and coarse rockfill material II (DB-C-P-21, DB-C-P-40) and have continuous and relatively regular measurements that are used for the back analysis of the permeability coefficients. The upstream and downstream time-varying reservoir water lever curves, filling process curves of some partitions, and measured time-varying hydraulic head curves of the selected osmometers are shown in Figure5.Figure 5
Filling process curves of some partitions and measured time-varying hydraulic head curves.As is shown in Figure5, the measurement points of the core-wall (DB-C-P-27, DB-C-P-35, and DB-C-P-43) have high excess pore-water pressures before the storage, while the measured hydraulic heads of the measurement points distributed in the filter material I and coarse rockfill material II have synchronized changes with the upstream reservoir level. This suggests that the excess pore-water pressure arose because of the compaction procedure during construction of the core-wall and not because of a head difference between the upstream and downstream. Consequently, the excess pore-water pressure measured when core-wall filling was completed, which is shown in Table 2, should be deducted from the measured value in the back analysis of the core-wall permeability coefficients. The correction time-varying hydraulic head curves are shown in Figure 6.Table 2
Excess pore-water pressure head of measure points distributed in the core-wall.
Osmometer number
DB-C-P-43
DB-C-P-35
DB-C-P-27
Elevation (m)
738.11
701.78
659.40
Measured hydraulic head on December 21, 2012 (m)
787.50
793.69
774.72
Excess pore water pressure head (m)
45.39
91.91
115.32Figure 6
Measured and correction time-varying hydraulic head curves.In this paper, the permeability coefficients of the dam materials are retrieved using the average hydraulic head (as shown in Table3) measured under the condition of steady seepage (2012.12.21~2013.6.20). Meanwhile, the parameters are tested by measuring the hydraulic heads (as shown in Table 4) for the condition of an unstable saturated-unsaturated seepage while the water lever rose (2013.6.21~2013.10.21).Table 3
Average hydraulic head under steady seepage condition.
Osmometer number
DB-C-P-27
DB-C-P-35
DB-C-P-43
DB-C-P-24
DB-C-P-33
DB-C-P-21
DB-C-P-40
Excess pore water pressure (m)
115.32
91.91
45.39
—
—
—
—
Measured hydraulic head (m)
763.35
783.97
777.49
772.91
773.12
773.80
773.39
Correction hydraulic head (m)
648.03
692.06
732.10
772.91
773.12
773.80
773.39Table 4
Measured hydraulic head under unsteady seepage condition.
Osmometer number
DB-C-P-27
DB-C-P-35
DB-C-P-43
DB-C-P-24
DB-C-P-33
DB-C-P-21
DB-C-P-40
Excess pore water pressure (m)
115.32
91.91
43.39
—
—
—
—
August 19, 2013
Measured hydraulic head (m)
763.15
790.68
787.84
788.35
789.63
789.06
792.35
Correction hydraulic head (m)
647.83
698.77
742.45
788.35
789.63
789.06
792.35
September 16, 2013
Measured hydraulic head (m)
770.66
798.14
796.32
801.79
803.53
801.20
804.85
Correction hydraulic head (m)
655.34
706.23
750.93
801.79
803.53
801.20
804.85


## 5.4. Back Analysis of the Permeability Coefficients

The flow of water through soil is one of the fundamental issues in geotechnical and geoenvironmental engineering. In engineering, the more important issue is the pore-water pressure. The pore-water pressure, whether positive or negative, has a direct bearing on the shear strength and volume change characteristics of the soil. This paper needs to simulate the saturated-unsaturated seepage field. It is no longer acceptable to take a simplified approach and ignore unsaturated flow above the phreatic surface. It is mandatory to deal with unsaturated flow in typical situations such as transient flow problems. Modeling the flow of water through soil with a numerical solution can be very complex. Natural soil deposits are generally highly heterogeneous and nonisotropic. In addition, boundary conditions often change with time and cannot always be defined with certainty at the beginning of an analysis. Furthermore, when a soil becomes unsaturated, the coefficient of permeability becomes a function of the negative pore-water pressure in the soil. The pore-water pressure is the primary unknown and needs to be determined, so iterative numerical techniques are required to match the computed pore-water pressure and the material property, which makes the solution highly nonlinear. These complexities make it necessary to use some form of numerical analysis to analyze seepage problems for all, but the simplest cases [36]. This paper uses commercial software SEEP/W which can consider saturated-unsaturated unsteady seepage in numerical modeling.Broadly speaking, there are three main parts to a finite element analysis. The first is discretization: dividing the domain into small areas called elements. The second part is specifying and assigning material properties. The third is specifying and applying boundary conditions [36]. The two-dimensional finite element (FEM) method is used to calculate and analyze the seepage field. The boundary of the calculation domain is as follows. In the x-direction, considering the upstream cofferdam bottom is very wide, a quarter of the height of the dam is chosen in the upstream dam foundation, while a distance equal to the height of the dam is chosen in the downstream dam foundation. In the z-direction, the depth of the dam foundation is the same as the elevation. According to the real condition of the material partitions of the dam body and its foundation, the finite element mesh is set up using the automatic meshing method for the dam. After discretization, the finite element mesh of the dam and its foundation is shown in Figure 7, and the solid black panes indicate the measuring points. The steps of permeability coefficients inversion are as follows.Figure 7
2D finite element meshes and measure points of Nuozhadu dam.Step 1 (generate 50 sets of samples of permeability coefficient).
In this paper, based on the good-point set theory, the number-theoretic method (NTM) is proposed to establish normalized permeability coefficients with good diversity by a uniform design [43]. The program of the NTM is written in a MATLAB environment, and 50 sets of normalized permeability coefficients are generated (as shown in Table 5).Table 5
Partial normalized permeability coefficient.
Samples
Partition
Core-wall (x-direction)
Core-wall (y-direction)
Filter material I
Coarse rockfill material II
1
0.979
0.654
0.915
0.811
2
0.957
0.322
0.830
0.622
3
0.936
0.102
0.746
0.433
4
0.915
0.615
0.661
0.244
5
0.893
0.304
0.576
0.054
6
0.872
0.169
0.491
0.865
7
0.851
0.575
0.406
0.676
8
0.829
0.285
0.322
0.487
9
0.808
0.091
0.237
0.298
10
0.787
0.614
0.152
0.109
11
0.766
0.266
0.067
0.920
12
0.744
0.118
0.982
0.731
13
0.723
0.494
0.898
0.542
14
0.702
0.247
0.813
0.353
15
0.680
0.082
0.728
0.163
⋮
⋮
⋮
⋮
⋮
47
0.147
0.106
0.608
0.436
48
0.126
0.049
0.523
0.247
49
0.104
0.016
0.438
0.058
50
0.084
0.060
0.353
0.869Combining the ranges of the inversed permeability coefficients (as shown in Table1) with the 50 sets of normalized data, 50 sets of permeability coefficients evenly distributed using the reverse-normalization equation are obtained, which are given as(19)k=kmin⁡+kmax⁡-kmin⁡x,where kmax⁡ and kmin⁡ are the maximum and minimum values of k, respectively, and x is a normalized value (as is shown in Table 2).Step 2 (generate samples).
As shown in Figure5, the upstream and downstream water level had almost no change from December 2012 to June 2013 and the average water level (upstream water level = 774.02 m; downstream water level = 601.42 m) is taken as the hydraulic boundary condition. Then, 50 sets of reverse-normalized permeability coefficients are entered into the finite element model and the seepage field is simulated. After that, the hydraulic heads at the measuring points are extracted and normalized. Then they are taken as input vectors of the PSO-RBFNN. Partial samples are shown in Table 6.Table 6
Partial samples from the numerical simulation (m).
Partition
Osmometer number
Hydraulic head (m)
1
2
3
4
⋯
48
49
50
Core-wall
DB-C-P-27
645.44
645.54
645.66
645.59
⋯
645.32
645.58
645.69
DB-C-P-35
687.86
687.84
687.84
687.82
⋯
687.81
687.82
687.85
DB-C-P-43
731.13
730.96
730.81
730.83
⋯
731.13
730.16
730.80
Filter material I
DB-C-P-24
773.98
773.99
774.00
774.02
⋯
774.01
774.02
773.99
DB-C-P-33
773.96
773.97
773.98
774.01
⋯
774.00
774.01
773.97
Coarse rockfill material II
DB-C-P-21
774.00
774.00
774.01
774.02
⋯
774.01
774.02
774.00
DB-C-P-40
774.01
774.01
774.02
774.02
⋯
774.02
774.02
774.01Step 3 (establish the RBF neural network optimized by the PSO algorithm).
(a) Partial parameters initialization: the number of particlesN=30; acceleration coefficients c1=c2=2; inertia weight ω=1; the mean squared error goal MSE=0.0001; and the maximum number of iterations nmax⁡=100. Particle swarm is firstly initialized randomly within the given scope of the parameters; that is, ci∈[0,1], σi∈[1,30], and wij∈[0,100].
(b) Generally, in order to test the generalization capacity of feed forward neural networks, the number of testing samples should not be less than 20% of the number of training samples. Thus, taking 40 sets of normalized hydraulic heads as input vectors and the associated normalized permeability coefficients for the hydraulic heads as the output vectors, the prediction model of the PSO-RBFNN, introduced in part 4, is established. The seven-dimensional input normalized hydraulic heads vectors are being mapped to four-dimensional output normalized permeability coefficients vectors. The number of neurons in the hidden layer is problem dependent. In this study, when the number of neurons is added to twenty eight, the mean squared error between actual and predicted normalized permeability coefficients falls beneath an error goal (1 × 10−4). The PSO algorithm fitness function is (18) and while the number of neurons is added to twenty eight, the best particle fitness value converges following the iteration shown in Figure 8. Using the trained PSO-RBFNN, take the last 10 sets of normalized hydraulic heads as input vectors to predict the normalized permeability coefficients, and calculate the mean squared error between the actual and predicted normalized permeability coefficients. The testing result shows the trained PSO-RBFNN hold a satisfactory generalization.Figure 8
The best particle fitness value converging following the iteration.Step 4 (inverse permeability coefficient).
Choosing the correction average hydraulic head that is shown in Table3 as input vector, the PSO-RBFNN is used to inverse the permeability coefficient. The predicted results are listed in Table 7.Table 7
Back analysis results and indoor testing results of the permeability coefficients.
Partition
Core-wall (x-direction)
Core-wall (y-direction)
Filter material I
Coarse rockfill material II
k
Prediction
1.85 × 10−6
3.92 × 10−7
1.71 × 10−2
3.79 × 10−2
k
experiment
5 × 10−6
1.51 × 10−2
2.97 × 10−2Step 5 (evaluate the accuracy of the inversed permeability coefficient).
Choose the steady seepage field (during the period from December 2012 to June 2013) to be the initial seepage field with the hydraulic boundary condition: upstream water level = 774.02 m and downstream water level = 601.42 m; the model simulated the unsteady saturated-unsaturated seepage field during the period from June 2013 to September 2013. The hydraulic boundary condition is shown in Figure5. The dam materials’ nonsaturated property should be considered, except the cofferdam, antiseepage curtain, and dam foundation, which are all underwater. Generally, the soil-water characteristic curve and permeability coefficient curve are representative of the materials’ nonsaturated property. The soil-water characteristic curve is estimated using the saturated volumetric water content and a spline function. The permeability coefficient curve is estimated using the saturated permeability coefficient, soil-water characteristic curve, and Fredlund-Xing model [44, 45]. Based on the required saturated volumetric water content saturated permeability coefficient listed in Table 8, the fitting soil-water characteristic curve and permeability coefficient curve are shown in Figures 9 and 10. The simulated and measured hydraulic head of the steady and unsteady seepage field (the selected moment: August 19, 2013, and September 16, 2013) at measuring points are listed in Table 9, which also contains absolute error and relative error of the simulated and measured hydraulic head. Figure 11 shows the comparison diagram of the measured and simulated hydraulic head for the condition of steady and unsteady seepage. Figures 12 and 13 are histograms of the absolute error and relative error of the simulated and measured hydraulic head.Table 8
Unsteady saturated-unsaturated seepage parameters of the Nuozhadu dam.
Partition
θ
saturated
K
saturated (cm/s)
Materials
θ
saturated
K
saturated (cm/s)
Core-wall (x-direction)
0.22
1.97 × 10−6
Fine rockfill material
0.32
0.2
Core-wall (y-direction)
0.22
3.81 × 10−7
Coarse rockfill material I
0.28
1
Filter material I
0.37
1.68 × 10−2
Coarse rockfill material II
0.25
3.85 × 10−2
Filter material II
0.41
0.2
Coffer dam
0.28
1Table 9
Measured and simulated hydraulic head under the condition of steady and unsteady seepage.
Partition
Coarse rockfill material II
Filter material I
Core-wall
Osmometer number
DB-C-P-21
DB-C-P-40
DB-C-P-24
DB-C-P-33
DB-C-P-43
DB-C-P-35
DB-C-P-27
Abscissa (m)
529.9
629.0
629.9
638.1
669.5
680.0
698.5
Steady seepage
Measured hydraulic head (m)
773.80
773.39
772.91
773.12
732.10
692.06
648.03
Simulated hydraulic head (m)
774.00
773.99
773.97
773.98
730.01
687.80
644.67
Absolute error (m)
0.20
0.60
1.06
0.86
−2.09
−4.26
−3.36
Relative error (%)
0.09
0.28
0.50
0.40
−1.22
−3.23
−3.82
Unsteady seepage (on August 19, 2013)
Measured hydraulic head (m)
789.06
792.35
788.35
789.63
742.45
698.77
647.83
Simulated hydraulic head (m)
790.26
790.25
790.21
790.22
746.17
695.72
648.73
Absolute error (m)
1.20
−2.10
1.86
0.59
3.72
−3.05
0.90
Relative error (%)
0.52
−0.90
0.82
0.26
2.03
−2.20
1.03
Unsteady Seepage (on September 16, 2013)
Measured hydraulic head (m)
801.20
804.85
801.79
803.53
750.93
706.23
655.34
Simulated hydraulic head (m)
802.86
802.85
802.80
802.81
752.97
701.27
650.86
Absolute error (m)
1.66
−2.00
1.10
−0.72
2.04
−4.96
−4.48
Relative error (%)
0.69
−0.82
0.42
−0.30
1.07
−3.39
−4.70
Note: absolute error = simulated hydraulic head − measured hydraulic head; relative error = absolute error/(measured hydraulic head − the lowest elevation of foundation surface, 560 m).Figure 9
Soil-water characteristic curve.Figure 10
Permeability coefficient curve.Figure 11
Measured and simulated hydraulic head under the condition of steady and unsteady seepage.Figure 12
Bar graph of absolute error.Figure 13
Bar graph of relative error.As shown in Figure11, the measured and simulated hydraulic head agree well and reflect basically the distribution and variation of the steady and unsteady seepage fields in the dam, which validate the proposed model and inversion method. Table 9 shows that, for the three different boundary conditions of the seepage field, the maximum absolute error between the simulated and measured hydraulic head is 4.96 m and more than 70% of the absolute error is less than 2.1 m; the maximum relative error is 4.7% and more than 70% of relative error is less than 1.1%; moreover, the average absolute relative error is 1.37%, which strongly indicates that the inversed permeability coefficients are reasonable.


## 6. Conclusions

The problem presented in this study is in essence one of estimating dam permeability coefficients while minimizing the error between the observed and computed hydraulic head values. An essential requirement of the solution to this optimization problem is the fulfillment of the implicit nonlinear relation that exists between the permeability coefficients and hydraulic head as described by the dam seepage flow equation. In this paper, a radial basis function (RBF) neural network is applied to solve the optimization problem and produce an estimate of the permeability coefficients. The fact that the RBF neural network can be trained to provide a good data-driven approximation of the dam seepage flow equation means that the solution (upon its inversion) is able to accurately predict the dam permeability coefficients in response to the present hydraulic head. This fact highlights a need to pay special attention to the structure parameters of the RBF neural network. The particle swarm optimization (PSO) algorithm, which is easy to program and has a quick rate of convergence and the powerful ability of global search, is used to train the neural network to obtain the optimal values of the RBF neural network structure parameters.The proposed approach (RBF neural network optimized using the PSO algorithm) is applied to inverse the permeability coefficients of the core-wall, filter material I, and coarse rockfill material II of the Nuozhadu Dam. The results show that the measured hydraulic head basically agrees well with calculated hydraulic head under both the condition of steady seepage and that of unsteady seepage, which indicates that the proposed method has high computational efficiency and high identification accuracy. The results illustrate the feasibility of the proposed method in the back analysis of permeability coefficients.



---

*Source: 124042-2015-11-17.xml*