Modal identification of civil structures via covariance-driven stochastic subspace method

1.   Introduction

Modal analysis aims to identify the modal parameters of a dynamic system that include natural frequency, damping ratio and modal shape. It usually involves as an essential part for studying structural dynamics, and provides prerequisite information for a wide branch of researches, such as assessment of dynamic responses, damage identification, vibration control and optimal design of dynamic structures.

The pioneering works on modal identification can be traced back to the 1930s, with the majority conducted in the fields of mechanical engineering and aerospace engineering. Since the mid-1970s, investigations have been extended gradually to other fields. In recent years, due to the fast development of large-scale civil structures, modal analysis has been also received increasing concerns in civil engineering. Unlike the cases in mechanical/aerospace engineering where both the input (i.e., excitation) and output (response) information of the dynamic system is often available, those in civil engineering are usually characterized by an absence of input records. Consequently, the modal identification problems cannot be solved via the traditional transform-function-based methods which rely on both the input and output information. For example, it would be considerably costly and/or difficult for high-rise buildings and large-span bridges to be excited by artificial devices with a controllable pattern. For such cases, one may have to identify the modal parameters under a working condition, with the dynamic system subjected to ambient excitation (e.g., wind load or traffic load). This is usually termed as Operational Modal Analysis (OMA).

Many OMA techniques have been developed. Clarkson et al. proposed a method by using cross function for lightly damped structures subjected to random excitation ^[1]. Akaike established the so-called auto-regressive (AR) model for modal identification of white-noise excited systems ^[2]. In 1973, Ibrahim proposed a time-domain method by using structural response information. Later, Ibrahim and his cooperators modified the above method and developed the famous ITD method ^[3]. In 1973, Cole proposed the Random Decrement Technique (RDT) ^[4]. This method was further studied by Ibrahim et al. so that it becomes more meaningful from a mathematical viewpoint. In 1993, Bendat et al. used the Peak Picking (PP) method ^[5,6] to determine the natural frequencies of a dynamic system. This method, however, tends to lose efficiency for cases with dense modes. Later, Brincker put forward the Frequency Domain Decomposition (FDD) method ^[7] which may be regarded as a modified version of the PP method. This method decomposes the response spectrum matrix into a set of single-degree-of-freedom power spectra by singular value decomposition (SVD), and then identifies the modal parameters of the system via the PP method. In 1995, James and Carne proposed the Natural Excitation Technology (NExT) method by using correlation function of output signal, and applied this method to identify the natural frequencies and damping ratios of steam turbine blades at a working state ^[8,9]. In 1998, Huang documented an Empirical Modal Decomposition (EMD) method, which forms the key step of the so-called Hilbert-Huang-Transform (HHT) ^[10]. The HHT method was soon used to deal with modal identification problems by peers ^[11,12]. Besides the above deterministic modal identification methods, techniques have been also developed in a Bayesian framework to quantify the uncertainty of estimation results. The Bayesian modal identification method was first documented by Beck and Katafygiotis ^[13], and further developed and applied in many following studies ^[14,15].

This study utilizes the Stochastic Subspace Identification (SSI) method which was firstly presented by Overschee and Moor in 1991. The SSI method identifies modal parameters from the so-called state space matrix based on either input-and-output records or output-only records. It involves a number of robust numerical techniques, such as QR decomposition, SVD and least-square fitting, which are of great importance for reducing computational amount and for suppressing noise contaminated in both the datasets and the system model. The SSI method is currently acknowledged as one of the most advanced modal identification methods ^[16,17,18]. It overcomes some typical shortcomings associated with a frequency domain method, such as inadequacy for identification of closely-spaced modes and insufficient resolution in the frequency domain. More attractively, this method is able to identify modal frequencies, mode shapes and damping ratios of multiple modes of the system simultaneously, which makes it of pretty efficiency. Thus, it has been adopted widely for modal analysis since its debut ^{[19,20,21,22,23,24]}.

