Practical Implementation of Recursive Subspace Identification on Seismically Excited Structures with Fixed Window

Abstract

As one of the most catastrophic natural disasters worldwide, earthquakes and their effect on structures are very important to structural health monitoring (SHM), particularly for the ones living around the Pacific ring of fire. In this regard, SHM techniques with real-time or online processing can be used to identify states of structures, track modal parameters, provide a warning message about damage, and help post-earthquake reconnaissance and rehabilitation. For instance, a recursive formulation based on subspace identification (SI) has been demonstrated that it is capable to track system changes. In this study, a recursive subspace identification (RSI) algorithm with a fixed window is proposed to investigate the time-varying dynamic characteristics under seismic excitations. Subsequently, some suggestion is described and discussed for practical implementation. For verifying the proposed algorithm, different datasets from full-scale experiments are applied to examine its applicability. In other words, the practicability of implementing RSI in real-time or online has been developed and examined in this paper.

1. Introduction

Structural health monitoring (SHM) is continuously emerging in the field of civil and structural engineering as it is an efficient way to monitor the current state of infrastructure. A well-developed SHM technique can provide useful information about the existence, location, and severity of damage [1], especially during and after a destructive input. For example, as one of the most catastrophic natural disasters worldwide, earthquakes often caused serious damage to structures in a moment. Therefore, the development of real-time or online processing is essential for SHM systems under seismic excitations.

Vibration-based SHM methods are the most common SHM techniques; they are based on the conceptual idea of a modal domain that changes in modal parameters, such as modal frequencies, damping ratios, and mode shapes, are a result of changes in the physical mass, damping, and stiffness properties of the structure [2]. To track the variation and detect the damage, vibration-based SHM methods need to be continuous and autonomous during seismic events for providing timely information and warnings. In the 1990s, Loh et al. introduced least-square estimation (LSE) to detect the changes in structural parameters under seismic excitations [3,4]. For nonlinear hysteretic structural systems, LSE approaches have been employed to identify cumulative damage [5,6]. Later, sequential non-linear LSE (SNLSE) has been developed based on acceleration responses, which is more practical in the field application [7]. Apart from LSE approaches, Hoshiya and Saito first adopted a filtering approach to iteratively identify structural systems [8]. In 2004, Yang et al. extended and used Kalman filter to adaptively track parameter changes [9,10]. Then, Kalman filter approach was applied to nonlinear systems [11,12] and verified via experiments [7,8]. Moreover, other filtering approaches, such as H_∞ filter [13] and Monte Carlo filter [14], have been proposed to track parameters over time.

During the past few decades, subspace identification (SI) and its relative approaches have attracted substantial interest in system identification, damage detection, and structural control since they can directly identify the system matrices from input and output data. Several well-known SI algorithms, including CVA (canonical variate analysis) [15], N4SID (numerical algorithms for subspace state space system identification) [16], PAST (projection approximation subspace tracking) [17], and so on, have been developed to identify structural systems and detect structural damage. To accommodate for real-time or online processing, various recursive forms have been proposed using either the matrix inversion lemma [18,19,20,21] or Yang’s criterion [22,23,24,25]. One of the major applications in the field of civil and structural engineering is to extract time-varying dynamic characteristics during seismic excitations. Compared to the LSE and filtering approaches, SI and its relative approaches have stable performance owing to the usage of geometry projection and less limitation on measurement because of the explicit input-output relationship. Consequently, SI based on the recursive formulation has been demonstrated to identify states of structures, track modal parameters, provide a warning message about damage, and help post-earthquake reconnaissance and rehabilitation [26,27,28,29,30].

In this study, a recursive subspace identification (RSI) algorithm originally developed in the field of control and automation is introduced to investigate the time-varying dynamic characteristics under seismic excitations. This algorithm is derived with a fixed window and, compared to the above-mentioned algorithms, the window length can be used to indicate extra time information. In Section 2, SI is first introduced to extract the system matrices and the modal parameters from input and output data. Furthermore, the proposed algorithm is carefully derived as a recursive form of SI to facilitate real-time or online processing. Then, some suggestions for practical implementation is described and discussed using numerical simulation in Section 3. To verify SHM applications based on RSI, two full-scale experiments using shaking tables were conducted and the proposed algorithm is applied to identify the modal parameters reclusively in Section 4. Last of all, brief summary and conclusion are drawn in Section 5 and some future research is illustrated.

2. Recursive Formulation

Considering a linear n degrees of freedom (DOFs) building structure subjected to seismic excitations, the equation of motion can be derived as the discrete-time state-space equations as

