Combining subspace approach and short time Fourier analysis for locating structural damage storeys

2.1. Subspace approach in time-frequency domain

The dynamic responses of a linear structure satisfy the equation of motion:

where $\mathbf{M}$, $\mathbf{C}$ and $\mathbf{K}$ are mass, damping and stiffness matrices, respectively; $\ddot{\mathbf{x}}$, $\dot{\mathbf{x}}$ and $\mathbf{x}$ are the acceleration, velocity, and displacement responses vectors of the system, and $\mathbf{f}$ is the input force vectors. Usually, not all degrees of freedoms of system are measured in a field experiment, for reasons of economy. Only some parts of $\ddot{\mathbf{x}}$ or $\dot{\mathbf{x}}$ are measured. Consequently, the measured response vector $\mathbf{y}$, which can be velocity or acceleration responses, can be described by the state-space model. Generally, the state-space model considered in the following will be presented as [5]:

where $\mathbf{A}$, $\mathbf{B}$, $\mathbf{E}$ and $\mathbf{D}$ are system matrices that relate to $\mathbf{M}$, $\mathbf{C}$ and $\mathbf{K}$; $\mathbf{z}$ is state-variable $\mathbf{z}={\left[\begin{array}{ll}{\mathbf{x}}^T& {\dot{\mathbf{x}}}^T\end{array}\right]}^T$.

From Eq. (2), one can construct:

In the frequency domain, one can arbitrarily remove uninteresting frequency-bands to reduce the effects of noise in accurately estimating the modal parameters of a linear time invariant system. Consequently, Fourier transform is further introduced into Eq. (3). Treating the columns of ${\mathbf{z}}_k$, ${\mathbf{y}}_k$ and ${\mathbf{f}}_k$ as vector functions and applying the Fourier transform to Eq. (3) yields:

Via Eq. (4), one can further construct:

where:

where the two sets $\mathcal{F}{\left\{\mathbf{y}\right\}}_{k,s}$ and $\mathcal{F}{\left\{\mathbf{f}\right\}}_{k,s}$ are the short-time Fourier transform coefficients of ${\mathbf{y}}_k$ and ${\mathbf{f}}_k$, respectively. From Eq. (5), the following relation can be established:

where:

For further details of subspace system identification approaches, see Huang and Lin [5], and the references therein. As a brief summary, the subspace approach procedure is as follows:

1) Define an orthogonal projection matrix ${\prod }_f^{\perp }=\mathbf{I}-{\widehat{\mathbf{F}}}_{k,N}^{\mathbf{T}}{\left({\widehat{\mathbf{F}}}_{k,N}{\widehat{\mathbf{F}}}_{k,N}^{\mathbf{T}}\right)}^{-1}{\widehat{\mathbf{F}}}_{k,N}$, onto the null-space of ${\widehat{\mathbf{F}}}_{k,N}$.

2) Introduce instrumental variables $\mathbf{P}={\left[\begin{array}{ll}{\widehat{\mathbf{F}}}_{p,N}^T& {\widehat{\mathbf{Y}}}_{p,N}^T\end{array}\right]}^T$.

3) Calculate weighting matrices ${\mathbf{W}}_r=\mathbf{I}$ and ${\mathbf{W}}_c={\left(\mathbf{P}{\mathrm{\Pi }}_f^{\perp }{\mathbf{P}}^T/N\right)}^{-1/2}$.

4) Define a matrix $\stackrel{-}{\mathbf{H}}={\mathbf{W}}_r{\widehat{\mathbf{Y}}}_{k,N}{\prod }_f^{\perp }\mathbf{P}{\mathbf{W}}_c/N$.

5) Apply singular value decomposition on $\stackrel{-}{\mathbf{H}}\cong {\mathbf{Q}}_{\stackrel{-}{n}}{\mathbf{D}}_{\stackrel{-}{n}}{\mathbf{V}}_{\stackrel{-}{n}}^T$.

6) One has ${\mathbf{\Gamma }}_{\alpha }={\stackrel{-}{\mathbf{\Gamma }}}_{\alpha }{\mathbf{T}}_{\stackrel{-}{n}}$, where ${\stackrel{-}{\mathbf{\Gamma }}}_{\alpha }={\mathbf{W}}_r^{-1}{\mathbf{Q}}_{\stackrel{-}{n}}{\mathbf{D}}_{\stackrel{-}{n}}^{-1/2}$; ${\mathbf{T}}_{\stackrel{-}{n}}={\mathbf{D}}_{\stackrel{-}{n}}^{-1/2}{\mathbf{V}}_{\stackrel{-}{n}}{\left({\widehat{\mathbf{Z}}}_{k,N}{\prod }_f^{\perp }\mathbf{P}{\mathbf{W}}_c/N\right)}^{-1/2}$.

