Philippe Mellinger
email: philippe.mellinger@epfedu.fr

Michael Döhler
email: michael.doehler@inria.fr

Laurent Mevel
email: laurent.mevel@inria.fr

Variance estimation of

Keywords: (Operational) modal analysis, Subspace algorithms, Variance estimation, Aircraft application

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction

The design and maintenance of mechanical, aeronautical or civil structures subject to noise and vibrations are an important topic in structural engineering. Laboratory and in-operation vibration tests are performed on structures for a modal analysis [START_REF] Heylen | Modal Analysis Theory and Testing[END_REF]. The desired modal parameters (frequencies, damping ratios, mode shapes) are related to the poles and observed eigenvectors of a linear time-invariant system. Subspace-based linear system identification algorithms have been proven efficient for the identification of such systems, fitting a linear model to either outputonly measurements or input and output measurements taken from a system [START_REF] Van Overschee | Subspace Identification for Linear Systems: Theory, Implementation, Applications[END_REF]. Usually, ambient excitation sources (e.g. wind, traffic) act on a structure during real tests. Besides these non-measurable excitation sources, in addition some measurable excitation sources (e.g. artificial excitation from hammer impact or shaker, ground motion) may be available. Thus, inputs to the considered system are in general unknown, but in some cases, some known inputs may be available in addition to the unknown inputs.

Operational Modal Analysis (OMA) refers to the identification of modal parameters from output-only measurements under unknown ambient excitation. In Operational Modal Analysis with eXogenous inputs (OMAX), some known or measured excitation is applied to the structure besides the ambient noise. Throughout this paper, the known input sequence is assumed to be a realization of a stochastic process. For example, control surface excitations can be used during either flight or ground vibration tests in the particular case of an aircraft. A stochastic framework can model the excitation of flaps actioned by a pilot during in-flight vibration tests, which are the known inputs in this case. Other unknown ambient excitation has to be considered, which is commonly assumed to be white noise.

The measured outputs are due to both known and unknown inputs, and may be corrupted by noise (i.e. sensor noise uncorrelated with known inputs).

Many methods for OMA and OMAX are available in the literature [START_REF] Reynders | System identification methods for (operational) modal analysis: review and comparison[END_REF]. In this work, we focus on subspace methods. By appropriate projections of the measurement data [START_REF] Van Overschee | Subspace Identification for Linear Systems: Theory, Implementation, Applications[END_REF], the state space representation of the corresponding linear system is properly estimated. Then, the modal parameters are retrieved from the estimated state space model.

There are a variety of subspace algorithms, which are defined by the way the measurements are projected. Data-driven subspace algorithms are characterized by large Hankel matrices filled with raw data. Orthogonal or oblique projections are carried out on their row space [START_REF] Van Overschee | Subspace Identification for Linear Systems: Theory, Implementation, Applications[END_REF][START_REF] Viberg | Subspace-based methods for the identification of linear time-invariant systems[END_REF], for which numerical schemes with the QR decomposition [START_REF] Van Overschee | Subspace Identification for Linear Systems: Theory, Implementation, Applications[END_REF] and iterative and fast approaches [START_REF] Mastronardi | A fast algorithm for subspace state-space system identification via exploitation of the displacement structure[END_REF][START_REF] Gandino | Identification of a duffing oscillator under different types of excitation[END_REF][START_REF] Döhler | Fast multi-order computation of system matrices in subspace-based system identification[END_REF] have been studied. In the stochastic framework, those projectors are essentially conditional expectations of the future output data with respect to some instrument depending on the algorithm [START_REF] Benveniste | Nonstationary consistency of subspace methods[END_REF]. Connections with Prediction Error Methods [START_REF] Ljung | System identification[END_REF] and closed loop identification have been investigated [START_REF] Jansson | A linear regression approach to state-space subspace system identification[END_REF][START_REF] Chou | Subspace algorithms for the identification of multivariable dynamic errors-in-variables models[END_REF][START_REF] Gustafsson | Subspace identification using instrumental variable techniques[END_REF]. It is well understood that the projections in the data-driven methods can be partly expressed in terms of auto-and cross-covariances of the respective time series [START_REF] Van Overschee | Subspace Identification for Linear Systems: Theory, Implementation, Applications[END_REF][START_REF] Gandino | Covariance-driven subspace identification: A complete input-output approach[END_REF]. The class of covariance-driven subspace algorithms is purely based on auto-and cross-covariance matrices of the measurements and has been used for modal identification [START_REF] Peeters | Reference-based stochastic subspace identification for output-only modal analysis[END_REF] or damage detection [START_REF] Basseville | Statistical model-based damage detection and localization: subspace-based residuals and damage-tonoise sensitivity ratios[END_REF][START_REF] Döhler | Subspace-based fault detection robust to changes in the noise covariances[END_REF] using vibration data for more than a decade. Like with data-driven methods, covariance-driven methods can be based on output-only data or input/output data [START_REF] Gandino | Covariance-driven subspace identification: A complete input-output approach[END_REF][START_REF] Mevel | Merging sensor data from multiple measurement set-ups for non-stationary subspacebased modal analysis[END_REF][START_REF] Guyader | Covariance-driven subspace methods: input/output vs. output-only[END_REF][START_REF] Mevel | Input/output versus output-only data processing for structural identificationapplication to in-flight data analysis[END_REF][START_REF] Döhler | Modular subspace-based system identification from multi-setup measurements[END_REF].

The purpose of this paper is the variance analysis of modal parameters obtained from subspace methods. Using noisy measurement data, subspace algorithms provide parameter estimates that are afflicted with statistical uncertainty due to finite data, unknown inputs and sensor noise properties. In-depth asymptotic variance analyses for subspace algorithms have been carried out in the automatic control literature [START_REF] Viberg | Analysis of state space system identification methods based on instrumental variables and subspace fitting[END_REF][START_REF] Bauer | Consistency and asymptotic normality of some subspace algorithms for systems without observed inputs[END_REF][START_REF] Chiuso | The asymptotic variance of subspace estimates[END_REF][START_REF] Bauer | Asymptotic properties of subspace estimators[END_REF], which are of important theoretical value.

In these works, analytical expressions for the asymptotic covariance of the system matrices are derived for several output-only and input/output subspace algorithms. However, explicit expressions for the actual covariance estimation from measurements are not available therein.

In the field of vibration analysis, such explicit expressions have been proposed for some subspace methods. For example, an application of the variance analysis in [START_REF] Viberg | Analysis of state space system identification methods based on instrumental variables and subspace fitting[END_REF] for input/output subspace identification was tempted in [START_REF] Raffy | Statistical asymptotic error on modal parameters in combined deterministic-stochastic identification algorithm[END_REF] in the context of modal analysis, however without entirely satisfying results according to the authors of the study.

A different and successful approach was developed in [START_REF] Pintelon | Uncertainty calculation in (operational) modal analysis[END_REF], where an estimated covariance on the measurements is propagated to the desired parameters based on a sensitivity analysis. The required sensitivities are derived analytically through the propagation of a first-order perturbation from the data to the identified parameters. This approach has the advantage of computational convenience: the sample covariance as the starting point is directly linked to the measurements and therefore easy to compute, and the sensitivities are computed using the system identification estimates. In [START_REF] Reynders | Uncertainty bounds on modal parameters obtained from stochastic subspace identification[END_REF], details of this scheme are developed for the covariance computation for output-only covariance-driven subspace identification. This approach offers advantages beyond technical correctness. The exploitation of the stabilization diagram, which is a basic tool for OMA(X), is a crucial step in modal analysis, but computationally costly for the covariance computation. A well-adapted fast and memory efficient implementation in an operational modal analysis context has been proposed in the covariance-driven output-only case [START_REF] Döhler | Efficient multi-order uncertainty computation for stochastic subspace identification[END_REF]. As a result, the computational cost for the variance computation of an entire stabilization diagram analysis is reduced significantly. This was mandatory for realistic applications on large structures like in [START_REF] Döhler | Uncertainty quantification for modal parameters from stochastic subspace identification on multi-setup measurements[END_REF][START_REF] Döhler | Structural health monitoring with statistical methods during progressive damage test of S101 Bridge[END_REF] and a prerequisite for implementation in modal analysis software.

While this approach has shown to be working both in theory and practice for the variance computation of modal parameters, it is only available for output-only covariance-driven subspace identification [START_REF] Reynders | Uncertainty bounds on modal parameters obtained from stochastic subspace identification[END_REF][START_REF] Döhler | Efficient multi-order uncertainty computation for stochastic subspace identification[END_REF]. In the present paper, this covariance computation approach is generalized to a wider class of subspace algorithms in both outputonly and input/output frameworks, and for covariance-driven and data-driven methods, while keeping the advantages of the previously published method in [START_REF] Döhler | Efficient multi-order uncertainty computation for stochastic subspace identification[END_REF] in terms of efficiency and precision. Notice that, while the computation of the modal mass is possible with the input/output algorithms [START_REF] Reynders | Reference-based combined deterministic-stochastic subspace identification for experimental and operational modal analysis[END_REF], we focus on the covariance computation for frequencies, damping ratios and mode shapes. In this paper, new covariance computation schemes are developed for four well-known subspace algorithms. One of them is an output-only algorithm:

• Output-only orthogonal projection data-driven, which is the UPC algorithm (Unweighted Principal Component) [START_REF] Van Overschee | Subspace Identification for Linear Systems: Theory, Implementation, Applications[END_REF]Sec. 3.3.2], [START_REF] Peeters | Reference-based stochastic subspace identification for output-only modal analysis[END_REF]. Furthermore, we consider three input/output algorithms:

• Input/output covariance-driven [START_REF] Ashari | Input-output subspace-based fault detection[END_REF],

• Input/output orthogonal projection data-driven, related to the projection algorithm [2, Sec. 2.3.2] and MOESP (Multivariable Output Error State Space) [START_REF] Verhaegen | Application of a subspace model identification technique to identify lti systems operating in closed-loop[END_REF], and

• Input/output oblique projection data-driven, which is the N4SID algorithm (Numerical algorithms for Subspace State Space System IDentification) [START_REF] Van Overschee | Subspace Identification for Linear Systems: Theory, Implementation, Applications[END_REF][START_REF] Van Overschee | N4SID: Subspace algorithms for the identification of combined deterministic-stochastic systems[END_REF].

The extension from the output-only covariance-driven algorithm to these four subspace algorithms is performed by a thorough sensitivity analysis of the underlying auto-and cross-covariance matrices in the formulation of the subspace algorithms. In particular, the extension to the data-driven algorithms is carried out by developing an equivalent formulation in a covariance-driven form, regarding the estimation of frequencies, damping ratios and mode shapes.

