Damage Detection of Structures Identified with Deterministic-Stochastic Models Using Seismic Data

A deterministic-stochastic subspace identification method is adopted and experimentally verified in this study to identify the equivalent single-input-multiple-output system parameters of the discrete-time state equation. The method of damage locating vector (DLV) is then considered for damage detection. A series of shaking table tests using a five-storey steel frame has been conducted. Both single and multiple damage conditions at various locations have been considered. In the system identification analysis, either full or partial observation conditions have been taken into account. It has been shown that the damaged stories can be identified from global responses of the structure to earthquakes if sufficiently observed. In addition to detecting damage(s) with respect to the intact structure, identification of new or extended damages of the as-damaged counterpart has also been studied. This study gives further insights into the scheme in terms of effectiveness, robustness, and limitation for damage localization of frame systems.


Introduction
Structural health monitoring of civil engineering structures has received considerable attention over the last two decades with the progress in digital signal processing and system identification techniques [1][2][3][4][5][6]. Damage to a structure alters its dynamic characteristics in terms of the modal frequencies, damping ratios, and mode shapes, so as the physical parameters in terms of the stiffness or flexibility matrices. Despite that damage of a structure might be revealed from the changes of frequencies, it is difficult to locate the damages only with this information. The mode shapes perhaps provide a preferable basis for damage detection as they are spatially specific and reflective of local structural behavior. However, information of mode shapes alone is not sensitive enough for damage localization from the global dynamic vibration of the structure. The physical parameters indeed are more useful as far as damage localization is concerned. The study by Zhao and DeWolf [7] comparing various vibrational parameters for damage detection of spring-mass systems concluded that the model flexibilities were more sensitive to structural damages than either the natural frequencies or mode shapes. Stiffness is intuitively the most direct physical parameter related to damages in structures. However, the sensitivity analysis developed for damage detection based on stiffness matrix requires an accurate analytical model of the intact structure as a reference. Unfortunately, obtaining such an analytical model is in itself a difficult task. Moreover, synthesis of the stiffness matrix requires contributions of the higher modes that are practically difficult to identify. In contrast, the flexibility matrix can be sufficiently synthesized with a limited number of the lower modes as the modal contribution to the flexibility matrix decreases with the square of the corresponding modal frequency. Therefore, flexibility-based techniques have been considered of great potential in damage localization of structures from the global dynamic responses.
Another branch on the development of SHM techniques based on wavelet transformation has also gained a great deal of progress recently. A review by Li et al. [8] gives an insight of the up-to-date innovations on the SHM of infrastructures. In this paper, related theories and methodologies including sensing technology, sensor placement, signal processing, data 2 The Scientific World Journal fusion, system identification, and damage detection have been thoroughly discussed. Li et al. [9] propose a novel wavelet approach integrated with an auto-regressive moving average (ARMA) model for the damage detection of structures with progressive damage in time. The methodology has been validated both numerically and experimentally through shaking table tests of a reinforced concrete frame. The concept of wavelet transformation has also been adopted by Li et al. [10] to develop a method for determining the critical incidence of the seismic wave interpreted via energy principles. Yi et al. [11] proposes an energy-based multistage structural damage diagnosis method by adopting the schemes of wavelet packet transform (WPT) and artificial neural network (ANN). Using the benchmark structure of the American Society of Civil Engineers (ASCE) as the target, the authors demonstrate via numerical simulations that the proposed method is able to detect structural damage of various extents. Yet, localization of damages requires other measures.
Among the damage localization techniques based on variation of the physical parameters as the structure deteriorates, the flexibility-based approaches have been shown to be very promising and computationally efficient. Pandey and Biswas [12] demonstrated that the damage locations of a wide-flange steel beam could be successfully identified by interrogating the change in the flexibility. This technique has been further extended for damage detection of the plane trusses [13]. The flexibility-based damage localization technique referred to as the method of damage locating vectors (DLV) has been advanced further by Bernal [14] and Bernal and Gunes [15]. This methodology has also been adopted for damage localization of space trusses or plates by Gao et al. [16,17] and Huynh et al. [18]. The concept of the DLV method is to identify the members with zero stress under some specific loading patterns, namely, the DLVs that span the null space of the change in flexibility matrix of the structure before and after the damage state. The loading vectors are obtained by performing the singular value decomposition (SVD) on the flexibility differential matrix. Those structural elements resulting with zero stresses (or internal forces) under the static loads of the DLVs are considered potentially damaged. The DLV technique is capable of identifying scenarios of multiple damages in the structure via a truncated modal basis without a predetermined reference model. This null space approach has recently been further extended by Bernal [19,20] as the dynamic damage locating vector (DDLV) theorem that connects changes in the transfer matrices to the spatial location of stiffness related damage. Bernal [21] revised the DLV method further in the context of linear statespace models to comply with the stochastic realization of the systems by acknowledging that the DLVs in the null space of the flexibility differential can be extracted without a priori knowledge of the flexibility matrix. This scheme is adopted in this study.
System identification problems are commonly concerned with determining, from the relation of input/output sequences, the parameters of either polynomial models (e.g., ARMA, auto-regressive model with external input (ARX)) using the prediction error method [1] or state-space models by subspace system identification (SSI) [22], numerical algorithms for subspace state space system identification (N4SID) [23] for deterministic/stochastic systems, or SRIM for deterministic systems [24]. In contrast to the prediction error method that identifies the system parameters by minimizing the errors between the observations and those predicted by the model in a least squares sense, the stochastic realization theory initiated by Akaike [25] and Faure [26] is not based on the concept of optimization and leads to noniterative, convergence-guaranteed and efficient numerical algorithms. The covariance matrix is first constructed from block Hankle matrix [27] of shifted process sequences. The state-space model is in turn realized from the observability and/or controllability matrix via singular value decomposition of the covariance matrix. The theory of covariance-driven subspace method has been unified by van Overschee and de Moor [23] for deterministic, stochastic, and combined systems by defining the estimated state sequences as the projection of input-output data. The projected state sequences turn out to be the outputs of non-steady-state Kalman filter banks. To facilitate implementation of the DLV method, system identification technique of subspace technique or SRIM that identifies the system matrix and output (observation) matrix of a state-space model is considered a preferable alternative. The method of system realization using information matrix (SRIM) proposed by Juang [24] is based on a deterministic state-space system. The equivalent system matrix is identified from the covariance matrix of the input and output signals. This simple and elegant approach works well for systems in response to transient excitation if the noise level is negligible. The performance degrades, however, as the noise becomes pronounced. On the contrary, the stochastic subspace identification (SSI) method based on a stochastic model [23] without knowing the input is less sensitive to noise and works well if the excitation is Gaussian white noise. This scheme is not sufficient, however, for systems under transient excitations such as earthquakes. Therefore, a mixed deterministic-stochastic model is considered more robust to transient systems with nonnegligible noise contamination.
In this study, the feasibility of DLV method for damage detection of frame structures using seismic response data is explored in association with the subspace identification algorithm developed by van Overschee and de Moor [23] for mixed deterministic-stochastic systems. As an effort to verify the robustness of the deterministic-stochastic model against noise, system identification analysis by the N4SID algorithm is conducted with various noise levels and compared with those obtained by the SRIM. The effectiveness of the DLV method for damage detection using system parameters identified from contaminated signals is further investigated. Moreover, a series of shaking table tests using a five-storey steel frame has been conducted in National Center for Research on Earthquake Engineering (NCREE), Taiwan. Damage condition is simulated by reducing the cross-sectional area of some of the columns at the bottom. Both single and combinations of multiple damage conditions at various locations have been considered. In the system identification analysis, either full or partial observation conditions have been taken into account. It has been shown that The Scientific World Journal 3 the damaged stories can be identified from global responses of the structure to earthquakes if sufficiently observed. In addition to detecting damage(s) with respect to the intact structure, identification of new or extended damages of the as-damaged (ill-conditioned) counterpart has also been studied. The proposed scheme proves to be effective. This study gives further insights into the scheme in terms of effectiveness, robustness, and limitation for damage localization of frame systems.

