Spectrum construction for non stationary vibration: Application to a moving flexible robot

This paper presents a method for constructing spectrum in non-stationary vibration using time series vector autoregressive model. A modal classification criterion and a modal amplification factor are introduced based on the eigendecomposition of the block data matrix. The classical spectrum computation can be therefore modified to amplify only the physical modes and lighten the computation induced modes. In combination with a sliding window technique, it is shown that the method provides a clear, lightweight spectral representation of non stationary vibration. Application on the vibration signals measured from a moving flexible robot along a given trajectory shows the effectiveness and applicability of the method. The technique is integrated in the online modal surveillance software call STAR (Short Time AutoRegressive).


Introduction
Experimental vibration has witnessed tremendous advancement in the recent decades. For industrial facilities and when excitation is unknown, the experimentalists are nowadays required to conduct dynamics and vibration tests on structures and machines in real working conditions. The methods based on operational vibration require signal processing developments. As the in-working measurement is continuous, the operational experimental vibration analysis need to be able to treat the signal with high level of perturbations, continuous, and most of the time with nonstationarities.
To date, the most convenient, fast and usable technique for processing and constructing the spectra of non stationary vibration signal is the Short Time Fourier Transform (STFT) [1]. However, STFT exhibits some drawbacks in dealing with short time windowing such as leakage and frequency resolution [2]. On the other hand by using autoregressive models, the spectra can be computed either from signal-noise separation, from the transfer function or from spectral decomposition. Akaike [3] early computed the power spectra through an autoregressive modelling. Since then, several techniques have been developed on the multichannel spectrum [4] or for improving the resolution of the spectra of the autoregressive models [5,6]. The most advanced developments on modal identification of non-stationary vibrations can be found in the recent collection on the identification of time varying structures and systems [7].
In this paper, a technique for modal decomposition and spectrum construction of non stationary vibration signal is presented. The vector autoregressive model is used for fitting the short time signal with the sliding window manner. Model parameters are estimated by the a Corresponding author: zhaoheng.liu@etsmtl.ca least squares via the construction of the QR factorization of the block data matrix. In the modal decomposition, a modal signal-to-noise based classification criterion is introduced in order to reorganize the modes. From that modal decomposition, a closed form of the spectral matrix is constructed with the insertion of an amplification in order to emphasize the physical modes and lighten the spurious modes.

Vector Auto-Regressive (VAR) for short time modelling
In operational modal analysis, the excitation is assumed unknown and can be modelled by a Gaussian white noise. As modal analysis is conducted by using several d channels of measurements, synchronized for data acquisition at a sampling period T s , the vector time series such as a Vector Auto-Regressive (VAR) or Vector Auto-Regressive Moving Average (VARMA) can be used for modelling the data. The modal parameters are therefore extracted from the autoregressive part of those models. Since the excitation is assumed Gaussian, the least squares estimation is assumed as non-biased and a Vector Auto-Regressive (VAR) model of p th order and of dimension d is sufficient for to fit the measured data [8]: where: is the matrix of parameters relating the output y(t − i) to y(t); z(t) = [y(t − 1) T y(t − 2) T . . . y (t − p) T ] T (size dp × 1) is the regressor for the output vector y(t); y(t − i) (size d × 1; i = 1 : p) is the output vector with delay time i × T s ; MATEC Web of Conferences e(t) (size d × 1) is the residual vector of all output channels considered as the error of model.
If the data are assumed to be measured in a white noise environment, the least squares estimation may be applied for estimating the model parameters. Taking into account N successive available output vectors of the responses from y(k) to y(k + N − 1) (k > p, N > dp + d), the model parameters matrix and the estimated covariance matrices of the deterministic partD and of the error partÊ (both of size d × d) can be estimated via the computation of the QR factorization as follows [8]: In these formulas, R 11 (size dp × dp), R 12 (size dp × d) and R 22 (size d × d) are sub-matrices of the upper triangular factor R (size N × dp + d) derived from the QR factorization of the data matrix as follows: where Q (size N × N ) is an orthogonal matrix (that is Q × Q T = I), R has the form: and data matrix K of size N × dp + d is constructed from N successive samples: Once the model's parameters matrix has been estimated, the modal parameters can directly be identified from the eigen-decomposition of the state matrix (size dp × dp) [9]: where λ i are discrete eigenvalues and L (size dp × dp) is the eigenvectors matrix, whose forms can be rewritten as follows for further usage: The continuous eigenvalues, natural frequencies, damping ratios and complex modes are derived as follow: Eigenvalues: Frequencies: Damping ratios: Complex modes:

Construction of the damped signal-to-noise ratio
In any model based modal analysis, the signal can be decomposed into a deterministic part and a stochastic part from the eigen-decomposition as follows [9]: where: y(t) is the measured sample data at time t, of dimension d × 1; d is the number of sensors, or vector dimension; p is the model order in the modelling; l i are taken from eigenvector matrix L; s i is the scale factor which can be computed from the initial data as: λ i are discrete eigenvalues; S i1 are extracted from the inverse eigenvector matrix; e(t) is the model error vector sampled at time t.

AVE2014
It is seen that the first term of Eq. (16) is the deterministic part and the second one is the stochastic part of the sampled signal. Those parts are the sum from the contribution of all eigenvalues and eigenvectors. If we consider the contribution of each mode to the deterministic part, we can take the modal power of the discrete signal over the N available sampled data: The contribution of each mode to the stochastic part can be computed in term of the modal variance, taken from Pandit [9]: whereÊ is the estimated error matrix of the underlying parametrical model.
Since the frequency and damping of each mode are identified, the factor called Damped Modal Signal to Noise (DMSN) ratio is proposed as follows [10]: The DMSN index is an effective criterion since it includes the stochastic participation in the denominator, and hence, the higher the DMSN, the more evident it characterizes the contribution to the deterministic part. The damping rate on the denominator penalizes the very high damped modes which usually belong to computational modes. The damping rate is more usually found ranging from 0.5% to 7% in industrial applications when no external damping is added. The modes are all sorted from the first highest DMSN index, with reasonable damping ratios. Consequently, the harmonic frequencies, if present, are first revealed, and can then be distinguished from the natural frequencies by their close-to-zero damping ratio, and hence a very high DMSN index. The available natural frequencies then follow and the computational frequencies end up with a very low DMSN.

Construction of an amplified spectrum.
The construction of a noise-free spectrum is of great interest for modal analysis and mechanical systems surveillance. If we assume that 2n modes, in conjugate pairs (eigenvalues and eigenvectors), are selected from the deterministic part, the spectral matrix function can be computed directly from the eigen-decomposition as follows: where: It is clear that the spectrum is composed of the sum over frequencies, which only exhibits the noise-free peaks corresponding to the natural frequencies and harmonic excitations, if present. This decomposition yields to a Hermitian multi-channels spectral matrix, and may be considered as a generalization of Pandit's formulations [9]. Further calculations on the multi-channels spectral matrix, such as the channels coherence function and phase can be found to be of interest. However, a difficulty persists when a low-amplitude peak is located close to a higher one, and is thus more difficult to identify. To exhibit all the frequency peaks in the spectrum representation, an innovation is introduced on an amplification factor for each frequency by dividing its participating amplitude by the real norm of the complex matrix d i (Eq. (24)). It is seen 01004-p.3   from this formulation that all the selected frequencies are revealed at equivalent amplitudes on the frequency plot:

Numerical simulations
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 at [M 1 , M 2 , M 3 ] = [3,1,2] kg, [C 1 , C 2 , C 3 , C 4 ] = [5, 10, 10, 5] Ns/m, K 1 , K 2 , K 3 , K 4 ] = [5000,10000,10000,5000] N/m. The system presents three natural frequencies at 6.4 Hz, 12.6 Hz and 25 Hz, and damping rates at 2%, 4% and 7.8%, respectively. In order to simulate the non-stationary vibration of the system, the stiffness K 1 is set to follow a linear variation as shown in Fig. 2. The stiffness K 1 is constant at 5000 N/m for the first 5 s and then increases linearly to 13000 N/m at 20 s. The simulated result is presented on Fig. 3 where the signal is stationary during the first 5 s and non-stationary during the last 15 s.
The 5 s stationary signal is considered first for proceeding. The vector autoregressive model order is estimated by using the NOF factor [8] which shown that any model order from 3 could be used for the modal identification. As seen even at the minimum model order 3, the identified modal parameters are given in Table 1 where a very good accordance is found with the simulated values.   The spectra of the AR model are first constructed on Fig. 4 at the order 3 by using some well-known methods such as the Yule, Burg, Covariance and Modified Covariance [11][12][13][14]. It is seen that at low model order, none of those methods are able to show all the three peaks. Figure 5 shows the spectrum performed by the proposed method at the same order 3 where the black curve is the raw spectrum as computed by Eq. (22) and the red curve is the spectrum as computed by the Eq. (25) with the introduction of the amplification factor. It is seen that even at the minimum model order 3, the method is able to capture all the three peaks. With the amplification factor, the spectrum is enhanced where the peaks are better displayed and separated from the neighbour peaks.
In this section, the proposed method is applied on the whole vibration signal including the non-stationary part when the stiffness K 1 varies linearly. The sliding window technique is combined with the method. The window length is chosen at 1 s (200 samples) with 50% overlapping in order to be able to identify the lowest frequency at 6.4 Hz, according to suggestion in [15] where the window length is suggested to be at least 4 times of the longest natural period. Because of the nonstationarities in the signal, the VAR model is chosen at order 10. Figure 6 displays the frequency identification and monitoring accordingly with the variation of the spring K 1 in comparison with the theoretical results (solid lines). It is seen that under random excitation, the method is able to identify and follow the change on the frequencies.
On the other hand, Fig. 7 shows the magnitude squared coherence (MSC) [12,16] of the three DOFs auto PSD. It is seen clearly that the spectrum displays all the 3 frequencies and reflects correctly their evolution in 01004-p.4

