Output-only subspace identification of structural properties and unknown ground excitation for shear-beam buildings

Abstract

A novel output-only subspace identification method is proposed to estimate the structural parameters of a shear-beam building under unknown ground excitation. The unique feature of the method is that both the structural parameters and the ground acceleration can be identified using only the absolute acceleration response, and no displacement and velocity responses are required. The method uses two subspace identification techniques sequentially and iteratively. The first technique extracts modal parameters of the building (natural frequencies, damping ratios, and mode shapes) using the absolute acceleration response and assumed or identified ground acceleration. These modal parameters further lead to the estimation of the building’s structural parameters when the mass properties of the building are known. The second technique then estimates the ground acceleration using the absolute acceleration and the obtained structural parameters from the first technique. The two techniques are performed iteratively until the estimated structural parameters converge. A numerical example and a laboratory test are used to illustrate the proposed identification method. Results show that the estimation of structural parameters is fairly robust to the presence of measurement noise while the unknown ground acceleration is not in high frequency range. Factors influencing the performance of the technique are also studied and discussed.

Introduction

Identifying modal and structural parameters for civil engineering structures is very important for their design verification, performance estimation, condition monitoring, damage detection, and safety assessment. This issue has been intensively researched and the solutions are well established when both the input and the output of a system are known.^1–3 However, it is quite common in practice that the input may not be available due to either unmeasured or immeasurable condition. Different methods^4–6 have been developed to estimate the unknown input with knowing all structural properties. However, in the field of structural identification, the structural properties are also unknown. Hence, it is necessary to develop techniques that can identify modal and structural parameters using output data only.

One common technique for output-only structural identification is to assume that the unknown inputs have a certain known characteristics such as white noise process.^7–9 Although this approach has seen many successful applications, significant errors seem unavoidable when the assumed characteristics of the inputs are not actually met. Another approach is to express the unknown input as a superposition of some orthogonal functions. ¹⁰ The structural parameters and the orthogonal coefficients are then identified simultaneously. However, the technique requires that the type and the number of the orthogonal functions be specified and might lead to an ill-posed problem if not considered carefully. The third approach estimates the structural parameters and the unknown input simultaneously without using any prior information of the input. Wang and Haldar^11,12 developed an iterative identification method for structures with viscous damping, in which the structural parameters were obtained using a least-squares method, while the unknown input was estimated from the identified structural parameters in each iteration. Note that the presence of proportional damping in this type of technique would lead to a nonlinear estimation problem which increases the complexity of computation. To address this issue, Ling and Haldar ¹³ proposed a modified iterative least-squares algorithm where the Taylor series approximation was used to transform a set of nonlinear equations to a set of linear ones. Zhao et al. ¹⁴ demonstrated that the structural parameters of a shear-beam building subject to ground excitation could not be uniquely identified using only the absolute structural responses. A hybrid identification method combining the time-domain information with the modal information was then proposed to identify the structural parameters and reconstruct the seismic input. Yang et al. ¹⁵ developed a recursive least-squares estimation technique with unknown input to identify the structural parameters as well as the unmeasured excitations. A comparison of the necessary information among the above different output-only approaches and the proposed approach in this article is summarized in Table 1.

Table 1.

Summary of output-only structural identification approaches.

Ref.	Category	Main features	Estimate structural properties	Estimate unknown input
Caicedo et al., 7 Alicioglu and Lus 8 and Mustafa and Necati 9	Ambient vibration analysis	Assumption of certain known statistical characteristics. Significant estimation errors when this assumption is not met	P	O
Zhang et al. 10	Basis function expansion	Unknown input expanded by orthogonal functions. Difficult to choose the type and number of orthogonal functions	P	P
Wang and Haldar11,12 and Ling and Haldar 13	Iterative least-squares approach	The response quantities including displacement, velocity, and acceleration required	P	P
Zhao et al. 14	Hybrid identification	Combining the time-domain information with the modal information to estimate both structural parameters and input	P	P
Yang et al. 15	Recursive least-squares approach	Real-time estimation of both structural parameters and input. Critical for the settings in the recursive algorithms	P	P
This article	Output-only subspace approach	Fast calculation based on QR and SVD. Only absolute responses needed	P	P

SVD: singular value decomposition.

QR: a decomposition of a matrix into a product of an orthogonal matrix (Q) and an upper triangular matrix (R).

Note that most of the least-square-based techniques mentioned above require that the acceleration, velocity, and displacement of the structure be measured either directly or indirectly via further data processing. Although the acceleration response of a structure can be easily and accurately obtained, the accurate measurement or determination of its velocity and displacement responses remains a challenging task in practice. Significant estimation errors seem unavoidable when these velocity and displacement responses are obtained via integration of acceleration due to the presence of measurement noise. ¹⁶

In this study, a novel identification method is proposed for a shear-beam building under unknown ground excitation. Modal, structural parameters, and the ground acceleration can be identified using only the absolute acceleration responses of the building, and no displacement and velocity responses are required. The method uses two subspace identification techniques sequentially and iteratively. The first technique extracts modal parameters of the building, including the natural frequencies, the damping ratios, and the mode shapes, using the absolute acceleration of the building and assumed/identified ground acceleration. When the mass properties of the building are known, these modal parameters further lead to the estimation of the building’s structural parameters. The second technique then estimates the ground acceleration based on the absolute acceleration and the obtained structural parameters of the building. The two techniques are performed iteratively until the estimated modal/structural parameters converge. A numerical example and a laboratory test are used to illustrate the proposed technique along with some factors influencing the identification performance.

Subspace identification of modal and structural parameters

The discrete state-space formulation for a general dynamic system with p inputs, p outputs, and N state variables can be expressed as follows

