# A Deformation Prediction Model of High Arch Dams in the Initial Operation Period Based on PSR-SVM-IGWO

Authors: Mingjun Li; Jiangyang Pan; Yaolai Liu; Hao Liu; Junxing Wang; Zhou Zhao
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

The deformation prediction of the dam in the initial stage of operation is very important for the safety of high dams. A hybrid model integrating chaos theory, support vector machine (SVM), and an improved Grey Wolf Optimization (IGWO) algorithm is developed for deformation prediction of dam in the initial operation period. Firstly, the chaotic characteristics of the dam deformation time series will be identified, mainly using the Lyapunov exponent method, the correlation dimension method, and the Kolmogorov entropy method. Secondly, the SVM-IGWO model based on phase space reconstruction (PSR) is established for deformation forecasting of the dam in the initial operation period. Taking SVM as the core, the deformation time series is reconstructed in phase space to determine the input variables of SVM and the GWO algorithm is improved to realize the optimization of SVM parameters. Finally, take the actual monitoring displacement of Xiluodu super-high arch dam as an example. The engineering application example shows that, compared with the existing models, the prediction accuracy of the PSR-SVM-IGWO model established in this paper is improved.


---

## Body

## 1. Introduction

The establishment of an accurate deformation prediction model based on the prototype observation data is of great significance to the long-term safe operation of the dam [1]. The uncertainty in the dam construction process reduces the prediction accuracy of the numerical model (also known as the deterministic model) based on the finite element method (FEM) [2]. Moreover, conventional statistical models are difficult to adapt to the complex nonlinear relationship between multiple factors and effect sizes, and the accuracy of model predictions is hard to guarantee [3]. In addition, the existing prediction models only consider the main influencing factors, such as water pressure, temperature, and aging, but do not consider the chaotic components that may be included in the deformed time series, which further limits the improvement of fitting accuracy [4].It is difficult to establish a stable and reliable deformation prediction model in the initial stage of high dam operation, because the initial stage of operation is a temporary stage and the most uncertain period. Studies [5, 6] have shown that nearly 50% of dam failures occur in the first 5 years of operation, which emphasizes the necessity and importance of early statistical modelling methods for dam safety monitoring. The purpose of this paper is to propose a hybrid model based on chaos theory, SVM, and IGWO algorithm and apply it to the deformation prediction of the dam in the initial operation stage.Chaos is an irregular random behavior with initial sensitivity and ergodicity [7–9]. Some studies have shown that there is chaos in the measured displacement data of dams. The efficient extraction of the chaotic component contained in the observation data is of great practical significance for improving the accuracy of prediction models. Wei [7] added a residual prediction item to the conventional statistical model of dam deformation prediction and applied chaos theory to discuss the calculation method of the residual prediction item and thus proposed a new hybrid prediction model. Zhang [10] proposed a new method based on empirical mode decomposition and phase space reconstruction theory to analyze the time-varying characteristics of dams. Su [11] combined SVM with phase space reconstruction, wavelet analysis, particle swarm optimization, and other methods to establish a dam deformation prediction model. Obviously, when the phase space reconstruction method is used to design the input variable, the machine learning algorithm shows high applicability in improving the prediction accuracy. However, there is little research on the chaotic characteristics in the deformation observation data during the initial operation period of the dam.SVMs are widely used in dam behavior prediction because of their obvious advantages in solving low-sample high-dimensional nonlinear problems [12]. The SVM can transform the original nonlinear problem into a high-latitude linear problem through the kernel parameters, so it is well adapted to the nonlinear deformation prediction [13–15]. However, the selection of the kernel parameters and penalty factors of SVM will directly affect its performance in dealing with nonlinear problems.Grey Wolf Optimizer (GWO) is a new intelligent optimization algorithm proposed by Mirjalili with reference to the social hierarchy and hunting behavior of grey wolves [16, 17]. The algorithm realizes intelligent optimization through the process of tracking, encircling, hunting, and attacking the grey wolf population [17, 18]. Studies have shown that the GWO algorithm is better than the other evolutionary algorithms [17, 19–21] in terms of quality, speed, and stability of the final solution [17, 22, 23]. However, the possibility of premature convergence reduces the probability of the algorithm finding the global optimum. The initial population is an important factor influencing the optimization performance of an intelligent algorithm, but it is difficult to guarantee the diversity of the initial population with conventional GWO algorithms.Many improved GWO (IGWO) algorithms have been proposed to overcome the shortcomings of conventional GWO and have been widely used to solve engineering problems. Li [24] proposed an IGWO algorithm combining tent chaos, Gaussian perturbation, and cosine control factor and applied it to the solution of four typical complex engineering optimization problems: pressure vessel design, tension spring design, welded beam design, and three-truss design. Zhao [25] proposed an IGWO algorithm based on the hybrid update equation of SM algorithm and GWO algorithm and verified the better performance of this IGWO algorithm through simulation experiments. Wang [26] introduced biological evolution and natural SOF principles into the standard GWO algorithm, while Liu [27] proposed an IGWO algorithm based on the differential evolution (DE) algorithm and the OTSU algorithm. In the IGWO algorithm proposed by Mohammad [28], the movement strategy is developed by introducing the hunting (DLH) search strategy based on dimensional learning. Whether it is the performance test result of the simulation test function or the actual engineering problem solving, the IGWO algorithm proposed by the above scholars shows better convergence speed and optimization performance than other intelligent optimization algorithms, such as the DE algorithm, the PSO algorithm, and the ABC algorithm. However, there are still shortcomings in some aspects. For example, it cannot be effectively applied to all problems, and the convergence accuracy can be further improved.To optimize the optimal kernel parameters and penalty factors of SVMs, an improved GWO (IGWO) algorithm is proposed in this paper. The IGWO algorithm enriches the diversity of the initial population through differential mutation and uses the global optimal search ability of the standard GWO algorithm to ensure convergence characteristics. On the basis of identifying the chaotic characteristics of the deformation observation data, a hybrid deformation prediction model of the high arch dam at the initial stage of operation is established.The structure of this article is as follows: Section2 introduces a variety of methods to identify chaotic characteristics of deformation time series. In Section 3, a new dam deformation prediction model is proposed on the basis of phase space reconstruction. An improved GWO algorithm based on the fusion of the DE algorithm and GWO algorithm is introduced to optimize the parameter settings of SVM. In Section 4, taking the measured displacement data of Xiluodu super-high arch dam in the initial operation stage as an example, the method proposed in this paper and other common methods are used to establish the deformation prediction models. The prediction accuracies of models are compared and analyzed. Section 5 summarizes the main conclusions reached here.


## 2. Identification of Chaotic Characteristics

Chaotic systems are deterministic and sensitive to initial conditions [8]. It is worth noting that it is neither random nor disordered. At present, the identification of chaotic characteristics of time series is mainly based on phase space reconstruction, which can obtain more hidden information by recovering the chaotic attractor in the so-called high-dimensional phase space. The Lyapunov exponent, correlation dimension, and Kolmogorov entropy of the singular attractor are calculated to correctly distinguish the chaotic system from the random system [29].


### 2.1. Phase Space Reconstruction

The reconstruction of the phase space is the basis for the quantitative analysis of chaotic time series, in which the embedding dimension and the delay coordinate are the two most critical parameters [30]. The time delay method is currently the most commonly used method. For univariate chaotic time series xi,i=1,2,…,n,(1)yi=xi,xi+τ,…,xi+m−1τ,i=1,2,…,n−m−1τ.According to Takens’ Theorem, the appropriate choice of the embedding dimensionm and the delay time τ can restore the dynamics properties of the original state space in the sense of topological equivalence.


#### 2.1.1. Delay Time

The mutual information method [31] is introduced to determine the delay time τ of the measured displacement sequence of the dam, shown as follows:(2)Iτ=∑xi,xi+τPxi,xi+τlog2Pxi,xi+τPxixi+τ,where Pxi is the normalized distribution of xi, Pxi+τ is the normalized distribution of xi+τ, and Pxi,xi+τ is the joint distribution of xi and xi+τ. The time at which the first minimum point appears in the τ∼Iτ curve is often selected as the optimal delay time.