Deterministic-Stochastic Subspace
System Identification A deterministic-stochastic linear time-invariant system is represented in a discrete-time state-space model as where z ∈ 2 ×1 and y ∈ ×1 are, respectively, the state and output vectors and u ∈ ×1 is the input vector at time instant with , , and being, respectively, the dimension of the dynamic system, observational outputs, and external inputs. A ∈ 2 ×2 is the system matrix, B ∈ 2 × is the input influence matrix, C ∈ ×2 is the observation matrix, and D ∈ × is the direct transmission matrix. w ∈ 2 ×1 and k ∈ ×1 are the unmeasurable vector signals assumed to be zero-mean, stationary white noise vector sequences. [ ] is the Kronecker .
The state and output vectors are split into deterministic and stochastic parts as z +1 = z +1 + z +1 and y = y + y satisfying, respectively, the deterministic subsystem and the stochastic subsystem By recursive substitution into the state-space equations of consecutive shifted processes, it leads to [28] where is the output block Hankle matrix [27] of the past where is the number of steps to be considered in the Markov process and represents the length of the data to be included in the analysis; is the input block Hankle matrix of the past; is the stochastic output block Hankle matrix of the past; is the deterministic state matrix; is the observability matrix; is the controllability matrix; is the triangular Toeplits matrix [27], and 4 The Scientific World Journal The stochastic covariance equations of the subspace can be defined as in which P and 0 represent, respectively, the auto-covariance matrix of the state (z ) and output (y ) vectors and G is the covariance matrix of the state (z +1 ) and output (y ) vectors.
The block Toeplitz covariance matrix Λ 0| −1 of the stochastic output is defined as where and the block Toeplitz cross covariance matrix Λ |2 −1 is defined similarly as where and represent the data length. Equation (18) is the stochastic controllability matrix. Projecting the future output state matrix Y |2 −1 onto the input matrix U 0|2 −1 and the past output state matrix Y 0| −1 , the output-input block Hankel matrices are constructed as where (⋅) * denotes the pseudoinverse of the matrix (⋅). To further simplify (19) and (20), the deterministic subspace and stochastic subspace are utilized to define what follows: Assuming that the deterministic input u and the deterministic state z are independent of the stochastic output y , then (19) and (20) can be simplified as The Scientific World Journal where in which It follows immediately from (22) that where Γ * or Γ * −1 denotes the pseudoinverse of its counterpart. The sequencesẐ andẐ +1 proved to be the states of a bank of non-steady-state Kalman filters [23]. The non-steady-state Kalman filter stateẐ can be defined in a recursive form as [23] If we collect consecutive columns of the Kalman states in parallel, (27) can then be extended aŝ where V | = (k k +1 k +2 ⋅ ⋅ ⋅ k + −1 ). Combining (30) and (31), one gets Substituting (26) forẐ andẐ +1 into (32), it leads to in which By forcing the noise terms in (33) to be zero, the coefficient matrices may be resolved as . A numerically stable and efficient algorithm referred to as the N4SID devised by van Overschee and de Moor [22] is adopted in this study to solve for the system matrix.