The new algorithms are validated in extensive Monte-Carlo simulations and in an application on an aircraft. This paper is organized as follows. In Section 2, the underlying vibration modeling is recalled. In Section 3, the general properties of subspace identification algorithms are described and the five considered subspace algorithms are detailed. An equivalent covariance-driven formulation of the data-driven algorithms is discussed in Section 4, which is the basis for the subsequent derivation of the covariance computation for the data-driven algorithms. Then variances for the considered algorithms are derived in Section 5. The validation of the methods in Monte-Carlo simulations is presented in Section 6 in an academic example. The application on an aircraft structure is presented in Sections 7 and 8 using simulations and real data, respectively, before discussing and concluding the paper in Sections 9 and 10. The vibration behavior of a linear time-invariant mechanical structure, which is observed at some sensor positions, can be described by the equations [1]

Nomenclature

⊗

         M q(t) + C q(t) + Kq(t) = F u(t) + ũ(t) y(t) = L a q(t) + L v q(t) + L d q(t) + v(t) (1)
where M, C and K ∈ R m×m are mass, stiffness and damping matrices, and m is the number of degrees of freedom. The vector q(t) ∈ R m×1 contains the displacements at the degrees of freedom generated by the known inputs

u(t) ∈ R N i ×1
and unknown inputs ũ(t) ∈ R m×1 , where N i is the number of known inputs, and F ∈ R m×N i represents how the known inputs are applied to the system. The vector y(t) ∈ R N o ×1 contains the observed outputs, with N o being the number of sensors. The matrices L a , L v and L d ∈ R N o ×m represent how accelerations, velocities and displacements are obtained from the model degrees of freedom. The vector v(t) ∈ R N o ×1 is the sensor noise. Both ũ(t) and v(t) are assumed to be white noise with finite fourth moments and uncorrelated with the known inputs [START_REF] Van Overschee | Subspace Identification for Linear Systems: Theory, Implementation, Applications[END_REF][START_REF] Verhaegen | Identification of the deterministic part of MIMO state space models given in innovations form from input-output data[END_REF].

From Eq. ( 1), the following continuous-time state space model can be obtained:

         ẋ(t) = A c x(t) + B c u(t) + Bc ũ(t) y(t) = C c x(t) + D c u(t) + Dc ũ(t) + v(t), (2)
where x(t) = [ q(t) T q(t) T ] T ∈ R n×1 is the state vector, n = 2m is the model order and

A c =           -M -1 C -M -1 K I m 0 m,m           ∈ R n×n , B c =           M -1 F 0 m,N i           ∈ R n×N i , C c = L v -L a M -1 C L d -L a M -1 K ∈ R N o ×n
and D c = L a M -1 F ∈ R N o ×N i are the state transition, input, output and feedthrough matrices, respectively. The matrices Bc and Dc are defined similarly with respect to the unknown inputs. The matrices I a ∈ R a×a and 0 a,b ∈ R a×b are the identity and the zero matrix, respectively.

Sampling Eq. ( 2) at rate τ and assuming zero-order hold for inputs u(t) and ũ(t) yield the discrete-time state space

representation [9]          x k+1 = Ax k + Bu k + Bũ k y k = Cx k + Du k + Dũ k + v k , (3)
where

A = e A c τ , B = (A -I)A -1 c B c , B = (A -I)A -1 c Bc , C = C c , D = D c
, and D = Dc . The modal parameters of system (1) are equivalently found in system (3) as follows. Let λ i and φ i be eigenvalues and eigenvectors of A, for i = 1, . . . , n. Then the poles µ i of systems (1) and (2), modal frequencies f i , damping ratios ξ i and mode shapes ϕ i are obtained by

µ i = log(λ i ) τ , f i = |µ i | 2π , ξ i = -Re(µ i ) |µ i | , ϕ i = Cφ i .

Stochastic nature of known inputs

Until now, no assumption has been stated on the known input {u k } k=1,...,N for the case of input/output algorithms.

Known inputs are used in what is called deterministic [START_REF] Van Overschee | Subspace Identification for Linear Systems: Theory, Implementation, Applications[END_REF] or combined deterministic/stochastic [START_REF] Van Overschee | Subspace Identification for Linear Systems: Theory, Implementation, Applications[END_REF][START_REF] Reynders | Reference-based combined deterministic-stochastic subspace identification for experimental and operational modal analysis[END_REF] system identification in the literature. In the designation of these methods, there is equivalence between "deterministic" and "the input sequence is known", and between "stochastic" and "there are unknown inputs and/or noise".

However, these designations should not be mistaken with the nature of the known inputs. The known input sequence itself may be a deterministic signal, or it may be a realization of a stochastic process. For the considered input/output methods in this paper, we assume that known inputs are realizations of a stochastic process satisfying mild ergodicity and mixing properties. This assumption is in particular coherent with the ergodic-algebraic framework of [START_REF] Verhaegen | Identification of the deterministic part of MIMO state space models given in innovations form from input-output data[END_REF] and cross-covariance estimates in [START_REF] Gandino | Covariance-driven subspace identification: A complete input-output approach[END_REF] for covariance-driven input/output algorithms. In combination with standard white noise assumptions for unknown inputs and measurement noise, this assumption allows the consistent computation of auto-covariances of the inputs, the outputs and of cross-covariances between inputs and outputs, which is a basic step in subspace identification.

These auto-and cross-covariances contribute to the modal parameter estimates. Then, the corresponding uncertainties for these auto-and cross-covariances play together in a non-trivial way to obtain the modal parameter uncertainties, for which the expressions are derived in this paper. Note that assuming known inputs as stochastic is no contradiction for the particular case of a system with no unknown inputs and no sensor noise. In this case, the modal parameters can be identified exactly from the known inputs and outputs, and their variance assuming stochastic inputs is indeed zero, as will be shown in Section 6.

Subspace identification

Subspace algorithms are based on projections of output or input/output data Hankel matrices with a "past" and "future" time horizon. These projections are designed in a way that the column space of the resulting matrix -we call it subspace matrix H -is defined by the observability matrix Γ p = [C T . . . (CA p-1 ) T ] T of system (3). All five subspace algorithms considered in this paper have the respective factorization property H = Γ p Z, where Z depends on the chosen algorithm [START_REF] Van Overschee | Subspace Identification for Linear Systems: Theory, Implementation, Applications[END_REF]. This property leads to the estimation of Γ p from measurement data and subsequently to the system matrices A and C, from where the desired modal parameters are identified.

After defining the underlying data Hankel matrices and the respective auto-and cross-covariances, the three steps of modal identification with subspace algorithms are detailed in this section: 1/ Estimation of auto-and crosscovariance matrices related to "past" and "future" data matrix products for the chosen subspace algorithm, 2/ Computation of the subspace matrix H for the five subspace algorithms considered in this paper, and finally 3/ Estimation of the modal parameters from H.

Definition 1 (Data Hankel matrix [START_REF] Van Overschee | Subspace Identification for Linear Systems: Theory, Implementation, Applications[END_REF]). Let N t samples a l ∈ R r×1 of a discrete signal a be given, and integers i and j with i ≤ j ≤ N t -N + 1. Then the data Hankel matrix A i| j is defined as

A i| j def =                   a i a i+1 . . . a i+N-1 . . . . . . . . . . . . a j a j+1 . . . a j+N-1                   ∈ R ( j-i+1)r×N .
Let N + p + q be the number of available samples for the known inputs u k and the outputs y k , where q and p are parameters that define a "past" and "future" time horizon. They are most often equal, and assumed to be large enough to satisfy the condition min{(p -1)N o , qN o } ≥ n. Under the above Definition 1, the "past" ( -) and "future" ( + ) data

Hankel matrices containing the known data are defined as in [START_REF] Van Overschee | Subspace Identification for Linear Systems: Theory, Implementation, Applications[END_REF],

U -def = 1 √ N U 0|q-1 , U + def = 1 √ N U q|q+p-1 , Y -def = 1 √ N Y 0|q-1 , Y + def = 1 √ N Y q|q+p-1 . (4)
Consider similarly the unknown signals, i.e. the unknown input ũ, the sensor noise v and the state x, defining the corresponding matrices

U -def = 1 √ N U 0|q-1 , U + def = 1 √ N U q|q+p-1 , V -def = 1 √ N V 0|q-1 , V + def = 1 √ N V q|q+p-1 and finally X -def = 1 √ N X 0|0 , X + def = 1
√ N X q|q , which are needed for deriving the equations of the model. Then Eq. ( 3) can be extended to [

2]                  Y -= Γ q X -+ H q U -+ Hq U -+ V - Y + = Γ p X + + H p U + + Hp U + + V + X + = A q X -+ ∆ q U -+ ∆q U -. (5)
The matrices ∆ q , ∆q , H p and Hp are defined in [START_REF] Van Overschee | Subspace Identification for Linear Systems: Theory, Implementation, Applications[END_REF].

First step: estimation of auto-and cross-covariances

Known data (outputs and, if available, inputs) are used to estimate auto-and cross-covariances in order to compute the subspace matrix H. Let a and b be placeholders for two of these signals (e.g. u and y, y and y or u and u). Then, their cross-covariance estimate at lag i is defined as

R [a,b] i = 1 N N k=1 a k+i b T k , (6)
considering that a and b are ergodic and zero mean random processes. The term auto-covariance is used for the crosscovariance of the same signal, i.e. a = b in [START_REF] Gandino | Identification of a duffing oscillator under different types of excitation[END_REF]. Note that due to the finite data length, these auto-and cross-covariance estimates are random variables with a covariance [START_REF] Hannan | Multiple time series[END_REF], which will be the starting point of the variance propagation to the modal parameter estimates in Section 5.

Besides the direct computation as in Eq. ( 6), an efficient way to compute auto-and cross-covariance functions at all lags i is to first estimate the PSD (Power Spectral Density) functions of signals by Welch's method [START_REF] Ljung | System identification[END_REF] and to compute the inverse discrete Fourier transform of the result, using the Fast Fourier Transform (FFT), as detailed in Appendix A.

Using "past" and "future" data Hankel matrices as defined in Eq. ( 4) brings some useful properties, since their products are related to the auto-and cross-covariance estimates. For example, it holds

U + Y -T =                           R [u,y] q R [u,y] q-1 . . . R [u,y] 1 R [u,y] q+1 R [u,y] q . . . R [u,y] 2 . . . . . . . . . . . . R [u,y] q+p-1 R [u,y] q+p-2 . . . R [u,y] p                           , (7)
and similarly

U -Y -T =                           R [u,y] 0 R [u,y] -1 . . . R [u,y] 1-q R [u,y] 1 R [u,y] 0 . . . R [u,y] 2-q . . . . . . . . . . . . R [u,y] q-1 R [u,y] q-2 . . . R [u,y] 0                           , U + Y + T =                           R [u,y] 0 R [u,y] -1 . . . R [u,y] 1-p R [u,y] 1 R [u,y] 0 . . . R [u,y] 2-p . . . . . . . . . . . . R [u,y] p-1 R [u,y] p-2 . . . R [u,y] 0                           (8)
The products Y -Y -T , Y + Y + T , Y + Y -T , U -U -T , U + U + T and U + U -T are computed analogously. These products based on auto-and cross-covariances are elementary blocks for the computation of the subspace matrices needed by the subspace algorithms of Section 3.2.