$${\dot{x}}_{k+1}=A{x}_k+B{u}_k+w_k$$

(1)

$$y_k=C{x}_k+D{u}_k+v_k$$

(2)

where x_k is state vector with 2n states; y_k is measured output vector with m measurement; u_k is input vector with l excitations; w_k and v_k are the process and measurement noise, respectively; the subscript k denotes k-th step which indicates $t=k\Delta t$ and ∆t is the sampling interval of measurement; A is linear elastic system matrix; B and D are excitation influence vector; C is output (or observer) matrix. Among those matrices, A and C are particularly important for the application of SHM and need to be frequently identified using SI.

2.1. Subspace Identification (SI)

To develop SI, the state-space equations need to be transformed into “Matrix Input-Output Equations” [31] as

$$Y_p={\Gamma }_i{X}_p+H_i{U}_p+G_i{W}_p+V_p$$

(3)

$$Y_f={\Gamma }_i{X}_f+H_i{U}_f+G_i{W}_f+V_f$$

(4)

$$X_f=A^i{X}_p+{\Delta }_i{U}_p$$

(5)

where Y_p and Y_f are past and future output data Hankel matrices with i rows and j columns (both in samples), respectively.

$$Y_p=\left[\begin{array}{cccc}{y}_{1|i}& y_{2|i+1}& \cdots & y_{j|i+j-1}\end{array}\right]$$

(6)

$$Y_f=\left[\begin{array}{cccc}{y}_{i+1|2i}& y_{i+2|2i+1}& \cdots & y_{i+j|2i+j-1}\end{array}\right]$$

(7)

$$y_{1|i}={\left[\begin{array}{cccc}{y}_1^T& y_2^T& \cdots & y_i^T\end{array}\right]}^T$$

(8)

Similarly, U_p and U_f are past and future input data Hankel matrices. As a result, the dimension of Y_p and Y_f is $im\times j$ and the dimension of U_p and U_f is $il\times j$. Δ_i is the reversed extended controllability matrix; H_i and G_i are the lower block triangular Toeplitz matrices. Γ_i is extended observability matrix which contains information of system matrices (A and C) and is the primary outcome of subspace identification. Details about those matrices can be also found in the literature [31]. So, the problem of identifying the system matrices is identifying the extended observability matrix.

To extract the key matrix, a geometric tool, called projection in the field of linear algebra, is utilized to obtain the extended observability matrix from the matrix input-output equations. Both orthogonal and oblique projections can be used to extract the extended observability matrix and the “Multivariable Output-Error State sPace (MOESP)” is especially efficient for performing the oblique projection [32]. Instead of a series of matrix multiplication, MOESP implements LQ decomposition to the Hankel matrices as

$$\left[\begin{array}{c}{U}_f\\ {\Xi }_p\\ Y_f\end{array}\right]=\left[\begin{array}{ccc}{L}_{11}& 0& 0\\ L_{21}& L_{22}& 0\\ L_{31}& L_{32}& L_{33}\end{array}\right]\left[\begin{array}{c}{Q}_1^T\\ Q_2^T\\ Q_3^T\end{array}\right]$$

(9)

where Ξ_p denotes an instrumental variable (IV) matrix consisting of both past input and output data Hankel matrices [33].

$${\Xi }_p=\left[\begin{array}{c}{U}_p\\ Y_p\end{array}\right]=\left[\begin{array}{cccc}{\xi }_{1|i}& {\xi }_{2|i+1}& \cdots & {\xi }_{j|i+j-1}\end{array}\right]$$

(10)

It is noted that the left-hand side matrix in Equation (9) is suggested to be a square matrix to secure a correct division after the LQ decomposition, meaning that the number of columns, j, has no choice but

$$j=2i(m+l)$$

(11)

where the dimension of L₃₂ is hence specified, $il\times i(m+l)$. Consequently, $L_{32}{Q}_2^T$ is equal to the oblique projection, O, and the column space of L₃₂ is the same with the column space of the extended observability matrix. It is the essence of why the algorithms are named “subspace” identification.

To retrieve the column space of the extended observability matrix, one can take advantage of singular value decomposition (SVD)

$$O=US{V}^T=\left[\begin{array}{cc}{U}_1& U_2\end{array}\right]\left[\begin{array}{cc}{S}_1& 0\\ 0& S_2\end{array}\right]\left[\begin{array}{c}{V}_1^T\\ V_2^T\end{array}\right]\approx U_1{S}_1{V}_1^T$$

(12)