7) Rewrite Eq. (5) as $\mathcal{F}{\left\{\mathbf{y}\right\}}_{k,s}={\stackrel{-}{\mathbf{\Gamma }}}_{\alpha }{\left\{\stackrel{-}{\mathbf{z}}\right\}}_k+{\mathbf{\Phi }}_{\alpha }\mathcal{F}{\left\{\mathbf{f}\right\}}_{k,s}$, where ${\left\{\stackrel{-}{\mathbf{z}}\right\}}_k={\mathbf{T}}_{\stackrel{-}{n}}\mathcal{F}{\left\{\mathbf{z}\right\}}_k$.

8) ${\stackrel{-}{\mathbf{\Gamma }}}_{\alpha }$ can be established by ${\stackrel{-}{\mathbf{\Gamma }}}_{\alpha }={\left[\begin{array}{llll}{\stackrel{-}{\mathbf{E}}}^T& {\left(\stackrel{-}{\mathbf{E}}\stackrel{-}{\mathbf{A}}\right)}^T& \cdots & {\left(\stackrel{-}{\mathbf{E}}{\stackrel{-}{\mathbf{A}}}^{s-1}\right)}^T\end{array}\right]}^T$.

9) Define:

10) Calculate $\stackrel{-}{\mathbf{A}}$ by solving linear equation ${\stackrel{-}{\mathbf{\Gamma }}}_{\alpha 2}={\stackrel{-}{\mathbf{\Gamma }}}_{\alpha 1}\stackrel{-}{\mathbf{A}}$.

2.2. Determination of modal parameters

When the equation of motion is expressed in terms of the space-state variable, it is well known that the dynamic characteristics of the structural system are determined by the eigenvalues and eigenvectors of $\mathbf{A}$ in Eq. (2a). However, $\mathbf{A}$ cannot be determined from the preceding derivation. Thus, to find the modal parameters, $\stackrel{-}{\mathbf{A}}$ is employed:

where ${\lambda }_j$ and ${\stackrel{-}{\lambda }}_j$ are the $j$th eigenvalue of $\mathbf{A}$ and $\stackrel{-}{\mathbf{A}}$, respectively, while ${\mathbf{\varphi }}_j$ and ${\stackrel{-}{\mathbf{\varphi }}}_j$ are the corresponding eigenvectors. Since ${\mathbf{T}}_{\stackrel{-}{n}}$ remains unknown, the eigenvectors of $\mathbf{A}$ from Eq. (9) cannot be evaluated. Nevertheless, usually, one only needs to determine the modal shape corresponding to the observed degrees of freedom, ${\mathbf{\varphi }}_{j,y}$. From the procedure in Section 2.1, the following relationship can be determined:

The eigenvalues are complex numbers. Let:

From Eq. (9), one has:

Then, the pseudo-undamped circular natural frequency and the modal damping ratio for the system are:

In a proportionally damped system, this is equivalent to its undamped circular natural frequency.

2.3. Damage assessment

Consider the response of a structure described by equation of motions such as Eq. (1). Assuming proportional damping, the orthogonality property of the mode shapes with respect to the mass and stiffness matrix leads to [24]:

in which ${\mathbf{M}}_D$ and ${\mathbf{K}}_D$ are the diagonal modal mass matrix and diagonal modal stiffness matrix, respectively. The equation used to express the square of modal frequencies as a matrix is:

Here, the unitary stiffness matrix were defined as:

The unitary stiffness matrix at the sensor locations constructed from measured data before and after damage are denoted as ${\widehat{\mathbf{K}}}_u$ and ${\widehat{\mathbf{K}}}_d$, respectively. Applying singular value decomposition (SVD) to ${\widehat{\mathbf{K}}}_u-{\widehat{\mathbf{K}}}_d$ yields:

The singular values $\mathbf{D}$ are the difference in potential energy between ${\widehat{\mathbf{K}}}_u$ and ${\widehat{\mathbf{K}}}_d$ under each state $\mathbf{V}$ and $\mathbf{W}$. The location of maximum relative displacement in each $\mathbf{V}$ and $\mathbf{W}$ that lead to the difference of potential energy implies the maximum difference between the damaged and undamaged systems. Therefore, the possible damaged floors can be identified by nonzero singular values and corresponding states $\mathbf{V}$ and $\mathbf{W}$. The damage index can be obtained from the following steps:

1) Find the location ($l_i$) with maximum relative displacement in ${\mathbf{V}}_i$ or ${\mathbf{W}}_i$_.

2) Define damage index ($DI$):

3.1. Three-storey steel frame