It should be noted that responses of a dynamic system subjected to ambient excitation are commonly characterized by small amplitude and high level of noise. Although great achievements have been made for modal analysis based on output-only records via the above mentioned OMA techniques, continuous efforts are still needed to improve the estimation performances, including accuracy, stability and efficiency, of these methods. In this study, main steps involved in the modal identification process via the covariance-driven SSI method are introduced first. A case study is then presented to demonstrate the accuracy and efficiency of this method, through comparing the corresponding results with those via an alternative method. The effects of noise contaminated in output signals on identification results are stressed. Special attention is also paid to how to determine the mode order accurately.

2.   State-space model of dynamic systems

2.1.   Continuous-time form

The dynamics of a multi-degree-of-freedom (MDOF) system can be expressed by:

in which, $M, {D_1}, K \in {\mathbb{R}^{{n_1} \times {n_1}}}$ respectively stand for the mass, damping and stiffness matrices; $U(t) \in {\mathbb{R}^{{n_{\rm{1}}} \times 1}}$ is the displacement vector at time instant t; $f(t){\rm{ = }}{H_1}u(t) \in {\mathbb{R}^{{n_1} \times 1}}$ denotes the external excitation which can be factorized in a matrix $u(t) \in {\mathbb{R}^{m \times 1}}$ and a matrix ${H_1} \in {\mathbb{R}^{{n_1} \times m}}$, with $u(t)$ describing the m inputs in time and ${H_1}$ describing the inputs in space, respectively.

For practical convenience, the above equation is conventionally converted into the state-space form ^[25,26]. Let's denote $x(t) = {\left[{\begin{array}{*{20}{c}} {U(t)} & {\dot U(t)} \end{array}} \right]^T}$ as the state-space vector of the dynamic system, Equation (1) can be then rewritten as:

where,

in which, ${A_c} \in {\mathbb{R}^{2{n_1} \times 2{n_1}}}$ is the state matrix, ${B_c} \in {\mathbb{R}^{2{n_1} \times m}}$ is the input matrix. Note that the dimension of the state-space vector equals to the number of independent variables required to describe the system.

On the other hand, as the output of dynamic systems is generally measured in terms of acceleration, velocity and displacement, the output vector can be expressed as:

where, $y(t) \in {\mathbb{R}^{l \times 1}}$ denotes the output, and ${C_a}, {C_v}, {C_d} \in {\mathbb{R}^{l \times {n_1}}}$ are the outputs matrices for acceleration, velocity and displacement, respectively.

Substituting (1) into (4) gives:

or:

with

where, ${C_c} \in {\mathbb{R}^{l \times {n_1}}}$ is the output matrix, and ${D_c} \in {\mathbb{R}^{l \times m}}$ is the direct transmission matrix.

Eqs 2 and 6 collectively refer to the state-space model of the continuous-time system. In reality, however, as the output signals of studied systems are usually provided in a discrete-time form, the state-space model of the continuous-time system needs to be converted into the corresponding discrete form.

2.2.   Discrete-time form

The integration form of Eq 2 can be readily derived as follows:

in which, t₀ denotes the initial time.

Let $t = (k + 1)\Delta t$, ${t_0} = k\Delta t$, $\tau '{\rm{ = }}t - \tau$, with k being a non-negative integer, and denote the state vector $x(k)$ at time $k\Delta t\;$ as $x(k) \equiv x(k\Delta t)$. By assuming that $u(\tau)$ remains as unchanged within each sampling interval $\Delta t$, Eq 8 turns to be:

Thus, Eqs 2 and 6 can be rewritten as the discrete form:

in which, $A{\rm{ = exp(}}{A_c}\Delta t{\rm{)}}$ is termed as the state matrix of the discrete-time system, and $B = (\int_0^{\Delta t} {} \exp ({A_c}\tau ')d\tau {\rm{'}})){B_c}$ is termed as the input matrix.

2.3.   Stochastic state-space model

Under real conditions, there are always external disturbances, modeling errors and measurement errors of the instrument. They are often referred to as noise. An overview of the linear time-invariant stochastic system is depicted in Figure 1. Denoting the noise caused by external environment as w_k and the one caused by measurement errors as v_k, and assuming that w_k and v_k are of mutually irrelative zero-mean white noise processes, one has:

where, E is a mathematical expectation operator, subscripts "p" and "q" represent two time instants, ${\delta _{pq}}$ is the Kronecker delta function whose value Equals to one for q = p and to zero for q≠ p, respectively.

To take into account the noise effects, the discrete-time form of the state-space model can be expressed as:

In the field of civil engineering, because it is usually difficult to quantify the external excitations through measurements, and more importantly, the external forces (e.g., wind load) commonly feature stochastic characteristics, for practical convenience, the input item u_k in Eqs 14 and 15 are conventionally regarded as a zero-mean white-noise process. Therefore, the last two items in the right-hand side of Eqs 14 and 15 can be respectively merged into a single zero-mean white-noise process:

where, for description convenience, the merged white-noise processes are expressed as the same form of the noise items w_k and v_k.

2.4.   Properties of stochastic state-space model

The SSI method is developed on the basis of the stochastic state-space model, in which the modal parameters of a dynamic system are determined via the state matrix A and the output matrix C. This section specifies the derivation process of the above two matrices by using output-only measurements.

Because both w_k and v_k belong to zero-mean white-noise processes, it is straightforward that$E[{w_k}\; {x_k}] = 0$, $E[{v_k}\; {x_k}] = 0$.

Denoting the auto-covariance matrix of the sate vector as $\Sigma {\rm{ = }}E\left[{{x_k}\; x_k^T} \right]$ and the covariance matrix between the state vector and output vector as $G{\rm{ = }}E\left[{{x_{k{\rm{ + 1}}}}\; y_k^T} \right]$, it can be deduced that:

Similarly, the auto-covariance matrix of the output vector, i.e., ${\Lambda _{\rm{0}}}$, is expressed as:

Meanwhile, it can be deduced from Eq 17 that:

Therefore,

Eqs 18–20 and 22 establish a close-form relationship between the output covariance and the state matrix A and output matrix C. Thus, the stochastic model can be identified.

3.   Modal identification via covariance-driven SSI method

This section introduces the basic theory of the covariance-driven SSI method. There are two basic steps involved in the covariance-driven SSI method: (1) to obtain the system matrix via SVD from the covariance matrix of output vector; (2) to identify the modal parameters of studied system from the system matrix.

Details of the methodology are shown in Figure 2. This method starts with constructing the so-called Hankel matrix from the output data. Through projection operations, the Hankel matrix can be used to obtain the Toeplitz matrix which contains sufficient information for determining the system matrix (i.e., A and C). Another attractive merit by using the Toeplitz matrix lies in that it can reduce the amount of analytical data significantly while maintaining the original information of the signal. Then the system matrix can be recognized via SVD manipulation. Finally, the modal parameters of the structure can be identified via eigenvalue decomposition.

3.1.   Hankel matrix

The Hankel matrix plays an important role in the SSI method ^[17,27]. It can be constructed from the output data directly. Define the output Hankel matrix as:

in which, the Hankel matrix Y_0/2i-1 is a matrix consisting of "2i" block rows and j columns, with each row composed of l sub-rows, l being the number of output channels (${y_k} \in {\mathbb{R}^{l \times 1}}$). Note that i and j are two key parameters involved in the SSI method which can influence the identification results markedly. This will be discussed in details in the case study section. Basically, to guarantee that the recognition algorithm is statistically meaningful, it is theoretically required that j should be infinitely large. In reality, however, the above requirement can be only roughly satisfied since the output signals are always finite in length. In Eq 23, Y_0/2i-1 is divided into two parts: "the past" and "the future", which are denoted as Y_p (past), Y_f (future) respectively. Both Y_p and Y_f consist of i block rows. The value of i can be determined somewhat artificially, given that it should be no less than the maximum computational mode-order of the system ^[12]. On the other hand, for computational efficiency, i should not be too large ^[28].

3.2.   Output covariance matrix

Define the output covariance as ${\Lambda _i} = E[{y_{k + i}}\; y_k^T]$. By assuming that the output data belongs to an ergodic process, ${\Lambda _i}$ can be estimated via the associated estimator ${\hat \Lambda _i}$ :

