Performance of tensor decomposition-based modal identification under nonstationary vibration

Natural hazards and extreme climatic events like earthquakes and strong wind introduce significant deterioration in aging civil structures. Human-induced vibration, on the other hand, is another form of nonstationary event that casues serviceability issues in flexible structures like pedestrian bridges. Structural health monitoring (Doebling et al 1996), a detailed diagnosis-based on input–output information of structure, plays a significant role in mitigating such catastrophic failures and predicting any future fatalities. Due to increasing size of modern civil engineering structures and cost of instrumentation, it is not feasible to measure inputs (i.e. earthquake or wind gusts) at all the locations. On the other hand, measuring dynamic response using inexpensive vibration sensors is proliferated. In the last few decades, operational modal analysis (also known as output-only modal identification) (Maia and Silva 2001) has been extensively utilized to estimate the modal parameters (i.e., natural frequencies, damping and modeshapes) based on vibration measurements. These methods are successfully used in a broad range of applications including model updating (Skolnik et al 2006), damage detection (Carden and Fanning 2004, Fan and Qiao 2011), serviceability assessment (Zivanovic et al 2005, Liu et al 2012), and vibration control (Moutinho et al 2011). In this paper, the authors have explored the performances of a relatively new tensor-decomposition based modal identification method under nonstationary response of structures.

Since early 1970s, various modal identification methods (Maia and Silva 2001) are developed for a multitude of structural and mechanical systems. Most of these methods operate either in the time or in frequency domain. The basic idea of the popular time-domain methods (Ibrahim and Mikulcik 1973, Juang and Pappa 1985, Allemang and Brown 1998, Perry and Koh 2000, Ma et al 2005) lies in extracting impulse response functions (IRFs) to develop state-space description of a dynamical system. There are other types of time-domain methods based on time-series models (Spiridonakos and Fassois 2014) focusing on identification of their model coefficients under both stationary and non-stationary environments. In fact, Stochastic subspace identification (SSI) method (VanOverschee and De Moor 1993, 1996) or Natural Excitation Technique (NExT) (James et al 1995) provides a powerful means of such parametric time-domain method. However, such techniques require the use of stability diagrams and modal assurance criteria (MAC) to select appropriate model order and stable modes. Unlike time-domain techniques, frequency-domain methods (Zhang et al 1985, Brincker et al 2001) are based on frequency response functions (FRFs) that are obtained using either peak picking (Zhang et al 1985) or frequency domain decomposition (FDD) (Brincker et al 2001). However, picking peaks is a subjective task, and the user-discrimination and issues associated with Fourier spectra are still present. Furthermore, most of the frequency and time-domain methods are based on the assumptions that the excitation is stationary. Therefore, these methods hinder their suitability towards modal identification under nonstationary excitations.

Recently several time-frequency and time-scale analysis methods employing wavelet transform (Amezquita-Sanchez and Adeli 2014, Dziedziech and Uhl 2015, Guo and Kareem 2015, Huang et al 2016), Laplace wavelet filtering and singular value decomposition (Guo and Kareem 2016), Hilbert transform (Huang et al 1998), empirical mode decomposition (EMD) (Darryll and Liming 2006, Sadhu 2015a) and blind source separation (BSS) techniques (Sadhu et al 2012, Sadhu 2013, 2015b, Sadhu and Narasimhan 2014) and time-series modeling based BSS methods (Musafere et al 2015) are explored towards modal identification under nonstationary excitations. These methods possess improved tracking, prediction and resolution capabilities to be readily used for nonstationary responses. However, each of these methods have their own limitations. For example, wavelet transform-based methods require a problem specific basis function, EMD is reliant on the intermittency criteria that is subjective in nature, whereas BSS methods (Sadhu et al 2012) require a proper selection of user-defined parameters depending on the extent of nonstationarity in the measurements. In parallel, tensor decomposition based methods (also known as parallel factor (PARAFAC) decomposition) have also garnered significant attention in the area of modal identification. In this paper, the performance of PARAFAC decomposition is explored under a wide range of nonstationary vibration.

In tensor decomposition methods, a covariance tensor is constructed from vibration measurements under multiple lags and is subsequently decomposed into covariances of hidden sources using multi-linear algebra tool such as alternating least squares (ALS) (Smilde et al 2004, Mokios et al 2006). This technique was introduced in two different independent works; canonical decomposition (Carroll and Chang 1970) and parallel factor analysis (Harshman 1970). Owing to the matricization operation among various lags, tensor decomposition-based methods are versatile and can be undertaken for both complete as well as partial measurement cases. Antoni and Chauhan (2011) used ALS to solve above tensor decomposition associated with modal analysis using a limited number of sensors (Papadimitriou et al 2000). The underdetermined source separation capability of PARAFAC decomposition was recently explored to identify the modal parameters of high-rise building (Mcneill 2012, Abazarsa et al 2013, Abazarsa et al 2015) and a structure equipped with tuned-mass damper (Sadhu et al 2014) using limited sensor measurements. In the case of complete measurement case, the modal responses can be obtained directly from the raw vibration measurements without any post-processing. Whereas, when only limited sensor measurements are available, PARAFAC decomposition results in the covariance of modal responses from which frequencies and damping are subsequently estimated.

Recently, PARAFAC decomposition was integrated with wavelet packet transform (WPT) to improve the source separation capability where mode-mixing in the WPT coefficients was alleviated using PARAFAC decomposition (Sadhu et al 2013, 2014). A pre-requisite in this method is the rank order selection a priori, which is an impediment for automated system identification. In order to alleviate this, a cluster diagram of modal frequencies is proposed under a suite of multiple rank orders (Sadhu et al 2015) from which the optimal rank order is chosen based on clustered densities of modal parameters. In this process, the effect of spurious modes is also circumvented and the modal frequencies are delineated from the excitation frequencies. However, none of the above studies investigates the performances of PARAFAC decomposition under nonstationary excitations. In this paper, the PARAFAC decomposition is studied under a wide range of nonstationary vibration due to earthquakes as well as human-induced excitations. The effect of lag parameters is investigated with respect to the extent of nonstationarity in the data. The performance of this method under limited sensor densities is also explored in a high-rise building model and a full-scale pedestrian bridge located in the campus of Lakehead University, Canada.