N4SID Algorithm
To facilitate implementation of the numerical algorithm, van Overschee and de Moor [22] proposed a numerically stable procedure N4SID utilizing the RQ decomposition of the block Hankle matrix whereQQ = I andR is a lower triangular matrix. It can be expressed in the following partitioned form as The Scientific World Journal which can be written in a more condensed manner as with Now if we express Y |2 −1 = (R 21 R 22 )Q and ( U 0|2 −1 Y 0| −1 ) = (R 11 0)Q = R 11 (Q 1 ) , then (19) can be written as By the same token, (20) may also be written as It has been shown by van Overschee and de Moor [22] that the observability matrix Γ can be obtained from SVD of The rank is determined from the dominant singular values of the decomposition, and Γ can be determined as and Γ −1 can be extracted directly from Γ without the last rows. With Γ and Γ −1 determined, one can calculate Γ * Y and Γ * −1 Y +1 without difficulty and in turn solve (35).

Review of the DLV Method
Bernal and Gunes [15] proposed that the structure subjected to the damage locating vectors, L, would undergo the same deformation before and after the damaged state. This statement immediately leads to where D is the flexibility differential of the structure before and after being damaged. When rank (D ) < ( is the degree of freedom of the structure), the basis corresponding to the null space of D is the damage locating vectors, L, which can be derived from singular value decomposition of the flexibility differential of the structure before and after the damage state. Members with nearly zero stress under the loadings of DLVs are considered potentially damaged.
The flexibility matrix of the structure can be expressed with the system matrices of the continuous-time state-space representation as where A = ln(A)/Δ ∈ 2 ×2 is the continuous-time system matrix; C 0 = [I 0] ∈ ×2 ; H = [ C 0 C 0 A ] ∈ 2 ×2 ;D = C 0 A B = −M −1 (M being the mass matrix of the system).
With (46), the flexibility differential D can be expressed as The weighted stress index WSI is defined as where nsi ,1 is the normalized stress index of the th member or d.o.f. subjected to the th DLV. Member (or storey ) is considered seriously damaged when the normalized weighted stress index nWSI ≤ 0.1 in which nWSI /WSI ,max , whereas it is considered moderately damaged as 0.1 < nWSI ≤ 0.2.