#### 2.1.2. Embedding Dimension

The embedding dimensionm is determined by the Cao method [32]. The distance at,m between the phase point and the nearest neighbor point is shown as follows:(3)at,m=ytm+1−yfm+1ytm−yfm,where ytm+1 is the t th vector in the reconstructed phase space with the embedding dimension m+1 and yfm+1 is the nearest neighbor to ytm+1. yfm is the nearest neighbor to ytm in the reconstructed phase space with the embedding dimension m.The average valueEm of at,m is calculated as follows:(4)Em=1n−mτ∑t=1n−mτat,m.Then, the change ofEm is as follows:(5)E1m=Em+1Em.Whenm>m0, if E1m no longer changes, m0 represents the minimum embedding dimension. In order to avoid the situation that the change of E1m is difficult to judge, the Cao method adds another definition.(6)E′m=1n−mτ∑t=1n−mτytm+1−yfm+1,E2m=E′m+1E′m.For chaotic time series, there will always be some value ofm, so that E2m≠1. Whether E1m tends to be stable and whether the value of m can achieve E2m≠1 are observed to determine the minimum m of the reconstructed phase space.


### 2.2. Lyapunov Exponent Method

The maximum Lyapunov exponent (λmax) is usually regarded as an indicator of chaotic motion [33]. λmax>0 indicates that the system is in a chaotic state. The specific process of calculating λmax by the small data amount method is as follows.Firstly, determine the nearest neighbor pointyj^ of the phase point yj, and its short-term distance is dj0=yj−yj⌢.Secondly, calculate the distancedji=yj+i−yj⌢+i between yj and yj^ after the i th discrete time steps.Thirdly, get the average of alllndji, which is yi:(7)yi=1qΔt∑t=1qlndji,where q represents the number of i satisfying the condition dji≠0 and Δt represents the time interval.Finally, calculateλmax. λmax is the slope of the regression line obtained by the least-squares method.


### 2.3. Correlation Dimension Method

The correlation dimensionD2 is mainly determined by the Grassberger-Procaccia algorithm [34]. Suppose that r is the radius of the sphere centered on yi and yj; then the correlation integral Cr is given by(8)Cr=limn⟶∞2nn−1∑i,j=1nθr−yi−yj,where θ⋅ is the Heaviside function:(9)θu=0,u≤0,1,u>0,limr⟶0Cr∞rDr⟶0,where D2 is the correlation dimension:(10)D2=logCnrlogr.Thus, drawlogCnr/logr curve, and then the value of D2 can be determined according to the slope of the curve. As the embedding dimension m gradually increases, the slope of the curve converges, and the limit of convergence is the correlation dimension D2. The slope of the curve of a stochastic system will continue to increase with the increase of the embedding dimension, and there will be no convergence phenomenon.


### 2.4. Kolmogorov Entropy Method

Kolmogorov entropy [35] describes the generation rate of chaotic orbital information over time. Kolmogorov entropy reflects the chaos level of nonlinear dynamic systems, and the Kolmogorov entropy proposed by Grassberger and Procaccia is most commonly used as its estimate. The Kolmogorov entropy could be estimated by the correlation integral method.(11)K2=−limr⟶0limd⟶∞1Δmτlog2CmrCm+Δmr.When the embedding dimensionm is continuously increasing at intervals of Δm, the stable estimation of Kolmogorov entropy can be obtained through the equal slope regression of equation (11). It should be noted that the minimum value of m must be an integer greater than the correlation dimension D2.In a phase space with embedding dimensioni, there is(12)xij=log2rij,yij=log2Crij,where(13)yij=axij+b.Leta=D2, and, for the embedding dimensions i and i+Δm, there is(14)K2=limi⟶∞ΔbiΔmτ,where Δb=bi−bi+m.The Kolmogorov entropy estimateK2 can be used to judge the motion properties of the nonlinear system: K2=0 means the nonlinear system performs regular motion, K2>0 means the nonlinear system performs chaotic motion, and K2<0 means the nonlinear system performs random motion.


## 2.1. Phase Space Reconstruction

The reconstruction of the phase space is the basis for the quantitative analysis of chaotic time series, in which the embedding dimension and the delay coordinate are the two most critical parameters [30]. The time delay method is currently the most commonly used method. For univariate chaotic time series xi,i=1,2,…,n,(1)yi=xi,xi+τ,…,xi+m−1τ,i=1,2,…,n−m−1τ.According to Takens’ Theorem, the appropriate choice of the embedding dimensionm and the delay time τ can restore the dynamics properties of the original state space in the sense of topological equivalence.


### 2.1.1. Delay Time

The mutual information method [31] is introduced to determine the delay time τ of the measured displacement sequence of the dam, shown as follows:(2)Iτ=∑xi,xi+τPxi,xi+τlog2Pxi,xi+τPxixi+τ,where Pxi is the normalized distribution of xi, Pxi+τ is the normalized distribution of xi+τ, and Pxi,xi+τ is the joint distribution of xi and xi+τ. The time at which the first minimum point appears in the τ∼Iτ curve is often selected as the optimal delay time.


### 2.1.2. Embedding Dimension

The embedding dimensionm is determined by the Cao method [32]. The distance at,m between the phase point and the nearest neighbor point is shown as follows:(3)at,m=ytm+1−yfm+1ytm−yfm,where ytm+1 is the t th vector in the reconstructed phase space with the embedding dimension m+1 and yfm+1 is the nearest neighbor to ytm+1. yfm is the nearest neighbor to ytm in the reconstructed phase space with the embedding dimension m.The average valueEm of at,m is calculated as follows:(4)Em=1n−mτ∑t=1n−mτat,m.Then, the change ofEm is as follows:(5)E1m=Em+1Em.Whenm>m0, if E1m no longer changes, m0 represents the minimum embedding dimension. In order to avoid the situation that the change of E1m is difficult to judge, the Cao method adds another definition.(6)E′m=1n−mτ∑t=1n−mτytm+1−yfm+1,E2m=E′m+1E′m.For chaotic time series, there will always be some value ofm, so that E2m≠1. Whether E1m tends to be stable and whether the value of m can achieve E2m≠1 are observed to determine the minimum m of the reconstructed phase space.


## 2.1.1. Delay Time

The mutual information method [31] is introduced to determine the delay time τ of the measured displacement sequence of the dam, shown as follows:(2)Iτ=∑xi,xi+τPxi,xi+τlog2Pxi,xi+τPxixi+τ,where Pxi is the normalized distribution of xi, Pxi+τ is the normalized distribution of xi+τ, and Pxi,xi+τ is the joint distribution of xi and xi+τ. The time at which the first minimum point appears in the τ∼Iτ curve is often selected as the optimal delay time.


## 2.1.2. Embedding Dimension

The embedding dimensionm is determined by the Cao method [32]. The distance at,m between the phase point and the nearest neighbor point is shown as follows:(3)at,m=ytm+1−yfm+1ytm−yfm,where ytm+1 is the t th vector in the reconstructed phase space with the embedding dimension m+1 and yfm+1 is the nearest neighbor to ytm+1. yfm is the nearest neighbor to ytm in the reconstructed phase space with the embedding dimension m.The average valueEm of at,m is calculated as follows:(4)Em=1n−mτ∑t=1n−mτat,m.Then, the change ofEm is as follows:(5)E1m=Em+1Em.Whenm>m0, if E1m no longer changes, m0 represents the minimum embedding dimension. In order to avoid the situation that the change of E1m is difficult to judge, the Cao method adds another definition.(6)E′m=1n−mτ∑t=1n−mτytm+1−yfm+1,E2m=E′m+1E′m.For chaotic time series, there will always be some value ofm, so that E2m≠1. Whether E1m tends to be stable and whether the value of m can achieve E2m≠1 are observed to determine the minimum m of the reconstructed phase space.


## 2.2. Lyapunov Exponent Method