where $S_2\approx 0$ is null space. The extended observability matrix can be obtained as

$${\Gamma }_i=U_1{S}_1^{1/2}$$

(13)

Furthermore, the discrete-time system matrices (A and C) shown in Equation (1) and Equation (2) can be evaluated as

$$A={\underset{\_}{{\Gamma }_i}}^{†}\overline{{\Gamma }_i}$$

(14)

where $\underset{\_}{{\Gamma }_i}$ denotes Γ_i without the last m rows, $\overline{{\Gamma }_i}$ denotes Γ_i without the first m rows, and the superscript $†$ is Moore-Penrose pseudoinverse. Moreover, C can be determined from the first m rows [31].

Finally, the continuous-time system matrices can be computed from discrete-time ones and the modal parameters, such as modal frequencies, damping ratios, and mode shapes, can be successfully estimated. Those modal parameters are one of the most important dynamic features in the application of SHM.

2.2. Recursive Subspace Identification (RSI)

There are several ways to implement SI recursively and they can be roughly divided into two categories. One uses the matrix inversion lemma to update the orthogonal or oblique projection matrix [18,19,20,21] and the other one uses an unconstrained criterion (known as Yang’s criterion) to approximate the projections [22,23,24,25]. Both of them introduce a forgetting factor, which is quite common in the recursive least squares (RLS) algorithm [34,35], to keep fading out the past measurement and updating the newest information of system matrices.

One of the significant results after introducing a forgetting factor is that the weighting of data is exponentially reduced as time goes on. There are pros and cons to using a forgetting factor; undoubtedly, the implementation is very easy, simple, and direct. On the other hand, the forgetting factors produce a time-varying weighting to the past data, as shown in Figure 1. To overcome the problem, constant weighting is preferable, such as a fixed moving window. Kameyama proposed a recursive formulation based on a swap of matrices, called bona fide recursive algorithm [36,37], and brought a moving window with fixed input and output data size, as shown in Figure 1. Definitely, the identified dynamic behaviors appeared only in the current window compared to the other algorithms using forgetting factors.

The bona fide recursive algorithm is based on MOESP; admittedly, MOESP needs to implement Equation (9) at each step. Given a set of new measurement of input data, u_k+1, and output data, y_k+1, the vectors $u_{k-i+2|k+1}$, $y_{k-i+2|k+1}$, and ${\xi }_{k-i+2|k+1}$ can be updated and the following relation holds

$$\left[\begin{array}{cc}{U}_{f,k}& u_{k-i+2|k+1}\\ {\Xi }_{p,k}& {\xi }_{k-i+2|k+1}\\ Y_{f,k}& y_{k-i+2|k+1}\end{array}\right]=\left[\begin{array}{cc}{u}_{k-i-j+2|k-j+1}& U_{f,k+1}\\ {\xi }_{k-i-j+2|k-j+1}& {\Xi }_{p,k+1}\\ y_{k-i-j+2|k-j+1}& Y_{f,k+1}\end{array}\right]$$

(15)

Regarding $\left[\begin{array}{cc}A& b\end{array}\right]=\left[\begin{array}{cc}c& D\end{array}\right]$, the equation $\left[\begin{array}{cc}A& b\end{array}\right]{\left[\begin{array}{cc}A& b\end{array}\right]}^T=\left[\begin{array}{cc}c& D\end{array}\right]{\left[\begin{array}{cc}c& D\end{array}\right]}^T$ can be held by $A{A}^T+b{b}^T=c{c}^T+D{D}^T$ and Equation (15) can be further expanded [36,37]. For simplicity, the following five covariance matrices are defined

$$K_{0,k}=L_{11,k}{L}_{11,k}^T$$

(16)

$$K_{1,k}=L_{21,k}{L}_{11,k}^T$$

(17)

$$K_{2,k}=L_{31,k}{L}_{11,k}^T$$

(18)

$$K_{3,k}=L_{21,k}{L}_{21,k}^T+L_{22,k}{L}_{22,k}^T$$

(19)

$$K_{4,k}=L_{31,k}{L}_{21,k}^T+L_{32,k}{L}_{22,k}^T$$

(20)

At next step, these covariance matrices are updated as

$$K_{0,k+1}=K_{0,k}+u_{k-i+2|k+1}{u}_{k-i+2|k+1}^T-u_{k-i-j+2|k-j+1}{u}_{k-i-j+2|k-j+1}^T$$

(21)

