Ordinary Kriging Interpolation Method Combined with FEM for Arch Dam Deformation Field Estimation

Abstract

The deformation characteristic of the arch dam can directly reflect its service performance, which can be analyzed on the basis of the dam deformation field. However, restricted by the limited number of dam monitoring points and the inhomogeneity of materials, an accurate measurement of arch dam deformation field is difficult to estimate by using the existing common methods, such as the spatial interpolation methods and the finite element method (FEM). With the aim of ensuring arch dam structure safety, the ordinary kriging interpolation method, combined with FEM, is proposed for arch dam deformation field estimation, in this study. Given the inversion of the computation parameters of the arch dam, FEM is used to calculate the basic arch dam deformation. Subsequently, the ordinary kriging interpolation method is introduced to estimate the spatial variance deformation at each point of the arch dam. One superhigh arch dam in China is selected as a case study; two additional methods are introduced as comparisons that are based on a numerical experiment and the actual monitoring data. The experimental results show that the proposed method considerably improves the accuracy and computational efficiency of the arch dam deformation field estimation and that it is of great practical importance for characterizing the deformation behavior of arch dams.

1. Introduction

An arch dam serves as the main component of a hydropower project, which substantially contributes to national, social, and economic development. Meanwhile, the service performances of arch dams have a significant impact on people’s lives and property [1,2,3]. Deformation can intuitively reflect the arch dam performance, and continuous deformation monitoring data in space can help to analyze and forecast the operational status of arch dam structure safety in real time [4,5]. Nevertheless, the number of dam monitoring points is limited, which is caused by some financial, logistical, and geological factors [6]. Therefore, it is impossible to directly acquire the spatially continuous deformation information of the arch dam only on the basis of limited monitoring data. For the purpose of accurately characterizing the arch dam deformation behavior, studying a proper deformation field estimation method is of great significance in such a case to estimate the deformation of unmonitored locations of the arch dam according to the available monitoring deformation values at surrounding locations, which has become a hot issue worthy of further study.

Recently, scholars have demonstrated that the methods for estimating the deformation field of dams include mainly methods based on monitoring data and methods based on physical mechanisms [7]. The methods based on monitoring data include some professional monitoring techniques and data interpolation methods. To characterize the overall dam spatial deformation, some engineering monitoring techniques have been applied in dam deformation field estimation [8,9,10,11]. However, the deformation field monitoring results of the above techniques always have the disadvantage of low precision. Data interpolation methods, such as spatial interpolation methods and temporal interpolation methods [12,13,14], are generally used to construct the dam deformation field. Spatial interpolation methods utilize the relationship between the deformation at the unmonitored locations and monitored locations to estimate the dam deformation field, such as the deterministic method (e.g., inverse distance weighting method (IDW)) [15] and the kriging interpolation method and its improved methods [16,17,18]. The kriging interpolation methods outperform the deterministic methods in that they account for the spatial correlation between the deformation of neighboring monitored locations and unmonitored locations [19,20]. However, the above interpolation methods all neglect the physical mechanism of dam deformation [21], so the accuracy of the deformation field estimation results encounter difficulties in trying to meet the requirements of dam safety monitoring.

Because the calculation results of the finite element method (FEM) are based on the physical mechanism, FEM has been broadly used to simulate dam deformation. In order to precisely evaluate the operation performance of arch dams, the arch dam deformation field can be simulated by an FEM that is based on the parameter inversion methods [22,23,24]. The objective spatial heterogeneity of dam material parameters [25,26,27] makes it difficult for the FEM to accurately characterize the actual deformation field of the dam.

In addition, some scholars have applied the Bayesian method to spatially varying parameter inverse analyses [28,29]. On the basis of the estimation results of parameter fields, the deformation field can be simulated. However, the results of the above studies indicate that the accuracy of estimation results of the Bayesian method depends mostly on the monitoring frequency and number of the monitoring data. In practical projects, the monitoring points and monitoring frequency of the arch dam are limited. Consequently, the Bayesian method is not well suited for arch dam deformation field estimation.

To compensate for the shortcomings of the above methods, it seems that the most reasonable method for estimating the dam deformation field should be a combination of spatial interpolation and FEM. Given that dam deformation has a certain spatial correlation at different positions [30], the ordinary kriging interpolation method, combined with FEM, is introduced in the method proposed, to simultaneously assess the reliability of the monitoring data and the distribution of the deformation in the spatial region of arch dams. FEM is used to calculate the basic dam deformation at each point of the arch dam in the ideal case, where the distribution of dam material is homogeneous. The ordinary kriging interpolation method is introduced to estimate the spatial random deformation at each point of the arch dam, which is caused by the discontinuity of the materials, pouring and damage to the dam, etc. Specifically, the estimation precision of the ordinary kriging interpolation method highly depends on the rationality of the spatial variogram model of the regionalization variable. To ensure that the estimation of the ordinary kriging interpolation method is as precise as possible, the calculation results of Monte Carlo stochastic FEM are applied to construct amounts of random deformation fields by accounting for the uncertainty of the numerical simulation, so that the appropriate variogram models can be determined. On the basis of the basic deformation and spatial variance deformation at each point of the arch dam, the reasonable and reliable arch dam deformation field can be effectively estimated with high accuracy.

The implementation process of the method proposed is introduced in Section 2. The Jinping superhigh arch dam is the engineering application case in Section 3, and the validity of the method is evaluated by comparing the arch dam deformation field estimation results of our method with the estimation results of two other methods. Finally, some conclusions are drawn in Section 4.

2. Method

2.1. FEM Simulation for Basic Deformation of Arch Dams

2.1.1. Hydrostatic-Seasonal-Time (HST) Statistical Model

The HST (hydrostatic-seasonal-time) model has been widely used to describe the three components of the arch dam deformation: the hydraulic component ${\delta }_H$, the temperature component ${\delta }_T$, the time effect component ${\delta }_{\theta }$, and the expressions of the models are as follows:

$$\delta ={\delta }_H+{\delta }_T+{\delta }_{\theta }$$

(1)

$${\delta }_H={\displaystyle \sum _{i=1}^4{a}_i}{H}^i$$

(2)

$${\delta }_T={\displaystyle \sum _{i=0}^1{b}_{1i}}\sin \frac{2\pi it}{365}+{\displaystyle \sum _{i=0}^1{b}_{2i}}\cos \frac{2\pi it}{365}$$

(3)

$${\delta }_{\theta }=c_1\theta +c_2\ln \theta$$

(4)

where $a_i,b_{1i},b_{2i},c_1 \mathrm{and} c_2$ demonstrate regression coefficients; H represents the upstream reservoir water depth; t denotes the number of cumulative days from the initial monitoring day, and $\theta =t/100$.

