Shock and Vibration

Volume 2014 (2014), Article ID 409298, 8 pages

http://dx.doi.org/10.1155/2014/409298

Research Article

A New Criterion for the Stabilization Diagram Used with Stochastic Subspace Identification Methods: An Application to an Aircraft Skeleton

E. Mrabet,¹ M. Abdelghani,² and N. Ben Kahla²

¹High Institute of Technological Studies of Kairouan, Campus Universitaire Raccada, 3100 Kairouan, Tunisia

²High Institute of Applied Sciences and Technologies of Sousse, University of Sousse, Cité Taffala (Ibn Khaldoun), 4003 Sousse, Tunisia

Received 5 July 2012; Accepted 20 November 2012; Published 17 March 2014

Academic Editor: Sami El-Borgi

Copyright © 2014 E. Mrabet et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

The modal parameters of a structure that is estimated from ambient vibration measurements are always subject to bias and variance errors. Accordingly the concept of the stabilization diagram is introduced to help users identify the correct model. One of the most important problems using this diagram is the appearance of spurious modes that should be discriminated to simplify modes selections. This study presents a new stabilization criterion obtained through a novel numerical implementation of the stabilization diagram and the discussion of model validation employing the power spectral density. As an application, an aircraft skeleton is used.

1. Introduction

The vibration and acoustical behaviors of a mechanical structure are determined by its dynamic characteristics.

This dynamic behavior is typically described with a linear system model. The procedure for the estimation of modal parameters of structures from measured data can be split into three distinct steps [1]: data collection, system identification, and determination of modal parameters from the identified system description.

Stochastic identification methods for systems with unknown input have been introduced decades ago. Among the most robust and accurate system identification methods for output-only modal analysis of mechanical structures is the stochastic subspace identification method. Two types of implementation are available: the covariance-driven (SSI-cov) [2] implementation and the data-driven (SSI-data) [3] implementation. For the first one (SSI-cov), three methods can be implemented: the balanced realization (BR), the principal component (PC), and the canonical variate analysis (CVA).

For dynamic structures such as the aircraft skeleton studied in this paper, the major setback in applying system identification for large-scale structures is the selection of the model order and the corresponding system poles.

To address this problem, the concept of the “stabilization diagram” is introduced, overestimating the structure model order. Therefore, spurious modes are going to surface out and we have to discriminate them. For this matter, many stabilization criteria have been implemented. The most recent one was the modal transform norm [4]. In this paper, a new stabilization criterion is implemented and a validation method is discussed. The stochastic subspace identification method used is the balanced realization.

2. Stochastic State Space Models for Vibrating Structures

The finite element method [4] is one of the most common tools for modeling mechanical structures. In the case of a linear dynamical model, one has the following system of ordinary differential equations:

where , and are the mass, damping, stiffness, and selection matrices, respectively; is the stochastic vector of nodal forces; is the vector of nodal displacement; is the sensors measurements vector; “” is the number of degree of freedom. In the case of nonproportional damping, (1) is written into a continuous time state space model. The classical form of this model is

where , , and .

The vector is called the state of the structure; is the measurements vector. The discrete time state model of the mechanical system is expressed as

where , and is the time sampling. For simplicity, the model given by (3) can be also written, when the output vector is noised, as follows:

where and are, respectively, the state transition and the output matrices; is the discrete state vector; is the discrete measured output vector; and are, respectively, the process and measurement noises.

These stochastic terms are unknown but it is assumed that they have a discrete white noise nature with an expected value equal to zero and that they have covariance matrices equal to

where denotes the Kronecker delta.

3. Identification of the System Matrices and Balanced Realization Algorithm

The starting point of the identification of the system matrices is based on the covariance matrices of the measured structural responses time series which are assumed to be realization of a stationary stochastic process. The covariance matrices are given by the following formula [2]:

An estimate of the covariance matrices is given [2] as follows:

where is the number of points of the time series, , and superscript means transpose. These covariance matrices can be organized in a Hankel matrix that can be factorized as follows:

where and denote the number of rows and columns, respectively, of the Hankel matrix. The first bloc matrix is called the observability matrix ; the second one is the controllability matrix . For all covariance-driven stochastic subspace identification, the estimation of and is based on the singular values decomposition of the Hankel matrix [5, 6].

For the balanced realization method, the decomposition gives

where , , and are the left singular vector, the singular values matrix, and the right singular vector, respectively; is the dominant singular bloc matrix . Afterwards, the observability and controllability matrices can be written as the following formulas:

The system matrix is immediately extracted from the observability matrix by taking the “” first rows; “” denotes the numbers of sensors used in the structure:

The system matrix is extracted from the observability matrix as follows:

where , , and the superscript (+) means the generalized inverse.

The system matrix is immediately extracted from the controllability matrix by taking the “” first columns:

4. Modal Analysis

With the assumption of low damping ratios and distinct eigenvalues, natural frequencies are obtained from as follows:

The damping ratios are also obtained from as follows:

From , the modes shapes are expressed as

where and are the eigenvalues and eigenvectors of and denotes the phase angle.

5. The Stabilization Diagram

The SSI-cov, as it happens with the balanced realization method, does not yield exact values for the parameters but only estimates with uncertainties. The origins of these uncertainties [4] can be described from two points of view:

(i) From the Operational/Experimental Point of View: the number of data samples is finite, the input might not be white noise and real structure can’t always be modeled into stationary linear system.(ii) From a Statistical Point of View: uncertainties can be induced either by the bias of the model or by the bias of the mode and the variance of the mode.

These uncertainties are responsible for the appearance of spurious modes. One of the primary challenges related to modal analysis is to remove these spurious modes. For this purpose, the concept of stabilization diagram is introduced.

The stabilization diagram is a graphical tool used to help in the manual selection of the modes that are more likely to represent the structure physical modes. The quality of a stabilization diagram [2] depends on the algorithm used in the identification, on the values of the input parameters of the algorithm, and also on noise ratio of the time series under analysis. The basic idea is that several runs of the complete pole identification process are made, by using models of increasing order. The stabilization diagram has frequencies in the horizontal axis and orders on the vertical one. Physical poles should readily appear into alignments. However, not only physical modes will come into view in this diagram but also spurious modes. Most of spurious modes can be removed by using the so-called stabilization criteria. These criteria can be split into two types.

(i) Preliminary Criteria. The frequencies and damping ratios are expected to be obtained in certain ranges and ; all modes having frequencies and damping ratios out of these ranges will be discriminated. Only modes verifying these criteria are plotted in the stabilization diagram.(ii) Stabilization Criteria. Modes, for which the differences in modal parameters between two consecutive model orders [4] are higher than certain threshold values, are not in the diagram.

In the classical implementation of the stabilization diagram, the typical stability criteria values [3] are as follows: for frequency, for damping, and the modal assurance criterion (MAC) for eigenvectors. Two modes identified in certain two orders, “” and “,” will be plotted in the stabilization diagram if

where subscript (*) denotes complex conjugate transpose.

This manner of plotting the stabilization diagram is not the only one. In [7], the authors present a new approach to establish the diagram. The method estimates model order in terms of component energy index (CEI); then the diagram is plotted in increasing Hankel matrix rows.

6. An Alternative Numerical Implementation for the Stabilization Diagram

When working with the stabilization diagram, the perfect situation is an utter disappearance of all spurious modes and only the alignments (corresponding to physical system modes) are plotted in the diagram. Using the classical numerical implementation and taking into account the bias and variance modes, if we minimize the thresholds values , , and , most of the spurious modes will disappear; nonetheless, the diagram is likely to lose some alignments.

This present study suggests an alternative numerical implementation of the stabilization diagram that intends not to calculate the difference between consecutive modal order parameters but to compare every modal order parameter with all the other modal orders (see Figure 1).

Figure 1: Flow chart, stabilization diagram implementation.

For a certain model orders “” and “”, , having “” and “” identified modes, the criteria will be transformed for all , as

The basic idea behind this algorithm is that in the simulated cases when the output is not noised, if a mode appears into an order of the diagram, it should show up for all the following orders too.

7. Application: The Aircraft Skeleton

To put this numerical implementation with the BR method into operation, an aircraft Skeleton is considered. The test rig and the model of the aircraft skeleton structure [8] are shown in Figure 2. The structure is excited by a white noise. Only seven sensors measurements and a sampling frequency  Hz are available for the identification.

(a)

(b)

Figure 2: (a) The test rig; (b) the schema of the aircraft skeleton structure [8].

Figure 3 shows the aircraft stabilization diagram with only preliminaries stabilization criteria. The diagram is plotted with frequencies and damping ranges (0, 120 Hz) and (0, 15%), respectively. The numbers of rows and columns of the Hankel matrix are , the maximal model order is taken as , and the time lags used are . In this diagram, 14 alignments are seen and seeming to be the physical system modes.

Figure 3: The aircraft stabilization diagram with only preliminaries criteria.

Figure 4 shows the stabilization diagram plotted with the following criteria: for frequencies, for damping, and for eigenvectors. The diagram is plotted by using the classical numerical implementation. The inspection of this diagram shows that there is discrimination of several spurious modes, certain alignments are affected (frequencies close to 3 and 80 Hz), and the modes selections become difficult. Using the suggested numerical implementation with the same values for the criteria, the stabilization diagram is shown in Figure 5. The alignments are already existent in spite of the fact that certain spurious modes are already there.

