Parameter identification of the dynamic Winkler soil–structure interaction model using a hybrid unscented Kalman filter

Abstract

Soil–structure interaction is a key aspect to take into account when simulating the response of civil engineering structures subjected to dynamic actions. To this end, and due to its simplicity and ease of implementation, the dynamic Winkler model has been widely used in practical engineering applications. In this model, soil–structure interaction is simulated by means of spring–damper elements. A crucial point to guarantee the adequate performance of the approach is to accurately estimate the constitutive parameters of these elements. To this aim, this article proposes the application of a recently developed parameter identification method to address such problem. In essence, the parameter identification problem is transformed into an optimization problem, so that the parameters of the dynamic Winkler model are estimated by minimizing the relative differences between the numerical and experimental modal properties of the overall soil–structure system. A recent and efficient hybrid algorithm, based on the combination of the unscented Kalman filter and multi-objective harmony search algorithms, is satisfactorily implemented to solve the optimization problem. The performance of this proposal is then validated via its implementation in a real case-study involving an integral footbridge.

Keywords

dynamic Winkler model finite element model updating hybrid local–global optimization algorithms operational modal analysis parameter identification soil–structure interaction

Introduction

Soil–structure interaction (SSI) is an important aspect to consider in civil engineering structures in order to mimic accurately their behaviour via numerical models (Tabatabaiefar et al., 2016). One of the most widely used SSI models is the static Winkler model (Hirai, 2011), where the SSI phenomenon is simulated via the implementation of spring elements, whose behaviour is characterized by a stiffness coefficient. When dynamic actions come into play, it is necessary to implement a dynamic Winkler model (Allotey and El Naggar, 2008), where the SSI is simulated via spring–dampers elements, whose behaviour is characterized by both stiffness and damping coefficients. However, both the static and dynamic Winkler models present a clear limitation for their practical engineering application: the need to accurately estimate the value of the above constitutive parameters (stiffness/damping coefficients) that resemble the actual SSI. In fact, this parameter identification is not a straightforward issue since it involves testing jointly both systems: the soil and the structure.

In order to address such parameter identification problem, two sets of methods have been traditionally employed. In the first type of methods, the SSI phenomenon is characterized via the use of numerical models whose constitutive laws are previously determined via either laboratory (Drnevich et al., 1978) or field tests (Badsar et al., 2010; Campanella and Stewart, 1992; Cheng and Leong, 2018; Karl et al., 2006; Lai and Rix, 1998; Lee et al., 2016). Among the different methods, the so-called mixed finite element (FE)–boundary element (BE) method has been widely implemented (Alamo et al., 2016; Yazdchi et al., 1999). In this approach, the FE method, used to simulate the structure, and the BE method, employed to model the soil, are coupled by imposing compatibility and equilibrium equations at the soil–structure interphase. Although these methods constitute a valid attempt to simulate numerically the SSI phenomenon, they present three main limitations that discourage their use for practical engineering applications, namely: (1) the constitutive laws of the soil are estimated without taking into account the SSI phenomenon; (2) the results rely strongly on such constitutive laws; and (3) they are somehow complex to be implemented by practitioner engineers.

In order to overcome these limitations, a second type of methods emerged that also combined numerical and experimental techniques, although with a different perspective: the parameters that characterize the SSI phenomenon are estimated via the implementation of some kind of parameter identification algorithm. Among these identification algorithms, the so-called estimators have been extensively employed. Within these estimators, the maximum likelihood method has been widely adopted, due to its ease of implementation and good balance between computational cost and accuracy. This method transforms the parameter identification problem into an optimization problem, in which the parameters to be estimated are determined via the minimization of an objective function, defined in terms of the relative differences between some numerical and experimental data that reflect the behaviour of the soil–structure system. According to the best of the authors’ knowledge, the first examples of this approach were reported by Maier and Gioda (1982) and Arai et al. (1984). In both studies, a gradient-based algorithm was considered as optimization method. Subsequently, Gens et al. (1996) used the Levenberg–Marquardt algorithm to solve the parameter identification problem during a tunnel excavation work. The main limitation of these methods is the local character of the optimization algorithm: local optimizers present a high risk of convergence to a local minimum and are highly dependent on the initial value of the design parameters.

Subsequently, in order to overcome these limitations, global optimization algorithms were considered to solve the parameter identification problem in geotechnical applications (Ledesma et al., 1996a, 1996b). Due to their good convergence and effectiveness to solve nonlinear problems, genetic algorithms have been widely used to this end (Srinivas and Deb, 1994). However, the use of genetic algorithms presents an important drawback: the high simulation time required to perform the optimization process for complex structures. For this reason, during the last 25 years, great efforts have been devoted to propose and implement alternative global optimizers that reduce the simulation time without comprising the accuracy of the solution.

However, when it becomes mandatory to take into account the stochastic character of either the soil or the structural system, another kind of identification methods is required, which involve estimators that filter the uncertainty associated with the stochastic nature of the constitutive parameters. Estimators based on the so-called Kalman filter (KF) have clearly imposed to this end (Kalman, 1960). The KF estimates the parameters of a linear system by considering statistical white noise of the experimental data. Thus, the mean and covariance of the parameters of the SSI model are estimated iteratively in two steps (prediction and correction) in order to minimize the estimation error, and, therefore, achieve an unbiased true estimation of the parameters of the model. As the use of the KF was limited to linear systems, a new algorithm – the extended Kalman filter (EKF) – was later developed to address the parameter identification problem for nonlinear systems (Jazwinski, 1970). The EKF algorithm is based on the local linearization of the nonlinear system and the subsequent application of the KF algorithm on the resulting linearized system. However, the performance of this algorithm is reduced as the nonlinear character of the system increases. In such case, a new enhanced version of the KF was proposed: the unscented Kalman filter (UKF). In this method, no linearization of the system is required, so that, instead, Gaussian behaviour of the system is assumed (Julier and Uhlmann, 1997). Thus, the covariance of the estimation process is determined from the propagation of a set of sample points through the nonlinear system. Thus, the UKF algorithm allows computing more accurately the mean and covariance of the estimated parameters than the EKF algorithm (up to the second order of the Taylor series expansion of the nonlinear function and up to the third order for Gaussian inputs) (Nguyen and Nestorović, 2015). However, these algorithms, based on KF, are local optimizers, and therefore may fail when finding the global optimum of functions with several extreme values.

For these reasons, the current trend to overcome all the above-mentioned limitations and reduce computational cost is to solve the parameter identification problem using hybrid algorithms that take advantage of the virtues of both local and global optimizers. For instance, a hybrid UKF-simulated annealing algorithm was proposed by Astroza et al. (2016) that allowed reducing the simulation time when compared to the sole implementation of the simulated annealing algorithm. More recently, Naranjo-Perez et al. (2020) proposed a hybrid UKF–harmony search (HS) algorithm for the FE model updating of complex civil engineering structures. This hybrid algorithm combined the virtues of two individual algorithms as follows: (1) the (local) UKF algorithm, which reduces the uncertainty associated with the parameter identification based on experimental data; and (2) the (global) HS algorithm (Geem et al., 2001), which reduces significantly simulation times when compared to other classical metaheuristic algorithms (Jiménez-Alonso et al., 2017).

In this article, the UKF–HS algorithm is extended to consider a more general multi-objective approach (UKF–MHS), and it is further implemented to solve the parameter identification problem associated with the experimental identification of the dynamic Winkler model parameters. The parameter identification problem is herein formulated based on the FE model updating method under the maximum likelihood approach (Mottershead et al., 2011). In this manner, the parameter identification problem is transformed into a multi-objective optimization problem (Marwala, 2010), with the objective functions defined in terms of the relative differences (residuals) between the numerical and experimental modal properties of the SSI system. In particular, three types of residuals are next considered as follows: (1) residuals based on the natural frequencies; (2) residuals based on the damping ratios; and (3) residuals based on the vibration modes. Therefore, the aim of the resulting optimization problem focuses on minimizing the value of the objective functions via the modification of some pre-selected physical parameters of the structure. The solution of this optimization problem provides a point estimation of the value of such physical parameters. In this process, the numerical modal properties of the soil–structure system are obtained via the FE method, while its experimental modal properties are identified via the signal processing of its response recorded during either a forced vibration (Maia et al., 1997) or an ambient vibration test (Magalhães and Cunha, 2011).