The maximum Lyapunov exponent (λmax) is usually regarded as an indicator of chaotic motion [33]. λmax>0 indicates that the system is in a chaotic state. The specific process of calculating λmax by the small data amount method is as follows.Firstly, determine the nearest neighbor pointyj^ of the phase point yj, and its short-term distance is dj0=yj−yj⌢.Secondly, calculate the distancedji=yj+i−yj⌢+i between yj and yj^ after the i th discrete time steps.Thirdly, get the average of alllndji, which is yi:(7)yi=1qΔt∑t=1qlndji,where q represents the number of i satisfying the condition dji≠0 and Δt represents the time interval.Finally, calculateλmax. λmax is the slope of the regression line obtained by the least-squares method.


## 2.3. Correlation Dimension Method

The correlation dimensionD2 is mainly determined by the Grassberger-Procaccia algorithm [34]. Suppose that r is the radius of the sphere centered on yi and yj; then the correlation integral Cr is given by(8)Cr=limn⟶∞2nn−1∑i,j=1nθr−yi−yj,where θ⋅ is the Heaviside function:(9)θu=0,u≤0,1,u>0,limr⟶0Cr∞rDr⟶0,where D2 is the correlation dimension:(10)D2=logCnrlogr.Thus, drawlogCnr/logr curve, and then the value of D2 can be determined according to the slope of the curve. As the embedding dimension m gradually increases, the slope of the curve converges, and the limit of convergence is the correlation dimension D2. The slope of the curve of a stochastic system will continue to increase with the increase of the embedding dimension, and there will be no convergence phenomenon.


## 2.4. Kolmogorov Entropy Method

Kolmogorov entropy [35] describes the generation rate of chaotic orbital information over time. Kolmogorov entropy reflects the chaos level of nonlinear dynamic systems, and the Kolmogorov entropy proposed by Grassberger and Procaccia is most commonly used as its estimate. The Kolmogorov entropy could be estimated by the correlation integral method.(11)K2=−limr⟶0limd⟶∞1Δmτlog2CmrCm+Δmr.When the embedding dimensionm is continuously increasing at intervals of Δm, the stable estimation of Kolmogorov entropy can be obtained through the equal slope regression of equation (11). It should be noted that the minimum value of m must be an integer greater than the correlation dimension D2.In a phase space with embedding dimensioni, there is(12)xij=log2rij,yij=log2Crij,where(13)yij=axij+b.Leta=D2, and, for the embedding dimensions i and i+Δm, there is(14)K2=limi⟶∞ΔbiΔmτ,where Δb=bi−bi+m.The Kolmogorov entropy estimateK2 can be used to judge the motion properties of the nonlinear system: K2=0 means the nonlinear system performs regular motion, K2>0 means the nonlinear system performs chaotic motion, and K2<0 means the nonlinear system performs random motion.


## 3. Chaotic Time Series Prediction

For the dam deformation time seriesxii=1,2,…,n, when the delay time τ and the embedding dimension m have been determined, the phase space reconstruction results of series xi are as follows:(15)X=x1x1+τx1+2τ⋯x1+m−1τx2x2+τx2+2τ⋯x2+m−1τx3x3+τx3+2τ⋯x3+m−1τ⋮⋮⋮⋮xn−m−1τxn−m−1τ+τxn−m−1τ+2τ⋯xn,Y=x2+m−1τx3+m−1τ⋮xn+1.After determining the input variables and output variables of the model, a novel model based on phase space reconstruction for dam deformation predicting is proposed.To clarify the influence of the reconstructed phase space as an input variable on the prediction performance of the model, a conventional model with dam deformation influence factors as input variables should be established. Early studies [3, 36, 37] have shown that water level, temperature, and aging are the main factors affecting dam deformation, shown as follows:(16)H−H0,H−H02,H−H03,H−H04,sin2πit365−sin2πit0365,cos2πit365−cos2πit0365,sin4πit365−sin4πit0365,cos4πit365−cos4πit0365,θ−θ0,lnθ−lnθ0.In order to match the consistency of the model and avoid the larger data information overwhelming the smaller data information, the input data of SVM is normalized. After completing the SVM training process, the output data of the SVM needs to be denormalized.Normalization equation is as follows:(17)xi′=xi−xminxmax−xmin.Antinormalization equation is as follows:(18)xi=xi′xmax−xmin+xmin,where xi is a sample data; xmin and xmax, respectively, represent the minimum and maximum sample data; and xi′ is the normalized data.The mean square error (MSE), mean absolute percentage error (MAPE), and square correlation coefficient (R2) are used to evaluate the performance of predictive models, as shown as follows [11]:(19)MSE=1N∑i=1NyDi−yi2,MAPE=1N∑i=1NyDi−yiyi,R2=∑i=1NyDi−y¯Dyi−y¯2∑i=1NyDi−y¯D2∑i=1Nyi−y¯2,where yD and y¯D represent predicted values and predicted average values, y and y¯ represent measured values and measured average values, and N represents the number of observed samples. The closer R2 is to 1, the smaller the MSE and the MAPE are and the better the prediction effect of the model is.


### 3.1. Support Vector Machine

SVM [13] usually needs to establish a suitable function fx to describe the nonlinear relationship between the characteristic value xi and the target value yi, shown as follows:(20)fxi=w⋅φxi+b,where w is the coefficient vector, φxi is the transformation function, and w and b represent the weight and bias, respectively. w and b are estimated by minimizing the regularized hazard function, shown as follows:(21)Rw=12w2+C∑i=1nLεyi,fxi,where 1/2w2 is the regularization term, C is the penalty coefficient, and Lεyi,fxi is the ε -insensitive loss function, calculated by the following equation:(22)Lεyi,fxi=max0,yi−fxi−ε,where ε represents the allowable error threshold. The optimization object can be deducted as follows:(23)minfw,ξ−,ξ+=12w2+C∑i=1nξ−,ξ+.It is subject to(24)yi−w⋅φxi−b≤ε+ξ−,ξ−≥0,w⋅φxi+b−yi≤ε+ξ+,ξ+≥0,where ξ+ and ξ− represent slack variables.The key is to establish the Lagrangian function.(25)maxH∂i−,∂i+=−12∑i=1n∑j=1n∂i−−∂i+∂j−−∂j+Kxi,xj+∑i=1nyi∂i−−∂i+−ε∑i=1nyi∂i−+∂i+.It is subject to(26)∑i=1n∂i−−∂i+=0,∂i−,∂i+∈0,C.Therefore,(27)fx=∑i=1n∂i−−∂i+Kxi,xj+b,where Kxi,xj represents kernel function, including polynomial, radial basis function, and sigmoid.


### 3.2. Improved Grey Wolf Optimizer (IGWO)

#### 3.2.1. Grey Wolf Optimizer

Grey wolves [17, 38] have a very strict social dominance hierarchy, which is mainly divided into four parts: α, β, δ, and ω. α is the best solution, followed by β and δ, and the remaining solutions belong to ω. The top three best wolves that are closest to their prey are α, β, and δ, and they guide ω to search for prey in promising search areas [17]. During the hunting process, the wolf will update its position around α, β, and δ, shown as follows [17]:(28)D⟶=C⟶⋅X⟶pt−X⟶t,X⟶t+1=X⟶pt−A⟶⋅D⟶,where t is the current iteration number, X⟶pt and X⟶t represent the current positions of the prey and the wolf, respectively; and D⟶ is the distance between the wolf and the prey.The coefficient vectorsA⟶ and C⟶ can be expressed as follows:(29)A⟶=2a⟶⋅r⟶1−a⟶,(30)C⟶=2r⟶2,where r⟶1 and r⟶2 are two vectors randomly generated between [0, 1] and the convergence factor a⟶ linearly decreases from 2 to 0 in the iterative process.Update the position of grey wolves.(31)D⟶α=C⟶1⋅X⟶α−X⟶,D⟶β=C⟶2⋅X⟶β−X⟶,D⟶δ=C⟶3⋅X⟶δ−X⟶,(32)X⟶1=X⟶α−A1⋅D⟶α,X⟶2=X⟶β−A2⋅D⟶β,X⟶3=X⟶δ−A3⋅D⟶δ,(33)X⟶t+1=X⟶1+X⟶2+X⟶33,where X⟶α, X⟶β, and X⟶δ represent the positions of α, β, and δ, respectively, X⟶ represents the location of the current solution, C⟶1, C⟶2, and C⟶3 represent randomly generated vectors, and D⟶α, D⟶β, and D⟶δ represent the distances of α, β, and δ from other grey wolves, respectively. A⟶1 , A⟶2, and A⟶3 are random vectors, and X⟶t+1 is the final position of ω.


