Regularization-Based Dual Adaptive Kalman Filter for Identification of Sudden Structural Damage Using Sparse Measurements

Featured Application

The dual adaptive filter proposed in this paper can identify sudden change in structural systems under dynamic excitations. The proposed filter method including the tuning process can be applied to a variety of engineering areas in which near-real-time tracking of system parameter is needed.

Abstract

This paper proposes a dual adaptive Kalman filter to identify parameters of a dynamic system that may experience sudden damage by a dynamic excitation such as earthquake ground motion. While various filter techniques have been utilized to estimate system’s states, parameters, input (force), or their combinations, the filter proposed in this paper focuses on tracking parameters that may change suddenly using sparse measurements. First, an advanced state-space model of parameter estimation employing a regularization technique is developed to overcome the lack of information in sparse measurements. To avoid inaccurate or biased estimation by conventional filters that use covariance matrices representing time-invariant artificial noises, this paper proposes a dual adaptive filtering, whose slave filter corrects the covariance of the artificial measurement noises in the master filter at every time-step. Since it is generally impossible to tune the proposed dual filter due to sensitivity with respect to parameters selected to describe artificial noises, particle swarm optimization (PSO) is adopted to facilitate optimal performance. Numerical investigations confirm the validity of the proposed method through comparison with other filters and emphasize the need for a thorough tuning process.

1. Introduction

In various engineering fields, identifying the uncertain status of structural systems is a critical task to avoid severe malfunctions due to degradation of system parameters caused by deterioration or sudden damage by natural or man-made disasters. Among various approaches to such system identification (SI), there have been notable research efforts for probabilistic SI to handle uncertainties in the problem. For instance, Tarantola (1987) formulated SI as a Bayesian inverse problem, in which a joint probability density function (PDF) was introduced for parameters of systems under static loads [1]. Gaussian assumption was used for both prior and posterior distributions. Recently, such approach was extended by several probabilistic SI methods utilizing Bayesian Network, Markov Chain Monte Carlo, and Branch-and-Bound method [2,3,4,5,6]. Although these methods can handle general distribution models for prior and posterior, they are still subject to fundamental limitations for dynamic systems, i.e., systems under dynamic loads.

For precise updates of knowledge of a dynamic system based on the given measurements, sequential Bayesian filters including Kalman filter have been utilized [7,8,9,10]. The main strength of Bayesian filters comes from their two complementary updates: time-update and measurement-update. In this approach, by employing models that are respectively describing time- and measurement-updates using common random variables, a given dynamic system does not need to be expressed as one highly complicated joint PDF. To apply a filter method to dynamic system identification, it is vital to build an accurate state-space model that represents the dynamic system by two numerical models respectively describing the time- and measurement-updates. Then, a suitable filter method should be selected especially when the state-space model is nonlinear or time-variant.

To track the status of nonlinear dynamic systems, the extended Kalman filter (EKF) has been utilized [7,11]. EKF uses first-order linearization to approximate underlying statistics of the state-space model, but this filter may suffer from some limitations including the requirement of sufficient differentiability of the dynamic model. To overcome this inherent problem, the unscented Kalman filter (UKF) was developed using so-called unscented transformation [10,12]. This transformation linearizes nonlinear dynamic model by using “sigma points” instead of the derivatives. Using the UKF as the basic filter method, this paper develops a near-real-time SI method for dynamic systems that may experience time-varying changes in system parameters caused by damage. To handle this target problem with sparsely measured accelerations, the state-space model of this study, which is generally used in parameter estimation only [13], utilizes the Newmark model to update modal-displacement, velocity, and acceleration with current system parameters for just one time-update. Then, regularization techniques with an assumed prior information about system are introduced to the state-space model to avoid divergence in estimations caused by sparse measurements.

Even if a proper state-space model is constructed for time-varying systems, conventional filters such as EKF and UKF using time-invariant artificial noises may not be able to achieve accurate estimations. This is because the covariance matrices of artificial noises, which are arbitrarily chosen to satisfy the optimality condition of filter during a certain time interval, may not be optimal as the statistical characteristics of measurements change over the excitation period. To address this issue, the UKF should be able to determine new artificial noises at every time step to satisfy the optimality condition. The authors thus explored several adaptive filtering techniques, which can adapt to time-varying statistical characteristics of measurements [11,13,14,15,16,17]. Through these investigations, the dual adaptive filtering for artificial measurement noise estimation is finally selected and used along with the state-space model.

It is also noted that the dual adaptive filtering, which adopts a slave filter to update artificial-measurement-noise of the master filter, requires additional selection of artificial noise covariance matrices for the slave filter. The increased number of covariance matrices of artificial noises makes it challenging to tune dual adaptive filter for various choices of those covariance matrices. It is difficult to select proper artificial noises of slave filter empirically due to the absence of accurate knowledge [13,15,18]. Therefore, it is important to construct a general framework for tuning process of dual adaptive filtering. For this purpose, particle swam optimization (PSO), which was utilized in a recent study to tune Kalman filter for mobile robot attitude estimation [14], is adopted for dual filtering with appropriate assumptions to decrease the number of tuning parameters. However, this tuning process should be performed thoroughly because of randomness in the PSO algorithm. This study thus aims to improve estimation of dual filtering by performing PSO again with different boundary conditions. In numerical investigations, the estimation results by the proposed approach will be presented along with several considerations including the use of iterative tuning process.

This paper first introduces the algorithm of UKF with its general formulation for dynamic state-space model in Section 2. Various applications of filter methods will be also mentioned together with practical issues. Then, Section 3 introduces the problem statement of dynamic system, which can be damaged during seismic excitations. A proper state-space model is proposed and explored to choose appropriate adaptive filtering techniques. In Section 4, for an actual earthquake excitation, several damage scenarios are assumed based on the calculated time histories of relative displacement between adjacent stories. For these scenarios, accuracy and robustness of the proposed method are tested and demonstrated in comparisons with normal UKF and other adaptive filtering techniques. Lastly, this paper concludes with several insights and future research topics.

2. Background: Identification of Dynamic Systems by Unscented Kalman Filter

2.1. A Brief Review of Unscented Kalman Filter

