Development of Recursive Subspace Identification for Real-Time Structural Health Monitoring under Seismic Loading

,


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 afect 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 beneft of SHM on structural life-cycle costs through the value of information (VoI) from Bayesian decision theory [1][2][3][4][5].Accordingly, SHM is a valuable tool for ensuring the safety, reliability, and longevity of structures and has increasingly been implemented to improve the efciency and efectiveness 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 benefts over periodic or manual inspection methods, and, therefore, it is often necessary for specifc 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 signifcant beneft 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 identifcation of time-varying dynamic characteristics and damage detection of the structural system [7,8].System identifcation 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 diferent 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 identifcation and damage detection methods that can efectively identify time-varying modal parameters.During the past two decades, timedomain methods based on subspace state-space system identifcation (4SID), or subspace identifcation (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) [15][16][17], PI/PO-MOESP (past input/past output multivariable output error state-space) [18][19][20], IV-4SID (instrumental variable 4SID) [21][22][23], 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 efort 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 identifcation (RSI) algorithms have been successfully applied to detect the timevarying dynamic characteristics under earthquake excitations [32][33][34].However, the algorithms still need extra processing, like down-sampling, to secure online and realtime implementation [35].
Te 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 frst 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.Terefore, a modifed algorithm is proposed by using the tools from signal processing techniques and the repeated system matrices in the extended observability matrix.Te 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 feld frame which is installed with an ofine SHM system.Trough 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 feld application.Last of all, a brief conclusion is drawn, and it may still need a cross-verifcation from a long-term monitoring system.

Introduction to Subspace Identification (SI)
SI is frst 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 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∆t and Δt is the sampling interval of measurement; A is linear elastic system matrix; B and D are excitation infuence vector; C is the output (or observer) matrix.Among those matrices, A and C are generally identifed by SI because they are particularly important for the applications of structural control or health monitoring.

Matrix Input-Output Equations.
To derive SI, the statespace equations can be transformed into matrix inputoutput equations [14] as follows: where y j|i+j � y T j y T j+1 Te 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].

Geometric Projections.
A geometric tool called projection in the feld of linear algebra is utilized to obtain the extended observability matrix from the matrix input-output equations.Hence, two diferent geometric projections can be implemented: (a) orthogonal projection, O orth , and (b) oblique projection, O obl .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 identifcation) [13,14], CVA (canonical variate analysis) [15][16][17], PI/PO-MOESP (past input/past output multivariable output error state-space) [18][19][20], IV-4SID (instrumental variable 4SID) [21][22][23], PC (principal component) algorithm [24,25], UPC (unweighted principal component) algorithm [25], and so on.Te general form of the projections is written as follows: where W L and W R 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 , where Π ⊥ U f denotes the operation that projects the row space of the matrix into the orthogonal complement of the row space of U f .
Te 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 V are left singular vectors, diagonal singular value matrix, and right singular vectors, respectively.S 1 ∈ R 2n×2n and S 2 ≈ 0 can be separated from the diagonal singular value matrix.Te extended observability matrix can be obtained as follows: Te small singular values in S 2 induced by noise data are meant to be neglected, and a truncated SVD is described, although DOF might be infnite or remain unknown for a feld structure and only a few modes can be identifed from the measurements.

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 discretetime version of the linear elastic system matrix, A, can be evaluated as follows: where ×2n denotes Γ i without the last m rows and ×2n denotes Γ i without the frst m rows.Moreover, output (or observer) matrix, C, can be retrieved from the frst m rows of Γ i , as shown in equation (7).
Ten, 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 shapes where Φ is the matrix of mode shapes.Te 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.