The paper is organized as follows. First, the basic motivation of the current research is presented in the introduction section. Background of PARAFAC decomposition is discussed next followed by the mathematical equivalence of tensor decomposition with modal identification. The proposed method is then validated using numerical and real-life case studies in the results section followed by key conclusions.

Any multi-dimensional signal can be expressed using a tensor (i.e., higher than two-dimensional) representation through multi-linear algebra tool which can be more effective than linear algebra tool (e.g. principal component analysis) (Bro 1997, Smilde et al 2004). A first-order tensor is essentially a vector ${\bf{s}}={s}_{i}\in {{\mathfrak{R}}}^{{n}_{1}}$, whereas a matrix ${\bf{S}}={s}_{{ij}}\in {{\mathfrak{R}}}^{{n}_{1}\times {n}_{2}}$ is a second-order tensor. Along the same line, a z-th order tensor is written as:

$$\begin{eqnarray}&&\bar{{\bf{S}}}={s}_{{ijk}...z}\in {{\mathfrak{R}}}^{{n}_{1}\times {n}_{2}\times {n}_{3}\times \cdots \times {n}_{z}}.\end{eqnarray} \tag{ 1 }$$

In general, a third-order tensor is primarily decomposed into a sum of outer products of triple vectors (Bro 1997):

$$\begin{eqnarray}&&\bar{{\bf{S}}}=\displaystyle \sum _{r=1}^{R}{{\boldsymbol{\theta }}}_{r}\circ {{\boldsymbol{\phi }}}_{r}\circ {{\boldsymbol{\psi }}}_{r}\,\Longleftrightarrow \,{\bar{S}}_{{ijk}}=\displaystyle \sum _{r=1}^{R}{\theta }_{{ir}}{\phi }_{{jr}}{\psi }_{{kr}},\end{eqnarray} \tag{ 2 }$$

where '$\circ$' denotes outer product with $i\in$ [1 I], $j\in$ [1 J] and $k\in$[1 K]. In equation (2), R is the number of rank-1 tensors present in $\bar{{\bf{S}}}$. This is also defined as trilinear model of $\bar{{\bf{S}}}$, $\bar{{\bf{S}}}=[{\boldsymbol{\Theta }},{\boldsymbol{\Phi }},{\boldsymbol{\Psi }}]$, where the matrices are given by ${\boldsymbol{\Theta }}=({{\boldsymbol{\theta }}}_{1},{{\boldsymbol{\theta }}}_{2},....,{{\boldsymbol{\theta }}}_{R})$, ${\boldsymbol{\Phi }}\ =({{\boldsymbol{\phi }}}_{1},{{\boldsymbol{\phi }}}_{2},....,{{\boldsymbol{\phi }}}_{R})$, and ${\boldsymbol{\Psi }}=({{\boldsymbol{\psi }}}_{1},{{\boldsymbol{\psi }}}_{2},....,{{\boldsymbol{\psi }}}_{R})$. As shown in equation (2), each triple vector product is a rank-1 tensor, namely PARAFAC component. Equation (2) represents the summation of R such PARAFAC components that fit the higher order tensor $\bar{{\bf{S}}}$ (Bro 1997, Lathauwer and Castaing 2008). The fundamental technique was developed by two different independent research: canonical decomposition (CANDECOMP) (Carroll and Chang 1970) and PARAllel FACtor (PARAFAC) analysis (Harshman 1970). The algorithms can be categorized in three main groups, namely (a) ALS (b) derivative-based methods and (c) direct or non-iterative approaches (Bro 1997). Out of three methods, the ALS method is the most popular method because of an easier implementation, smooth convergence and robust handling for higher order tensors. The key steps of the ALS are briefly presented in appendix A.1. The details of tensor decomposition can be found in the literature (Bro 1997) and are not repeated for brevity.

In order to demonstrate the signal separation capability of PARAFAC decomposition, following mixtures of three harmonics (sources, s) with frequencies 1.0, 2.5 and 1.2 Hz containing a measurement noise of $20 \%$ are considered as shown in the first two rows of figure 1:

$$\begin{eqnarray}&&\left[\begin{array}{c}{x}_{1}\\ {x}_{2}\\ {x}_{3}\end{array}\right] & = & \left[\begin{array}{ccc}1 & 1 & 1\\ 1 & 1.5 & 0.8\\ 2 & 2.2 & 1\end{array}\right]\,\left[\begin{array}{c}{s}_{1}\\ {s}_{2}\\ {s}_{3}\end{array}\right]\\ x & = & \theta s.\end{eqnarray} \tag{ 3 }$$

Since the mixtures contain three sources, rank-3 PARAFAC decomposition (i.e., R = 3) is performed over these mixtures to extract three hidden sources. The PARAFAC decomposition yields ${\hat{{\boldsymbol{\theta }}}}_{{\bf{1}}}$, ${\hat{{\boldsymbol{\theta }}}}_{{\bf{2}}}$ and ${\hat{{\boldsymbol{\theta }}}}_{{\bf{3}}}$ as shown in the last row of figure 1 from which the mixing matrix is estimated by concatenating successive normalized (say with respect to x₁) $\hat{{\boldsymbol{\theta }}}$'s that is very close to equation (3):

$$\begin{eqnarray}\hat{{\boldsymbol{\theta }}} & = & \left[\begin{array}{ccc}{\hat{\theta }}_{1} & {\hat{\theta }}_{2} & {\hat{\theta }}_{3}\end{array}\right]=\left[\begin{array}{ccc}0.93 & 0.81 & -0.97\\ 0.92 & 1.22 & -0.77\\ 1.82 & 1.79 & -0.92\end{array}\right]\\ & & \iff \hat{{\boldsymbol{\theta }}}=\left[\begin{array}{ccc}1.0 & 1.0 & 1.0\\ 0.99 & 1.5 & 0.79\\ 1.95 & 2.2 & 0.95\end{array}\right].\end{eqnarray} \tag{ 4 }$$

