A Multi-Strategy Improved Sooty Tern Optimization Algorithm for Concrete Dam Parameter Inversion

Abstract

A original strategy for optimizing the inversion of concrete dam parameters based on the multi-strategy improved Sooty Tern Optimization algorithm (MSSTOA) is proposed to address the issues of low efficiency, low accuracy, and poor optimizing performance. First, computational strategies to improve the traditional Sooty tern algorithm, such as chaos mapping to improve the initial position of the population, a new nonlinear convergence factor, the LIMIT threshold method, and Gaussian perturbation to update the optimal individual position, are adopted to enhance its algorithmic optimization seeking ability. Then, the measured and finite element data are combined to create the optimization inversion fitness function. Based on the MSSTOA, the intelligent optimization inversion model is constructed, the inversion efficiency is improved by parallel strategy, and the optimal parameter inversion is searched. The inversion strategy is validated through test functions, hypothetical arithmetic examples, and concrete dam engineering examples and compared with the inversion results of the traditional STOA and other optimization algorithms. The results show that the MSSTOA is feasible and practical, the test function optimization results and computational time are better than the STOA and other algorithms, the example inversion of the elastic modulus is more accurate than the traditional STOA calculation, and the results of the MSSSTOA inversion are reasonable in the engineering example. Compared with other algorithms, the local extremes are skipped, and the time consumption is reduced by at least 48%. The finite element hydrostatic components calculated from the inversion results are well-fitted to the statistical model with minor errors. The intelligent inversion strategy has good application in concrete dam inverse analysis.

1. Introduction

Dams are an essential part of national economic development and play a significant role in flood control, power generation, and other fields [1,2]. According to recent statistics, there are about 61,988 dams, which are spread in 166 countries [3]. In China, there are over 98,000 dams in service, of which concrete dams account for 56% of the high dams over 200 m [4,5]. Dams are subjected to complex loading from internal and external factors during operation. Forward analysis of concrete dam monitoring data can be used to understand the operational status of dams, thus effectively evaluating the real-time safety status of dams, of which displacement field monitoring is the most widespread and accessible concrete dam measurement [6,7,8,9,10].

The accuracy of the dam displacement calculation depends on the reasonableness of the selected mechanical parameters of the dam construction materials and foundation rock [11]. The mechanical parameters of dams are constantly changing with time, especially after many years of service. The actual parameter values differ significantly from the design values [12]. Some material parameters, such as the modulus of elasticity, can be directly estimated using a list of in situ tests. The common tests are the core drilling method, rebound method, stress wave propagation method, ultrasonic pulse method, and the acoustic emission method. However, in practice, it is impractical to obtain accurate mechanical parameters of concrete dams determined entirely by experimental methods [13]. And in fact, measurement is point wise and instruments can only acquire measurements within a small area around the measurement point, usually displacement-based. Inverse analysis, based on engineering-measured information to invert the physical and mechanical parameters or measured loads of materials under the assumption of known material ontological relationships, has been applied in several fields [14,15,16,17,18]. Therefore, to improve the accuracy of the values of the dam material parameters, thus laying the foundation for the accurate analysis of the deformation state of the dam, attention should be paid to the inverse analysis of the mechanical parameters of the dam foundation based on the measured data.

Dam inversion analysis is a complex nonlinear space search problem with a multi-parameter combination [19], whose objective function is often a nonlinear multi-peak function, and numerous scholars have used traditional optimization algorithms to achieve the global optimization search of dam material parameters. Kang et al. [20] proposed a displacement inversion analysis model based on Gaussian Process Regression and applied the improved sparrow search algorithm to the identification of parameters of arch dams. Gu et al. [21] used a chaotic genetic optimization algorithm to invert the zonal modulus of concrete dams. Zhu et al. [22] conducted a study on applying the Quantum Genetic algorithm in the inversion of zonal mechanical parameters of extra-high dams. Dou et al. [23] introduced adaptive factors to improve the fireworks algorithm for parameter identification of concrete dams and improved the accuracy and efficiency of inverse analysis of concrete dams. Chen et al. [24] combined improved cuckoo search and improved particle swarm optimization to construct a hybrid optimization algorithm for identifying the material parameters of partitioned concrete dams and achieved good results. Han et al. [25] constructed a particle swarm algorithm-based inversion model for dam foundation rock parameters and verified the reasonableness of the method. Lin et al. [26] used the GWO optimization algorithm to search and identify the viscoelastic parameters of concrete dam foundations, which is very competitive in parameter inversion and a real-time problem approach. Chen et al. [27] and Liu et al. [28] both built a machine learning-based inversion model of a deformation modulus for high-arch dams and identified the mechanical parameters of the concrete dams and their foundation rock body by GSA and WOA, respectively, and they all achieved reasonable inversion results.