Shaking table tests are often employed in a laboratory to investigate the behaviors of structures under earthquake conditions. To generate a set of earthquake response data for a benchmark model of a 1/8-scaled three-storey steel structure, CV-NCTU (Department of Civil Engineering, National Chiao Tung University) performed a series of shaking table tests on three frames (Fig. 1). One steel frame was represented as “std”, where each floor weighed about 3.8 kg and each column had a cross-sectional area of 80 mm² and was 440 mm in height. The second frame and third frame, which were structurally irregular with respect to stiffness and denoted “R_k1”, were identical to “std” except that a stiffening central column was removed from the first story. The third frame, denoted as “R_k12”, was identical to “R_k1” except that the stiffening central column was removed from the second story. These frames were subjected to base excitations specified by records of real earthquakes such as the 1995 Kobe earthquake. The responses were sampled at 100 Hz. Fig. 1 depicts the acceleration responses of each floors in the weak axis direction for frame “std”, subject to 20 % of the strength of the Kobe earthquake.

Fig. 1A photo of frame “std” on shaking table and its responses

Table 1 summarizes the identified dynamic characteristics of the three frames obtained using the acceleration responses of all floors and the input excitation in the weak axis direction. The results were obtained from the responses of all floors at $t=$ 5-25 seconds. The natural frequencies of “R_k1” and “R_k12” differ significantly from those of “std”. As expected, frame “R_k1” has lower natural frequencies than frame “std”, whereas frame “R_k12” has the lowest natural frequencies in all frames. To evaluate the correlation of mode shapes obtained from different methods and frames, the index of modal assurance criterion (MAC) [15] and coordinate modal assurance criterion (COMAC) [16] was computed to indicate the correlation between any two mode shapes of interest, and which are defined as:

To compare the identified mode shapes of undamaged frames in “std” with damaged frames in “R_k1” and “R_k12”, ${\phi }_{ij,s}$ and ${\phi }_{ij,R}$ respesent the deformation of the $i$th degree of freedom in the $j\mathrm{t}\mathrm{h}$ identified mode shapes, respectively.

Fig. 2 compares the results obtained using the various methods. The frame “std” is treated as a reference structure. Since the MAC and COMAC values are always between zero and unity, the maximum values of the results from proposed approach were normalized to unity for comparison. The results denoted by “present” in Fig. 2 normalize the difference of potential energy between damaged and undamaged systems. Fig. 2 clearly demonstrates that the proposed approach is superior to MAC and COMAC in identifying the floor of the structure whose properties differ from those of the corresponding floor in the reference structure.

Fig. 2Normalized index

a) R_k1

b) R_k12

3.2. Eight-storey steel frame

To generate a set of earthquake response data for a benchmark model of an eight-storey steel frame, the National Center for Research on Earthquake Engineering (NCREE) in Taiwan performed a series of shaking table tests on this model (Fig. 3). The eight-storey steel frame analyzed in this test was 1.8 m long, 1.2 m wide and 8.5 m high. Lead plates were piled on each floor such that the total mass of the steel frame was approximately 4.519 tons. In the shaking table tests, two accelerometers and two linear displacement transducers were installed in the long-span direction on the two edges of each floor to measure acceleration and displacement responses of the floor, respectively. The frame was subjected to base excitations equivalent to the Chi-chi earthquake at different reduced levels, including 200 Gal and 1200 Gal of base excitation levels. Data were recorded at a sampling rate of 200 Hz. The acceleration responses of the base and all floors at $t=$ 5-35 seconds were used in calculating modal parameters of the frame. Fig. 4 shows the base excitations and the responses of all floors in the long-span direction of the frame, subjected to 1200 Gal, equivalent to the Chi-hi earthquake. The strain records shown in Fig. 5 for different levels of earthquake input confirm that the columns of the first to third floors yielded when the frame was subjected to 1200 Gal base excitation.

Fig. 3A photo of frame on shaking table

Fig. 6 shows the results of the identified modal parameters of eight modes for the steel frame under 200 Gal and 1200 Gal of the Chi-chi earthquake force. Since the frame is symmetric with eight storeys, it is expected to have eight modes in the long-span direction of the frame. Damage that is caused by a strong earthquake will change the natural frequencies of the structure.

Fig. 7 compares the results obtained using the various methods. When the frame was subjected to 200 Gal of the Chi-chi earthquake, no nonlinear behaviors were observed, and it was treated as a reference structure. Again, the maximum values of the results from proposed approach were normalized to unity for comparison. The results denoted by “present” in Fig. 7 normalize the difference of potential energy between damaged and undamaged systems. Fig. 7 clearly indicates that the stiffness of the first to third floors significantly changed when the frame was under 1200 Gal of the Chi-chi earthquake, and these floors were somewhat damaged. Notably, Fig. 5 shows the measurements obtained by strain gauges installed on the columns of the first to third stories, which also indicated that the columns yielded during such base excitation. The success of the proposed procedure, as revealed in capturing the differences in properties of experimental frames, demonstrates the practical applicability of this procedure on an actual building with symmetry.

Fig. 4Responses of frame under 1,200 Gal of the Chi-chi earthquake

Fig. 5Strain responses at first, second and third floor columns