Numerical Illustration of the Proposed Schemes
As an effort to verify the robustness of the deterministicstochastic model against noise, system identification analysis by the N4SID algorithm is conducted with various noise levels and compared with those obtained by the SRIM. A five-story diagonally braced shear frame ( Figure 1) with parameters shown in Table 1 is considered in the numerical example. The 1940 El Centro earthquake with peak ground acceleration (PGA) scaled to 0.1 g is considered as the input.
To simulate noisy conditions, white noises of various noiseto-signal ratios (NSR) are added to the dynamic response signals. The noise-to-signal ratio at the th storey, NSR , is defined as where RMS , and RMS , represent, respectively, the root mean squares (RMS) of the added white noise and the original response signals at the th storey. In the simulation, the noise levels at all floors are assumed to be of the same amount in each case with NSR = 0%, 5% and 10%, respectively.
To assess the accuracy of the identification results, an error index ( ) of the th mode is defined as where is the total number of modes, , = (̂, − , ) 2 /̂2 , in whicĥ, and , are, respectively, the th element of the th mode shape derived from the eigen analysis and system identification. Figure 2(a) shows the comparison of mode shapes identified by the N4SID algorithm at various NSR levels with the analytical solutions. When the signals are noise-free (NSR = 0%), all the identified mode shapes coincide with their corresponding analytical counterparts. Deviation of the mode shapes increases with the noise level and becomes more pronounced for higher modes as observed from Figure 2(a). This is also revealed quantitatively from the error index ( ) bar charts as illustrated in Figure 3(a). When the noise level achieves NSR = 10%, only the first three modes can be identified with acceptable accuracy. The results by the SRIM exhibit similar trends, as shown in Figures 2(b) and 3(b), but deviate in a larger extent from the analytical solutions. When the noise level achieves NSR = 10%, moreover, only the first mode is identified with fidelity. The N4SID algorithm proves more robust to noises than the SRIM, as expected.
In order to examine if the damage detection scheme works sufficiently when the dynamic response signals are contaminated with noises, numerical simulations are next conducted using the same analytical model. Damage conditions of the structure are simulated by removing some of the diagonal bracings. Both single and multiple damages have been considered, which includes Cases Dl, D2, D3, D4, D5, D13, D15, and D135 where the number(s) denote the damaged stories. Only the N4SID algorithm will be adopted and a noise level of NSR = 10% is considered. Analysis for the noise-free condition has also been carried out as a reference. The results are summarized and compared in Table 2. The shaded area in the table with the normalized weighted stress index nWSI < 0.2 indicates those identified   as potentially damaged. It is observed that in both cases with NSR = 0% and 10%, the damaged stories can be successfully identified without exception, despite that results of the noisefree condition are shown to be more sensitive in the sense of getting smaller nWSI values at the damaged locations. The proposed scheme proves effective even if there exists a certain level of noise and only the first three modes are accurately identified.