By using Eq 23, the Toeplitz matrix T_l/i can be readily obtained as follows:

Using the Toeplitz matrix of covariance sequence can reduce the computational amount during the modal analysis significantly, and thus improve the efficiency. It can be seen that after the above manipulation, the output data is changed from an original form of the Hankel matrix with a dimension of 2li × j to that of the Toeplitz matrix whose dimension is only li × li.

3.3.   Toeplitz matrix decomposition

Substituting Eq 22 into Eq 25 leads to:

which can be equally written as:

with:

where, ${\Gamma _i}$, ${\Delta _i}$ respectively represent the so-called extended observable matrix and the inversely extended controllable matrix.

Through SVD manipulation, the Toeplitz matrix can be decomposed as follows:

where, U and V are orthogonal matrix with ${U_1}, {V_1} \in {\mathbb{R}^{li \times 2{n_1}}}$, while S is a diagonal matrix${S_1} \in {\mathbb{R}^{2{n_1} \times 2{n_1}}}$, with the number of nonzero elements indexing the rank of the matrix, or the computational order of the studied dynamic system.

Reforming Eq 30 as follows:

in which, T is a non-singular matrix with an order equal to n. Particularly, one can use the unit matrix I instead of T.

Replacing T with I in the above equation, and comparing this equation with Eq 27, one has:

Based on Eqs 30 and 32–34, T_2/i+1 can be deduced:

It is clear that T_2/i+1 can be calculated directly based on the output records. Substituting Eqs 33–34 into the above equation results in:

From Eqs 28–29, the matrix C is a sub-matrix, or the first l rows, of ${\Gamma _i}$; while the matrix G is the last l columns of ${\Delta _i}$. Thus, both of the matrices A and C can be determined.

3.4.   Modal parameters

On obtaining the state matrix A and output matrix C, the modal parameters of the dynamic system can be identified according to the following analysis.

Through eigenvalue decomposition, A can be reformed as:

where, $\Lambda$ is a diagonal matrix which contains complex eigenvalues of the system, and Ф is a matrix composed of eigenvectors.

From Eq 10, the state matrix of discrete-time system A correlates with the one of the continuous-time system A_c via the following relationship:

According to Eqs 37–38, it is clear that:

in which, ${\Lambda _c}, {\Phi _c}$ are the counterparts of $\Lambda$ and Ф for the continuous-time system.

Based on the above equation, the relationship of eigenvalues between the continuous-time system and the discrete-time system can be deduced:

On the other hand, the eigenvalues of matrix A_c correlate with the natural frequency ${\omega _i}$ and damping ratio ${\xi _i}$ of the system via:

The values of f_i and ${\xi _i}$ can be then computed as follows:

where, a and b are the real part and imaginary part of the eigenvalue:

The mode shape of the system is:

So far, all the modal parameters of natural frequencies, damping ratios and mode shapes associated with multiple modes of the system can be determined.

3.5.   Determine system order via stability diagram method

A crucial step involved in the modal analysis via the SSI method lies in the determination of computational order of the system n. This parameter should be pre-assigned before identifying the modal parameters from the system matrices.

Some methods exist for determining the value of this parameter. Theoretically, it can be determined as the number of non-zero singular values of the Toeplitz matrix (refer to Eq 31). To take into account the noise effect in practice, it may be determined as the number of the first several descending singular values after which a distinct jump of the amplitude of singular values appears. Unfortunately, the above method tends to lose efficacy for the cases of ambient-excited civil structures, since the associated output records are usually contaminated with much more severe noise.

In reference to the above problem, Reynders et al. studied the uncertainty bounds of modal parameters identified via the SSI method, and proposed the stability diagram method ^[29]. The basic ideal of this technique is that the true modes correspond to more stable values of modal parameters than those of false modes which are usually caused by noise. Thus, one can estimate modal parameters via the SSI method for multiple times, with the computational order for each time varying at an acceptable interval within a certain range. Then, the identified values of a nominal modal parameter can be plotted against the computational order. Statistically, the false modes can only be identified occasionally, while the true modes should exist more stably. Therefore, a reasonable computational order can be recognized according to the above stability diagram where only the stable (or true) modes exist in the estimation results.

