Abstract

Structural health monitoring (SHM) can continuously and nondestructively evaluate the state and performance of structures using the structural responses to external loads or environmental conditions. Moreover, online or real-time SHM of civil structures provides significant advantages over periodic or manual inspection methods, especially under disaster loadings, where the consequences of failure can be severe. To achieve it, performing system identification and damage detection recursively, said recursive subspace identification (RSI), is a promising solution, and SHM based on the algorithms can evaluate damage or deterioration of civil structures, give insight into the health and performance of a structural system, and provide valuable information for decision-making on maintenance and repair. However, the time-consuming decompositions frustrate these algorithms. As a compromise, additional processing is required to implement online and real-time applications. This study demonstrates a modified algorithm that takes advantage of the projection approximation subspace tracking (PAST) algorithm and the repeated system matrices in the extended observability matrix. The modification can reduce numerical decompositions and improve important timeliness for online or real-time SHM of civil structures. Both the numerical simulation and experimental investigation have been used to verify the proposed method, and the results show its capability to determine the changes in the dynamic characteristics of a structure in either the laboratory experiment or in the field application. In the last place, the discussion and some conclusions are also drawn in this paper.

1. Introduction

Structural health monitoring (SHM) is a process to continuously and nondestructively evaluate the state and performance of structures, such as bridges and buildings, and detect any damage and deterioration that may affect their safety, reliability, or lifespan. It utilizes various types of sensors to collect the structural responses to external loads or environmental conditions and then analyzed them using advanced algorithms to identify any changes in the structural behavior that may indicate damage or degradation. For example, many research works tried to quantify the benefit of SHM on structural life-cycle costs through the value of information (VoI) from Bayesian decision theory [15]. Accordingly, SHM is a valuable tool for ensuring the safety, reliability, and longevity of structures and has increasingly been implemented to improve the efficiency and effectiveness of structural monitoring and maintenance, especially for seismic hazard areas.

With advancements in technology over the past decade, including sensing networks, data acquisition, communication, signal processing, information management, the internet of things (IoT), and intelligent algorithms, there has been a growing interest among scientists and engineers in applying real-time SHM to civil structures [6]. It can provide several benefits over periodic or manual inspection methods, and, therefore, it is often necessary for specific applications where safety and reliability are critical, especially under disaster loadings. Other advantages include early detection of damage, continuous monitoring, cost savings, improved safety, etc. Overall, real-time SHM has a significant benefit in constantly monitoring the structural performance and continuously tracking the structural state during natural disasters, especially in critical applications where the consequences of failure can be severe.

One promising way to provide online or real-time SHM is to perform system identification of time-varying dynamic characteristics and damage detection of the structural system [7, 8]. System identification involves analyzing the dynamic behavior of a structure under various loading and environmental conditions and using this information to develop models that can accurately predict its responses to future conditions. In addition, damage detection involves determining the changes in the dynamic characteristics of a structure, such as natural frequencies and damping ratios, and using this information to detect, localize, and quantify the structural damage. Both of them can be performed in different domains; for example, Civera et al. recently proposed frequency-domain techniques [9, 10], and Dessena et al. recently developed time-domain techniques for SHM [11, 12]. Most importantly, online or real-time SHM combining these approaches can evaluate the damage or deterioration of civil structures, give insight into the health and performance of a structural system, and provide valuable information for decision-making on maintenance and repair.

To accurately assess the dynamic behavior of a structure during natural disasters, such as strong earthquakes, considering the time-dependent dynamic characteristics, the focus should be placed on system identification and damage detection methods that can effectively identify time-varying modal parameters. During the past two decades, time-domain methods based on subspace state-space system identification (4SID), or subspace identification (SI) in short, have attracted a great deal of interest in the control and monitoring community because they can identify the system matrices of a state-space model directly from the input and output data. Except for the original 4SID derived by Van Overschee and De Moor [13, 14], a lot of well-developed algorithms were proposed in the 1990s, such as CVA (canonical variate analysis) [1517], PI/PO-MOESP (past input/past output multivariable output error state-space) [1820], IV-4SID (instrumental variable 4SID) [2123], PC (principal component) algorithm [24, 25], UPC (unweighted principal component) algorithm [25], and so on. Unfortunately, those algorithms are not suitable for online or real-time implementation due to the large computational effort of some numerical tools, such as decompositions, and the ability to detect the time-varying modal parameters.

In the 2000s, several recursive algorithms have been developed to provide tracking capability. For example, Oku et al. and Tamaoki et al. try to recursively update the projection matrix used in SI with a similar lemma [26, 27]. Mercère et al. focus on IV-based recursive algorithms and develop various variants [28, 29]. Kameyam et al. rederives the decomposition to a recursive form [30, 31]. Moreover, in the past decade, recursive subspace identification (RSI) algorithms have been successfully applied to detect the time-varying dynamic characteristics under earthquake excitations [3234]. However, the algorithms still need extra processing, like down-sampling, to secure online and real-time implementation [35].

