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

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: 𝑀̈̇(𝑡),𝑌(𝑡)(𝑡),𝑢+𝐷𝑢+𝐾𝑢=𝑓=𝐿𝑢 (1) 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 𝑥̇(𝑡)=𝐴_𝑐𝑥(𝑡)+𝐵_𝑐𝑓(𝑡),𝑦(𝑡)=𝐶_𝑐𝑥(𝑡), (2) where 𝐴_𝑐=^0𝐼_{−𝑀−1𝐾−𝑀−1𝐷}, 𝐵_𝑐=⁰_𝑀−1, 𝐶_𝑐=𝐿0 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 𝑥_𝑘+1=𝐴𝑥_𝑘+𝐵𝑓_𝑘,𝑦_𝑘=𝐶𝑥_𝑘, (3) 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: 𝑥_𝑘+1=𝐴𝑥_𝑘+𝑤_𝑘𝑦_𝑘=𝐶𝑥_𝑘+𝑣_𝑘, (4) 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 𝐸𝑤_𝑝𝑣_𝑝𝑤_𝑞𝑣_𝑞^𝑇=𝑆𝑄𝑆_𝑇𝑅𝛿_𝑝𝑞;∀𝑝,𝑞∈𝐙, (5) 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]: 𝑅_𝑚=𝐸𝑦_𝑘+𝑚𝑦^𝑇_𝑘. (6) An estimate of the covariance matrices is given [2] as follows: 𝑅_𝑚=1𝑁−𝑚^𝑁−𝑚_𝑘=1𝑦_𝑘+𝑚𝑦^𝑇_𝑘=𝐶𝐴^𝑚𝐺, (7) 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: 𝐻=𝐶⋮𝑝,𝑞𝐶𝐴𝐶𝐴^𝑝−1⋅𝐺𝐴𝐺𝐴²𝐺⋯𝐴^𝑞−1𝐺, (8) 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 𝐻𝑝,𝑞=𝑈⋅𝑆⋅𝑉^𝑇=𝑈₁𝑈₂𝑆₁00𝑆₂𝑉^𝑇₁𝑉^𝑇₂≅𝑈₁⋅𝑆₁⋅𝑉^𝑇₁=𝑈₁⋅𝑆₁^1/2𝑆₁^1/2⋅𝑉^𝑇₁=𝑂⋅𝐶, (9) 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: 𝑂=𝑈₁⋅𝑆₁^1/2,𝐶=𝑆₁^1/2⋅𝑉^𝑇₁. (10)

The system matrix 𝐶 is immediately extracted from the observability matrix by taking the “𝑟” first rows; “𝑟” denotes the numbers of sensors used in the structure: 𝐶=𝑂(1∶𝑟,∶). (11) The system matrix 𝐴 is extracted from the observability matrix as follows: 𝐴=𝑂_↓⁺𝑂_↑, (12) where 𝑂_↓=𝐶𝐶𝐴𝐶𝐴²⋯𝐶𝐴^𝑞−2𝑇, 𝑂_↑=𝐶𝐴𝐶𝐴²⋯𝐶𝐴^𝑞−1𝑇, and the superscript (+) means the generalized inverse.

The system matrix 𝐺 is immediately extracted from the controllability matrix by taking the “𝑟” first columns: 𝐺=𝐶(∶,1∶𝑟). (13)

4. Modal Analysis

With the assumption of low damping ratios and distinct eigenvalues, natural frequencies are obtained from 𝐴 as follows: 𝑓_𝑝=12𝜋ln𝜇_𝑝𝑇²+arg(𝜇_𝑝)𝑇². (14) The damping ratios are also obtained from 𝐴 as follows: 𝜉_𝑝=−ln𝜇_𝑝ln𝜇_𝑝²+arg(𝜇_𝑝)². (15) From 𝐶, the modes shapes are expressed as 𝜙_𝑝=C𝜓_𝑝, (16) where 𝜇_𝑝 and 𝜓_𝑝 are the eigenvalues and eigenvectors of 𝐴 and arg(⋅) 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 𝑓_min≤𝑓_𝑖≤𝑓_max and 𝜉_min≤𝜉_𝑖≤𝜉_max; 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: 𝜀_𝑓=1% for frequency, 𝜀_𝜉=5% for damping, and the modal assurance criterion (MAC) MAC_𝜙=98% for eigenvectors. Two modes identified in certain two orders, “𝑖” and “𝑖+1,” will be plotted in the stabilization diagram if 𝑓_𝑖−𝑓_𝑖+1𝑓_𝑖≤𝜀_𝑓,𝜉_𝑖−𝜉_𝑖+1𝜉_𝑖≤𝜀_𝜉,MAC𝜙_𝑖,𝜙_𝑖+1=𝜙^∗_𝑖⋅𝜙_𝑖+1²𝜙^∗_𝑖⋅𝜙_𝑖⋅𝜙^∗_𝑖+1⋅𝜙_𝑖+1≥MAC_𝜙, (17) 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 (1−MAC_𝜙), 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 𝑖∈{1,…,𝑝_𝑚}, 𝑗∈{1,…,𝑝_𝑛} as 𝑓_𝑖−𝑓_𝑗𝑓_𝑖≤𝜀_𝑓,𝜉_𝑖−𝜉_𝑗𝜉_𝑖≤𝜀_𝜉,MAC𝜙_𝑖,𝜙_𝑗=𝜙^∗_𝑖⋅𝜙_𝑗²𝜙^∗_𝑖⋅𝜙_𝑖⋅𝜙^∗_𝑗⋅𝜙_𝑗1≥MAC_𝜙. (18)

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 𝑓_𝑠=256 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 𝑝=𝑞=25, the maximal model order is taken as 𝑂_max=53, and the time lags used are 𝑁_lag=128. 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: 𝜀_𝑓=1% for frequencies, 𝜀_𝜉=5% for damping, and MAC_𝜙=98% 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: 𝑅_𝑚=𝐶𝐴^𝑚𝐺=(𝐶Ψ)Λ^𝑚Ψ⁻¹𝐺=ΦΛ^𝑚𝐺=Φ₁Φ₁Λ₁00Λ₁𝐺^𝑇₁𝐺^𝑇₁, (19) 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 𝐺𝐺=^𝑇1𝐺^𝑇₁=𝑔^𝑇₁𝑔^𝑇₂⋯𝑔^𝑇_𝑛𝑔^𝑇₁𝑔^𝑇₂⋯𝑔^𝑇_𝑛^𝑇. (20)