Thus, the following three are the main aspects that characterize the proposed parameter identification method:

It is a model-based identification method, so that the solution of the parameter identification problem involves the solution of an FE model updating problem.

The identification of the parameters defining the dynamic Winkler model is performed indirectly, by inferring them from the analysis of the overall experimental behaviour of the soil–structure system.

A hybrid UKF–MHS algorithm is implemented to efficiently perform the parameter identification process.

Finally, the performance of this parameter identification method has been validated via its implementation in a real case study. Concretely, a real steel–concrete composite integral footbridge has been considered to conduct the parameter identification of the dynamic Winkler model. The selection of this structural type is based on the high influence of the SSI phenomenon on the dynamic behaviour of integral bridges (Rodriguez et al., 2011).

This article is organized as follows. First, some basics about the use of dynamic Winkler models to numerically simulate the SSI phenomenon in civil engineering structures are presented. Subsequently, the parameter identification problem of the dynamic Winkler model is formulated in terms of a multi-objective optimization problem. Next, the performance of the proposed method is validated via its implementation for a real steel–concrete composite integral footbridge. Finally, some concluding remarks are drawn to close this article.

SSI under the dynamic Winkler model

An adequate simulation of the SSI is a basic aspect for the accurate estimation of the response of the civil engineering structure under dynamic actions (Tabatabaiefar et al., 2016). The SSI phenomenon depends clearly on the type of foundation. Two extreme cases are normally considered to simulate the behaviour of these foundations, the so-called: (1) direct and (2) deep foundations.

Direct foundations are usually rigid massive structures. Direct foundations may be modelled as an equivalent spring–damper element (Gazetas, 1991). The dynamic stiffness and damping coefficient of this element may be determined in terms of the amplitude of the vibration level of the structure. For low vibration levels, which are assumed herein, the SSI phenomenon on direct foundation may be simulated adequately considering only the stiffness component.

Thus, the value of the equivalent stiffness of an embedded direct foundation, either $k_{i}$ for translational degrees of freedom (N/m) or $k_{ri}$ for rotational degrees of freedom (N/rad) (i is the direction of the considered degree of freedom), may be determined (Gazetas, 1991) in terms of the equivalent stiffness of a surface direct foundation, either $k_{i, surf}$ for translational degrees of freedom (N/m) or $k_{ri, surf}$ for rotational degrees of freedom (N/rad), as follows (x is the longitudinal direction, y is the lateral direction and z is the vertical direction)

k_{z} = k_{z, surf} [1 + \frac{1}{21} \frac{D}{B} (1 + 1.3 χ)] [1 + 0.2 {(\frac{A_{w}}{A_{b}})}^{\frac{2}{3}}]

(1)

k_{x} = k_{x, surf} [1 + 0.15 \sqrt{\frac{D}{B}}] [1 + 0.52 {(\frac{d A_{w}}{B L^{2}})}^{0.4}]

(2)

k_{y} = k_{y, surf} [[1 + 0.15 \sqrt{\frac{D}{B}}] [1 + 0.52 {(\frac{d A_{w}}{B L^{2}})}^{0.4}]]

(3)

k_{rx} = k_{rx, surf} [1 + 1.26 \frac{d}{B} [1 + \frac{d}{B} {(\frac{d}{D})}^{- 0.2} \sqrt{\frac{B}{L}}]]

(4)

k_{ry} = k_{ry, surf} [1 + 0.92 {(\frac{d}{L})}^{0.6} [1.5 + {(\frac{d}{L})}^{1.9} {(\frac{d}{D})}^{- 0.6}]]

(5)

k_{rz} = k_{rz, surf} [1 + 1.4 (1 + \frac{B}{L}) {(\frac{d}{B})}^{0.9}]

(6)

k_{z, surf} = \frac{2 G_{s} L}{1 - ν_{s}} {0.73 + 1.54 χ^{0.75}}

(7)

k_{y, surf} = \frac{2 G_{s} L}{2 - ν_{s}} {2 + 2.5 χ^{0.85}}

(8)

k_{x, surf} = k_{y, surf} - \frac{0.2}{0.75 - ν_{s}} G_{s} L (1 - \frac{B}{L})

(9)

k_{rx, surf} = \frac{G_{s}}{1 - ν_{s}} {I_{bx}}^{0.75} {(\frac{L}{B})}^{0.25} (2.4 + 0.5 \frac{B}{L})

(10)

k_{ry, surf} = \frac{G_{s}}{1 - ν_{s}} {I_{by}}^{0.75} [3 {(\frac{L}{B})}^{0.15}]

(11)

k_{rz, surf} = G {J_{b}}^{0.75} [4 + 11 {(1 - \frac{B}{L})}^{10}]

(12)

where $G_{s} = E_{s} / 2 (1 + ν_{s})$ is the shear modulus of the soil (N/m²); $E_{s}$ is the longitudinal modulus of the soil (N/m²); $ν_{s}$ is the Poisson’s ratio of the soil (–), D is the total depth of the direct foundation (m), d is the depth of the direct foundation in contact with the soil (m); B is the half of the width of the direct foundation (m); L is the half of the length of the direct foundation (m); $A_{b} = 2 B \cdot 2 L$ is the lower surface of the direct foundation in contact with the soil (m²); $A_{w} = 2 L \cdot d$ is the lateral surface of the direct foundation in contact with the soil (m²); $I_{bx}$ is the moment of inertia of the lower surface of the direct foundation with respect to the x-axis (m⁴); $I_{by}$ is the moment of inertia of the lower surface of the direct foundation with respect to the y-axis (m⁴); $J_{b} = I_{bx} + I_{by}$ is the torsional moment of inertia of the lower surface of the direct foundation (m⁴); and $χ$ is the shape parameter defined as $χ = A_{b} / 4 L^{2}$ (–).

Deep foundations are slender flexible structures whose structural behaviour depends mainly on their length. The following two types of models have been usually used to analyse the SSI phenomenon in the case of deep foundations: (1) spring models and (2) spring–dampers models.

On one hand, the spring models present as main assumption a fully coupling between the soil and the deep foundation. According to this model, the effect of the soil adjacent to the deep foundation is simulated by several spring elements located along the deep foundation. The response of these spring models have been usually obtained via either approximate solutions (Novak, 1974; Tajimi, 1966) or numerical methods, such as the FE or the BE method (Alamo et al., 2016).

On the other hand, the spring–damper models, in which the effect of the soil adjacent to the deep foundation is simulated by several spring–damper elements located along its length. These second models have been widely used to simulate the SSI phenomenon of civil engineering structure for practical engineering applications due to the well-balanced equilibrium between their complexity and the accurate of the results provided. Among the different spring–damper models, the dynamic Winkler model has been considered herein.

According to this model, the stiffness of each spring–damper element, $k = k (z)$ , may be determined as follows

k (z) = k_{b} (z) \cdot L_{k} \cdot D_{p}

(13)

where $k (z)$ is the stiffness of the spring element (N/m) at depth z (m), $k_{b} (z)$ is the ballast coefficient (N/m³) at depth z (m), $L_{k}$ is the length between springs (m) and $D_{p}$ is the diameter of the deep foundation (m).

In addition, the damping coefficient of each spring–damper element, c (sN/m), may be determined via the sum of two energy dissipation mechanisms (Gazetas and Dobry, 1984) as follows: (1) hysteretic damping and (2) radiation damping

c = c_{h} + c_{r}

(14)

The hysteretic damping, $c_{h}$ (sN/m), also called the material damping, may be expressed as (Thavaraj, 2000)

c_{h} = \frac{2 β k_{b}}{ω}

(15)

where $β$ is the damping ratio of the soil (–); $k_{b}$ is the ballast coefficient of the soil (N/m³) and $ω$ is the angular frequency of the applied force (rad/s).

Similarly, the radiation damping coefficient, $c_{r}$ (sN/m), may be determined through an experimental relationship (Thavaraj, 2000). Among the different proposals, the relationship proposed by Berger et al. (1977) has been considered herein because it is an angular frequency independent expression

c_{r} = 4 \frac{D_{p}}{2} ρ_{s} v_{s} (1 + \frac{v_{c}}{v_{s}})

(16)