Application to a flexible robot in movement
A flexible robot named SCOMPI has been developed at Hydro Quebec's Research Institute (IREQ) for grinding process. Since the robot is used for performing the grinding processes, it is subjected to various dynamic excitations. The dynamic characterization and modal analysis on the structure are thus required for analysing its dynamic behaviour in the operational working conditions. The proposed method is applied to the operational modal analysis steps of the robot, while considering the dynamic motion of the robot as a random excitation. The signals are acquired from the flexible robot moving along a predefined trajectory without grinding operation interfaced at the endeffector. A translation displacement scheme is selected from its break position (initial position) to its deployment position (last position) as shown in Fig. 8.
A constant deployment speed was selected at approximately 1 cm/s. A 3-axial accelerometer PCB-352C21 is mounted at the robot's end-effector and 5 uniaxial accelerometers PCB-352C34 are mounted at the joint and mid-links (8 channels). The LMS Test.lab system is used for the data acquisition at the sampling frequency of 512 Hz.
The signal is seen non-stationary on Fig. 8 because of the continuous change on the geometric configuration of the robot during the process. The signal starting from the 60 th second to the end is seen as free vibrations with the damped amplitudes when the robot arrived to its final position.
The proposed method is applied on the 8 channels accelerations signal. A constant sliding window length at 2 s is applied with 50% overlapping. A VAR at model   order 20 is used throughout the signal. Figure 9 shows the magnitude squared coherence (MSC) of the 8 channels.
It is seen that the spectrum reveals clearly the variation of the change on frequencies during the non stationary process of the robot. We can also observe the electric frequencies at 60 Hz, 120 Hz and 180 Hz. As for comparison, Fig. 10 shows the Short Time Fourier Transform (STFT) of the first channel signal at the same window length 2 s and 50% overlapping.
It is seen that the proposed spectrum is better on the frequency identification and resolution. Parallel with the spectrum where the amplitudes of the auto power can be of interest, the natural frequencies can be also precisely identified and monitored during the process, as shown on Fig. 11. As expected, it is found that the deployment process of the robot arm implies the continuous variation of the natural frequencies of the structure. The first frequency decreases from 8 Hz to 6 Hz while the 2 nd and the 3 rd frequencies increase from 25 Hz to 33 Hz and from 27 Hz to 37 Hz respectively, which are well matched with the spectrum evolution.
This present algorithm is integrated into a software program called STAR for modal surveillance on mechanical systems in operation [17][18][19].

Conclusions
The ongoing research is aimed at developing a numerical model which is able to simulate the working process of the robot during operation. The stiffness of joints and the flexibility of the robot links will be taken into account in order to evaluate their effect on the modal parameters of the robot and therefore on its dynamic behaviour. A method for the computation and construction of the spectrum from the vector autoregressive model applying in operational modal analysis of non stationary vibration was presented. With the introduction of the modal classification criterion and the modal spectrum amplification factor, the spectrum of the multivariate autoregressive can be applied to the multi-channels vibration measurement during the non-stationary process. The method reveals a fast, peaks distinct, free of noise and harmonized time-frequency representation which outperforms the conventional SFFT method and the modified covariance method.
The current research work is being focused on the development and integration of the finite element model updating approach with the method proposed in this paper in order to identify time-varying modal parameters.