Paper Ref: S1101_P0212 3rd International Conference on Integrity, Reliability and Failure, Porto/Portugal, 20-24 July 2009 OPERATIONAL MODAL ANALYSIS OF NON-STATIONARY MECHANICAL SYSTEMS BY SHORT-TIME AUTOREGRESSIVE (STAR) MODELLING V.H. Vu1, M. Thomas1, A.A. Lakis2 and L. Marcouiller3 1 Department of Mechanical Engineering, École de technologie supérieure, 1100 Notre Dame West, Montreal, QC H3C 1K3, Canada 2 Department of Mechanical Engineering, École Polytechnique, Montreal, QC H3C 3A7, Canada 3 Hydro-Québec’s research institute, Varennes, QC J3X 1S1, Canada Corresponding author: Email: marc.thomas@etsmtl.ca ABSTRACT In this paper, a method is developed for the modal analysis of vibrating structures whose properties may vary with time; it is based on an autoregressive model in a short- time scheme, and is called Short-Time Auto Regressive (STAR). This new method allows the successful modeling and identification of an output-only modal analysis of non-stationary systems. The originality of the proposed method lies in its specific handling of non-stationary vibrations, which allows the tracking of modal parameter changes in time. The paper presents an update of the model used with respect to order as well as a new criterion for the selection of an optimal model order, based on the noise-to-signal ratio. It is shown that this criterion works stably with respect to block data length. As for the block size itself, a length equal to four times the period of the lowest natural frequency has been found to be efficient. To validate the proposed method, a system with three degrees of freedom is first numerically simulated under a random excitation, and both linear and non-stationary vibrations are considered. In both cases, only a time-invariant autoregressive model is used in each sliding window, and the model parameters are estimated using least squares implemented via a QR factorization of the data matrix. The method is finally experimentally applied on the real multichannel data measured on a steel plate emerging from water, and compared to the conventional Short-Time Fourier Transform (STFT) method. It is shown that the proposed method outperforms in terms of frequency identification, whatever the non-stationary behaviour (either slow or abrupt change) due to the added mass effect of the fluid. Keywords: Non-stationary vibration; operating modal analysis; autoregressive model; shorttime modeling; structural change. 1. INTRODUCTION Operational modal analysis [Vu (2006); Mohanty (2004)] is a necessary method for modal parameters identification in structural testing and monitoring, especially when the structures operate under severe working conditions with a non linear behavior. In modern modal analysis, it is important to monitor changes in the modal properties of the structures over time [Basseville (1988); Uhl (2005)]. Non-stationary vibration in systems with time-dependent 1 properties may be analyzed through non-parametric and time-frequency methods [Bellizzi, 2001] such as the Short-Time Fourier transform (STFT), Wavelet or Wigner-Ville [Safizadeh (2000)]. However, parametric methods offer a number of advantages, such as improved accuracy and resolution [Poulimenos (2004), Owen (2001], which is why time domain methods are generally preferred. The last two decades have witnessed a trend toward the use of the time series model for non stationary modal analysis. Notable industrial applications [Hermans, 1999] include traffic-excited bridge vibration [Vu, 2007], seismic vibration, robotic devices, wind turbines [Gagnon, 2006], fluid-structure interaction [Thomas, 2005], etc. In non-stationary cases, most of these methods involve a nonlinear model [Fassois (2001), Hammond (1996), Basseville (1993)]. In this paper, we use the linear autoregressive model, which is sequentially computed in a short-time scheme [Risanen, 1978], by using a sliding window. The multivariate autoregressive model is used, and the model parameters are estimated by least squares using the QR factorization. Rather than varying the parameters on a sample by sample basis, we propose the use of a sliding window [Mahon (1993), Strobach (1993)], and keep the parameters constant inside each window. Consequently, the parameters are adjusted window by window. The proposed method may thus be considered as the time domain counterpart of the Short-Time Fourier Transform method, and can be used with multi-channel measurements. The model order is robustly selected from the noise-to-signal separation via the computation of a new criterion. The block size automatically adjusted from the lowest natural frequency. The method has shown its efficiency in tracking the modal parameters and monitoring their evolution on both linear and nonlinear structures. 2. VECTOR AUTOREGRESSIVE (VAR) MODELING Assuming a random environment, the excitation may be ignored, and since the modal analysis required multiple measurement locations, a vector autoregressive model, with d sensors, can be expressed as follows:  y (t )   a1  y (t  1)   a2  y (t  2)  ...   a p   y(t  p)  e(t ) (1) The model, as rewritten here below, is a convenient form of a dimension d vector at order p [Vu, 2009].  y (t )d 1   Ad dp  (t )dp1  e(t )d 1 (2) where  Ad dp     a1    a2  ...   ai  ...   a p   is the model parameter matrix,  ai d d  y(t ) , is the matrix of autoregressive parameters relating the output  (t )dp1   y(t  1) ;  y(t  2); vector  y (t ) , ...;  y (t  p)  y (t  i ) to is the regressor for the output 2  y(t  i )d 1 (i=1:p) is the output vector with delay time i  T , T is the sampling period (s), e(t )d 1 is the residual vector of all output channels, and considered as the error of the model. If the data are assumed to be measured in a white noise environment, the least squares estimation is applicable. If N successive output vectors of the responses from  y (t ) to  y (t  N  1) are considered (N ≥ dp+d), then the model parameters matrix  Ad dp and the estimated covariance matrices of the unnoised part  Dˆ  and of the error part  Eˆ  can be d d d d estimated via the computation of the QR factorization as follows [Vu, 2009]:  A  (  R12T   R11 ).(  R11T   R11 ) 1  (  R111   R12 )T (3) 1  Dˆ    R12T R12     N (4) 1 T  Eˆ    R22    N  R22  (5) In these formulas,  R11  ,  R12  and  R22  are sub-matrices of the upper triangular factor  R  derived from the QR factorization of the data matrix as follows:  K   Q . R  where  Q N  N is an orthogonal matrix (that is  Q  . QT    I  ), and  R  has the form  R N ( dp d )  R11 dpdp   0  0   R12 dpd    R22 dxd  0   (6) (7) and data matrix  K  is constructed from N successive samples:  K N  dp d    (t )Tdp    (t  1)T dp   ...  T  (t  N  1)dp   y (t )d T  y (t  1)d      ...  T y t N (   1)  d  T (8) 3 Once the model parameters matrix has been estimated, modal parameters can be directly identified from the eigen-decomposition of the state matrix    [Pandit, 1991].   dpdp    a1  d d   I   0   ...   0   a2 d d   ai d d 0 0 ... 0 0 ... 0 I  ... ... 0 0 ...   a p  I  ... d d ... 0         (9) 3. THE STAR METHOD In operational modal analysis, the dynamic parameters of the system are unknown, and thus, a priori knowledge about the model order is not available. Since we are concerned with shorttime modelling, we propose that the data be processed in block-wise Gabor expansion [Gabor, 1946]. From the above modelling, it is found that the number of samples in each block N must satisfy N ≥ dp+d, where p is the computing model order, and thus can be variable. It is also clear that the block size must be long enough to allow an exhibition of the vibratory features of the system and to cover the largest period in the signal. The block length of the sliding window must be adjusted from the greatest period. The optimal model must be selected from the order 2 to the maximum available order fitting data of all block sizes. Because it is time-consuming to repeat the computation for each order value, this procedure should be avoided. Below, we present an algorithm allowing an effective updating of the solution by model order, and it requires only the triangularization on a sub-matrix of the equations. 3.1 Order updating and new criterion for optimal model order selection The data matrix K ( p ) at order p can be rewritten as: K ( p ) N ( dp  d )   T (k )   T ( k  1)   ...  T  ( k  N  1)   y ( k  1)    K ( p)   1 N dp ...  y T (k  N  1)  y T (k ) T K2 N d  (10) If the model order is updated to p  1 , the data matrix has the form: K ( p1) N ( d ( p 1)  d )   K1( p ) N dp where K ' Nxd K ' N d K2 N d  (11) are the added d columns 4  yT (k  ( p  1))   T  y (k  1  ( p  1))      ...  T   y (k  N  1  ( p  1))  K ' N d (12) We can then compute the following matrix: Q ( p )T K ( p1)  Q( p )T K1( p ) N dp dp d where T1 and T2 ( N  dp ) d Q( p )T K ' N d Q ( p )T K 2 N d  R( p ) T1    11  0 T2 T  are extracted from Q ( p )T K '   1  . T2  R12( p )   R22( p )  (13) We must now triangularize the right term matrix in equation (13). This can be done with a set of Householder transformations or Givens rotations. If we decompose only the small submatrix T2 , it easily yields: R  T2  QT  T  0 (14) where RT d d is an upper diagonal matrix and QT Householder transformations or Givens rotations. Equation (13) then becomes: Q ( p )T K  I dpdp  0  ' where R22 ( p 1) N ( d ( p 1)  d )  I dpdp   0  R( p) 0  ( p )T ( p 1)  11 Q K  0 QTT   0  '' d  d and R22 ( N dpd )d  R( p ) 0   11 0 QT    0 T1 RT 0 T1 RT 0 R12( p )  '  R22  ''  R22  ( N  dp )( N  dp ) is the product of the   QTT  R22( p )    R12( p ) (15) (16)  R'  are obtained from multiplication  22  QTT R22 . "  R  22  It can be seen that the first d  p rows of the right hand side in equation (16) are not affected by the above transformations, and the factor matrix R ( p1) at order p  1 is thus updated: ( p 1) 11 R  R( p)   11  0 T1  ; ( p 1)  R12  ; R ( p 1)  R"  '  22 22  R12 RT   R22  ( p) (17) as is the Q matrix: 5  I dpxdp Q ( p 1)  Q ( p )   0 0 QT  (18) as well as the two covariance matrices from (4) and (5): T T T T ' ' '  R22  R22'  Dˆ ( p )   R22 Dˆ ( p 1)   R12( p 1)  R12( p1)   R12( p )  R12( p )   R22 (19) T T T ( p 1) T ( p 1) '' '' ( p) T ( p) ' ' '  R22   R22    R22  R22   Eˆ ( p )   R22   R22  (20) Eˆ ( p 1)   R22   R22   R22'   R22 Finding the optimal model order pmin is crucial in parametric model-based methods [Smail (1999), Gang (1993), Hannan (1980)]. From a statistical point of view, AIC and MDL criteria can be used to select the optimal model order [Lutkepohl, 1993]. It is seen from equations (19) and (20) that as the model order increases, the norm of the deterministic covariance matrix increases as well, while one of the error parts decreases by the same amount. The global noise-to-signal ratio is:  NSR  eˆ(t ) Tr ( Eˆ )  yˆ(t ) Tr ( Dˆ ) (21) It is therefore monotonically decreased in terms of the model order. The Noise-rate Order Factor (NOF) defines the change in the NSR within two successive order values:   NOF ( p )  NSR ( p )  NSR ( p 1) (22) It is seen that the NOF is always positive, falling quickly at low orders and converging at high orders, which is why it serves as a criterion for selecting the model order. The optimal order should therefore be chosen after a significant change in the NOF. The NOF converges to a minimum value after the optimum order is reached. Since the solution is effectively updated, the selection of optimal order can be applied for time-varying systems. 3.2 Working procedure Since there is no leakage present in parametric model-based modal analyses, it is therefore not necessary to use a window, and so we propose that the data be processed in combination with a progressive search for the model order, as follows. Firstly, the above VAR model is initially applied to a block of data with a reasonable low order value. The length of the first block size could be specified if the smallest natural frequency of interest of the structure is known. Modal parameters are identified and the natural frequencies are estimated by using the signalto-noise ratio of each eigen-value (MSN) [Pandit, 1991] in order to find the smallest frequency to use to specify the length of the next block. Once the block size is chosen, the optimal model order is selected by the NOF and an order equal to or higher than this minimum value is used to get the modal parameters. The overlapping process can also be employed by changing the sliding step, which can vary from only one sample to the whole of the block window. 6 4. NUMERICAL SIMULATION 4.1 Simulation on mechanical system A numerical simulation of the proposed method was applied on a 3 degrees of freedom (d.o.f.) system, as shown in Fig. 1, under a random force excitation. The mechanical properties of the system are first considered constant, and present three natural frequencies at 6.4 Hz, 12.6 Hz and 25 Hz, and damping rates at 2%, 4% and 7.8%, respectively. Fig. 1 Three d.o.f mechanical system 4.2 Discussion on block data length Fig. 2 shows the frequencies and damping rates of the above structure identified by the VAR method, with varying data lengths. a) Frequencies b) Damping rates Fig. 2 Modal parameter identification with block size 7 We can observe that the block size must be larger than 3 times the longest period Tmax in order to produce the smallest natural frequency value. For that reason, the block size was chosen to be equal to 4 times the period. The results found were obtained from multiple simulations. In conclusion, when a moderate damped system is subjected to a random excitation, its modal parameters can be monitored at block sizes equal to 4 times the period of the smallest natural frequency. Fig. 3 plots the NOF curves of the 3 d.o.f. system at various data length sizes. It is seen that the structural model order is found accurately at 3 regardless the data length showing the stability of the NOF criterion with respect to data length. It also confirms that if the properties of the structure are subjected to change, the optimal order can still be tracked, and does not depend on the block size once this latter is long enough to show the smallest frequency. Fig. 3 Optimal model order at different block sizes 4.3 Simulation on mechanical system with time-dependent parameters The 3 d.o.f. system above has been modified to vary its mechanical properties, and is always subjected to a random excitation. The mass M2 is now a time variant factor which changes following the function, as shown in Fig. 4. Fig. 4 Simulated time-varying mass 8 Since the data is non-stationary, the optimal model and modal parameters can vary with time, and therefore, the block size should be changed. If the smallest frequency is available either through prior knowledge of the system or from the frequency range of interest, the initial block size is chosen to be 4 times its period. When more data are acquired, the block size is adjusted based on the smallest frequency identified in the previous step. Figures 5 and 6 show the optimal model and the computing block size to use to track the change in the system properties. It is shown that the optimal order is primarily monitored at 3, except for some outer values. We can observe an adjustment on the block length when the change appears at the half end of the monitoring period. Fig. 5 Monitoring of optimal order on simulation Fig. 6 Monitoring of block size on simulation The changes in frequencies and damping rates are plotted on Figures 7 and 8. As is logically seen when the mass is increasing, all frequencies tend to decrease, and can be well tracked. On the other hand, we see a high variance for the damping rates. The damping identification only gives an approximate range of values. 9 a) Mode 1 b) Mode 2 c) Mode 3 Fig. 7 Monitoring of frequencies on simulation 10 a) Mode 1 b) Mode 2 c) Mode 3 Fig. 8 Monitoring of damping rates on simulation 11 5. EXPERIMENTAL APPLICATION ON AN EMERGING STEEL PLATE The method is applied to the monitoring of vibration seen in a submerged steel plate. The plate measures 500 mm x 200 mm x 2 mm, and emerges from water while, and is always excited by a random turbulent flow. Fig. 9 shows the configuration of the test while Fig. 10 presents a temporal response data, with the low amplitude portion corresponding to the submerging period and the high amplitude portion attributable to the emergence of the plate from the water to the air. Before the plate rises, its modal parameters are both calculated and identified using analytical and experimental methods [Vu, 2007], as shown in Table 1. Tab. 1 Modal identification of the plate Frequency (Hz) in submerging conditions (Depth/plate length ratio) 0.6 (totally submerged) 11.9 0.4 0.2 0.1 12.0 12.2 12.7 0 (totally in air) 39.4 nd 34.1 34.1 34.2 35.0 75.0 rd 77.7 77.9 78.1 79.5 108.6 4th 135.3 135.4 135.6 137.5 164.0 5th 151.3 151.3 151.4 152.6 210.0 Mode 1st 2 3 Sensor 1 100 Sensor 2 Sensor 5 Sensor 3 150 150 Sensor 4 100 Sensor 6 Fig. 9 Plate test configuration 12 Fig. 10 Plate temporal response Fig. 11 plots the optimal model order applied to the data over time. It is seen that the optimal order is found primarily from 3 to 6. Fig. 12 shows the block size over time for computing purposes, which follows the same pattern than the frequency variation. Fig. 11 Optimal model order Fig. 12 Computing block size Fig. 13 shows the monitoring of frequencies computed from order pcom = pmin+ 10, with the variations clearly revealed. The changes in frequencies correspond to the emergence of the 13 plate. Both the slow change when the plate was still in water and the abrupt change when it appears on the surface are monitored. The natural frequencies highly match the calculated values in Table 1, showing that the effect of added mass on the plate is accurately monitored [Vu, 2007]. Compared to the STFT in Fig. 14, which is computed on the first channel with the same parameters, it is seen that the proposed STAR method outperforms in terms of revealing the natural frequencies. Fig. 13 Monitoring of frequencies Fig. 14 Short-time Fourier transform (STFT) Unfortunately, the results on the damping rate identification were not significant, and research is still ongoing on this subject. 14 6. CONCLUSIONS An application of the multivariate autoregressive model in the short-time scheme (STAR) was presented, and covered the monitoring of changes in the modal parameters of the structure under unknown random excitation. In order to track the frequency variations, it is not suitable to use a stabilisation diagram and it is preferred to select an optimal order for computing the modal parameters. It was found that the minimum structural model order does not depend on the block size if this latter is long enough. The block size was minimally found to be equal to four times the wavelength of the first natural frequency, and consequently, block sizes vary with variations in the first natural frequency. Numerical simulations and experiments show that the proposed method can be used to track a slow change as well as a sudden change in the modal properties of the structure, and that it outperforms the STFT method. While the monitoring of natural frequencies has been successfully dealt with, those of damping rates are not enough precise. Research is still ongoing on damping identification of time-varying system as it is found in fluid-structure interaction. 7. ACKNOWLEDGEMENTS The support of NSERC (Natural Sciences and Engineering Research Council of Canada) through Research Cooperative grants is gratefully acknowledged. The authors would like to thank HydroQuebec’s Research Institute for its collaboration. 8. REFERENCES 1. Basseville M., Benveniste A., Gach-Devauchelle B., Goursat M., Bonnecase D., Dorey P., Prevosto M., Olagnon M., 1993, Damage monitoring in vibration mechanics: issues in diagnostics and predictive maintenance, Mechanical Systems and Signal Processing, Vol. 7, No. 5, Sept., pp. 401-423. 2. Basseville M., May 1988, Detecting changes in signals and systems - A survey, Automatica, Vol. 24, No. 3, pp. 309-326. 3. Bellizzi S., Guillemain P., Kronland-Martinet R., 2001. Identification of Coupled Nonlinear Modes from Free Vibrations using Time-Frequency Representations. Journal of Sound and Vibrations, Vol. 243 (2), pp. 191-213. 4. Fassois S.D., 2001. Parametric identification of vibrating structures. In: S.G. Braun, D.J. Ewins and S.S. Rao, Ed., Encyclopedia of Vibration, Acad. Press, N. York, pp. 673-685. 5. Gabor D., 1946. “Theory of Communication,” J. IEEE (London), Vol. 93, pp. 429-457. 6. Gagnon M., Tahan S.-A., Coutu A. and Thomas M., 2006. Operational modal analysis with harmonic excitations: application to a wind turbine. Proceedings of the 24th Seminar on machinery vibration, Canadian Machinery Vibration Association, ISBN 2-921145-618, Montreal, pp. 320-329. 7. Gang Liang, Wilkes D. M. & Cadzow J. A., 1993. ARMA Model Order Estimation Based on the Eigenvalues of Covariance Matrix. Transactions on Signal Processing, Vol. 41, No 10, pp. 3003-3009 8. Hammond J.K. and White P.R., 1996, The analysis of non stationary signals using timefrequency methods, Journal of Sound and Vibration, Vol. 190, pp. 419-447. 9. Hannan E. J., 1980. The estimation of the order of an ARMA process. The Annals of Statistics, Vol. 8 (5), pp. 1071-1081. 10. Hermans L. and Van der Auweraer H., 1999. Modal testing and analysis of structures under operational conditions: Industrial applications. Mechanical Systems and Signal Processing 13(2), pp. 193-216. 15 11. Lutkepohl H., Introduction to Multiple Time Series Analysis (2nd ed.). Springer-Verlag, Berlin, 1993, 545 p. 12. Mahon M., L. Sibul, and H. Valenzuela, 1993. A sliding window update for the basis matrix of the QR-decomposition, IEEE Trans. Signal Processing, Vol. 41, pp. 1951-1953. 13. Mohanty P. and Rixen D. J., February 2004. Operational modal analysis in the presence of harmonic excitation. Journal of Sound and Vibration, 270 (1-2, 6), pp. 93-109. 14. Owen J.S., B.J. Eccles, Choo B.S. and Woodings M.A., 2001. The application of autoregressive time series modelling for the time-frequency analysis of civil engineering structures, Engineering Structures 23, pp. 521–536. 15. Pandit S.M., 1991, Modal and spectrum analysis: data dependent systems in state space. New York, N.Y., J. Wiley and Sons, 415 p. 16. Poulimenos A.G. and Fassois S.D., 2004. Non stationary vibration modelling and analysis via functional series TARMA models. 5th International conference on acoustical and vibratory surveillance methods and diagnostic techniques, Surveillance5, Senlis, 10 p. 17. Rissanen, J. 1978. Modeling by shortest data description. Automatica, 14, pp. 465-471. 18. Safizadeh, M.-S., Lakis, A.A. et Thomas, M, 2000, Using Short Time Fourier Transform in Machinery Fault Diagnosis, International Jour. of Condition Monitoring and Diagnosis Engineering Management , 3 (1), pp 5-16. 19. Smail M., Thomas M. and Lakis A.A., 1999. Assessment of optimal ARMA model orders for modal analysis, Mechanical systems and Signal Processing journal, 13 (5), 803-819. 20. Strobach P. and D. Goryn, 1993. A computation of the sliding window recursive QR decomposition, Proc. ICASSP, pp. 29-32. 21. Thomas M., Abassi K., Lakis A. A. and Marcouiller L., October 2005. Operational modal analysis of a structure subjected to a turbulent flow. Proceedings of the 23rd Seminar on machinery vibration, Canadian Machinery Vibration Association, Edmonton, AB, 10 p. 22. Uhl, T., 2005. Identification of modal parameters for non-stationary mechanical systems. In: Arch. Appl. Mech. 74, pp. 878–889. 23. Vu V.H., M. Thomas, A.A. Lakis and L. Marcouiller, May 2009. Online monitoring of varying modal parameters by operating modal analysis and model updating. CIRI2009, Reims, France, 18 p. 24. Vu V.H., Thomas M., Lakis A.A. and Marcouiller L., April 2007. Identification of modal parameters by experimental operational modal analysis for the assessment of bridge rehabilitation. Proceedings of International operational modal analysis conference (IOMAC 2007). Copenhagen, Denmark, pp. 133-142. 25. Vu V.H., Thomas M. and Lakis A.A., 2006. Operational modal analysis in time domain. Proceedings of the 24th Seminar on machinery vibration, Canadian Machinery Vibration Association, ISBN 2-921145-61-8, Montreal, pp. 330-343. 16