where $ρ_{s}$ is the density of the soil (kg/m³); $v_{s} = \sqrt{G_{s} / ρ_{s}}$ is the shear wave velocity (m/s) and $v_{c}$ is the velocity defined as $v_{c} = \sqrt{(2 / 1 - ν_{s})} v_{s}$ (m/s), where $ν_{s}$ is Poisson’s ratio of the soil (–).

Finally, it is necessary to remark that herein the contribution of the hysteretic damping coefficient, $c_{h}$ , has been neglected as assumption because it is two orders of magnitude lower than the radiation damping coefficient, $c_{r}$ , within the frequency range of interest, 0–20 Hz.

Parameter identification method of dynamic systems using a hybridUKF–MHS algorithm

Parameter identification problem consists in estimating the constitutive parameters of a system by considering both the actual measurements of its real response and the predictions provided by a dynamic numerical model which simulates its behaviour. Thus, the objective of the parameter identification problem is to infer the value of these constitutive parameters which minimizes the difference between the numerical and experimental behaviour of the system. For this reason, the parameter identification problem may be also referred as an inverse problem.

There are several methods to simulate numerically the response of the system during the parameter identification process. Among these methods, FE analysis has been widely used to model numerically the response of civil engineering structures. If the parameter identification method is based on an FE model, the solution of this parameter identification problem implies the updating of the corresponding FE model. Thus, this parameter identification problem may be transformed into an FE model updating problem (Marwala, 2010).

FE model updating of civil engineering structures for practical engineering applications is usually solved using some type of estimators. Estimators can be classified into two general groups (Chen, 2003; Rao et al., 2007) as follows: (1) point estimators and (2) interval estimators. The first type returns a single value of each considered parameter and the second type returns an interval of possible values of each considered parameter. The maximum likelihood method, a point estimator, has been considered herein to perform this problem due to its efficiency and accuracy in the solution of the updating problem (Jiménez-Alonso et al., 2017).

According to this method, the FE model updating problem may be formulated as a multi-objective optimization problem (Mottershead et al., 2011). Thus, the objective of this problem is to minimize the value of the different terms of the objective function via the modification of some pre-selected physical parameters of the model (considered as design variables). This objective function is defined in terms of the relative differences between the numerical and experimental modal properties of the structure. Computational intelligence algorithms are normally considered to solve this optimization problem (Marwala, 2010). However, these computational algorithms present the following two clear limitations: (1) they elapse a high simulation time when the complexity of the FE model increases; and (2) they are not able to deal with the uncertainty associated with the experimental modal properties of the system.

In order to overcome these limitations, hybrid algorithms are normally employed. Among these algorithms, a recent proposal – the UKF–MHS algorithm (Naranjo-Perez et al., 2020) – has been considered herein due to its great efficiency and accuracy to solve the FE model updating problem of complex civil engineering structures.

As result of the optimization process, the so-called Pareto front is obtained. The Pareto front shows the set of possible solutions of the optimization problem. Subsequently, a decision-making problem must be solved to select the best solution among all the elements of the Pareto front. The criterion provided by Jin et al. (2014) has been considered herein to cope with this decision-making problem.

In next sections, the proposed parameter identification method is described in detail. First, the formulation of the parameter identification problem as a multi-objective optimization problem has been presented. Subsequently, a theoretical background of the MHS and the UKF algorithms has been included. Finally, the proposed hybrid UKF–MHS algorithm has been described briefly.

Formulation of the parameter identification problem

As it has been mentioned above, the parameter identification problem may be formulated as a multi-objective optimization problem. The formulation of this problem may be expressed as follows

\begin{matrix} min (\begin{matrix} f_{1} (θ) & \begin{matrix} f_{2} (θ) & f_{3} (θ) \end{matrix} \end{matrix}) = \min \\ (\frac{1}{2} \begin{matrix} {[\sum_{j = 1}^{n_{f}} r_{f, j} {(θ)}^{2}]}^{\frac{1}{2}} & \begin{matrix} \frac{1}{2} {[\sum_{j = 1}^{n_{f}} r_{m, j} {(θ)}^{2}]}^{\frac{1}{2}} & \frac{1}{2} {[\sum_{j = 1}^{n_{f}} r_{d, j} {(θ)}^{2}]}^{\frac{1}{2}} \end{matrix} \end{matrix}) \\ s . t . θ_{l} \leq θ \leq θ_{u} \end{matrix}

(17)

where $r_{f, j} (θ)$ , $r_{m, j} (θ)$ and $r_{d, j} (θ)$ are the residuals of the natural frequency, vibration mode and damping ratio, respectively, of the j considered vibration mode; $n_{f}$ is the number of vibration modes considered in the analysis; $θ$ is the vector containing the pre-selected physical parameters of the FE model; and $θ_{l}$ and $θ_{u}$ are the lower and the upper bounds of the search domain of these physical parameters.

The three residuals may be determined in terms of the relative differences between the numerical and experimental modal properties of the structure as follows

r_{f, j} (θ) = \frac{f_{num, j} (θ) - f_{\exp, j}}{f_{\exp, j}}

(18)

r_{m, j} (θ) = \sqrt{\frac{{(1 - \sqrt{{MAC}_{\exp, j}^{num, j} (θ)})}^{2}}{{MAC}_{\exp, j}^{num, j} (θ)}}

(19)

r_{d, j} (θ) = \frac{ξ_{num, j} (θ) - ξ_{\exp, j}}{ξ_{\exp, j}}

(20)

where $f_{num, j} (θ)$ and $f_{\exp, j}$ are the numerical and experimental natural frequencies, respectively; $ξ_{num, j} (θ)$ and $ξ_{\exp, j}$ are the numerical and experimental damping ratios, respectively; and ${MAC}_{\exp, j}^{num, j} (θ)$ is the Modal Assurance Criterion ratio defined as (Allemang and Brown, 1982)

{MAC}_{\exp, j}^{num, j} (θ) = \frac{{(ϕ_{num, j} (θ) \cdot ϕ_{\exp, j})}^{2}}{(ϕ_{num, j}^{T} (θ) \cdot ϕ_{num, j} (θ)) \cdot (ϕ_{\exp, j}^{T} \cdot ϕ_{\exp, j})}

(21)

where $ϕ_{num, j} (θ)$ and $ϕ_{\exp, j}$ are the numerical and experimental vibration modes, respectively, and T denotes the transpose.

MHS algorithm

The HS algorithm is a global metaheuristic algorithm inspired by the mental process for the creation of musical harmony (Geem et al., 2001). The aim of the algorithm is to find the global minimum of a previously defined objective function by modifying a set of design variables (physical parameters) of a numerical model. The HS algorithm has been implemented successfully for several practical engineering applications (Yang and Koziel, 2011). This algorithm has shown a great effectiveness (a reduced simulation time without compromising the accuracy) in comparison with other conventional metaheuristic algorithms when it is implemented to solve nonlinear optimization problems (Wang et al., 2015). The operating rules of this algorithm are similar to the ones which govern the performance of other conventional metaheuristic algorithms, such as genetic algorithms. However, in this case, the complexity of the numerical operations employed to control the evolution of the algorithm has been reduced to improve its performance. Despite the great effectiveness of this algorithm, it has been rarely applied, to the best of the authors’ knowledge, for the parameter identification of systems based on the FE model updating method. The MHS algorithm is an extension of the previous algorithm which allows minimizing more than one objective function.

The MHS algorithm consists of the following steps:

The harmony matrix, $H$ , is initialized. This matrix contains the initial HMS generated solutions. Each solution represents a vector with the physical parameters to be updated. Subsequently, the different terms of the objective functions are evaluated for the different solutions in the harmony matrix, $H$ .

New harmonies are improvised. This step consists in generating a new set of harmonies (parameter vectors) based on: (a) memory considerations, (b) randomization and (c) pitch adjustments. According to this, each parameter can adopt a value from a previous value in the harmony matrix, $H$ , or be assigned randomly. The probability that a parameter is selected from the harmony matrix, $H$ , is controlled by the variable HMCR. Thus, for each parameter, a random number, between 0 and 1, is generated and if this number is lower than or equal to HMCR, the parameter is chosen randomly from the matrix $H$ . Otherwise, a random value among the possible values of the search domain is assigned to the parameter. When the parameter is defined from $H$ , it can be modified additionally. The probability – that this parameter can be modified – is determined by the variable PAR. A random number is generated again between 0 and 1 and if it is lower than or equal to PAR, the parameter value is adjusted by subtracting or adding a previously defined value, the so-called bandwidth, bw. Otherwise, the parameter value remains unchanged.