Remark 2 (Cross-covariance with unknown inputs and sensor noise). Products between data Hankel matrices containing outputs or known inputs and noise Hankel matrices are established analogously as above. These crosscovariance based products will appear in the derivation of subspace algorithms in the following section. Some of these cross-covariances tend to zero due to the assumption that unknown inputs and sensor noise are white and uncorrelated with the known inputs. In particular it holds U + U +T = 0 and V + U +T = 0 for an infinite amount of data.

Under these assumptions, past outputs are uncorrelated with future unknown inputs and with future sensor noise, and it holds similarly U + Y -T = 0 and V + Y -T = 0. These properties are required for the development of the subspace methods in the next section.

Second step: computation of subspace matrix H

Each subspace algorithm is related to a subspace matrix H with factorization property H = Γ p Z, where Z depends on the algorithm. Subspace algorithms may be divided into two families: covariance-driven and data-driven.

Covariance-driven subspace matrices are products of auto-or cross-covariance matrices. Data-driven matrices can be described as a product of covariance-driven matrices with a data Hankel matrix [START_REF] Benveniste | Nonstationary consistency of subspace methods[END_REF]. The main difference between both kinds of algorithms is the size of the associated subspace matrix. When the number of samples increases, covariancedriven matrices converge to their expected value, while their column dimension does not increase. Concerning datadriven matrices, their covariance components converge as the number of samples increases, but also their column dimension increases due to the increasing dimension of the data-matrix part. In both cases, the column space of the respective subspace matrices converges to the subspace defined by the observability matrix Γ p .

Orthogonal and oblique projections

Orthogonal and oblique projections of the data Hankel matrices are essential in the definition of subspace methods [START_REF] Van Overschee | Subspace Identification for Linear Systems: Theory, Implementation, Applications[END_REF][START_REF] Benveniste | Nonstationary consistency of subspace methods[END_REF][START_REF] Gandino | Covariance-driven subspace identification: A complete input-output approach[END_REF]. In the following, they are defined and some of their properties are recalled, before they are used in the subsequent sections.

Definition 3. Let M, P and Q be three matrices with the same number of columns. The orthogonal projection M/P of the row space of M onto the row space of P and the oblique projection M/ Q P of the row space of M on the row space of P along the row space of Q are defined by

M/P def = Mπ P def = MP T PP T † P, M/ Q P def = Mπ Q ⊥ Pπ Q ⊥ † P,
where π P is the orthogonal projector on P, Q ⊥ is the orthogonal complement of Q and † denotes the Moore-Penrose pseudoinverse.

Recall the properties of an orthogonal projector: It holds π P ⊥ = I -π P , π T P = π P , π P π P = π P and Pπ P ⊥ = 0. These properties are always exact and no convergence assumption is required [START_REF] Golub | Matrix computations[END_REF].

These geometric projections have a stochastic interpretation, since their definition coincides with the definition of the projection of Gaussian random vectors. Indeed, the projection of such vectors is defined through the conditional expectation. For the particular case of centered Gaussian random vectors X, Y, and Z it holds

Y/X def = E(Y/X) = Yπ X def = Cov(Y, X)Cov(X) † X, Y/ Z X def = Yπ Z ⊥ (Xπ Z ⊥ ) † X. (9)
Finally, let A and B be placeholders for any of the data Hankel matrices U -, U + , Y -or Y + (see Eq. ( 4)) that are filled with the discrete-time samples of the known input or output data, or any of the Hankel matrices in Section 2 related to the unknown signals. Then, AB T tends to a matrix filled with auto-or cross-covariances as in Eqs. ( 7)-( 8) and the orthogonal projection A/B = AB T (BB T ) † B relates to the conditional expectation in Eq. ( 9) when the number of data samples goes to infinity. In particular, if A and B are uncorrelated and covariances have converged, then AB T = 0 and A/B = 0.

Output-only algorithms

These algorithms are based on the assumptions that the unknown white noise input ũ(t) excites the system, and no known inputs are available (u(t) = 0). The observations may be corrupted by an unknown white sensor noise v(t).

Then Eq. ( 5) becomes:

                 Y -= Γ q X -+ Hq U -+ V - Y + = Γ p X + + Hp U + + V + X + = A q X -+ ∆q U -

3.2.2.a. OOcov: Output-only covariance-driven

This algorithm [START_REF] Benveniste | Nonstationary consistency of subspace methods[END_REF][START_REF] Peeters | Reference-based stochastic subspace identification for output-only modal analysis[END_REF][START_REF] Benveniste | Single sample modal identification of a nonstationary stochastic process[END_REF] is based on the auto-covariances between the N o outputs as in Eq. ( 7), which can be computed as a product of the future and past output data matrices,

H OOcov def = Y + Y -T (10) = (Γ p X + + Hp U + + V + )Y -T = Γ p X + Y -T = Γ p Z OOcov 3.2.2.b. OOort: Output-only orthogonal projection data-driven
It is also known as UPC algorithm [START_REF] Van Overschee | Subspace Identification for Linear Systems: Theory, Implementation, Applications[END_REF][START_REF] Benveniste | Nonstationary consistency of subspace methods[END_REF][START_REF] Peeters | Reference-based stochastic subspace identification for output-only modal analysis[END_REF], and is based on the projection of the row space of the future output data onto the row space of the past output data,

H OOort def = Y + /Y -= Y + Y -T Y -Y -T † Y - (11)
= Γ p Z OOcov Y -Y -T † Y -= Γ p Z OOort

Input/output algorithms

These algorithms are based on the assumption that the system is excited by known inputs u(t) and unknown white noise inputs ũ(t), and measurements are corrupted by white sensor noise v(t). It holds

                 Y -= Γ q X -+ H q U -+ Hq U -+ V - Y + = Γ p X + + H p U + + Hp U + + V + X + = A q X -+ ∆ q U -+ ∆q U - 3.2.3.a. IOcov: Input/output covariance-driven
This algorithm is based on the data-driven projection algorithm [2, Sec. 2.3.2] right multiplied by the past output data matrix Y -T in order to obtain a covariance-driven algorithm [START_REF] Ashari | Input-output subspace-based fault detection[END_REF],

H IOcov def = (Y + /U +⊥ )Y -T (12) = (Γ p X + + H p U + + Hp U + + V + )/U +⊥ Y -T = Γ p X + /U +⊥ + Hp U + + V + Y -T = Γ p X + /U +⊥ Y -T = Γ p Z IOcov ,
using the assumptions that known inputs are uncorrelated with unknown inputs and with sensor noise, and that past outputs are uncorrelated with future unknown inputs and with future sensor noise (see Remark 2).

3.2.3.b. IOort: Input/output orthogonal projection data-driven

This algorithm is described as the projection algorithm in [2, Sec. 2.3.2] for the "no noise case" when there are no unknown inputs and no sensor noise,

H IOort def = Y + /U +⊥ (13) = (Γ p X + + H p U + + Hp U + + V + )/U +⊥ = Γ p X + /U +⊥ + Hp U + + V + Indeed,
H IOobl def = Y + / U + W -= Y + /U + ⊥ W -/U + ⊥ † W - (14)
= Γ p Z IOobl where W -def =           U - Y -           and Z IOobl is defined in [2].
3.2.4. Some remarks about the considered algorithms

3.2.4.a. H OOcov and H IOcov

Assuming U + does not depend on the past, Z OOcov and Z IOcov are both product of future states and past observations. Both algorithms differ by the knowledge of a component of the excitation. For the output-only covariance-driven algorithm the whole excitation is considered to be unknown. The unknown input term in Eq. ( 10) is canceled by the product with Y -T as future unknown inputs are not correlated to past outputs (see Remark 2). The product converges to zero when the number of samples tends to infinity. The difference between Eq. ( 10) and Eq. ( 12) is how the input terms are canceled. In both Eq. ( 10) and Eq. ( 12), canceling U + is based on cross-covariance convergence. The known part U + in Eq. ( 12) is canceled from the orthogonal projection as U + /U +⊥ is zero. The cancellation through the orthogonal projection is satisfied all the time, whereas cancellation by the cross-covariance product is satisfied only when the number of samples tends to infinity.

H IOcov = H IOort Y -T ,
where the multiplication by Y -T removes the effect of both unknown inputs and sensor noise assuming convergence for the cross-covariance estimates (see Remark 2).

Third step: estimation of the modal parameters from H

Assume H is computed according to some subspace algorithm of interest. The final step is to retrieve modal parameters from H. First notice that all presented subspace methods have been defined such that H factorizes as H = Γ p Z at convergence. The objective is then to obtain the observability matrix

Γ p from H, since Γ p = [C T . . . (CA p-1
) T ] T contains the system matrices A and C. Thanks to the factorization property, H and Γ p share the same column space [START_REF] Van Overschee | Subspace Identification for Linear Systems: Theory, Implementation, Applications[END_REF]. Thus, Γ p can be obtained by the singular value decomposition (SVD) [START_REF] Golub | Matrix computations[END_REF]

H = US V T = U 1 U 2           S 1 0 0 S 2                     V T 1 V T 2           (15)
as

Γ p = U 1 S
from the shift invariance property of Γ p as the least squares solution A = Γ † up Γ dw [START_REF] Peeters | Reference-based stochastic subspace identification for output-only modal analysis[END_REF], where Γ dw and Γ up are obtained from Γ p by removing the first and last block row, respectively, as

Γ dw =                   CA . . . CA p-1                   , Γ up =                   C . . . CA p-2                   .
Finally, the modal parameters are obtained from matrices A and C as stated in Section 2. Hence, the computation of the modal parameters is straightforward when the subspace matrix H is computed.

From data-driven to covariance-driven subspace algorithms

In contrast to the covariance-driven algorithms, the number of columns of matrix H depends on the number of data samples N for the data-driven algorithms. Hence, the matrix H does not converge to a fixed limit for data-driven algorithms for N → ∞, which is a problem for the subsequent covariance analysis. In this section, we propose a simple way to transform data-driven algorithms to covariance-driven algorithms for the sake of covariance computation, where the resulting covariance-driven algorithm yields identical estimates of the observability matrix Γ p as the original data-driven algorithm. Hence, the variance of the resulting modal parameter estimates is also identical between both algorithms. Note that only Γ p is required for the estimation of the modal parameters (frequencies, damping ratios and mode shapes) from the system matrices (A, C), and the identification of (B, D) and the modal mass is beyond the scope of this paper.

Let H dat be the subspace matrix of a data-driven algorithm. Then, the corresponding covariance-driven algorithm is defined through the squared subspace matrix

H cov = H dat H T dat
Notice that this was already exploited in [START_REF] Benveniste | Nonstationary consistency of subspace methods[END_REF] to relate the consistency of both data-driven and covariance-driven algorithms together. It is easy to see that both matrices H cov and H dat indeed have the same column space: let the thin SVD of H dat = US V T be given, then H cov = H dat H T dat = US2 U T , and hence both H cov and H dat have the same left singular vectors (up to a change of basis). Thus, subspace-based system identification using H cov or H dat leads to identical results (see Section 3.3), and the covariance computation for data-driven algorithms can be based on H cov .

Based on the data-driven algorithms described in Sections 3.2.2.b, 3.2.3.b, 3.2.3.c, the respective covariance-driven algorithms can be formulated as follows.