$$K_{1,k+1}=K_{1,k}+{\xi }_{k-i+2|k+1}{u}_{k-i+2|k+1}^T-{\xi }_{k-i-j+2|k-j+1}{u}_{k-i-j+2|k-j+1}^T$$

(22)

$$K_{2,k+1}=K_{2,k}+y_{k-i+2|k+1}{u}_{k-i+2|k+1}^T-y_{k-i-j+2|k-j+1}{u}_{k-i-j+2|k-j+1}^T$$

(23)

$$K_{3,k+1}=K_{3,k}+{\xi }_{k-i+2|k+1}{\xi }_{k-i+2|k+1}^T-{\xi }_{k-i-j+2|k-j+1}{\xi }_{k-i-j+2|k-j+1}^T$$

(24)

$$K_{4,k+1}=K_{4,k}+y_{k-i+2|k+1}{\xi }_{k-i+2|k+1}^T-y_{k-i-j+2|k-j+1}{\xi }_{k-i-j+2|k-j+1}^T$$

(25)

Obviously, from Equation (21) to Equation (25), the matrices K₀, K₁, K₂, K₃, and K₄ are the covariance of input data, output data, and IV.

However, to get the column space of the extended observability matrix, L₃₂ needs to be extracted from Equation (20)

$$L_{32,k}=(K_{4,k}-L_{31,k}{L}_{21,k}^T){(L_{22,k}^T)}^{-1}$$

(26)

where $L_{31,k}{L}_{21,k}^T$ can be derived from Equations (16)–(18); $L_{22,k}{L}_{22,k}^T$ can be derived from Equations (16), (17) and (19). Therefore, the product of $L_{32,k}{L}_{32,k}^T$ can be constructed as

$$L_{32,k}{L}_{32,k}^T=(K_{4,k}-K_{2,k}{K}_{0,k}^{-1}{K}_{1,k}^T){(K_{3,k}-K_{1,k}{K}_{0,k}^{-1}{K}_{1,k}^T)}^{-1}{(K_{4,k}-K_{2,k}{K}_{0,k}^{-1}{K}_{1,k}^T)}^T$$

(27)

Consequently, the column space of $L_{32,k}{L}_{32,k}^T=O$ is the same with the column space of the extended observability matrix.

Generally, RSI implementation using the bona fide recursive algorithm consists of four steps, as shown in Figure 2. The first step carries out the initial window and directly uses MOESP with LQ decomposition. Moreover, the five covariance matrices are computed using the equations from Equation (16) to Equation (20). The second step constructs $L_{32,k}{L}_{32,k}^T=O$ and, afterward, the extended observability matrix can be retrieved using SVD. In the third step, the discrete-time system matrices can be evaluated using Equation (16), the continuous-time system matrices can be computed using Equation (14), and the modal parameters can therefore be estimated. Finally, the new measurement is acquired for the next step and the covariance matrices are updated with those equations from Equation (21) to Equation (25). Then, repeating from the second step to the fourth step produces a close-loop recursion.

3. Numerical Simulation

Different from the forgetting factor adopted by other algorithms, the proposed algorithm brings a moving window with fixed input and output data size. To demonstrate the effectiveness of the proposed method for seismically excited structures, an 8-story shear-type frame is numerically simulated with a seismic excitation of the El Centro Earthquake (18 May 1940) in the north-south direction. The mass of each story is lumped at the center and assumed to be 10 tons. The stiffness at each story is supposed to be 1200 kN/m. The damping coefficients are calculated via the modal damping and the damping ratio is assumed to be 2% for each mode. Accordingly, the 8 modal frequencies from 1.02 Hz to 10.84 Hz are listed in

Table 1. The 8 modal frequencies of numerically simulated structure.

Mode Number	1st	2nd	3rd	4th
Frequency (Hz)	1.0174	3.0176	4.9150	6.6450
Mode Number	5th	6th	7th	8th
Frequency (Hz)	8.1487	9.3750	10.2820	10.8388

. The measurement is assigned to be the 8 acceleration responses of the frame and an additional measurement is the ground acceleration. Although the frame is simulated using a 1 kHz sampling rate, the one in measurement can be variated. All the simulation is done using the software, MATLAB, issued by MathWorks.

3.1. User-Defined Parameter

Before applying SI or RIS, the first and inevitable question is the size of Hankel matrices, such as U_f, Ξ_p, and Y_f. The only user-defined parameter is i because the left-hand side matrix in Equation (9) is a square matrix. Once i is assigned, j can be determined using Equation (11), and the total data size within the current window, l_w, called window length, is resolved as

$$l_w=2i(m+l+1)-1$$

