Mechanical Systems and Signal Processing 25 (2011) 1028–1044 Contents lists available at ScienceDirect Mechanical Systems and Signal Processing journal homepage: www.elsevier.com/locate/jnlabr/ymssp Operational modal analysis by updating autoregressive model V.H. Vu a, M. Thomas a,n, A.A. Lakis b, L. Marcouiller c a Department of Mechanical Engineering, École de Technologie Supérieure, Montréal, Quebec, Canada H3C 1K3 Department of Mechanical Engineering, École Polytechnique, Montréal, Quebec, Canada H3C 3A7 c Hydro-Québec’s Research Institute, Varennes, Quebec, Canada J3X 1S1 b a r t i c l e i n f o abstract Article history: Received 16 July 2008 Received in revised form 25 August 2010 Accepted 27 August 2010 Available online 16 September 2010 This paper presents improvements of a multivariable autoregressive (AR) model for applications in operational modal analysis considering simultaneously the temporal response data of multi-channel measurements. The parameters are estimated by using the least squares method via the implementation of the QR factorization. A new noise ratebased factor called the Noise rate Order Factor (NOF) is introduced for use in the effective selection of model order and noise rate estimation. For the selection of structural modes, an orderwise criterion called the Order Modal Assurance Criterion (OMAC) is used, based on the correlation of mode shapes computed from two successive orders. Specifically, the algorithm is updated with respect to model order from a small value to produce a costeffective computation. Furthermore, the confidence intervals of each natural frequency, damping ratio and mode shapes are also computed and evaluated with respect to model order and noise rate. This method is thus very effective for identifying the modal parameters in case of ambient vibrations dealing with modern output-only modal analysis. Simulations and discussions on a steel plate structure are presented, and the experimental results show good agreement with the finite element analysis. & 2010 Elsevier Ltd. All rights reserved. Keywords: Modal identification Multivariate-autoregressive model Model orders selection QR factorization updating Parameters estimation Confidence intervals 1. Introduction Modal analysis identification techniques give useful information on modal parameters to understand the dynamic behavior of a structure [1]. However, in many cases which present nonlinear behaviors, the modal parameters have to be estimated under operating conditions and very often the excitations cannot be measured [2]. We must thus consider the operational modal analysis. Although the identification technique can be conducted both in the frequency [3] or time [4–6] domains, it is seen that the time domain is more suitable for operational modal analysis and can be classified into two groups. The first group lies in the fitting of response correlation functions, including the Ibrahim Time Domain (ITD) method [7], the Least Squares Complex Exponential (LSCE) method [8], the Covariance-driven Stochastic Subspace Identification (SSI-COV) method [9], and several modified versions of these methods for more suitable applications, particularly under harmonic excitations [10–13]. Other methods, based on parametric models, involve choosing a mathematical model to idealize the structural dynamic responses, including AutoRegressive Moving Average (ARMA) and AutoRegressive (AR) models [14–18]. For these autoregressive methods, a system identification algorithm is needed for estimating the model parameters. Among them, the Prediction Error Method (PEM) [19] is a common technique based on either the least squares estimate or on the Gauss–Newton iterative search. When the mode identification is necessary, multiples sensors have to be simultaneously recorded, and several applications of the multivariate AR model can be found using the ordinary least squares in the form of normal equations [20–22] or the Levinson algorithm [23]. The PEM iterative n Corresponding author. E-mail address: marc.thomas@etsmtl.ca (M. Thomas). 0888-3270/$ - see front matter & 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.ymssp.2010.08.014 V.H. Vu et al. / Mechanical Systems and Signal Processing 25 (2011) 1028–1044 d ^ D e(t) E^ fi I K K1,2 Kn L N p peff(com) Q Qa R Ra submatrices of R factorized matrices of R22 factor of the order updating t time index Ts sampling period R factor corresponding to added data T1,2,3 ui discrete complex eigenvalue y(t  i) the output vector with time delay i  Ts z(t) the regressor for the output vector y(t) zi damping ratio H real mode shapes matrix li continuous complex eigenvalue K model parameters matrix p Pi number P state matrix W complex mode shapes matrix (p) parameter at order p ^ estimated value T matrix/vector transpose conjugated transpose ¯ || absolute value Im(y) imaginary part Re(y) real part Trace(y) trace norm of a matrix Rij RðÞ 22 Nomenclature Ai 1029 matrix of parameters relating the output y(t i) to y(t) vector dimension or number of sensors estimated covariance matrix of the deterministic part the residual vector of all output channels estimated covariance matrix of the error part natural frequency unity matrix data matrix subdivided data matrices added data columns matrix complex eigenvectors matrix number of available data samples model order efficient (computing) model order orthogonal factor matrix of the QR factorization Q factor of the added columns factorization upper-diagonal factor matrix of the QR factorization R factor of the added columns factorization method, generally used to search for minimization, requires intensive computation and initial start-up values which are normally calculated using the least squares method. Furthermore, in some cases, the local minimization problem poses a big challenge [19]. In this paper, the multivariate autoregressive model is expressed in a convenient fashion for computation. The QR factorization gives an easy, fast and well-conditioned formulation for the least squares estimate of model parameters, and can be effectively updated with respect to the model order. A new factor based on the separation and evolution of the signal and noise is developed for the model order selection and noise rate estimation. The modal parameters are derived using the eigen-decomposition of the state matrix. An orderwise version of the correlation criterion called the Order Modal Assurance Criterion (OMAC) is defined for a user-friendly selection of modes. For interest on uncertainty in the parameter estimates, the confidence intervals for each modal parameter are computed. Finally, the method is applied both on simulated and experimental data of a steel plate in comparison to finite element method. 2. Vector-autoregressive model and its order updating In operational modal analysis, we assume that the excitation is unknown. Since the modal analysis is conducted by using several d channels of measurements synchronized for data acquisition at sampling period Ts, a multivariate autoregressive model of pth order and dimension d can be utilized to fit the measured data [14,19]. yðtÞ ¼ KzðtÞ þeðtÞ where K ¼ ½ A1 A2 ð1Þ  Ap  of size d  dp is the parameters matrix, Ai of size d  d is the matrix of parameters relating the output y(t  i) to y(t), z(t) of size dp  1 is the regressor for the output vector y(t), zðtÞT ¼ ½ yðt1ÞT yðt2ÞT  yðtpÞT , y(t  i) of size d  1 (i= 1:p) is the output vector with delayed time i  Ts, e(t) of size d  1 is the residual vector of all output channels considered as the error of model. If N consecutive output vectors of the responses from y(t) to y(t +N  1) are taken into account, the model parameters can be obviously estimated with the least squares method. The following section reports the solution of this least squares problem by using of the well-known QR factorization [24–26]. It is summarized for better understanding as follows: The N  dp+ d data matrix is first constructed from available data: 3 2 zðtÞT yðtÞT 7 6 6 zðt þ1ÞT yðt þ1ÞT 7 7 K¼6 ð2Þ 7 6 ... ... 5 4 T T yðt þN1Þ zðt þ N1Þ 1030 V.H. Vu et al. / Mechanical Systems and Signal Processing 25 (2011) 1028–1044 The QR factorization of the data matrix K =Q  R gives an orthogonal matrix Q (size N  N) and R (size N  dp+ d) an upper triangular matrix with the form 2 3 R11 R12 6 7 R22 5 ð3Þ R¼4 0 0 0 where R11, R12 and R22 are of size dp  dp, dp  d and d  d, respectively. ^ and of The parameters matrix is finally derived, as well as the estimated covariance matrices of the deterministic part D the error E^ (both of size d  d): T K ¼ ðRT12 R11 ÞUðRT11 R11 Þ1 ¼ ðR1 11 R12 Þ ð4Þ ^ ¼ 1 RT R12 D N 12 ð5Þ 1 E^ ¼ RT22 R22 N ð6Þ The solution thus yields to the computation of R submatrices. The computation of those submatrices should be updated with respect to model order to avoid the problem of other algorithms which need a prior evaluation of model order and require a repetitive computation for a set of model orders. Although algorithms for the updating of QR factorization had been developed in [25] by using the Givens rotations, these transformations need to be repetitively performed on half of matrix elements which can be large for multi-sensors modeling. In this paper, we present a computation for the updating of the model when the order is increasing and it is shown that if only submatrices of R are required, the problem results in applying the QR factorization or Givens rotations to only a submatrix of the previous solution. Suppose that we have the data matrix K(p) at the order p and its factored matrices Q(p), R(p) computed from Eqs. (2) and (3): 3 2 yðtÞT zðtÞT 7 6 i 6 zðt þ 1ÞT yðt þ 1ÞT 7 h ðpÞ 7 ¼ K1 K2 ð7Þ KðpÞ ¼ 6 7 6 ... ... 5 4 T T yðt þ N1Þ zðt þ N1Þ R ðpÞ RðpÞ 11 6 ¼4 0 0 2 RðpÞ 12 RðpÞ 22 0 3 7 5 At the order p + 1, the data matrix can be subdivided as follows [25]: h i ðpÞ K K2 Kðp þ 1Þ ¼ K1 ð8Þ n where K of size N  d comprises the added d columns: 3 2 yðkðp þ 1ÞÞT 6 T 7 6 yðk þ 1ðp þ1ÞÞ 7 7 K ¼ 6 7 6 ... 5 4 yðk þN1ðp þ 1ÞÞT ð9Þ Then we can compute the following matrix: h ðpÞT ðpÞ Q ðpÞT Kðp þ 1Þ ¼ Q K1 Q ðpÞT K Q ðpÞT K2 i 2 RðpÞ 6 11 6 ¼4 0 0 3 T1 RðpÞ 12 T2 7 7 RðpÞ 22 5 T3 ð10Þ 0 where T1 of size dp  d, T2 of size d  d and T3 of size N  dp  d  d are extracted from 2 3 T1 6 7 ðpÞT  Q K ¼ 4 T2 5 T3 " # T2 If we apply the QR factorization on the submatrix , it yields to T 3 " #   T2 Ra ¼ Qa T3 0 ð11Þ V.H. Vu et al. / Mechanical Systems and Signal Processing 25 (2011) 1028–1044 1031 where Ra of size d  d is an upper diagonal matrix and Qa of size N dp  N  dp is the product of the Householder transformations or Givens rotations. Then Eq. (10) becomes 2 ðpÞ 3 " # R11 T1 RðpÞ 12 I 0 6  7 Ra R22 7 Q ðpÞT Kðp þ 1Þ ¼ 6 ð12Þ 4 0 5 0 Q Ta R 0 0 22 where R22 of size d  d and R 22 of size N dp  d  d are obtained from multiplication "  #   R22 R22 ¼ Q Ta R 0 22 It can be noticed that the submatrix R 22 in the right hand side of Eq. (12) is not an upper diagonal matrix and must be triangularized by a small orthogonal transformation to yield the QR decomposition of the data matrix K(p + 1). This ðp þ 1Þ ðp þ 1Þ ðp þ 1Þ modification evidently does not affect other components of Eq. (12). Hence the submatrices R11 , R12 and R22 of model at order p+ 1 are exactly updated: " ðpÞ # " ðpÞ # R12 R11 T1 ðp þ 1Þ ðp þ 1Þ ðp þ 1Þ ¼ ¼ ¼ R ð13Þ R11 ; R12 ; R22 22 R22 0 Ra ^ ðp þ 1Þ are updated, as shown in Eqs. (5) and (6): Next, the model parameters matrix K(p + 1) and covariance matrix D h iT ðp þ 1Þ 1 ðp þ 1Þ Kðp þ 1Þ ¼ ½R11  R12 ð14Þ ^ ðpÞ þ R T R ðÞ ^ ðp þ 1Þ ¼ Rðp þ 1ÞT Rðp þ 1Þ ¼ RðpÞT RðpÞ þ R T R ðÞ ¼ D D 22 22 22 22 12 12 12 12 ð15Þ ðp þ 1Þ is also updated: Since all transformations are orthogonal, covariance matrix E^ ðp þ 1Þ T ðp þ 1Þ ðpÞ T  T  ^ ðpÞ E^ ðp þ 1Þ ¼ R22 R 22 ¼ RðpÞ 22 R22 R 22 R 22 ¼ E R22 R 22 ð16Þ Once the model parameters are updated, the calculation of modal parameters is straightforward [14]. The state matrix of the system is first established from autoregressive parameters: 3 2 A1 A2 . . . Ap1 Ap 6 I 0 ... 0 0 7 7 6 7 6 6 I ... 0 0 7 P¼6 0 ð17Þ 7 7 6 ... ... ... ... 5 4 ... 0 0 ... I 0 The eigen-decomposition of the state matrix can be obtained: 2 3 u1 0 0 0 6 0 u 0 0 7 6 7 1 2 7L P ¼ L6 6 0 7 0 & ^ 4 5 0 0 . . . udp ð18Þ The continuous eigenvalues, natural frequencies, damping rates and complex modes of the system can be computed: Eigenvalues : Frequencies : li ¼ lnðui Þ Ts fi ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Re2 ðli Þ þ Im2 ðli Þ Damping rates : Complex modes : 2p xi ¼  Reðli Þ 2pfi h W ¼ W1 W2 ... i Wdp ¼ I 0 . . . 0 L   ð19Þ In modal analysis, the selection of modes may be realized by using either a modal signal-to-noise classification MSN [14] or the MAC (Modal Assurance Criterion [14,27]) which defines the correlation between the identified mode shape vector with a reference vector (extracted from numerical or experimental results). When conducting an operational modal analysis under a noisy environment, the MSN index is still reliable but is difficult to apply because it needs a rigorous attention on the number of modes selected from the spurious ones and an apriori knowledge about frequency and damping [14] which in fact is not available in output-only modal analysis. Since our algorithm is updated with respect to 1032 V.H. Vu et al. / Mechanical Systems and Signal Processing 25 (2011) 1028–1044 model order with availability of stabilized diagrams, we propose a new version of correlation criterion, which we have called OMAC (Order Modal Assurance Criterion). In this index, the MAC [14] of the ith mode is replaced by the correlation of the identified mode shapes given by two successive model orders p and p  1, formulated as follows:  h iT    WðpÞ Wðp1Þ    i i ðpÞ ffi ffirffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð20Þ OMACi ¼ rhffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi iT h iT WiðpÞ WiðpÞ Wiðp1Þ Wiðp1Þ where WðpÞ and Wiðp1Þ are the ith identified complex mode shape vectors at orders p and p  1, respectively, and i ðpÞ W i signifies the conjugated transpose of WðpÞ i . This OMAC parameter can be thus considered as a complementary information of the usual MAC (which compares the computed mode with a reference), by analysing the stability of the mode with the order of computation, since it allows to eliminate perturbation modes, before to compute the MAC. It is evident that the measurement with unobserved perturbations produces an uncertainty in the parameter estimation and hence on the modal parameters. Therefore, an analysis of confidence intervals of modal parameters should be taken into account. Uncertainty in structural dynamics has been generally introduced and bibliographically reviewed in [31] where all experimental modal analysis contributed on non-parametric methods. A recent derivation on variances of modal parameters has been presented in [32] which supposed an availability of the transfer function in the frequency domain. It is seen that the uncertainty of autoregressive model parameters can be transferred from the estimation of the least square estimate by differentiation [33]. In this paper, the derivatives of modal parameters are derived from Eq. (19) and their confidence intervals are straightforward transferred as described in [33]: Derivative of frequencies : Derivative of damping rates : Reðli ÞReðl_ i Þ þImðli ÞImðl_ i Þ f_ i ¼ 4p2 fi Reðl_ i Þ f_ i þ z_ i ¼ zi Reðli Þ fi ! Fig. 1. NOF evolution and efficient order selection. Fig. 2. Configuration of testing plate. V.H. Vu et al. / Mechanical Systems and Signal Processing 25 (2011) 1028–1044 Fig. 3. NOF on simulated data. Fig. 4. Frequency stabilization diagram on simulated data. 1033 1034 V.H. Vu et al. / Mechanical Systems and Signal Processing 25 (2011) 1028–1044 Derivative of Hl,i, the lth component of the ith real mode shape: ðHl,i Þ2 ¼ ðWl,i Þ2 ¼ Re2 ðLl,i Þ þIm2 ðLl,i Þ, l ¼ 1,2,. . .,d _ ¼ 9ReðL Þ9ReðL_ Þ þ9ImðL Þ9ImðL_ Þ 9Wl,i 9H l,i l,i l,i l,i l,i ð21Þ Several observations should now be pronounced:  With only a small factorization of R submatrices, the model parameters and covariance matrices can be exactly updated to a higher order. This technique is much preferred than the repetitive QR factorization for each order value.  From the update of covariance matrices, it is noticed that, as the model order increases, the estimated signal part D^   increases and the estimated noise part E^ decreases monotonically. The changing amount is significant at low orders and negligible at high orders. This feature is thus an idea for a new orders selection criterion as discussed later. The introduction of the OMAC is very convenient with the order updating algorithm since the modes shapes are successively constructed. It becomes easy to check if a stabilized frequency comes from the structural properties when its corresponding OMAC is closed to unity within all considered model orders. The updating of the solution with respect to model order allows an evaluation of uncertainty of the modal parameters with respect to model order and to the noise rate, as seen in the experiment section. 3. A new index for model order selection and noise estimation The selection of the optimal model orders is the first step in the model based identification process. Among the better ^ and a penalty function, such as known methods are the criteria based on the statistical properties of prediction errors eðtÞ the Final Prediction Error (FPE) and the Minimum Description Length (MDL) [34]. However some remarks should be tackled:  These criteria are primarily based on the evolution of the error covariance which monotonically decreases with respect to the model order. Table 1 Results on simulations. Mode 1 2 3 4 5 6 7 Simulated frequencies (Hz) Simulated damping rate (%) Identified frequencies (Hz) Impact force, order 5 Identified damping rate (%) Impact force, order 5 Identified frequencies (Hz) Random force, order 5 Identified damping rate (%) Random force, order 5 40.8 1 40.78 77.5 1 77.5 112.9 1 112.8 172.3 1 171.6 222.9 1 221.5 290.4 1 287.4 294.4 1 291.3 1 1 1 1 1.1 1.1 1.1 39.7 75.7 112.3 169.7 218.8 278.1 286.5 1.9 0.7 0.6 0.6 0.6 Fig. 5. Real testing data. 1 0.2 V.H. Vu et al. / Mechanical Systems and Signal Processing 25 (2011) 1028–1044 1035  FPE and its variants asymptotically choose the correct order model if the underlying multiple time series has high dimensions d [35], but tends to overestimate the model order as the data length increases [36].  MDL and its variants outperform with long recorded data and are strongly consistent when the data length tends to infinite [28]. However, the application of these criteria first requires the selection of a possible interval for the model order to be used, and then an evaluation of model parameters. The selection of the order is made on the basis of Fig. 6. NOF evolution at different noise rates. Fig. 7. Comparison NOF and MDL. 1036 V.H. Vu et al. / Mechanical Systems and Signal Processing 25 (2011) 1028–1044 minimum variances. It is necessary, therefore, to have a prior evaluation with different orders, which results in a significant computation time. Other authors [29,30] have proposed new approaches derived from the MDL criterion, and based on the minimum eigenvalue of the covariance matrix without prior evaluation. To determine the model order, it is necessary to determine a limiting upper value in order to establish the order of the matrix. However, it is difficult to experimentally determine that value when significant noise is present. In modal analysis, we must consider not only the error sequence, but also the deterministic part of the signal, since this latter contains the modal information of the system. Based on these facts and to compare to the well-known optimal model criterion MDL, this paper proposes an innovative factor for the selection of an efficient model order peff based on an analysis of the noise-to-signal ratio (NSR). This order is defined as the smallest order value which can be used to fit the data with a negligible discrepancy and hence can be effectively used for modal analysis. While the MDL criterion considers only the prediction error, the estimated noise-to-signal ratio (NSR) can be defined ^ as previously derived in from the trace norm part of the estimated deterministic and error covariance matrices E^ and D, Eqs. (15) and (16): ^ ¼ NSR ^ TraceðEÞ ^ TraceðDÞ ð%Þ or ^ ¼ 10 log NSR 10 ^ TraceðEÞ ^ TraceðDÞ ðdBÞ ð22Þ Then a Noise-ratio Order Factor (NOF) is defined as being the variation of the NSR between two successive orders: NOFðpÞ ¼ NSRðpÞ NSRðp þ 1Þ ð23Þ It is seen that the NSR decreases monotonically with respect to model order, and contains properties of both stochastic (on numerator) and deterministic (on denominator) norms. The NOF is therefore always positive and is a representative factor for the convergence of the NSR, which changes significantly at low orders and converges at high orders. Since this factor is positive and close to zero, the convergent to zero property can be used to set the selection of efficient model order which is thus easier to do on the curve evolution (Fig. 1). 4. Application on a plate structure and discussions The method presented above is applied to a clamped-free–clamped-free rectangular steel plate [37] with dimensions of 500  200  1.9 mm3 (Fig. 2). Material properties are as follows: elastic modulus E= 200 GPa, the Poisson coefficient n = 0.29 Table 2 Selection of orders and noise rate estimation. Simulated noise rate (NSR) Estimated noise rate (%) Efficient order peff from NOF Optimal order from MDL Computational order pcom 0% 0.03 6 21 9 1% 1.20 6 11 12 10% 10.52 6 6 12 Fig. 8. Frequency stabilization diagram at NSR= 100%. 100% 102.1 6 6 12 200% 201.2 6 6 14 400% 386.5 6 4 14 V.H. Vu et al. / Mechanical Systems and Signal Processing 25 (2011) 1028–1044 1037 and density r = 7872 kg/m3. Six accelerometers are mounted on the plates to simultaneously record the responses at the measurement locations. 4.1. Numerical simulations The plate has been numerically excited by assumed impact and random forces. The responses have been computed for the first seven modes till 350 Hz by the modal superposition method and by considering a 1% damping ratios for all modes. Fig. 9. Modal parameters identification with confidence intervals (NSR= 100%). 1038 V.H. Vu et al. / Mechanical Systems and Signal Processing 25 (2011) 1028–1044 Fig. 9. (Continued) A sampling frequency of 1280 Hz has been applied and no noise has been added on the responses. It can be noticed that the sixth and seventh frequencies are close. Fig. 3 shows the NOF factor under the two excitations and it is seen that in both cases, the efficient model order can be clearly identified as 5 under an impact force (Fig. 3a), while it is less clear under a random force (Fig. 3b). V.H. Vu et al. / Mechanical Systems and Signal Processing 25 (2011) 1028–1044 1039 Fig. 9. (Continued) The stabilization diagrams in Fig. 4 show the first seven natural frequencies, with good accuracy for both excitations even if some frequencies are close (6th and 7th). However, we must notice that a random excitation produces a higher variance for the identification, particularly at these two close frequencies (Fig. 4b and Table 1). The identification of damping is good for both excitations except at the 7th mode, where the result is erroneous due to the proximity of the 6th frequency. Table 1 compares the first seven modal parameters of the plate as simulated by Ansys and those identified at efficient order 5, when no noise is added. A good agreement is found between simulated and identified modal parameters. However, a larger variance of the identified damping ratios is observed when the excitation is random. 4.2. Real structure testing The same plate has been experimentally excited with an uncontrolled shock, and the excitation forces (amplitudes and locations) have not been considered. The responses at the six locations have been simultaneously acquired (Fig. 5) and sampled at a frequency of 1280 Hz. We can see from the figure that the transient response has been perturbed by a second excitation. Despite this unexpected perturbation, the entire signal of 4000 samples has been considered. Furthermore, several additive random white noises (1–400% of the signal r.m.s.) have been numerically added to each measurement channel for the assessment of noise effect. 4.2.1. Model order selection and noise rate estimation In order to select an efficient model order for the system at different noise rates, the method described earlier has been used. Fig. 6 presents the evolution of the NOF factor at different noise rates. It can be observed that the efficient order of the system can be set close to 6 whatever the noise rate. 1040 V.H. Vu et al. / Mechanical Systems and Signal Processing 25 (2011) 1028–1044 To check the efficiency of the NOF method, the results have been compared to the MDL criterion [28,34]. Fig. 7 exhibits the comparison of NOF and MDL at very low noise rate (Fig. 7a, NSR = 1%) and very high noise rate (Fig. 7b, NSR =400%). It can be noticed that MDL gives the order between 4 and 11, respectively, with the noise rate while the NOF is found at order 6, regardless the noise rate. In fact, there is a good agreement between NOF and MDL under a moderate noise environment, but the MDL may overestimate the order with noisy-free data and underestimate the order at very high noise levels. Since we are looking for a criterion which is stable regardless the noise level, the proposed NOF factor revealed thus better than MDL for order selection in a wide range of noise environment. Table 2 shows the comparison of model orders, computed from MDL and NOF methods and the estimation of noise at different additive noises levels, as computed from Eq. (23). The computational order pcom shown in Table 2 may be selected greater or equal to the efficient order peff. The pcom has been selected from a threshold of the uncertainty on modal parameters as it is described later. As expected, it is shown that the computational order increases with the noise level. 4.2.2. Selection of modes and modal parameters identification Consider the case with a noise rate of 100%. For an assessment of model orders, a frequency stabilization diagram is constructed from orders 2 to 30 (Fig. 8). It is seen in agreement with NOF that a model order of 6 is the smallest value Fig. 10. OMAC diagram (NSR= 100%). Table 3 Identified modal parameters with their 95% confidence intervals. Mode Mode 1 Mode 2 Mode 3 Mode 4 NSR= 0% (pcom = 9) Frequency/confidence Damping/confidence 38.6/0.04 1.29/0.11 74.6/0.03 0.27/0.05 107.7/0.14 0.40/0.13 162.9/0.78 2.19/0.47 NSR= 1% (pcom = 12) Frequency/confidence Damping/confidence 38.6/0.09 1.47/0.24 74.6/0.2 0.27/0.03 107.8/0.14 0.4/0.13 NSR= 10% (pcom =12) Frequency/confidence Damping/confidence 38.6/0.11 1.54/0.3 74.6/0.03 0.27/0.04 NSR= 100%(pcom = 12) Frequency/confidence Damping/confidence 38.6/0.12 1.44/0.3 NSR= 200% (pcom = 14) Frequency/confidence Damping/confidence NSR= 400% (pcom = 14) Frequency/confidence Damping/confidence Mode 5 Mode 6 Mode 7 208.2/0.02 0.36/0.01 273.0/0.06 0.26/0.02 288.2/0.01 0.18/0.005 163.9/0.51 1.24/0.31 208.2/0.018 0.37/0.01 273.0/0.06 0.26/0.02 288.2/0.13 0.21/0.05 107.8/0.14 0.4/0.13 163.9/0.52 1.26/0.32 208.2/0.018 0.37/0.01 273.0/0.07 0.27/0.03 288.2/0.18 0.23/0.06 74.6/0.02 0.27/0.03 107.8/0.1 0.36/0.09 163.9/0.4 0.89/0.25 208.2/0.02 0.37/0.01 273.0/0.07 0.27/0.03 288.2/0.19 0.24/0.06 38.6/0.12 1.44/0.31 74.6/0.02 0.27/0.03 107.8/0.1 0.36/0.09 163.9/0.4 0.89/0.25 208.2/0.02 0.37/0.01 273.0/0.05 0.27/0.02 288.2/0.13 0.21/0.04 38.6/0.12 1.44/0.31 74.6/0.02 0.27/0.03 107.8/0.1 0.36/0.1 163.9/0.4 0.89/0.25 208.2/0.02 0.37/0.01 273.0/0.05 0.27/0.02 288.2/0.13 0.21/0.04 V.H. Vu et al. / Mechanical Systems and Signal Processing 25 (2011) 1028–1044 1041 required for the stabilization of all natural frequencies. Seven natural frequencies of Table 1 can distinctly be identified on the stabilization diagram. However, one can wonder if the frequency at 120 Hz is a real frequency or not. For the effective selection of computational order and structural modes, stabilized diagrams of the corresponding prospective frequencies and damping rates with their 95% confidence interval are constructed in Fig. 9. The OMAC for those frequency candidates are plotted in Fig. 10. Several discussions may be made:  On the frequency stabilization diagram, a natural frequency must be stable from the minimum order, since the identified part of the signal does not change significantly above this order value.  A stable frequency belongs to a mode if its damping rate takes a non-zero positive stable value on the diagram and if  the OMAC is closed to unity above the efficient order. Otherwise, it must belong to a harmonic excitation (zero damping) or to a computational frequency. The destabilizations of OMAC, frequency and damping reveal that the doubtful frequency at 120 Hz (Fig. 9d) was not a natural frequency (and can be due to electrical perturbation or computational mode). Even by considering an NSR of 100%, the results show that the confidence intervals converge to the good value at high orders, and give us the uncertainty of the estimation (such as at the 5th frequency (Fig. 9e) that presents a relatively large uncertainty). By considering an acceptable uncertainty (threshold), a computational order pcom may thus be selected for the modal parameter identification. Similar conclusions were found for other simulated noise rate cases (Table 1). Fig. 11. Mode shape identification at NSR= 100% and by FEA. 1042 V.H. Vu et al. / Mechanical Systems and Signal Processing 25 (2011) 1028–1044 Fig. 11. (Continued) Table 3 shows the results of modal parameter identification (frequency in Hz and damping rate in %) of the seven structural modes selected at the corresponding computational order with their 95% confidence interval, for different noise rates ranging from 0% to 400%. It is found from Table 3 that:  Whatever the noise rate, the modal parameters are accurately identified at the computational order.  When the noise rate changes, the modal parameters and their confidence intervals are considered to be stable which proves the performance of the proposed method even in a noisy environment.  The confidence intervals of frequencies are relatively much smaller than those of damping ratios. This explains the reliability on identification of frequency compared to the damping. Fig. 11 plots the seven scaled mode shapes identified with a noise rate NSR = 100% versus the mode shapes numerically computed by finite element method and the corresponding MAC factors (after having removed the perturbation mode at V.H. Vu et al. / Mechanical Systems and Signal Processing 25 (2011) 1028–1044 1043 120 Hz). In the left figures, the mode shapes are represented and interpolated from the deformation of the six sensor locations (on grid lines) while the confidence intervals are displayed by the darkness. The MAC is computed for each mode in order to compare with numerical results. It is found that even if the signal is perturbed by a very high noise rate and some modes are close, each identified mode shape is corresponding to an analytical one. As expected, the seventh and higher mode shapes are not well identified because of the limitation of the sensor number. 5. Conclusions An operational modal analysis based on a multivariable autoregressive model is presented with the ability of computing the uncertainty on modal parameters identification of selected modes, after updating and selecting efficient model orders. The least squares implemented in the form of QR factorization of the data matrix produces a fast, conditioned and convenient algorithm for updating. The Noise rate Order Factor (NOF) is effective for selecting the efficient model order from which the modal parameters converge regardless the noise rate. It is seen that, compared to the MDL method, the NOF is more stable with respect to noise and hence the modal parameters start to be stable from this order in the stabilization diagrams. The new orderwise version of the correlation factor called OMAC derived between two successive model orders is useful for the selection of structural modes. The confidence intervals of each natural frequency and damping ratio are added to the modal parameters on the stabilization diagram, which decrease as the model order increases. Computational model order can be chosen from acceptable confidence intervals for automatic identification of modal parameters, and hence varies with the noise rate. Simulations and experiments on a steel plate show that this new method is a good technique for operational modal analysis even in noisy environments ranging from 0% to 400%. Further studies are conducted to develop the algorithm for online identification in the time domain and for nonstationary systems. Acknowledgements The support of NSERC (Natural Sciences and Engineering Research Council of Canada), through Research Cooperative grants was gratefully acknowledged. The authors would like to thank Hydro-Quebec Research Institute for their collaboration. References [1] N.M.M. Maia, J.M.M. Silva, Modal analysis identification techniques, Royal Society, No. 359-2001, 2001, pp. 29–40. [2] M. Thomas, K. Abassi, A.A. Lakis, L. Marcouiller, Operational modal analysis of a structure subjected to a turbulent flow, in: Proceedings of the 23rd Seminar on Machinery Vibration, Canadian Machinery Vibration Association, Edmonton, AB, October 2005, 10 pp. [3] N.-J. Jacobsen, P. Andersen, R. Brincker, Using enhanced frequency domain decomposition as a robust technique to handle deterministic excitation in operational modal analysis, in: Proceedings of the International Operational Modal Analysis Conference (IOMAC 2007), Copenhagen, Denmark, vol. 4, 2007, pp. 193–200. [4] L. Hermans, H. Van der Auweraer, Modal testing and analysis of structures under operational conditions: industrial applications, Mechanical Systems and Signal Processing 13 (2) (1999) 193–216. [5] V.H. Vu, M. Thomas, A.A. Lakis, Operational modal analysis in time domain, in: Proceedings of the 24th Seminar on Machinery Vibration, Canadian Machinery Vibration Association, ISBN:2-921145-61-8, Montréal, 2006, pp. 330–343. [6] V.H. Vu, M. Thomas, A.A. Lakis, L. Marcouiller, A time domain method for modal identification of vibratory signal, in: Proceedings of the First International Conference on Industrial Risk Engineering CIRI, Montreal, ISBN:978-2-921145-65-7, December 2007, pp. 202–218. [7] S.R. Ibrahim, E.C. Mikulcik, Method for the direct identification of vibration parameters from the free responses, Shock and Vibration Bulletin (47) (1977) 183–198. [8] D.L. Brown, R.J. Allemang, R.D. Zimmerman, M. Mergeay, Parameter estimation techniques for modal analysis, SAE Paper No. 790221, SAE Transactions, vol. 88, 1979, pp. 828–846. [9] B. Peeters, System identification and damage detection in civil engineering, Ph.D. Thesis, K.U Leuven, Belgium, 2000, 256 pp. [10] P. Mohanty, D.J. Rixen, Operational modal analysis in the presence of harmonic excitation, Journal of Sound and Vibration 270 (1–2) (2004) 93–109 (6 February). [11] P. Mohanty, D.J. Rixen, Modified SSTD method to account for harmonic excitations during operational modal analysis, Mechanism and Machine Theory 39 (12) (2004) 1247–1255 (December). [12] P. Mohanty, D.J. Rixen, A modified Ibrahim time domain algorithm for operational modal analysis including harmonic excitation, Journal of Sound and Vibration 275 (1–2) (2004) 375–390 (6 August). [13] M. Gagnon, S.-A. Tahan, A. Coutu, M. Thomas, Operational modal analysis with harmonic excitations: application to a wind turbine, in: Proceedings of the 24th Seminar on Machinery Vibration, Canadian Machinery Vibration Association, ISBN:2-921145-61-8, Montréal, 2006, pp. 320–329. [14] S.M. Pandit, Modal and Spectrum Analysis: Data Dependent Systems in State Space, J. Wiley and Sons, New York, NY, 1991, 415 pp. [15] F. Gonthier, M. Smail, P.E. Gauthier, A time domain method for identification of dynamic parameters of structures, Mechanical Systems and Signal Processing 7 (1) (1993) 45–56. [16] P. Andersen, Identification of civil engineering structures using vector ARMA models, Ph.D. Thesis, Aalborg University, 1997, 244 pp. [17] M. Smail, M. Thomas, A.A. Lakis, ARMA models for modal analysis: effect of model order and sampling frequency, Mechanical Systems and Signal Processing 13 (6) (1999) 925–941. [18] V.H. Vu, M. Thomas, A.A. Lakis, L. Marcouiller, Identification of modal parameters by experimental operational modal analysis for the assessment of bridge rehabilitation, in: Proceedings of the International Operational Modal Analysis Conference (IOMAC 2007), Copenhagen, Denmark, April 2007, pp. 133–142. [19] L. Ljung, System Identification—Theory for the User, Prentice-Hall, Upper Saddle River, NJ, 1999, 609 pp. [20] U. Kadakal, O. Yuzugullu, A comparative study on the identification methods for the autoregressive modelling from the ambient vibration records, Soil Dynamics and Earthquake Engineering 15 (1) (1996) 45–49. [21] X. He, G. De Roeck, System identification of mechanical structures by a high-order multivariate autoregressive model, Computers and Structures 64 (1–4) (1997) 341–351. 1044 V.H. Vu et al. / Mechanical Systems and Signal Processing 25 (2011) 1028–1044 [22] C.S. Huang, Structural identification from ambient vibration measurement using the multivariate AR model, Journal of Sound and Vibration 241 (3) (2001) 337–359. [23] C.S. Li, W.J. Ko, H.T. Lin, R.J. Shyu, Vector autoregressive modal analysis with application to ship structures, Journal of Sound and Vibration 167 (1) (1993) 1–15. [24] A. Bjorck, Numerical Methods for Least Squares Problems, Society for Industrial and Applied Mathematics, Philadelphia, PA, 1996, 408 pp. [25] G. Golub, C. Van Loan, Matrix Computations, third ed., The Johns Hopkins University Press, London, 1996, 694 pp. [26] B.A. Cipra, The best of the 20th century: editors name top 10 algorithms, SIAM News 33 (4) (2000) 1–2. [27] R.J. Allemang, D.L. Brown, A correlation coefficient for modal vector analysis, in: Proceedings of the First International Modal Analysis Conference, Orlando, November 1982, pp. 110–116. [28] B.R. Mace, K. Worden, G. Manson, Uncertainty in structural dynamics, Journal of Sound and Vibration 288 (3) (2005) 423–429. [29] R. Pintelon, P. Guillaume, J. Schoukens, Uncertainty calculation in (operational) modal analysis, Mechanical Systems and Signal Processing 21 (6) (2007) 2359–2373. [30] A. Neumaier, T. Schneider, Estimation of parameters and eigenmodes of multivariate autoregressive models, ACM Transactions on Mathematical Software 27 (2001) 27–57. [31] H. Lutkepohl, Introduction to Multiple Time Series Analysis, second ed., Springer-Verlag, Berlin, 1993, 545 pp. [32] J. Paulsen, D. Tjostheim, On the estimation of residual variance and order in autoregressive time series, Journal of the Royal Statistical Society B 47 (1985) 216–228. [33] R.L. Kashyap, Inconsistency of the AIC rule for estimating the order of autoregressive models, IEEE Transactions on Automatic Control AC-25 (1980) 996–998. [34] E.J. Hannan, The estimation of the order of an ARMA process, Annals of Statistics 8 (5) (1980) 1071–1081. [35] G. Liang, D.M. Wilkes, ARMA model order estimation based on the eigenvalues of covariance matrix, Transactions on Signal Processing 41 (10) (1993) 3003–3009. [36] M. Smail, M. Thomas, A.A. Lakis, Assessment of optimal ARMA model orders for modal analysis, Mechanical Systems and Signal Processing 13 (5) (1999) 803–819. [37] V.H. Vu, M. Thomas, A.A. Lakis, L. Marcouiller, Multi-autoregressive model for structural output only modal analysis, in: Proceedings of the 25th Seminar on Machinery Vibration, Canadian Machinery Vibration Association, St. John, Canada, October 2007, pp. 41.1–41.20.