H OOort 2 = Y + π Y -Y + T = H OOort Y + T = Y + Y -T Y -Y -T † Y + Y -T T . ( 16
H IOort 2 = Y + π U + ⊥ Y + T = H IOort Y + T . ( 17
)
Note that this right multiplication with Y + T does not remove the effect of unknown inputs and sensor noise, whereas this would be the case with Y - = H IOobl H T IOobl . Using the pseudoinverse property A † = A T (AA T ) † for an arbitrary matrix A, the pseudoinverse in Eq. ( 14) can be replaced, leading to

H IOobl = Y + /U + ⊥ W -/U + ⊥ T W -/U + ⊥ W -/U + ⊥ T † W -.
Squaring this matrix and using projector properties leads finally to the covariance-driven matrix

H IOobl 2 = Y + π U + ⊥ W -T W -π U + ⊥ W -T † W -W -T W -π U + ⊥ W -T † W -π U + ⊥ Y + T . (18)

Variance estimation

In this section, the notations and basic principles of the covariance computation are introduced based on [START_REF] Reynders | Uncertainty bounds on modal parameters obtained from stochastic subspace identification[END_REF][START_REF] Döhler | Efficient multi-order uncertainty computation for stochastic subspace identification[END_REF].

Until now this method was developed only for the output-only covariance-driven subspace algorithm. Here, the variance estimation algorithms are extended to the considered data-driven and input/output subspace algorithms by developing expressions for the covariance of the corresponding subspace matrices.

Theory and definitions

The uncertainty of the modal parameter estimates comes from the fact that only output or input/output sequences of finite length N are known. Expressions for their variance based on an input/output data sequence of length N are hence also estimates, and the central limit theorem (CLT) ensures that these expressions are asymptotically correct (i.e. for N → ∞) [9, Ch. 9].

The computation of the modal parameter covariance results from the propagation of sample covariances on autoor cross-covariance estimates (involving measured outputs, and for the case of input/output algorithms also the known inputs) through all steps of the modal identification algorithm. These sample covariances reflect in particular the unknown inputs due to non-measurable excitation sources and the sensor noise, and they contribute in a non-trivial way to the covariance of the modal parameter estimates. The propagation to the modal parameter estimates is based on the delta method [START_REF] Casella | Statistical inference[END_REF], where the analytical sensitivity matrices are obtained from perturbation theory [START_REF] Pintelon | Uncertainty calculation in (operational) modal analysis[END_REF][START_REF] Reynders | Uncertainty bounds on modal parameters obtained from stochastic subspace identification[END_REF].

Definition 4. Let θ(0) be the expected value of a parameter θ and θ(ε) its estimate, where ε denotes a small perturbation. A first order Taylor approximation of θ(ε) is given by θ(ε) ≈ θ(0)

+ ∂θ ∂ε ε=0 ε. Define ∆θ def = θ(ε) -θ(0) ≈ ∂θ ∂ε ε=0 ε.
Lemma 5. Let g be a continuous and differentiable function of θ. A first order Taylor approximation of g

(θ(ε)) is given by g(θ(ε)) ≈ g(θ(0)) + ∂g(θ) ∂θ θ=θ(ε) (θ(ε) -θ(0)), then ∆g(θ) = g(θ(ε)) -g(θ(0)) ≈ ∂g(θ) ∂θ θ=θ(ε) (θ(ε) -θ(0)). Then ∆g(θ) ≈ ∂g(θ) ∂θ θ=θ(ε) ∆θ,
and the covariance yields

Cov (g(θ)) = E ∆g(θ)∆g(θ) T ≈ ∂g(θ) ∂θ θ=θ(ε) Cov(θ) ∂g(θ) ∂θ T θ=θ(ε) . (19)
A sample covariance is easily computed from the measured data (see Section 5.2). With Eq. ( 19), it can be propagated through the whole computation chain of the subspace algorithms (see Sections 5.3 and 5.4). Since the auto-and cross-covariances of the data are asymptotically Gaussian (for N → ∞) [START_REF] Hannan | Multiple time series[END_REF], the delta method ensures that each variable in this computation chain -down to the modal parameters -is also asymptotically Gaussian [START_REF] Casella | Statistical inference[END_REF]. Hence the respective covariance estimate in Eq. ( 19) is consistent.

The challenge is to obtain the analytical sensitivity matrices for each step of the algorithm by perturbation theory that are required in Eq. [START_REF] Mevel | Input/output versus output-only data processing for structural identificationapplication to in-flight data analysis[END_REF]. To illustrate with an example, let M = PQ be a matrix product of P ∈ R k×l and Q ∈ R l×m .

The perturbation theory consists then on relating the covariance of M to the covariances of P and Q. Let vec(.) be the vectorization operator which stacks the columns of a matrix. With ∆M = ∆PQ + P∆Q, it follows using the vec properties vec (∆M) = (Q T ⊗ I k )vec (∆P) + (I m ⊗ P)vec (∆Q), and finally

vec (∆M) = J           vec (∆P) vec (∆Q)           where J = (Q T ⊗ I k ) (I m ⊗ P) , Cov(vec(M)) = E vec (∆M) vec (∆M) T = J Cov                     vec(P) vec(Q)                     J T .
In the following, some further properties of the involved matrix operations are needed. Define the permutation

matrix P a,b def = a k=1 b l=1 E a,b k,l ⊗ E b,a l,k
, where E a,b k,l ∈ R a×b are matrices which are equal to 1 at entry (k, l) and zero elsewhere, and ⊗ denotes the Kronecker product [START_REF] Brewer | Kronecker products and matrix calculus in system theory[END_REF]. This matrix satisfies P a,b vec(G) = vec(G T ), for G ∈ R a×b .

Finally, we need to derive the sensitivities of components containing the pseudo-inverse in the projection terms within the definition of the subspace methods. While the relation ∆M -1 = -M -1 ∆MM -1 holds for invertible matrices M, a similar relationship for the pseudo-inverse requires some stronger assumptions and is established for our case in the following lemma.

∆(R a R † b R c ) = ∆R a R † b R c -R a R † b ∆R b R † b R c + R a R † b ∆R c .
Proof. See Appendix C.

First step: covariance estimation of the auto-and cross-covariance matrices

The starting point of the variance propagation to the modal parameters is the covariance of the auto-and crosscovariance matrices that are involved in a chosen subspace algorithm. For the considered algorithms in this paper, the auto-and cross-covariance matrices are the following:

R 1 def = Y + Y + T ∈ R pN o ×pN o , R 2 def = Y + Y -T ∈ R pN o ×qN o , R 3 def = Y -Y -T ∈ R qN o ×qN o R 4 def = Y + U + T ∈ R pN o ×pN i , R 5 def = Y -U + T ∈ R qN o ×pN i , R 6 def = Y + U -T ∈ R pN o ×qN i , R 7 def = Y -U -T ∈ R qN o ×qN i , R 8 def = U + U + T ∈ R pN i ×pN i , R 9 def = Y + W -T ∈ R pN o ×q(N o +N i ) , R 10 def = W -U + T ∈ R pN i ×q(N o +N i ) , R 11 def = W -W -T ∈ R q(N o +N i )×q(N o +N i )
Notice that the convergence of these matrices depends on the data length. The central limit theorem ensures that they are asymptotically Gaussian and afflicted with uncertainty for a finite data length [START_REF] Hannan | Multiple time series[END_REF].

Hence, the computation of Cov(vec(R i ), vec(R j )) between any of the matrices R i and R j for the considered subspace algorithm (including i = j) is required. This computation follows the lines of the OOcov algorithm, i.e. for Cov(vec(R 2 )), as described in detail in [START_REF] Reynders | Uncertainty bounds on modal parameters obtained from stochastic subspace identification[END_REF][START_REF] Döhler | Efficient multi-order uncertainty computation for stochastic subspace identification[END_REF], and is generalized as follows.

Let A, B, C and D be placeholders for any of the data Hankel matrices U -, U + , Y -or Y + (see Eq. ( 4)), such that R i = AB T and R j = CD T are the respective auto-or cross-covariance matrices. Divide the data matrices into n b blocks and normalize them with respect to their length, such that

√ N A = √ N b A 1 A 2 . . . A n b , √ N B = √ N b B 1 B 2 . . . B n b , √ N C = √ N b C 1 C 2 . . . C n b , √ N D = √ N b D 1 D 2 . . . D n b ,
where each block A k , B k , C k and D k may have the same length N b , with n b • N b = N for simplicity. Each block may be long enough to assume statistical independence between the blocks. On each of these blocks, the respective auto-or cross-covariance estimate can be computed as R k i = A k B kT and R k j = C k D kT , which can be assumed to be i.i.d., yielding

R i = 1 n b n b k=1 R k i , R j = 1 n b n b k=1 R k j . It follows Cov (vec(R * )) = 1 n b Cov vec(R k * )
, and the covariance between the auto-or cross-covariance matrices can be computed from the usual sample covariance as

Cov(vec(R i ), vec(R j )) = 1 n b (n b -1) n b k=1 vec(R k i ) -vec(R i ) vec(R k j ) -vec(R j ) T . ( 20
)
With this expression, covariance estimates for covariance-driven subspace identification are achievable. It is not directly applicable for data-driven algorithms, where the size of the matrices changes with the number of available samples, requiring the squaring of the subspace matrix as shown in Section 4. Further subspace algorithms besides the ones investigated in the current paper can fit to this approach as long as the corresponding subspace matrix H can be expressed with respect to similar auto-or cross-covariance matrices as R i , i = 1, . . . , 11.

Second step: sensitivity of H for each subspace algorithm

In this section, the sensitivities of the subspace matrices H are derived with respect to the underlying auto-and cross-covariance matrices in their computation, depending on each considered subspace algorithm. With the sample covariance of the auto-and cross-covariance matrices from the previous section, the covariance of the subspace matrix may then be computed based on Eq. ( 19), and propagated to the system matrices A and C and the modal parameters in the next step. From its definition in Eq. ( 16), it follows

H OOort 2 = R 2 R † 3 R T 2
and thus, together with Lemma 6, From its definition in Eq. ( 12) it follows

∆H OOort 2 = ∆R 2 R † 3 R T 2 -R 2 R † 3 ∆R 3 R † 3 R T 2 + R 2 R † 3 ∆R T 2 , vec ∆H OOort 2 = J H OOort 2 ,R           vec (∆R 2 ) vec (∆R 3 )           , where J H OOort 2 ,R = R 2 R † 3 ⊗ I pN o + (I pN o ⊗ R 2 R † 3 )P pN o ,qN o -R 2 R † 3 ⊗ R 2 R † 3 .

Input/output algorithms

H IOcov = Y + Y -T -Y + U + T U + U + T † U + Y -T = R 2 -R 4 R † 8 R T 5
, and thus From its definition in Eq. ( 17), it follows

∆H IOcov = ∆R 2 -∆R 4 R † 8 R T 5 + R 4 R † 8 ∆R 8 R † 8 R T 5 -R 4 R † 8 ∆R T 5 , vec(∆H IOcov ) = J H IOcov ,R                          vec(∆R 2 ) vec(∆R 4 ) vec(∆R 8 ) vec(∆R 5 )                          , where J H IOcov ,R = I pN o qN o -R 5 R † 8 ⊗ I pN o R 5 R † 8 ⊗ R 4 R † 8 -(I qN o ⊗ R 4 R † 8 )P qN o ,pN i .
H IOort 2 = Y + Y + T -Y + U + T U + U + T † U + Y + T = R 1 -R 4 R † 8 R T 4

and thus

∆H IOort 2 = ∆R 1 -∆R 4 R † 8 R T 4 + R 4 R † 8 ∆R 8 R † 8 R T 4 -R 4 R † 8 ∆R T 4 , vec(∆H IOort 2 ) = J H IOort 2 ,R                  vec(∆R 1 ) vec(∆R 4 ) vec(∆R 8 )                  , where J H IOort 2 ,R = I pN o pN o -R 4 R † 8 ⊗ I pN o -(I pN o ⊗ R 4 R † 8 )P pN o ,pN i R 4 R † 8 ⊗ R 4 R † 8 .
5.3.2.c. IOobl 2 : Input/output oblique projection data-driven squared

From its definition in Eq. ( 18), it follows

H IOobl 2 = Ω 1 Ω † 2 R 11 Ω † 2 Ω T 1 ,
where

Ω 1 = Y + π U + ⊥ W -T = R 9 -R 4 R † 8 R T 10 ∈ R pN o ×q(N o +N i ) and Ω 2 = W -π U + ⊥ W -T = R 11 -R 10 R † 8 R T 10 ∈ R q(N o +N i )×q(N o +N i ) .
After perturbation and vectorization of

H IOobl 2 , vec(∆H IOobl 2 ) = J H IOobl 2 ,R                                  vec (∆R 9 ) vec (∆R 4 ) vec (∆R 8 ) vec (∆R 10 ) vec (∆R 11 )                                  , (21)
where

J H IOobl 2 ,R = J (1) H,R J (2)
H,R J

H,R J

H,R J (5) H,R is given in Appendix D.

Third step: covariance of H

In the previous section, the sensitivity of the subspace matrix was related to the sensitivities of the underlying autoand cross-covariance matrices for each considered subspace algorithm, and the covariance Cov(vec(R i ), vec(R j )) between the auto-and cross-covariance matrices was obtained in Section 5.2. To finally compute the covariance of the subspace matrix, Eq. ( 19) is applied. For example, the OOcov algorithm requires the computation of Cov(vec(H OOcov )) =

Cov(vec(R 2 )) (see Section 5.3.1.a), or the OOort 2 algorithm yields Cov(vec

H OOort 2 ) = E            J H OOort 2 ,R           vec (∆R 2 ) vec (∆R 3 )                     vec (∆R 2 ) vec (∆R 3 )           T J T H OOort 2 ,R            = J H OOort 2 ,R           Cov(vec(R 2 )) Cov(vec(R 2 ), vec(R 3 )) Cov(vec(R 2 ), vec(R 3 )) T Cov(vec(R 3 ))           J T H OOort 2 ,R , (22)
requiring thus the computation of Cov(vec(R 2 )), Cov(vec(R 3 )) and Cov(vec(R 2 ), vec(R 3 )) (see Section 5.3.1.b). With the covariance of the subspace matrix, the covariance of the modal parameters is obtained as described in Appendix B.

Algorithmic summary

In this section, a summary of the variance computation for the modal parameters from the subspace methods of this paper is given. Starting point of the computation is the subspace matrix H computed on the output dataset or input/output dataset for the chosen subspace method and the resulting modal parameters (see Section 3). Then the following steps are taken:

1. Covariance of auto-and cross-covariance matrices: For the chosen subspace method, H depends on a number of auto-or cross-covariance matrices R i (see Sections 5.3.1 and 5.3.2). The covariance of and between all involved auto-or cross-covariance matrices is computed on n b data blocks of the available output or input/output data as shown in Eq. (20).

Covariance of subspace matrix:

The sensitivity matrix J H,R for the chosen subspace method is computed in Sections 5.3.1 and 5.3.2. Then, the covariance of the subspace matrix Cov(vec (H)) is computed as in Eq. ( 22).

Covariance of modal parameters:

The covariance of the identified modal parameters is then computed as outlined in Appendix B. A detailed computation is described in [START_REF] Reynders | Uncertainty bounds on modal parameters obtained from stochastic subspace identification[END_REF][START_REF] Döhler | Efficient multi-order uncertainty computation for stochastic subspace identification[END_REF].

Note that an efficient numerical implementation of these steps can be formulated based on [START_REF] Döhler | Efficient multi-order uncertainty computation for stochastic subspace identification[END_REF][START_REF] Döhler | Multi-order covariance computation for estimates in stochastic subspace identification using QR decompositions[END_REF].

Numerical validation: mass-spring chain

A simulation study was made on a damped mass-spring chain model of six degrees of freedom (see Figure 1)

to validate the variance computation for the different subspace methods in an academic example. The mass of all elements is m i = 1/20, stiffnesses are k 1 = k 3 = k 5 = 100, k 2 = k 4 = k 6 = 200 and damping is defined such that all modes have a damping ratio of ξ i = 0.020.

Monte Carlo simulations including unknown excitation and sensor noise

Gaussian white noise excitation force at element 1 is simulated, which is considered as the known input for the input/output algorithms, plus additional (unknown) Gaussian white noise excitation force at all six elements with 10% of the standard deviation of the known input. From these inputs, acceleration output time series of length N = 10,000 are simulated with time step τ = 0.02 s at elements 1, 2 and 5, where Gaussian white sensor noise with 5% of the standard deviation of each output is added, respectively.

These simulations are repeated 10,000 times. Each time the modal parameters and their variances are estimated with the five output-only and input/output subspace algorithms from this paper. The parameters used for the identification and variance computation in the subspace algorithms are:

• system order n = 12 (see Section 2 and Eq. ( 15)),

• future and past time horizon p = q = 16 (see Eq. ( 4)),

• number of data blocks for variance computation n b = 50 (see Section 5.2).

The means of the identified modal parameters with the different subspace algorithms are shown together with the theoretical modal parameters from the mass-spring chain model in Table 1. The identified values are very close to the theoretical ones, except for the IOort 2 method that is not adapted to noisy data (see Section 3.2.3.b), and except for the first mode in the OOcov method that would have required longer datasets.

In each simulation, the standard deviations of the modal parameters are computed for the different subspace algorithms with the methods in Section 5. They are obtained from the square root of the respective modal parameter point of view. Note that the estimates of the standard deviations show good agreement even for the cases where the identification results were not satisfying (in particular for the first mode in the OOcov method and all modes in the IOort 2 method). Note also that the standard deviations are significantly lower with the input/output methods, compared to the output-only methods, roughly by factor 4 in this application.

Table 2 also shows the one sigma sample standard deviation of the estimated standard deviations during the Monte Carlo simulations. It is mostly smaller than 15% of their mean value (except for the badly estimated modes), both for frequencies and damping ratios. This means that the estimation error of the standard deviation on a single dataset is in most cases less than 15% with the proposed methods, which is very reasonable for a second-order statistics. As examples, the histograms of the identified frequencies and damping ratios of the first mode, and the histograms of their estimated standard deviations with the methods from this paper are shown in Figure 2 for the Monte Carlo simulations with the IOcov method. The histograms validate empirically the Gaussian assumption for the distribution of the frequency and damping estimates, and they show the sample distribution of the their estimated standard deviations that also seems to be approximately Gaussian. The frequencies seem to fit better to a Gaussian distribution than the damping ratios, which is expected considering frequencies are easier to estimate [START_REF] Gersch | On the achievable accuracy of structural system parameter estimates[END_REF]. Still, the distribution for the damping ratios is satisfactory and well centered.

Identification without unknown inputs and without sensor noise

For the particular case where all the inputs are known (and not corrupted by noise) and the measured outputs are not corrupted by sensor noise, the modal identification with input/output algorithms should be perfect if the true model order is known, and if (even short) data of some minimum length is available. The identified modal parameters should match exactly the theoretical values in this case, and the uncertainty on their estimation should thus be zero. This is verified in this section.

For this purpose, Gaussian white noise excitation at element 1 of length N = 1,000 is simulated with time step τ = 0.02 s, which is considered as the known input of the system. From this input, one acceleration output dataset at elements 1, 2 and 5 is simulated without adding sensor noise. From this input/output dataset, modal identification with the estimation of their standard deviations is performed with the same parameters as in the previous section.

In Table 3, the errors between the identified modal parameters and the theoretical parameters from the model, as well as the estimated standard deviation σ f and σ ξ from the methods of this paper are shown. It can be seen that both the errors and the standard deviations are zero as expected, up to numerical errors, validating the uncertainty quantification with the input/output algorithms for this "no noise" case.

Note that in the present variance computation methods it is assumed that both the input and the output signals are (known) stochastic sequences. The standard deviations (σ f and σ ξ ) in Table 3 are correctly estimated at zero in this "no noise" case, which shows that assuming at the same time known and stochastic inputs is no contradiction in the derivation of the uncertainty quantification methods. On the contrary, the derivation of the terms related to the auto-or cross-covariances involving the inputs as it is done in the present methods is necessary to correctly obtain this result.

σ d f 6.3 • 10 -3 7.4 • 10 -3 6.0 • 10 -3 2.7 • 10 -3 4.8 • 10 -3 1.9 • 10 -2 IOort 2 | f -f model | 3.3 • 10 -
To show the effect of neglecting the stochastic nature of the known inputs in the derivation of the uncertainty quantification methods, let us assume that the terms related to uncertainties of the input covariances are zero, i.e. ∆R 8 = 0, and ∆(U -U +T ) = 0 and ∆(U -U -T ) = 0 in ∆R 10 and ∆R 11 (see definitions in Section 5.2). The corresponding standard deviations are denoted by σ d f and σ d ξ and shown in Table 3. Note that these values are distinctively different from zero in our example, which again supports the assumption on using the stochastic nature of the inputs in the derivation.

Numerical validation: aircraft model

In this section, the five considered subspace algorithms and the new methods for their variance computation are applied using a complex numerical model of an unnamed aircraft, simulating a ground vibration test. The performance of the methods is evaluated in an OMA(X) context. In particular, robustness to noise on sensors is tested.

Benchmark

The simulation of time series was performed on a linear time invariant system, from a finite element model of an unnamed aircraft with around 200,000 degrees of freedom using the Newmark numerical scheme. A constant damping ratio of 15 was chosen for each mode. The inputs of the system are simulated as white noise excitation triggered by control surfaces of the aircraft. The responses of the system are observed at 46 sensors (outputs of the system) along the aircraft. Simulations of a length of 128 s were made and sampled at 2048 Hz. Such a high sampling frequency was necessary to reach scheme convergence until 40 Hz. The 262,144 samples of each time signal are divided into ten consecutive blocks. Power spectral densities are estimated as described in Appendix A for frequencies from 0 Hz to 40 Hz. Finally, auto-and cross-covariances are estimated using 10 blocks. The simulations are performed in order to respect ground vibration test conditions.

Strategy for uncertainty validation

The purpose of this section is to validate the variance computation based on one realization of each input and output signals in a more realistic context than in Section 6. Similarly, if the variance computation procedure is valid, it will coincide with the variance computed by the much more computationally intensive Monte Carlo approach. Notice that the Monte Carlo approach is not feasible in an operational context. Notice also that the proposed sensitivity-based approach will only need one sequence of each input and output signals since they are ergodic signals. This section focuses on the uncertainty validation for frequencies and damping ratios.

To validate the developed algorithms, vibration simulations of N s = 200 datasets are performed. Each simulation provides responses of the system under white noise excitation triggered by control surfaces. The input signals are different realizations for all simulations but all have the same statistical properties. From these simulations, N s = 200 modal identifications are carried out. Since the system order is rarely known (never when using real-life measured data), the modal parameters are estimated at several orders in a stabilization diagram. Considering a specific order, each identification provides a set of estimated modes ( f , ξ and ϕ) and their estimated standard deviation (σ f , σ ξ and σ ϕ ). From these 200 sets of modal parameters, the Monte Carlo validation is performed. Note that the different estimates are i.i.d. random variables for each mode.

To perform the validation of the variance computation, we compare the sample covariance σ 2 MC of the 200 Monte Carlo estimates of a quantity with the mean of their computed variances σ 2 est using the new methods of this paper, similar as in Section 6.1. Then the following ratio is computed for each estimated quantity m (frequency or damping ratio) at each model order for validation, to compare both variances:

C m = min σ 2 MC,m , σ 2 est,m max σ 2 MC,m , σ 2 est,m ∈ [0, 1]
The criterion to validate the variance computation for a quantity is set at C m ≥ 0.7, i.e. the computed variance with the methods of this paper deviates in mean less than 30 % from the Monte Carlo sample variance, and the respective standard deviations (as the square root of the variance) deviate less than 16 %. Note that these thresholds are empirical.

Plotting frequencies (or damping ratios) for each order yields so-called frequency (or damping) stabilization diagrams [START_REF] Döhler | Fast multi-order computation of system matrices in subspace-based system identification[END_REF][START_REF] Peeters | Reference-based stochastic subspace identification for output-only modal analysis[END_REF]. The solid lines in the following plots correspond to the exact values from the model. Stable modes are modes that appear at multiple model orders unlike spurious modes. Frequencies are normalized by the highest displayed frequency value. Axes are adjusted from 0 to 1 for the frequency range and for orders 1 to 100. In our uncertainty validation framework, results are presented in particular stabilization diagrams corresponding to the Monte Carlo simulations: instead of plotting the identified values, their means obtained from the 200 identifications are plotted. To find the same mode in the different identifications, a criterion was used based on the frequency ratio and the MAC value with the theoretical values from the model. An empirical identification threshold is set up, which needs to be exceeded for mode pairing. If a mode is not found in each identification, i.e., it is always below the identification threshold, it is removed from the analysis. Empty circles correspond to values with C m < 0.7, for which the variance computation is not validated. Circles filled with stars correspond to values with C m ≥ 0.7, where the variance computation is validated.

Validation results

7.3.1. OOcov: Output-only covariance-driven -Figure 3 In terms of identification: Seven stable modes are identified, and one more mode seems to appear around order 90. Both frequencies and damping ratios stabilize well. Several spurious modes appear when model order becomes too high. Sensor noise has no significant effect on stable modes. This algorithm, in this specific case, is robust to noise on sensors. Note that the number of spurious modes seems to decrease when adding noise on sensors. However, recall that the averages of 200 identifications are shown in the figures. The number of spurious modes does not really decrease at each identification. However these modes become much more unstable comparing all identifications. The automated selection removes them, when they are not found in every identification.

In terms of uncertainty estimations: With or without noise on sensors, uncertainty estimations on spurious modes are rarely validated. It seems that the uncertainty estimation is more accurate when the identification is more stable. Noise on sensors reduces the number of validated uncertainty estimations. This algorithm, in this specific case, does not seem to be robust to noise on sensors in terms of uncertainty estimations for both frequencies and damping ratios. In terms of identification: In this specific case, the identification is not well stabilized and numerous spurious modes appear. In particular, the damping identification is quite biased. This may be due to ill-conditioning of Y -Y -T .

Noise on sensors seems to help with stabilization. This effect can be explained by a sensors decorrelation and then an improvement of the conditioning of Y -Y -T . As with the OOcov algorithm, noise on sensors reduces stability of spurious modes comparing all identifications and then reduces the number of spurious modes that are kept. Notice that even if this algorithm provides badly stabilized identifications, it is robust to and even improved by noise on sensors both for frequencies and damping ratios.

In terms of uncertainty estimations: As with the OOcov algorithm, the accuracy of the uncertainty estimations increases with the stability of the modes. However, numerous modes can be unstable with respect to model order but still validated in terms of uncertainty estimations. Concerning noise effects, adding noise on sensors improves uncertainty estimations, like for the identification. These results show that uncertainty estimations may be accurate even if identification is not stable with respect to the model order. In terms of identification: In this specific case, ten modes are identified, three more than in output-only algorithms, and modes are found "sooner" (i.e. for lower orders). This may be due to the use of input information which eases modal identification. Noise on sensors seems to have no impact at all on the identification. We will see that it is the only algorithm which is completely robust to noise on sensors.

In terms of uncertainty estimations: Except for some spurious modes, uncertainty estimations on the modal identification are almost all validated for both frequencies and damping ratios. As for identification, noise on sensors has no effect on the validity of uncertainty estimations. In this specific case, the IOcov algorithm gives much better results than the output-only algorithms. In terms of identification: In this specific case, the same ten modes are identified as with the IOcov algorithm, plus one additional mode (the eleventh physical mode of the system) which becomes stable around order 90. Two other modes (the seventh and thirteenth physical modes of the system) seem to appear around order 95. However, the frequency of the seventh mode is not well identified and what looks like the thirteenth mode might only be a spurious mode. Depending on the mode, identification using the IOort 2 algorithm is either equivalent or better than using the IOcov algorithm. However, when noise is applied on sensors, the number of identified modes decreases.

Identification becomes close to the performance of the OOcov algorithm but remains better in terms of stability with respect to model order, especially for damping ratio estimations. This algorithm is not as robust to noisy sensor conditions as the IOcov algorithm due its theoretical properties (see Section 3.2.3.b) but it remains accurate enough to identify well excited modes.

In terms of uncertainty estimations: This algorithm is not really accurate in terms of uncertainty estimations.

It barely validates half the identified modes. However, conclusions stated about previous algorithms are still valid.

Uncertainty estimation is more accurate on physical modes than on spurious modes. Noise on sensors seems to have no significant effect on uncertainty estimation. This algorithm can be considered robust to noise on sensors but in this specific case less modes are validated than using the IOcov algorithm.

7.3.5. IOobl 2 : Input/output oblique projection data-driven squared -Figure 7 In terms of identification: Identification results are similar to the IOort 2 algorithm in terms of both number and stability of identified modes. Even the sensor noise effect is equivalent for these two algorithms. It significantly reduces the number of identified modes. However, when noise is applied on sensors the IOobl 2 algorithm gives better identification results than the IOort 2 algorithm, especially for the fifth and sixth physical modes of the system.

In terms of uncertainty estimations: In this specific case, this algorithm is the worst proposed algorithm for uncertainty estimations. As with IOort 2 , less than half of the identified modes are validated in terms of uncertainty estimations. However, when noise is applied on sensors, it becomes one of the best algorithm to estimate uncertainty both for frequencies and damping ratios. Results become equivalent to results obtained using the IOcov algorithm, except that less (physical and spurious) modes are identified with IOobl 2 than with IOcov.

Global conclusions on the simulations

In this section five modal identification algorithms have been tested in terms of both identification and uncertainty estimations. Note that these conclusions should not be generalized and may be specific to the studied case.

In terms of identification: Based on Figure 3-7 it is clear that using input information eases modal identification.

Modes are found sooner (i.e. for lower orders) and are more numerous in the input/output framework than in the output-only framework. In the output-only framework, the OOcov algorithm provides a more stable identification in general and better damping ratios than the OOort 2 algorithm. This may be due to ill conditioning of Y -Y -T (see Eq. ( 11)). In the input/output framework, the three studied algorithms are equivalent. However, a slight distinction can be made concerning the seventh and the eleventh physical modes. Only the IOort 2 algorithm identifies both of them. Until order 100, the IOobl 2 algorithm identifies only the seventh one and the IOcov algorithm none of them. Concerning noise on sensors effect, it was shown that it strongly depends on the algorithm. Remind that the presented results are an average of stable mode estimates for different identifications. Based on that, noise has no effect on the IOcov algorithm and almost no effect on the OOcov algorithm. It reduces the number of modes and improves their stability when using OOort 2 , IOort 2 and IOobl 2 .

In terms of uncertainty estimations: Concerning the use of input information, it seems to have no clear effect on uncertainty estimations. Validations are not better in the input/output framework than in the output-only framework.

Both the OOort • OOort 2 , Figure 9: As with the OOcov algorithm, seven modes are clearly identified. However, they are much less stable with respect to model order. Moreover, this identification contains more spurious modes than the identification obtained with OOcov. This example shows how useful the uncertainty-based clean up may be.

• IOcov, Figure 10: It is the first presented algorithm using input information. On can observe that modes appear "sooner" (i.e. for lower orders) when input information is used. Here, modes are also more stable and have smaller uncertainty values. Only a few spurious modes are identified (only 7 points in the whole stabilization charts). Once again, as with previous algorithms, one can observe that uncertainty estimates are smaller for more stable modes. shows as well that the more stable a mode, the smaller its uncertainty estimation. • IOobl 2 , Figure 12: Identification is similar to the IOcov algorithm. In this specific case, these two algorithms provide the best identification and small uncertainty values.

Discussion

In this paper, variance estimates were developed for output-only and input/output system identification schemes generalizing [START_REF] Reynders | Uncertainty bounds on modal parameters obtained from stochastic subspace identification[END_REF][START_REF] Döhler | Efficient multi-order uncertainty computation for stochastic subspace identification[END_REF] previously developed for the OOcov (output-only covariance-driven) subspace algorithm. Based on this method, variance expressions have been developed for OOort 2 (output-only orthogonal projection data-driven squared), IOcov (input/output covariance-driven), IOort 2 (input/output orthogonal projection data-driven squared) and IOobl 2 (input/output oblique projection data-driven squared) algorithms. Providing also a computation scheme for the mode shapes, the paper focusses on the estimation of the modal frequencies and damping ratios, which are parameters of interest for operational modal analysis during in flight tests.

These five modal identification algorithms have been validated in terms of both identification and variance estimations. In this validation, their robustness to sensor noise has also been tested. Finally they have been applied on real measurements. Besides an academic test case on a mass-spring chain, simulations and real data were provided on an unnamed aircraft from industrial cases and results have been normalized to allow publication.

The validation of the modal identification algorithms and their variance computation were based on Monte Carlo simulations. Output-only and input/output algorithms have been compared. In the investigated examples, the number of identified modes increases and modes appear at lower orders when using input information. Furthermore, the variance decreases with input information when using the appropriate algorithms. The validation of the variance computation was equally comparable with or without using inputs. In the academic test case, the results of the variance computation were close to the respective variances from the Monte Carlo simulations for all algorithms, including the use of sensor noise. In the industrial test cases, the effect of sensor noise has been shown strongly dependent on the respective algorithms. For the OOcov algorithm a degradation of both identification and variance estimations could be observed with sensor noise. For the OOort 2 algorithm, it was a slight improvement of both identification in terms of stability with respect to model order and variance estimations. Sensor noise did not change validation results for the IOcov algorithm. In terms of identification, impact of noise was similar with IOort 2 and IOobl 2 algorithms. One could observe a reduction of the number of identified modes. Only the more stable ones were found again with sensor noise. However, in terms of uncertainty validation, a huge improvement was observed for the IOobl 2 algorithm.

After validation, the considered algorithms have been tested on real-life measured data obtained during in-flight vibration tests. It appeared that non-physical poles (i.e., spurious modes) have often much higher variance values than physical ones (real modes of the structure). This allows to use variance estimation to clean up stabilization charts and then to ease their interpretation during identification. Furthermore, the mode selection from the stabilization diagram could be supported by the actual computed uncertainty bounds. While different model orders indeed help to reliably identify the modes, the results show that the model order should not be too high in same cases. There, the results deteriorate slightly or the uncertainty estimation is less validated for high orders, which might be due to overfitting or worse conditioning of the involved matrices for (pseudo-)inverting in the subspace methods at high orders.

Concerning the efficiency of the algorithms, it seems that input information eases identification and reduces uncertainty values. However, if only outputs are available, output-only algorithms give relevant results. If the five algorithms had to be classified in terms of both identification and uncertainty estimations, their ranking could be as follows:

IOcov > IOobl 2 > OOcov > IOort 2 > OOort 2
This classification has to be taken carefully. It cannot be generalized and it could be specific to the test case.

Recall that the chosen approach for variance computation cannot be used directly on data-driven algorithms. Indeed, the column dimension of data-driven matrices for subspace identification increases with the number of samples.

Hence, they are not compatible with a sample covariance formulation and the perturbation-based approach from [START_REF] Pintelon | Uncertainty calculation in (operational) modal analysis[END_REF][START_REF] Reynders | Uncertainty bounds on modal parameters obtained from stochastic subspace identification[END_REF][START_REF] Döhler | Efficient multi-order uncertainty computation for stochastic subspace identification[END_REF] cannot be applied directly. Fortunately, as shown in the current paper, the singular value decomposition properties ensure that the proposed variance computation scheme for the squared subspace algorithms provides equivalent uncertainty bounds for the pair (A, C) obtained from the data-driven algorithms. Thus, frequency, damping ratio and mode shape estimates are exactly the same using either data-driven algorithms or their squared counterpart. This allows to use the proposed scheme on the squared subspace matrices to estimate variances of modal parameters from data-driven algorithms. Still, the actual implementation of the data-driven algorithms may differ from the direct computation of the subspace matrices and implies usually computation procedures such as QR decomposition [START_REF] Van Overschee | Subspace Identification for Linear Systems: Theory, Implementation, Applications[END_REF], which is a supplementary numerical step before the SVD to obtain the observability matrix. While such a step does not change the variance expressions thanks to the equivalence of the algorithms, it should be incorporated as in [START_REF] Döhler | Multi-order covariance computation for estimates in stochastic subspace identification using QR decompositions[END_REF]. Furthermore, the efficient numerical framework for the variance computation of [START_REF] Döhler | Efficient multi-order uncertainty computation for stochastic subspace identification[END_REF] should be applied to the new methods from this paper. With this framework, computation times for all the methods are less than a minute for a whole stabilization diagram containing model orders until 100 on a standard desktop computer and hardly differ between the different subspace methods.

Finally, we recall that the input/output methods in this paper have been derived under the premise that the known inputs are the realization of a stochastic process. The Monte Carlo simulations have shown that this assumption gives coherent variance estimates with the methods of this paper. The stochastic properties of the inputs contribute to the variance estimates of the modal parameters. They do not necessarily mean an increase of the variance estimates of the modal parameters. Taking them into account is necessary for a correct variance estimation. In particular, it has been demonstrated in Section 6.2 that the variance is correctly estimated at zero in the case of no unknown inputs and no sensor noise when taking the stochastic properties of the inputs into account. Neglecting the stochastic properties of the inputs has led to a variance estimate different from zero in this case.

Deterministic inputs are not considered in this paper. Under a few assumptions for the known inputs such as persistent excitation [START_REF] Van Overschee | Subspace Identification for Linear Systems: Theory, Implementation, Applications[END_REF], quasi-stationarity for deterministic signals such as in [START_REF] Viberg | Analysis of state space system identification methods based on instrumental variables and subspace fitting[END_REF], and under the assumption of both known deterministic input and unknown stochastic noise, the current approach would be applicable provided the uncertainty is derived with respect to the noise only, and as long as auto-and cross-covariances can be defined for the known input as in [START_REF] Ljung | System identification[END_REF]. The corresponding variance computation could be handled by the current approach as a special case, where the variance with respect to the functions of the deterministic and known inputs is taken as zero.

Appendix B. From covariance of H to the modal parameters

Let H be the subspace matrix corresponding to a chosen subspace algorithm, which is given together with an estimate of Cov(vec(H)). This covariance is propagated to the observability matrix Γ p , then to the system matrices A and C and finally to the modal parameters, as outlined as follows. Matrix Γ p is obtained from an SVD of H in Eq. ( 15) as Γ p = U 1 S 1/2 1 . The sensitivity of Γ p is obtained by perturbation as ∆Γ p = ∆U 1 S 1/2 1 + U 1 ∆S 1/2 1 and finally, using J Γ,H as detailed in [START_REF] Reynders | Uncertainty bounds on modal parameters obtained from stochastic subspace identification[END_REF][START_REF] Döhler | Efficient multi-order uncertainty computation for stochastic subspace identification[END_REF],

Cov(vec(Γ p )) = E vec(∆Γ p )vec(∆Γ p ) T = J Γ,H Cov(vec(H))J T

Γ,H

Then, A and C are obtained from Γ p , and perturbations on A an C are obtained by perturbations of the equations in Section 3.3. Then vec(∆A) = J A,Γ vec(∆Γ p ) and vec(∆C) = J C,Γ vec(∆Γ p ). The covariances of A and C write as

Σ A,C = Cov                     vec(A) vec(C)                     =           J A,Γ J C,Γ           J Γ,H Σ H J T Γ,H J T A,Γ J T

4. 1 .

 1 OOort 2 : Output-only orthogonal projection data-driven squared With the definition of the underlying UPC algorithm in Section 3.2.2.b, set H OOort 2 def = H OOort H T OOort . It follows from projector properties

Lemma 6 .

 6 Let A, B and C be placeholders for any of the data Hankel matrices U -, U + , Y -or Y + , and let R a = AB T , R b = BB T and R c = BC T be the respective auto-or cross-covariance matrices. Then,

5. 3 . 1 .

 31 Output-only algorithms 5.3.1.a. OOcov: Output-only covariance-driven Due to Eq. (10), it holds H OOcov = R 2 and thus vec(∆H OOcov ) = vec(∆R 2 ).5.3.1.b. OOort 2 : Output-only orthogonal projection data-driven

5. 3

 3 .2.b. IOort 2 : Input/output orthogonal projection data-driven squared

Figure 1 :

 1 Figure 1: Mass-spring chain (with modal damping), three sensors, for simulations.

Figure 2 :

 2 Figure 2: Histograms for first mode identified with IOcov from Monte Carlo simulations: estimated frequencies, estimated standard deviations on frequencies, estimated damping ratios and estimated standard deviations on damping ratios.

15 1 . 8 • 15 σ f 5 . 5 •d f 6 . 1 • 10 -2 9 . 8 • 10 -3 2 . 8 • 10 -2 1 . 2 • 10 -1 3 . 4 • 10 -2 1 . 5 • 10 - 1 14 σ f 1 . 4 • 13 σ ξ 2 . 8 • 13 σ ξ 2 . 7 • 13 σ ξ 9 . 3 •

 18155561109810281012103410151011414132813271393 10 -15 3.6 • 10 -15 1.8 • 10 -15 1.8 • 10 -15 3.6 • 10 -10 -16 4.6 • 10 -16 3.8 • 10 -16 6.3 • 10 -16 6.0 • 10 -16 8.2 • 10 -16 σ IOobl 2 | ff model | 1.4 • 10 -14 4.4 • 10 -15 0 7.1 • 10 -15 1.8 • 10 -15 1.4 • 10 -10 -15 4.9 • 10 -16 4.0 • 10 -16 2.1 • 10 -15 1.8 • 10 -15 5.5 • 10 -15 IOcov |ξ -ξ model | 1.9 • 10 -12 1.2 • 10 -12 3.0 • 10 -13 1.2 • 10 -12 1.1 • 10 -12 5.9 • 10 -10 -13 9.7 • 10 -14 8.6 • 10 -13 3.7 • 10 -13 5.7 • 10 -13 3.9 • 10 -13 model | 1.4 • 10 -12 5.1 • 10 -14 5.2 • 10 -13 2.6 • 10 -13 8.9 • 10 -14 6.3 • 10 -10 -13 6.3 • 10 -14 5.2 • 10 -14 1.2 • 10 -13 3.2 • 10 -14 1.5 • 10 -13 model | 5.9 • 10 -12 1.0 • 10 -12 9.3 • 10 -14 6.8 • 10 -13 5.5 • 10 -13 9.8 • 10 -10 -13 2.1 • 10 -13 6.9 • 10 -14 1.8 • 10 -13 6.9 • 10 -14 2.3 • 10 -13

Figure 3 :

 3 Figure 3: OOcov -Monte Carlo stabilization charts with variance validation values for orders 1 to 100

7. 3 . 2 .Figure 4 :

 324 Figure 4: OOort 2 -Monte Carlo stabilization charts with success test values for orders 1 to 100

7. 3 . 3 .Figure 5 :

 335 Figure 5: IOcov -Monte Carlo stabilization charts with success test values for orders 1 to 100

7. 3 . 4 .Figure 6 :

 346 Figure 6: IOort 2 -Monte Carlo stabilization charts with success test values for orders 1 to 100

Figure 7 :

 7 Figure 7: IOobl 2 -Monte Carlo stabilization charts with success test values for orders 1 to 100

8 .Figure 9 :Figure 10 :

 8910 figures, the original uncleaned stabilization diagrams are shown on the left, whereas the cleaned stabilization diagrams using the uncertainty-based threshold are shown. Thresholds on uncertainties have been set to 0.25 Hz for frequencies and 100 for damping ratios. Results are as follows:• OOcov, Figure8: Seven modes are clearly identified. An eighth one appears at order 58. There are few spurious modes. Uncertainty-based clean up is relevant. It removes only spurious or unstable modes. One can observe that more stable modes have a smaller uncertainty estimation. Moreover, uncertainties are coherent with operationally

Figure 12 :

 12 Figure 12: IOobl 2 : Input/output oblique projection data-driven squared -stabilization charts of normalized frequencies (left: full, right: cleaned)

  there is no factorization property H IOort = Γ p Z IOort in the presence of unknown inputs or sensor noise. Hence, the theory of this algorithm is only adapted to Experimental Modal Analysis framework. Note that the formulation of this algorithm is equivalent to the MOESP algorithm[START_REF] Van Overschee | Subspace Identification for Linear Systems: Theory, Implementation, Applications[END_REF] Sec. 4.3.2] if there are no unknown inputs and no sensor noise.