Figure 3 shows the flowchart of using the stability diagram method to recognize the mode order. One can compare the estimation results of each modal parameter for any two computational times that are associated with two adjacent computational mode orders. If the difference between the two kinds of results is less than a critical level, the identified modal parameters can be regarded as stable.

For the natural frequency f, the above critical level may be treated as the frequency resolution. For example, if a critical level of 0.01 Hz is adopted to judge the stability of identified natural frequency, and the difference of f between two adjacent computational mode orders is less than 0.01 Hz, then 0.01 Hz is the frequency resolution of identified f results.

Compared to f, the damping ratio $\xi$ becomes more sensitive to noise effects. Thus, the identified damping values for different computational times may vary significantly from one another. It is, therefore, less reasonable to set the critical level in an absolute form. Instead, a relative form can be adopted:

where, $\bar \xi$ is the average of identified damping ratios, and 0.05 is the critical level in the relative form.

For the mode shape, the similarity between two mode shape vectors${\rm{\{ }}\phi {\rm{\} }}$ and${\rm{\{ }}\varphi {\rm{\} }}$ can be quantified according to the Modal Assurance Criterion (MAC) ^[30]:

For MAC = 1.0, it means that the two mode shape vectors are identical, while for MAC = 0, it means that the two vectors are completely different. Conventionally, the mode shape can be considered as acceptably consistent for MAC > 0.8–0.9.

4.   Alternative modal identification method

In this study, an alternative modal identification method is also exploited for comparison purpose. This method is actually a combination of power spectral analysis and the Random Decrement Technique (RDT). Two specific methods are utilized to estimate the Power Spectral Density (PSD): the Welch method and the Yule-Walker method. One may resort to ^[31,32,33] for details about these methods. Based on the PSD estimation results, the natural frequencies can be identified by using the PP method, given that the modes are not closely spaced in the frequency domain. The single degree of freedom (SDOF) response component can be then separated from the output signal through filtering. The mode shape for each separated mode can be further determined via the cross-spectrum method based on SDOF responses at varied test points, or simply, by comparing the RMS (i.e., root mean square) values of the SDOF responses at different test points (refer to as RMS method, hereafter) ^[34].

The RDT method is a widely adopted time-domain technique to investigate the amplitude-dependence of natural frequency and damping ratio ^[35,36,37]. This method starts with the idea that the total response of a single degree of freedom (SDOF) dynamical system is superimposed by three components which are respectively associated with an initial displacement x₀, an initial velocity ${\dot x_0}$, and the external input

force at time instant t:

Compose a signature as:

where, sgn represents the sign operation, $E\{ x\left| C \right.\}$ is the conditional expectation under the triggering condition C, $\tau$ is the length of each sample involved in the signature, with$\hat x$ being the initial amplitude of the sample.

By assuming that the external load can be treated as a zero-mean white noise process, it is clear that both the velocity-related and force-related response components tend to vanish in the signature, leaving only the x₀-related component. By choosing different levels of displacement threshold $\hat x$, a series of signatures can be calculated. The damping ratio $\xi$ can be then determined by fitting the envelope of each signature:

where, the radian natural frequency ${\omega _0}$ can be estimated through PSD analysis from the output signal.

Note that although the RDT method is developed based on the displacement output records, this method can be practically operated on acceleration signals as well ^[38].

5.   Case study

A case study is presented in this section to demonstrate the model identification process via the SSI method. As shown in Figure 4a, the studied civil structure is a pedestrian bridge. The experiment was conducted under weak wind conditions, and consequently, the structure should be mainly subjected to vehicle-induced environmental excitation. Six tri-axial wireless accelerometers (Product model: JM3873) were utilized and fixed at different positions above the bridge deck. All the accelerometers had been calibrated and synchronized carefully before test. The experiment lasted for 1 hour. During the experiment, acceleration responses of the bridge at each testing point were first low-pass filtered and then sampled respectively at 32 Hz and 64 Hz. Figure 5 shows the time history of collected acceleration response at point No. 3 (Figure 4b).