The new solutions and the previous solutions in the harmony matrix $H$ are then classified using both the non-dominating sorting technique (Deb et al., 2002) and the crowding distance. Hence, the so-called Pareto front is defined in terms of the non-dominated solutions. Finally, in order to restore the initial size of the harmony matrix, $H$ , the worst solutions, according to a crowding distance criterion, are removed.

These steps are repeated iteratively until a stop criterion is reached.

Thus, Figure 1 shows the flow chart of the MHS algorithm.

Figure 1.

Flow chart of the MHS algorithm.

UKF

The UKF algorithm is a derivative free estimator for nonlinear systems. This algorithm belongs to the class of filters named sigma points Kalman filters (Julier and Uhlmann, 1997). Its implementation as parameter estimator was proposed by Wan and Van Der Merwe (2000) and Wan et al. (1999).

The parameter identification can be addressed as

θ_{k} = θ_{k - 1} + w_{k - 1}

(22)

z_{k} = h (θ_{k}) + v_{k}

(23)

where $θ$ is a vector with the physical parameters to be estimated; $z$ is a vector with the outputs of the modelling function $h$ ; $w$ is a vector which takes into account the process noise; and $v$ is a vector which takes into account the observation noise.

Both types of noise are assumed to be white Gaussian with zero mean and uncorrelated covariance matrices $Q$ and $R$ , respectively. The matrix $R$ can be divided into two components (Tarantola, 2005) as follows: (1) the modelling noise covariance and (2) the measurement noise covariance. If the same FE model is considered for each step, the first component of the matrix $R$ can be neglected.

The UKF algorithm addresses the nonlinear estimation by considering $2 n_{d} + 1$ ( $n_{d}$ is the number of parameters) deterministic sampling points (so-called sigma points) to derive the posterior mean and covariance matrix, $P$ , via their propagation (the so-called unscented transformation) through the nonlinear function, $h$ (Julier and Uhlmann, 2004). Thus, the main computational effort of this algorithm is the computation of the new sigma points. In addition, it is necessary to remark as limitation that the posterior covariance matrix, $P$ , must be positive semidefinite, as these sigma points are obtained by the square root decomposition of this matrix. Even though the square root decomposition can be efficiently derived using the Cholesky factorization ( $A = \sqrt{P} = chol (P)$ where $P = A A^{T}$ and $T$ is the transpose function), the matrix $P$ is still updated at each iteration and numerical errors can give a non-positive-semidefinite matrix $P$ . The square root UKF algorithm, proposed by (Van Der Merwe and Wan, 2001), overcomes this problem as it avoids factorizing at each step. This algorithm propagates directly the matrix $A$ and guarantees that the covariance matrix, $P$ , is positive semidefinite.

The $2 n_{d} + 1$ sigma points are defined as

{(χ_{k - 1})}_{0} = {\hat{θ}}_{k - 1 | k - 1}

(24)

{(χ_{k - 1})}_{i} = {\hat{θ}}_{k - 1 | k - 1} + η {(A_{k - 1 | k - 1})}_{i} i = 1, 2, \dots, n_{d}

(25)

{(χ_{k - 1})}_{i + n_{d}} = {\hat{θ}}_{k - 1 | k - 1} - η {(A_{k - 1 | k - 1})}_{i} i = 1, 2, \dots, n_{d}

(26)

where $\hat{θ}$ are the posterior parameters estimated at the previous step and $η = \sqrt{(n_{d} + λ)}$ , where $λ$ is a scaling parameter. The sigma points have their associated weights

W_{0} = \frac{λ}{n_{d} + λ}

(27)

W_{i} = W_{i + n_{d}} = \frac{1}{2 (n_{d} + λ)} i = 1, 2, \dots, n_{d}

(28)

As a particular case of the KF, the square root UKF consists of two steps as follows: (1) the prediction step and (2) the correction (update) step.

First, the prediction step is conducted based on the prior model to assess the sigma points and predict the estimates of the estimation error covariance, $A^{θ}$ , the model output error covariance, $S^{z}$ , and the cross covariance, $P^{θ z}$ . The estimation error covariance is given by $A_{k} = γ^{- 1 / 2} A_{k - 1}$ , where $γ$ is a scalar factor slightly less than 1. Since $W_{i} > 0$ for all $i \geq 1$ , the model output error covariance, $S^{z}$ , can be expressed as

\begin{matrix} S_{k}^{z} = \sum_{0}^{2 n_{d}} W_{i} [(z_{k | k - 1}^{i} - {\hat{z}}_{k | k - 1}) \cdot {(z_{k | k - 1}^{i} - {\hat{z}}_{k | k - 1})}^{T}] \\ + R = [\sqrt{W_{i}} (z_{k | k - 1}^{i} - {\hat{z}}_{k | k - 1}), \sqrt{R}] \cdot \\ {[\sqrt{W_{i}} {(z_{k | k - 1}^{i} - {\hat{z}}_{k | k - 1})}^{T}, {\sqrt{R}}^{T}]}^{T} \\ + W_{0} [(z_{k | k - 1}^{0} - {\hat{z}}_{k | k - 1}) \cdot {(z_{k | k - 1}^{0} - {\hat{z}}_{k | k - 1})}^{T}] \\ for i = 1 : 2 n_{d} \end{matrix}

(29)

where $z_{k | k - 1}^{i}$ indicates that the variable $z^{i}$ was computed at the step k considering the data from the step $k - 1$ .

A $QR$ factorization,¹qr, can be used to express the first term as the product of an orthogonal matrix $O_{k} \in R^{2 n_{d} + r \times r}$ and an upper triangular matrix $S_{k}^{z} \in R^{r \times r}$ , where r is the number of measurements. The last term can be taken into account through a rank 1 update to the Cholesky factorization,2cholupdate (MATLAB R2019a, n.d.). Therefore, the matrix $S_{k}^{z}$ can be calculated as (Terejanu, 2011)

S_{k}^{z} = qr ([\sqrt{W_{1 : 2 n_{d}}} \cdot [{(z_{k | k - 1})}_{1 : 2 n_{d}} - {\hat{z}}_{k | k - 1}], \sqrt{R}])

(30)

S_{k}^{z} = cholupdate (S_{k}^{z}, {(z_{k | k - 1})}_{0} - {\hat{z}}_{k | k - 1}, sgn (W_{0}) W_{0})

(31)

Subsequently, the correction step uses both the measurements, $z^{obs}$ , and the model outputs to determine the posterior mean and covariance of the parameter estimation by considering the Kalman’s gain matrix, $K$ . The following algorithm is shown below.

Initial step:

{\hat{θ}}_{0} = θ_{prior}

(32)

A_{0}^{θ} = chol (P_{0}^{θ})

(33)

Main loop: $for k = 1 : N_{UKF}$ (number of iterations of the UKF)

Prediction step:

Calculate the $2 n_{d} + 1$ sigma points: ${(χ_{k - 1})}_{i}$

{(χ_{k | k - 1})}_{i} = {(χ_{k - 1})}_{i}

θ_{k | k - 1} = \sum_{0}^{2 n_{d}} W_{i} \cdot {(χ_{k | k - 1})}_{i}

A_{k | k - 1}^{θ} = γ^{- \frac{1}{2}} A_{k - 1}^{θ}

{(z_{k | k - 1})}_{i} = h ({(χ_{k | k - 1})}_{i})

{\hat{z}}_{k | k - 1} = \sum_{0}^{2 n_{d}} W_{i} \cdot {(z_{k | k - 1})}_{i}

S_{k | k - 1}^{z} = qr ([\sqrt{W_{1 : 2 n}} \cdot [{(z_{k | k - 1})}_{1 : 2 n} - {\hat{z}}_{k | k - 1}] \sqrt{R}])

S_{k | k - 1}^{z} = cholupdate (S_{k | k - 1}^{z}, {(z_{k | k - 1})}_{0} - {\hat{z}}_{k | k - 1}, sgn (W_{0}))