To overcome the aforementioned challenges in using EKF for estimation of nonlinear state-space models, Julier and Uhlmann (2004) developed Unscented Kalman Filter (UKF) [12], also known as Sigma-Point Kalman Filter (SPKF). The unscented transformation of UKF is based on a statistical technique using sigma points (SPs), which provides derivative-free state estimation for a nonlinear state-space model. UKF is known to achieve the second-order accuracy with the same computational complexity as that of EKF.

For details of UKF, let us consider a dynamic state-space model governed by a set of discrete-time nonlinear equations [19]

$$\begin{gathered} x_{k+1}=f\left(x_k\right)+w_k \\ {\mathit{z}}_{k+1}=h\left(x_{k+1}, {\mathit{u}}_{k+1}\right)+v_{k+1} \end{gathered}$$

(1)

where $k$ is the time index; $x_k\in {\mathbb{R}}^{n_x\times 1}$ is the state vector representing displacement, velocity, or system parameters such as elasticity and damping coefficient at time $t_k$; $f(·)$ and $h(·)$ respectively denote nonlinear system-process model and measurement-output model; $w_k\in {\mathbb{R}}^{n_x\times 1}$ represents the artificial process noise modeled by a Gaussian random vector with zero means and covariance matrix $Q_k\in {\mathbb{R}}^{n_x\times n_x}$, i.e., $w_k~N\left(\mathbf{0},Q_k\right)$; $z_k\in {\mathbb{R}}^{n_z\times 1}$ is the vector of the measured quantities such as acceleration; ${\mathit{u}}_k$ is the vector of inputs; and $v_k\in {\mathbb{R}}^{n_z\times 1}$ denotes artificial measurement noise modelled as $v_k~N\left(\mathbf{0},R_k\right)$. The noises $w_k$ and $v_k$ are assumed to be uncorrelated with each other.

Figure 1 illustrates the general algorithm of UKF [20]. In the figure, ${\widehat{x}}_{k|k}$ and $P_{k|k}^{x\mathrm{x}}$ respectively denote the mean (column) vector and covariance matrix of the states at the k-th time step ($k=0,1,\dots$) before the updates are performed. At each time step, sigma points (SPs) ${\chi }_{k|k}^{\left[i\right]}$, $i=0,\dots ,2n_x$, are first generated based on ${\widehat{x}}_{k|k}$ and $P_{k|k}^{x\mathrm{x}}$ as

$$\begin{gathered} {\mathit{\chi }}_{k|k}^{\left[0\right]}={\widehat{x}}_{k|k} \\ {\mathit{\chi }}_{k|k}^{\left[i\right]}={\widehat{x}}_{k|k}+{\left(\sqrt{\left(n_X+\kappa \right)P_{k|k}^{xx}}\right)}^i, i=1,\dots ,n_X \\ {\mathit{\chi }}_{k|k}^{\left[i\right]}={\widehat{x}}_{k|k}-{\left(\sqrt{\left(n_X+\kappa \right)P_{k|k}^{xx}}\right)}^{i-n_X},i=n_X+1,\dots ,2n_X \end{gathered}$$

(2)

where the square-root symbol denotes the “matrix square root”, which can be computed by Cholesky factorization [20], and ${( )}^i$ indicates the i-th column vector of the matrix. The tuning parameter $\kappa$ is often used to consider high order moments in the approximation and assumed to be zero in this paper. Next, using the system-process model $f(\cdot )$, ${\widehat{x}}_{k|k}$ and $P_{k|k}^{x\mathrm{x}}$ are updated to ${\widehat{x}}_{k+1|k}$ and $P_{k+1|k}^{x\mathrm{x}}$, respectively, via the sigma points in Equation (2). In detail, this “time-update” is performed as

$$\begin{gathered} {\mathit{\chi }}_{k+1|k}^{\left[i\right]}=f\left({\mathit{\chi }}_{k|k}^{\left[i\right]}\right) \\ {\widehat{x}}_{k+1|k}=\sum _{i=0}^{2n_X}{W}_m^i{\mathit{\chi }}_{k+1|k}^{\left[i\right]} \\ {P}_{k+1|k}^{xx}=\sum _{i=0}^{2n_X}{W}_c^i\left[{\mathit{\chi }}_{k+1|k}^{\left[i\right]}-{\widehat{x}}_{k+1|k}\right]{\left[{\mathit{\chi }}_{k+1|k}^{\left[i\right]}-{\widehat{x}}_{k+1|k}\right]}^{\mathrm{T}}+Q_k \end{gathered}$$

(3)

where $W_m^i$ and $W_C^i$ are the weighting coefficients of SPs for which $W_m^0=W_C^0=\kappa /\left(n_X+\kappa \right)$ and $W_m^i=W_C^i=1/2\left(n_X+\kappa \right)$ are used [19]. ${\widehat{x}}_{k+1|k}$ and $P_{k+1|k}^{xx}$ in Equation (3) are considered prior estimates of the states at the k-th time step.

For the “measurement-update” of the prior estimates in Figure 1, SPs are first generated as

$$\begin{gathered} {\mathit{\chi }}_{k+1|k}^{\left[0\right]}={\widehat{x}}_{k+1|k} \\ {\mathit{\chi }}_{k+1|k}^{\left[i\right]}={\widehat{x}}_{k+1|k}+{\left(\sqrt{P_{k+1|k}^{xx}}\right)}^i, i=1,\dots ,n_X \\ {\mathit{\chi }}_{k+1|k}^{\left[i\right]}={\widehat{x}}_{k+1|k}-{\left(\sqrt{P_{k+1|k}^{xx}}\right)}^{i-n_X},i=n_X+1,\dots ,2n_X \end{gathered}$$

(4)

Then, the measurement-updates are performed to obtain the following sequentially:

$$\begin{gathered} {\mathit{\zeta }}_{k+1|k}^{\left[i\right]}=h\left({\mathit{\chi }}_{k+1|k}^{\left[i\right]}\right), i=0,\dots ,2n_X \\ {\widehat{z}}_{k+1|k}={\displaystyle \sum }_{i=0}^{2n_X}{W}_m^i{\zeta }_{k+1|k}^{\left[i\right]} \\ {\mathit{v}}_{k+1}=z_{k+1}-{\widehat{z}}_{k+1|k} \\ {P}_{k+1|k}^{zz}={\displaystyle \sum }_{i=0}^{2n_X}{W}_C^i\left[{\zeta }_{k+1|k}^{\left[i\right]}-{\widehat{z}}_{k+1|k}\right]{\left[{\zeta }_{k+1|k}^{\left[i\right]}-{\widehat{z}}_{k+1|k}\right]}^{\mathrm{T}}+R_{k+1} \\ {P}_{k+1|k}^{xz}={\displaystyle \sum }_{i=0}^{2n_X}{W}_C^i\left[{\mathit{\chi }}_{k+1|k}^{\left[i\right]}-{\widehat{x}}_{k+1|k}\right]{\left[{\zeta }_{k+1|k}^{\left[i\right]}-{\widehat{z}}_{k+1|k}\right]}^{\mathrm{T}} \\ K=P_{k+1|k}^{xz}{\left(P_{k+1|k}^{zz}\right)}^{-1} \end{gathered}$$

(5)

where ${\zeta }_{k+1|k}^{\left[i\right]}$ is an SP calculated by substituting ${\mathit{\chi }}_{k+1|k}^{\left[i\right]}$ to the function $h(·)$; $v_k$ is an innovation vector, i.e., theoretical residuals; $P_{k+1|k}^{zz}$ and $P_{k+1|k}^{xz}$ are covariance of the measurements and cross covariance of state and measurements; and $K$ is the calculated Kalman gain.

Finally, the results of the two updates are used to obtain the posterior estimates of the states and covariance matrix by the “state-update” step, i.e.,

$$\begin{gathered} {\widehat{x}}_{k+1|k+1}={\widehat{x}}_{k+1|k}+K{\mathit{v}}_{k+1} \\ {P}_{k+1|k+1}^{xx}=P_{k+1|k}^{xx}-K{P}_{k+1|k}^{zz}{K}^{\mathrm{T}} \end{gathered}$$

(6)

These time-, measurement- and state-updates are performed at each time step $k$ using the results of the previous step as illustrated in Figure 1.

2.2. Applications of Filter Methods and Practical Issues

In applications of filter methods to structural systems, it is generally assumed that the accelerations can be measured. For the purpose of state estimation, parameter estimation, or input estimation, the corresponding state vector is respectively defined in terms of the dynamic responses (e.g., displacement and velocity), system parameters (e.g., Young’s modulus, damping coefficient, and parameters of Bouc-Wen hysteretic models), or input forces. If system parameters are estimated concurrently with dynamic responses, the filter is termed as joint state-parameter estimation. In the so-called joint input-state estimation, input forces are estimated together with the responses.

For example, as a joint state-parameter estimation, Yang et al. (2006) proposed an adaptive EKF-based tracking technique to detect sudden damage in a five-story linear structure [11]. Wu and Smyth (2007) applied UKF to a nonlinear system having 2 degrees of freedom (DOFs) for the purpose of SI to find that the UKF provided better estimation than EKF for highly nonlinear systems [21]. Chatzi and Smyth (2009) also utilized UKF for SI of nonlinear 3-DOF structural system assuming non-collocated heterogeneous sensing [22]. Most previous studies on joint state-parameter estimation have focused on improving performance of filters in tracking parameters of highly nonlinear structural system. These studies employed the system equation of motion under the assumption that either acceleration or displacement is measured at every DOF to avoid rank-deficiency. However, measuring dynamic responses at all DOFs is infeasible especially for structures consisting of many components. This may limit the applicability of joint state-parameter estimation to SI in practice.

Recently, many research efforts have been made regarding the joint input-state estimation to track the input force whose exact information is not available in many situations. Gillijns and De Moor (2007) proposed the recursive-three-step filter (RTSF), which updates responses and input force at the different step, for unbiased minimum-variance estimation of input and state [23]. Since RTSF suffered from rank-deficiency caused by a limited number of measured acceleration or reduced-order model from modal analysis, Lourens et al. (2012) proposed an advanced version of RTSF by using truncated covariance matrices through singular value decomposition [24]. Azam et al. (2015) developed a dual Kalman filter consisting of two stages to prevent numerical issues such as un-observability and rank-deficiency of RTSF when using sparse acceleration measurements [25]. Aucejo et al. (2019) proposed an augmented Kalman filter, whose state vector consists of responses and input together, with the state-space model developed based on the generalized-α method [18]. This was to address the issue that the dual Kalman filter is sensitive to selection of covariance matrices of artificial noises. These studies focused on joint input-state estimation of linear structural systems with time-invariant system parameters. Although full-measurements are not required because of the linear model assumption, rank-deficiency and divergence are still critical issues when SI relies on sparse measurements or a reduced-order model is used.

Meanwhile, as a parameter estimation example, Astroza et al. (2019) applied UKF to Finite-Element (FE) updating with adaptive techniques to consider uncertainty of nonlinear FE models [13]. The state-space model based on Newmark method was able to utilize sparse measurements whereas the use of equilibrium equations in state-space models of previous joint state-parameter estimations needs full measurements. Then, dual adaptive filtering was performed based on the proposed state-space model to consider uncertainty in FE model. The applications using three different approaches, i.e., joint state-parameter estimation, input-state estimation, and parameter estimation confirmed that the state-space model must be developed appropriately based on target estimations, and it is essential to use a proper filtering method to obtain optimal estimations.

It is also noted that, in all filter-related research, it is critical to tune the filter by selecting covariance matrices $Q_k$ and $R_k$ which describe uncertainties of the artificial noises in system-process and measurement-output model respectively. Although the tuning process of those covariance matrices is vital for good filter performance, there is no general technique available. In efforts to handle this, many studies assumed that the artificial measurement noise can be prescribed using the available specifications of the sensors. Then, the artificial process noise is determined by minimizing the root mean square (RMS) of time-history “innovations,” which are defined as follows, to quantify the differences between measurements and responses corresponding to the system parameters predicted during the time-update part [13,16]:

$$f_I\left(Q_k\right)=\sqrt{\frac{1}{n_T}\left({\sum }_{k=0}^T{\Vert {\mathit{v}}_k\left(Q_k\right)\Vert }^2\right)}$$

(7)