3.2.3.c. IOobl: Input/output oblique projection data-driven This data-driven algorithm is known as N4SID [2, Sec. 4.3.1] and is again adapted to the input/output framework with unknown inputs and sensor noise. It yields

  ) 4.2. IOort 2 : Input/output orthogonal projection data-driven squared With the definition of the underlying projection algorithm in Section 3.2.3.b, set H IOort 2

	def = H IOort H T IOort , and
	projector properties yield

Table 1 :

 1 Means for frequencies and damping ratios from modal identification in 10,000 Monte Carlo simulations.

				Frequency (Hz)				Damping ratio ( )
	Algorithm	1	2	3	4	5	6	1	2	3	4	5	6
	Actual model 1.930 5.618 8.682 14.49 15.85 17.01 20.0 20.0 20.0 20.0 20.0 20.0
	OOcov 1.937 5.618 8.683 14.49 15.85 17.01 29.4 20.1 20.1 20.0 20.0 20.0
	OOort 2 1.933 5.618 8.683 14.49 15.84 17.01 21.7 20.1 20.0 20.0 19.7 19.9
	IOcov 1.930 5.618 8.682 14.49 15.85 17.01 20.2 20.0 20.0 20.0 20.0 20.0
	IOort 2 2.016 5.623 8.684 14.52 15.85 17.00 36.1 19.0 19.2 18.0 18.3 17.1
	IOobl 2 1.930 5.618 8.682 14.49 15.85 17.01 20.1 20.0 20.0 20.0 20.0 20.0