2.1.2. Parameter Inversion of Arch Dams

The computation parameters of arch dams are commonly affected by complicated and diverse service conditions, and thereby, the parameters usually change at time. To improve the accuracy of the FEM simulation, it is necessary to inverse the computation parameters of arch dams under the actual situation by using the monitoring data.

For an arch dam, the stiffness matrix of the overall structure $\left[K\right]$ consists of the stiffness matrixes of each unit ${\left[K\right]}_{e_j}$. An equation is established between the column vector of the nodal displacement $\left\{\delta \right\}$ and the column vector of the nodal load $\left\{R\right\}$:

$$\left[K\right]\left\{\delta \right\}=\left\{R\right\}$$

(5)

The expression of the stiffness matrix of the overall structure $\left[K\right]$ can be represented as

$$\begin{array}{l}\left[K\right]={\displaystyle \sum _{e_j\in {\Omega }_1}{\left[C\right]}_{e_j}^T{\left[K\right]}_{e_j}{[C]}_{e_j}}+{\displaystyle \sum _{e_j\in {\Omega }_2}{\left[C\right]}_{e_j}^T{\left[K\right]}_{e_j}{[C]}_{e_j}}\\ \begin{array}{cc}& \end{array}+{\displaystyle \sum _{e_j\in {\Omega }_3}{\left[C\right]}_{e_j}^T{\left[K\right]}_{e_j}{[C]}_{e_j}}\end{array}$$

(6)

where ${\Omega }_1$, ${\Omega }_2$, and ${\Omega }_3$ represent the computational domains of the arch dam body, the dam foundation, and the reservoir basin, respectively; ${\left[C\right]}_{e_j}$ is the stiffness transformation matrix of the element $e_j$; and ${\left[K\right]}_{e_j}$ is the stiffness matrix of the unit e_j. The latter can be expressed as

$${\left[K\right]}_{e_j}={\displaystyle {\iiint }_{{\Omega }_j}{\left[B\right]}^T\left[D\right]\left[B\right]d\Omega }$$

(7)

where [D] represents the elastic matrix, which depends on the elasticity modulus E and Poisson ratio $\mu$, and [B] represents the unit of the geometric characteristics.

Consequently, given the influence of water load P at different times, the relationship between the hydraulic component and the comprehensive elastic modulus of the arch dam can be described as follows:

$${\delta }_H=f\left(L,E_r,E_c,{\mu }_r,{\mu }_e,P\right)$$

(8)

where E_r is the comprehensive elastic modulus of the arch dam body; E_c is the comprehensive elastic modulus of the rock dam foundation; ${\mu }_r$ and ${\mu }_e$ represent the Poison ratio of the arch dam body and the foundation, respectively; and L is the size of the arch dam. For a known structure, L can be determined. Because the Poisson ratio slightly impacts on [K] and ${\delta }_H$, the formula can be simplified:

$${\delta }_H=f\left(E_r,E_c,P\right)$$

(9)

The arch dam deformation caused by the temperature load $\Delta T$ is as follows:

$${\delta }_T=J\left(\alpha ,\Delta T\right)$$

(10)

where $\alpha$ is the linear expansion coefficient. If temperature load $\Delta T$ is known, Equation (10) can be simplified as follows:

$${\delta }_T=J\left(\alpha \right)$$

(11)

The creep deformation at each point i of the arch dam can be described as

$${\delta }^i=G\left(P,\dot{P},t\right)$$

(12)

where P is the undetermined parameter; $\dot{P}$ is the changing rate of P; t is the time; and P and t reflect the spatial distribution characteristics of the materials and the temporal characteristics of the materials, respectively. In addition, $P={\left[p_1,p_2,\cdots ,p_l\right]}^T$, l is the number of undetermined parameters.

To describe the creep properties of the materials of arch dams, the generalized Kelvin modal is introduced. The general Kelvin model is illustrated in Figure 1.

In this work, in particular, the comprehensive elastic modulus inversion of arch dams can be summarized as the following fitting optimization problems:

$$Var\left[F\left(E\right)\right]=\min {\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=1}^n{\omega }_{ij}{\left[{\delta }_{Hij}^c-{\delta }_{Hij}^m\right]}^2}}$$

(13)

where $F\left(E\right)$ is the objective function of the inversion analysis of parameter $E_r$ and $E_c$; m is the total number of monitoring points; n is the number of monitoring times; ${\omega }_{ij}$ is the weight of every monitoring value, which is often set to 1; and ${\delta }_{Hij}^c$ and ${\delta }_{Hij}^m$ are the j-th calculation value and the monitoring value of the hydraulic component separated from the monitoring deformation of the monitoring point, respectively.

The linear expansion coefficient inversion of arch dams can be summarized as the following fitting optimization problems:

$$Var\left[F\left(\alpha \right)\right]=\min {\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=1}^n{\omega }_{ij}{\left[{\delta }_{Tij}^c-{\delta }_{Tij}^m\right]}^2}}$$

(14)

where $F\left(\alpha \right)$ is the objective function of inversion analysis of parameter $\alpha$; ${\delta }_{Tij}^c$ and ${\delta }_{Tij}^m$ are the j-th calculation value and the monitoring value of the temporal component separated from the monitoring deformation of the monitoring point, respectively.

The creep parameter inversion of arch dams can be summarized as the following fitting optimization problems:

$$\begin{array}{ll}Var\left[F\left(P\right)\right]& =\min {\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=1}^n{\omega }_{ij}{\left({\delta }_{\theta ij}^c-{\delta }_{\theta ij}^m\right)}^2}}\hfill \\ & =\min {\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=1}^n{\omega }_{ij}{\left({\delta }_{\theta ij}^c-\left({\delta }_{ij}^m-{\delta }_{Hij}^m-{\delta }_{Tij}^m\right)\right)}^2}}P\in D^N\hfill \end{array}$$

(15)

where $F\left(P\right)$ is the objective function of inversion analysis of parameter $P$; ${\delta }_{\theta ij}^c$ and ${\delta }_{\theta ij}^m$ are the j-th calculation value and the monitoring value of the time effect component separated from the monitoring deformation of the monitoring point, respectively; and $D^N$ is the feasible computational domain.

The above optimization problems can be solved by applying the optimization algorithm, and the parameters calculated corresponding to the smallest value of the above optimization objective functions are the inversion results.

2.1.3. Basic Deformation of Arch Dam Calculation

According to the inversion results, the basic deformation of arch dams can be calculated by using an FEM simulation. According to [K] and Equation (5), the hydraulic component of the deformation of arch dam elemental nodes can be simulated. On account of the temperature load, the temperature deformation of each node is equal to the deformation caused by the body force and the surface force. The equilibrium equation used in temperature deformation calculation is