Experimental Verifications
As a further step in verifying the feasibility of the proposed scheme experimentally, a series of shaking table tests has been carried out in the National Center for Research on Earthquake Engineering (NCREE), Taiwan, using a benchmark model (Figure 4). Weighing 22.5 tons, the benchmark model consisting of steel I-beams ( 150×150×7×10) is a five-storey frame of 7.5 m in height and 3 m × 2 m in plane, as illustrated in Figure 5. Accelerometers have been implemented at each floor and the base to monitor the dynamic responses of the model to ground motion which together serve as the basis for system identification. In order to get more insight of the structural behavior, additional velocity meters, LVDT, and load cells have also been implemented in the tests, as illustrated in Figure 6. Only the acceleration information is utilized, however, in this study. Damage of the structure is simulated by cutting out a small portion of the flange near the bottom of the column(s), as shown in Figure 7. In order to sufficiently examine the damage at various extents, the damages were progressively enforced on one side of the frame from the first storey up to the third storey from one test to another. It is meant to represent a moderate damage condition as a single column is damaged in the same storey while representing a serious damage condition as two columns are damaged. In total, six damage conditions have been considered in the test program. The 1940 El Centro earthquake has been adopted as the input with the PGA scaled to 0.1 g.
The test results are analyzed under considerations of full observation (utilizing acceleration responses of all stories), partial observation (ignoring accelerations at some stories), and the ill-conditioned condition where the reference structure has been wounded in previous tests.

Full Observation.
The six damaged conditions simulated in the tests are designated as: M1: single column damaged at the first storey, representing a moderate damage condition of the first storey; S1: two columns damaged at the first storey, representing a serious damage condition of the first storey; S1M2: two columns damaged at the first storey and single column damaged at the second storey; S12: two columns damaged at both the first and second stories; S12M3: two columns damaged at both the first and second stories and single column damaged at the third storey; S123: two columns damaged from storey 1 to 3.
The assessment results of various damage conditions based on full observation data of the structure are summarized in Table 3 where the shaded area corresponds to those being screened out as potentially damaged stories. Take case M1 for example, only the normalized weighted stress index of 1F is below 0.1, indicating damage of the first storey. This agrees with the actual damage condition. And in case S1M2, the index is 0.04 for 1F and 0.12 for 2F, agreeing with the actual condition that the first storey is seriously damaged and the second storey is moderately damaged. It is evident that, under a full observation condition, the damaged location(s) are successfully identified, regardless of single or multiple damage conditions.

Partial
Observation. The 12 cases analyzed in a partial observation condition are designated below as follows: M1/135: single column damaged at the first storey, representing a moderate damage condition of the first storey with only the first, third, and fifth floors observed; S1/135: two columns damaged at the first storey, representing a serious damage condition of the first storey with only the first, third, and fifth floors observed; M1/124: single column damaged at the first storey with only the first, second, and fourth floors observed; S1/124: two columns damaged at the first storey with only the first, second, and fourth floors observed; S1M2/124: two columns damaged at the first storey and single column damaged at the second storey with only the first, second, and fourth floors observed; S1M2/125: two columns damaged at the first storey and single column damaged at the second storey with only the first, second, and fifth floors observed; S1M2/1235: two columns damaged at the first storey and single column damaged at the second storey without observing the fourth floor; S12/124: two columns damaged at both the first and second stories with only the first, second, and fourth floors observed; S12/125: two columns damaged at both the first and second stories with only the first, second, and fifth floors observed; S12/1235: two columns damaged at both the first and second stories without observing the fourth floor; The Scientific World Journal  S12M3/1235: two columns damaged at both the first and second stories and single column damaged at the third storey without observing the fourth floor; S123/1235: two columns damaged at the first, second and third stories without observing the fourth floor.
The assessment results of various damage conditions based on partial observation data of the structure are summarized in Table 4 where the shaded area corresponds to those being screened out as potentially damaged stories. The performance index (PMC) is similarly defined as in Table 3. The results show that the scheme with partial observation remains effective for single damage conditions, as in the cases of M1/135, S1/135, M1/124 and S1/124 where the index corresponding to the first storey is below 0.2. The scheme fails, however, to locate the damaged stories in multiple damage conditions except for case S12/1235 where the first two stories are seriously damaged and 4 out of 5 stories are observed.