$$\{\begin{array}{c}\mathbf{x}(k+1)=\mathbf{Ax}(k)+\mathbf{Bf}(k)+\mathbf{w}(k)\hfill \\ \multicolumn{1}{c}{\mathbf{y}(k)=\mathbf{Cx}(k)+\mathbf{Df}(k)+\mathbf{v}(k)}\end{array}$$

(1)

where $\mathbf{x}\in {\mathbb{R}}^{\mathit{N}\times 1}$ is the state vector; $\mathbf{f}\in {\mathbb{R}}^{\mathit{p}\times 1}$ is the deterministic input vector; $\mathbf{y}\in {\mathbb{R}}^{\mathit{q}\times 1}$ is the observed output vector; $\mathbf{A}\in {\mathbb{R}}^{\mathit{N}\times \mathit{N}}$ , $\mathbf{B}\in {\mathbb{R}}^{\mathit{N}\times \mathit{p}}$ , $\mathbf{C}\in {\mathbb{R}}^{\mathit{q}\times \mathit{N}}$ , and $\mathbf{D}\in {\mathbb{R}}^{\mathit{q}\times \mathit{p}}$ are the system matrices; $\mathbf{w}\in {\mathbb{R}}^{\mathit{N}\times 1}$ is the input disturbance vector and $\mathbf{v}\in {\mathbb{R}}^{q\times 1}$ is the measurement noise vector; and k is the instant in the multiple of the sampling interval $\mathrm{\Delta }t$ . Both w and v are assumed to be zero-mean white noise processes in this study. The subspace identification problem can be postulated as the determination of the four system matrices (A, B, C, and D) up to a similarity transform given that y and f are available. In this study, the multivariable output error state space (MOESP) method is adopted to solve this problem. ¹⁷ Formulate the following two input block Hankel matrices

$$\begin{array}{c}{\mathbf{F}}_{\mathit{p}}=\left[\begin{array}{cccc}\mathbf{f}(0)\hfill & \mathbf{f}(1)\hfill & \cdots \hfill & \mathbf{f}(j-1)\hfill \\ \multicolumn{1}{c}{\mathbf{f}(1)}& \mathbf{f}(2)\hfill & \cdots \hfill & \mathbf{f}(j)\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{\mathbf{f}(i-1)}& \mathbf{f}(i)\hfill & \cdots \hfill & \mathbf{f}(i+j-2)\hfill \end{array}\right]\hfill \\ {\mathbf{F}}_f=\left[\begin{array}{cccc}\mathbf{f}(i)\hfill & \mathbf{f}(i+1)\hfill & \cdots \hfill & \mathbf{f}(i+j-1)\hfill \\ \multicolumn{1}{c}{\mathbf{f}(i+1)}& \mathbf{f}(i+2)\hfill & \cdots \hfill & \mathbf{f}(i+j)\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{\mathbf{f}(2i-1)}& \mathbf{f}(2i)\hfill & \cdots \hfill & \mathbf{f}(2i+j-2)\hfill \end{array}\right]\hfill \end{array}$$

(2)

The number of block rows i is a user-defined quantity which should be at least larger than the maximum order of the system. ¹⁷ Theoretically, a larger i gives more stable and robust estimation of system information from the measured data if the data length is infinitely long. However, very large i may result in the over-smoothing issue when there are disturbance and measurement noises for the situation of data with finite length. A procedure for the choice of i is given in the section of illustrative examples. The parameter j is typically set to $s-2i+1$ (s is the total length of sampled data) which implies that all the given data are used. The two output block Hankel matrices ${\mathbf{Y}}_{\mathit{p}}$ and ${\mathbf{Y}}_{\mathit{f}}$ are defined similarly by replacing f with y in equation (2).

The extended observability matrix O is defined as follows

$$\mathbf{O}=[\begin{array}{cccc}{\mathbf{C}}^T\hfill & {(\mathbf{CA})}^T\hfill & \cdots \hfill & {(\mathbf{C}{\mathbf{A}}^{i-1})}^T\hfill \end{array}{]}^T$$

(3)

where $(.{)}^T$ denotes the transpose of a matrix. The column space of O can be retrieved by an oblique projection of the row space of ${\mathbf{Y}}_{\mathit{f}}$ along the row space of ${\mathbf{F}}_{\mathit{f}}$ on the row space of ${\mathbf{W}}_{\mathit{p}}=[\begin{array}{cc}{\mathit{F}}_{\mathit{p}}^{\mathit{T}}\hfill & {\mathit{Y}}_{\mathit{p}}^{\mathit{T}}\hfill \end{array}{]}^{\mathit{T}}$ . This oblique projection can be calculated by the LQ decomposition as follows

$$\left[\begin{array}{c}{\mathbf{F}}_f\hfill \\ \multicolumn{1}{c}{{\mathbf{W}}_p}\\ \multicolumn{1}{c}{{\mathbf{Y}}_f}\end{array}\right]=\left[\begin{array}{ccc}{\mathbf{L}}_{11}\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{{\mathbf{L}}_{21}}& {\mathbf{L}}_{22}\hfill & 0\hfill \\ \multicolumn{1}{c}{{\mathbf{L}}_{31}}& {\mathbf{L}}_{32}\hfill & {\mathbf{L}}_{33}\hfill \end{array}\right]\left[\begin{array}{c}{\mathit{Q}}_1^T\hfill \\ \multicolumn{1}{c}{{\mathit{Q}}_2^T}\\ \multicolumn{1}{c}{{\mathit{Q}}_3^T}\end{array}\right]$$

(4)

The column space of O is then only associated with ${\mathbf{L}}_{32}$ in equation (4). Since the system in equation (1) is generally assumed to be a minimal realization and O is of full column rank, the observation matrix O can be estimated by taking the singular value decomposition (SVD) of ${\mathbf{L}}_{32}$ and retaining the first N dominant singular vectors. Once O is obtained, A and C could be estimated based on the shifted structure of O as follows

$$\begin{array}{c}\mathbf{C}=\mathbf{O}(1:q,:)\hfill \\ \mathbf{A}=\mathbf{O}(1:M-q,:{)}^{†}\mathbf{O}(q+1:M,:)\hfill \end{array}$$

(5)

For simplicity, the MATLAB syntax is adopted in the above equations. Also, $(.{)}^{†}$ denotes the Moore–Penrose pseudo-inverse of a matrix. $M=q\times i$ is the row number of the matrix O. Note i should be chosen sufficiently large such that A can be estimated stably and accurately from equation (5). Once A and C are obtained, the modal parameters of the system, including its natural frequencies ${\omega }_r$ , damping ratios ${\zeta }_r$ , and mode shapes ${\mathit{\phi }}_r$ , could be obtained as follows ¹⁸

$${\omega }_r=\frac{|\ln ({\mu }_r)|}{\mathrm{\Delta }t};{\zeta }_r=-\frac{Re(\ln ({\mu }_r))}{|\ln ({\mu }_r)|};{\phi }_r=C{\psi }_r;{\phi }_r^{*}=C{\psi }_r^{*};r=1,2,\dots ,N/2$$

(6)

where ${\mu }_r$ and ${\mathit{\psi }}_r$ are the rth pair of eigenvalues and eigenvectors of the state transition matrix A, $\ln (·)$ and $\mathrm{Re}(·)$ are the natural logarithm and the real part of the argument, respectively, ${\mathit{\phi }}_r^{*}$ and ${\mathit{\psi }}_r^{*}$ are the corresponding complex conjugates of ${\mathit{\phi }}_r$ and ${\mathit{\psi }}_r$ .

As mentioned above, the subspace identification method can be used to extract the system matrices up to a similarity transform from inputs and outputs. These system matrices cannot be directly connected to the physical parameters of the structure unless additional information such as the structural model is provided. In this study, the following assumptions are adopted: (1) the structure is an n-degree-of-freedom (n-DOF) shear-beam building with story mass and stiffness coefficient of $m_r$ and $k_r(r=1,2,\dots ,n)$ , respectively, and the stiffness matrix is as follows

$$\mathbf{k}=\left[\begin{array}{ccccc}{k}_1+k_2\hfill & -k_2\hfill & 0\hfill & \cdots \hfill & 0\hfill \\ \multicolumn{1}{c}{-k_2}& k_2+k_3\hfill & -k_3\hfill & \cdots \hfill & 0\hfill \\ \multicolumn{1}{c}{0}& -k_3\hfill & k_3+k_4\hfill & \cdots \hfill & 0\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & ⋮\hfill & \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & \cdots \hfill & k_n\hfill \end{array}\right]$$

(7)

(2) The building mass matrix is assumed to be known. Typically, for shear-beam building, it can take the lumped diagonal form $\mathbf{m}=\mathrm{diag}([\begin{array}{cc}\begin{array}{cc}{m}_1\hfill & m_2\hfill \end{array}\hfill & \begin{array}{cc}\cdots \hfill & m_n\hfill \end{array}\hfill \end{array}])$ which is used in this article. However, consistent mass matrix should be used if the structure has specified distributed inertial mass. (3) The outputs are assumed to be the absolute acceleration responses measured at all n DOFs of the building. (4) The damping matrix c is assumed to be of Rayleigh type: $\mathbf{c}=a_1\mathbf{m}+a_2\mathbf{k}$ where $a_1$ and $a_2$ are the Rayleigh damping coefficients. The objective of structural identification is to determine the stiffness coefficients $k_r(r=1,2,\dots ,n)$ and the Rayleigh damping coefficients $a_1$ , $a_2$ based on the modal parameters obtained from the subspace technique and the known mass m. Note that the number of identified modes from the subspace technique could be equal to or smaller than the building’s DOFs $(N/2\le n)$ depending on whether all modes of the building are excited by the input. Since proportional damping is assumed here, the theoretical mode shapes are normal modes or real mode shapes. However, the eigenvectors of the identified A from the subspace approach is generally complex pairs as shown in equation (6) for general damping. In order to identify the structural parameters here, real mode shapes are needed and approximated from equation (6) as follows

$${\varphi }_r\approx \left|{\mathbf{\phi }}_{\mathit{r}}\right|\mathrm{sign}(\mathrm{Re}({\mathbf{\phi }}_r))$$

(8)

Consider first the following Eigen equation for typical structural dynamics problem

$$\mathbf{k}{\varphi }_{\mathit{r}}={\omega }_r^2\mathbf{m}{\varphi }_{\mathit{r}}$$

(9)

Equation (9) can be cast into a regression form ⁷

$${\mathbf{\Phi }}_{\mathit{r}}\mathbf{\theta }={\mathrm{Z}}_r$$

(10)

in which

$${\Phi }_r=\left[\begin{array}{cc}\begin{array}{cc}\begin{array}{c}\begin{array}{c}{\varphi }_{r,1}\\ 0\end{array}\\ \begin{array}{c}⋮\\ 0\end{array}\end{array}& \begin{array}{c}\begin{array}{c}{\varphi }_{r,1}-{\varphi }_{r,2}\\ {\varphi }_{r,2}-{\varphi }_{r,1}\end{array}\\ \begin{array}{c}⋮\\ 0\end{array}\end{array}\end{array}& \begin{array}{cc}\begin{array}{c}\begin{array}{c}\cdots \\ \cdots \end{array}\\ \begin{array}{c}⋮\\ \cdots \end{array}\end{array}& \begin{array}{c}\begin{array}{c}0\\ 0\end{array}\\ \begin{array}{c}⋮\\ {\varphi }_{r,n}-{\varphi }_{r,n-1}\end{array}\end{array}\end{array}\end{array}\right];\theta =\left[\begin{array}{c}\begin{array}{c}{k}_1\\ k_2\end{array}\\ \begin{array}{c}⋮\\ k_n\end{array}\end{array}\right];{{\rm Z}}_r=\left[\begin{array}{c}\begin{array}{c}{\omega }_r^2{m}_1{\varphi }_{r,1}\\ {\omega }_r^2{m}_2{\varphi }_{r,2}\end{array}\\ \begin{array}{c}⋮\\ {\omega }_r^2{m}_n{\varphi }_{r,n}\end{array}\end{array}\right]$$

(11)

The term ${\varphi }_{r,i}$ is the ith element of ${\varphi }_r$ . Assembling h $(h\le N/2)$ pairs of natural frequencies and mode shapes leads to the following equations

$$\begin{array}{c}\mathbf{\Phi }\mathbf{\theta }=\mathbf{Z}\hfill \\ \mathbf{\Phi }={\left[\begin{array}{cccc}{\mathbf{\Phi }}_1^T\hfill & {\mathbf{\Phi }}_2^T\hfill & \cdots \hfill & {\mathbf{\Phi }}_h^T\hfill \end{array}\right]}^T;\hfill \\ \mathbf{Z}={\left[\begin{array}{cccc}{\mathbf{Z}}_1^T\hfill & {\mathbf{Z}}_2^T\hfill & \cdots \hfill & {\mathbf{Z}}_h^T\hfill \end{array}\right]}^T\hfill \end{array}$$

(12)

Note that equation (12) contains $h\times n$ equations and a least-squares technique can be used to estimate the n stiffness parameters.

It should be pointed out that if all the modes of the structure can be identified, the full stiffness matrix $\mathbf{k}$ can be calculated directly as follows

$$\mathbf{k}=\mathbf{m}\mathbf{\Phi }\mathbf{\Lambda }{\mathbf{\Phi }}^{-1}$$

(13)

where $\mathbf{\Lambda }=\mathrm{diag}([\begin{array}{cc}\begin{array}{cc}{\omega }_1^2\hfill & {\omega }_2^2\hfill \end{array}\hfill & \begin{array}{cc}\cdots \hfill & {\omega }_n^2\hfill \end{array}\hfill \end{array}])$ . If it is the case, the assumption of shear-beam building aforementioned can be relieved.

It should be pointed out that since both equations (9) and (13) are function of mass m, error in the mass can have direct influence on the accuracy of the stiffness estimation.

The two Rayleigh damping coefficients $a_1$ and $a_2$ can be estimated based on the natural frequencies and the damping ratios. For the Rayleigh damping, the damping ratios can be expressed as follows ¹⁹

$$\frac{1}{2}\left[\begin{array}{cc}\frac{1}{{\omega }_r}\hfill & {\omega }_r\hfill \end{array}\right]\left[\begin{array}{c}{a}_1\hfill \\ \multicolumn{1}{c}{a_2}\end{array}\right]={\zeta }_r$$

(14)

Similar to the identification of stiffness coefficients, h groups of natural frequencies and damping ratios can be assembled, and the damping coefficients are then obtained by the least-squares technique.

Subspace identification of input ground acceleration

After obtaining the stiffness and the damping matrices in previous section, the system matrices for the building in state-space form like equation (1) take on the following explicit expressions

$$\begin{array}{c}\mathbf{A}=e^{{\mathbf{A}}_c\Delta t};\mathbf{B}=(\mathbf{A}-\mathbf{I}){\mathbf{A}}_c^{-1}{\mathbf{B}}_c;\hfill \\ \mathbf{C}=[\begin{array}{cc}-{\mathbf{m}}^{-1}\mathbf{k}\hfill & -{\mathbf{m}}^{-1}\mathbf{c}\hfill \end{array}];\mathbf{D}=[0]\hfill \\ {\mathbf{A}}_c=\left[\begin{array}{cc}0\hfill & \mathbf{I}\hfill \\ \multicolumn{1}{c}{-{\mathbf{m}}^{-1}\mathbf{k}}& -{\mathbf{m}}^{-1}\mathbf{c}\hfill \end{array}\right];{\mathbf{B}}_c=\left[\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{\mathbf{-}1}\end{array}\right]\hfill \end{array}$$

(15)

where 0 and 1 $\in {\mathbb{R}}^{\mathit{n}\times 1}$ are vectors containing elements of 0 and 1, respectively, I is an identity matrix of appropriate sizes, and $(.{)}^{-1}$ is the inverse of a matrix. Note that the sizes of system matrices of the physical building now are as follows: $\mathbf{A}\in {\mathbb{R}}^{2\mathit{n}\times 2\mathit{n}}$ , $\mathbf{B}\in {\mathbb{R}}^{2\mathit{n}\times 1}$ , $\mathbf{C}\in {\mathbb{R}}^{\mathit{n}\times 2\mathit{n}}$ , and $\mathbf{D}\in {\mathbb{R}}^{\mathit{n}\times 1}$ . Here, same symbols in equations (1) and (15) are used; however, the difference in values from the physical modeling and input–output data modeling should be noticed. Wang and Haldar ¹¹ showed that when the mass, stiffness, and damping matrices of a structure are available, the input ground acceleration can be estimated if acceleration, velocity, and displacement are all measured. In this section, it is to show that the input ground acceleration can be obtained using the absolute acceleration measurement only by the subspace formulation.

After substituting equation (15) into equation (1), it can formulate the block input–output relation as follows

$$\mathbf{Y}(k)=\mathbf{Ox}(k)+\mathbf{S}\mathbf{F}(k)+\mathbf{V}(k)$$

(16)

where

$$\begin{array}{c}\mathbf{Y}(k)={\left[\begin{array}{cccc}{\mathbf{y}}^T(k)\hfill & {\mathbf{y}}^T(k+1)\hfill & \cdots \hfill & {\mathbf{y}}^T(k+l-1)\hfill \end{array}\right]}^T\hfill \\ \mathbf{F}(k)={\left[\begin{array}{cccc}{\mathbf{f}}^T(k)\hfill & {\mathbf{f}}^T(k+1)\hfill & \cdots \hfill & {\mathbf{f}}^T(k+l-2)\hfill \end{array}\right]}^T\hfill \\ \mathbf{O}=\left[\begin{array}{c}\mathbf{C}\hfill \\ \multicolumn{1}{c}{\mathbf{CA}}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{\mathbf{C}{\mathbf{A}}^{\mathit{l}-1}}\end{array}\right];\mathbf{S}=\left[\begin{array}{cccc}0\hfill & 0\hfill & \cdots \hfill & 0\hfill \\ \multicolumn{1}{c}{\mathbf{CB}}& 0\hfill & \cdots \hfill & 0\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{\mathbf{C}{\mathbf{A}}^{\mathit{l}-2}\mathbf{B}}& \mathbf{C}{\mathbf{A}}^{\mathit{l}-3}\mathbf{B}\hfill & \cdots \hfill & \mathbf{CB}\hfill \end{array}\right]\hfill \end{array}$$

(17)

The rth block component of $\mathbf{V}(k)$ in equation (16) is as follows

$$V_r(k)=C{A}^{r-2}w(k)+C{A}^{r-3}w(k+1)+\cdots +Cw(k+r-2)+v(k+r-1)$$

(18)

Note i in equation (3) is replaced by l (the row number of the output vector) in equation (17) to highlight that they are not necessarily the same.

Assume matrices $\mathbf{O}\in {\mathbb{R}}^{\mathit{nl}\times 2\mathit{n}}$ and $\mathbf{S}\in {\mathbb{R}}^{\mathit{nl}\times (\mathit{l}-1)}$ meet the following condition

$$\mathbf{T}=[\begin{array}{cc}\mathbf{O}\hfill & \mathbf{S}\hfill \end{array}];\mathrm{Rank}(\mathbf{T})=2n+l-1$$

(19)

where $\mathrm{Rank}(·)$ denotes the rank of a matrix. Equation (19) suggests that the matrix $\mathbf{T}$ has a full column rank, which implies that $\mathbf{CB}$ must also have a full column rank. The assumption also implies the degree of freedom n must be greater than the number of unknown inputs, that is, $n>1$ when the input is the ground acceleration. Since the matrix $\mathbf{O}$ has a full column rank, a projection matrix can be defined which projects the column space of a matrix onto the orthogonal complement of the column space of the matrix $\mathbf{O}$

$$\mathbf{P}=\mathbf{I}-\mathbf{O}({\mathbf{O}}^T\mathbf{O}{)}^{-1}{\mathbf{O}}^T$$

(20)

Multiplying $\mathbf{P}$ on both sides of equation (16) yields

$$\mathbf{P}\mathbf{Y}(k)={\mathbf{S}}_p\mathbf{F}(k)+{\mathbf{V}}_p(k)$$

(21)

where

$${\mathbf{S}}_p=\mathbf{PS};{\mathbf{V}}_p(\mathit{k})=\mathbf{PV}(\mathit{k});\mathbf{PO}=0$$

(22)

Equation (19) ensures that the matrix ${\mathbf{S}}_p$ has a full column rank; hence, equation (21) can be solved by a least-square method when l is sufficiently large. In equation (16), assembling all the outputs into a block Hankel matrix as in equation (2), the input $\mathbf{f}(k)$ except the last sample can be obtained. Theoretically speaking, the identified input matrix must also be a Hankel matrix with the same skew-diagonals as in equation (2), which however is not the case in practice due to the presence of estimation error. Nonetheless, the ground acceleration can be estimated as the average of the skew diagonals of the estimated input matrix. Note that the above derivations are valid only when the input is earthquake ground excitation $(p=1)$ . It is shown in Appendix 1 that the absolute acceleration responses alone are not sufficient to determine a more general type of loading when $p>1$ . For these more general cases, it is necessary to use more measured information for the identification.

Sequential estimation of structural parameters and ground acceleration

The two subspace techniques shown above can be used sequentially to identify the structural parameters of a shear-beam building and the unknown ground acceleration. Note that for identifying the structural parameters, it is necessary to know the input and the output of the building as well as its mass matrix, while the output and the building’s system matrices are needed to estimate the ground acceleration. A sequential identification procedure with iteration is shown in Figure 1 to estimate both the structural parameters and the ground acceleration using only the absolute acceleration measurement of the building. To initialize the algorithm, the following parameters or values are needed: the system order N, the dimensions of the input/output Hankel matrices i, j, and l, and the initial guess of input ground acceleration (zero or a white noise much smaller than the response intensity). The system order N can be inferred from either the power spectral density (PSD) of the observed output or the singular values of the matrix ${\mathbf{L}}_{32}$ or the dimension of mass m. Note that a subjective judgment is often needed to determine N especially when the output contains unknown noise. ²⁰ As aforementioned, the row number of Hankel matrices i should be set to be larger than N. The column number of Hankel matrices j is set to s − 2i herein, as the last sample of the ground motion cannot be obtained as seen from equation (17) when using all the response output data y. The row number of the output vector l should satisfy equation (19). The initial guess of input ground acceleration can be taken to be zero. The first subspace technique is used to extract the modal parameters and then the structural parameters of the building when the mass matrix is known. The obtained structural parameters can be used to formulate the stiffness and damping matrices, hence the system matrices as shown in equation (15). Last, the input ground acceleration can then be determined using the second subspace technique. The two subspace techniques are implemented sequentially and iteratively until the identified modal/structural parameters converge.

Figure 1.

Flowchart of the proposed subspace identification method with iteration.

Numerical example

First, a numerical example of a three-story shear-beam building ¹¹ was used to demonstrate the proposed technique. The structural parameters of the building were $m_1=m_2=m_3=10\mathrm{kg}$ ; $k_1=3000\mathrm{N}/\mathrm{m}$ , $k_2=2000\mathrm{N}/\mathrm{m}$ , and $k_3=1000\mathrm{N}/\mathrm{m}$ . The Rayleigh damping coefficients were $a_1=0.586/\mathrm{s}$ and $a_2=1.41\times {10}^{-3}\mathrm{s}$ . The corresponding natural frequencies of the building were 1.03, 2.41, and 3.99 Hz and the damping ratios were 5.0%, 3.0%, and 2.9%, respectively. White noise ground excitation was considered first since the traditional subspace identification method was developed assuming the whiteness of the unknown input and measurement noise. A historical earthquake record (non-stationary input) was then applied on the building to testify the flexibility of the proposed subspace approach.

White noise ground excitation

A Gaussian white noise process with an intensity of 0.1 m²/s⁴ was used to excite the building. Absolute acceleration responses at all the three floors were assumed to be measured at a sampling frequency of 50 Hz for 100 s. Thus, the number of data for each absolute acceleration response is s = 5000. These measured responses were superimposed with uncorrelated white noises. Different noise levels in terms of root mean square (RMS) were considered in this study. One dataset of measured responses which were used for identification is shown in Figure 2. The average PSD of the measured responses is also plotted which shows that all three modes of the building were well excited. This result suggested that the system order N should be set as 6.

Figure 2.

Absolute acceleration responses (with 1% RMS noise) and their averaged PSD.

Known input case

The subspace identification algorithm was first tested with known ground excitation. Hence, only the subspace identification for modal and structural parameters was performed without the iteration process for this case. It is instructive to see how the row number of the Hankel matrix i could be chosen to obtain acceptable identification performance. As for the column number j, it was so determined that all the 100-s data were used for identification. Figure 3 shows the identified stiffness and Rayleigh damping coefficients for different values of i ranging from 3 to 200 using the measured data in Figure 2. It can be found that all the identified parameters are almost the same as their exact values while there are only slight fluctuations for the two Rayleigh damping coefficients for small values of i. These results indicate that the estimation of modal and structural parameters is insensitive to i as long as it is larger than the minimum value ( $i_{\min }=3$ ) for the known input case.

Figure 3.

Effect of i on the identification performance under known input (1% RMS noise).

Statistical analysis of the identified parameters using 5000 simulations of different randomly generated ground motions for $i=3,10,20,\mathrm{and}50$ each was further performed and the results under 1% RMS noise are listed in Table 2. It shows that all the estimated parameters are unbiased and the standard deviations (stds) of the corresponding stiffness coefficients for different i are quite close. However, a larger i gives smaller stds for the two damping coefficients, indicating that the identified damping results are more reliable for larger i. Since the subspace identification algorithm is linear on processing the data, the statistics for different noise levels could be easily deduced by linear scaling. The statistical results in Table 2 suggest that i should be chosen as large as possible for better estimation, which unavoidably would demand more computational efforts. One possible way for choosing i is based on the required confidence level. As shown in Table 2, $i=20$ could be chosen after taking into account both the estimation accuracy and computational efforts. Note that a rule of thumb approach was proposed in Reynders et al. ¹⁸ to choose the value of i by

Table 2.

Statistics of estimated parameters under known input for different i (1% RMS noise).

Quantity	Exact	$i=3$		$i=10$		$i=20$		$i=50$
		Mean	Std	Mean	Std	Mean	Std	Mean	Std
$k_1$	3000	3000	1.17	3000	0.92	3000	0.91	3000	0.91
$k_2$	2000	2000	0.62	2000	0.60	2000	0.60	2000	0.61
$k_3$	1000	1000	0.28	1000	0.24	1000	0.24	1000	0.25
$a_1(\%)$	58.64	58.62	1.33	58.64	0.096	58.64	0.037	58.64	0.027
a 2 (‰)	1.405	1.406	0.033	1.405	0.0027	1.405	0.0014	1.405	0.0011

RMS: root mean square.

$$i\ge \frac{1}{2f_0\Delta t}$$

(23)

where $f_0$ is the lowest frequency of interest. The value is 25 for the example here according to Reynders et al. ¹⁸ when using $f_0=1\mathrm{Hz}$ .

Unknown input case

It is expected that the identification for the unknown input case will have larger estimation error statistically than the known input case. Based on the known input case, $i=20$ was chosen for the unknown input case. The initial value of the input ground acceleration was assumed to be zero. Since the unknown input was also estimated using the subspace technique, the effect of the parameter l on the identification performance was studied. As aforementioned, increasing the value of l means using more output data to estimate the input, resulting in better response fitting. The response fitting performance is defined by the error ratio (ER) as

$$\mathrm{ER}\ge \frac{||\hat{\mathbf{y}}-\mathbf{y}||}{||\mathbf{y}||}$$

(24)

where $\hat{\mathbf{y}}$ is the predicted responses using the estimated ground acceleration for each iteration. ER could be used as an index to determine the overall convergence rate and response fitting performance during the identification process for the proper selection of parameter l. Figure 4 shows the response fitting performance (ER) after 100 iterations for l = 10, 20, 50, and 100, respectively. It can be seen that a larger l gives a better fitting performance and quicker convergence rate. The identification converges quickly with only a small number of iterations for $l=50$ and 100. Actually, $l=100$ gives the best fitting results because its converged prediction error is close to that of the measurement noise level (1% RMS). However, a larger l will also increase the computational effort. Therefore, l was chosen based on the fitting performance, convergence rate, and computational load. Here, $l=50$ was used in the following analyses for the identification of unknown input.

Figure 4.

Effect of l on the response fitting performance under unknown input (1% RMS noise).

The identified modal and structural parameters using the measured data in Figure 2 are shown in Figures 5 and 6 after 20 iterations. Identified results of both 1% and 5% RMS noises are plotted for comparison. It can be seen that the overall convergence among different parameters could be clearly observed and most parameters converge quickly with around 10 iterations and closely approach to their exact values. It is seen that even in the first iteration where the unknown excitation was assumed to be white noise and initially taken to be zero in the identification algorithm, the modal frequencies were estimated quite accurately since all the vibrational modes of building were well excited. Furthermore, although the estimation of the damping ratios is less accurate in the first iteration, they approach their exact values very quickly. Using the identified modal parameters and the mass matrix, the structural parameters of the building were calculated and shown in Figure 6. The estimation errors after 20 iteration for stiffness coefficients are 1.41%, 0.07%, and 0.05% for 1% RMS noise and 2.42%, 0.40%, and 0.24% for 5% RMS noise, respectively. The estimation errors for the Rayleigh coefficients are 0.61% and 2.67% for 1% RMS noise, and 1.32% and 5.39% for 5% RMS noise, respectively. These results indicate that the proposed technique is not significantly affected by the increase in measurement noise intensity from 1% to 5%.

Figure 5.

Iteratively identified natural frequencies and damping ratios for the three-story shear-beam building under white noise ground excitation.

Figure 6.

Iteratively identified structural parameters for the three-story shear-beam building under white noise ground excitation.

The identified structural parameters were iteratively used to identify the ground acceleration. The identified ground acceleration time history after 20 iterations and its PSD are shown in Figures 7 and 8 for 1% RMS and 5% RMS noises, respectively. Note that only 5 s of data is plotted instead of the entire 100 s such that a clear comparison can be made. It is seen that the identified ground acceleration in time domain matches very well with the exact ground acceleration for 1% RMS noise, while large deviations are observed for 5% RMS. Figure 7 also shows that the PSDs of identified and exact ground acceleration are close in the low-frequency region and different slightly in the high-frequency region. This indicates that the input estimation is less accurate in the high-frequency region. This phenomenon can be observed more vividly for 5% RMS in Figure 8. It is seen that the measurement noise intensity can have a significant influence on the accuracy of ground acceleration estimation especially in the high-frequency region. However, it is found that by estimating also the input, the accuracy of structural properties could be significantly improved as seen from the iteratively estimated parameters in Figures 5 and 6.

Figure 7.

Identified ground acceleration and its PSD (1% RMS noise).

Figure 8.

Identified ground acceleration and its PSD (5% RMS noise).

Statistical analysis of the identified parameters for the unknown input case was further performed using 5000 simulations of different randomly generated ground motions. Table 3 shows the statistics of the identified parameters under 1% RMS and 5% RMS noises, respectively. In the statistical analysis, more number of iterations, 200 for each simulation was chosen to study the different convergence trends of the parameters. It can be seen that all the parameters are approximately unbiased. However, it should be pointed out that the stds and the coefficients of variation (covs) of the stiffness coefficient $k_1$ are the largest among the three stiffness coefficients for both noise levels. This shows that $k_1$ is less accurate when the input ground motion is unknown. The reason is that the first story restoring force from $k_1$ is due to the response difference between the ground and first floor. Thus, $k_1$ has a coupling effect with the ground motion during estimation, resulting in higher uncertainty than other stiffness coefficients. Also, the uncertainty of parameter ${\alpha }_2$ is the highest among all the parameters with cov of 2.48% and 12.60% for 1% RMS and 5% RMS, respectively. The result of high uncertainty in damping estimation under high noise level is quite typical for structural identification algorithms. In order to study the statistical trends of the identified parameters, the means and the 1 − σ bounds of the stiffness coefficients versus the iteration number for the two noise levels are plotted in Figure 9. It is seen that $k_2$ and $k_3$ converge faster than $k_1$ , Figure 9 shows that the estimation of $k_3$ is the most reliable, while $k_1$ is the worst. However, the estimation of $k_1$ is still reasonably accurate with covs under 2% and 6% for 1% RMS and 5% RMS, respectively. Figure 9 also shows that the convergence curve of std of $k_1$ is narrowing for 1% RMS noise, while expanding for 5% RMS noise. This indicates that the estimation of $k_1$ is highly affected due to the unknown input comparing with other parameters.

Table 3.

Statistics of estimated parameters after 200 iterations under unknown input.

Quantity	Exact	1% RMS			5% RMS
		Mean	Std	cov %	Mean	Std	cov %
$k_1$	3000	3004.7	34.85	1.16	3020.5	170.02	5.63
$k_2$	2000	2000.0	0.86	0.043	1999.1	4.43	0.22
$k_3$	1000	1000.0	0.24	0.012	1000.0	1.22	0.12
$a_1(\%)$	58.64	58.68	0.45	0.77	58.95	2.33	3.95
a 2 (‰)	1.405	1.403	0.035	2.48	1.402	0.18	12.60

RMS: root mean square; cov: coefficient of variation.

Figure 9.

Mean values and 1 − σ bounds of the estimated structural parameters under (a) 1% RMS noise and (b) 5% RMS noise.

Historical earthquake ground excitation

As demonstrated above, the proposed subspace identification technique is quite effective when the ground motion is white noise. In the following, the applicability of the proposed identification technique under non-white historical earthquake records is studied. The 1940 El Centro earthquake was used to excite the building and the absolute acceleration responses with a sampling frequency of 50 Hz were assumed to be measured output. High levels of measurement noises with 5% and 10% RMS were added to study the robustness of the technique. Following the previous analysis results, the parameters for the subspace identification technique were set as $N=6$ , $i=20$ , and $l=50$ . Figures 10 and 11 show the iterative results of the identified modal parameters and structural parameters, respectively. It is seen that the subspace technique provided fairly good initial estimates although the input was not white noise process. Furthermore, converged results could be obtained in a few iterations. The errors of modal frequencies and damping ratios are within 0.5% and 5%, respectively, for 10% RMS noise. The identified results after 20 iterations are $k_1=2938.8\mathrm{N}/\mathrm{m}$ , $k_2=2001.2\mathrm{N}/\mathrm{m}$ , $k_3=999.5\mathrm{N}/\mathrm{m}$ , $a_1=0.565/\mathrm{s}$ , and $a_2=1.44\times {10}^{-3}\mathrm{s}$ , and $k_1=2915.9\mathrm{N}/\mathrm{m}$ , $k_2=2000.6\mathrm{N}/\mathrm{m}$ , $k_3=998.9\mathrm{N}/\mathrm{m}$ , $a_1=0.549/\mathrm{s}$ , and $a_2=1.47\times {10}^{-3}\mathrm{s}$ , for 5% RMS and 10% RMS, respectively. The maximum errors for the first story stiffness are 2.0% and 2.8% and for the damping coefficient $a_1$ are 3.6% and 6.4%, under 5% RMS and 10% RMS, respectively. Note that the estimated stiffness coefficients for $k_2$ and $k_3$ are almost the same as their respective exact values. The slightly higher error in the first story is again due to the coupling effect between $k_1$ and the ground excitation. Figure 12 shows the PSDs of the identified input ground acceleration time history. From the PSD plots of the identified and the exact ground acceleration, it is noted that the low-frequency regions of the PSDs match quite well and the discrepancies are mainly seen in the high-frequency region which might be due to the high measurement noise level. In spite of this, the good accuracy of the identified results demonstrates that with the inclusion of input estimation, much better estimation could be achieved and the result is not sensitive to the measurement noise.

Figure 10.

Iteratively identified natural frequencies and damping ratios for the three-story shear-beam building under the El Centro earthquake excitation.

Figure 11.

Iteratively identified structural parameters for the three-story shear-beam building under the El Centro earthquake excitation.

Figure 12.

PSD of the identified El Centro ground acceleration.

Experimental study

An experimental study was performed on a three-story shear-beam building model as shown in Figure 13 to further validate the proposed technique. The building was made of aluminum and its dimensions are $401\times 314\times 1158\mathrm{mm}$ $(\mathrm{width}\times \mathrm{depth}\times \mathrm{height})$ . The mass coefficients were $m_1=5.63\mathrm{kg}$ , $m_2=6.03\mathrm{kg}$ , and $m_3=4.66\mathrm{kg}$ . The stiffness coefficients were obtained from a static test as $k_1=20.88\mathrm{kN}/\mathrm{m}$ , $k_2=22.37\mathrm{kN}/\mathrm{m}$ , and $k_3=24.21\mathrm{kN}/\mathrm{m}$ . The natural frequencies were calculated as 4.54, 13.02, and 18.21 Hz. The damping ratios for the first two modes were obtained from processing the free vibration response through the Hilbert–Huang transform and estimated as 2.39% and 0.87%, respectively. ²¹ Assuming Rayleigh damping, the two damping coefficients were calculated as 1.36/s and $9.0\times {10}^{-6}\mathrm{s}$ , respectively. The building was fixed on a shake table and subjected to the 1999 Chi-Chi earthquake ground acceleration (TCU075 station). The absolute acceleration responses at all three stories were measured with a sampling frequency of 250 Hz. An accelerometer was also installed on the shake table to measure the input ground acceleration for comparison. Before identification, the data were decimated by 5 through resampling at 50 Hz. This operation is to reduce the number of data points for fast computation and make the identification more accurate in the considered frequency range 0–25 Hz. ²² Figure 14 shows the resampled acceleration responses and their averaged PSD. It is seen that all three modes were well excited under the Chi-Chi earthquake. The parameters for the technique were set as $N=6$ , $i=20$ , and $l=50$ . Same as in the numerical example, the initial estimate of the ground acceleration was assumed to be zero.

Figure 13.

The three-story building model and its simplified three-DOF shear-beam model.

Figure 14.

Absolute acceleration measurements and their averaged PSD for the three-story laboratory model.

Figure 15 presents the iterative results of the identified natural frequencies and stiffness parameters. The results from the static test are also presented for comparison. It is seen that the initial estimates of both the natural frequencies and the stiffness coefficients assuming white noise input are close to their respective values obtained from the static test. The three natural frequencies converged to 4.15, 12.75, and 18.64 Hz, and the three stiffness coefficients to 19.51, 25.39, and 23.42 kN/m very quickly. The maximum difference in the estimated and the static stiffness coefficients is 13.5% in $k_2$ . It is seen that the estimated frequencies are closer to the excited modes from the PSD in Figure 14. Hence, the estimated stiffness coefficients are more accurate than those from the static test. As for the damping, the three damping ratios were 1.17%, 0.50%, and 0.54% as shown in Figure 16 and the Rayleigh damping coefficients were estimated to be $a_1=0.58/\mathrm{s}$ and $a_2=4.42\times {10}^{-5}\mathrm{s}$ after 10 iterations, which differ significantly from those obtained from the Hilbert–Huang transform. ²¹ This might be due to the reason that the actual damping is not proportional damping or is very different from the Rayleigh damping. The results show that the damping coefficients are much more difficult to estimate comparing to the stiffness coefficients. This conclusion is consistent to the reported results in Bernal and Ussia ⁵ and Ceravolo and Abbiati. ²³ Finally, the identified ground acceleration is plotted in Figure 17 together with the measured ground acceleration. Also shown are their respective PSD curves. It is observed that although the magnitude of the identified ground acceleration matches with that of the actual Chi-Chi earthquake ground acceleration in the time domain; however, there is a large discrepancy in the frequency domain. This discrepancy might be from the coupling estimation effect between the ground excitation and the first story, also the high measurement noise level, which is already discussed in the numerical studies of previous section.

Figure 15.

Iteratively identified natural frequencies and stiffness parameters for the three-story laboratory model.

Figure 16.

Iteratively identified damping ratios for the three-story laboratory model.

Figure 17.

Identified ground acceleration and its PSD for the three-story laboratory model.

Concluding remarks

In this article, a novel approach to identify modal and structural parameters for a shear-beam building under unknown ground excitation is presented. The unique feature of the method is that both the structural parameters and the ground acceleration can be identified using only the absolute acceleration responses of the building. This method consists of two subspace techniques which are used sequentially and iteratively. The first technique extracts modal parameters of the building, including the natural frequencies, the damping ratios, and the mode shapes, using the absolute acceleration of the building and assumed/identified input forces. When the mass properties of the building are known, these modal parameters further lead to the estimation of the building’s structural parameters. The second technique then estimates the input ground acceleration based on the absolute acceleration and the obtained structural parameters of the building. The two techniques are performed iteratively until the estimated structural parameters converge through the response fitting error curve. A numerical example and a laboratory test are used to illustrate the proposed technique. Results show that the method could accurately identify most of the structural parameters and reasonably estimate the input ground acceleration in low frequency range using only the absolute acceleration of the building. Also, the identification performance of the structural parameters does not seem to be significantly affected by the presence of measurement noise and the type of input ground acceleration. Note that the performance of the method is affected by some pre-assigned parameters such as the dimension of Hankel matrix and the row number of the output vector. Some suggestions on the selection of these parameters are discussed and presented in this study.