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

Themodal parameters of a structure that is estimated from ambient vibrationmeasurements 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.


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 (SSIcov) [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.

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 , ,  ∈ R × , 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 ).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  =     ,  = (−) −1    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 V  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.

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 C. 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;  1 is the dominant singular bloc matrix  2 ≺≺  1 .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:  = C (:, 1 : ) . (13)

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 arg(⋅) denotes the phase angle.

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 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.

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).
For a certain model orders "" and "",  ̸ = , having "  " and "  " identified modes, the criteria will be transformed for all  ∈ {1, . . .,   },  ∈ {1, . . .,   } 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.Plot poles "i" and "j"

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.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 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.

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 Ψ −1  = G 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 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 (  ,   ): 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%.
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.

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.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.
Consequently, the results comparison of the BR algorithm identification with the others confirms that our results are satisfactory.condition [3], in case, the inequality (28).Therefore, the extended covariance matrices [10] might not be positive.Thus, the spectrum divergences are well justified.

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.

Figure 4 :Figure 5 :
Figure 4: The aircraft stabilization diagram with classical implementation.

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

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