The aim of this paper is to develop an RSI algorithm with a limited computational complexity for online and real-time SHM. In the following sections, SI is first introduced, as are the overall procedures for extracting the dynamic characteristics of a system from the dominant subspace of geometric projection. Subsequently, RSI is reviewed to reveal the computational bottleneck, which is the large number of input/output data and the need for numerical decompositions. Therefore, a modified algorithm is proposed by using the tools from signal processing techniques and the repeated system matrices in the extended observability matrix. The proposed method is elaborately derived and preliminarily exanimated using a numerical simulation. Following the numerical study, two datasets are collected to verify the method; one is from a full-scale specimen which is experimentally tested using a shaking table, and the other one is from a field frame which is installed with an offline SHM system. Through the experimental study, the proposed method shows its capability to determine the changes in the dynamic characteristics of a structure in either the laboratory experiment or in the field application. Last of all, a brief conclusion is drawn, and it may still need a cross-verification from a long-term monitoring system.

2. Introduction to Subspace Identification (SI)

SI is first described as the foundation of the proposed method. Considering that a linear n degree-of-freedom (DOF) structure is subjected to an earthquake excitation, the motion equation can be expressed as a discrete-time state-space equation as follows:where xk is state vector with 2n states; yk is measured output vector with m measurement; uk is input vector with l excitations; wk and vk are the process and measurement noise, respectively; the subscript k denotes k-th step which indicates and Δt is the sampling interval of measurement; A is linear elastic system matrix; B and D are excitation influence vector; C is the output (or observer) matrix. Among those matrices, A and C are generally identified by SI because they are particularly important for the applications of structural control or health monitoring.

2.1. Matrix Input-Output Equations

To derive SI, the state-space equations can be transformed into matrix input-output equations [14] as follows:where Yp and Up are past output and input data Hankel matrices; Yf and Uf are future output and input data Hankel matrices, respectively. All those Hankel matrices have i rows and j columns (both in samples) aswhere ; Up and Uf are defined similarly as Yp and Yf, respectively. As a result, the dimension of Yp and Yf is im by j, and the dimension of Up and Uf is il by j. Similarly, Wp, Wf, Vp, and Vf are the Hankel matrices of the process and measurement noise in the matrix input-output equations. Furthermore, Xp and Xf are past and future state sequences (in column); Δi is the reversed extended controllability matrix; Hi and Gi are the lower block triangular Toeplitz matrices; and Γi is the extended observability matrix.

The extended observability matrix is composed of the information of system matrices (A and C) which are the primary outcome of SI. Details about these matrices can also be found in the literature [14].

2.2. Geometric Projections

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. Hence, two different geometric projections can be implemented: (a) orthogonal projection, Oorth, and (b) oblique projection, Oobl. One can choose either orthogonal or oblique projections for extracting the extended observability matrix. Actually, there are many methods to achieve the projections such as N4SID (numerical algorithms for subspace state-space system identification) [13, 14], CVA (canonical variate analysis) [1517], PI/PO-MOESP (past input/past output multivariable output error state-space) [1820], IV-4SID (instrumental variable 4SID) [2123], PC (principal component) algorithm [24, 25], UPC (unweighted principal component) algorithm [25], and so on. The general form of the projections is written as follows:where and are the left and right nonsingular Hermitian weighting matrices that are possibly data dependent. For example, PO-MOESP computes the projection using QR decomposition and eventually produces the results with and , where denotes the operation that projects the row space of the matrix into the orthogonal complement of the row space of Uf.

The most important observation about the projection is that the extended observability matrix lies in the column space of the projection matrix. One can use singular value decomposition (SVD) to decompose the projection matrix as follows:where U, S, and are left singular vectors, diagonal singular value matrix, and right singular vectors, respectively. and can be separated from the diagonal singular value matrix. The extended observability matrix can be obtained as follows:

The small singular values in S2 induced by noise data are meant to be neglected, and a truncated SVD is described, although DOF might be infinite or remain unknown for a field structure and only a few modes can be identified from the measurements.

2.3. Modal Parameters

Once the extended observability matrix is obtained, modal parameters such as modal frequencies, damping ratios, and mode shapes can be determined by the following procedures. First, the discrete-time version of the linear elastic system matrix, A, can be evaluated as follows:where denotes Γi without the last m rows and denotes Γi without the first m rows. Moreover, output (or observer) matrix, C, can be retrieved from the first m rows of Γi, as shown in equation (7).