#### 3.2.2. Differential Evolution

Differential evolution (DE) algorithm is an efficient group-based heuristic search algorithm [39]. It mainly obtains the optimal solution through three operations of mutation, crossover, and selection. The population size is N, D represents the dimensionality in the solution space, xi=x1,i,x2,i,…,xD,L is the individual vector, and G=0,1,…,Gmax is the iteration time. xiG=x1,iG,x2,iG,…,xD,iG is the i-th individual in the G-th generation. xL=x1,L,x2,L,…,xD,L and xU=x1,U,x2,U,…,xD,U record the lower limit and upper limit, respectively, of search space.The initial populationP0 is mainly randomly generated within the upper limit and lower limit xL,xU, and the j-th index of the i-th individual is obtained by the following equation:(34)xj,i0=xj,L+rand0,1×x1,U−xj,L.The mutation operation generates a new mutation vectorviG, shown as follows:(35)viG=xr1G+Fxr2G−xr3G,where r1,r2,r3∈1,2,…,N are randomly generated integers and F is the magnification ratio of the control difference vector, which is a real number with a varying range between [0, 2].The crossover operation is as follows:(36)uiG=vj,iG,ifrandj,i≤CR orj=jrand,xj,iG,otherwise,where CR is the cross probability, which takes a value between [0, 1].Finally, the selection operation is performed. For specific problems, all mutation crossover individuals are evaluated. For specific problems, all mutation crossover individuals are evaluated. If the fitness of the current individual exceeds the previous generation, it means that the mutation crossover operation is successful, and the current individual is retained; if the current individual's fitness is not as good as the previous generation, the better individual is retained. The individual with the optimal fitness will become the optimal value of this generation of individuals. When the termination condition is met, the evolution will stop; otherwise, the next round will continue.


#### 3.2.3. Improved Grey Wolf Optimization (IGWO) Algorithm

The conventional GWO algorithm randomly generates the initial population, which may fall into the dilemma of local optimization [17, 40]. The DE algorithm generates group intelligence through mutual cooperation and competition between individuals. Based on the respective advantages and disadvantages of the GWO and DE algorithms, a more efficient improved algorithm (DE-GWO, IGWO) is proposed. Specific steps are as follows:(1)
Set the relevant parameters of the IGWO algorithm, such as the population sizeN, the maximum number of iterations tmax, and the upper and lower limits of the search range, ub and lb.(2)
Initialize parametersa, A, and C. Generate intermediates (variant populations) through evolutionary mutation operations, and then generate initial population individuals through competitive selection operations, and set the iteration time t=1.(3)
Calculate the objective function value of a single grey wolf individual and determine the best three individuals asXα, Xβ, and Xδ, respectively.(4)
Calculate the distance between other grey wolf individuals and optimalXα, Xβ, and Xδ according to equation (31), and update the position of each grey wolf according to equation (32).(5)
Update the values ofa, A, and C. Crossover and competition operations are applied to retain better individual positions and generate new individuals, respectively.(6)
Update the position of the first three grey wolvesXα, Xβ, and Xδ.(7)
Determine whether the maximum number of iterationstmax has been reached. If yes, exit and output the current objective function value of Xα; otherwise, t=t+1, and move to the third step to continue.


### 3.3. Construction Process of PSR-SVM-IGWO Model Predicting Dam Deformation

A hybrid model combining chaos theory and SVM is proposed, and the improved GWO algorithm is used to select optimal parameters for concrete dam deformation analysis and prediction. Use the method mentioned in Section2 to identify the chaotic characteristics of the deformed time series, and then determine the input variables of the support vector machine. According to the IGWO algorithm proposed in Section 3.2, the parameter optimization of SVM is realized. By performing training operations, a PSR-SVM-IGWO prediction model can be established. The specific process is shown in Figure 1.(1)
Obtain dam deformation data.(2)
Reconstruct the phase space of the observed displacement time series. The mutual information method and the Cao method are adopted to determine the optimal delay timeτ and the minimum embedding dimension m, and the phase space is reconstructed accordingly.(3)
Identification of chaotic characteristics of the observed displacement time series. Estimate the maximum Lyapunov exponentλmax, the correlation dimension D2, and the Kolmogorov entropy K.λmax>0. The Kolmogorov entropy estimate K2 is a finite positive value and the saturation of the correlation dimension indicates that the observed displacement time series has the chaotic characteristics.(4)
Determine the input and output variables according to equation (15).(5)
Divide the data set into training data set and test data set, and normalize the data set.(6)
Initialize the parameters of the grey wolf population. Set the parameters such as the population size and the maximum number of iterations to determine the scope of the search space, and at the same time encode the penalty coefficientC and the kernel parameter γ.(7)
Use DE algorithm to initialize grey wolf population.(8)
Initializea, A, and C. Use equations (29) and (30) to initialize coefficients a and A, and the convergence factor a = 2.(9)
Calculate the fitness value of the grey wolf individual. The combination parameters (C and γ) are determined by the grey wolf's position, and, based on this combination of parameters, the SVM model is trained with training data, and the calculated MSE is used to determine the fitness value of the grey wolf individual.(10)
Keep the three grey wolves with the best fitness value, namely,α wolf, β wolf, and δ wolf.(11)
Refer to the positions ofα wolf, β wolf, and δ wolf, and use equations (31)–(34) to update the positions of individual grey wolves.(12)
Calculate the fitness values of all grey wolves. Determine the combination parameters (C and γ) based on the grey wolf's position, and use the training data to train the SVM model based on C and γ, and determine the individual fitness value of the grey wolf through the calculated MSE.(13)
Update the positions and fitness values ofα wolf, β wolf, and δ wolf.(14)
Updatea, A, and C. Use equations (29) and (30) to initialize coefficients A and C, and the convergence factor a here starts from a = 2 and decreases linearly until a = 0.(15)
Determine whether the maximum number of iterations has been reached. If it is reached, output the best parameter combination (C, γ); if it is not reached, return to step 11, and increase the number of iterations by 1.(16)
Output the optimal parameter combination (C, γ), and retrain the SVM based on the training set.(17)
Use the IGWO-SVM model to predict the sample to be tested, and evaluate the performance of the model's prediction results.Figure 1
Construction process of dam deformation prediction model based on PSR-SVM-IGWO.


## 3.1. Support Vector Machine

SVM [13] usually needs to establish a suitable function fx to describe the nonlinear relationship between the characteristic value xi and the target value yi, shown as follows:(20)fxi=w⋅φxi+b,where w is the coefficient vector, φxi is the transformation function, and w and b represent the weight and bias, respectively. w and b are estimated by minimizing the regularized hazard function, shown as follows:(21)Rw=12w2+C∑i=1nLεyi,fxi,where 1/2w2 is the regularization term, C is the penalty coefficient, and Lεyi,fxi is the ε -insensitive loss function, calculated by the following equation:(22)Lεyi,fxi=max0,yi−fxi−ε,where ε represents the allowable error threshold. The optimization object can be deducted as follows:(23)minfw,ξ−,ξ+=12w2+C∑i=1nξ−,ξ+.It is subject to(24)yi−w⋅φxi−b≤ε+ξ−,ξ−≥0,w⋅φxi+b−yi≤ε+ξ+,ξ+≥0,where ξ+ and ξ− represent slack variables.The key is to establish the Lagrangian function.(25)maxH∂i−,∂i+=−12∑i=1n∑j=1n∂i−−∂i+∂j−−∂j+Kxi,xj+∑i=1nyi∂i−−∂i+−ε∑i=1nyi∂i−+∂i+.It is subject to(26)∑i=1n∂i−−∂i+=0,∂i−,∂i+∈0,C.Therefore,(27)fx=∑i=1n∂i−−∂i+Kxi,xj+b,where Kxi,xj represents kernel function, including polynomial, radial basis function, and sigmoid.


## 3.2. Improved Grey Wolf Optimizer (IGWO)

### 3.2.1. Grey Wolf Optimizer