In studying the above traditional optimization algorithm for the inversion of dam parameters, there are still two areas for improvement. First, the inversion accuracy of traditional optimization algorithms is still to be improved, and problems of premature convergence and poor robustness will occur when solving complex problems; second, the use of the traditional serial computation strategy leads to problems of a time-consuming, low computational efficiency. There is an urgent need to introduce novel intelligent optimization algorithms and improve the finite element orthogonal computational architecture to solve the above problems. The Sooty tern algorithm is a new type of optimization algorithm of the meta-inspired method, with solid optimization ability and fast convergence speed, which has achieved specific applications in several fields [29,30]. In order to improve the accuracy and efficiency of solving complex engineering practical problems such as dam parameter inversion, a multi-core parallel computing strategy is introduced into the traditional Sooty Tern Optimization algorithm in this paper, and several strategies are used to improve the traditional Sooty Tern Optimization and propose the multi-strategy-based improved Sooty Tern Optimization (MSSTOA). On this basis, the proposed improved algorithm is applied to test functions and the inversion of mechanical parameters in arithmetic and engineering examples, and the inversion results are compared with those obtained from some of the aforementioned intelligent algorithms. It is found that the time cost for inversion is reduced dramatically, and the accuracy of inversion is superior to those of the aforementioned intelligent algorithms, which then verifies the feasibility of the inversion strategy.

The structure of this paper is as follows. Section 2 introduces the theory of parameter inversion for concrete gravity dams. Section 3 presents the inverse analysis framework proposed in this paper. Section 4 and Section 5 provide method arithmetic examples and example validation results, respectively. Conclusions are provided in Section 6.

2. Basic Principles of Parametric Inversion for Gravity Dams

The displacement of any point of the dam body can be regarded as the structural response under multiple factors such as water pressure, ambient temperature, etc. Électricité De France (EDF) proposed the hydrostatic seasonal time (HST) model in the 1960s. By establishing the HST model for the displacement of each measurement point, the displacement can be classified into three components, namely, hydrostatic pressure, temperature, and time effect [31], as shown in the following equation:

$$\begin{array}{ll}\delta & ={\delta }_H+{\delta }_T+{\delta }_{\theta }\\ & =a_0+{\displaystyle \sum _{i=1}^3{a}_i{H}^i+}{\displaystyle \sum _{i=1}^2(b_{i1}\sin \frac{2\pi it}{365}+b_{i2}\cos \frac{2\pi it}{365})+}{c}_1\theta +c_2\ln \theta \end{array}$$

(1)

where δ_H is the hydrostatic component, δ_T is the temperature component, and δ_θ is the aging component. a_i, b_i, c_i are the coefficients to be determined for the statistical model, H is the water depth upstream of the dam with a unit of m, t is the number of days from the initial monitoring date, θ = t/100. c₁θ + c₂lnθ is the time effect component, ${\displaystyle \sum _{i=1}^2(b_{i1}\sin \frac{2\pi it}{365}+b_{i2}\cos \frac{2\pi it}{365})}$ is the temperature component, and $a_0+{\displaystyle \sum _{i=1}^3{a}_i{H}^i}$ is the required hydrostatic component.

The hydrostatic pressure component is related to the values of the material mechanical parameters of the concrete dam and foundation rock [32]. In addition, the hydrostatic pressure component captures the elastic displacement behavior due to the water level variation, which coincides with the water level variation in the reservoir, and the hydrostatic pressure component can be easily obtained from the finite element model. It is only to compute the displacement response between different water levels, so the separated hydrostatic component δ_H is used to invert the mechanical parameters. It is assumed that the material parameters to be inverted are Z, which contains M parameters. The principle of inversion is to take the parameters to be inverted as input values and consider the displacement δ as the best approximation between the finite element calculation value and the measured material, that is, the inversion objective function J(Z), which transforms the inversion problem into an optimization problem and establishes the objective function of inversion for multiple measurement points as follows:

$$\left\{\begin{array}{l}\mathit{K}\mathit{u}=\mathit{R}\\ \min J(\mathit{Z})=\frac{{\displaystyle {\sum }_{i=1}^N{\displaystyle {\sum }_{j=1}^M{\left(\mathsf{\Delta }{\delta }_{ij}^{\mathrm{measure}}-\mathsf{\Delta }{\delta }_{ij}^{\mathrm{inverse}}\right)}^2}}}{N\times M}\\ \mathit{Z}=[Z_1,Z_2,\dots ,Z_M{]}^T\\ \mathrm{s}.\mathrm{t}.\begin{array}{cc}& \end{array}{Z}_{min}\le Z_k\le Z_{max}\end{array}\right.$$

(2)

where K is the stiffness matrix, u is the displacement vector, R is the load vector, J(·) is the objective function, Z is the parameter vector to be inverted, Z_k is the i-th material parameter to be inverted, Z_max and Z_min are the upper and lower limits of Z_k, N is the number of displacement measurement points of the dam body, and M is the number of the selected typical water level. $\mathsf{\Delta }{\delta }_{ij}^{\mathrm{measure}}$ is the difference of the separated hydrostatic component of the measured data of the i-th point, and $\mathsf{\Delta }{\delta }_{ij}^{\mathrm{inverse}}$ is the difference of the calculated component of hydraulic pressure of the same point. The calculation expression is as follows:

$$\left\{\begin{array}{c}\mathsf{\Delta }{\delta }_{ij}^{\mathrm{measure}}={\delta }_{ij}^{\mathrm{measure}}-{\delta }_{i0}^{\mathrm{measure}}\\ \mathsf{\Delta }{\delta }_{ij}^{\mathrm{inverse}}={\delta }_{ij}^{\mathrm{inverse}}-{\delta }_{i0}^{\mathrm{inverse}}\\ {\delta }_{ij}^{\mathrm{inverse}}=f(H_j,\mathit{Z})\end{array}\right.$$

(3)

where ${\delta }_{ij}^{\mathrm{measure}}$ is the measured water pressure component at point i under the typical water level H_j, ${\delta }_{i0}^{\mathrm{measure}}$ is the measured water pressure component at point i under the initial water level H₀, ${\delta }_{ij}^{\mathrm{inverse}}$ is the calculated water pressure component at point i under the characteristic water level H_j, and ${\delta }_{i0}^{\mathrm{inverse}}$ is the calculated water pressure component at point i under the initial water level H₀.

3. A Multi-Strategy Improved STOA Inversion Framework

3.1. STOA

The Sooty Tern Optimization algorithm is a search optimization algorithm proposed by Dhiman et al. 2019 by simulating the feeding behavior of sooty terns [33] which consists of two phases: a global search phase that simulates terns’ migration and a local search phase in which the terns circle and attack their prey.

As with other population algorithms, each sooty tern as represents a search agent, and all search agents together construct the algorithm’s solving matrix X. In the initialization phase, the algorithm starts by constructing an initial solving matrix in the pre-set search space.

$$\mathit{X}={\left[\begin{array}{c}{X}_1\\ ⋮\\ X_i\\ ⋮\\ X_N\end{array}\right]}_{N\times m}={\left[\begin{array}{ccccc}{x}_{11}& \cdots & x_{1j}& \cdots & x_{1m}\\ ⋮& \ddots & ⋮& \ddots & ⋮\\ x_{i1}& \cdots & x_{ij}& \cdots & x_{im}\\ ⋮& \ddots & ⋮& \ddots & ⋮\\ x_{N1}& \cdots & x_{Nj}& \cdots & x_{Nm}\end{array}\right]}_{N\times m}$$

(4)

where X is the population solving matrix of the Sooty Tern algorithm; X_i is the solving vector of a search agent, which represents a random initial solution in the initial stage and is updated during the iterative computation to the iterative solution at current iteration times, i is the index of search agent; x_ij is the value of a search agent in a dimension of the search space; j is the dimension index; N is the number of populations of the algorithm; and m is the dimension of the search space.

3.1.1. Migration Behavior (Global Search)

Sooty terns engage in seasonal migratory behavior in search of abundant food sources. The purpose of this phase is to quickly identify the optimal area by randomly searching in the search space for the global search phase. This phase consists of three parts: conflict avoidance, aggregation, and position update.

• Conflict avoidance: To avoid collisions between individuals during migration, additional mass S_A is introduced in the iterative computation to update individual positions.

  $$C=S_A\times X\left(t\right)$$

  (5)

  where X(t) is the position of the search agent at a particular iteration, t is the current iteration index, C is the position in the case of no collision with other individuals, and S_A is an additional variable used to avoid collision, which is calculated as follows:

  $$S_A=C_f\times (1-t/T)$$

  (6)

  where C_f is the control constant used to adjust S_A, which is generally set to 2 [33], and T is the maximum number of iterations. Therefore, S_A will linearly decrease from 2 to 0.

2.

Aggregation: After avoiding collisions between neighboring agents, the search agent will move towards the best position among the neighboring agents, that is, towards the position of the optimal solution, as expressed in the following equation:

$$M=C_B\times \left(X_{best}(t)-X(t)\right)$$

(7)

where M is the process of moving X at different locations towards the location X_best of the optimal solution, and C_B is a random number used to make the exploration more comprehensive. The equation is as follows:

$$C_B=0.5\times Rand$$

(8)

where Rand is a random value between (0, 1).

3.

Position update: The search agent updates its position based on the best position. The equation is as follows:

$$D=C+M$$

(9)

where D is the distance between the current individual position and the global optimal position.

3.1.2. Migration Behavior (Local Search)

When sooty terns need to attack their prey during migration, they will hover in a spiral formation in the air. The description is as follows:

$$\left\{\begin{array}{l}{x}^{\prime }=R\times \sin (i)\\ y^{\prime }=R\times \cos (i)\\ z^{\prime }=R\times i\\ R=u\times e^{kv}\end{array}\right.$$

(10)

where R is the radius of each spiral and i is a random variable which lies between the range of (0, 2π). u and v are constants defining the shape of the spirals and take the value of 1 [33]. The search agent position update equation is as follows:

$$X(t)=\left(D\times (x^{\prime }+y^{\prime }+z^{\prime })\right)\times X_{best}\left(t\right)$$

(11)

The pseudo-code for the STOA is shown in Figure 1a.

3.2. An Improved STOA with Multiple Strategies

For the shortcomings of the traditional STOA, this paper introduces circle mapping to improve the distribution quality of the initial solution, improves the additional quality to balance the global and local search ability of the algorithm, introduces limit threshold and Gaussian variation to avoid the algorithm from jumping into the trap of local minima when searching for the optimal value in the local search, and improves the computational efficiency of the algorithm based on parallel computation. Based on the above improvement strategies, the MSSTOA is proposed for inverse computation.

3.2.1. Population Initialization by Circle Mapping

The spatial distribution of the initial solution set in metaheuristic algorithms is an essential factor that affects the global search speed and optimization results of the algorithm. The initialized solution set generated in the STOA based on the random number strategy is difficult to uniformly distribute in the search space, which can lead to a reduction in the search efficiency of the algorithm [34]. Chaotic mapping has the characteristics of randomness, traversal, and regularity, which can be used to initialize the positions of STOA individuals using chaotic series to avoid the algorithm falling into local extremes. In this paper, circle chaotic mapping is used for population initialization. Circle mapping produces chaotic sequence expression, as shown in Equation (12).

$$nu{m}_{i+1}=\mathrm{mod}(nu{m}_i+0.2-0.5/2\pi \times \sin (2\pi \times nu{m}_i),1)$$

(12)

where num_i denotes the value of the i-th chaotic series, and mod denotes the residual operation.

Assuming that the search space domain of the optimization objective problem is [Z_min, Z_max] and num_ij is the value of the chaotic series of the circle mapping, the initial solution vector X is expressed as follows:

$$x_{ij}=Z_{\min }+(Z_{max}-Z_{min})\times nu{m}_{ij}$$

(13)

3.2.2. Non-Linear Additional Mass S_A

It is crucial to balance the two phases of global and local search of the meta-heuristic optimization algorithm. As can be seen from Equation (6), S_A is linearly decreasing and is unable to cope with complex nonlinear problems in the actual optimization process, which will reduce the ability of the algorithm to find the optimal search. For this reason, a nonlinear additional mass based on sinusoidal function is adopted, and the improved S_Anew is defined as

$$S_{Anew}=1-\sqrt{2}\sin (\frac{\pi }{2}\times {(\frac{t}{T})}^4-\frac{\pi }{4})$$

(14)

where t is the current number of iterations of the algorithm, and T is the set maximum number of iterations.

As shown in Figure 2, compared with the traditional STOA, the improved nonlinear added mass shows a nonlinear trend of a slight decrease in the initial period and sudden decrease in the later period as the number of iterations increases, which indicates that the individual explores the whole search space powerfully in the early period, which improves the global searching ability of the algorithm; the later period implies a faster convergence of the algorithm. Compared with the original algorithm, the nonlinear additional mass better balances the global and local search ability.

3.2.3. Limit Threshold and Gaussian Variation

In the iteration of the STOA, agents will rapidly assimilate due to aggregation to the optimal individual, thus shrinking the population distribution and making it easy to appear as the local optimum. To deal with this problem, the limit threshold mechanism is adopted, a limited number of values is set before the iteration starts, and it is observed whether the current optimal solution is stagnant and unchanged within the limited number of times and if the optimal solution does not change, the Gaussian mutation strategy is applied. This strategy is to perturb the location of the optimal solution individual to obtain a new location and compare the fitness of the two locations before and after. Furthermore, the location of the updated individual with smaller fitness is selected, which aims to improve the ability of the algorithm to jump out of the local minimum. The mathematical model of this strategy is expressed as follows:

$$x_m=x_{ij}^{}+x_{ij}^{}\times \mathrm{Gauss}(0,1)$$

(15)

where X_m is the position of the population after mutation. Gauss(0, 1) is the Gaussian operator.

The pseudo-code of the MSSTOA is shown in Figure 1b.

3.2.4. Algorithm Parallelism Improvement

The traditional STOA uses a continuous computation strategy to simulate population feeding behavior. That is, for each iteration step, when one individual finishes feeding, the next one can start foraging. For the parameter inversion problem of large-scale engineering, the analysis process needs to call the finite element operation repeatedly, and the use of continuous computation means that the computer can only carry out finite element analysis once at the same moment, which is highly time-consuming. The parallel computing strategy has the same number of computational threads as the number of computer kernels by dividing the original population into sub-populations with a comparable number of kernels and assigning them to each kernel to perform the finite element analysis computation at the same time to be merged to discover the optimal solution of this iteration after the end of the sub-population computation and then enter into the next iteration step. Through parallel improvement, the space complexity of the algorithm can be increased and the time complexity can be reduced, thus improving the algorithm’s efficiency.

3.3. The Process of Inverse Optimization Implementation

Combining the above basic principle of optimized inversion of concrete dam parameters and the multi-strategy improved STOA, the optimization inversion framework of the mechanical parameters of concrete dams based on the MSSTOA is proposed. Its implementation process is shown in Figure 3, and the specific steps are as follows:

Step 1: Combine the measured data, separate the water pressure component according to the statistical model, and establish the finite element model according to the engineering data.

Step 2: Select the mechanical parameter Z to be inverted, construct N sooty tern individuals, set the iteration number T, termination conditions, and parameters, including the population size N, the upper and lower bounds of the search lb and ub, and initialize the population by circle mapping chaotic series within the search domain of parameter Z.

Step 3: Allocate sub-populations according to the principle of parallel computing, call the finite element software to calculate the corresponding objective function value of each sooty tern according to Equation (2), respectively, and compare them to discover the optimal value of the population.

Step 4: Calculate and update the position of the sooty tern according to Equations (5)–(15).

Step 5: Determine whether the iteration termination condition is reached; if so, output the global optimal value, namely the inversion result, and end the procedure; otherwise, update the individual sooty tern position and repeat Steps 3–4.

4. Algorithm Verification

4.1. The Performance on the Test Functions

Before applying the MSSTOA for concrete dam inversion, the performance of the MSSTOA needs to be tested. The proposed algorithms are compared: the Sooty Tern Optimization algorithm (STOA), Grey Wolf Optimization (GWO), Gravitational Search Algorithm (GSA), Particle Swarm Optimization (PSO), and Whale Optimization algorithm (WOA), and the test functions are the Sphere, Rosenbrock, Rastrigin, Ackley, Griewank, and Schwefel functions from the benchmarking function set of the Congress on Evolutionary Computation (CEC) [35] (see

Table 1. Test function.

Test Function	Search Range	Formula	fmin
Sphere	(−100, 100)	$f\left(x\right)={\displaystyle \sum _{i=1}^n{x}_i^2}$	0
Rosenbrock	(−30, 30)	$f\left(x\right)={\displaystyle \sum _{i=1}^{n-1}\left[100\left(x_{i+1}-x_i^2\right)+{(x_i-1)}^2\right]}$	0
Rastrigin	(−5.12, 5.12)	$f\left(x\right)={\displaystyle \sum _{i=1}^n\left[x_i^2-10\cos \left(2\pi x_i\right)+10\right]}$	0
Ackley	(−32, 32)	$f\left(x\right)=-20\exp (-0.2\sqrt{\frac{1}{2}{\displaystyle \sum _{i=1}^n{x}_i^2}})-\exp \left(\frac{1}{2}{\displaystyle \sum _{i=1}^n\cos \left(2\pi x_i\right)}\right)+20+e$	0
Griewank	(−600, 600)	$f\left(x\right)=\frac{1}{4000}{\displaystyle \sum _{i=1}^n{x}_i^2-{\displaystyle \prod _{i=1}^n\cos \left(\frac{x_i}{\sqrt{i}}\right)}}+1$	0
Schwefel 1.2	(−100, 100)	$f\left(x\right)={\displaystyle \sum _{i=1}^n{\left({\displaystyle \sum _{j-1}^i{x}_j}\right)}^2}$	0

). The maximum number of iterations and population size of all algorithms are 500 and 30, respectively. Each algorithm is run independently 30 times to avoid algorithmic chance. All algorithms shown in this subsection are the serial version. The mean and standard deviation are taken to judge the solving accuracy and robustness.

Table 2. Comparison results of different algorithms on CEC test functions.

Test Function		MSSTOA	STOA	GWO	GSA	PSO	WOA
Sphere	M	2.89 × 10−249	6.42 × 10−90	1.67 × 10−27	1.16 × 102	9.17	1.5 × 10−72
	SD	3.77 × 10−250	2.80 × 10−89	1.59 × 10−27	1.77 × 102	3.78	6.49 × 10−72
	TC	0.09	0.11	0.18	0.36	0.14	0.12
Rosenbrock	M	9.32 × 10−3	1.59 × 10−2	27.20	3.46 × 102	1.46 × 103	27.9
	SD	1.31 × 10−2	1.05 × 10−2	6.52 × 10−1	2.48 × 102	9.12 × 102	4.84 × 10−1
	TC	0.10	0.16	0.18	0.36	0.19	0.12
Rastrigin	M	0.00	0.00	3.72	38.10	99.20	2.84 × 10−15
	SD	0.00	0.00	4.45	10.90	18.30	4.45
	TC	0.09	0.13	0.19	0.34	0.17	0.11
Ackley	M	8.88 × 10−16	1.07 × 10−15	1.02 × 10−13	4.43 × 10−1	6.16	4.80 × 10−15
	SD	3.92 × 10−20	7.74 × 10−15	1.98 × 10−14	6.47 × 10−1	1.11	1.91 × 10−15
	TC	0.09	0.14	0.18	0.35	0.16	0.12
Griewank	M	0.00	1.33 × 10−2	2.76 × 10−3	1.05 × 102	4.35 × 10−1	5.55 × 10−18
	SD	0.00	3.47 × 10−18	5.76 × 10−3	14.6	1.54 × 10−1	2.42 × 10−17
	TC	0.10	0.17	0.17	0.38	0.19	0.13
Schwefel 1.2	M	3.76 × 10−261	1.43 × 10−80	2.39 × 10−5	1.49 × 103	7.87 × 102	4.48 × 104
	SD	5.87 × 10−280	2.90 × 10−80	4.15 × 10−5	5.65 × 102	3.33 × 102	1.02 × 104
	TC	0.15	0.39	0.28	0.45	0.23	0.19

Note: Abbreviations: M, Mean; SD, Standard Deviation, TC, Time Cost (s).

and Figure 4 show the test results and average fitness convergence curves, respectively. As seen in Table 2, the results of the MSSTOA computation outperform that of the initial STOA and the other algorithms. In the same test function, the proposed algorithm consumes the least time compared to all the other algorithms, the GSA algorithm consumes the longest time, and except for the GSA algorithm, the difference in time consumption between the serial versions of the algorithms is not very significant. In Figure 4, the MSSTOA shows excellent global search ability and the fastest convergence speed. The test results prove that the defects of the original algorithm have been improved, and the proposed algorithm effectively improves the global search ability and accelerates the convergence speed.

4.2. The Performance on the Assumed Example

An ideal two-dimensional concrete gravity dam is then used as an example to verify the performance of the proposed inversion algorithm. The height of this concrete gravity dam is 113.0 m. The finite element model of its right bank dam section consists of 841 elements and 921 nodes. There are two points in the vertical direction of the cross-section, which are point O within the dam foundation and the inversed pump line measurement point P1, as shown in Figure 5. The foundation extends two times the height of the dam both upstream and downstream, as well as at the foundation depth. The X-direction is the downstream direction, and the Y-direction is the vertical direction. Normal constraints are taken at the sides and fixed constraints at the bottom.

Assuming that the dam is in the elastic change stage, two sets of elastic moduli to be inverted are set. Poisson’s ratio has a negligible effect on dam deformation [36] and is set as the default value in the inversion. The specific inversion parameters are designed as shown in

Table 3. Model material parameters.

Sets	Materials	Density (kg·m−3)	Assumed E (GPa)	Search Range of E (GPa)	Poisson’s Ratio
1	concrete	2400	30.0	15~40	0.167
	rock	2790	15.0	10~30	0.2
2	concrete	2400	24.0	15~40	0.167
	rock	2790	12.0	10~30	0.2

. Water load was considered for loading, and the relative displacement of measurement point P1 concerning point O inside the dam foundation was taken as the calculated displacement of P1. The upstream water depth of 93 m was taken as the initial water level and 103 m as the characteristic water level, the calculated displacements of point P1 under two sets of hypothetical elastic modes and two characteristic water levels are obtained, respectively, and the specific calculation results are shown in

Table 4. Calculated displacements at point P1 for different parameter sets.

Sets	Water Level (m)	X-Direction Displacement (mm)	ΔX (mm)	Y-Direction Displacement (mm)	ΔY (mm)
1	93	0.97	/	0.73	/
	103	1.67	0.70	1.25	0.52
2	93	1.21	/	0.91	/
	103	2.08	0.87	1.56	0.65

.

All calculations after this section were completed in MATLAB R2023a and ABAQUS 2022 via the same station computer manufactured by DELL, Round Rock, TX, America, with detailed technical parameters shown as two CPUs of E5-2678V3 and twelve core processors, RAM of 64 GB, and 64-bit system type, and the number of kernels of the adopted computers for parallelism computing was 4.

Parameter inversion was carried out using the MSSTOA and STOA, respectively. And Figure 6 and Figure 7 show the inversion iterative process and the elastic modulus change process, respectively. From Figure 6 and Figure 7, it can be seen that the fitness of the MSSTOA is close to 0 when the number of iterations reaches 20 and 37, after which the fitness is further reduced due to the LIMIT jump rule. On the contrary, in the STOA, after an initial rapid decrease in fitness, the computational fitness does not change or need more iterations to obtain a better fitness, and the corresponding inversion results are also at odds with the assumed parameters. The two sets of inversion results of the MSSTOA are (1) Ec = 30.0 GPa, Er = 15.0 Gpa and (2) Ec = 24.0 Gpa, Er = 12.0 Gpa, which are consistent with the preset parameters, indicating that the algorithm can be used for the inversion of concrete dam parameters.

5. Case Study

An actual project is selected for case analysis to further verify the practicality of the inverse analysis method. The actual project is a raised concrete gravity dam (Figure 8a). The original dam crest elevation is 162 m, the dam height is 97 m, and the dam is raised to increase the reservoir capacity. Construction began in September 2005 and was completed in March 2010. The height of the dam after raising is 176.6 m, the maximum height of the dam is 111.6 m, and the foundation rock is a metamorphic magma rock [37]. The evolution of the material properties of old and new concrete dams has an essential effect on the performance of the raised dams. It is unreasonable to treat concrete dams as a whole to calculate the integrated modulus. It is necessary to calculate the material mechanical parameters of concrete dams in a partitioned manner. Taking the 13# spillway dam section in the middle of the riverbed as an example, a plumb line is arranged in the dam body to monitor the horizontal displacement of the dam section (Figure 8b), point O is the suspension point of the plumb line, and P1~P5 are the measurement points of the plumb line. The inversion of the material mechanical parameters of the old and new dams and the dam foundation during the operation period is carried out using the actual measured horizontal displacements.

The monitoring of the plumb line started on 1 June 2013, and the stepwise regression analysis of the data from P1 to P5 was carried out to separate δ_H, δ_T, and δ_θ. Considering the sensitivity of the inversion, using the data from the measurement points with a high fitting degree can be more conducive to the accurate inversion of the dam material parameters. The horizontal displacement of the gravity dam in the cross-river direction has a small measurement value with a large error relatively, and at the same time, points P4 and P5 have a low fitting degree. Therefore, the horizontal displacements along the river from measurement points P1 to P3, which are located in the middle and upper part of the dam body, are selected for the parameter inversion. The fitting degrees and component coefficients of the three points are shown in

Table 5. Component coefficients and fitting of statistical model for measurement points.

Points	R	a0	a1	a2	a3	b1	b2	b3	b4	c1	c2
P1	0.9360	−0.5126	9.8443	−0.1452	7.16 × 10−4	−2.6507	−0.8143	0	0	0	−0.4797
P2	0.9553	−0.3813	11.6993	−0.1756	8.75 × 10−4	−2.5646	−1.3217	0	0.3649	0	−0.5857
P3	0.9578	0.076	6.9122	−0.1014	4.95 × 10−4	−2.486	−1.4652	0	0.2188	0	−0.4452

, and the results of the separation of the displacement components are shown in Figure 9. The monitoring data of the starting measurement day are selected as the initial condition of the finite element. The recent two days with a high fit in time series are set as the typical days. The corresponding hydrostatic components’ displacement difference between them and the initial condition is the target of the inversion analysis, shown in

Table 6. Target values for dam inversion analyses.

Points	Typical Day	Upstream Water Level (m)	Measured Value (mm)	Fitted Value (mm)	Water Load Component (mm)	Component Difference (mm)
P1	1 June 2013	138.43	0	0	0	/
	26 November 2017	166.03	2.58	2.52	5.07	5.07
	27 December 2017	165.52	3.04	3.14	4.45	4.45
P2	1 June 2013	138.43	0	0	0	/
	26 November 2017	166.03	1.69	1.81	4.00	4.00
	27 December 2017	165.52	2.34	2.44	3.64	3.64
P3	1 June 2013	138.43	0	0	0	/
	26 November 2017	166.03	1.40	1.41	1.77	1.77
	27 December 2017	165.52	2.02	1.98	1.62	1.62

.

A finite element analysis model of the 13# dam section was established. As shown in Figure 10, the whole model consists of 135,232 elements and 146,122 nodes, the foundation range extends upstream and downstream, and the depth of the foundation is about two times the maximum height of the dam, respectively, with the X direction directed to the downstream, the Y direction directed to the left bank, and the Z direction vertically upward. Zone A is the raised concrete dam, zone B is the original dam section, and zone C is the dam foundation. Normal constraints are taken at the side and fixed constraints at the bottom. The dam concrete and foundation materials are taken as design values except for the parameters to be inverted; the concrete’s Poisson’s ratio is 0.167 and the foundation rock’s Poisson’s ratio is 0.22. The loads consider upstream water pressure and lift pressure. The displacements of points P1~P3 relative to the suspension point O at the top of the dam are taken as the calculated displacement values of P1~P3 in the numerical calculations, which can be used to simulate the hydrostatic component δ_H in the displacements of the concrete dam since there is only the effect of hydraulic pressure in the external loads.

According to the related information, the values of the concrete and bedrock elastic modulus are taken in the range of Ec = 15~40 GPa and Er = 10~30 GPa. Taking Equation (2) as the objective function, the material parameters of concrete dams at the optimum fitness are solved by the GWO, GSA, PSO, STOA, and MSSTOA, respectively. The benchmark algorithm uses serial computation, and the MSSTOA additionally uses parallel computation to compare the computation time. The maximum number of iterations is 100. The optimal fitness change curves and the material parameter identification process for the MSSTOA parallel computation are shown in Figure 11. Additionally, the results of the benchmark algorithm for comparison are all shown in Figure 11. The inversion results of the five algorithms are shown in

Table 7. Comparison of elastic modulus inversion results for raised gravity dams.

		MSSTOA	STOA	GWO	GSA	PSO
EA (GPa)		36.16	36.17	33.77	35.04	37.71
EB (GPa)		30.18	30.49	29.26	29.32	30.30
EC (GPa)		24.45	24.12	26.29	25.86	24.00
Optimal fitness		5.37 × 10−31	9.99 × 10−7	1.23 × 10−6	7.99 × 10−6	7.91 × 10−6
Cost time (s)	parallel	34,014	/	/	/	/
	serial	61,225	65,758	91,359	94,093	70,047

.

From Figure 11 and Table 7, it can be seen that the proposed algorithm has a faster convergence speed relative to the benchmark algorithm. At the same time, when the algorithm is at the local optimal solution, it will jump out according to the LIMIT rule, and the fitness fluctuates. It continues to converge until it obtains the global optimal solution. In contrast, other algorithms will fall into the local optimal trap and keep searching for the optimal solution in the local range. It is proved that the inverse analysis strategy proposed in this paper is feasible and effective, and the inverse fitness is better than other algorithms. Meanwhile, due to the parallel computing strategy, the computation time cost is reduced by 48%~64% relative to other benchmark optimization algorithms.

The inversion identification parameters are substituted into the finite element for calculation. The obtained finite element water pressure components are compared with the statistical water pressure components. The results of the water pressure component and finite element water pressure component calculations for 2017 are shown in Figure 12 and the comparison results for some typical days are shown in

Table 8. Comparison of positive analyses of inversion results for raised gravity dams.

Points	Date	Upstream Water Level (m)	δHmeasured (mm)	δHFEM (mm)	Relative Error (%)
P1	1 January 2017	154.54	1.56	1.568	0.51
	26 September 2017	162.47	3.28	3.296	0.49
	26 November 2017	166.03	5.07	5.098	0.55
	27 December 2017	165.52	4.45	4.437	0.29
P2	1 January 2017	154.54	1.30	1.305	0.39
	26 September 2017	162.47	1.95	1.956	0.31
	26 November 2017	166.03	4.00	3.989	0.28
	27 December 2017	165.52	3.64	3.633	0.19
P3	1 January 2017	154.54	0.62	0.622	0.33
	26 September 2017	162.47	1.39	1.395	0.36
	26 November 2017	166.03	1.77	1.776	0.34
	27 December 2017	165.52	1.62	1.627	0.43

. As can be seen in Figure 12 and Table 7, the relative error between the measured displacement water pressure components and the finite element calculation results is within 0.6%, and the finite element results in the figure match well with the measured displacement water pressure components with small errors, which indicates that the inversion results are reasonable.

6. Conclusions

To improve the accuracy and efficiency of the inversion algorithm, the MSSTOA is proposed in this paper. Combining the statistical model of concrete dam displacement measured data and the MSSTOA, the inversion analysis method of concrete dam mechanical parameters is presented to calculate the integrated deformation modulus and zonal deformation modulus of concrete dams. The objective function was established using horizontal incremental displacement and examples verified the method. The main conclusions are as follows:

(1)

The proposed improved algorithm shows good global search capability and convergence speed by enhancing the STOA with multiple strategies and setting the jump-out rule. It provides the possibility to eliminate local minima.

(2)

Combined with finite element computation, the inversion framework of concrete dam mechanics parameters based on the MSSTOA is constructed. Validated by two examples, the concrete dam material mechanical parameters can be effectively identified, and its inversion results are better than other benchmark algorithms, indicating that the inversion strategy has high search accuracy and fast inversion speed. Meanwhile, based on multi-core CPUs to subdivide populations in a sub-population manner for computation, it dramatically improves the solution rate of complex inversion problem computations.

(3)

The identified material parameters are used for finite element prediction of displacements, and the results are in good agreement with the elastic hydrostatic component separated by the statistical model, indicating that the method can identify the mechanical parameters related to the hydrostatic component when the dam is in an elastic state, while the proposed method can be adopted for the inversion of most mechanical responses.