Ill-Conditioned Structures.
At times the structural health monitoring system might be introduced after the target building has been previously damaged. It is of interest to verify if the scheme is able to identify new or extended damage(s) in an earthquake event of an ill-conditioned structure that has been previously damaged. The system identification analysis will be based on full observation data as it provides more  reliable structural information for damage assessment. The 10 cases considered for damage detection analysis of illconditioned structures are designated below as follows: S1/M1: structure seriously damaged at the first storey of the current state versus the one with its first storey moderately damaged earlier; S1M2/S1: structure damaged seriously at the first storey and moderately at the second storey of the current state versus the one with its first storey seriously damaged earlier; S12/S1M2: structure damaged seriously at both the first and second stories of the current state versus the one seriously damaged at its first storey and moderately at the second in advance; S12M3/S12: structure damaged seriously at both the first and second stories and moderately at the third storey of the current state versus the one with its both first and second stories seriously damaged earlier; S123/S12M3: structure seriously damaged at the first, second, and third stories of the current state versus the one seriously damaged at its both first and second stories and moderately at the third in advance; S12/M1: structure seriously damaged at both the first and second stories of the current state versus the one with its first storey moderately damaged earlier; S12M3/M1: structure seriously damaged at both the first and second stories and moderately at the third of the current state versus the one with its first storey moderately damaged earlier; S123/M1: structure seriously damaged at the first, second, and third stories of the current state versus the one with its first storey moderately damaged earlier; S12/S1: structure seriously damaged at both the first and second stories of the current state versus the one with its first storey seriously damaged earlier; S123/S1: structure seriously damaged at the first, second, and third stories of the current state versus the one with its first storey seriously damaged earlier.
The assessment results of various damage conditions in respect to an ill-conditioned structure damaged earlier are summarized in Table 5. The performance index (PMC) is similarly defined as in Table 3. In all the cases considered, the scheme proves sufficient in identifying new or extended damage(s) without exception under a full observation condition in the system identification process, as indicated from the indices.

Summary and Conclusion
A scheme integrated with deterministic-stochastic subspace system identification and the method of damage localization vector for damage detection of structures is considered and verified based on seismic response data in this study. Both numerical simulation and experimental verification have been carried out to verify the feasibility of the proposed scheme. A series of shaking table tests using a five-storey steel frame has been conducted in National Center for Research on Earthquake Engineering (NCREE), Taiwan. Both single and combinations of multiple damage conditions at various locations have been considered. Either a full or partial observation conditions have been taken into account in the system identification process. This study gives further insights into the scheme in terms of effectiveness, robustness, and limitation for damage localization of frame systems. Based on the numerical simulation and test results, the conclusions are drawn as the following.
(a) The N4SID algorithm for deterministic-stochastic models is more robust to noise contamination than the SRIM algorithm for purely deterministic models. With the first three modes of the shear-type frame system accurately identified from noise-contaminated signals with 10% NSR by using the N4SID and the flexibility-based DLV method, all damaged conditions considered can be successfully predicted. (e) The scheme proves to be sufficient in identifying new or extended damage(s) without exception under a full observation condition.