P_{k | k - 1}^{θ z} = \sum_{0}^{2 n_{d}} (W_{i} [{(χ_{k | k - 1})}_{i} - {\hat{θ}}_{k | k - 1}] \cdot {[{(z_{k | k - 1})}_{i} - {\hat{z}}_{k | k - 1}]}^{T})

Correction (update) step:

K_{k} = \frac{(P_{k | k - 1}^{θ z} / {S_{k | k - 1}^{z}}^{T})}{S_{k | k - 1}^{z}}

{\hat{θ}}_{k | k} = {\hat{θ}}_{k | k - 1} + K_{k} (z^{obs} - {\hat{z}}_{k | k - 1})

U = K_{k} S_{k | k - 1}^{z}

A_{k | k}^{θ} = cholupdate (A_{k | k - 1}^{θ}, U, - 1)

end

The hybrid UKF–MHS algorithm

The hybrid UKF–MHS algorithm is a local–global optimization algorithm which takes advantage of the virtues of the two well-known optimization algorithms to reduce the simulation time of the optimization process without compromising the accuracy of the solutions provided (Naranjo-Pérez et al., 2020).

The general layout of this algorithm is shown in Figure 2 and it may be summarized as follows: (1) the MHS algorithm proposes a set of candidate values for the parameters; (2) for each candidate value, the UKF algorithm obtains the value of $θ_{k | k}$ ; (3) the different terms of the objective function are evaluated for each value of $θ_{k | k}$ ; (4) the solutions are classified according to the non-dominating sorting technique; (5) the process is repeated iteratively until some convergence criterion is met; and (6) finally, the Pareto front is obtained.

Figure 2.

Flow chart of the hybrid UKF–MHS algorithm (Naranjo-Perez et al., 2020).

Application example: parameter identification of the dynamic Winkler model of a real integral footbridge

Description of the structure and preliminary finite element model

As benchmark structure, a real integral footbridge located at Paradas (Seville, Spain) has been considered herein (Figure 3(a)). The footbridge is a steel–concrete composite integral structure with a single span of 12 m of length. The cross-section of the structure consists in a steel–concrete composite box girder formed by a steel box girder of 0.40 m of depth and a concrete slab of 0.2 × 2.5 m². The two abutments are supported on four micropiles with a diameter, $D_{p}$ , of 0.2 m and a length, $L_{p}$ , of 15 m (Figure 3(b)).

Figure 3.

Lateral view of Paradas footbridge (Seville, Spain): (a) graphical representation and (b) footbridge in its current state.

A preliminary FE model of the footbridge was built using the FE analysis package Ansys (2019). The concrete slab, the steel box girder, the handrailing and the micropiles were modelled by three-dimensional (3D) beam elements (BEAM188). The mesh consists of 0.25 m long elements for the concrete slab and the steel box girder and 0.50 m long for the handrailing and the micropiles. The SSI has been simulated by one-dimensional (1D) spring–damper elements (COMBIN14) according to the model proposed by Krizek (2011) for this type of structure. A spring–damper element was placed at each node of the different micropiles. The tip of each micropile was assumed to be simply supported. The FE model of the overall structure is shown in Figure 4(c).

Figure 4.

Numerical models considered in this study: (a) geometric definition of the direct foundation of each abutment, (b) dynamic Winkler model for modelling the soil–pile interaction, (c) overall view of the FE model of the Paradas footbridge and (d) detail of the FE model showing the concrete–slab, steel box girder, abutment (as a massive node), handrailing and micropiles.

The mechanical properties of both the concrete and the steel were established according to the Eurocode 4 (1994). For the steel, density, $ρ_{st} = 7850 kg / m^{3}$ , Young’s modulus, $E_{st} = 2.1 \times 10^{11} N / m^{2}$ and Poisson’s ratio, $ν_{st} = 0.3$ ; and for the concrete, density, $ρ_{c} = 2500 kg / m^{3}$ , Young’s modulus, $E_{c} = 3 \times 10^{10} N / m^{2}$ and Poisson’s ratio, $ν_{c} = 0.2$ . A global structural damping ratio was considered for all the vibration modes. The value of this damping ratio was taken from the Synpex guidelines (Butz et al., 2007), which proposes $ξ_{num} = 0.6 %$ for steel–concrete composite footbridges.

According to the geotechnical report of the design project, the vertical profile of the soil was characterized by: (1) a first backfill material layer about 5 m of thickness and (2) a subsequent clay layer. The abutments were embedded around 0.9 m in the soil. The following mechanical properties were assumed for the simulation of the SSI phenomenon (according to the recommendations of the geotechnical report): (1) for the backfill material, density, $ρ_{s, 1} = 1800 kg / m^{3}$ , Young’s modulus, $E_{s, 1} = 7.2 \times 10^{8} N / m^{2}$ , Poisson’s ratio, $ν_{s, 1} = 0.35$ and the ballast coefficient, $k_{b 1} = 4 \times 10^{7} N / m^{3}$ ; (2) for the clay layer, density, $ρ_{s, 2} = 1800 kg / m^{3}$ , Young’s modulus, $E_{s, 2} = 4 \times 10^{7} N / m^{2}$ , Poisson’s ratio, $ν_{s, 2} = 0.3$ and the ballast coefficient, $k_{b 2} = 1.2 \times 10^{8} N / m^{3}$ .

The SSI phenomenon was simulated considering the following two interaction submodels: (1) the soil–abutment interaction (direct foundation) and (2) the soil–pile interaction (deep foundation). The soil–abutment interaction was simulated by six spring elements (three translations and three rotations) in each abutment. The stiffness of these springs has been calculated using equations (1) to (12) according to the proposal of Gazetas (1991). In addition, each abutment was modelled by a massive node which simulated its mechanical properties (Figure 4(a) and (d)).

The soil–pile interaction was simulated via the dynamic Winkler model (Figure 4(b)). The equivalent spring–damper elements were placed along the length of the pile. Thus, the stiffness, k, of each spring–damper element was computed using equation (13) and the damping, c, of each spring–damper element was computed using equation (16) which relates this magnitude with the mechanical properties of the soil.

Subsequently, a numerical FE damped modal analysis was performed in order to obtain the modal properties of the footbridge. As result of this analysis, Table 1 shows the first three numerical natural frequencies of the footbridge, $f_{num, j} (θ)$ (j is the considered vibration mode). In addition, Figure 5(b) illustrates the first three numerical vibration modes of the considered footbridge (labelled as $Num .$ ).

Table 1.

Numerical, $f_{num, j} (θ)$ , and experimental, $f_{\exp, j}$ , natural frequencies, numerical, $ζ_{num, j} (θ)$ , and experimental, $ζ_{\exp, j}$ , damping ratios, relative differences between natural frequencies, $Δ f_{\exp, j}^{num, j} (θ)$ , and damping ratios, ${Δ ζ}_{\exp, j}^{num, j} (θ)$ , and the ${MAC}_{\exp, j}^{num, j} (θ)$ ratios.

Mode	$f_{num, j} (θ) (Hz)$	$f_{\exp, j} (Hz)$	$ζ_{num, j} (θ) (%)$	$ζ_{\exp, j} (%)$	$Δ f_{e x p, j}^{n u m, j} (θ) (%)$	$Δ ζ_{e x p, j}^{n u m, j} (θ) (%)$	${MAC}_{\exp, j}^{num, j} (θ) (-)$
1	15.322	15.602	3.103	4.650	1.795	33.268	0.975
2	32.179	35.971	1.121	1.049	10.542	6.863	0.692
3	38.614	40.017	0.700	4.846	3.506	85.555	0.907

Figure 5.

(a) Layout of the ambient vibration test and (b) Numerical ( $Num .$ ), experimental ( $Exp .$ ) and updated ( $Upd .$ ) vibration modes.

Ambient vibration test, operational modal analysis and parameter identification

Subsequently, the modal properties of the footbridge were identified experimentally via the signal processing of the measures recorded during an ambient vibration test. In this experimental test, the response of the structure under ambient vibration was recorded at different points. Concretely, a gridline with 2 × 7 points was instrumented. The points were separated longitudinally 2 m and transversally 2.5 m. Three high sensitive triaxial force-balanced accelerometers (Kinemetrics EpiSensor ES-T) were used during this ambient test. Two of these devices were considered as references, at two fixed locations, and the other accelerometer was moving to the defined points (Figure 5(a)). Twelve set-ups, with a duration of 600 s and a sampling frequency of 100 Hz, were performed.