Grey wolves [17, 38] have a very strict social dominance hierarchy, which is mainly divided into four parts: α, β, δ, and ω. α is the best solution, followed by β and δ, and the remaining solutions belong to ω. The top three best wolves that are closest to their prey are α, β, and δ, and they guide ω to search for prey in promising search areas [17]. During the hunting process, the wolf will update its position around α, β, and δ, shown as follows [17]:(28)D⟶=C⟶⋅X⟶pt−X⟶t,X⟶t+1=X⟶pt−A⟶⋅D⟶,where t is the current iteration number, X⟶pt and X⟶t represent the current positions of the prey and the wolf, respectively; and D⟶ is the distance between the wolf and the prey.The coefficient vectorsA⟶ and C⟶ can be expressed as follows:(29)A⟶=2a⟶⋅r⟶1−a⟶,(30)C⟶=2r⟶2,where r⟶1 and r⟶2 are two vectors randomly generated between [0, 1] and the convergence factor a⟶ linearly decreases from 2 to 0 in the iterative process.Update the position of grey wolves.(31)D⟶α=C⟶1⋅X⟶α−X⟶,D⟶β=C⟶2⋅X⟶β−X⟶,D⟶δ=C⟶3⋅X⟶δ−X⟶,(32)X⟶1=X⟶α−A1⋅D⟶α,X⟶2=X⟶β−A2⋅D⟶β,X⟶3=X⟶δ−A3⋅D⟶δ,(33)X⟶t+1=X⟶1+X⟶2+X⟶33,where X⟶α, X⟶β, and X⟶δ represent the positions of α, β, and δ, respectively, X⟶ represents the location of the current solution, C⟶1, C⟶2, and C⟶3 represent randomly generated vectors, and D⟶α, D⟶β, and D⟶δ represent the distances of α, β, and δ from other grey wolves, respectively. A⟶1 , A⟶2, and A⟶3 are random vectors, and X⟶t+1 is the final position of ω.


### 3.2.2. Differential Evolution

Differential evolution (DE) algorithm is an efficient group-based heuristic search algorithm [39]. It mainly obtains the optimal solution through three operations of mutation, crossover, and selection. The population size is N, D represents the dimensionality in the solution space, xi=x1,i,x2,i,…,xD,L is the individual vector, and G=0,1,…,Gmax is the iteration time. xiG=x1,iG,x2,iG,…,xD,iG is the i-th individual in the G-th generation. xL=x1,L,x2,L,…,xD,L and xU=x1,U,x2,U,…,xD,U record the lower limit and upper limit, respectively, of search space.The initial populationP0 is mainly randomly generated within the upper limit and lower limit xL,xU, and the j-th index of the i-th individual is obtained by the following equation:(34)xj,i0=xj,L+rand0,1×x1,U−xj,L.The mutation operation generates a new mutation vectorviG, shown as follows:(35)viG=xr1G+Fxr2G−xr3G,where r1,r2,r3∈1,2,…,N are randomly generated integers and F is the magnification ratio of the control difference vector, which is a real number with a varying range between [0, 2].The crossover operation is as follows:(36)uiG=vj,iG,ifrandj,i≤CR orj=jrand,xj,iG,otherwise,where CR is the cross probability, which takes a value between [0, 1].Finally, the selection operation is performed. For specific problems, all mutation crossover individuals are evaluated. For specific problems, all mutation crossover individuals are evaluated. If the fitness of the current individual exceeds the previous generation, it means that the mutation crossover operation is successful, and the current individual is retained; if the current individual's fitness is not as good as the previous generation, the better individual is retained. The individual with the optimal fitness will become the optimal value of this generation of individuals. When the termination condition is met, the evolution will stop; otherwise, the next round will continue.


### 3.2.3. Improved Grey Wolf Optimization (IGWO) Algorithm

The conventional GWO algorithm randomly generates the initial population, which may fall into the dilemma of local optimization [17, 40]. The DE algorithm generates group intelligence through mutual cooperation and competition between individuals. Based on the respective advantages and disadvantages of the GWO and DE algorithms, a more efficient improved algorithm (DE-GWO, IGWO) is proposed. Specific steps are as follows:(1)
Set the relevant parameters of the IGWO algorithm, such as the population sizeN, the maximum number of iterations tmax, and the upper and lower limits of the search range, ub and lb.(2)
Initialize parametersa, A, and C. Generate intermediates (variant populations) through evolutionary mutation operations, and then generate initial population individuals through competitive selection operations, and set the iteration time t=1.(3)
Calculate the objective function value of a single grey wolf individual and determine the best three individuals asXα, Xβ, and Xδ, respectively.(4)
Calculate the distance between other grey wolf individuals and optimalXα, Xβ, and Xδ according to equation (31), and update the position of each grey wolf according to equation (32).(5)
Update the values ofa, A, and C. Crossover and competition operations are applied to retain better individual positions and generate new individuals, respectively.(6)
Update the position of the first three grey wolvesXα, Xβ, and Xδ.(7)
Determine whether the maximum number of iterationstmax has been reached. If yes, exit and output the current objective function value of Xα; otherwise, t=t+1, and move to the third step to continue.


## 3.2.1. Grey Wolf Optimizer

Grey wolves [17, 38] have a very strict social dominance hierarchy, which is mainly divided into four parts: α, β, δ, and ω. α is the best solution, followed by β and δ, and the remaining solutions belong to ω. The top three best wolves that are closest to their prey are α, β, and δ, and they guide ω to search for prey in promising search areas [17]. During the hunting process, the wolf will update its position around α, β, and δ, shown as follows [17]:(28)D⟶=C⟶⋅X⟶pt−X⟶t,X⟶t+1=X⟶pt−A⟶⋅D⟶,where t is the current iteration number, X⟶pt and X⟶t represent the current positions of the prey and the wolf, respectively; and D⟶ is the distance between the wolf and the prey.The coefficient vectorsA⟶ and C⟶ can be expressed as follows:(29)A⟶=2a⟶⋅r⟶1−a⟶,(30)C⟶=2r⟶2,where r⟶1 and r⟶2 are two vectors randomly generated between [0, 1] and the convergence factor a⟶ linearly decreases from 2 to 0 in the iterative process.Update the position of grey wolves.(31)D⟶α=C⟶1⋅X⟶α−X⟶,D⟶β=C⟶2⋅X⟶β−X⟶,D⟶δ=C⟶3⋅X⟶δ−X⟶,(32)X⟶1=X⟶α−A1⋅D⟶α,X⟶2=X⟶β−A2⋅D⟶β,X⟶3=X⟶δ−A3⋅D⟶δ,(33)X⟶t+1=X⟶1+X⟶2+X⟶33,where X⟶α, X⟶β, and X⟶δ represent the positions of α, β, and δ, respectively, X⟶ represents the location of the current solution, C⟶1, C⟶2, and C⟶3 represent randomly generated vectors, and D⟶α, D⟶β, and D⟶δ represent the distances of α, β, and δ from other grey wolves, respectively. A⟶1 , A⟶2, and A⟶3 are random vectors, and X⟶t+1 is the final position of ω.


## 3.2.2. Differential Evolution

Differential evolution (DE) algorithm is an efficient group-based heuristic search algorithm [39]. It mainly obtains the optimal solution through three operations of mutation, crossover, and selection. The population size is N, D represents the dimensionality in the solution space, xi=x1,i,x2,i,…,xD,L is the individual vector, and G=0,1,…,Gmax is the iteration time. xiG=x1,iG,x2,iG,…,xD,iG is the i-th individual in the G-th generation. xL=x1,L,x2,L,…,xD,L and xU=x1,U,x2,U,…,xD,U record the lower limit and upper limit, respectively, of search space.The initial populationP0 is mainly randomly generated within the upper limit and lower limit xL,xU, and the j-th index of the i-th individual is obtained by the following equation:(34)xj,i0=xj,L+rand0,1×x1,U−xj,L.The mutation operation generates a new mutation vectorviG, shown as follows:(35)viG=xr1G+Fxr2G−xr3G,where r1,r2,r3∈1,2,…,N are randomly generated integers and F is the magnification ratio of the control difference vector, which is a real number with a varying range between [0, 2].The crossover operation is as follows:(36)uiG=vj,iG,ifrandj,i≤CR orj=jrand,xj,iG,otherwise,where CR is the cross probability, which takes a value between [0, 1].Finally, the selection operation is performed. For specific problems, all mutation crossover individuals are evaluated. For specific problems, all mutation crossover individuals are evaluated. If the fitness of the current individual exceeds the previous generation, it means that the mutation crossover operation is successful, and the current individual is retained; if the current individual's fitness is not as good as the previous generation, the better individual is retained. The individual with the optimal fitness will become the optimal value of this generation of individuals. When the termination condition is met, the evolution will stop; otherwise, the next round will continue.