Zoom In Zoom Out Reset image size

A linear, proportionally damped, and discrete lumped-mass system with n_d degrees-of-freedom (DOF) when subjected to an excitation force, ${\bf{P}}(t)$ can be expressed as:

$$\begin{eqnarray}&&{\bf{M}}\ddot{{\bf{y}}}(t)+{\bf{C}}\dot{{\bf{y}}}(t)+{\bf{Ky}}(t)={\bf{P}}(t),\end{eqnarray} \tag{ 5 }$$

where, ${\bf{y}}(t)$ is a vector of displacement at the DOFs. ${\bf{M}}$, ${\bf{C}}$, and ${\bf{K}}$ are the mass, damping and stiffness matrices of the multi-DOF system. The solution of equation (5) can be expressed in terms of modal superposition of vibration modes with the following matrix form:

$$\begin{eqnarray}&&{\bf{y}}={\boldsymbol{\Gamma }}{\bf{q}},\end{eqnarray} \tag{ 6 }$$

where, ${\bf{y}}\in {{\mathfrak{R}}}^{{n}_{y}\times N}$ is a matrix consisting of measurements ${\bf{y}}$, ${\bf{q}}\in {{\mathfrak{R}}}^{{n}_{d}\times N}$ is a matrix of the corresponding modal responses, ${{\boldsymbol{\Gamma }}}_{{n}_{y}\times {n}_{d}}$ is the mode shape matrix, and N is the number of data points in the measurement. n_y is the number of measurement channels.

The covariance matrix ${{\bf{Z}}}_{{\bf{y}}}({\tau }_{k})$ of vibration measurements (${\bf{y}}$) evaluated at time-lag ${\tau }_{k}$ can be written as:

$$\begin{eqnarray}&&{{\bf{Z}}}_{{\bf{y}}}({\tau }_{k})=E\{{\bf{y}}(n){{\bf{y}}}^{{\mathsf{T}}}(n-{\tau }_{k})\}={\boldsymbol{\Gamma }}{{\bf{Z}}}_{{\bf{q}}}({\tau }_{k}){{\boldsymbol{\Gamma }}}^{{\mathsf{T}}},\end{eqnarray} \tag{ 7 }$$

where

$$\begin{eqnarray}&&{{\bf{Z}}}_{{\bf{q}}}({\tau }_{k})=E\{{\bf{q}}(n){{\bf{q}}}^{{\mathsf{T}}}(n-{\tau }_{k})\}\end{eqnarray} \tag{ 8 }$$

and q are the hidden sources. Let us consider the following annotations to simplify the mathematical notations of successive derivations

$$\begin{eqnarray}{Z}_{{y}_{1}{y}_{2}}({\tau }_{k}) & = & {Z}_{12k}^{y}\,\Longleftrightarrow \,{Z}_{{y}_{i}{y}_{j}}({\tau }_{k})={Z}_{{ijk}}^{y}\\ {Z}_{{q}_{1}}({\tau }_{k}) & = & {Z}_{k1}^{q}\,\Longleftrightarrow \,{Z}_{{q}_{l}}({\tau }_{k})={Z}_{{kl}}^{q}.\end{eqnarray} \tag{ 9 }$$

Considering a case with three available measurements where ${\bf{y}}=\{{y}_{1},{y}_{2},{y}_{3}\}$, equation (7) with above simplified notations can be represented as follows:

$$\begin{eqnarray}\left[\begin{array}{ccc}{Z}_{11k}^{y} & {Z}_{12k}^{y} & {Z}_{13k}^{y}\\ {Z}_{21k}^{y} & {Z}_{22k}^{y} & {Z}_{23k}^{y}\\ {Z}_{31k}^{y} & {Z}_{32k}^{y} & {Z}_{33k}^{y}\end{array}\right] & = & \left[\begin{array}{ccc}{{\rm{\Gamma }}}_{11} & {{\rm{\Gamma }}}_{12} & {{\rm{\Gamma }}}_{13}\\ {{\rm{\Gamma }}}_{21} & {{\rm{\Gamma }}}_{22} & {{\rm{\Gamma }}}_{23}\\ {{\rm{\Gamma }}}_{31} & {{\rm{\Gamma }}}_{32} & {{\rm{\Gamma }}}_{33}\end{array}\right]\\ & & \times \text{}\left[\begin{array}{ccc}{Z}_{k1}^{q} & 0 & 0\\ 0 & {Z}_{k2}^{q} & 0\\ 0 & 0 & {Z}_{k3}^{q}\end{array}\right]\,\left[\begin{array}{ccc}{{\rm{\Gamma }}}_{11} & {{\rm{\Gamma }}}_{21} & {{\rm{\Gamma }}}_{31}\\ {{\rm{\Gamma }}}_{12} & {{\rm{\Gamma }}}_{22} & {{\rm{\Gamma }}}_{32}\\ {{\rm{\Gamma }}}_{13} & {{\rm{\Gamma }}}_{23} & {{\rm{\Gamma }}}_{33}\end{array}\right].\end{eqnarray} \tag{ 10 }$$

Equation (10) can now be expressed as:

$$\begin{eqnarray}&&{Z}_{12k}^{y}={{\rm{\Gamma }}}_{11}{{\rm{\Gamma }}}_{21}{Z}_{k1}^{q}+{{\rm{\Gamma }}}_{12}{{\rm{\Gamma }}}_{22}{Z}_{k2}^{q}+{{\rm{\Gamma }}}_{13}{{\rm{\Gamma }}}_{23}{Z}_{k3}^{q}\end{eqnarray} \tag{ 11 }$$

which can be generalized for ${Z}_{{ijk}}^{y}$ of equation (10) as:

$$\begin{eqnarray}{Z}_{{ijk}}^{y} & = & {{\rm{\Gamma }}}_{i1}{{\rm{\Gamma }}}_{j1}{Z}_{k1}^{q}+{{\rm{\Gamma }}}_{i2}{{\rm{\Gamma }}}_{j2}{Z}_{k2}^{q}+{{\rm{\Gamma }}}_{i3}{{\rm{\Gamma }}}_{j3}{Z}_{k3}^{q}\\ & = & \displaystyle \sum _{r=1}^{3}{{\rm{\Gamma }}}_{{ir}}{{\rm{\Gamma }}}_{{jr}}{Z}_{{kr}}^{q}\quad i,j=1:3;k=1:K.\end{eqnarray} \tag{ 12 }$$

For any general n_d-DOF dynamical system, equation (12) can be simplified as:

$$\begin{eqnarray}&&{Z}_{{ijk}}^{y}=\displaystyle \sum _{r=1}^{{n}_{d}}{{\rm{\Gamma }}}_{{ir}}{{\rm{\Gamma }}}_{{jr}}{Z}_{{kr}}^{q}\,\Longleftrightarrow \,{{\bf{Z}}}^{{\bf{y}}}=\displaystyle \sum _{r=1}^{{n}_{d}}{{\boldsymbol{\Gamma }}}_{{\bf{r}}}\circ {{\boldsymbol{\Gamma }}}_{{\bf{r}}}\circ {{\bf{Z}}}_{{\bf{r}}}^{{\bf{q}}}.\end{eqnarray} \tag{ 13 }$$

Considering the similarity between equation (2) and (13), it is observed that by decomposing the third order tensor ${{\bf{Z}}}^{{\bf{y}}}$ into n_d number of PARAFAC components (i.e., modal responses), the mixing matrix (i.e., Γ) can be estimated. By using PARAFAC decomposition of ${{\bf{Z}}}^{{\bf{y}}}$, the resulting solutions yield the mixing matrix (i.e., mode shape matrix) ${\boldsymbol{\Gamma }}=[{{\boldsymbol{\Gamma }}}_{{\bf{1}}},{{\boldsymbol{\Gamma }}}_{{\bf{2}}},{{\boldsymbol{\Gamma }}}_{{\bf{3}}},\ldots ,{{\boldsymbol{\Gamma }}}_{{{\bf{n}}}_{{\bf{d}}}}]$ and the auto-correlation function ${{\bf{Z}}}_{{\bf{r}}}^{{\bf{q}}}$ for $r=1,2,3,\ldots ,{n}_{d}$ from which the natural frequency and damping of the individual modal responses can be estimated. In this way, the mathematical equivalence of PARAFAC decomposition and modal identification in the context of operational modal analysis is established (Sadhu et al 2013). However to the authors knowledge, the performance of PARAFAC decomposition has not been yet studied under nonstationary vibration. It may be noted that the lag parameter (i.e., K) plays an important role in separating the PARAFAC components. In this paper, the effect of K is studied under a wide range of nonstationary measurements obtained using more than 1400 ground motions and human-induced vibrations.

One of the salient feature of the PARAFAC decomposition is that, it can yield unique decomposition even if the rank is greater than the smallest dimension of the tensor (Lathauwer and Castaing 2008). This property of the PARAFAC decomposition can be utilized to solve underdetermined modal identification problems in structural health monitoring. Comparing equations (2) and (13), it is observed that equation (13) is a special case of PARAFAC tensor model with ${\boldsymbol{\phi }}={\boldsymbol{\theta }}$. Hence, a more relaxed uniqueness condition is proposed where the following inequality is satisfied (Lathauwer and Castaing 2008):

$$\begin{eqnarray}\displaystyle \frac{{n}_{d}({n}_{d}-1)}{2} & = & \displaystyle \frac{{n}_{y}({n}_{y}-1)}{4}[\displaystyle \frac{{n}_{y}({n}_{y}-1)}{2}+1]\\ & & -\displaystyle \frac{{n}_{y}!}{({n}_{y}-4)!4!}{({n}_{y})}_{\{{n}_{y}\geqslant 4\},}\end{eqnarray} \tag{ 14 }$$

where

$$\begin{eqnarray}{({n}_{y})}_{\{{n}_{y}\geqslant 4\}} & = & 0\quad {n}_{y}\lt 4\\ {({n}_{y})}_{\{{n}_{y}\geqslant 4\}} & = & 1\quad {n}_{y}\geqslant 4.\end{eqnarray} \tag{ 15 }$$

For a given number of measurements (n_y), an upper bound of source separability for PARAFAC decomposition can be computed using equation (14) which is tabulated in table 1 where ${n}_{d}^{u}$ is the highest number of PARAFAC components (i.e., modal responses) that can be extracted from n_y measurements. Table 1 shows that with rank-n_d PARAFAC decompositions of equation (13), n_d number of sources can be extracted from n_y measurements, when $2\leqslant {n}_{d}\leqslant {n}_{d}^{u}$. In this way, one can undertake a straight-forward approach to solve underdetermined modal identification problems, where n_d sources are identified from n_y vibration measurements even when ${n}_{y}\leqslant {n}_{d}$. In this paper, the effect of fewer number sensors (i.e n_y) and their locations are also studied under a wide range of nonstationary measurements.

Table 1.  Upper bound of source separation capability for PARAFAC decomposition.

ny	2	3	4	5	6	7	8	9	10
${n}_{d}^{u}$	2	4	6	10	15	20	26	33	41

The main contribution of this paper lies in investigating the performance of PARAFAC decomposition under nonstationary vibration first time. The PARAFAC decomposition technique is validated using a suite of more than 1400 ground motions and human induced vibration characterized by stationary duration. The effect of lag parameter is studied with respect to the severity of nonstationarity and the performance of the proposed method is evaluated using fewer number of sensors. Subsequently, the uncertainties associated with the accuracy of the PARAFAC decomposition using optimal sensor locations are also investigated.

In order to investigate the effects of lag parameter (i.e., K) and sensor densities (i.e., n_y) on modal identification, two different building models (5-storey and 10-storey models) and a wide range of earthquakes are considered.