The experimental identification of the modal properties was conducted in time domain using the unweighted principal component-Merge algorithm, as it is implemented in the software Artemis (2016). This identification method, framed within the stochastic subspace identification method (Magalhães and Cunha, 2011), merges the records of the different set-ups which allows improving the estimates of the experimental modal properties of the footbridge (Döhler et al., 2010). The three first identified natural frequencies, $f_{\exp, j}$ , and associated damping ratios, $ζ_{\exp, j}$ (Magalhães et al., 2010) are shown in Table 1. Figure 5(b) also illustrates the first three experimental vibration modes of the considered footbridge (labelled as $Exp .$ ).

In order to determine the correlation between the numerical and experimental modal properties of this footbridge, the relative differences between natural frequencies, $Δ f_{\exp, j}^{num, j} (θ)$ , and the ${MAC}_{\exp, j}^{num, j} (θ)$ ratios are also shown in Table 1.

Thus, Table 1 shows a bad correlation between the numerical and experimental modal properties of the footbridge: (1) the relative difference between the second numerical and experimental natural frequency, $Δ f_{\exp, 2}^{num, 2} (θ)$ , is greater than 5%; (2) the ${MAC}_{\exp, 2}^{num, 2} (θ)$ ratio of the second vibration mode is clearly lower than 0.90 (Živanović et al., 2007); and (3) the relative differences of the damping ratios among the first and third vibrations modes are greater than 30%.

Hence, the value of the physical parameters of the FE model originally considered to simulate numerically the dynamic behaviour of the footbridge is not adequate, and a parameter identification method must be used to improve the performance of this numerical model.

In this manner, the above-mentioned parameter identification method was implemented in the mathematical package MATLAB (R2019a, n.d.) to perform the identification process. As one of the main factors, which may have a greater influence on the dynamic behaviour of the structure, is the SSI phenomenon, the parameter identification method focused on the parameter identification of the dynamic Winkler model.

Thus, the main objective of this identification process is to estimate the value of the most relevant physical parameters of the numerical model which minimizes the relative differences between the numerical and experimental modal properties of the footbridge. In order to select these parameters, two criteria were considered as follows: (1) to consider those physical parameters which mobilize more modal strain energy (Fox and Kapoor, 1968); and (2) to take into account those parameters which characterize the SSI phenomenon. As a result of this selection process, a vector of six parameters, $θ$ , was considered (Figure 4) as follows: (1) Young’s modulus of the concrete, $E_{c}$ ; (2) Young’s modulus of the steel, $E_{st}$ ; (3) Young’s modulus of the backfill, $E_{s, 1}$ ; (4) Young’s modulus of the clay, $E_{s, 2}$ ; (5) the ballast coefficient of the backfill, $k_{b, 1}$ ; and (6) the ballast coefficient of the clay, $k_{b, 2}$ . In relation to these parameters, it is necessary to remark that there is a clear relationship between Young’s modulus of the soil and the damping coefficient of the dynamic Winkler model (as it has been expressed by equation (16)). Therefore, it may be assumed that the first two parameters characterize the behaviour of the structure and the remaining four parameters characterize the SSI phenomenon.

In addition, a search domain has been established in order to avoid ill-conditioning problems and to ensure that the estimative of these parameters maintains an adequate physical meaning. In this sense, the following three facts must be remarked: (1) the upper bound of Young’s modulus of the steel is strangely high because the stiffeners and diaphragms have been modelled in a simplified way, only an equivalent added mass has been included for this purpose; (2) the lower and the upper bounds of the search domain of the parameters which characterize the behaviour of the soil have been established according to the recommendations of the European standards (Eurocode 7, 2004); and (3) the upper bound of Young’s modulus of the backfill was increased in order to simulate the effect of the dynamic compaction on this layer.

The lower, $θ_{l}$ , and upper, $θ_{u}$ , bounds of the search domain of each considered physical parameter are shown in Table 2.

Table 2.

Bounds (lower, $θ_{l}$ , and upper, $θ_{u}$ , bounds) of the search domain for each considered physical parameter and the corresponding updated value, $θ_{upd}$ , after the identification process.

$θ$	Definition	$θ_{l}$	$θ_{upd}$	$θ_{u}$
$θ_{1}$	Concrete Young’s modulus $E_{c}$ (N/m²)	$2.7 \times 10^{10}$	$3.67 \times 10^{10}$	$4 \times 10^{10}$
$θ_{2}$	Steel Young’s modulus, $E_{st}$ (N/m²)	$1.9 \times 10^{11}$	$2.28 \times 10^{11}$	$2.4 \times 10^{11}$
$θ_{3}$	Backfill Young’s modulus, $E_{s, 1}$ (N/m²)	$3 \times 10^{6}$	$7.77 \times 10^{8}$	$4 \times 10^{9}$
$θ_{4}$	Clay Young’s modulus, $E_{s, 2}$ (N/m²)	$3 \times 10^{6}$	$3.5 \times 10^{7}$	$4.8 \times 10^{7}$
$θ_{5}$	Backfill ballast coefficient, $k_{b, 1}$ (N/m³)	$1.3 \times 10^{7}$	$4.16 \times 10^{7}$	$2.1 \times 10^{8}$
$θ_{6}$	Clay ballast coefficient, $k_{b, 2}$ (N/m³)	$1.3 \times 10^{7}$	$1.27 \times 10^{8}$	$2.1 \times 10^{8}$

An initial population of 40 harmonies (parameter vectors) were generated. For each harmony, $θ_{k}$ , the UKF algorithm obtains the value of $θ_{k | k}$ . For this purpose, the number of iterations of the UKF algorithm is set to 5. The measurement noise covariance matrix is derived from the standard deviations provided by the experimental identification of the modal properties. The vector of observed measurements, $d_{obs}$ , was defined as follows: (1) the first three terms and the last three terms corresponded to the three experimental natural frequencies and modal damping ratios, respectively; (2) the second three terms were associated with the mode shapes; and (3) as the residual of the mode shapes is defined in terms of the MAC ratios and therefore the objective is to make these values closer to one, then these three second terms were equal to one. The parameter, $λ$ , which defines the sigma points, was defined as $λ = 0.0001$ according to the recommendations provided by Astroza et al. (2016) and the scalar factor was set as, $γ = 0.99$ , following the recommendations provided by Van Der Merwe and Wan (2001). The non-dominated solutions are stored in the matrix $H$ . Then, an iterative process with a total number of iterations of 20 was run.

In each iteration, a new set of 20 individuals (harmonies) was created. For each harmony, the UKF algorithm calculates $θ_{k | k}$ . Then, the above-mentioned FE model was updated considering the parameters $θ_{k | k}$ as inputs. Subsequently, the numerical natural frequencies, mode shapes and damping ratios were obtained considered this updated FE model. Next, the three objective functions were computed and the dominated solutions were neglected. When the stop criterion (maximum number of iterations) was reached, the optimal Pareto front was obtained.

Finally, a decision-making problem was solved to select the best solution among the different elements of the Pareto front. As it was mentioned above, the criterion provided by Jin et al. (2014) has been considered herein. Thus, the point estimation, $θ_{upd}$ , of the physical parameters of the model after the identification process is shown in Table 2. In addition, Figure 5(b) illustrates the first three numerical vibration modes of the considered footbridge after the parameter identification process (labelled as $Upd .$ ).

Once the identification process was finished, it was possible to estimate the value of the parameters of the dynamic Winkler model. Thus, the parameter of the spring–damper elements for each layer may be estimated as follows: (1) for the backfill layer, $k_{1} = 4.16 \times 10^{6} N / m$ and $c_{1} = 7.93 \times 10^{5} sN / m$ ; and (2) for the clay layer, $k_{2} = 1.27 \times 10^{7} N / m$ and $c_{2} = 1.68 \times 10^{5} sN / m$ .

Thus, the adequate performance of the proposed identification method is supported by the following three facts:

The good correlation between the numerically updated and experimental modal properties of the footbridge that follows the identification process: the relative differences, $Δ f_{\exp, j}^{upd, j} (θ)$ , between the numerical and experimental natural frequencies were lower than 5%, the relative differences, ${Δ ξ}_{\exp, j}^{upd, j} (θ)$ , between the numerical and experimental damping ratios were lower than 3%, and the ${MAC}_{\exp, j}^{upd, j} (θ)$ ratios were greater than 0.90.