(28)

In the field application, the window length is very important since the product of the window length and the sampling interval is the time of providing the first identification result after an earthquake starts. If the window length is longer than the duration of the mainshock, RSI can only identify the modal parameters belonging to the damaged structure. The results may be useless without comparing with the original state. By setting i as 10, the modal parameters identified by RSI are shown in Figure 3.

Obviously, RSI implementation using the bona fide recursive algorithm gives an excellent estimation of the modal parameters. The grids in Figure 3b are drawn according to the 8 modal frequencies listed in Table 1 and, clearly, the different modes can be successfully identified at each step. Note that the identified damping ratios, ζ, in Figure 3c are multiplied by the number of modes to avoid overlap surrounding 0.02 (2%). Therefore, the identified results evenly distributed from 0.02 to 0.16 illustrated that the damping ratios over 8 modes are all 2%. The identified mode shapes are also correct but not shown here due to limited space. Overall, RSI is capable to provide correct modal parameters during the seismic excitation.

As a key parameter, i is crucial to the performance of SI and RSI and some literature gives suggestions about how to decide it [28,29,30]. The concept is that i shall be larger than the half-cycle of the fundamental period and the fundamental period is usually the first modal period of the target structure. Hence, it can be computed by multiplying the Nyquist frequency, N_f (in Hz), and the fundamental period, T_f (in seconds). However, in the above example, i is 25 if the measurement is taken via 50 Hz and the window length is 499, which indicates almost 10 s, for the first identification result. Apparently, the suggested criterion is too conservative and inadequate for the practical implementation of building structures or infrastructures under seismic excitations.

Considering the projection is done by the whole part of the matrix, the criterion shall be taken on the data number rather than the row number of the Hankel matrices. To be specific, the proposed criterion is

$$i\ge C_E\frac{2N_f{T}_f+1}{2(m+l+1)}$$

(29)

where C_E is the extension coefficient. C_E shall be at least larger than 1, it can be more than 3 or 4 for stable identified results, and it needs to be further increased while dealing with a noisy measurement. With a 50 Hz measurement, Figure 4 exhibits the identified modal frequencies across various C_E. It can be observed that the identification fails when i is equal to 2 ($C_E=0.78$) and provide stable results when i is around 9 or 10 ($C_E=3.53\mathrm{or}3.92$). Moreover, Figure 4b,c indicate that C_E still works between 1 and 3 after a few recursive steps. Definitely, the proposed criterion brings a reasonable and effective selection for the practical implementation of RSI on seismically excited structures.

3.2. Noise Effect

In the field application, the measurement is always perturbed by noise so it is important to examine the performance using noisy measurement. Signal-to-noise ratio (SNR) is commonly used to represent the level of background noise was adopted in this study. By adding different levels of noise into the measurement, the identified system’s poles are displayed in Figure 5. In the noise-free case shown in Figure 5a, both RSI via the bona fide recursive algorithm and SI via the moving window method can extract the identical 8 system’s poles during the seismic excitation. When noise is blended into the measure, spurious poles appear and dramatically increased (in number) as the noise level. It is observed that the spurious poles only come out under the small row number of the Hankel matrices, implying a larger extension coefficient, C_E, can somehow deal with the noise problem. However, even an enlarged row number (i = 20) can get rid of the spurious poles if SNR keeps reducing, as shown in Figure 5d.

The modal assurance criterion (MAC) is usually used to discriminate between the spurious poles and the true system’s poles while implementing SI [38,39,40,41]. It is defined as the coherence between two vectors, ${\overline{\mathsf{\phi }}}_k$ and ${\tilde{\mathsf{\phi }}}_k$, as

$$\mathrm{MAC}=\frac{{\overline{\mathsf{\phi }}}_k^H{\tilde{\mathsf{\phi }}}_k}{\sqrt{{\overline{\mathsf{\phi }}}_k^H{\overline{\mathsf{\phi }}}_k}\sqrt{{\tilde{\mathsf{\phi }}}_k^H{\tilde{\mathsf{\phi }}}_k}}$$

(30)

where the superscript H is Hermitian transpose. Generally, those vectors are the mode shapes taken out from a stabilization diagram generated over various values of the user-defined parameters, but RSI uses a consistent value and has no stabilization diagram. In the proposed method, one of the vectors can be the column space of the extended observability matrix and the other one can be the reconstructed column space from the eigenvalues and eigenvectors. Sieving these vectors with MAC and an adequate threshold, C_MAC, can remove the spurious modes during the recursive steps. Figure 6 shows the identified system poles with and without MAC criterion under the high noise level. Clearly, all the spurious modes are removed after C_MAC is set to 0.98 and the effectiveness of the proposed method is demonstrated.