4.1. 5-DOF building model and example ground motions

A dynamical system with 5 DOFs (Lin et al 1994) as shown in figure 2 is used to demonstrate the performance of PARAFAC method under base excitation. The natural frequencies of the model are 0.9, 3.4, 7.1, 10.7 and 12.7 Hz, respectively. The details of the model parameters can be found in appendix A.2. The model is subjected to a suite of ground motions and the resulting vibration responses are processed through PARAFAC method to extract the modal parameters under a wide range of nonstationary excitations.

Zoom In Zoom Out Reset image size

Table 2 shows six typical ground motions selected as input base excitations in the 5-DOF model along with the detailed information of peak ground acceleration (PGA) and ground motion duration (T). The sampling frequencies of all earthquakes are 50 Hz. The extent of nonstationarity is characterized by the ratio of stationary duration (T_s) (Trifunac and Brady 1975) and T. T_s (Trifunac and Brady 1975) may be computed using the time interval containing the energy envelope between 5% and 95% of the total energy of an earthquake. An earthquake is considered stationary when the fraction of time between 5% and 95% energy is close to 1. On the other hand, when this ratio attains lower value say, $\leqslant 0.3$, the earthquake (i.e., NR and PF) can be considered to be nonstationary. Figure 3 shows the Fourier spectra of the excitations revealing both wideband and narrowband characteristics with respect to the modal frequencies of the model (i.e., 0.9–12.7 Hz). For example, other than EC and KC earthquakes, the energies of the ground motions are distributed in a very narrow frequency range. Furthermore, the example excitations cover a broad range of PGA values with 0.01–0.37 g. Therefore, these ground motions form a perfect test bed to validate the PARAFAC method.

Zoom In Zoom Out Reset image size

Table 2.  Details of example ground motions.

Earthquake	PGA (g)	T (s)	$\tfrac{{T}_{{\rm{s}}}}{T}$
El Centro (EC), 1940	0.004	50.0	0.50
Northridge (NR), 1994	0.009	60.0	0.12
Imperial Valley (IV), 1940	0.36	53.8	0.47
Kern County (KC), 1952	0.16	54.4	0.63
Parkfield (PF), 1966	0.37	44.0	0.21
San Fernando (SF), 1971	0.02	68.7	0.67

Table 3.  Details of the suite of ground motion records used.

No.	Event	Magnitude	% of records
1	Imperial Valley, 1940	6.9	8.4
2	Santa Barbara, 1941	5.0	6.2
3	Eureka, 1954	4.5	6.7
4	San Francisco, 1957	5.3	6.2
5	Borrego Mountain, 1968	6.2	5.3
6	San Fernando, 1971	6.6	7.6
7	Oroville, 1975	5.7	6.7
8	Northern California, 1975	5.2	6.2
9	Coyote Lake, 1979	5.9	6.2
10	Imperial Valley, 1979	6.0	9.3
11	Mammoth Lakes, 1980	5.3	5.3
12	Livermore, 1980	5.5	7.6
13	Coalinga, 1983	6.5	4.4
14	Morgan Hill, 1984	6.1	7.1
15	Northridge, 1994	6.7	6.7

Figure 4 shows the Fourier spectra of top floor measurements of the building under example ground motions. As evident from figure 3 that most of the energies of the earthquakes are distributed within 0−8 Hz, the first three modes (i.e., 0.91, 3.4 and 7.1 Hz) of the 5-DOF building model are mostly excited and showed up in the vibration spectra. Therefore, the modal identification of 5-DOF model is restricted to only first three modes for brevity.

Zoom In Zoom Out Reset image size

4.2. Application of PARAFAC method

The floor vibration responses resulting from each of the earthquakes are processed through PARAFAC decomposition. In this section, all the floor measurements are utilized and covariance matrices with lag (K) up to 15 s are utilized. The resulting modal responses as obtained from PARAFAC method are shown in figures 5–7 for six ground motions. It may be noted except under SF earthquake, PARAFAC separates out the target modes under all ground motions. In fact, as shown in figures 3 and 4, SF earthquake can be considered as a very narrowband excitation (with a frequency content less than 3 Hz) that excites only few modal frequencies of the model and leads to poor performance in figure 7(b).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 8 shows the resulting modeshapes as obtained using the vibration data under EC earthquake resulting a MAC of more than 0.99 in all the target modes. Similar accuracy is obtained in all other five earthquakes and is not repeated. The results reveal that PARAFAC method can be considered as a versatile modal identification method for both stationary and nonstationary earthquake excitations.

Zoom In Zoom Out Reset image size

The effectiveness of PARAFAC method depends on the lag parameter and the number of sensors used. In order to evaluate the performance of PARAFAC under different lag parameters, a suite of different lags is undertaken. Figure 9 shows the performance of the identified modal parameters using three different lags which reveals its insensitivity towards its choice of lag parameters even under nonstationary vibration. Similar exercise is repeated for NR and KC earthquakes and the results are shown in figures 10 and 11 respectively. As shown in table 2, NR and KC have extreme $\tfrac{{T}_{{\rm{s}}}}{T}$ values ranging between 0.12 and 0.63, even though the PARAFAC method is successful in separating key modal frequencies under any choice of K. Due to tensor representation of covariance matrices under multiple discrete lags and its successive decomposition, the PARAFAC method shows excellent signal separation capabilities using nonstationary measurements. It further corroborates its suitability as a possible modal identification tool to any nonstationary response.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Performance of PARAFAC with limited sensors is shown in figures 12 and 13 for NR and KC earthquakes respectively. Results show that removing a sensor does not have a significant impact to the performance of PARAFAC method for low-rise buildings.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.3. Performance under a larger suite of ground motions