$$\left[K\right]\left\{{\delta }_T\right\}=\left\{R_T\right\}$$

(16)

where $\left\{{\delta }_T\right\}$ is the column vector of the temperature deformation and $\left\{R_T\right\}$ is the column vector of the temperature load of each node.

The column vector of the temperature load of each node $\left\{R_T\right\}$ can be demonstrated as

$$\left\{R_T\right\}=\frac{\alpha E}{1-2\mu }\left({\displaystyle \sum _{i=1}^n{N}_i{T}_i}\right){\left[\begin{array}{cccccc}1& 1& 1& 0& 0& 0\end{array}\right]}^T$$

(17)

where $E$ is the elastic modulus of the arch dam; $\alpha$ is the linear expansion factor; $\mu$ is the Poisson ratio; and $N$ is the element shape function. Consequently, the temperature deformation of each node can be simulated on the basis of the temperature variation field obtained by the FEM simulation. Because the temperature load generates only positive strain, the matrix ${\left[\begin{array}{cccc}1& 1& 1& \begin{array}{cc}\begin{array}{cc}0& 0\end{array}& 0\end{array}\end{array}\right]}^T$ in Equation (17) is used to guarantee only the generation of normal stress.

The FEM simulation of time effect component of dam deformation should be considered from the construction period, and because of the complex loads, it is necessary to introduce the initial strain increment method to calculate the time effect deformation.

When grading the multiple external loads and loading each level of loads, the upper level of the viscous strain of the dam and the rock foundation can be regarded as the initial strain. The nonlinear incremental constitutive relation can be expressed as

$$\left\{\Delta \sigma \right\}=\left[D\right]\left(\left\{\Delta \epsilon \right\}-\left\{\Delta {\epsilon }_V\right\}\right)$$

(18)

where $\left\{\Delta \sigma \right\}$ is the stress increment, $\left\{\Delta \sigma \right\}={\left[\Delta {\sigma }_x,\Delta {\sigma }_y,\Delta {\sigma }_z,\Delta {\sigma }_{xy},\Delta {\sigma }_{yz},\Delta {\sigma }_{zx}\right]}^T$; $\Delta \epsilon$ is strain increment, $\left\{\Delta \epsilon \right\}={\left[\Delta {\epsilon }_x,\Delta {\epsilon }_y,\Delta {\epsilon }_z,\Delta {\epsilon }_{xy},\Delta {\epsilon }_{yz},\Delta {\epsilon }_{zx}\right]}^T$; $\left\{\Delta {\epsilon }_V\right\}$ is viscous strain increment; and $\left[D\right]$ is the elastic matrix.

The basic equilibrium equation of the incremental initial strain method can be described as

$$\left[K\right]\left\{\Delta u\right\}=\left(\left\{\Delta R\right\}-\left\{\Delta R_V\right\}\right)$$

(19)

where $\left\{\Delta R\right\}$ is the column vector of the increment of load of one grade; $\left\{\Delta u\right\}$ is the column vector of the deformation of one element node; and $\left\{\Delta R_V\right\}$ is the column vector of the equivalent node load of the initial strain.

Given the little variation of the Poisson ratio of the dam body and its foundation concrete materials over time, the Poisson ratio of viscous deformation could be replaced by the Poisson ratio of elastic deformation. For $t=t_0+\Delta t$, it is assumed that stress remained virtually unchanged in $\Delta t$. Therefore, the increment of viscous strain of generalized Kelvin modal is

$${\left\{\Delta {\epsilon }_V\right\}}_t={\displaystyle \sum _{i=1}^n\left[1-\exp \left(-\frac{E_K}{{\eta }_K}\Delta t\right)\right]\left\{\frac{1}{E_K}\left[C\right]\left\{\sigma \right\}-{\left\{{\epsilon }_{V,K}\right\}}_{t0}\right\}}$$

(20)

where $E_K$ and ${\eta }_K$ are the elastic shear modulus and the viscosity coefficient of materials of the generalized Kelvin modal, respectively, and [C] is the matrix of the Poisson ratio $\mu$.

According to the inversion results of the main parameters, the basic arch dam deformation m(x) of each point at each time can be calculated by using the FEM simulation. Obviously, the arch dam deformation field can be constructed as a calculation result.

2.2. Ordinary Kriging Interpolation Method for Spatial Variance Deformation of Arch Dams

In order to take the relationship between the points in adjacent area into account, the spatial interpolation method, which factors in the spatial correlations, can be introduced to estimate the spatial variance deformation of arch dams.

However, the deformation magnitude is different in various regions of arch dams, which is caused by the inhomogeneity of the dam materials and the boundary. The spatial variance deformation of arch dams fails to meet the second-order stationarity assumption. Therefore, the ordinary kriging interpolation method cannot be directly used, and the spatial variation characteristic of dam deformation needs to be analyzed further [31].

2.2.1. Spatial Variation Characteristic of Arch Dam Deformation Analysis

The dam material inhomogeneity is discussed below as an example to study the perturbation of deformation caused by material differences. To simplify the problem, the arch dam is approximated as consisting of a series of cantilever beams with variable sections, as shown in Figure 2.

As shown in Figure 2a, the load on the vertical cantilever beam of the arch dam can be expressed as a power series polynomial in head H. According to the differential relationship, the variable-section cantilever beam subjected to the load can be simplified to a uniform-section cantilever beam subjected to two concentrated forces, which is shown in Figure 2b. The deflection curve equation can be demonstrated as