3.3. Computation Time

Time consumption for computation is another issue in the field application; actually, it is even more crucial if SHM systems are designed for a real-time or online application. To ensure the computation time, the frame responses are successively down-sampled to several sampling rates, such as 500, 250, 200, 100, 50, and 40 Hz, as the measurement. Two scenarios, fixed i and fixed C_E, are studied and RSI is individually conducted using the same input, output, and parameters in each scenario. Moreover, all the analyses are carried out in the same computer environment for comparison. For the hardware environment, the CPU is Intel^(R) Core^(TM) i7-6700K and the RAM is 32.0 GB. For the software environment, OS is 64 bit Windows and the platform is MATLAB R2021a [42]. The detailed parameters and the results of computation time are listed in

Table 2. The computation time over various sampling rates and a fixed i.

Measurement	2nd, 4th, 6th, and 8th Floor
i	20
CMAC	0.98
Sampling Rate (Hz)	500	250	200	100	50	40
Total Computation Time (seconds)	214.77	100.89	80.58	36.23	15.97	12.03
Averaged Computation Time (ms)	10.96	10.51	10.6	10.06	9.97	10.01

and

Table 3. The computation time over various sampling rates and a fixed α.

Measurement	2nd, 4th, 6th, and 8th Floor
CE	around 4
CMAC	0.98
i	100	50	40	20	10	8
Sampling Rate (Hz)	500	250	200	100	50	40
Total Computation Time (seconds)	5763.64	574.62	268.68	36.23	7.79	4.98
Averaged Computation Time (ms)	320.17	63.83	37.31	10.06	4.32	3.45

.

In Table 2, as expected, the averaged computation time, t_c, per step is quite similar; however, it exponentially increases if C_E is considered in Table 3. Doubtless, too frequent identification may bring no advantage to SHM systems for infrastructures. A ratio is introduced to represent the percentage spent on computation in each step as

$$r_c=\frac{t_c}{\Delta t}\times 100\%$$

(31)

The ratios, r_c, are shown in Figure 7 and must be smaller than 100% to guarantee real-time or online processing. Figure 7a clearly reveals the exponential increase in the second scenario and the averaged computation time surpasses the result from the fixed i around 112 Hz. Furthermore, according to Figure 7b, the sampling ratio should be smaller than 92 Hz so that the proposed recursive formulation can be finished within each sampling interval. Although the computation time varies with different computer environments, one can conclude that down-sampling processing is generally needed for a real-time or online application. Otherwise, the proposed method may not complete the computation in each step because SHM systems commonly have a sampling rate from 100 Hz to 200 Hz [43,44,45]. One can note that on the other hand, the sampling rate needs to be large enough to see those significant frequencies. In some literature, 50 Hz or 60 Hz is a good suggestion for the sampling rate [46,47,48,49].

4. Experimental Validation

The proposed method utilizes the bona fide recursive algorithm with fixed input and output data size to extract the dynamic behaviors. The effectiveness has been demonstrated through a numerically simulated 8-story frame under seismic excitation. Meanwhile, the criterion for selecting the row number of Hankel matrices, the measure for re-moving the spurious modes contaminated by noise, and the down-sampling processing for a real-time or online SHM have been considered and studied with plain suggestions. To further verify the reliability of the proposed method for seismically excited structures, two full-scale specimens are experimentally tested with shaking tables. One is a 3-story steel frame subjected to the El Centro Earthquake (18 May 1940) in the north-south direction. A brace with a special mechanism is designed so that the inter-story stiffness can be reduced during seismic excitations. The other one is a 10-story concrete frame subjected to the Kobe Earthquake (17 January 1995) in three directions. The intensity of excitations is escalated until the peak ground acceleration (PGA) of the original time history. In each experiment, the ground shaking and the structural responses were measured from the accelerometers. Unfortunately, in both experiments RSI was not implemented in real-time because of some technical problems. Instead, the measurement was first recorded and then processed with the bona fide recursive algorithm.

4.1. Experiment of 3-Story Steel Frame