## 3.2.3. Improved Grey Wolf Optimization (IGWO) Algorithm

The conventional GWO algorithm randomly generates the initial population, which may fall into the dilemma of local optimization [17, 40]. The DE algorithm generates group intelligence through mutual cooperation and competition between individuals. Based on the respective advantages and disadvantages of the GWO and DE algorithms, a more efficient improved algorithm (DE-GWO, IGWO) is proposed. Specific steps are as follows:(1)
Set the relevant parameters of the IGWO algorithm, such as the population sizeN, the maximum number of iterations tmax, and the upper and lower limits of the search range, ub and lb.(2)
Initialize parametersa, A, and C. Generate intermediates (variant populations) through evolutionary mutation operations, and then generate initial population individuals through competitive selection operations, and set the iteration time t=1.(3)
Calculate the objective function value of a single grey wolf individual and determine the best three individuals asXα, Xβ, and Xδ, respectively.(4)
Calculate the distance between other grey wolf individuals and optimalXα, Xβ, and Xδ according to equation (31), and update the position of each grey wolf according to equation (32).(5)
Update the values ofa, A, and C. Crossover and competition operations are applied to retain better individual positions and generate new individuals, respectively.(6)
Update the position of the first three grey wolvesXα, Xβ, and Xδ.(7)
Determine whether the maximum number of iterationstmax has been reached. If yes, exit and output the current objective function value of Xα; otherwise, t=t+1, and move to the third step to continue.


## 3.3. Construction Process of PSR-SVM-IGWO Model Predicting Dam Deformation

A hybrid model combining chaos theory and SVM is proposed, and the improved GWO algorithm is used to select optimal parameters for concrete dam deformation analysis and prediction. Use the method mentioned in Section2 to identify the chaotic characteristics of the deformed time series, and then determine the input variables of the support vector machine. According to the IGWO algorithm proposed in Section 3.2, the parameter optimization of SVM is realized. By performing training operations, a PSR-SVM-IGWO prediction model can be established. The specific process is shown in Figure 1.(1)
Obtain dam deformation data.(2)
Reconstruct the phase space of the observed displacement time series. The mutual information method and the Cao method are adopted to determine the optimal delay timeτ and the minimum embedding dimension m, and the phase space is reconstructed accordingly.(3)
Identification of chaotic characteristics of the observed displacement time series. Estimate the maximum Lyapunov exponentλmax, the correlation dimension D2, and the Kolmogorov entropy K.λmax>0. The Kolmogorov entropy estimate K2 is a finite positive value and the saturation of the correlation dimension indicates that the observed displacement time series has the chaotic characteristics.(4)
Determine the input and output variables according to equation (15).(5)
Divide the data set into training data set and test data set, and normalize the data set.(6)
Initialize the parameters of the grey wolf population. Set the parameters such as the population size and the maximum number of iterations to determine the scope of the search space, and at the same time encode the penalty coefficientC and the kernel parameter γ.(7)
Use DE algorithm to initialize grey wolf population.(8)
Initializea, A, and C. Use equations (29) and (30) to initialize coefficients a and A, and the convergence factor a = 2.(9)
Calculate the fitness value of the grey wolf individual. The combination parameters (C and γ) are determined by the grey wolf's position, and, based on this combination of parameters, the SVM model is trained with training data, and the calculated MSE is used to determine the fitness value of the grey wolf individual.(10)
Keep the three grey wolves with the best fitness value, namely,α wolf, β wolf, and δ wolf.(11)
Refer to the positions ofα wolf, β wolf, and δ wolf, and use equations (31)–(34) to update the positions of individual grey wolves.(12)
Calculate the fitness values of all grey wolves. Determine the combination parameters (C and γ) based on the grey wolf's position, and use the training data to train the SVM model based on C and γ, and determine the individual fitness value of the grey wolf through the calculated MSE.(13)
Update the positions and fitness values ofα wolf, β wolf, and δ wolf.(14)
Updatea, A, and C. Use equations (29) and (30) to initialize coefficients A and C, and the convergence factor a here starts from a = 2 and decreases linearly until a = 0.(15)
Determine whether the maximum number of iterations has been reached. If it is reached, output the best parameter combination (C, γ); if it is not reached, return to step 11, and increase the number of iterations by 1.(16)
Output the optimal parameter combination (C, γ), and retrain the SVM based on the training set.(17)
Use the IGWO-SVM model to predict the sample to be tested, and evaluate the performance of the model's prediction results.Figure 1
Construction process of dam deformation prediction model based on PSR-SVM-IGWO.


## 4. Case Study

### 4.1. Engineering Overview

The Xiluodu dam is a concrete arch dam on the Jinsha River in China, and its geographic location is shown in Figure2. The dam crest elevation is 610.0 m, and the maximum dam height is 285.5 m. The normal water storage level is 600.0 m, and the dead water level is 540.0 m. The layout of the dam vertical line monitoring system is shown in Figure 3.Figure 2
Geographical location of Xiluodu arch dam.Figure 3
Vertical arrangement of the dam body.The monitoring point PL15-1 is taken as an example. From June 2014 to June 2018, the radial displacement and water level changes of the measuring point are shown in Figure4. Take the observation data from June 2014 to December 2015 as the training set, and take the rest as the prediction set. As shown in Figure 4, the symbol “-” indicates downstream displacement, and the symbol “+” indicates upstream displacement. During the study period, the water level change range is [540 m, 610 m], and the displacement showed obvious periodic changes with the water level change, and the change range is [-60 mm, 20 mm]. The displacement of the measuring point presents an obvious periodic change with time, but the peak displacement gradually decreases with the change of time. The main reason is that the Xiluodu arch dam was in the primary impoundment stage during the study period, and the displacement was affected by creep.Figure 4
Time series of the radial displacement and water level of the measuring point.


### 4.2. Prediction Model Construction