In previous section, 5-DOF model is investigated using six example ground motions. In this section, a database of 1443 accelerograms is considered to validate the performance of PARAFAC under a wide range of nonstationarity and sensor combinations. The accelerograms considered in the database have been recorded during 14 earthquake events in western USA between 1931 and 1984 (Trifunac and Brady 1975) and during the 1994 Northridge earthquake. Table 3 gives the names of these events along with their magnitudes and the number of records chosen from each event. The 1443 records have been chosen in such a way that the suite has a balanced distribution of records in terms of magnitude, epicentral distance, strong motion duration, and geologic site conditions. All records have PGAs greater than 0.1 g, and have magnitude M ranging from 4.5 to 6.9, epicentral distance R from 4 to 62 km, strong motion duration T_s from 1.8 to 42 s, and site conditions from alluvium to rock.

Table 4.  Number of combinations for different number of limited sensors.

No. of sensors	6	7	8	9	10
No. of sensor combinations	210	120	45	10	1

Figures 14 and 15 show the distribution of PGA and $\tfrac{{T}_{{\rm{s}}}}{T}$ ratios of the suite of ground motions. It reveals that the ground motions form a perfect test bed for the validation purpose with respect to two nonstationary measures (i.e., PGA and $\tfrac{{T}_{{\rm{s}}}}{T}$). Figure 16 shows the identification results of the PARAFAC method with the PGA values of earthquake under full sensor (i.e., ${n}_{y}=5$) and fewer sensor (i.e., ${n}_{y}=4$) case using five different sensor combinations (i.e., resulting ${C}_{4}^{5}=5$ sensor combinations), while figure 17 shows similar results w.r.t. $\tfrac{{T}_{{\rm{s}}}}{T}$ ratio. In both the figures, solid lines represent exact values of natural frequencies. It is seen that the identified frequencies are nearly clustered with the actual values even under multiple sensor combinations using fewer sensors. These results reveal that the performance of PARAFAC method is insensitive to the location of partial sensors. However it is observed that when $\tfrac{{T}_{{\rm{s}}}}{T}$ is more than 0.4 (i.e, with increasing stationarity of earthquake), the clusters of identification results are relatively sparse and associated with less uncertainties. Similar conclusions can be drawn under partial measurement combinations. Therefore, the PARAFAC method can be treated as a robust method under wide range of nonstationary excitations even with complete and partial measurements.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.4. 10-DOF building model

In order to show the effectiveness of limited sensors in the PARAFAC method, a 10 DOF model (Villaverde and Koyama 1993) is now utilized. The natural frequencies of the 10-DOF modal are 0.8, 1.8, 2.8, 3.9, 5.0, 6.1, 7.2, 8.5, 9.9 and 11.5 Hz respectively. Using EC, KC, PF and NR earthquakes as the base excitations, the performance of PARAFAC method in 10-DOF model is studied under a wide range of nonstationarities. As shown in table 1, the minimum number of sensors that can be used for a 10-DOF modal is 6. Therefore, unlike 5-DOF model (as in section 4.1), 10-DOF model can be effectively used to demonstrate the proposed method under fewer sensors cases owing to their various combinations.

Table 4 shows the number of sensor combinations for the number of limited sensors used. For example, the total number of sensor combinations using 6 sensors is 210. For a specific earthquake and given number of sensors, PARAFAC decomposition is performed for all possible sensor combinations and then the total number of target frequencies are identified. The statistics of the number of identified frequencies are shown in figures 18 and 19 under above mentioned four different earthquakes, respectively. The circle represents the mean number of frequencies detected. The top and bottom solid lines indicate the standard deviations ($\pm \sigma$), whereas the dotted lines show the range of the number of frequencies identified. The results reveal that the mean number of identified frequencies remain approximately the same irrespective of the number of fewer sensors, however the uncertainties associated with the identification become larger with reduced number of sensor cases. Figure 20 shows the coefficient of variation of the mean number of frequencies of figures 18 and 19 revealing significant accuracy even with fewer number of sensors. With these results, it can be concluded that the selection of optimal number of sensors will play a key role to identify the maximum number of target frequencies which however is beyond the scope of present study. Finally, figure 21 shows the Fourier spectra of identified frequencies using nine sensors under EC earthquake. The results of the identified frequencies under EC earthquake using 9 sensors are compared with the FE frequencies as well as identified frequencies obtained from the SSI method in table 5. It can be observed that the performance of PARAFAC method is relatively better than the SSI method.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 5.  Identification results of the 10-DOF model under EC earthquake.

	${\omega }_{1}$	${\omega }_{2}$	${\omega }_{3}$	${\omega }_{4}$	${\omega }_{5}$	${\omega }_{6}$	${\omega }_{7}$	${\omega }_{8}$	${\omega }_{9}$	${\omega }_{10}$
Exact	0.8	1.8	2.8	3.9	5.0	6.1	7.2	8.5	9.9	11.5
PARAFAC	0.8	1.8	2.8	3.8	5.0	6.0	7.3	—	9.9	—
SSI	0.65	1.8	3.2	4.1	5.0	5.5	—	—	8.5	—

5.1. Bridge description and instrumentations

In order to illustrate the proposed method, a footbridge (as shown in figure 22) crossing the Mcintyre River located in the campus of Lakehead University is utilized under a wide range of pedestrian-induced nonstationary excitations. Designed in 1967, it is composed of two main girders fixed into concrete abutments on both ends with steel struts spaced evenly along the length of the bridge and wooden lumber for decking as shown in figure 23.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

This bridge is instrumented with the accelerometers along the deck to measure the pedestrian-induced vibration. Figure 24(a) shows the layout of the sensor nodes used in this test. A total of eight sensors with a sensitivity of 1 V g^–1 are used with four distributed evenly on each side of the bridge. The sensors were attached to a data acquisition system operated in a portable computer as shown in figure 24(b). The sampling frequency was set to 200 Hz.

Zoom In Zoom Out Reset image size

5.2. Vibration testing and modal identification