The 3-story specimen is a full-scale steel frame; it was designed and constructed at the National Center for Research on Earthquake Engineering (NCREE) in Taipei, Taiwan. In each story, the specimen is a single bay with 2.15 m by 3.15 m in floor area and 3.0 m in height, as shown in Figure 8. The beams and columns are made by wide flange H-beams (H150 × 150) and each beam-column joint is designed as a bolted connection. Additional mass blocks are located on the floors so that the final weight of each story reaches 6 tons. Then, the 9.0 m height and 18 tons weight specimen are bolted on the shaking table and excited using the pre-specified earthquake. It should be noted that an extra brace is located in the first story. Due to the brace, the approximate stiffness of the first story is 2692 kN/m while the stiffness of the other two stories is 1346 kN/m. As a result, the modes can be identified using the white-noise excitation and the three modal frequencies are 1.23, 3.66, and 5.33 Hz, respectively.

The extra brace is used to produce the stiffness change and, therefore, to mimic the structural damage. In order to achieve that, a lock-up system is designed to connect the brace and the base floor in the first story, as shown in Figure 8. A trigger is manually held by the technician and can be arbitrarily released during the shaking table test. Once the brace is released (un-connected), the inter-story stiffness is immediately reduced to simulate the abrupt change in stiffness. In the shaking table test, the El Centro earthquake is normalized to 100 gal and then inputted to the actuator system to excite the specimen. Here, the unit, gal, is an acceleration unit (1 gal is equal to 1 cm/s²). The lock-up system is assigned to release the brace around 15 s after the earthquake starts. To confirm the true releasing time, an additional displacement meter is installed on the lock-up system in such a way that the relative displacement between the brace and the base floor can be determined.

For the application of SHM, RSI implementation using the bona fide recursive algorithm can detect the abrupt change in stiffness at any instant of time. The measurement was taken using a 200 Hz sampling rate in the test and is down-sampled to 50 Hz for efficient computation. Based on the suggestions given in Section 3, the row number of Hankel matrices, i, is set as 20 and the window length, l_w, is 199 (3.98 s). The threshold for MAC, C_MAC, is set as 0.95 to discriminate between the physical and spurious modes. Hence, the modal parameters identified by RSI and the relative displacement measured by the displacement meter are shown in Figure 9.

The structural damage represented by the stiffness change can be easily observed as the three modal frequencies are lower after 15 s in Figure 9a. The first mode close-up displays two distinct modal frequencies before 15 s and after 19 s, which points out that the dynamic behaviors of the frame have been changed during these 4 s. Considering the fixed window (produce by the bona fide recursive algorithm) and the window length (3.98 s in this test), the frame probably has a very rapid change at around 15 s. The relative displacement measured by the displacement meter can provide a more definite indicator. In Figure 9c, the brace is firmly locked until 14.74 s, generating no displacement at all. After this instant, the trigger is set and the brace is released with the result that the relative vibration is created between the brace and the base floor.

To further confirm the stiffness change produced by the brace and the lock-up system, the hysteresis loop of each story before and after the releasing time is reconstructed using the measurement. Figure 10a–c show the hysteretic behaviors of the first, second, and third stories, respectively. It is noted that the force is reconstructed from the measured acceleration and the floor weights. Unquestionably, the inter-story stiffness in the last two stories is the same before 14.5 s and after 15.0 s, while the one in the first story becomes significantly soft. Moreover, the slope of the hysteresis loop in the first story is similar to the other two after 15.0 s because the inter-story stiffness is supported only by four columns. It can be obviously found that the change of stiffness in the first story can be clearly identified, which is consistent with the original test arrangement. Finally, the computation time for the 46.06 s shaking table test is only 7.94 s, showing the capability for real-time or online processing.

4.2. Experiment of 10-Story Concrete Frame

The second specimen is a 10-story full-scale reinforced concrete (RC) frame, as shown in Figure 11. It was conducted on the world’s largest three-dimensional (3D) earthquake shaking table test facility, E-Defense, at the National Research Institute for Earth Science and Disaster Resilience (NIED) in Hyogo, Japan. The dimensions of the base are 15.7 and 9.7 m in the longitudinal and transverse directions, respectively. On the standard floors, the dimensions are 13.5 and 9.5 m. The longitudinal direction has 3 spans with 4.0 m each and the transverse direction has 3 spans with 3.1, 1.8, and 3.1 m in each floor. The floor heights are 2.8 m for the first floor, 2.6 m for the second to the fourth floors, 2.55 m for the fifth to the seventh floors, and 2.5 m for the eighth to tenth floors. Considering the base, the total height reaches 27.45 m from the shaking table. Furthermore, the longitudinal direction is a simple frame structure composed of beams and columns, whereas the transverse direction is a frame structure with a multi-story shear wall from the first to the seventh floor.