The delay timeτ of the observed displacement time series is estimated by the mutual information method, as shown in Figure 5. According to the principle of mutual information, the optimal delay time τ=2.Figure 5
Iτ∼τ curve for the observed displacement data series of PL15-1.According to the determinedτ, the Cao method is adopted to calculate m. The change curves of E1m and E2m as m increases are shown in Figure 6. When m>4, the change of E1m starts to become smaller, and thus the optimal embedding dimension is determined; m=4. It can also be found that there are some values of m such that E2m≠1, which can prove that the observed displacement data series comes from a deterministic process.Figure 6
E1m∼m and E2m∼m curves for the observed displacement data series of PL15-1.According to the determinedτ, the G-P method is adopted to calculate the correlation dimension of the observed displacement data series. When the value of the embedding dimension d ranges from 1 to 10, the corresponding double logarithmic curve of lnCr∼lnr is shown in Figure 7. It is shown in Figure 7 that the correlation index is saturated to a constant. According to the definition, it is concluded that the observed displacement data series is a chaotic series and the correlation dimension D2 = 0.9043.Figure 7
Double logarithmic curve of correlation dimension.The small data set method is introduced to calculate the maximum Lyapunov exponentλmax of the observed displacement data series, as shown in Figure 8. λmax can be calculated from the fitting straight line of the corresponding curve: λmax=0.001511. λmax>0 indicates that the observed displacement data series has chaotic characteristics.Figure 8
The calculation diagram of the largest Lyapunov exponent of the observed displacement data series of PL15-1.The correlation integral method is used to calculate the Kolmogorov entropy estimateK2 of the observed displacement data series. According to the double logarithmic curve of lnCr∼lnr, as shown in Figure 7, the double logarithmic curve of log2r∼log2Cr can be derived, so that the Kolmogorov entropy estimate K2 of the observed displacement data series can be obtained when the delay time is determined, and K2 = 0.0306. K2 takes a finite positive value, which indicates the chaotic characteristics of the observed displacement data series at the measuring point PL15-1.Based on the abovementioned results, it can be concluded that the observed displacement data series of the measuring point PL15-1 has chaotic characteristics, and a chaotic prediction model of dam displacement can be established.The proposed IGWO algorithm is adopted to realize the parameter optimization of SVM. The phase space of reconstructed observation displacement and the influencing factors of dam deformation are used as input variables, respectively, to evaluate the predicting performance of the IGWO and GWO algorithms. The population sizeN=20, the maximum number of iterations tmax=200, ε=0.01, the penalty factor C∈0.01,100, and the kernel parameter γ∈0.01,1000. The calculation is terminated when the number of iterations reaches 200.The fitness curves of the GWO and IGWO algorithms are shown in Figure9. For the conventional model with dam deformation influencing factors as input variables, the optimal parameters [C, γ] of SVM obtained by the IGWO algorithm are [0.1088, 24.2515]. For the chaotic model with the reconstructed observed deformation phase space as the input variable, the optimal parameters [C, γ] of the SVM obtained by the IGWO algorithm are [3.5801, 29.4460].Figure 9
Relation curves between fitness value and iteration number: conventional model (a) and chaotic model (b).It can be seen from Figure9 that, whether it is a conventional model with dam deformation influencing factors as input variables or a chaotic model with reconstructed observation displacement data series phase space as the input variable, the IGWO algorithm reduces the number of iterations and can find the solution closest to the best goal faster. The main explanation is that the IGWO algorithm enriches the diversity of the initial population and improves the global search capability, thereby accelerating the convergence speed and convergence accuracy.


### 4.3. Results

For the PSR-SVM-IGWO-based dam observation displacement prediction model, the relevant information is introduced as follows. The SVM is at the heart of this innovative combination model. The input variable is the reconstructed phase space of the measured displacement data sequence, and the IGWO algorithm is used to realize the parameter optimization of SVM. In addition, the kernel function of SVM is a radial basis function.In order to better analyze the predictive performance of the PSR-SVM-IGWO model, the PSR-SVM and PSR-SVM-GWO models that take the reconstructed phase space of the observed displacement data sequence as input variables are established, respectively. In addition, the SVM, SVM-GWO, and SVM-IGWO models with the factors affecting dam deformation as input variables are established to explore the influence of different input variables on the accuracy of model prediction.The prediction performances of the SVM, SVM-GWO, SVM-IGWO, PSR-SVM, PSR-SVM-GWO, and PSR-SVM-IGWO models are shown in Table1 and Figure 10. From Table 1, we can see that when the influencing factors of dam deformation are used as the input variables of the model, the square correlation coefficient (R2) is ranked from large to small as SVM-IGWO model > SVM-GWO model > SVM model, mean absolute percentage error (MAPE) is ranked from small to large as SVM-IGWO model < SVM model < SVM-GWO model, and mean square error (MSE) is ranked from small to large as SVM-IGWO model < SVM-GWO model < SVM model. The SVM-IGWO model has the largest R2 of 0.9747, and the smallest MAPE and MSE are 0.0569 and 6.1547, respectively. When the reconstructed phase space of the dam deformation time series is used as the input variable of the model, R2 is sorted from large to small as PSR-SVM-IGWO model > PSR-SVM-GWO model > PSR-SVM model. The ranking of MAPE and MSE from small to large is PSR-SVM-IGWO model < PSR-SVM-GWO model < PSR-SVM model. The PSR-SVM-IGWO model has the largest R2 of 0.9902, and the smallest MAPE and MSE are 0.0423 and 2.6086, respectively.Table 1
Predictive performances of SVM, SVM-GWO, and SVM-IGWO models.
Prediction modelMSEMAPER2SVM33.12950.11480.9517SVM-GWO21.05680.13950.9616SVM-IGWO6.15470.05690.9747PSR-SVM6.66190.06180.9802PSR-SVM-GWO2.95860.04440.9842PSR-SVM-IGWO2.60860.04230.9902Figure 10
The predicted results of all models.As shown in Figure10, the measured displacement shows an obvious periodic change, but, due to the influence of creep, the peak displacement presents a gradual decrease trend. SVM, SVM-GWO, SVM-IGWO, PSR-SVM, PSR-SVM-GWO, and PSR-SVM-IGWO models can effectively predict the trend of dam displacement; and the PSR-IGWO-SVR model has the highest prediction accuracy and the smallest fluctuation range of the prediction error. For the conventional model, as time goes by, the deviation between the predicted value of the model and the observed value becomes larger and larger, while for the PSR model, the deviation is greatly improved. The calculation results show the following: ① The GWO optimized SVM algorithm is applicable in dam deformation prediction. ② The IGWO algorithm proposed in this paper has more outstanding optimization ability in optimizing SVM parameters than the conventional GWO. ③ For ultrahigh arch dams, the PSR model with the reconstructed phase space as the input variable has higher prediction accuracy and smaller prediction error than its corresponding conventional model with dam deformation influence factors as input variables.In general, SVM, SVM-GWO, and SVM-IGWO models can predict dam deformation trends well. The prediction accuracy of the PSR model with the phase space of the reconstructed observation displacement data sequence as the input variable is higher, which provides a strong support for the discussion of the chaotic phenomenon in the dam deformation time series. The calculation results also prove the applicability of the GWO algorithm in the field of dam deformation prediction. IGWO has more outstanding optimization capabilities than GWO. However, the algorithm proposed in this paper has only obtained good application results on dam deformation prediction. The application prospects of dam seepage and stress prediction are temporarily unknown, which will be tried in future research work.


## 4.1. Engineering Overview

The Xiluodu dam is a concrete arch dam on the Jinsha River in China, and its geographic location is shown in Figure2. The dam crest elevation is 610.0 m, and the maximum dam height is 285.5 m. The normal water storage level is 600.0 m, and the dead water level is 540.0 m. The layout of the dam vertical line monitoring system is shown in Figure 3.Figure 2
Geographical location of Xiluodu arch dam.Figure 3
Vertical arrangement of the dam body.The monitoring point PL15-1 is taken as an example. From June 2014 to June 2018, the radial displacement and water level changes of the measuring point are shown in Figure4. Take the observation data from June 2014 to December 2015 as the training set, and take the rest as the prediction set. As shown in Figure 4, the symbol “-” indicates downstream displacement, and the symbol “+” indicates upstream displacement. During the study period, the water level change range is [540 m, 610 m], and the displacement showed obvious periodic changes with the water level change, and the change range is [-60 mm, 20 mm]. The displacement of the measuring point presents an obvious periodic change with time, but the peak displacement gradually decreases with the change of time. The main reason is that the Xiluodu arch dam was in the primary impoundment stage during the study period, and the displacement was affected by creep.Figure 4
Time series of the radial displacement and water level of the measuring point.


## 4.2. Prediction Model Construction