Then, the eigenvalues and eigenvectors can be generated using the eigenvalue decomposition (ED) as follows:where Λ is the diagonal eigenvalue matrix which is composed of eigenvalues, and Ψ is the eigenvector matrix. Finally, the natural frequencies and damping ratios can be extracted from the eigenvalues. It is noticed that the eigenvalues and eigenvectors appear in complex conjugated pairs, and a pair of conjugated eigenvalues is associated with the same natural frequency and damping ratio.

Furthermore, the mode shapeswhere Φ is the matrix of mode shapes. The above-mentioned procedures clearly demonstrate how the dynamic characteristics of the structural system can be determined from the extended observability matrix and the extended observability matrix is estimated by dominant subspace computed from projection matrix.

3. Implementation of Recursive Subspace Identification (RSI)

There are several ways to implement SI recursively, and most of them maintain the same procedure as the one in SI. To move the time step forward, these methods try to avoid the construction and projection of Hankel matrices. To be specific, RSI updates the future output and input data vectors instead of constructing Hankel matrices. For example, the future output data vector in k-th time step is . Furthermore, RSI renews the projection with successive computations instead of multiplying the Hankel matrices. Three main streams are available: one is sophisticatedly derived via matrix inversion lemma [26, 27, 36, 37], the other is accomplished by the cross-multiplication of matrices from QR decomposition [30, 31, 38, 39], and the last one rotates the newest vector in QR decomposition to update the projection [28, 29, 40, 41]. Some comparative studies of these two methods can be found in the literature [3235, 42].

Despite the success brought by those methods, the same procedure with SI produces inevitable decomposition steps. Every method needs both SVD and ED for the estimation of the extended observability matrix and the modal parameters, respectively. The last two streams even need QR decomposition to get the initial matrices at the 1-st time step. Admittedly, the larger the input and output data, the more computational complexity is required to perform these decompositions. For online or real-time SHM, the computation of decompositions has been a bottleneck in RSI. To reduce laborious computation, save precious time, and secure timeliness, one can avoid estimating the extended observability matrix using SVD by introducing additional algorithms from different fields. Fortunately, from a signal processing point of view, the projection approximation subspace tracking (PAST) algorithm is another way to implement subspace identification recursively. Furthermore, taking advantage of the repeated system matrices (A and C) in the extended observability matrix can also reduce the computational effort.

3.1. Transformation Matrix for Extended Observability Matrix

First, an alternative to the extended observability matrix is proposed to avoid the decomposition. Considering the projection illustrated by PO-MOESP in equation (8), the orthogonal projection can be expressed as follows:

The orthogonal projection keeps its row number but increases its column number as time moves. For example, the projection in the k-th time step is

The orthogonal projection vector, ok, is reclusively computed in the proposed method instead of the projection matrix to avert increasing matrix size and computational effort. According to equations (14) and (15), the most updated vector can be computed as follows:where , , , and is recursive matrices. These matrices can be sophisticatedly derived using matrix inversion lemma aswhere λ is the forgetting factor used to eliminate the influence of previous steps so as to identify the latest state of the system [32, 42]. From the SI point of view, the input and noise Hankel matrices are canceled out after the orthogonal projection, meaning that the extended observability matrix is still lying in the column space of the orthogonal projection vector as .

By observing the extended observability matrix in equation (7), one can see that system matrices; A and C is repeatedly occurring in each row, so it is unnecessary to collect the full column space of the orthogonal projection vector. Therefore, the extended observability matrix is divided into the upper part and the lower part to reduce the dimension of the vector as follows:and a transformation matrix, T, can be used to transfer between two parts.

Notably, the row number of the upper part (as well as the lower part) cannot be an arbitrary number. Considering the repeatability, the row number of the upper part should be where ; in other words, the row number of the lower part is where α is denoted as a dividing factor. Details about the dividing factor are discussed in the numerical study. According to equations (20) and (21), the following equation can be derivedand serves as the target equation for the PAST algorithm.

3.2. Projection Approximation Subspace Tracking (PAST)

PAST algorithm is one of the successful examples for subspace tracking [43, 44]. It uses a projection-like unconstrained criterion (known as Yang’s criterion) as follows:where z is a signal vector and a matrix argument H; moreover, H is expected to be full rank without loss of generality. The global minimum of equation (23) is attained if and only if H is full of the dominating eigenvectors (with similarity transformation) of . Similar RSI algorithms have been developed based on the criterion, including IV-PAST (instrumental variable PAST), 2IV-PAST (second-order IV-PAST), and EIV-PAST [4547].

Furthermore, taking equation (22) as a target equation and implementing PAST recursively, the criterion is replaced by a series summation with the exponentially forgetting as follows:

By using the matrix inversion lemma, the minimization of equation (24) approximates the criterion, which leads to a recursive least squares (RLS)-like algorithm for tracking the signal subspace [43].

Note that the recursive computation keeps updating the window from 1-st step to k-th step while forgetting memory so that the latest subspace can be tracked.

After the transformation matrix, T, is calculated recursively, the modal parameters can be extracted from equations (11) to (13) by replacing the extended observability matrix with the transformation matrix. As a result, the dynamic characteristics of the structural system can be determined from the dominant subspace. The above-mentioned method has several novelties and can be highlighted herein. First, the projection illustrated by PO-MOESP is approximated by the PAST algorithm to save computation in this study. Then, a dividing factor and a transformation matrix are proposed to perform system identification. The advantages of using dividing factors include increasing the identified accuracy and enhancing the robustness against noise, which will be examined in the following sections. The equivalence of extracting modal parameters from the extended observability matrix and transformation matrix is also demonstrated in this study.

Generally, implementing RSI with the proposed method consists of four steps, as shown in Figure 1. The first step carries out the initial window; all those Hankel matrices are constructed similarly to SI as well as the user-defined parameters [14]. Meanwhile, the other matrices are assumed according to their definitions. The second step constructs the projection matrix, and, afterward, the projection vector, ok, is reclusively computed from equations (16) to (19) once yf,k and uf,k are measured. Then, the projection vector can be divided according to equations (20) and (22). In the third step, the transformation matrix can be computed or updated from equations (25) to (27). Finally, the continuous-time and discrete-time system matrices can be evaluated via equation (11) and, therefore, the modal parameters can be estimated. After the new measurement is acquired for the next step, reiterating from the second step to the fourth step produces a closed-loop recursion. Although ED is still inevitable for the extraction of the linear elastic system matrix, no SVD or QR decomposition is needed in the above-mentioned procedures, and the bottleneck in RSI is therefore eased. Furthermore, the computational complexity is reduced, and the computation time is saved for the application of online or real-time identification.

4. Numerical Study for the Proposed Method

Different from other RSI algorithms, the proposed method avoids using SVD and solves the time-consuming computation involved in online or real-time identification. In the numerical study, a 10-story shear-type frame is excited by the earthquake, El Centro Earthquake (May 18, 1940) in the north-south direction to illustrate the effectiveness of the proposed method. The schematic diagram of the simulated frame is shown in Figure 2(a). The numerical simulation can be done by constructing the motion equation with the following structural parameters. The mass is assumed to be 1 ton and lumped at the center of each story. The stiffness is specified as 20,00 kN/m for each story. The damping ratio is assumed to be 2% for each mode, and the damping matrix is calculated using modal damping. With this configuration, the structural responses are easily evaluated under the 1-dimensional earthquake excitation. Accordingly, the 10 modal frequencies range from 1.06 Hz to 14.08 Hz, as listed in Table 1. The measurement is assigned to be the 10 acceleration responses of all stories, and an additional measurement is taken from the ground acceleration. All the simulation is done with 1 kHz sampling rate via the software, MATLAB, from Math-Works.

4.1. Preliminary Verification Using Numerical Simulation

For the numerical simulation, it is assumed that the response of each story and the ground excitation are measured by the sensors. To study the robustness of the proposed method, a white Gaussian noise with a 120 signal-to-noise ratio (SNR) is artificially added to each sensor to introduce measurement contamination. Hence, the number of measured responses, m, is 10, and the number of measured excitations, l, is 1. Although the numerical simulation is taken using a high sampling rate, the measurement is down-sampled to 100 Hz for efficient computation. Considering the highest modal frequency is less than 15 Hz, a 100 Hz sampling rate (with a 50 Hz Nyquist frequency) is quite sufficient to extract all modes.

Before implementing SI or RSI, the size of Hankel matrices, including i and j, needs to be determined. Once i and j are assigned, the length of initial window, , can be uniquely determined as follows

The length of the initial window is very important because the product of the length of the initial window and the sampling interval is the time for providing the first identification result. It should be small enough, meaning that the size of Hankel matrices should be small enough, to provide a timely warning in the field application [35]. In the meantime, i needs to be larger enough to consider the fundamental period, which is usually the first modal period of the target structure. The minimum requirement can be estimated by the fundamental period, the number of measurements (both input and output), and Nyquist frequency [35, 42]. In the numerical verification, 100 is adopted for both i and j.