Implementation of Recursive Subspace Identification (RSI)
Tere 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 Structural Control and Health Monitoring construction and projection of Hankel matrices.To be specifc, 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.Tree 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 [32][33][34][35]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.Te 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 diferent felds.Fortunately, from a signal processing point of view, the projection approximation subspace tracking (PAST) algorithm is another way to implement subspace identifcation recursively.Furthermore, taking advantage of the repeated system matrices (A and C) in the extended observability matrix can also reduce the computational efort.

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: Te 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 Te orthogonal projection vector, o k , is reclusively computed in the proposed method instead of the projection matrix to avert increasing matrix size and computational efort.According to equations ( 14) and ( 15), the most updated vector can be computed as follows: where is recursive matrices.Tese matrices can be sophisticatedly derived using matrix inversion lemma as where λ is the forgetting factor used to eliminate the infuence 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.Terefore, 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 αm where α ∈ N; in other words, the row number of the lower part is (i − α)m 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 derived and serves as the target equation for the PAST algorithm.

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.Te global minimum of equation ( 23 Furthermore, taking equation ( 22) as a target equation and implementing PASTrecursively, 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.Te abovementioned 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.Ten, a dividing factor and a transformation matrix are proposed to perform system identifcation.Te advantages of using dividing factors include increasing the identifed accuracy and enhancing the robustness against noise, which will be examined in the following sections.Te 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.Te frst step carries out the initial window; all those Hankel matrices are constructed similarly to SI as well as the user-defned parameters [14].Meanwhile, the other matrices are assumed according to their defnitions.Te second step constructs the projection matrix, and, afterward, the projection vector, o k , is reclusively computed from equations ( 16) to (19) once y f,k and u f,k are measured.Ten, 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 abovementioned 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 identifcation.

Numerical Study for the Proposed Method
Diferent from other RSI algorithms, the proposed method avoids using SVD and solves the time-consuming computation involved in online or real-time identifcation.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 efectiveness of the proposed method.Te schematic diagram of the simulated frame is shown in Figure 2(a).Te numerical simulation can be done by constructing the motion equation with the following structural parameters.Te mass is assumed to be 1 ton and lumped at the center of each story.Te stifness is specifed as 20,00 kN/m for each story.Te damping ratio is assumed to be 2% for each mode, and the damping matrix is calculated using modal damping.With this confguration, the structural responses are easily evaluated under the 1dimensional earthquake excitation.Accordingly, the 10 modal frequencies range from 1.06 Hz to 14.08 Hz, as listed in Table 1.Te 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.

Preliminary Verifcation 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 artifcially 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 downsampled to 100 Hz for efcient computation.Considering the highest modal frequency is less than 15 Hz, a 100 Hz sampling rate (with a 50 Hz Nyquist frequency) is quite sufcient 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, l w , can be uniquely determined as follows   Te 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 frst identifcation result.It should be small enough, meaning that the size of Hankel matrices should be small enough, to provide a timely warning in the feld application [35].In the meantime, i needs to be larger enough to consider the fundamental period, which is usually the frst modal period of the target structure.Te minimum requirement can be estimated by the fundamental period, the number of  Structural Control and Health Monitoring measurements (both input and output), and Nyquist frequency [35,42].In the numerical verifcation, 100 is adopted for both i and j.
To eliminate the spurious modes identifed by subspace methods, the modal assurance criterion (MAC) is needed while implementing SI or RSI [48][49][50][51].It is defned as the coherence between two vectors, φ and  φ, as follows: where the superscript H 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 identifcation.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.Te 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 defnitions 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, C MAC , MAC can sieve these vectors and remove the spurious modes during the recursive steps.In the numerical verifcation, 0.98 is adopted for C MAC .
For the other user-defned parameters, the forgetting factor, λ, and the dividing factor, α, are set as 0.95 and 2, respectively, and the identifcation result is shown in Figure 2. Obviously, RSI with the proposed method gives an excellent estimation of the modal parameters.Te grids in Figure 2 are drawn according to the correct modal frequencies and damping ratio.It is noteworthy that the identifed damping ratio is multiplied by the mode number for concise separation, avoiding overlap surrounding 0.02.Terefore, the values evenly distributed from 0.02 to 0.2 illustrate that the damping ratios over 10 modes are all 0.02 (2%).Tis result well demonstrates that the diferent modes can be successfully identifed at each step once the measurement is noise-free, which is represented by the black circles.Certainly, the identifed 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 identifed by decreasing the MAC threshold, C MAC .For example, adopting a value of 0.95 for C MAC retrieves the eighth mode.However, this choice creates a trade-of between noise and spurious modes, suggesting that the user's experiences may infuence the fnal performance.Besides, the identifed 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.

Study on Dividing
Factor.In the proposed method, the extended observability matrix is divided into upper and lower parts by the dividing factor, α.Te row numbers of the upper and lower parts is αm and (i − α)m, respectively, and the dividing factor can be assigned by any positive integer smaller than i.Terefore, 6 diferent 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 identifed, and the errors are then calculated, as shown in Figure 3. Te horizontal axis in the bar chart represents the errors in percentage and is limited only to ±3%.Te vertical axis is the cumulative numbers of the identifed frequencies, and the number are listed on the top of the bar.
Clearly, the identifed modes demonstrate high accuracy in Figure 3; most of the results are located within ±1% error.To be specifc, 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.Te 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 frst few recursive steps.Consequently, the frst results only appeared after 5 seconds in Figure 2. If the gap before the frst identifcation 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 frst modal frequency can be correctly identifed (with ±1% error), but only 96.8% to 97.73% of the tenth modal frequency can be found with the same deviation.Te diference comes from the vibration energy; lower vibration modes carry more energy so that it is easy to be identifed.
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 identifed values are unrelated to noise.Te 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.Te percentage signifcantly increases as α increases since the transformation matrix has enough column and row spaces to against the noise efect.Fortunately, the frst three fundamental modes are more robust to noise and can be correctly identifed with various α.Te identifcation 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.Te 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 α � 1 produces the worst performance and the diference is signifcant.Te frst mode can only be identifed 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 (T ∈ R (i− α)m×mα ), and the linear elastic system matrix (A ∈ R 2n×2n ), the dynamic characteristics of the structural system can be fully extracted if and only if 2n ≥ αm and 2n ≥ (i − α)m.In other words, the transformation matrix   10 Structural Control and Health Monitoring 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 verifcation studied with α � 1, αm is equal to 10 which is smaller than 2n (20), leading the poor identifcation.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 identifed at each step.Admittedly, even though the dividing factor is set to 1, the mode can be perfectly identifed once the transformation matrix has enough ranks.

Study on Time Efciency.
Te consumption of computation time is an important issue in feld applications, especially if SHM systems are designed for an online or realtime application.To study the computation time, RSI is conducted using the proposed method with diferent dividing factors, α, and sampling rates.To examine the efciency of the proposed method, diferent methods are also conducted with similar user-defned 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 (Te Math Works, Inc., 2022).Te detailed parameters and the results of the computation time are listed in Tables 2 and 3.
To study efciency, the efect of diferent dividing factors on the time consumption is frst investigated, and it has been set as 1, 2, 3, 4, 6, and 8 while the other confgurations and the other user-defned parameters are identical, as shown in Table 2. Apparently, because of the matrix sizes in the recursive steps, the diference 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 α � 1 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 identifcation, 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-defned parameters the same.Te details and the computation time are listed in Table 3. Doubtless, a high sampling rate increases computation time and may be disadvantageous to SHM.Te ftting result is shown in Figure 5 where an overhead time can be estimated at 3.19 milliseconds.Te curve clearly reveals the exponential increase in the computation time and eventually surpasses the sampling interval before the 150 Hz sampling    rate.Terefore, the sampling rate should be equal to or smaller than the 100 Hz numerical study so that the proposed method can be fnished within each sample.Although the computation time varies with diferent 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 signifcant frequencies.
To compare the efciency of the proposed method, three convention methods are implemented using similar userdefned parameters.Te frst 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.Te second one is the RSI method proposed by Tamamoki et al., and the recursive formulation is also developed based on the orthogonal projection [27].Te 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 confgurations, and the results are displayed in Table 4. Obviously, all these three methods need extensive time to complete the computation in each step.Te 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.Te computation time of the three methods confrms the efciency 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.

Experimental Study and Verification
Te efectiveness 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 efcient identifcation 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 feld frame which is installed with an ofine 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 realtime applications.

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.Te 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 fange 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 fnal weight is increased to 6 tons per story.Ten, the full-scale 3-story steel frame was bolted to the shaking table and tested using the whitenoise and earthquake excitation.Hence, the modes identifed 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 frst story has an extra brace, as shown in Figure 6(a).Te approximate stifness of the frst story is doubled due to the brace, and the original stifness of all the stories is 1346 kN/m.Tis extra brace is used to generate a stifness change as well as structural damage.A lock-up system is attached to the brace and the base foor 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 stifness is immediately simulated by reducing the interstory stifness.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 foor can be determined to confrm 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 foor.An additional accelerometer is mounted on the shaking table to measure the ground acceleration.Te above measurement is used to implement RSI.For the user-defned parameters, the size of Hankel matrices, i and j, is set as 100, the threshold of MAC, C MAC , is set as 0.9, the dividing factor, α, is set as 8, and the forgetting factor, λ, is set as 0.95.Te total computation time is 20.46 seconds, which is only half of the 46 seconds long measurement.If the identifcation 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 identifed by RSI are shown in Figure 7. Te structural damage represented by the stifness 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 fxed below the foor of the second story.
To investigate the efect of the forgetting factor, the modal frequencies identifed using diferent forgetting factors and the relative displacement measured by the displacement meter are shown in Figure 8. Defnitely, the relative displacement provides an accurate releasing time.In Figure 8(c), the brace is frmly 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 foor 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 frst modal frequency identifed by λ � 0.99 keeps time-delayed tracking before 19 seconds; however, the ones identifed by λ � 0.95 and λ � 0.93 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 identifed by a larger forgetting factor are more stable.

Field Measurement under Seismic Event.
To further verify the proposed method, a dataset collected from an ofine SHM system during Chi-chi Earthquake (September 21, 1999) is utilized.Te SHM system is installed on the 7story reinforced concrete (RC) frame of the National Chung Hsing University in Taichung, Taiwan.Te frame was constructed in 1992 and was strongly impacted by the  Structural Control and Health Monitoring earthquake in 1999.From the Chi-Chi earthquake reconnaissance report, it is identifed 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.Considering that the feld seismograph is unaligned with the building, the ground acceleration measured from the basement is taken as input data, while acceleration responses on the fourth foor 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, C MAC , is set as 0.9, the dividing factor, α, is set as 2, and forgetting factor, λ, is set as 0.98.Te total computation time is 29.41 seconds, around one-third of the 90 seconds long measurement.If the identifcation 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 confrm the timeliness for online and real-time  Structural Control and Health Monitoring applications.Only longitudinal modes are identifed by RSI, and the modal frequencies are displayed in Figure 10.Two modes can be identifed, and they have apparently decreased from 40 seconds to 70 seconds, which corresponds to the major part of the seismic event.Te frst 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 identifed 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.Structural Control and Health Monitoring

Potential Challenges and Limitations of the Proposed
Method.Te proposed method has been used to demonstrate the system identifcation of time-varying dynamic characteristics and the damage detection of the structural system.Te 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 identifcation results could be worse if signifcant 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.Tose identifcation results may sometimes be confused with the structural damage.Te hysteretic behaviors could be helpful in diferentiating 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 signifcantly in a small area, the use of free-feld 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 identifcation and damage detection of building structures.Te proposed method, however, can be applied to various civil infrastructures such as bridges, tunnels, dams, and more.Tis is because the derivation does not require any specifc geometry or input form.Te investigation can be continued for further verifcation.

Conclusion
Online or real-time SHM of civil structures provides signifcant advantages under disaster loadings.As a promising solution, RSI can perform system identifcation and damage detection recursively to detect the deterioration of a structural system.In this study, a modifed algorithm is proposed that takes advantage of the PAST algorithm and the repeated system matrices in the extended observability matrix.Te recursive formulation is frst derived with great detail, and the user-defned 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) Te identifed 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) Te proposed method with a small forgetting factor can immediately detect the damage although the results identifed 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 identifcation with the original sampling rate and is suitable for the practical implementation of the SHM system.(v) Te computation time of the convention methods confrms the efciency of the proposed method, and the diference demonstrates that timeliness can be secured even with 100 sampling rate measurement, which is common in the SHM system.

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