The delay timeτ of the observed displacement time series is estimated by the mutual information method, as shown in Figure 5. According to the principle of mutual information, the optimal delay time τ=2.Figure 5
Iτ∼τ curve for the observed displacement data series of PL15-1.According to the determinedτ, the Cao method is adopted to calculate m. The change curves of E1m and E2m as m increases are shown in Figure 6. When m>4, the change of E1m starts to become smaller, and thus the optimal embedding dimension is determined; m=4. It can also be found that there are some values of m such that E2m≠1, which can prove that the observed displacement data series comes from a deterministic process.Figure 6
E1m∼m and E2m∼m curves for the observed displacement data series of PL15-1.According to the determinedτ, the G-P method is adopted to calculate the correlation dimension of the observed displacement data series. When the value of the embedding dimension d ranges from 1 to 10, the corresponding double logarithmic curve of lnCr∼lnr is shown in Figure 7. It is shown in Figure 7 that the correlation index is saturated to a constant. According to the definition, it is concluded that the observed displacement data series is a chaotic series and the correlation dimension D2 = 0.9043.Figure 7
Double logarithmic curve of correlation dimension.The small data set method is introduced to calculate the maximum Lyapunov exponentλmax of the observed displacement data series, as shown in Figure 8. λmax can be calculated from the fitting straight line of the corresponding curve: λmax=0.001511. λmax>0 indicates that the observed displacement data series has chaotic characteristics.Figure 8
The calculation diagram of the largest Lyapunov exponent of the observed displacement data series of PL15-1.The correlation integral method is used to calculate the Kolmogorov entropy estimateK2 of the observed displacement data series. According to the double logarithmic curve of lnCr∼lnr, as shown in Figure 7, the double logarithmic curve of log2r∼log2Cr can be derived, so that the Kolmogorov entropy estimate K2 of the observed displacement data series can be obtained when the delay time is determined, and K2 = 0.0306. K2 takes a finite positive value, which indicates the chaotic characteristics of the observed displacement data series at the measuring point PL15-1.Based on the abovementioned results, it can be concluded that the observed displacement data series of the measuring point PL15-1 has chaotic characteristics, and a chaotic prediction model of dam displacement can be established.The proposed IGWO algorithm is adopted to realize the parameter optimization of SVM. The phase space of reconstructed observation displacement and the influencing factors of dam deformation are used as input variables, respectively, to evaluate the predicting performance of the IGWO and GWO algorithms. The population sizeN=20, the maximum number of iterations tmax=200, ε=0.01, the penalty factor C∈0.01,100, and the kernel parameter γ∈0.01,1000. The calculation is terminated when the number of iterations reaches 200.The fitness curves of the GWO and IGWO algorithms are shown in Figure9. For the conventional model with dam deformation influencing factors as input variables, the optimal parameters [C, γ] of SVM obtained by the IGWO algorithm are [0.1088, 24.2515]. For the chaotic model with the reconstructed observed deformation phase space as the input variable, the optimal parameters [C, γ] of the SVM obtained by the IGWO algorithm are [3.5801, 29.4460].Figure 9
Relation curves between fitness value and iteration number: conventional model (a) and chaotic model (b).It can be seen from Figure9 that, whether it is a conventional model with dam deformation influencing factors as input variables or a chaotic model with reconstructed observation displacement data series phase space as the input variable, the IGWO algorithm reduces the number of iterations and can find the solution closest to the best goal faster. The main explanation is that the IGWO algorithm enriches the diversity of the initial population and improves the global search capability, thereby accelerating the convergence speed and convergence accuracy.


## 4.3. Results

For the PSR-SVM-IGWO-based dam observation displacement prediction model, the relevant information is introduced as follows. The SVM is at the heart of this innovative combination model. The input variable is the reconstructed phase space of the measured displacement data sequence, and the IGWO algorithm is used to realize the parameter optimization of SVM. In addition, the kernel function of SVM is a radial basis function.In order to better analyze the predictive performance of the PSR-SVM-IGWO model, the PSR-SVM and PSR-SVM-GWO models that take the reconstructed phase space of the observed displacement data sequence as input variables are established, respectively. In addition, the SVM, SVM-GWO, and SVM-IGWO models with the factors affecting dam deformation as input variables are established to explore the influence of different input variables on the accuracy of model prediction.The prediction performances of the SVM, SVM-GWO, SVM-IGWO, PSR-SVM, PSR-SVM-GWO, and PSR-SVM-IGWO models are shown in Table1 and Figure 10. From Table 1, we can see that when the influencing factors of dam deformation are used as the input variables of the model, the square correlation coefficient (R2) is ranked from large to small as SVM-IGWO model > SVM-GWO model > SVM model, mean absolute percentage error (MAPE) is ranked from small to large as SVM-IGWO model < SVM model < SVM-GWO model, and mean square error (MSE) is ranked from small to large as SVM-IGWO model < SVM-GWO model < SVM model. The SVM-IGWO model has the largest R2 of 0.9747, and the smallest MAPE and MSE are 0.0569 and 6.1547, respectively. When the reconstructed phase space of the dam deformation time series is used as the input variable of the model, R2 is sorted from large to small as PSR-SVM-IGWO model > PSR-SVM-GWO model > PSR-SVM model. The ranking of MAPE and MSE from small to large is PSR-SVM-IGWO model < PSR-SVM-GWO model < PSR-SVM model. The PSR-SVM-IGWO model has the largest R2 of 0.9902, and the smallest MAPE and MSE are 0.0423 and 2.6086, respectively.Table 1
Predictive performances of SVM, SVM-GWO, and SVM-IGWO models.
Prediction modelMSEMAPER2SVM33.12950.11480.9517SVM-GWO21.05680.13950.9616SVM-IGWO6.15470.05690.9747PSR-SVM6.66190.06180.9802PSR-SVM-GWO2.95860.04440.9842PSR-SVM-IGWO2.60860.04230.9902Figure 10
The predicted results of all models.As shown in Figure10, the measured displacement shows an obvious periodic change, but, due to the influence of creep, the peak displacement presents a gradual decrease trend. SVM, SVM-GWO, SVM-IGWO, PSR-SVM, PSR-SVM-GWO, and PSR-SVM-IGWO models can effectively predict the trend of dam displacement; and the PSR-IGWO-SVR model has the highest prediction accuracy and the smallest fluctuation range of the prediction error. For the conventional model, as time goes by, the deviation between the predicted value of the model and the observed value becomes larger and larger, while for the PSR model, the deviation is greatly improved. The calculation results show the following: ① The GWO optimized SVM algorithm is applicable in dam deformation prediction. ② The IGWO algorithm proposed in this paper has more outstanding optimization ability in optimizing SVM parameters than the conventional GWO. ③ For ultrahigh arch dams, the PSR model with the reconstructed phase space as the input variable has higher prediction accuracy and smaller prediction error than its corresponding conventional model with dam deformation influence factors as input variables.In general, SVM, SVM-GWO, and SVM-IGWO models can predict dam deformation trends well. The prediction accuracy of the PSR model with the phase space of the reconstructed observation displacement data sequence as the input variable is higher, which provides a strong support for the discussion of the chaotic phenomenon in the dam deformation time series. The calculation results also prove the applicability of the GWO algorithm in the field of dam deformation prediction. IGWO has more outstanding optimization capabilities than GWO. However, the algorithm proposed in this paper has only obtained good application results on dam deformation prediction. The application prospects of dam seepage and stress prediction are temporarily unknown, which will be tried in future research work.


## 5. Conclusions

This research proposes an innovative model combining chaos theory, support vector machine, difference algorithm, and grey wolf algorithm, namely, the PSR-SVM-IGWO model, to predict dam deformation. Taking the measured displacement data of the Xiluodu super-high arch dam during the initial operation as an example, the prediction effect of the PSR-SVM-IGWO model is compared and verified. The main conclusions are as follows:(1)
As the correlation dimension of the deformation time series tends to be saturated (D2 = 0.900), the largest Lyapunov exponent (λmax=0.001511) is greater than 0 and the Kolmogorov entropy estimate (K2 = 0.0306) is a finite positive value; it can be seen that there is chaos in the deformation observation data during the initial operation period of the dam.(2)
The optimization performance of the IGWO algorithm is superior to that of the GWO algorithm. Using the DE algorithm to ensure the initial population diversity can effectively improve the grey wolf optimization algorithm's ability to find high-quality solutions.(3)
The predictive performances of the PSR-SVM, PSR-SVM-GWO, and PSR-SVM-IGWO models with the reconstructed observation data sequence phase space as the input variable are superior to those of the corresponding SVM, SVM-GWO, and SVM-IGWO conventional models with deformation influence factors as input variables. When the conventional model predicts the deformation of an ultrahigh arch dam, although the accuracy meets the requirements, the predicted value will gradually deviate from the measured value as time goes by. On the contrary, it is difficult to observe such large deviations in models that adopt the reconstructed phase space of the measured displacement sequence as the input variable. Among all the models calculated in this paper, the PSR-SVM-IGWO model has the best prediction performance. In the next step, we will try to apply the method proposed in this paper to the prediction and analysis of seepage, stress, and strain of concrete dams.



---

*Source: 8487997-2021-11-05.xml*