where ${\mathit{v}}_k\left(Q_k\right)$ is an innovation vector defined for the time duration $k=0,\dots ,T$; and $n_T=T+1$. In several references [25,26], it has been shown that $f_I\left(Q_k\right)$ in Equation (7) tends to decrease following an L-shaped curve as $Q_k$ increases. Generally, the inflection point of the L-shaped curve is recommended as a near-optimal solution in terms of the filter performance. The opposite case, in which the artificial process noise is known, and the corresponding artificial measurement noise is to be determined, can be solved in a similar manner. If both types of artificial noises are unknown, the objective function can feature both $Q_k$ and $R_k$. To address this, for instance, Odry et al. (2018) utilized particle swarm optimization (PSO) [14]. In this paper, PSO is used to obtain the best tuning parameters of the proposed dual adaptive filtering.

3. Proposed Method: Adaptive Unscented Kalman Filter Using Sparse Measurements of Accelerations

3.1. Regularization-Based State-Space Model

This paper aims to develop a near-real-time system identification method for structures under disastrous dynamic loads such as earthquake excitations. In particular, the method focuses on detecting sudden changes in the system parameters such as Young’s modulus caused by structural damage during excitations. It is assumed that input data such as an earthquake (EQ) ground motion has been acquired accurately. Conventional filter methods such as EKF and UKF show limitations in tracking the status-history having discontinuity points because those methods using constant artificial noises cannot adapt to time-varying statistics in the measurements. To address this issue, Yang et al. (2006) proposed an adaptive joint state-parameter estimation method requiring full measurement condition, which is however not available in general practice [11]. The adaptive rule was tested with simple events where just one component experiences sudden damage. On the other hand, state-space models in joint input-state estimations do not accord with the objective of this paper because those models are constructed based on modal analysis for which the system parameters are assumed to be time-invariant.

To facilitate SI under practical measurement conditions, we aim to develop a filter that performs well even with sparse measurements. To this end, let us first consider the following state-space model proposed by Astroza et al. (2019) for the purpose of parameter estimation [13]:

$$x_{k+1}=x_k+w_k$$

(8)

$${\mathit{z}}_{k+1}=h\left(x_{k+1}, {\mathit{u}}_0, {\dot{\mathit{u}}}_0, {\mathit{p}}_{1:k+1}\right)+v_{k+1}$$

(9)

where $x_k$ is the state vector representing parameters to estimate, e.g., Young’s moduli at time $t_k$; ${\mathit{u}}_0$ and ${\dot{\mathit{u}}}_0$ are the initial displacement and velocity respectively (hereafter, the dot means time-derivative); and ${\mathit{p}}_{1:k+1}$ denotes the time-history of excitation up to $t_{k+1}$, e.g., seismic loads. The function $h(·)$ represents a response function from an FE model, e.g., the Newmark method ($\gamma =1/2, \beta =1/4$). The system-process model in Equation (8) is linear, thus Equation (3) can be simplified into ${\widehat{x}}_{k+1|k}=f\left({\widehat{x}}_{k|k}\right)$ and $P_{k+1|k}^{x\mathrm{x}}=P_{k|k}^{xx}+Q_k$ in this case, i.e., without deriving SPs as in Equation (2). The measurement-output model in Equation (9) predicts responses, which will be compared with the measured response. In this state-space model, full-measurement is not required because the measurement-output model adopts the response function rather than the equilibrium equation.

It is noted that the measurement-output model in Equation (9) was originally designed for time-invariant system parameters. Consequently, a transient period can appear as the estimations of conventional filter converge into the changed status of damaged components. If that transient period is considerably long, it may hamper accurate estimations. To address this limitation of the state-space model of Equations (8) and (9) in SI of time-variant systems, Equation (9) is modified as follows to track sudden changes:

$${\tilde{q}}_{k+1}=h\left(x_{k+1}, {\tilde{q}}_k, {\mathit{p}}_{k+1}\right)$$

(10)

$${\mathit{z}}_{k+1}=S\mathsf{\Phi }\left(x_{k+1}\right){\tilde{q}}_{k+1}^{\mathrm{T}}+{\mathit{v}}_{k+1}$$

(11)

where ${\tilde{q}}_k$ is the matrix consisting of modal displacement, velocity, and acceleration, i.e., ${\tilde{q}}_k=\left[{\mathit{q}}_k {\dot{\mathit{q}}}_k {\ddot{\mathit{q}}}_k\right]$; ${\tilde{q}}_0=\Phi {\left(x_0\right)}^{-1}\left[{\mathit{u}}_0 {\dot{\mathit{u}}}_0 {\ddot{\mathit{u}}}_0\right]$; $\Phi \in {\mathbb{R}}^{n_{\mathrm{DOF}}\times n_{\mathrm{DOF}}}$ is the matrix of the mode shapes consisting of the eigenvectors ${\Phi }_j$, $j=1,\dots ,n_{DOF}$ (calculated from undamped eigenvalue problem $K\Phi =M\Phi {\mathsf{\Omega }}^2$, where $K$, $M$, and $\mathsf{\Omega }$ are stiffness, dumped mass matrix, and a diagonal matrix of the eigen frequencies, respectively); and $S_{n_{\mathrm{DOF}}\times 3}$ is the selection matrix for the measured responses. In Equation (9), at every time step, the current response such as acceleration is calculated from initial responses, which becomes a computational burden as the amount of calculation increases linearly as time goes on. By contrast, Equation (10) uses the Newmark method to predict the next modal responses based on the current information, i.e., system parameters estimated from time-update part and modal responses of the previous time-step. This modified measurement-output model facilitates handling time-variant system parameters. In Equation (11), the vector of the predicted modal responses is transformed to nodal responses by the matrix of calculated mode shapes. Each element in the selection matrix is set as 1 if measurements are made at the corresponding DOF. Finally, the nodal responses predicted by the changed system parameters are matched to the measured responses to calculate the Kalman gain.