Figure 4: The aircraft stabilization diagram with classical implementation.

Figure 5: The stabilization diagram with new implementation.

8. How to Improve the Stabilization Diagram Quality?

In order to make the stabilization diagram cleaner, a criterion is introduced in this study. It is based on the expression of the covariance matrix that can be written in the identified base as follows:

where are the modes shapes expressed in the identified base and is the covariance matrix between the state and the system output expressed also in the same base; it takes the form

The stabilization criteria described up to now are derived from the vectors if the is used and from the eigenvalue matrix if frequencies and damping criteria are employed.

The new criterion is based on the calculation of the MAC (modal assurance criterion) between two identified covariance vectors ():

For all , , these covariance vectors should be the same, at least theoretically, for every identified mode onto an alignment in the stabilization diagram. For the aircraft diagram shown in Figure 6, the criteria were , , , and .

Figure 6: The stabilization diagram with the new stabilization criterion.

In Figure 6, it is clear that 12 alignments are already stable in spite of the fact that the criteria are considered severe and most of the spurious modes have been discriminated. The structure has likely 12 vibration modes.

9. Balanced Realization (BR) Results and Comparison with Other Stochastic Subspace Identification Algorithms

For the model order 24, stable modes, belonging to alignments, are taken and the identification is presented in Table 1. In order to validate this system identification, the balanced realization results are compared with others obtained by two different algorithms implemented into a commercial modal analysis software, the ARTeMIS Extractor pro 2009 [9]. The algorithms are the canonical variate analysis (CVA) and the unweighted principal component (UPC), which are two data-SSI methods. Table 1 shows that both BR and CVA methods present the same identified parameters. Among the 12 modes identified by these two algorithms, 9 are also identified in frequencies by the UPC method in which 7 are also identified in damping.

Table 1: The identified frequencies and damping using several algorithms.

Modes
Frequencies: (Hz)
Damping: (%)
BR/UPC (%)
BR/CVA (%)

BR
CVA
UPC
BR
CVA
UPC

1
3.71
3.69
***
4.90
6.57
***
***
***
0.54
34.08

2
6.13
6.18
***
3.60
3.17
***
***
***
0.82
11.94

3
16.06
16.02
16.07
4.05
4.07
3.94
0.06
2.54
0.25
0.49

4
21.71
21.77
***
5.73
5.40
***
***
***
0.28
5.76

5
39.87
39.87
39.86
0.38
0.38
0.40
0.02
5.26
0.00
0.00

6
44.33
44.34
44.37
0.03
0.02
0.06
0.09
106.66
0.02
33.33

7
54.35
54.36
54.37
0.11
0.10
0.09
0.03
16.36
0.02
9.09

8
55.42
55.42
55.39
0.22
0.20
0.31
0.05
40.90
0.00
9.09

9
58.44
58.47
58.49
0.38
0.39
0.47
0.08
23.68
0.05
2.63

10
***
***
63.83
***
***
3.93
***
***
***
***

11
69.41
69.37
69.42
1.12
1.18
1.25
0.01
11.60
0.06
5.36

12
80.06
80.15
79.8
0.09
0.09
0.72
0.32
700
0.11
0.00

13
***
***
86.18
***
***
3.00
***
***
***
***

14
100.25
100.3
100.3
0.11
0.11
0.12
0.04
9.09
0.05
0.00

15
***
***
107.6
***
***
4.51
***
***
***
***

16
***
***
114.5
***
***
1.63
***
***
***
***

***Unidentified mode.

These differences between the first two algorithms and the third are understood given the fact that all of them can generate uncertainties. Moreover, it is known that the identification of the damping ratio is difficult even in simulated cases.

Figure 7 shows the repartition of the damping ratios corresponding to the stabilization diagram of Figure 6 and the dashed lines represent the first four identified modes. The inspection of Figure 7 shows that the damping identification presents a large dispersion along the model order range, mainly, for the first identified mode (3.71 Hz), in which the damping ratio varies between, about, 4% and 14% which is considered as large dispersion. These facts justify the relative damping error (BR/CVA) for the first mode (34.08%) which is considered as acceptable.

Figure 7: Repartition of damping ratios.

Consequently, the results comparison of the BR algorithm identification with the others confirms that our results are satisfactory.

10. Spectral Analysis

10.1. Expression of the Power Spectral Densities Matrix

Once a state space is identified, it is possible to compute the power spectral densities which are written into the identified base for as follows:

Proof. Consider

Let the half PSD , ;

The matrix is stochastic, so for all eigenvalue, ; therefore, we have

Finally, after substitution of the expression of the matrices into the identified base, we have

10.2. The Positive Real Sequences Conditions

A sequence cannot be always considered as a valid output covariance sequence [3]. The sequence has to satisfy the positive real sequence condition. The following statements are equivalent:

(i) is a positive real sequence.(ii)The double infinite matrix is positive definite:

(iii)The power spectral densities are positive definite matrix for all :

11. Validation of the Identification by Spectral Analysis

A second procedure to validate the system identification is the spectral analysis. Four spectra are presented in Figures 8, 9, 10, and 11. They are compared with those obtained by applying FFT to measurements.

Figure 8: Cross spectrum between sensors 3 and 4.

Figure 9: Autospectrum sensor 5.

Figure 10: Cross spectrum between sensors 1 and 5.

Figure 11: Cross spectrum between sensors 3 and 7.

From these comparisons, it is clear that the structure model parameters are well identified because the model and data spectra fit in spite of the divergence in the autospectrum of sensor 5 and the cross spectrum between sensors 3 and 7. Figures 8 and 10 show that both peaks and valleys are well superposed which implies a perfect identification of the model. The inspection of Figure 9 shows that all peaks and valleys are well superposed except the valley close to the Nyquist frequency (120 Hz). The same remark can be done for Figure 11 where model and data peaks are well superposed but not all the valleys over the frequencies range. In effect, model valleys on the frequencies range [40 Hz, 55 Hz] and close to the Nyquist frequency do not fit the data spectrum valley.

From these facts, we can conclude that peaks in all data and model spectra are better identified than the valleys. This can be understood from the implementation of the SSI-cov method that does not guarantee the positive real sequence condition [3], in case, the inequality (28). Therefore, the extended covariance matrices [10] might not be positive. Thus, the spectrum divergences are well justified.

12. Conclusions

A new stabilization criterion associated with a new numerical implementation of the stabilization diagram has been presented using the BR (SSI-cov) method.

The new stabilization implementation makes the alignments in the diagram more robust and only spurious modes were removed. However, some other spurious modes are not removed and modes selection is still difficult. To remedy to this fact, the new stabilization criterion, based on the calculation of the modal assurance criterion (MAC) between the identified covariance vectors, was associated with the new numerical implementation. The obtained diagram is cleaner, most of the spurious modes were removed, and only alignments corresponding to physical modes are already there; then the modes selections are easier.

The validation by comparison with other results derived through the CVA and UPC algorithms gave satisfactory results in spite of some large errors in the damping identification. This fact has been justified in the previous section and we can conclude that the BR algorithm was a robust identification method when used with the new numerical implementation associated with the new criterion.

The validation by spectral analysis also confirms the quality of the model identification and the divergences into some spectra have been well justified.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgment

The authors are grateful to S. D. Fassois of the Stochastic Mechanical Systems (SMS) Group of the University of Patras (Greece) for providing the experimental measurements.

References

E. Reynders and G. D. Roeck, “Reference-based combined deterministic-stochastic subspace identification for experimental and operational modal analysis,” Mechanical Systems and Signal Processing, vol. 22, no. 3, pp. 617–637, 2008.

F. Magalhães, Á. Cunha, and E. Caetano, “Online automatic identification of the modal parameters of a long span arch bridge,” Mechanical Systems and Signal Processing, vol. 23, no. 2, pp. 316–329, 2009.

P. Van Overschee and B. De Moor, Subspace Identification for Linear Systems: Theory, Implementation and Applications, Kluwer Academic Publishers, 1996.

E. Reynders, R. Pintelon, and G. De Roeck, “Uncertainty bounds on modal parameters obtained from stochastic subspace identification,” Mechanical Systems and Signal Processing, vol. 22, no. 4, pp. 948–969, 2008.

F. Magalhães, E. Caetano, and Á. Cunha, “Operational modal analysis and finite element model correlation of the Braga Stadium suspended roof,” Engineering Structures, vol. 30, no. 6, pp. 1688–1698, 2008.

D. Montalvão, N. M. M. Maia, and A. M. R. Ribeiro, “A review of vibration-based structural health monitoring with special emphasis on composite materials,” Shock and Vibration Digest, vol. 38, no. 4, pp. 295–324, 2006.

Y. Zhang, Z. Zhang, X. Xu, and H. Hua, “Modal parameter identification using response data only,” Journal of Sound and Vibration, vol. 282, no. 1-2, pp. 367–380, 2005.

V. Papakos and S. D. Fassois, “Multichannel identification of aircraft skeleton structures under unobservable excitation: a vector AR/ARMA framework,” Mechanical Systems and Signal Processing, vol. 17, no. 6, pp. 1271–1290, 2003.

ARTeMIS, “ARTeMIS Extractor pro 2009. s.l.,” Denmark, SVIBS, 2009, http://www.svibs.com/

A. Lindquist and G. Picci, “Canonical correlation analysis, approximate covariance extension, and identification of stationary time series,” Automatica, vol. 32, no. 5, pp. 709–733, 1996.