To eliminate the spurious modes identified by subspace methods, the modal assurance criterion (MAC) is needed while implementing SI or RSI [4851]. It is defined as the coherence between two vectors, and , as follows:where the superscriptH is Hermitian transpose. Generally, those vectors are the mode shapes taken out from a stabilization diagram generated over various i, but RSI has no stabilization diagram and uses a consistent value during online or real-time identification. 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 [35]. For example, at each step, each column of the extended observability matrix obtained from equation (21) can be one vector for calculating MAC. The eigenvalues and eigenvectors can also be used to reconstruct the system matrices (A and C) as well as the extended observability matrix according to their definitions shown in equation (7); each column can be the other vector, too. Most spurious modes cannot keep their shapes after reconstruction and, therefore, have a very small coherence. By choosing an adequate threshold, CMAC, MAC can sieve these vectors and remove the spurious modes during the recursive steps. In the numerical verification, 0.98 is adopted for CMAC.

For the other user-defined parameters, the forgetting factor, λ, and the dividing factor, α, are set as 0.95 and 2, respectively, and the identification result is shown in Figure 2. Obviously, RSI with the proposed method gives an excellent estimation of the modal parameters. The grids in Figure 2 are drawn according to the correct modal frequencies and damping ratio. It is noteworthy that the identified damping ratio is multiplied by the mode number for concise separation, avoiding overlap surrounding 0.02. Therefore, the values evenly distributed from 0.02 to 0.2 illustrate that the damping ratios over 10 modes are all 0.02 (2%). This result well demonstrates that the different modes can be successfully identified at each step once the measurement is noise-free, which is represented by the black circles. Certainly, the identified result is deteriorated by contaminated measurements; the higher modes cannot be extracted under these noise conditions, as shown by red crosses. Furthermore, those modes could be again identified by decreasing the MAC threshold, CMAC. For example, adopting a value of 0.95 for CMAC retrieves the eighth mode. However, this choice creates a trade-off between noise and spurious modes, suggesting that the user’s experiences may influence the final performance. Besides, the identified mode shapes are also correct but not shown here due to limited space. Overall, RSI with the proposed method is capable of providing correct modal parameters during the seismic excitation.

4.2. Study on Dividing Factor

In the proposed method, the extended observability matrix is divided into upper and lower parts by the dividing factor, α. The row numbers of the upper and lower parts is and , respectively, and the dividing factor can be assigned by any positive integer smaller than i. Therefore, 6 different values are selected to study the performance numerically. By setting the dividing factors as 1, 2, 3, 4, 6, and 8, the 10 modal frequencies are identified, and the errors are then calculated, as shown in Figure 3. The horizontal axis in the bar chart represents the errors in percentage and is limited only to ±3%. The vertical axis is the cumulative numbers of the identified frequencies, and the number are listed on the top of the bar.

Clearly, the identified modes demonstrate high accuracy in Figure 3; most of the results are located within ±1% error. To be specific, the overall accuracy within ±1% error is 91.06%, 91.25%, 91.38%, 92.16%, and 92.36% when α set as 2, 3, 4, 6, and 8, respectively. The percentage slightly increases as α increases although the variation is minute. It should be noted that PAST is a RLS-like algorithm so it is hard to obtain reliable solutions in the first few recursive steps. Consequently, the first results only appeared after 5 seconds in Figure 2. If the gap before the first identification result is ignored, the overall accuracy within ±1% error is boosted to 97.97%, 98.17%, 98.31%, 99.15%, and 99.37%. Moreover, the accuracy of lower modes is higher than that of high modes. For example, around 99.45% to 100% of the first modal frequency can be correctly identified (with ±1% error), but only 96.8% to 97.73% of the tenth modal frequency can be found with the same deviation. The difference comes from the vibration energy; lower vibration modes carry more energy so that it is easy to be identified.

To further study the dividing factors, contaminated measurements are again considered, and the errors of the modal frequencies are shown in Figure 4. Most of the results are still located within ±1% error, indicating that the identified values are unrelated to noise. The overall accuracy within ±1% error is 37.95%, 46.36%, 50.42%, 53.91%, and 55.13% when α set at 2, 3, 4, 6, and 8, respectively. The percentage significantly increases as α increases since the transformation matrix has enough column and row spaces to against the noise effect. Fortunately, the first three fundamental modes are more robust to noise and can be correctly identified with various α. The identification from the fourth to eighth mode may sometimes be lost under noise conditions, and the modes are relatively hard to extract in the higher modes. The ninth and tenth modes are barely reported because of weak energy compared to noise.