The state-space model consisting of Equations (8), (10) and (11) may achieve desirable accuracy of estimation in some cases. However, when using sparse measurements, it has been reported that even if the predicted response well matches the measured responses, unobserved response quantities may be considerably different from the actual ones due to the biased estimation of parameters [27]. Noise in sparse measurement may exacerbate the bias and instability of estimation. To alleviate this undesirable situation, this paper proposes to utilize a regularization technique to avoid rapid and wrong transitions during the estimation. In particular, adopting a well-known regularization technique, Tikhonov term is included in the objective function with a regularization factor to obtain stabilized solutions based on a priori information and the measurements [28,29,30,31,32]. From the Bayesian perspective, the Tikhonov term corresponds to a prior (as a normal distribution in general) to inform a possible solution space [1]. In detail, a regularization term is added to measurement-output model as constraints as follows:

$$\begin{gathered} x_{k+1}={\widehat{x}}^P+r_x \\ \mathrm{where} {\widehat{x}}^P={\sum }_{i=k-N_w}^k{\widehat{x}}_i \end{gathered}$$

(12)

where ${\widehat{x}}^P$ is a vector representing prior information about system parameters, which is updated at every time-step in Equation (12); $r_x$ denotes the standard deviation of prior as regularization factor to constrain the estimation, i.e., $x_k~N\left({\widehat{x}}^P, \Sigma \right)$ and $diag\left(\Sigma \right)=r_x^2$; and $N_w$ is set up as 5% of $n_T$. Actually, it is a critical issue to determine the regularization factor in existing studies of regularization. In this paper, the regularization factor in Equation (12) is selected intuitively not to increase the number of tuning parameters. Because the following dual adaptive filtering techniques requires more than two tuning parameters, it is not desirable to further increase the number of those parameters. The assumption about intuitive selection will be explained in the part regarding numerical investigations.

By including the regularization term in the measurement-output model, stable estimation is achieved while avoiding divergence due to noisy sparse measurements. Therefore, this regularization-based state-space model is utilized in this study to track sudden damage in system parameters under seismic loads.

3.2. Dual Adaptive Filtering for Measurement Noise Estimation

For accurate estimation of time-variant system parameters, the advanced state-space model needs to be used in conjunction with a proper adaptive filtering technique. There are two reasons for the limited tracking ability of conventional filters when statistics in measurements change. First, the predicted state vector in the time-update part is influenced by the size of covariance of artificial process noise. Second, both covariance matrices of artificial process and measurement noises affect the calculation of the Kalman gain. Therefore, if both covariance matrices are constant during estimation, the tracking performance of conventional filters will be limited for time-variant dynamic system.

While developing the proposed filter, a few other adaptive techniques were explored. For example, Liu et al. (2005) proposed an adaptive tracking technique for the purpose of leak detection and location in the management of a pipeline system [17]. This tracking technique was performed based on particle filter, which is also known as a sequential Monte Carlo method. However, it is noted that the idea of adaptive technique is utilized also in the UKF-based approaches. The details and limitations of this technique (with modification for UKF) are described in Appendix A. It is especially noted that this adaptive technique controls the artificial process noise, which causes severe bias and instability of estimations. The other adaptive technique explored in this study is a master-slave Kalman filter (MS-KF), which is also called as a dual adaptive filtering. The main concept of MS-KF is that, when a master filter (MF) estimates the given random variables, the slave filter (SF) will estimate the covariance matrix of either artificial process or measurement noises in MF at every time-step. There are two types of MS-KF based on process noise estimation (MS-KF_PE) and measurement noise estimation (MS-KF_ME). Among these, details of MS-KF_PE are briefly presented in Appendix B because MS-KF_PE is not suitable for the goal of this study for the same reason as the adaptive technique described in Appendix A.

Unlike the two techniques mentioned above, MS-KF_ME controls artificial measurements noise and can be utilized with the advanced state-space model. The general algorithm of MS-KF is illustrated by a block diagram in Figure 2. In the slave filter (SF) part, two random variables ${\omega }_k\in {\mathbb{R}}^{n_x\times 1}$ and ${\upsilon }_k\in {\mathbb{R}}^{n_z\times 1}$ represent the artificial noises for process and measurements respectively modeled as ${\omega }_k~N(0,T)$ and ${\upsilon }_k~N(0,U)$ in which $T\in {\mathbb{R}}^{n_x\times n_x}$ and $U\in {\mathbb{R}}^{n_z\times n_z}$ are time-invariant covariance matrices; and the superscript S indicates a variable related to the SF. A well-known drawback of MS-KF is its sensitivity with respect to the choice of two covariance matrices of artificial noises in SF. Therefore, these artificial noises should be selected carefully with those of MF [15,27]. Linear KF (LKF) is utilized for SF in most cases because there is no accurate knowledge available for the state-space model [15]. In this paper, a covariance-matching technique is utilized in the state-space model of SF because this intuitive technique can be implemented easily with low computational costs [13,15,27,33].

In Figure 2, using the equations related to the number (2) indicates that the MS-KF estimates the covariance of artificial measurements noise at every time-step of MF [13,15] (if the equations noted by the number (1) are used, MS-KF_PE is performed). MS-KF_ME can be conducted based on two types: “the innovation and its covariance” or “the residual and its covariance”. In this paper, the latter one is selected because its calculations are less complex.

The goal of SF in MS-KF_ME is to make the innovations, which is employed to provide a measure of optimality of the filter, compatible with their expected covariance matrix. For this purpose, the covariance-matching technique is conducted, i.e.,

$$z_{k+1}^s=diag\left(\mathrm{E}\left[{\mathit{v}}_{k+1}{\mathit{v}}_{k+1}{}^{\mathrm{T}}\right]\right)$$

(13)

$$I_{k+1}=diag\left(\sum _{i=0}^{2n_x}{W}_c^i\left[{\zeta }_{k+1|k}^i-{\widehat{\mathit{z}}}_{k+1|k}\right]{\left[{\zeta }_{k+1|k}^i-{\widehat{\mathit{z}}}_{k+1|k}\right]}^{\mathrm{T}}\right)$$

(14)

where Equation (13) represents the actual covariance based on innovations, and the Equation (14) is the predicted covariance of innovations from MF. If two matrices from Equation (13) and (14) are similar to each other, the current covariance of artificial measurement noise of MF is considered appropriate. If the calculated value of Equation (13) is larger than that of Equation (14), $R_{k+1}$ needs to be increased to avoid divergence of MF, and consequently this procedure increases the covariance of states, $P_{k+1}^{xx}$.