The footbridge was subjected to excitation through a variety of activities such as walking, running, jogging and cycling. The different excitations used during the test are intended to represent the normal operational conditions of the bridge. This bridge mostly sees only light traffic but occasionally is subjected to periods of higher traffic specially during class hours. The excitation tests conducted are described in table 6 including the test duration and ${T}_{{\rm{s}}}/T$ at the mid span. As listed in the table, ${T}_{{\rm{s}}}/T$ indicates the extent of nonstationarity present in the vibration data with ${T}_{{\rm{s}}}/T$ ratio ranging between 0.5–0.8.

Table 6.  Details of vibration testing in the footbridge.

Test	Details	T (s)	$\tfrac{{T}_{{\rm{s}}}}{T}$
Single walking	A single person walked at a normal pace across the bridge	10.9	0.62
Group walking	A group of four people walked at a normal pace across the bridge	20.1	0.61
Single running	A single person ran at a jogging pace across the bridge	7.2	0.53
Group running	A group of four people ran at a jogging pace across the bridge	12.9	0.54
Single jumping	A single person jumped at the center of the bridge	12.5	0.52
Group jumping	A group of four people jumped at the center of the bridge	7.7	0.6
Biking	A single person rode a bicycle at a normal pace over the bridge	10.1	0.79

Figures 25 and 26 show the vibration response of the bridge under single and group walking as listed in table 6. Due to the nature of transient and spatial excitation over the bridge, it can be seen that the responses are nonstationary which are also reflected through the $\tfrac{{T}_{{\rm{s}}}}{T}$ ratio of table 6. For example, sudden jumping is always more nonstatioanry than the gentle walking which is revealed by smaller $\tfrac{{T}_{{\rm{S}}}}{T}$ ratio in the jumping data. Therefore, the entire data forms a perfect test bed of the proposed algorithm. In this paper, the vibration data under single running is utilized to validate the performance of the proposed algorithm. The acceleration data is first processed through the PARAFAC method and the resulting modal responses for single run and group walk data are shown in figure 27. In order to validate the accuracy of the results, a finite element (FE) model of the pedestrian bridge is developed using S-frame software as shown in figure 28. The identified frequencies (4.2, 11.8 and 28.5 Hz) are reasonably matching with the FE frequencies. Slight discrepancies may be observed due to inexact modeling and lack of structural details of the pedestrian bridge.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

This paper investigates the performance of the PARAFAC method in modal identification of structures under nonstationary excitations. A suite of ground motions under a wide range of nonstationarity is utilized to validate PARAFAC method in scalable building models and real-life footbridge structure. The results show the performance of lag parameters is insensitive to the extent of nonstationarity which also confirm the suitability of the PARAFAC method under wide range of nonstationary excitation. The effect of sensor densities is also investigated and the results indicate the number of mean target frequencies is insensitive to the choice of sensor combinations under limited sensor case which is an attractive feature of PARAFAC method. However the uncertainties associated with the identification becomes larger for limited sensor case which leads to further studies related to optimal sensor configurations. The validation using real-life footbridge data also corroborates with the accuracy of the PARAFAC decomposition method under human-induced nonstationary excitations.

The authors gratefully acknowledge the funding provided by the Natural Sciences Engineering Research Council of Canada (NSERC) through the second author's Discovery Grant program to undertake this work. The authors would also like to thank a group of final-year degree project students (Hon Cheng, Ian Dornbush, Jason Guenther and Landon Little) to conduct the vibration testing and developing the finite element modeling of the footbridge used in this paper. Special thanks are reserved to Hugh Briggs (the director of physical plant of Lakehead University) for his kind technical assistance to conduct the testing on the footbridge.

A.1. Alternating least square

ALS is mainly comprised of the following key steps to undertake simultaneous unfolding of three model matrices (Bro 1997):

• 1.  
  Keeping ${\boldsymbol{\theta }}$ and ${\boldsymbol{\phi }}$ same, ${\boldsymbol{\psi }}$ is solved using:

  $$\begin{eqnarray}&&\mathop{\min }\limits_{{\boldsymbol{\psi }}}\parallel \bar{{\bf{S}}}-[[{\boldsymbol{\theta }},{\boldsymbol{\phi }},{\boldsymbol{\psi }}]]{\parallel }^{2}\equiv \mathop{\min }\limits_{{\boldsymbol{\psi }}}\parallel {\bar{{\bf{S}}}}_{(3)}-{\boldsymbol{\psi }}{({\boldsymbol{\phi }}\odot {\boldsymbol{\theta }})}^{{\mathsf{T}}}{\parallel }^{2},\end{eqnarray} \tag{ 16 }$$

  where $\odot$ represents Khatri–Rao product. Given ${\boldsymbol{\theta }}\epsilon {{\mathfrak{R}}}^{I\times R}$ and ${\boldsymbol{\phi }}\epsilon {{\mathfrak{R}}}^{J\times R}$, then ${\boldsymbol{\theta }}\odot {\boldsymbol{\phi }}$ is a matrix with IJ rows and R columns and is expressed as:

  $$\begin{eqnarray*}\,\begin{array}{l}{\boldsymbol{\theta }}\odot {\boldsymbol{\phi }}=\left[\begin{array}{cccc}{\theta }_{11} & {\theta }_{12} & \cdots \cdots & {\theta }_{1R}\\ {\theta }_{21} & {\theta }_{22} & \cdots \cdots & {\theta }_{2R}\\ \vdots & \vdots & \cdots \cdots & \vdots \\ {\theta }_{I1} & {\theta }_{I2} & \cdots \cdots & {\theta }_{{IR}}\end{array}\right]\odot \left[\begin{array}{cccc}{\phi }_{11} & {\phi }_{12} & \cdots \cdots & {\phi }_{1R}\\ {\phi }_{21} & {\phi }_{22} & \cdots \cdots & {\phi }_{2R}\\ \vdots & \vdots & \vdots & \vdots \\ {\phi }_{J1} & {\phi }_{J2} & \cdots \cdots & {\phi }_{{JR}}\end{array}\right]\\ \quad \quad \,\,\,\,=\text{}\left[\begin{array}{cccc}{\theta }_{11}{\phi }_{:1} & {\theta }_{12}{\phi }_{:2} & \cdots \cdots & {\theta }_{1R}{\phi }_{:R}\\ {\theta }_{21}{\phi }_{:1} & {\theta }_{22}{\phi }_{:2} & \cdots \cdots & {\theta }_{2R}{\phi }_{:R}\\ \vdots & \vdots & \vdots & \vdots \\ {\theta }_{I1}{\phi }_{:1} & {\theta }_{I2}{\phi }_{:2} & \cdots \cdots & {\theta }_{{IR}}{\phi }_{:R}\end{array}\right],\end{array}\end{eqnarray*}$$

  where ${\phi }_{:k}$ represents kth column of ${\boldsymbol{\phi }}$.