The stabilization criteria described up to now are derived from the vectors Φ if the MAC_Φ 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 (𝑔_𝑖,𝑔_𝑗): 𝑔^∗_𝑖⋅𝑔_𝑗²𝑔_𝑖^∗⋅𝑔_𝑖⋅𝑔^∗_𝑗⋅𝑔_𝑗≤1−MAC_𝑔. (21)

For all 𝑖∈{1,…,𝑝_𝑚}, 𝑗∈{1,…,𝑝_𝑛}, 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 𝜀_𝑓=1%, 𝜀_𝜉=3%, MAC_𝜙=99%, and MAC_𝑔=99.95%.

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.

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 𝑧=exp(𝑗𝜔𝑇) as follows: 𝐺(𝑧)=^+∞_−∞𝑅_𝑖𝑧^−𝑖=12Φ𝑔+Φ(𝑧𝐼−Λ)⁻¹+1Λ𝑔21Φ𝑔+Φ𝑧𝐼−Λ⁻¹Λ𝑔^𝑇. (22)

Proof. Consider 𝐺(𝑧)=^+∞_−∞𝑅_𝑖𝑧^−𝑖=⁰_−∞𝑅_𝑖𝑧^−𝑖+^+∞₀𝑅_𝑖𝑧^−𝑖=12𝑅₀+¹_−∞𝑅_𝑖𝑧^−𝑖+12𝑅₀+^+∞₁𝑅_𝑖𝑧^−𝑖. (23) Let the half PSD 𝐺₊(𝑧)=((1/2)𝑅₀+∑¹_−∞𝑅_𝑖𝑧^−𝑖), 𝐺(𝑧)=𝐺₊(𝑧)+(𝐺₊(1/𝑧))^𝑇; ^+∞₁𝑅_𝑖𝑧^−𝑖=^+∞₁𝐶⋅𝐴^𝑖⋅𝐺⋅𝑧^−𝑖=𝐶^+∞₁𝐴^𝑖−1𝑧^𝑖−1𝑧⁻¹𝐴⋅𝐺=𝐶^+∞₁𝐴𝑧^𝑖−1𝑧⁻¹⋅𝐴⋅𝐺. (24) The matrix 𝐴 is stochastic, so for all eigenvalue, |𝜆|≤1; therefore, we have ^+∞₁𝐴𝑧^𝑖−1=𝐴𝐼−𝑧⁻¹,𝐺₊(𝑧)=12𝑅₀+𝐶(𝑧𝐼−𝐴)⁻¹⋅𝐴⋅𝐺. (25) Finally, after substitution of the expression of the matrices into the identified base, we have 𝐺₊(𝑧)=12Φ𝑔+Φ(𝑧𝐼−Λ)⁻¹Λ𝑔. (26)

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: 𝐿_∞=𝑅₀𝑅₋₁𝑅₋₂⋯𝑅₁𝑅₀𝑅₋₁…𝑅₂𝑅₁𝑅₀⋯⋯⋯⋯⋱≻0. (27)(iii)The power spectral densities 𝐺(𝑧) are positive definite matrix for all 𝑧=exp(𝑗𝜔𝑇): 𝐺(𝑧)=12Φ𝑔+Φ(𝑧𝐼−Λ)⁻¹+1Λ𝑔21Φ𝑔+Φ𝑧𝐼−Λ⁻¹Λ𝑔^𝑇≻0. (28)

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.