There are three advantages in using MS-KF_ME: First, no additional parameter is introduced for the algorithm (as for MS-KF_PE, the time-window is required). Second, while the other two adaptive techniques shift the solution space rapidly by increasing $Q_k$, MS-KF_ME widens the solution space by increasing $P_{k+1}^{xx}$. Indeed, the biased estimation may be avoided by increasing $P_{k+1}^{xx}$ instead of increasing the size of $Q_k$ rapidly, and stable estimation will be achieved. Lastly, this correction process for artificial measurements noise is implemented in measurement-update of MF before calculating the Kalman gain. Therefore, the corrected covariance of artificial measurements noise by SF can be utilized immediately. These advantages enable MS-KF_ME to track the sudden decrease of system parameter in a stable manner during the seismic excitations.

The algorithm of MS-KF_ME with the newly proposed state-space model is summarized in Figure 3. In the time-update part, calculations of SPs are not needed because of the linear system-process model in Equation (8). The part corresponding to the slave LKF is represented in gray background along with detailed descriptions. It is noteworthy that the regularization term in Equation (12) is considered in the form of dummy variables [34] when the auto- and cross- covariance matrices are calculated in the measurement-update part of MF. In these calculations for two covariance matrices, the prior information of regularization is required. When the time-step is shorter than the selected time-window, the values of initial state are used for this prior information. After that time point, this prior is updated in the state-update part.

4. Numerical Investigations

4.1. SI under Simulated Damage Scenarios

To test the combination of the selected dual adaptive filtering techniques and the proposed state-space model based on regularization, a shear building having 6 DOFs in Figure 4 is investigated. The lumped mass of each floor is assumed as 1 kg, and Young’s modulus of all inter-story is equal to 1 Pa (i.e., inter-story stiffness of each floor is 24 N/m). Rayleigh damping, which is proportional to mass and stiffness matrices, is utilized by assuming a damping ratio of 5% for the first and second modes. For the linear system, Young’s modulus of each component is estimated using an EQ ground motion. The time-history of the 180-degree component of ground acceleration recorded at El Centro station during the earthquake event in 1940 with M_w = 6.95 is used in this numerical investigation. The length of this ground acceleration history is 42 s, and discretized by 0.01 s. To consider potential damage scenarios, the EQ time-history is first applied to the shear building without damage events. The relative displacements between adjacent stories are calculated to find the peak values at each story and the occurrence-time as shown in

Table 1. Calculated peak relative displacements and occurrence-times without considering damage.

	Peak Relative Displacement	Occurrence-Time of Peak Response
1st story	0.0559	8.58 s
2nd story	0.0454	4.94 s
3rd story	0.0393	5.23 s
4th story	0.0500	7.50 s
5th story	0.0568	9.92 s
6th story	0.0490	6.15 s

.

Based on the simulated sequence of peak displacements, five damage scenarios are created in terms of sudden decrease of Young’s moduli as summarized in

Table 2. Damage scenarios based on simulated peak relative displacements in Table 1.

Scenario No.	Description of Scenario
Scenario 1	E2 decreases by 25% at 4.94 s ➔ E1 decreases by 33% at 8.58 s
Scenario 2	E4 at 7.50 s (25%) ➔ E1 at 8.58 s (33%)
Scenario 3	E5 at 9.92 s (28%) ➔ E2 at 14.89 s (24%)
Scenario 4	E3 at 5.23 s (24%) ➔ E2 at 8.58 s (24%) ➔ E2 at 16.00 s (41%)
Scenario 5	E6 at 6.15 s (25%) ➔ E4 at 7.50 s (22%) ➔ E1 at 8.58 s (33%)

. Scenarios 1 and 2 depict sequence of sudden decreases of story stiffness. Scenario 3 simulates a sequential failure initiated by the maximum-peak response at 5th story. Scenarios 4 and 5 are assumed to depict more complicated sequential scenarios. In Scenario 4, the second story experiences the additional sudden decrease 7.42 s after the initial decrease of corresponding story’s stiffness. Scenario 5 consists of three sudden decreases of different stories. After simulating the above scenarios, the condition of sparse measurements is defined such that accelerations at the first, second, fourth, and sixth story can be measured. From this assumption, the selection matrix S will be set up as $S_{13}=S_{23}=S_{43}=S_{63}=1$ and $S_{ij}=0$ otherwise. Then, the noise-history, which is simulated from a uniform distribution between 95% and 105% of RMS of measured accelerations, is added to simulate noise in the measurements.

4.2. Tuning Process Using Particle Swarm Optimization

As stated in Section 2.2, all filter methods should be tuned by selecting proper covariance matrices of artificial noises to achieve an optimal estimation. In using LKF, EKF, and UKF, just one parameter is tuned. Although each diagonal term of those covariance matrices may require individual tunings, it is assumed that all elements in vectors of states and measurements have same statistical characteristics for simplicity. Generally, the covariance of artificial measurements noise is determined based on sensor specification. The last tuning parameter is then the covariance of artificial process noise and can be selected from the derived L-curve. However, this process does not work for adaptive filtering techniques, because more than two tuning parameters are required even under the assumption that sensor specification can determine the covariance of the artificial measurement noise. The tuning process for MS-KF_ME with the proposed state-space model handles three parameters, i.e., two parameters corresponding to covariance matrices of master and slave filters’ artificial process noise, and the regularization factor in the proposed state-space model. In this case, it is impossible to derive the L-curve for two tuning parameters.

This paper utilizes the particle swarm optimization (PSO) to tune the dual adaptive technique with the proposed state-space model. PSO was originally inspired by the information circulation and social behavior observed in bird flocks and fish schools and is suited to solve problems where the optimal solution is a point in a multi-dimensional space of the parameter [35]. Odry et al. (2018) showed that PSO can find proper tuning parameters for optimal estimation performance [14]. However, if there are many variables in global optimization problem, this optimization process may give a computational burden. This is because PSO uses many particles to probe the best in the solution space, and the required number of particles increases exponentially as the number of tuning parameters increases. Another issue is that the performance of MS-KF is sensitive to choices of covariance matrices of artificial noises for master and slave filter. Actually, the positive definiteness condition for calculating SPs related with Equations (2) and (4) may be violated by subtle changes in those covariance matrices, and this numerical issue will interrupt the procedure of UKF as master filter. Therefore, it is desirable to decrease the number of tuning parameters in PSO for stable tuning process, and thus, the following assumptions will be introduced.