Despite the accuracy of the various dividing factors, Figure 3 (as well as Figure 4) also reveals that produces the worst performance and the difference is significant. The first mode can only be identified over 823 steps (23.92% accuracy) with this factor, and no frequencies can be extracted from the third mode to the tenth mode. Considering the dimensions of the transformation matrix (), and the linear elastic system matrix (), the dynamic characteristics of the structural system can be fully extracted if and only if and . In other words, the transformation matrix needs to provide enough column and row spaces (as well as ranks) to solve every mode in a dynamic system otherwise, some of them would be lost. In the numerical verification studied with , is equal to 10 which is smaller than 2n (20), leading the poor identification. To further exanimate the inference, a 1 DOF simulation is conducted, both the acceleration and displacement responses are collected, and the proposed method is again used to identify the modal parameters. In this simulation, the modal frequency and damping ratio can be successfully identified at each step. Admittedly, even though the dividing factor is set to 1, the mode can be perfectly identified once the transformation matrix has enough ranks.

4.3. Study on Time Efficiency

The consumption of computation time is an important issue in field applications, especially if SHM systems are designed for an online or real-time application. To study the computation time, RSI is conducted using the proposed method with different dividing factors, α, and sampling rates. To examine the efficiency of the proposed method, different methods are also conducted with similar user-defined parameters for comparison. Moreover, all the study is done in the same computer environment and platform. For the hardware environment, the CPU is Intel(R) Core(TM) i7-4790K and the RAM is 32.0 GB. For the software environment, the OS is 64bit Windows, and the platform is MATLAB R2022b (The Math Works, Inc., 2022). The detailed parameters and the results of the computation time are listed in Tables 2 and 3.

To study efficiency, the effect of different dividing factors on the time consumption is first investigated, and it has been set as 1, 2, 3, 4, 6, and 8 while the other configurations and the other user-defined parameters are identical, as shown in Table 2. Apparently, because of the matrix sizes in the recursive steps, the difference is observable, and the computation time is increased as the dividing factors. In this study, it shall be smaller than 4 since the 100 Hz sampling rate bring us only 10 milliseconds sampling interval. However, it is still free to use a larger factor if the sampling rate is lower in another case. Noteworthily, the computation time of the case is conducted using only acceleration responses, meaning that the results are in accurate (as shown in the above subsection). Once the 10 displacement responses are included for better identification, the computation time for each step is escalated to 14.53 milliseconds and eventually goes beyond the sampling interval. To end up, α is set as 2 for all analyses in the numerical study.

For the study of sampling rate, the acceleration responses are successively down-sampled to several sampling rates, such as 250, 200, 100, 50, 40, and 20 Hz, as the measurement. Corresponding i and j are selected, and RSI is individually conducted with other user-defined parameters the same. The details and the computation time are listed in Table 3. Doubtless, a high sampling rate increases computation time and may be disadvantageous to SHM. The fitting result is shown in Figure 5 where an overhead time can be estimated at 3.19 milliseconds. The curve clearly reveals the exponential increase in the computation time and eventually surpasses the sampling interval before the 150 Hz sampling rate. Therefore, the sampling rate should be equal to or smaller than the 100 Hz numerical study so that the proposed method can be finished within each sample. Although the computation time varies with different computer environments, the down-sampling processing could be a common solution for an online or real-time application. It can be noted that, on the other hand, the sampling rate needs to be large enough to see those significant frequencies.

To compare the efficiency of the proposed method, three convention methods are implemented using similar user-defined parameters. The first one is moving window SI which successively performs SI in each time window. For comparison, the orthogonal projection is selected to extract the extended observability matrix. The second one is the RSI method proposed by Tamamoki et al., and the recursive formulation is also developed based on the orthogonal projection [27]. The last RSI method was proposed by Oku et al. and is derived from oblique projection [26, 36, 37]. All the cases are applied with the same configurations, and the results are displayed in Table 4. Obviously, all these three methods need extensive time to complete the computation in each step. The second method proposed by Tamamoki et al. uses less time compared to the moving window SI, but the last method needs even more time to complete the computation because its formulation for oblique projection is built on the orthogonal projection. The computation time of the three methods confirms the efficiency of the proposed method (which needs only 5.1 milliseconds in each step); the orthogonal projection vector, the transformation matrix, and the PAST algorithm dramatically reduce laborious computation, save precious time, and secure timeliness.

5. Experimental Study and Verification

The effectiveness of the proposed method has been investigated through a numerically simulated 10-story frame under earthquake excitation. Meanwhile, the correctness of various dividing factors and the potential for efficient identification have been studied to support practical usage. To further demonstrate the proposed method, two datasets are collected and analyzed: one is from a full-scale specimen which is experimentally tested using a shaking table, and the other is from a field frame which is installed with an offline SHM system. Both datasets contain acceleration responses with a 200 Hz sampling rate, and they are down-sampled to 100 Hz to consider the processing time for online or real-time applications.

5.1. Experiment with a 3-Story Steel Frame