Figure 1 :
Figure 1: Te fowchart for implementing RSI with the proposed method.

Figure 2 :
Figure 2: Te schematic diagram of the 10-story shear-type frame and the identifed modal parameters under the seismic excitation.(a) 10-story shear-type frame for numerical simulation.(b) Identifed the frequency over time.(c) Identifed the damping ratio (multiplied by mode number) over time.

Figure 5 :
Figure 5: Te computation time with various sampling rates.

Figure 9 (
b) demonstrates the side view and (frst foor) plan view of this building and the locations of all the accelerometers with their channel numbers; meanwhile, Figures9(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.

Figure 6 :Figure 7 :Figure 8 :
Figure 6: Te photograph of the full-scale specimen and the close-up of the lock-up system.(a) Full-scale specimen.(b) Lock-up system.

Figure 9 :
Figure 9: Te 7-story RC building excited by the Chi-chi earthquake on September 21, 1999.(a) Measured acceleration responses from the roof in the longitudinal direction.(b) Side view and (frst foor) plan view of the feld frame installed with an ofine SHM system.(c) Measured acceleration responses from the basement in the longitudinal direction.

Figure 10 :
Figure 10: Te identifcation result of the RC frame during the earthquake excitation.

2
Structural Control and Health Monitoring where Y p and U p are past output and input data Hankel matrices; Y f and U f are future output and input data Hankel matrices, respectively.All those Hankel matrices have i rows and j columns (both in samples) as , respectively.As a result, the dimension of Y p and Y f is im by j, and the dimension of U p and U f is il by j.Similarly, W p , W f , V p , and V f are the Hankel matrices of the process and measurement noise in the matrix input-output equations.Furthermore, X p and X f are past and future state sequences (in column); Δ i is the reversed extended controllability matrix; H i and G i are the lower block triangular Toeplitz matrices; and Γ i is the extended observability matrix.
• • • y T i+j   T ; U p and U f are defned similarly as Y p and Y f

)
Structural Control and Health Monitoring

Table 1 :
Te modal frequencies and damping ratios of the numerically simulated structure.

Table 2 :
Te computation time over various dividing factors.

Table 3 :
Te computation time over various sampling rates.

Table 4 :
Te computation time of conventional methods.