$$\left\{\begin{array}{c}y=\frac{F_1{l}^3}{6EI}\left(3\frac{x^2}{l^2}-\frac{x^3}{l^3}\right)-\frac{F_2{l}^2}{6EI}\left(\frac{3l}{2}-x\right)\left(0\le x\le \frac{l}{2}\right)\\ y=\frac{F_1{l}^3}{6EI}\left(3\frac{x^2}{l^2}-\frac{x^3}{l^3}\right)-\frac{F_2{l}^2}{24EI}\left(3x-\frac{l}{2}\right)\left(\frac{l}{2}\le x\le l\right)\end{array}\right.$$

(21)

where $F_1$ is the left concentrated force of the dam body; $F_2$ is the right concentrated force of the dam body; E denotes the elastic modulus of the dam body; I demonstrates the moment of the inertia of the dam section for the neutral axis; l is the height of the dam body; and x is the vertical coordinate.

In the case of estimating the continuous materials of arch dams, the elastic modulus can be described as E_r, the average value of variable $\beta$ is 0, and the deformation of point $x_1\left(0\le x_1\le \frac{l}{2}\right)$ and point $x_2=x_1+\frac{l}{2}\left(\frac{l}{2}\le x_1\le l\right)$ can be described as

$$m\left(x_1\right)=\frac{F_1{l}^3}{6E_{r}I}\left(3\frac{x_1{}^2}{l^2}-\frac{x_1{}^3}{l^3}\right)-\frac{F_2{l}^2}{6E_{r}I}\left(\frac{3l}{2}-x_1\right)\left(0\le x_1\le \frac{l}{2}\right)$$

(22)

$$m\left(x_2\right)=\frac{F_1{l}^3}{6E_{r}I}\left(3\frac{x_2{}^2}{l^2}-\frac{x_2{}^3}{l^3}\right)-\frac{F_2{l}^2}{24E_{r}I}\left(3x_2-\frac{l}{2}\right)\left(\frac{l}{2}\le x_2\le l\right)$$

(23)

In the case of estimating the discontinuous materials of arch dams, the elastic modulus can be described as ${E^{\prime }}_r=\left(1\pm \beta \right)E_r$, the average value of variable $\beta$ is 0, and the deformation of point $x_1\left(0\le x_1\le \frac{l}{2}\right)$ and point $x_2\left(\frac{l}{2}\le x_2\le l\right)$ can be demonstrated as

$$Z\left(x_1\right)=\frac{F_1{l}^3}{6{E^{\prime }}_{r}I}\left(3\frac{x_1{}^2}{l^2}-\frac{x_1{}^3}{l^3}\right)-\frac{F_2{l}^2}{6{E^{\prime }}_{r}I}\left(\frac{3l}{2}-x_1\right)\left(0\le x_1\le \frac{l}{2}\right)$$

(24)

$$Z\left(x_2\right)=\frac{F_1{l}^3}{6{E^{\prime }}_{r}I}\left(3\frac{x_2{}^2}{l^2}-\frac{x_2{}^3}{l^3}\right)-\frac{F_2{l}^2}{24{E^{\prime }}_{r}I}\left(3x_2-\frac{l}{2}\right)\left(\frac{l}{2}\le x_2\le l\right)$$

(25)

If $x_2=x_1+\frac{l}{2}$, then

$$\frac{Z\left(x_1\right)-m\left(x_1\right)}{m\left(x_1\right)}=\frac{r\left(x_1\right)}{m\left(x_1\right)}=\frac{\beta }{1\mp \beta }$$

(26)

$$\frac{Z\left(x_2\right)-m\left(x_2\right)}{m\left(x_2\right)}=\frac{r\left(x_2\right)}{m\left(x_2\right)}=\frac{\beta }{1\mp \beta }$$

(27)

If $\beta \left(x\right)=r\left(x\right)/m\left(x\right)$, according to the above analysis results, then it can be recognized that the variable $\beta \left(x\right)$ meets the requirement of the second-order stationarity assumption, which means that the distribution laws of two random variables with k-dimensional vectors $\left\{\beta \left(x_1\right),\beta \left(x_2\right),\cdots ,\beta \left(x_k\right)\right\}$ and $\left\{\beta \left(x_1+h\right),\beta \left(x_2+h\right),\cdots ,\beta \left(x_k+h\right)\right\}$ are the same.

In addition, the increment $\beta \left(x\right)-\beta \left(x+h\right)$ of the variable $\beta \left(x\right)$ meets the following two conditions:

(1)

In the whole study domain, the mathematical expectation of increment $\beta \left(x\right)-\beta \left(x+h\right)$ of random function $\beta \left(x\right)$ is 0:

$$E\left[\beta \left(x\right)-\beta \left(x+h\right)\right]=0\hspace{1em}\forall x,\forall h$$

(28)

(2)

The variance function of increment $\beta \left(x\right)-\beta \left(x+h\right)$ of all the vectors exists and is stable:

$$\begin{array}{ll}Var\left[\beta \left(x\right)-\beta \left(x+h\right)\right]& =E{\left[\beta \left(x\right)-\beta \left(x+h\right)\right]}^2\hfill \\ & =2\gamma \left(x,h\right)\hfill \\ & =2\gamma \left(h\right)\hspace{1em}\forall x,\forall h\hfill \end{array}$$

(29)

where $\gamma \left(h\right)$ is variation function of $\beta \left(x\right)$.

The estimation value of $\gamma \left(h\right)$ is ${\gamma }^{*}\left(h\right)$, which can be expressed as

$${\gamma }^{*}\left(h\right)=\frac{1}{2N\left(h\right)}{\displaystyle \sum _{i=1}^{N\left(h\right)}{\left[\beta \left(x_i\right)-\beta \left(x_i+h\right)\right]}^2}$$

(30)

where $N\left(h\right)$ represents the number of the experimental data pair, which is separated by vector h.

Therefore, $\beta \left(x\right)$ meets the requirement of the intrinsic hypothesis. The ordinary kriging interpolation method can be introduced to estimate the $\beta \left(x\right)$ of each known point of arch dams. After estimating the ${\beta }^{*}\left(x\right)$ of the unmonitored points according to the $\beta \left(x\right)$ of the known points and calculating the basic deformation m(x), the random spatial deformation $r^{*}\left(x\right)$ of the unmonitored points can be described as follows:

$$r^{*}\left(x\right)={\beta }^{*}\left(x\right)\cdot m\left(x\right)$$

(31)

Consequently, it is necessary to estimate $\beta \left(x\right)$ by using the ordinary kriging interpolation method.

2.2.2. Estimation of Variable $\beta \left(x\right)$

In order to estimate the $\beta \left(x\right)$ of each point of arch dams, the Monte Carlo stochastic FEM is used to obtain several groups of samples of $\beta \left(x\right)$,which can be applied in calculating the theoretical variogram model of $\beta \left(x\right)$.

(1)

The calculation of the samples of variable $\beta \left(x\right)$ is based on the Monte Carlo stochastic FEM [32].

The reasons behind the results of the spatial variance deformation of arch dams mainly include the discontinuous materials of the dam body and the dam foundation, the initial strain caused by pouring, the arch closure, the temperature control, the random microcracks of arch dam and uneven uplift pressure, and so on. The above reasons can be regarded as variables that follow definite distribution law and can be used to generate an arch deformation field under different situations on the basis of using the Monte Carlo stochastic FEM. First, the interval can be divided into N nonoverlapping subintervals by using the Latin hypercube sampling method. On this basis, in interval [0, 1], N random numbers can be generated after separately completing equal probability sampling in every subinterval. The random number in the n-th subinterval is

$$U_n=\frac{U}{N}+\frac{n-1}{N}\left(n=1,2,\cdots ,N\right)$$

(32)

where U is a uniform random number in interval [0, 1].

After the above computational analysis, the random numbers in each subinterval are transformed into sample values by using the transformation method. In the sampling simulation, N values of each input sample x_k are associated with a random sequence of integers representing the integers 1, 2, …, N, and each input sample is then substituted into an FEM simulation to obtain a number of the output samples of the arch dam deformation field corresponding to each input sample.

On the basis of one of the sets of deformation field calculation results, the deformation of each point is sampled, and the corresponding basic deformation of the point is calculated. Next, a group sample of $\beta \left(x\right)$ can be calculated as follows:

$$\left\{\begin{array}{c}\beta \left(x_1\right)=\frac{r\left(x_1\right)}{m\left(x_1\right)}=\frac{Z^{\prime }\left(x_1\right)-m\left(x_1\right)}{m\left(x_1\right)}\\ \beta \left(x_2\right)=\frac{r\left(x_2\right)}{m\left(x_2\right)}=\frac{Z^{\prime }\left(x_2\right)-m\left(x_2\right)}{m\left(x_2\right)}\\ \begin{array}{c}⋮\hspace{1em}\hspace{1em}=\hspace{1em}\hspace{1em}⋮\hspace{1em}=\hspace{1em}\hspace{1em}⋮\\ \beta \left(x_n\right)=\frac{r\left(x_n\right)}{m\left(x_n\right)}=\frac{Z^{\prime }\left(x_n\right)-m\left(x_n\right)}{m\left(x_n\right)}\end{array}\end{array}\right.$$

(33)

Repeating the above steps will eventually result in several sets of samples of $\beta \left(x\right)$.

(2)

The variation function calculation of variable $\beta \left(x\right)$ when the spatial point varies on the x-axis, half of the variance of the difference between the value of the variable $\beta \left(x\right)$ at x and x + h is defined as the variance function of $\beta \left(x\right)$ in the direction of the x-axis. This is denoted as $\gamma \left(x,h\right)$, and the expression is as follow:

$$\begin{array}{ll}\gamma \left(x,h\right)& =\frac{1}{2}Var\left[\beta \left(x\right)-\beta \left(x+h\right)\right]\hfill \\ & =\frac{1}{2}E{\left[\beta \left(x\right)-\beta \left(x+h\right)\right]}^2-\frac{1}{2}{\left\{E\left[\beta \left(x\right)\right]-E\left[\beta \left(x+h\right)\right]\right\}}^2\hfill \end{array}$$

(34)

Under the hypothesis of the second-order stationarity, the following equation stands:

$$E\left[\beta \left(x+h\right)\right]=E\left[\beta \left(x\right)\right]\forall h$$

(35)

When the variance function is independent of x but only related to the distance h between two sample points, the following equation stands:

$$\gamma \left(h\right)=\frac{1}{2}E{\left[\beta \left(x\right)-\beta \left(x+h\right)\right]}^2$$

(36)

On the basis of the solved distance of samples and semivariance, the variation function of the variation curve can be fitted, which reflects the relationship between the value of the variation function ${\gamma }^{*}\left(h\right)$ and the certain lag distance h. In Figure 3, $C_1$ is the ‘nugget’ of $\gamma \left(h\right)$, which means the monitoring errors and the short distance are unknown or constitute random spatial variation. $C_0$ is the block gold effect, which indicates the fluctuation of the variables between two points when h is very small. a is the variation range. When $h\le a$, the correlation of the variable values between any two points decreases as h becomes larger. In addition, when $h>a$, the correlation no longer appears. a reflects the range of influence. $C\left(C=C_0+C_1\right)$ is the total abutment value, which reflects the intensity of the variation of a regionalized variable in the study area. When $h\to \infty$, the value of $\gamma \left(h\right)$ can be calculated as follows:

$$\gamma \left(\infty \right)=C_0+C_1=C$$

(37)

Commonly used theoretical models include the spherical model, which can be expressed as follows [33]:

$$\gamma \left(h\right)=\left\{\begin{array}{c}{C}_0\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} h=0\\ C_0+C_1\left(\frac{3}{2}\frac{h}{a}-\frac{1}{2}\frac{h^3}{a^3}\right)\hspace{1em}0<h\le a\\ C_0+C_1\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}h>a\end{array}\right.$$

(38)

In practice, $\beta \left(x\right)$ is more complex and has a certain probability of being anisotropic or containing multiple levels of variability at different scales in the same direction. This anisotropy of regionalized variables can be characterized by using a nested structure, which is a combination of variability that appears to simultaneously act at different distances h and in different directions $\eta$.

The nested structure can be expressed as the sum of multiple variance functions, each of which represents a variability at a particular scale, and the expression of the nested structure is

$$\gamma \left(h\right)={\gamma }_0\left(h\right)+{\gamma }_1\left(h\right)+\cdots +{\gamma }_n\left(h\right)+\cdots$$

(39)

Consequently, the best linear and unbiased estimate of $\beta \left(x\right)$ at a set of unmonitored points can be performed by using a known set of data $\left\{\beta \left(x_i\right),i=1,2,\cdots ,n\right\}$ that is based on the ordinary kriging interpolation algorithm.

Assume that the estimation of any block V is ${\beta }_k^{*}$. This estimation is a linear combination of n information values ${\beta }_i$ in the estimated neighborhood, which is demonstrated as follows:

$${\beta }_k^{*}={\displaystyle \sum _{i=1}^n{\lambda }_i{\beta }_i}$$

(40)

In order to obtain an unbiased estimation under the second-order stationarity assumption, it is required that ${\displaystyle \sum _{i=1}^n{\lambda }_i=1}$.

The equation group of the ordinary kriging interpolation method is obtained as follows:

$$\left\{\begin{array}{c}{\displaystyle \sum _{j=1}^n{\lambda }_j\overline{\gamma }\left(v_i,v_j\right)+\mu =\overline{\gamma }\left(v_i,v_0\right)}\\ {\displaystyle \sum _{j=1}^n{\lambda }_j=1}\end{array}\right.\left(i=1,2,\cdots ,n\right)$$

(41)

where $\overline{\gamma }\left(v_i,v_j\right)$ is the semivariogram of $\beta \left(x_i\right)$ and $\beta \left(x_j\right)$; $\mu$ is the Lagrangian multiplier with the minimum variance; $\overline{\gamma }\left(v_i,v_0\right)$ is the semivariogram of $\beta \left(x_i\right)$ and $\beta \left(x_0\right)$; and $x_0$ is the location of the point that needs to be estimated.

The matrix form of the above equation group is written as

$$\left[K^{\prime }\right]\left[{\lambda }^{\prime }\right]=\left[{M^{\prime }}_2\right]$$

(42)

$${\sigma }_K^2={\left[{\lambda }^{\prime }\right]}^T\left[M_2^1\right]-\gamma \left(v_0,v_0\right)$$

(43)

where

$$\left[K^{\prime }\right]=\left[\begin{array}{ccccc}\overline{\gamma }\left(v_1,v_1\right)& \overline{\gamma }\left(v_1,v_2\right)& \cdots & \overline{\gamma }\left(v_1,v_n\right)& 1\\ \overline{\gamma }\left(v_2,v_1\right)& \overline{\gamma }\left(v_2,v_2\right)& \cdots & \overline{\gamma }\left(v_2,v_n\right)& 1\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ \overline{\gamma }\left(v_n,v_1\right)& \overline{\gamma }\left(v_n,v_2\right)& \cdots & \overline{\gamma }\left(v_n,v_n\right)& 1\\ 1& 1& \cdots & 1& 0\end{array}\right]$$

(44)

$$\left[{\lambda }^{\prime }\right]=\left[\begin{array}{c}{\lambda }_1\\ {\lambda }_2\\ ⋮\\ {\lambda }_n\\ \mu \end{array}\right]$$

(45)

$$\left[{M^{\prime }}_2\right]=\left[\begin{array}{c}\overline{\gamma }\left(v_1,v_0\right)\\ \overline{\gamma }\left(v_2,v_0\right)\\ ⋮\\ \overline{\gamma }\left(v_n,v_0\right)\\ 1\end{array}\right]$$

(46)

According to Equations (42) and (43), the weight coefficients can be calculated, then the estimation value ${\beta }_k^{*}$ of any block segment V can be calculated by using Equation (40), and the spatial variance deformation at any point can be estimated by using Equation (31). After calculating the value of the basic deformation m_k at the corresponding point, the estimation value of the dam deformation at each point can finally be obtained.

2.3. Process of the Arch Dam Deformation Field Estimation

Figure 4 demonstrates the process of the arch dam deformation field estimation on the basis of using the proposed method. The arch dam deformation field estimation steps of applying the ordinary kriging interpolation method, combined with FEM, can be implemented as follows:

(1)

Preprocessing the original deformation monitoring data of arch dam and obtaining valid deformation monitoring data by using outlier detection and missing data processing.

(2)

Inverse analyzing the arch dam materials and mechanic parameters by using valid data.

(3)

Calculating arch dam basic deformation field on the basis of using the FEM.

(4)

Simulating deformation field samples by factoring in the discontinuity of materials, pouring, and damage to the dam on the basis of using the Monte Carlo stochastic FEM.

(5)

Constructing the variation model of variable $\beta \left(x\right)$ and estimating the value of $\beta \left(x\right)$ at each point of arch dam by using the ordinary kriging interpolation method.

(6)

Estimating spatial variance deformation at each point of the arch dam.

(7)

Adding the basic deformation $m\left(x\right)$ and spatial variance deformation $r\left(x\right)$ at each point of arch dam, then finishing the arch dam deformation field estimation.

3. Case Study

3.1. Brief Introduction to a Superhigh Arch Dam and Settlement of Its Deformation Monitoring Points

In recent decades, China has built a number of arch dams in the southwest of the country. Among them, the height of the Jinping superhigh arch dam is 305 m. It is currently the highest double-curvature arch dam in the world. The location of the Jinping arch dam is shown in the topographic map in Figure 5b.

The automatic dam structural health monitoring (SHM) system of the Jinping arch dam was completed at the end of 2015. Dam deformation is an essential and intuitive monitoring item to reflect dam structural service characteristics, especially for super-high arch dams. A plumb line (PL) system is the most popular dam deformation monitoring technique for arch dams. The automatic structural monitoring system of the Jinping super-high arch dam was completed and began to operate at the end of 2015. Dam deformation is a key monitoring item that reflects the dam structural operating characteristics, especially for super-high arch dams.

At the elevations of 1601 m, 1664 m, 1730 m, 1778 m, and 1829 m, 34 deformation monitoring points of PL systems were arranged, which are used to monitor the deformation of the dam body. Figure 6 shows the layout of the PL systems of the Jinping arch dams.

3.2. Numerical Case Study

3.2.1. The Design of Numerical Case

Because of the limited quality of the monitoring points of the arch dam, the spatially continuous deformation field of the arch dam cannot be directly monitored. To verify the estimation accuracy of the method proposed, it is necessary to generate several spatially continuous artificial deformation fields as sample fields. In this case, the sample deformation field is constructed by using the Monte Carlo stochastic FEM and the deterministic FEM by including the spatial inhomogeneity of the materials. Consequently, a numerical case is introduced in this study, and the sample deformation field generation process can be described as follows:

(1)

Simulating several deformation field samples of the arch dam by using the Monte Carlo stochastic FEM according to the statistical characteristics of the material parameter of the arch dam.

(2)

Calculating the variance distribution characteristic of the deformation field samples.

(3)

Generating more sample deformation fields of arch dams on the basis of the variance distribution characteristic.

(4)

Simulating the arch dam deformation field under the normal loads by using the deterministic FEM.

(5)

Superposing the deformation field generated in (3) and the deformation field simulated in (4) and then generating sample deformation fields.

The 3D finite element model of the Jinping arch dam is selected as the generalized model in this case, which is shown in Figure 7. The dam body and the surrounded mountain are divided into 38,537 elements. The average size of the element is 9 m, and the cell types consist of the linear hexahedral cell and trigonal cells. The model can be divided into three main parts: the dam body, the dam foundation, and the nearby mountains. The model was calculated to include 800 m from the axis of the arch dam on the left and right banks, 800 m upstream and downstream of the origin of the dam crest, 600 m below the footprint, and from the rocky slope above the crest to the height of the natural boundary. The hills to the left and right of the dam crest extend 2.5 times the upstream and downstream quadrants, respectively. The rock conditions and mound shapes on both sides of the dam were simulated after excavation.

The fixed constraints were imposed at the bottom of the model, and normal constraints were imposed around the perimeter. The water load was applied to the upstream dam face and the upstream reservoir area considered in the model, and the dam body and foundation adopted the linear elasticity principal structure. The actual parameters for the concrete material and the inversion parameters used in the case study are listed in

Table 1. Real parameters of arch dam materials and the inversion parameters.

Zone	$\mathit{\gamma }$ (KN/m3)	$\mathit{\mu }$	$\mathit{\alpha }$	${\mathit{E}}^{\mathit{i}\mathit{n}\mathit{v}\mathit{e}\mathit{r}\mathit{s}\mathit{e}}\left(\mathbf{GPa}\right)$	${\mathit{E}}_{\mathit{M}}^{\mathit{i}\mathit{n}\mathit{v}\mathit{e}\mathit{r}\mathit{s}\mathit{e}}\left(\mathbf{GPa}\right)$	${\mathit{\eta }}_{\mathit{M}}^{\mathit{i}\mathit{n}\mathit{v}\mathit{e}\mathit{r}\mathit{s}\mathit{e}}(\mathbf{GPa}.\mathbf{S})$
Dam body	2400	0.17	1.0 × 10−5	30	75.5	1.36 × 107
Dam foundation	2700	0.20		25

.

The locations of the monitoring points in PL systems were selected as the locations of the sampling points. According to the results of the finite element numerical simulation, the basic deformation field of the arch dam can be calculated so that the basic deformation values of the monitoring points can be obtained: $m\left(x_1\right),m\left(x_2\right),\cdots ,m\left(x_l\right)$. On this basis, $Z^{\prime }\left(x_1\right),Z^{\prime }\left(x_2\right),\cdots ,Z^{\prime }\left(x_l\right)$ can be calculated by using the Monte Carlo stochastic FEM. Consequently, according to Equation (33), $\beta \left(x_1\right),\beta \left(x_2\right),\cdots ,\beta \left(x_l\right)$ can be calculated. For one of the generated sample deformation fields, the distribution of the sample distance and semivariance of 100 samples of $\beta \left(x\right)$ are shown in Figure 8.

After fitting the variance function five times, the function model can be obtained as

$$\gamma \left(h\right)=22.31+\mathrm{99,080.13}sph\left(347.04\right)$$

(47)

where the anisotropy ratio is $r=1.06$ and the geometric anisotropy angle is $\phi ={10.95}^{\circ }$.

Once the variation function model has been established, the weight coefficients can be calculated according to the kriging equation in Equation (41), and $\beta \left(x_i\right)\left(i=1,2,\cdots ,l\right)$ at the unmonitored points are estimated from the known values of ${\beta }^{*}\left(x_p\right)\left(p=1,2,\cdots ,P\right)$ at the monitoring points. After obtaining the basic deformation values $m\left(x_p\right)\left(p=1,2,\cdots ,P\right)$ at the unmonitored points from the FEM simulation results, the estimated values of $r^{*}\left(x_p\right)\left(p=1,2,\cdots ,P\right)$ at the unmonitored points can be solved according to Equation (31). The estimated values of the deformation at the unmonitored points can be estimated on the basis of the estimated values of $m\left(x_p\right)\left(p=1,2,\cdots ,P\right)$ and $r^{*}\left(x_p\right)\left(p=1,2,\cdots ,P\right)$ at the unmonitored points and according to the deformation of other known points.

3.2.2. Analysis of the Estimation Results of Arch Dam Sample Deformation Fields

The simulation of the FEM selects the mean value of the parameter samples as parameters, and the simulation results of arch dam deformation fields are demonstrated in Figure 9. The deformation fields of the arch dam are estimated on the basis of 10 groups of the arch dam radial and tangential sample deformation fields and by using the ordinary kriging method, the universal kriging interpolation, and the proposed method. The radial deformation fields and tangential deformation fields of one group are illustrated in Figure 10 and Figure 11, respectively. It can be seen that simulation results of arch dam deformation fields based on FEM apparently vary from the radial and tangential sample deformation fields illustrated in Figure 10a and Figure 11a. The rest of the groups have similar results.

It can be observed that the variation of the deformation field simulated by using the deterministic FEM is pretty uniform, but the variation of the sample deformation field is fluctuant because of the spatial inhomogeneity of the materials. However, the deformation fields calculated by using the deterministic FEM and the sample deformation field have fairly similar basic spatial distribution characteristics. Therefore, given the spatial inhomogeneity of the arch dam materials, the deformation field of the deterministic FEM cannot completely coincide with the sample deformation fields.

3.2.3. Comparison Analysis of Numerical Deformation Field Estimation Results

As can be observed from Figure 10 and Figure 11, the deformation fields estimated by the proposed method are extremely similar to those of the sample deformation field, and the variation of the deformation fields estimated according to the proposed method is essentially the same as that of the sample deformation fields. Although the ordinary kriging interpolation method has relatively stable estimation accuracy at points near the sampling points, it does not work well when it is used to estimate the deformation at points far from the sampling points. Owing to the limited number of sampling points, the trend surface of the universal kriging interpolation does not describe the actual trend surface of the sample deformation field. As a result, the deformation field estimates from the ordinary kriging interpolation method significantly differ from the sample deformation fields and cannot describe the actual deformation characteristics of the arch dam.

In order to quantify the estimation precisions of the above methods, the mean absolute error (MAE) is introduced to analyze the errors at all points of the sample deformation fields $j\left(j=1,2,\cdots ,10\right)$. The MAE of one group of the deformation estimation results $MAE\left(A_j\right)$ can be expressed as

$$MAE\left(A_j\right)=\frac{1}{n_p}{\displaystyle {\sum }_{i=1}^{n_p}\left|Z_j\left(p_i\right)-{Z^{\prime }}_j\left(p_i\right)\right|}j=1,2,\cdots ,10$$

(48)

where $Z_j\left(p_i\right)$ represents the deformation at $p_i$ of sample deformation field j; ${Z^{\prime }}_j\left(p_i\right)$ indicates the estimation result of the arch dam deformation at node $p_i$ of sample deformation field j; and $n_p$ represents the number of points, $n_p=31941$ in this case. Therefore, a comprehensive index that factors in the 10 groups of deformation field estimation results can be defined as

$$AMAE=\frac{1}{10}{\displaystyle \sum _{j=1}^{10}MAE(A_j)}$$

(49)

where AMAE means the average value of MAE(A_j)(j = 1,2,…,10).

The MAE(A_j) and AMAE of the above methods are illustrated in Figure 12. The comparison results show that the proposed method achieves a higher estimation precision, and the AMAE values of the 10 groups of radial deformation fields and those of the tangential deformation fields are 1.977 mm and 1.654 mm, respectively. Meanwhile, the MAE(A_j) of the estimation results of radial and tangential deformation fields are much lower than the estimation results of the ordinary kriging method and the universal kriging interpolation method.

3.3. Accuracy Verification Based on the Actual Monitoring Data

The monitoring deformation values of all the monitoring points of 3 typical days in 2018 were selected to verify the estimation accuracy and robustness of the proposed method. Here, 3 typical days consist of a day with a low water level and a high temperature (12 June 2018, 1801.19 m, 15.60 °C), a day with a high water level and a high temperature (4 August 2018, 1879.70 m, 17.60 °C), and a day with a high water level and a low temperature (16 December 2018, 1878.83 m, 5.4 °C). Given that the actual monitoring data of the monitoring points are discrete in space, this provides an idea to introduce the cross-validation method to compare the estimation accuracy of the proposed method, the ordinary kriging method, and the universal kriging algorithm. The basic ideas are as follows:

(1)

Individually remove the monitoring values of all the monitoring points $Z\left(x_i\right)\left(i=1,2,\cdots ,n\right)$ in order.

(2)

Estimate the value of $Z^{*}\left(x_i\right)\left(i=1,2,\cdots ,n\right)$ at each removed monitoring point on the basis of the remaining monitoring values by using the above three methods (n is the number of monitoring points).

(3)

Repeat the above operations for each monitoring point by using different methods.

(4)

Calculate the estimation error between the original data and the estimated value of each monitoring point.

The RMSE can be calculated according to Equation (50) to verify the precision of the interpolation. At the same time, the correlation degree can be assessed by using the coefficient R². The calculation function is shown in Equation (51). The stronger the correlation between the two values is, the more accurate the value of R² is.

$$RMSE=\sqrt{\frac{{\displaystyle \sum _{i=1}^n{\left(Z\left(x_i\right)-Z^{*}\left(x_i\right)\right)}^2}}{n}}$$

(50)

$$R^2=1-\frac{S_{res}}{S_{total}}$$

(51)

In Equation (51), $S_{res}$ and $S_{total}$ represent the sum of the squared residuals and the total sum of squares, respectively. The calculation equations are as follows:

$$S_{res}={\displaystyle \sum _{i=1}^n{\left(Z\left(x_i\right)-Z^{*}\left(x_i\right)\right)}^2}$$

(52)

$$S_{total}={\displaystyle \sum _{i=1}^n{\left(Z\left(x_i\right)-\overline{Z}\right)}^2}$$

(53)

where $\overline{Z}$ is the average value of the monitoring deformation values of the sample points.

The cross-validation results of the three methods are shown in

Table 2. Cross-validation results of arch dam radial deformation fields.

NO.	Time	The Ordinary Kriging Interpolation		The Universal Kriging Interpolation		The Proposed Method
		RMSE (mm)	R2	RMSE (mm)	R2	RMSE (mm)	R2
1	12 June 2018	6.302	0.773	5.321	0.826	1.843	0.918
2	4 August 2018	5.914	0.764	4.997	0.742	1.352	0.901
3	16 December 2018	6.014	0.811	5.738	0.826	1.547	0.925

and

Table 3. Cross-validation results of arch dam tangential deformation fields.

NO.	Time	The Ordinary Kriging Interpolation		The Universal Kriging Interpolation		The Proposed Method
		RMSE (mm)	R2	RMSE (mm)	R2	RMSE (mm)	R2
1	12 June 2018	3.894	0.671	2.595	0.715	1.029	0.879
2	4 August 2018	4.025	0.749	2.706	0.763	1.133	0.853
3	16 December 2018	3.998	0.792	2.615	0.807	1.154	0.898

. The following conclusions can be drawn:

(1)

The values of RMSE of the proposed method at different times are less than the value of other methods, and the R² of the proposed method is higher than the R² of other methods.

(2)

Regardless of the time, the deformation field estimation results of the proposed method are more accurate than the estimation results of other methods. It is more useful for the analysis of the whole deformation trend of the arch dam.

(3)

According to the results of the methods, R² is basically consistent at different times, which illustrates that the proposed method has high stability and that the estimation results are accurate.

It can be seen that the RMSE of the proposed method at any time is smaller than the other two methods, and the R² is larger than that of the other two methods. The MAE values of the radial deformation and that of the tangential deformation of the proposed method are 1.582 mm and 0.963 mm, respectively. They are both much lower than those of the ordinary kriging method (6.161 mm and 3.635 mm) and those of the universal kriging (5.162 mm and 2.401 mm) interpolation method. For monitoring points with other sampling points in the neighborhood (e.g., PL11-2 and PL19-2), the precision of the proposed method and that of the universal kriging interpolation method are higher than that of the ordinary kriging interpolation method. This is mostly because the trend surface of the universal kriging interpolation method is to a certain extent able to reflect the basic deformation characteristic of the arch dam. On the contrary, for the monitoring points where there are no monitoring points near them (e.g., PL1-1, PL23-1), the proposed method performs better than the other two methods, which primarily results from the inclusion of the physical significance of the dam deformation in the proposed method.

Figure 13 and Figure 14 show the estimation results of radial and tangential arch dam deformation fields on 12 June 2018 based on the above methods, respectively. It can be found that the deformation fields estimated on the basis of the ordinary kriging method and the universal kriging interpolation method do not agree with those simulated by using the FEM in some regions of the arch dam. In summary, the practical application shows that the proposed method accurately estimates the arch dam deformation fields.

However, the estimation errors of the ordinary kriging method and universal kriging interpolation method are small at the location close to the monitoring points, and the errors are large at the location away from the monitoring points. The error range of the radial deformation of the universal kriging method and the ordinary kriging interpolation method are both pretty large. Because the proposed method factors in the advantages of the FEM with the ordinary kriging interpolation method, the error range of the radial deformation estimation results is from 0.755 mm to 3.458 mm. The error range of the tangential deformation estimation results is from 0.560 mm to 1.531 mm.

4. Conclusions

Current data-based deformation field estimation methods and the FEM encounter difficulty in trying to accurately estimate arch dam deformation field. This study proposes combining the ordinary kriging interpolation method with the FEM to create the arch dam deformation field estimation method, which is supposed to be successfully applied to the determination of arch dam spatial deformation characteristics. The performance of the proposed method is tested through a numerical case study and the engineering application of the Jinping superhigh arch dam. Some conclusions can be drawn, as follows:

(1)

The estimation result of the proposed method in this paper is based on the FEM simulation result, which is rich with data and has low precision. After merging the ordinary kriging interpolation results, the precision of the arch dam deformation estimation has been effectively improved. Compared with the kriging-based method, this method is more stable and has higher estimation accuracy.

(2)

Sufficient and valid monitoring data help to ensure the validity of the arch deformation field estimation results of the proposed method.

(3)

The effect of inhomogeneous material properties on arch dam deformation is reflected by 10 sample deformation fields; the uneven dam deformation caused by material inhomogeneity leads to errors in the arch dam deformation field simulation results of the FEM. The errors of the ordinary kriging method and the universal kriging interpolation method come from their failing to sufficiently account for and inaccurately calculating the physical insignificance of the arch dam deformation field.

The method proposed in this paper is able to also estimate the temperature or stress field of the arch dam on the basis of using valid monitoring data. The estimation effects of the above fields based on the proposed method also need to be further verified.