To experimentally study the proposed method, a full-scale steel specimen 3 stories was designed and constructed at the National Center for Research on Earthquake Engineering (NCREE) in Taipei, Taiwan. The 3-story steel frame is a single-bay structure with 2.15 meters wide, 3.15 meters long, and 3.0 meters height. In each story, both the beams and columns are made by wide flange H-beams (H150 × 150), and they are all connected using bolts. Mass blocks are added to each story to adjust the dynamic characteristics so that the final weight is increased to 6 tons per story. Then, the full-scale 3-story steel frame was bolted to the shaking table and tested using the white-noise and earthquake excitation. Hence, the modes identified using the white noise excitation are 1.23, 3.66, and 5.33 Hz initially [35]. Moreover, the El Centro earthquake is normalized to 100 gal and then inputted 468 into the actuator system for earthquake excitation.

It is important to note that the first story has an extra brace, as shown in Figure 6(a). The approximate stiffness of the first story is doubled due to the brace, and the original stiffness of all the stories is 1346 kN/m. This extra brace is used to generate a stiffness change as well as structural damage. A lock-up system is attached to the brace and the base floor to connect or release the brace, as shown in Figure 6(b). A trigger signal can be sent to the system by the technician at any instant during the shaking table test. Once the signal is sent, the lock-up system releases the brace, and the abrupt change in stiffness is immediately simulated by reducing the interstory stiffness. In the El Centro earthquake, the technician released the lock-up system (and the brace) around 15 seconds after the earthquake starts. In addition, the lock-up system is installed with an additional displacement meter, and the relative displacement between the brace and the base floor can be determined to confirm the true releasing time.

To verify that the proposed method is capable of tracking the modal parameters, the specimen is excited by the shaking table, and the acceleration responses are measured by accelerometers installed on each floor. An additional accelerometer is mounted on the shaking table to measure the ground acceleration. The above measurement is used to implement RSI. For the user-defined parameters, the size of Hankel matrices, i and j, is set as 100, the threshold of MAC, CMAC, is set as 0.9, the dividing factor, α, is set as 8, and the forgetting factor, λ, is set as 0.95. The total computation time is 20.46 seconds, which is only half of the 46 seconds long measurement. If the identification is done by moving window SI, recursive orthogonal projection, or oblique projection, the total computation time is increased to 56.23, 73.14, and 77.75 seconds, respectively. As a result, the modal frequencies identified by RSI are shown in Figure 7. The structural damage represented by the stiffness change can be easily observed as the three modal frequencies decline after 15 seconds in Figure 7(b). Notably, the proposed method captures an additional frequency (16.25 Hz) that has not been observed before 15 seconds. After conducting a thorough investigation, the mode belongs to the brace fixed below the floor of the second story.

To investigate the effect of the forgetting factor, the modal frequencies identified using different forgetting factors and the relative displacement measured by the displacement meter are shown in Figure 8. Definitely, the relative displacement provides an accurate releasing time. In Figure 8(c), the brace is firmly locked until 14.74 seconds, generating no displacement at all. After this moment, the trigger is set, the brace is released, and the relative vibration is generated between the brace and the base floor at a 16.25 Hz natural frequency. All modal frequencies are clearly decreased from 15 to 19 seconds, which points out that the dynamic behaviors of the frame have been changed during these 4 seconds. Four forgetting factors are used to compare the results. Considering the changing trend, a smaller forgetting factor can detect the damage more rapidly. For example, the first modal frequency identified by keeps time-delayed tracking before 19 seconds; however, the ones identified by and is already stabilized after 17 seconds in Figure 8(b). Admittedly, the proposed method with a small forgetting factor can immediately detect the damage although the results identified by a larger forgetting factor are more stable.

5.2. Field Measurement under Seismic Event

To further verify the proposed method, a dataset collected from an offline SHM system during Chi-chi Earthquake (September 21, 1999) is utilized. The SHM system is installed on the 7-story reinforced concrete (RC) frame of the National Chung Hsing University in Taichung, Taiwan. The frame was constructed in 1992 and was strongly impacted by the earthquake in 1999. From the Chi-Chi earthquake reconnaissance report, it is identified as moderate damage after the seismic event. Fortunately, the building is equipped with 29 accelerometers, and the sensors recorded all the acceleration responses during this earthquake excitation. Figure 9(b) demonstrates the side view and (first floor) plan view of this building and the locations of all the accelerometers with their channel numbers; meanwhile, Figures 9(a) and 9(c) show the ground acceleration measured from the basement in the longitudinal direction, said channels 4 and 8, and the acceleration responses measured from the roof, said channels 18 and 21. It can be observed that the peak ground acceleration (PGA) is over 250 gal and the peak acceleration response is almost 650 gal.