The good agreement observed between the parameters of the dynamic Winkler model provided by this proposal and the values recommended by European standards (Eurocode 7, 2004).

The good correlation between the damping coefficient of the soil obtained by this identification method and the values reported by other studies in the literature (Hashash and Park, 2001).

Finally, one additional conclusion can be drawn from this study: the importance of considering the SSI phenomenon in the modelling of integral bridges. As Table 3 illustrates, the experimental damping ratios associated with the three identified vibration modes are greater than the values recommended by international guidelines (Butz et al., 2007). Therefore, the adequate characterization of the dynamic Winkler model plays a key role when simulating the actual behaviour of this type of civil engineering structures.

Table 3.

Updated, $f_{upd, j} (θ)$ , and experimental, $f_{\exp, j}$ , natural frequencies, updated, $ζ_{upd, j} (θ)$ , and experimental, $ζ_{\exp, j}$ , damping ratios, relative differences between natural frequencies, $Δ f_{\exp, j}^{upd, j} (θ)$ , and damping ratios, ${Δ ζ}_{\exp, j}^{upd, j} (θ)$ , and the ${MAC}_{\exp, j}^{upd, j} (θ)$ ratios.

Mode	$f_{upd, j} (θ) (Hz)$	$f_{\exp, j} (Hz)$	$ζ_{upd, j} (θ) (%)$	$ζ_{\exp, j} (%)$	$Δ f_{e x p, j}^{u p d, j} (θ) (%)$	$Δ ξ_{e x p, j}^{u p d, j} (θ) (%)$	${MAC}_{\exp, j}^{upd, j} (θ) (-)$
1	16.380	15.602	4.588	4.650	4.987	1.333	0.975
2	35.377	35.971	1.058	1.049	1.651	0.858	0.907
3	40.599	40.017	4.701	4.846	1.454	2.992	0.964

Conclusion

In this article, the experimental identification of the parameters of the dynamic Winkler SSI model (stiffness and damping coefficients) has been addressed based on the extension of a recently developed hybrid UKF–MHS algorithm. Implementation details of the proposed approach have been discussed in detail, and its performance has been illustrated via a real engineering case study. According to this proposal, the parameter identification problem is transformed successively into two problems as follows: (1) an FE model updating problem and (2) a multi-objective optimization problem.

In this manner, the parameter identification problem has been formulated in terms of a multi-objective optimization problem, where the objective function to be minimized is defined in terms of the relative differences (residuals) between the numerical and experimental modal properties of the soil–structure system. To this end, the following three types of residuals have been considered: (1) residuals associated with the natural frequencies; (2) residuals associated with the damping ratios; and (3) residuals associated with the vibration modes. The numerical behaviour of the system has been simulated using the FE method, while its experimental modal properties have been identified via the signal processing of the dynamic response of the structure recorded during an experimental vibration test. As design variables, the most relevant physical parameters of the structure were adopted. The optimization has been conducted using a hybrid UKF–MHS algorithm.

The performance of the method has been validated via its application to a real steel–concrete composite integral footbridge, for which the parameters of the dynamic Winkler model have been identified experimentally. The obtained results illustrate the following: (1) the good agreement between the numerical and experimental modal properties of the structure after the identification process; and (2) the good correlation between the values of the Winkler parameters obtained according to this proposal and the values recommended by the European standards (Eurocode 7, 2004).

The proposed approach exhibits the following three main advantages when addressing the parameter identification of the dynamic Winkler SSI model: (1) the parameter identification can be performed indirectly based on the solution of an inverse problem defined in terms of the overall numerical and experimental behaviour of the structure; (2) no complex geotechnical test is needed to perform the parameter identification; and (3) the method is quite simple and may be straightforwardly implemented for practical engineering applications.

In conclusion, the proposed method becomes a valid tool to characterize the parameters of the dynamic SSI Winkler model for practical engineering applications. Nevertheless, additional studies are recommended in order to better characterize the parameters of the dynamic Winkler model for different types of civil engineering structures supported on different types of soils.

Footnotes

Acknowledgements

The authors acknowledge the Planning Department of Paradas (Spain) for the support during the development of the experimental tests.

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship and/or publication of this article: This work was partially funded by the Ministerio de Economía y Competitividad of Spain and the European Regional Development Fund under project RTI2018-094945-B-C21. In addition, the co-author, J.N.-P., has been supported by the research contract, USE-17047-G, provided by the Universidad de Sevilla.

Notes

ORCID iDs

Javier Naranjo-Pérez

Javier Fernando Jiménez-Alonso

References

Alamo

Martínez-Castro

Padrón

, et al. (2016) Efficient numerical model for the computation of impedance functions of inclined pile groups in layered soils. Engineering Structures 126: 379–390.

Allemang

Brown

(1982) A correlation coefficient for modal vector analysis. In: Proceedings of the 1st International Modal Analysis Conference, Vol. 1, Orlando, FL, 8–10 November, pp. 110–116. Orlando, FL: Union College Press.

Allotey

El Naggar

(2008) Generalized dynamic Winkler model for nonlinear soil–structure interaction analysis. Canadian Geotechnical Journal 45(4): 560–573.

Ansys (2019) Ansys Mechanical 19.0. Available at: https://www.ansys.com/

Arai

Ohta

Kojima

(1984) Estimation of soil parameters based on monitored movement of subsoil under consolidation. Soils and Foundations 24(4): 95–108.

Artemis (2016) Artemis Modal. Available at: http://www.svibs.com/

Astroza

Nguyen

Nestorović

(2016) Finite element model updating using simulated annealing hybridized with unscented Kalman filter. Computers and Structures 177: 176–191.

Badsar

Schevenels

Haegeman

, et al. (2010) Determination of the material damping ratio in the soil from SASW tests using the half-power bandwidth method. Geophysical Journal International 182: 1493–1508.

Berger

Mahi

Pyke

(1977) Simplified method for evaluating soil-pile-structure interaction effects. In: Proceedings of the 9th Offshore Technology Conference, Houston, TX, 2–5 May.

10.

Butz

Heinemeyer

Goldack

, et al. (2007) Advanced Load Models for Synchronous Pedestrian Excitation and Optimised Design Guidelines for Steel Footbridges (SYNPEX). RFCS-Research Project RFS-CR-03019. Available at: https://op.europa.eu/en/publication-detail/-/publication/432f1d42-ed71-47db-9087-c5b74b4bb0fe

11.

Campanella

Stewart

(1992) Seismic cone analysis using digital signal processing for dynamic site characterization. Canadian Geotechnical Journal 29: 477–486.

12.

Chen

(2003) Bayesian filtering: from Kalman filters to particle filters, and beyond. Statistics 182: 1–69.

13.

Cheng

Leong

(2018) Determination of damping ratios for soils using bender element tests. Soil Dynamic and Earthquake Engineering 111: 8–13.

14.

Deb

Pratap

Agarwal

, et al. (2002) A fast and elitist multi-objective genetic algorithm: NSGA-II. IEEE Transactions on Evolutionary Computation 6: 182–197.

15.

Döhler

Andersen

Mevel

(2010) Data merging for multi-setup operational modal analysis with data driven SSI. In: 28th International Modal Analysis Conference, Jacksonville, FL, 1–4 February.

16.

Drnevich

Hardin

Shippy

(1978) Modulus and damping of soils by the resonant-column method. Dynamic Geotechnical Testing 654: 91–125.

17.

Eurocode 4 (1994) Design of Composite Steel and Concrete Structures. Brussels: CEN.

18.

Eurocode 7 (2004) Geotechnical Design – Part 1: General Rules. Brussels: CEN.

19.

Fox

Kapoor

(1968) Rates of change of eigenvalues and eigenvectors. AIAA Journal 6: 2426–2429.

20.

Gazetas

(1991) Formulas and charts for impedances of surface and embedded foundations. Journal of Geotechnical Engineering 117: 1363–1381.

21.

Gazetas

Dobry

(1984) Horizontal response of piles in layered soils. Journal of Geotechnical Engineering 110: 20–40.

22.

Geem

Kim

Loganathan

(2001) A new heuristic optimization algorithm: harmony search. Simulation 76: 60–68.