Table 2 :

 2 Mean of estimated standard deviations (σ est ) and sample standard deviations (σ MC ) for frequencies and damping ratios of the six modes in 10,000 Monte Carlo simulations. The mean σ est of these 10,000 estimated standard deviations from each of the Monte Carlo simulations is shown for each modal parameter and algorithm in Table2. These values are compared to the sample standard deviations σ MC of the 10,000 estimated modal parameters in Table2. An excellent agreement between the estimated standard deviations with the methods from this paper and the sample statistics from the Monte Carlo simulations can be found for all subspace methods, validating the new methods for variance computation from a theoretical

					Standard deviation of frequency (Hz •10 -2 )
	Algorithm		1	2	3	4	5	6
	OOcov	σ est	2.02 ± 1.80 1.10 ± 0.15 1.55 ± 0.21 2.60 ± 0.36 2.64 ± 0.38 3.51 ± 0.48
		σ MC 1.99	1.14	1.59	2.63	2.65	3.48
	OOort 2	σ est	0.74 ± 0.12 1.05 ± 0.14 1.41 ± 0.18 2.45 ± 0.33 2.42 ± 0.33 3.01 ± 0.42
		σ MC 0.69	1.08	1.45	2.46	2.44	2.99
	IOcov	σ est	0.42 ± 0.06 0.24 ± 0.03 0.29 ± 0.03 0.66 ± 0.08 0.49 ± 0.05 0.87 ± 0.10
		σ MC 0.41	0.25	0.29	0.68	0.50	0.87
	IOort 2	σ est	2.00 ± 1.17 0.23 ± 0.03 0.27 ± 0.03 0.62 ± 0.07 0.45 ± 0.05 0.75 ± 0.08
		σ MC 2.74	0.24	0.27	0.63	0.46	0.75
	IOobl 2	σ est	0.34 ± 0.05 0.23 ± 0.03 0.27 ± 0.03 0.62 ± 0.07 0.45 ± 0.05 0.77 ± 0.09
		σ MC 0.33	0.23	0.27	0.64	0.46	0.78
					Standard deviation of damping ratio ( )
	Algorithm		1	2	3	4	5	6
	OOcov	σ est	16.46 ± 26.64 1.98 ± 0.28 1.80 ± 0.25 1.70 ± 0.22 1.67 ± 0.23 2.09 ± 0.28
		σ MC 17.12	2.07	1.85	1.71	1.71	2.09
	OOort 2	σ est	4.08 ± 0.81	1.83 ± 0.26 1.60 ± 0.21 1.61 ± 0.21 1.50 ± 0.20 1.81 ± 0.23
		σ MC	3.91	1.91	1.64	1.62	1.53	1.81
	IOcov	σ est	2.13 ± 0.34	0.43 ± 0.05 0.33 ± 0.04 0.45 ± 0.05 0.31 ± 0.04 0.52 ± 0.06
		σ MC	2.07	0.44	0.34	0.45	0.31	0.53
	IOort 2	σ est	10.46 ± 12.20 0.39 ± 0.05 0.30 ± 0.03 0.39 ± 0.04 0.27 ± 0.03 0.42 ± 0.05
		σ MC 15.69	0.40	0.30	0.39	0.27	0.42
	IOobl 2	σ est	1.72 ± 0.25	0.40 ± 0.05 0.31 ± 0.04 0.42 ± 0.05 0.28 ± 0.03 0.47 ± 0.05
		σ MC	1.65	0.41	0.31	0.43	0.29	0.48
	variances.