Considering that the field seismograph is unaligned with the building, the ground acceleration measured from the basement is taken as input data, while acceleration responses on the fourth floor and roof are used as output data for implementing the proposed method. In this example, the size of Hankel matrices, i and j, is set as 100, the threshold of MAC, CMAC, is set as 0.9, the dividing factor, α, is set as 2, and the forgetting factor, λ, is set as 0.98. The total computation time is 29.41 seconds, around one-third of the 90 seconds long measurement. If the identification is done by moving window SI, recursive orthogonal projection, or oblique projection, the total computation time is increased to 94.94, 122.31, and 122.01 seconds, respectively; those results again confirm the timeliness for online and real-time applications. Only longitudinal modes are identified by RSI, and the modal frequencies are displayed in Figure 10. Two modes can be identified, and they have apparently decreased from 40 seconds to 70 seconds, which corresponds to the major part of the seismic event. The first modal frequency is reduced from 2.98 Hz to 2.09 Hz and, despite the scattering, the second mode changes from 8.52 Hz to 6.39 Hz. It is believed that moderate damage lowers the modal frequency by at least 25%. Moreover, the time-varying dynamic characteristics just happened in between the major seismic waves, and the identified results are relatively stable either before 40 seconds or after 70 seconds. Consequently, the proposed method is able to identify the changes in structural behavior and track the time-varying modal parameters.

5.3. Potential Challenges and Limitations of the Proposed Method

The proposed method has been used to demonstrate the system identification of time-varying dynamic characteristics and the damage detection of the structural system. The results have shown that the developed RSI algorithm is suitable for online and real-time SHM. However, there are some potential challenges and limitations that should be highlighted and addressed before practical implementation.(i)Although the performance of the proposed method under noise conditions has been roughly investigated through numerical simulation, the identification results could be worse if significant noise is present. As a potential extension, the signal vector in equation (20) could be replaced by the cross-covariance, known as IV-PAST, 2IV-PAST, or EIV-PAST, to enhance the robustness of the proposed method.(ii)Another limitation of the proposed method is the assumption of linear behaviors shown in the derivation. When civil structures exhibit non-linear behaviors, the proposed method generates an equivalent linear model that averages the modal parameters over a short period of time. Those identification results may sometimes be confused with the structural damage. The hysteretic behaviors could be helpful in differentiating them if hysteresis loops are available.(iii)Obviously, the proposed method extracts the modal parameters from the input-output relationship under seismic events. It indicates that ground excitations must be measured during earthquakes. Considering that ground motions do not vary significantly in a small area, the use of free-field seismometers near structures or accelerometers installed in lower-level basements could be considered as an alternative if SHM systems are unable to measure earthquake excitations.(iv)So far, this study has demonstrated the system identification and damage detection of building structures. The proposed method, however, can be applied to various civil infrastructures such as bridges, tunnels, dams, and more. This is because the derivation does not require any specific geometry or input form. The investigation can be continued for further verification.

6. Conclusion

Online or real-time SHM of civil structures provides significant advantages under disaster loadings. As a promising solution, RSI can perform system identification and damage detection recursively to detect the deterioration of a structural system. In this study, a modified algorithm is proposed that takes advantage of the PAST algorithm and the repeated system matrices in the extended observability matrix. The recursive formulation is first derived with great detail, and the user-defined parameters are studied to provide selection guidance while online or real-time implementation. Additionally, both the numerical simulation and experimental investigation have been presented, and the following conclusion can be obtained accordingly:(i)Either with numerically simulated or experimentally measured responses, RSI with the proposed method is capable of extracting modal parameters during the seismic excitation, even if the parameters are time-varying ones.(ii)The identified accuracy from the proposed method is high, and it slightly increases as the dividing factor increases; however, more computation time is required with a larger factor.(iii)The proposed method with a small forgetting factor can immediately detect the damage although the results identified by a larger forgetting factor are more stable.(iv)Although down-sampling processing is a common solution for online or real-time applications, the proposed method shows its potential for providing identification with the original sampling rate and is suitable for the practical implementation of the SHM system.(v)The computation time of the convention methods confirms the efficiency of the proposed method, and the difference demonstrates that timeliness can be secured even with 100 sampling rate measurement, which is common in the SHM system.

The RSI with the proposed method is able to continuously identify changes in the structural behavior, track the structural state, and monitor the structural performance. The evaluation of damage or deterioration of civil structures gives valuable information for decision-making on maintenance and repair. However, some further verification is still recommended, especially a real-time example in field applications to demonstrate the timely performance of SHM in civil structures.

Data Availability

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

The study was supported by the National Science and Technology Council under contract No. NSTC 112-2221-E-005-033- and the National Center for Research on Earthquake Engineering (NCREE, Taiwan).