23.

Gens

Ledesma

Alonso

(1996) Estimation of parameters in geotechnical backanalysis: II. Application to a tunnel excavation problem. Computers and Geotechnics 18(1): 29–46.

24.

Hashash

YMA

Park

(2001) Non-linear one-dimensional seismic ground motion propagation in the Mississippi embayment. Engineering Geology 62: 185–206.

25.

Hirai

(2011) A Winkler model approach for vertically and laterally loaded piles in nonhomogeneous soil. International Journal for Numerical and Analytical Method in Geomechanics 36(17): 1869–1897.

26.

Jazwinski

(1970) Stochastic Processes and Filtering Theory. New York: Academic Press.

27.

Jiménez-Alonso

Hudson

Pavic

, et al. (2017) Maximum likelihood methods for finite element model updating of civil engineering structures: a comparative study. In: Proceedings of 4th International Conference on Mechanical Models in Structural Engineering, Madrid, November–December.

28.

Jin

Cho

Jung

, et al. (2014) A new multi-objective approach to finite element model updating. Journal of Sound and Vibration 333(11): 2323–2338.

29.

Julier

Uhlmann

(1997) New extension of the Kalman filter to nonlinear systems. In: Signal Processing, Sensor Fusion, and Target Recognition VI, Vol. 3068, pp. 182–194. International Society for Optics and Photonics. DOI: https://doi.org/10.1117/12.280797.

30.

Julier

Uhlmann

(2004) Corrections to ‘unscented filtering and nonlinear estimation’. Proceedings of the IEEE 92(12): 1958–1958.

31.

Kalman

(1960) A new approach to linear filtering and prediction problems. Journal of Basic Engineering 82: 35–45.

32.

Karl

Haegeman

Degrande

(2006) Determination of the material damping ratio and the shear wave velocity with the seismic cone penetration test. Soil Dynamic and Earthquake Engineering 26: 1111–1126.

33.

Krizek

(2011) Soil–structure interaction of integral bridges. Structural Engineering International 21: 169–174.

34.

Lai

Rix

(1998) Simultaneous Inversion of Rayleigh Phase Velocity and Attenuation for Near-Surface Site Characterization. Atlanta, GA: Georgia Institute of Technology, National Science Foundation, and US Geological Survey.

35.

Ledesma

Gens

Alonso

(1996a) Estimation of parameters in geotechnical backanalysis – I. Maximum likelihood approach. Computers and Geotechnics 18: 1–27.

36.

Ledesma

Gens

Alonso

(1996b) Parameter and variance estimation in geotechnical backanalysis using prior information. International Journal for Numerical and Analytical Methods in Geomechanics 20: 119–141.

37.

Lee

J-H

Lee

J-J

Choi

J-S

, et al. (2016) A semi-analytical approach to predict ground vibration by identification of soil properties and train-transit loads. Advances in Structural Engineering 15(6): 1013–1029.

38.

Magalhães

Cunha

(2011) Explaining operational modal analysis with data from an arch bridge. Mechanical Systems and Signal Processing 25: 1431–1450.

39.

Magalhães

Cunha

Caetano

, et al. (2010) Damping estimation using free decays and ambient vibration tests. Mechanical Systems and Signal Processing 24: 1274–1290.

40.

Maia

NMM

Silva

JMM

, et al. (1997) Theoretical and Experimental Modal Analysis. Somerset: Research Studies Press Ltd.

41.

Maier

Gioda

(1982) Optimization methods for parametric identification of geotechnical systems. In: Martins

(ed.) Numerical Methods in Geomechanics (NATO Advanced Study Institutes Series (Series C – Mathematical and Physical Sciences Proceedings of the Advanced Study Institute Programme), Vol. 92. Dordrecht: Springer, pp. 273–304.

42.

Marwala

(2010) Finite-Element-Model Updating Using Computational Intelligence Techniques: Applications to Structural Dynamics. London: Springer, p. 231.

43.

MATLAB R2019a (n.d.) Available at: http://www.mathworks.com/

44.

Mottershead

Link

Friswell

(2011) The sensitivity method in finite element model updating: a tutorial. Mechanical Systems and Signal Processing 25: 2275–2296.

45.

Naranjo-Pérez

Jiménez-Alonso

Pavic

, et al. (2020) Finite-element-model updating of civil engineering structures using a hybrid UKF-HS algorithm. Structure and Infrastructure Engineering. Epub ahead of print. DOI:10.1080/15732479.2020.1760317.

46.

Nguyen

Nestorović

(2015) Nonlinear Kalman filters for model calibration of soil parameters for geomechanical modeling in mechanized tunneling. Journal of Computing in Civil Engineering 30(2): 04015025.

47.

Novak

(1974) Dynamic stiffness and damping of piles. Canadian Geotechnical Journal 11: 574–598.

48.

Rao

Toutenburg

Fieger

(2007) Linear Models and Generalizations: Least Squares and Alternatives. 2nd ed. New York: Springer.

49.

Rodriguez

Martinez

Marti

(2011) Integral bridge for high-speed railway. Structural Engineering International 21: 297–303.

50.

Srinivas

Deb

(1994) Multiobjective optimization using nondominated sorting in genetic algorithms. Evolutionary Computation 2: 221–248.

51.

Tabatabaiefar

SHR

Fatahi

Samali

(2016) Numerical and experimental investigations on seismic response of building frames under influence of soil-structure interaction. Advances in Structural Engineering 17(1): 109–130.

52.

Tajimi

(1966) Earthquake response of foundation structures. Report of Faculty of Science and Engineering, Nihon University, Tokyo, Japan. (In Japanese)

53.

Tarantola

(2005) Inverse Problem Theory and Methods for Model Parameter Estimation, Vol. 89. Philadelphia, PA: Society for Industrial and Applied Mathematics (SIAM).

54.

Terejanu

(2011) Unscented Kalman Filter Tutorial. Buffalo, NY: University at Buffalo.

55.

Thavaraj

(2000) Seismic analysis of pile foundations for bridges. PhD thesis, University of British Columbia, Vancouver, BC, Canada.

56.

Van Der Merwe

Wan

(2001) The square-root unscented Kalman filter for state and parameter-estimation. In: 2001 IEEE International Conference on Acoustics, Speech, and Signal Processing – Proceedings, Salt Lake City, UT, 7–11 May, pp. 3461–3464. New York: IEEE.

57.

Wan

Van Der Merwe

(2000) The unscented Kalman filter for nonlinear estimation. In: Proceedings of the IEEE 2000 Adaptive Systems for Signal Processing, Communications, and Control Symposium, Lake Louise, AB, Canada, 4 October, pp. 153–158. New York: IEEE.

58.

Wan

Van Der Merwe

Nelson

(1999) Dual estimation and the unscented transformation. In: Proceedings of the 12th International Conference on Neural Information Processing Systems, pp. 666–672. MA, USA: MIT Press.

59.

Wang

Gao

X-Z

Zenger

(2015) An Introduction to Harmony Search Optimization Method. Cham: Springer International Publishing.

60.

Yang

X-S

Koziel

(2011) Computational Optimization and Applications in Engineering and Industry, Vol. 359. Berlin: Springer Science + Business Media.

61.

Yazdchi

Khalili

Valliappan

(1999) Dynamic soil–structure interaction analysis via coupled finite-element–boundary-element method. Soil Dynamic and Earthquake Engineering 18: 499–517.

62.

Živanović

Pavic

Reynolds

(2007) Finite element modelling and updating of a lively footbridge: the complete process. Journal of Sound and Vibration 301(1–2): 126–145.

Parameter identification of the dynamic Winkler soil–structure interaction model using a hybrid unscented Kalman filter–multi-objective harmony search algorithm

Abstract

Keywords

Introduction

SSI under the dynamic Winkler model

Parameter identification method of dynamic systems using a hybridUKF–MHS algorithm

Formulation of the parameter identification problem

MHS algorithm

UKF

The hybrid UKF–MHS algorithm

Application example: parameter identification of the dynamic Winkler model of a real integral footbridge

Description of the structure and preliminary finite element model

Ambient vibration test, operational modal analysis and parameter identification

Conclusion

Footnotes

Acknowledgements

Declaration of Conflicting Interests

Funding

Notes

ORCID iDs

References