First, the regularization factor in the proposed state-space model is assumed as 10% of initial states, which means that the solution space will be constrained by the normal probability density, i.e., $x~\mathit{N}({\widehat{x}}^P,\Sigma )$ and $diag\left(\Sigma \right)={\left(0.1·|x_0|\right)}^2$. The reason for adopting the regularization was to support sequential estimations of filter to avoid divergence caused by the lack of information due to sparse measurements. Therefore, there is no need to select the accurate regularization factor necessarily because the regularization term can confine solution space weakly. In the following numerical examples, moderate damage, which cannot be inspected visually, is assumed. Thus, in this paper, 10% of initial states’ values is used as the regularization factor, which gave rapid convergence to the exact value. Other values may be also appropriate, but too small values such as 1% of initial values are likely to result in long transient period before converging to an accurate solution. Second, each diagonal term of covariance of artificial measurements noise is assumed as 5% of the RMS of measurements for master filter. This assumption of artificial measurements noise can be improved by accurate information about actual sensor specification. Similarly, the covariance of artificial measurement noise for slave filter is determined as 1% of the identity matrix for simplicity.

Through the above assumptions, two tuning parameters that control the size of covariance matrices of artificial process noise for master and slave filter are used as variables for PSO. As in Equation (7), these two parameters affect the innovation in the following objective function:

$$\begin{gathered} f_I\left(x_1,x_2\right)=\sqrt{\frac{1}{n_T}\left({\sum }_{k=0}^T{\Vert {\mathit{v}}_k\left(Q(x_1\right),T(x_2))\Vert }^2\right)} \\ \mathrm{where} Q=diag\left(|x_0|·{10}^{-x_1}\right) \mathrm{and} T=I_{n_x\times n_x}·{10}^{-x_2} \end{gathered}$$

(15)

in which $x_1$ and $x_2$ are tuning parameters, and $I_{n_x\times n_x}$ is an identity matrix. The objective function in Equation (15) is constructed to find an optimal estimation having the least RMS of innovation-history. This paper used the PSO algorithm in the global optimization toolbox of 2018b-version of MATLAB^®.

When performing PSO, there are two important points to achieve efficiency: (1) setting boundary conditions in PSO for random variables and (2) determining the number of particles representing initial locations. The boundary conditions constrain the solution space, and thus a wrong setting may lead to finding no optimal solution. In some cases, the results of PSO are located near boundaries. In this situation, PSO should be re-performed, and it is recommended to shift the relevant variables’ boundary to include the corresponding end. On the other hand, to determine the number of particles properly, the generated particles should be scattered uniformly around the whole solution space. For example, one could discretize uniformly in the range of each variable such as [0.1, 0.2, 0.3, …, 0.9]. In this case, the total number of particles is assumed as 81 for 2-dimentional problem.

There may be more efficient alternatives than using these recommendations. However, it is essential to perform PSO iteratively to find a near-optimal solution because of intrinsic randomness in PSO.

Table 3. Tuning parameters identified by particle swarm optimization.

Scenario	Results of PSO $\left({\mathit{x}}_1,{\mathit{x}}_2\right)$	Initial Boundary		Changed Boundary
		for x1	for x2	for x1	for x2
Scenario 1	(3.51, 2.33)	[2, 4]	[1, 5]	-	-
Scenario 2	(3.81, 2.19)	[2, 4]	[1, 5]	-	-
Scenario 3	(3.98, 5) ➔ (3.88, 6.99)	[2, 4]	[1, 5]	-	[3, 7]
Scenario 4	(4, 1.98) ➔ (4.21, 1.73)	[2, 4]	[1, 5]	[3, 5]	-
Scenario 5	(3.95, 2.68)	[2, 4]	[1, 5]

presents results of PSO using 81 particles about the damage scenarios with the selected boundary conditions given as [lower bound, upper bound]. In Scenarios 3 and 4, PSO was re-performed because the results were located at one of the boundaries.

With the tuning parameters obtained by PSO, the MS-KF_ME is conducted for the scenarios. For normal-UKF (n-UKF) and adaptive tracking technique-based UKF (AUKF) in Appendix A, PSO using Equation (15) (with the variable x₂ excluded) is performed to determine x₁ as 3.23 and 3.88, respectively. MS-KF_PE (in Appendix B) was not performed due to the difficulty in selecting the size of time-window. Figure 5a–c shows the estimation results by the proposed method, n-UKF and AUKF, respectively, under Scenario 1. From the results, it is confirmed that, while the results of n-UKF and AUKF show biased estimations, the MS-KF_ME derives accurate estimations about all parameters. Although divergence in n-UKF and AUKF was avoided by adopting the regularization term into state-space model, those filters still converged into wrong values. These results show that filter methods, which control artificial process noise, cannot adapt to measurements, whose statistics change due to damage.

Figure 6 reports estimation results by MS-KF_ME for Scenario 2. After a short range of unstable estimation at around 11 s, the estimation converges to accurate values without bias. When initial boundaries by PSO (Table 3) are used for Scenario 3 (Figure 7a), the estimation of the fifth story’s parameter, where acceleration-history is missing, shows divergence after around 38 s. The estimations of first and third stories’ parameters are inaccurate and biased as well. Figure 7b shows that, when the newly derived PSO’s result is used, the bias in estimation of fifth story’s parameter disappears, and the accuracy of first and third stories’ estimations is improved. From these results, it is noted that the tuning process using PSO should be performed iteratively for finding the near-optimal solution.

Next, more complicated sequential damage scenarios are considered to test stable tracking-capability of the proposed method. The PSO’s result for Scenario 4 (the second story experiences the additional sudden damage after the initial damage) was derived at one of the boundaries. Thus, PSO was re-performed to find new tuning parameters reported in Table 3. The estimations with the newly derived tuning parameters are presented in Figure 8. In this figure, the first story’s parameter is not assumed to be damaged. By the way, the range from around 10 to 25 s shows slightly biased estimation of the first story, and this bias is being severe before the estimation approaches to the exact value at around 35 s. This is probably because the corresponding acceleration for the damaged third story is not measured, and the damage occurred at adjacent stories. It may be a burden for slave filter to infer accurate information based on only innovations from master filter without structural knowledge. However, as mentioned in Section 3.2, the MS-KF_ME aims to obtain stable estimations from sparse measurements by searching solution space gradually rather than rapid tracking. The results confirm that the estimations about adjacent stories are distinguishable and converged to accurate values respectively after the transient period.