Table 3 :

 3 Identification error and estimated standard deviations on input/output data without noise in one simulation.

	Frequency (Hz)

model | 7.8 • 10 -16 3.6 • 10 -16 1.1 • 10 -15 1.8 • 10 -16 7.1 • 10 -16 1.1 • 10 -15 σ f 3.9 • 10 -16 6.8 • 10 -16 1.8 • 10 -15 6.1 • 10 -15 6.9 • 10 -15 3.2 • 10 -15

  2 and IOcov algorithms provide much more validated uncertainty estimations than OOcov, IOort 2 and IOobl 2 . The OOcov algorithm is the only studied algorithm which is not robust to noise on sensors. Noise has absolutely no influence on uncertainty estimated by the IOcov algorithm. It has no clear impact when using OOcov: Output-only covariance-driven -stabilization charts of normalized frequencies (left: full, right: cleaned) either OOort 2 or IOort 2 . Finally, noise on sensors makes IOobl 2 the best proposed algorithm in terms of uncertainty estimations.

	50											50
	45											45
	40											40
	25 30 35 model order											25 30 35 model order
	20											20
	15											15
	10											10
	5											5
	0.5	0.55	0.6	0.65	0.7	0.75	0.8	0.85	0.9	0.95	1	0.5	0.55	0.6	0.65	0.7	0.75	0.8	0.85	0.9	0.95	1
					normalized frequency								normalized frequency
			Figure 8:

1 , where the SVD is truncated at the desired model order n [2]. The output matrix C is then obtained by direct extraction of the first block row of Γ p as C = I N o 0 N o ,(p-1)N o Γ p . The state transition matrix A is obtained

Conclusion

Modal identification is commonly used in industry in the design process for a structure. The uncertainty information related to the identified modal parameters is at least as valuable as the modal parameters themselves in order to judge their quality and confidence. Moreover, in some highly dangerous and expensive vibration tests (e.g. flutter tests [START_REF] Pickrel | Flight flutter testing of transport aircraft: in-flight modal analysis[END_REF]), uncertainty estimations may be very helpful in terms of both safety and cost reduction. In this paper, the main results are new expressions for variance estimation for four output-only and input/output subspace-based modal identification algorithms. It has been shown that considering the stochastic properties of the known input leads to a correct formulation of the variance estimation scheme. Concerning the considered subspace algorithms, it has been shown how data-driven algorithms as found in [START_REF] Van Overschee | Subspace Identification for Linear Systems: Theory, Implementation, Applications[END_REF][START_REF] Ashari | Input-output subspace-based fault detection[END_REF] can be related to an equivalent covariance-driven form yielding identical estimates for the pair (A, C). From this remark, covariance estimates for (A, C) and subsequently for frequencies, damping ratios and mode shapes were derived. The validation of this approach has been carried out with Monte Carlo simulations, and it has been applied on real in-flight aircraft data. Further work should include the variance estimation for the pair (B, D) as well as for the modal mass.

Appendix A. Computation of auto-and cross-covariance functions

Let F be the discrete Fourier transform operator, and let a and b be placeholders for two discrete-time signals (e.g. u and y, or y and y) whose cross-covariance shall be computed. The signals are divided into n b segments a (k) and b (k) , k = 1, . . . , n b , of equal size (with or without overlap). An apodization window is often applied to these segments to reduce the leakage effect. Fourier transformations of a (k) and b (k) are given by a (k)

where J A,Γ and J C,Γ are detailed in [START_REF] Reynders | Uncertainty bounds on modal parameters obtained from stochastic subspace identification[END_REF][START_REF] Döhler | Efficient multi-order uncertainty computation for stochastic subspace identification[END_REF]. Covariances of modal parameters are finally as in [START_REF] Reynders | Uncertainty bounds on modal parameters obtained from stochastic subspace identification[END_REF][START_REF] Döhler | Efficient multi-order uncertainty computation for stochastic subspace identification[END_REF].

, and it remains to show the identity

to conclude the proof. Let the SVD of B be given,

where S 1 contains all non-zero singular values, and

Both sides of (C.1) are now developed. It holds To derive vec(∆H IOobl 2 ) in Eq. ( 21) with

and

we first show

First, consider the "noisy" case with the presence of non-zero unmeasured inputs (process noise) ũ and sensor noise v. For this case, it has been shown in [START_REF] Van Overschee | A unifying theorem for three subspace system identification algorithms[END_REF] that the term W -/U + ⊥ is of full row rank. Hence, Ω 2 is invertible and Eq. (D.1) holds since Ω † 2 = Ω -1 2 and due to the property ∆(Ω -1 2 ) = -Ω -1 2 ∆(Ω 2 )Ω -1 2 . Second, consider the "no noise" case with no unmeasured inputs and no sensor noise. For this case, it has been shown in [START_REF] Van Overschee | Subspace Identification for Linear Systems: Theory, Implementation, Applications[END_REF]Appendix A.1] that the term W -/U + ⊥ is indeed rank deficient, but that W -/U + ⊥ and W -have the same column space. Due to this property, Lemma 6 can be applied to show Eq. (D.1). Hence, Eq. (D.1) holds for both the "noisy" and "no noise" cases. Then, a subsequent relationship for

) is straightforward to show, and the relationship for vec(∆H IOobl 2 ) in Eq. ( 21) follows after vectorization and re-ordering, where