The PSD of the acceleration response of the test bridge is first studied. The block length involved in the Welch method is set as 2¹³ = 8192 so that a frequency resolution no rougher than 0.01 Hz can be achieved. The order of the Yule-Walker method is taken as 500. Figure 6 shows the results for the vertical response from sensors No. 3 and No. 4 (Figure 4). As can be seen, the two kinds of results show good agreement, which indicates the effectiveness of both methods. There are four distinct peaks identified in the range below 10 Hz, i.e., 2.06 Hz, 3.12 Hz, 6.81 Hz, and 8.11 Hz. These peaks should correspond to the first four modes of the bridge.

Figure 7 shows the some typical results of the damping ratio estimated via the RDT method. It is clear that the random decrement signature can be well regarded as a free decay curve of a damped SDOF system. The estimated damping ratios at the two test points for the 2^nd mode under the concerned amplitude condition (i.e., ~2 mm/s²) are consistently equal to 0.52%.

Figure 8 shows the stability diagram obtained via the SSI method which is overlaid by the PSD estimation results. During the modal analysis process, three values of the block line number i involved in the Hankel Matrix are considered, i.e., 50, 75, and 100; while the computational mode order n is set to range from 2 to 100. The critical levels for natural frequency, damping ratio (in a relative form) and mode shape (i.e., MAC) are set as 0.01, 0.05 and 0.95, respectively.

Through comparison, it is seen that the modal identification results via the stability diagram method agree basically well with those obtained through PSD analysis, if the order of the system matrix involved in the SSI method is set in a range of 8-18, or equally the computational mode order is set in a range of 4-9. The results also demonstrate that the stability diagram method fails to provide accurate modal identification results when the values of i are too small. On the other hand, if i is set as a too large value, the obtained results would be contaminated with some false modes. It should be noted that the first mode component occupies only a considerably low proportion in the total response for this study. Due to the severe noise effects, the stable poles associated with this mode cannot be well identified via the stability diagram method.

Table 1 lists the modal parameters for the first four modes of the pedestrian bridge identified via both the covariance-driven SSI method and the alternative method, while Figure 9 compares the mode shapes determined via the above two methods. Basically, the two kinds of results agree well with each other, which demonstrates the effectiveness of both kinds of methods. Note that the SSI method can provide identification results of natural frequency, damping ratio and mode shape for multiple modes simultaneously. Thus, it is much more efficient than the alternative method which focuses either on one modal parameter or on a single mode.

6.   Conclusion

During the last decades, the SSI method has been developed as an advanced modal identification technique which is driven by output-only records. In this study, main steps involved in the modal identification process via the covariance-driven SSI method are introduced first. A case study on a real bridge is presented to illustrate the applicability of this method. The natural frequencies, damping ratios and mode shapes identified with the SSI method are found to agree well with the results obtained via an alternative method.

There are some overwhelming advantages of the SSI method compared to other traditional modal identification methods. Typically, the SSI method can identify not only the natural frequencies but also the modal shapes and damping ratios associated with multiple modes of the system simultaneously, making it of particular efficiency. Meanwhile, since the SSI method uses a series of powerful mathematical tools (SVD, QR decomposition, etc.), it possesses attractively good anti-noise ability. Moreover, the SSI method is able to work effectively for cases with closely spaced modes. These merits facilitate the SSI method to be exploited widely in dealing with OMA problems.

It is stressed that great attention should be needed for determining the computational mode order when one uses the SSI method for modal analysis. It is found in this study that when the energy contained in a SDOF response component is too limited, the SSI method may fail to work effectively due to the severe noise effects. Thus, continuous efforts are required to improve the working performance of the SSI method for such cases.

Acknowledgements

This study is supported by two grants from the National Natural Science Foundation of China (Project Nos: 51878194, 51578169).

Conflict of interest

All authors declare no conflict of interest in this paper.