Finally, Scenario 5, which consists of three sudden damages in three different stories, are investigated as shown in Figure 9. When compared with Scenario 4, this scenario assumes that the damage occurs at separate stories, and corresponding accelerations of these stories are measured. In this case, faster convergence to accurate values is achieved.

4.3. Further Investigations Using Different Ground Motions

To test general applicability of the proposed method, another real EQ ground motion is applied to the same shear building in Figure 4. The time-history of the 65-degree component of ground acceleration from the 1971 San Fernando earthquake (M_w = 6.61 at the ‘2516 Via Tejon PV’ station) is used. The 56-second-long time history is discretized by 0.005 s. This EQ ground motion matches an assumption that damage happens at the middle of EQ time-history. That is, when calculating the relative displacements between adjacent stories, the peak values were located around the middle of EQ time-history. In

Table 4. Damage scenarios for San Fernando EQ and tuning parameters from PSO.

Scenario	Description of Scenario	Results of PSO
Scenario 1’	E2 decreases by 25% at 21.78 s ➔ E1 decreases by 33% at 25.45	(4.12, 8.41)
Scenario 2’	E4 at 22.45 s (25%) ➔ E3 at 22.94 s (25%) ➔ E1 at 25.45 s (33%)	(3.27, 9.05)
Scenario 3’	E5 at 22.73 s (18%) ➔ E3 at 22.94 s (25%) ➔ E1 at 25.45 s (33%)	(3.40, 8.50)
Scenario 4’	E5 at 22.73 s (19%) ➔ E5 at 45.47 s (38%)	(3.67, 8.53)

, the scenarios based on the calculated peak values are presented. In the previous numerical tests, the n-UKF and AUKF showed early convergence before approaching accurate values. Scenarios in this section require that MS-KF_ME should not be converged too early during the filtering-process and adapt to time-varying statistics of the measurements. From this further investigation, it is proved that the proposed method enables us to estimate status of system parameters without early convergence.

To derive tuning parameters in Equation (15), PSO using 81 particles was performed with boundary conditions as $x_1\in \left[3, 5\right]$ and $x_2\in \left[5, 10\right]$. These boundary conditions, which are different from those in Table 3, are obtained by performing PSO several times. This result proves that PSO requires boundary conditions appropriate for the given problem, and PSO-based tuning process should be re-performed iteratively to obtain proper boundary conditions.

The results in Figure 10, Figure 11, Figure 12 and Figure 13 confirm that MS-KF_ME with parameters tuned by PSO can accurately detect the sudden damage of systems despite sparse and noisy measurements. These also demonstrate that the proposed method helps avoid early convergence and enables us to track status of system accurately throughout the excitation period. In Figure 11 about Scenario 2’, it is noted that the time-history of the fifth story’s target parameter decreases along with the damaged parameter of the fourth story, but the estimation of the fifth story’s parameter returns back to the accurate value. It should be noted that the proposed method achieves this result even though the acceleration of the fifth story was missing. A similar observation is made in Scenario 3’ where the fifth story is damaged instead of the fourth story. The results of Scenario 3’ in Figure 12 show better separation between the estimations of fourth and fifth stories when compared with the results of Scenario 2’. This phenomenon indicates that when acceleration-history, whose statistics change due to the damage of corresponding parameter, is measured, estimations are less stable than the opposite case. In Figure 13, the estimations in the range after successive damage are more accurate than those in the range corresponding to initial damage.

5. Conclusions and Topics for Future Study

In this paper, a novel methodology was proposed to track sudden damage of structural system under dynamic excitations using sparse measurements. First, the state-space model often used for the purpose of parameter estimation was modified to handle time-variant system parameters. The regularization term was also introduced to avoid divergence of estimation caused by sparse noisy measurements. Second, the master-slave Kalman filter for measurement noise estimation (MS-KF_ME) was utilized to adapt to time-varying status of the target system. The unscented Kalman filter (UKF) as master filter estimates the status of target system parameters while a linear Kalman filter (LKF) as slave filter estimates the covariance of artificial measurement noise in UKF through a covariance-matching technique. Lastly, particle swarm optimization (PSO) was exploited to tune the artificial noises of MS-KF_ME with practical considerations. To test the proposed methodology, sets of structural response data were simulated for a shear-building system subjected to two real earthquake ground motion time histories. Several damage scenarios were introduced based on the simulated data such as peak values of inter-story drifts. For comparison purposes, the normal-UKF and adaptive UKF, which control artificial process noise, were conducted with sparse noisy measurements. Unlike the other two methods, the proposed method provided accurate estimation of sudden damage in system parameters even if the corresponding parameter’s response was missing. These results demonstrated the potential of the proposed method for practical applications.

In general, it is critical to tune filter methods to obtain near-optimal estimation results. Especially, the tuning process of MS-KF_ME is known to be sensitive to choices of covariance matrices of artificial noises for master and slave filter. Even if the number of tuning parameters is reduced to just two parameters through the recommended assumptions, it is impossible to determine the tuning parameters empirically. To this end, this paper utilizes particle swarm optimization (PSO). The PSO-based tuning process is affected by choice of boundary conditions, and thus it is essential to perform PSO repeatedly to find the best tuning parameters.

It is recommended to further explore the use of an efficient global optimization method to avoid assumptions that were introduced to reduce the number of tuning parameters. If so, a tuning process of dual adaptive filter will be more flexible to adapt to practical applications. In addition, the slave filter in this paper used the covariance-matching technique as a state-space model. The matching technique depends on differences between actual and theoretical covariance while structural knowledge is not exploited. Because of this limited capability of covariance-matching technique from SI point of view, some important structural information may not be revealed accurately, which causes slow convergence into accurate values or biased estimations. To infer meaningful information appropriately from sparse noisy measurements, the latest machine-learning techniques and Bayesian methods could provide a suitable state-space model for the slave filter.