• 2.  
  Optimal ${\boldsymbol{\psi }}$ is the least square solution which can be obtained using:

  $$\begin{eqnarray}&&{\boldsymbol{\psi }}={\bar{{\bf{S}}}}_{(3)}({\boldsymbol{\phi }}\odot {\boldsymbol{\theta }}){({{\boldsymbol{\phi }}}^{{\mathsf{T}}}{\boldsymbol{\phi }}* \text{}{{\boldsymbol{\theta }}}^{{\mathsf{T}}}{\boldsymbol{\theta }})}^{\dagger }.\end{eqnarray} \tag{ 17 }$$
• 3.  
  Then each component of ${\boldsymbol{\theta }}$, ${\boldsymbol{\phi }}$ and ${\boldsymbol{\psi }}$ are solved until the desired convergence is achieved:

  $$\begin{eqnarray}{\boldsymbol{\theta }} & \longleftarrow & {\bar{{\bf{S}}}}_{(1)}({\boldsymbol{\psi }}\odot {\boldsymbol{\phi }}){({{\boldsymbol{\psi }}}^{{\mathsf{T}}}{\boldsymbol{\psi }}\ast {{\boldsymbol{\phi }}}^{{\mathsf{T}}}{\boldsymbol{\phi }})}^{\dagger }\\ {\boldsymbol{\phi }} & \longleftarrow & {\bar{{\bf{S}}}}_{(2)}({\boldsymbol{\psi }}\odot {\boldsymbol{\theta }}){({{\boldsymbol{\psi }}}^{{\mathsf{T}}}{\boldsymbol{\psi }}\ast {{\boldsymbol{\theta }}}^{{\mathsf{T}}}{\boldsymbol{\theta }})}^{\dagger }\\ {\boldsymbol{\psi }} & \longleftarrow & {\bar{{\bf{S}}}}_{(3)}({\boldsymbol{\phi }}\odot {\boldsymbol{\theta }}){({{\boldsymbol{\phi }}}^{{\mathsf{T}}}{\boldsymbol{\phi }}\ast {{\boldsymbol{\theta }}}^{{\mathsf{T}}}{\boldsymbol{\theta }})}^{\dagger }.\end{eqnarray} \tag{ 18 }$$

Finally, the ALS estimates a tensor $\hat{\bar{{\bf{S}}}}={\sum }_{{\bf{r}}={\bf{1}}}^{{\bf{R}}}{{\boldsymbol{\theta }}}_{{\bf{r}}}\circ {{\boldsymbol{\phi }}}_{{\bf{r}}}\circ {{\boldsymbol{\psi }}}_{{\bf{r}}}$ such that the following objective function is minimized:

$$\begin{eqnarray}&&f({\boldsymbol{\theta }},{\boldsymbol{\phi }},{\boldsymbol{\psi }})=\parallel \bar{{\bf{S}}}-\hat{\bar{{\bf{S}}}}{\parallel }^{2}.\end{eqnarray} \tag{ 19 }$$

A unique decomposition can be obtained if the Kruskal condition (Kruskal 1977) is satisfied:

$$\begin{eqnarray}&&{l}_{{\boldsymbol{\theta }}}+{l}_{{\boldsymbol{\phi }}}+{l}_{{\boldsymbol{\psi }}}\geqslant 2R+2,\end{eqnarray} \tag{ 20 }$$

where ${l}_{{\boldsymbol{\theta }}}$, ${l}_{{\boldsymbol{\phi }}}$ and ${l}_{{\boldsymbol{\psi }}}$ are l-rank of the matrices ${\boldsymbol{\theta }}$, ${\boldsymbol{\phi }}$ and ${\boldsymbol{\psi }}$ respectively, where l-rank is defined as maximum number l such that each set of l columns of the matrix is linearly independent.

A.2. Details of the 5-DOF model

The 5-DOF model used in this paper is taken from (Lin et al 1994). Following are the mass, damping, and stiffness matrices of the 5-DOF model

$$\begin{eqnarray}M({{\rm{Ns}}}^{2}\,{{\rm{cm}}}^{-1})=\left[\begin{array}{ccccc}19.57 & 0 & 0 & 0 & 0\\ 0 & 19.57 & 0 & 0 & 0\\ 0 & 0 & 19.57 & 0 & 0\\ 0 & 0 & 0 & 19.57 & 0\\ 0 & 0 & 0 & 0 & 19.57\end{array}\right],\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&\ C({\rm{Ns}}\,{{\rm{cm}}}^{-1})\\ \ \ =\,\,\left[\begin{array}{ccccc}47.19 & -13.67 & -0.79 & 0.30 & 0.06\\ & 37.46 & -15.61 & -1.04 & 0.46\\ & & 36.22 & -16.46 & 0.11\\ & {\rm{sym}}. & & 34.26 & -14.28\\ & & & & 15.93\end{array}\right],\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&K({\rm{N}}\,{{\rm{cm}}}^{-1})\\ \ \ =\text{}\left[\begin{array}{ccccc}77\,108 & -36\,564 & 4549 & 1612 & -211\\ & 58\,596 & -35\,825 & 5481 & 1169\\ & & 58\,\,344 & -36587 & 7463\\ & {\rm{sym}}. & & 52\,688 & -22\,962\\ & & & & 14\,621\end{array}\right].\end{eqnarray} \tag{ 23 }$$