Because of some restrictions, the specimen is divided into upper and lower portions. They were constructed separately and steel beams in the sixth floor are designed to facilitate transportation and assembly. The approximate weight of the steam beams is 16 tons and taken into account in the specimen’s design live load. The design load is around 6 tons for the second and the tenth floors and is around 3 tons for other floors. Moreover, the dead load variates from 70 tons to 97 tons for each floor. Consequently, the modes can be identified using the white-noise excitation; the first three modal frequencies in the longitudinal direction are 1.37, 4.05, and 7.31 Hz, and the three in the transverse direction are 1.49, 5.94, and 12.30 Hz, respectively.

The Kobe earthquake is designated as ground excitation in the shaking table test. The excitation is simultaneously applied in three directions and the intensity of excitations is sequentially increased with 10%, 25%, 50%, and 100% PGA of the original time history. Moreover, white-noise excitations are individually applied in each direction before and after seismic excitations to extract the dynamic features. The sampling rate of measurement was 1000 Hz in all the tests and it is down-sampled to 50 Hz before RSI implementation. According to the suggestions, the row number of Hankel matrices, i, is set as 18 and the window length, l_w, is 431 (8.62 s). The threshold for MAC, C_MAC, is set as 0.95 to discriminate between the physical and spurious modes. Due to limited space, only the identification result in the transverse direction of the 100% PGA’s one is shown here.

From Figure 12a, 0.761 g was achieved by the shaking table, which is 93.27% of the original time history, and the roof response reaches 1.379 g. It is noted that the specimen has been slightly damaged during the test of 50% PGA, so the first three modal frequencies in the transverse direction were reduced to 1.34, 5.25, and 10.96 Hz, respectively. Subsequently, the modal parameters identified by RSI illustrate severe damage, as shown in Figure 12b. The modal frequencies continuously decrease and, finally, converge to 0.87, 3.02, and 6.5 Hz, respectively. According to the literature, the beam-column joint failure further progressed until the concrete cover spalled and the compression side concrete at the base of the multi-story shear wall has crushed in the transverse direction [50].

The hysteresis loop of the specimen’s first story during the shaking table test is reconstructed using the measurement to study the nonlinear behaviors produced by the seismic excitation. Figure 13 shows the comparison of hysteretic behaviors of the specimen over time and the dash line represents the initial stiffness fitted around 10 s after an earthquake starts. The stiffness degradation can be clearly observed after 11 s, indicating that the damage is increasing gradually. The RC frame exhibits significant hysteretic and hardening behaviors to dissipate the huge input energy from the earthquake between 14 and 30 s. In consequence, the identified modal frequencies, which describes equivalent ones, have dramatic reduction during the period with the structural nonlinearity. On the other side, the hysteresis loop is similar between 30 and 40 s, as shown in Figure 13e, so that the identified three mods are quite steady. After 50 s, the excitation is about to cease, the hysteretic behaviors are barely seen, and then the identified modal frequencies equivalently gain a small increase. It can be clearly found that the change of modal frequencies obtained from RSI is consistent with the observation of the hysteresis loop in the first story.

5. Conclusions

Continuous system evaluation and autonomous damage assessment under seismic excitations are critical to structural safety evaluation. Generally, the recursive algorithms are exploited to achieve real-time or online processing for SHM. In this study, RSI implementation using the bona fide recursive algorithm is introduced to identify the time-varying dynamic characteristics under the modal domain. To demonstrate and verify the effectiveness for seismically excited structures, both numerical simulation and experimental investigation have been presented and the results are discussed.

In the numerical simulation, the performance of the proposed method is clarified through an 8-story shear-type frame. The only user-defined parameter, i, is suggested to account for the data number rather than the row number of the Hankel matrices. The extension coefficient, C_E, brings a reasonable and effective selection for the practical implementation of RSI. For the spurious modes, MAC is taken between the column spaces instead of the stability diagram. The sampling rate shall be reduced for a real-time or online application, but it needs to be larger enough to see those significant frequencies. In the experimental investigation, the capability of the proposed method is verified through two full-scale frames. In the test of the 10-story RC frame, the results show that RSI can track the slow time-varying parameter changes and provide the equivalent modal parameters for the nonlinear behaviors. Next, in the test of the 3-story steel frame, the RSI is capable to detect this kind of sudden damage and estimate the event time by leveraging the fixed window shown in the proposed method.

Overall, the bona fide recursive algorithm is demonstrated and verified in this study. Future works will be focused on the implementation of field application and the long-term performance under the effect of the environment, such as temperature, moisture, and so on.