Subspace-Based Algorithms for Structural Identiﬁcation, Damage Detection, and Sensor Data Fusion

This paper reports on the theory and practice of covariance-driven output-only and input/output subspace-based identiﬁcation and detection algorithms. The motivating and investigated application domain is vibration-based structural analysis and health monitoring of mechanical, civil, and aeronautic structures.


Framework
Detecting and localizing damages for monitoring the integrity of structural and mechanical systems is a topic of growing interest, due to the aging of many engineering constructions and machines and to increased safety norms. Many current approaches still rely on visual inspections or local nondestructive evaluations performed manually, for example acoustic, ultrasonic, radiographic or eddy-current methods. These experimental approaches assume an a priori knowledge and the accessibility of a neighborhood of the damage location. Automatic global vibration-based monitoring techniques have been recognized to be useful alternatives to those local evaluations [1][2][3][4][5].
Many structures to be monitored (e.g., civil engineering structures subject to wind and earthquakes, aircrafts subject to turbulence) are subject to both fast and unmeasured variations in their environment and small slow variations in their modal (vibrating) properties. While any change in the excitation is meaningless, damages or fatigues on the structure are of interest. But the available measurements (e.g., from accelerometers) do not separate the effects of the external forces from the effect of the structure. Moreover the changes of interest (1% in eigenfrequencies) neither are visible on the signals nor on their spectra. A global health monitoring method must rather rely on a model which will help in discriminating between the two mixed causes of the changes that are contained in the data.
Most classical modal analysis and vibration monitoring methods basically process data registered either on test beds or under specific excitation or rotation speed conditions. However a need has been recognized for vibration monitoring algorithms devoted to the processing of data recorded inoperation, namely, during the actual functioning of the considered structure or machine, without artificial excitation, speeding down or stopping [6,7].
In this framework, covariance-driven input/output and output-only subspace-based algorithms have been developed for the purpose of structural identification, damage detection and diagnosis, and merging sensor data from multiple measurements setups registered at different periods of time. The purpose of this paper is to present an overview of the theory and practice of these algorithms.

Paper outline
The paper is organized as follows. In Section 2 we recall the main features of the output-only covariance-driven subspace-based identification algorithm. We exhibit a key factorization property and introduce the parameter estimating function associated with this algorithm. We elaborate further on the factorization property in Sections 3, 4, and 5 and on the estimating function in Sections 6 and 7.

EURASIP Journal on Advances in Signal Processing
We explain in Section 3 how the joint use of the factorization property and various projections helps handling both known inputs and unknown excitations; in Section 4 how the factorization property helps extending the covariance-driven subspace algorithm to the joint processing of data from multiple measurements setups recorded under nonstationary excitation; and in Section 5 the role of the factorization property in the analysis of the consistency of these identification algorithms under nonstationary excitation. Section 6 is devoted to a batch-wise change detection algorithm built on the covariance-driven subspace-based estimating function and the statistical local approach to the design of change/fault/damage detection algorithms. In Section 7, another asymptotic for this estimating function is used for designing a sample-wise recursive CUSUM detection algorithm.
Some typical application examples are discussed in Section 8. Finally comments on ongoing research and conclusions are drawn in Section 9.

OUTPUT-ONLY EIGENSTRUCTURE IDENTIFICATION
It is well established [8][9][10][11] that the vibration-based structural analysis and health monitoring problems translate into the identification and monitoring of the eigenstructure of the state transition matrix F of a linear dynamical statespace system excited by a zero-mean Gaussian white noise sequence (V k ): namely, the (λ, ϕ λ ) defined by The associated parameter vector is where Λ is the vector whose elements are the eigenvalues λ, Φ is the matrix whose columns are the ϕ λ 's, and vec is the column stacking operator. This parameter is canonical, that is invariant with respect to a change in the state space basis. Subspace-based methods are the generic name for linear systems identification algorithms based on either time domain measurements or output covariance matrices, in which different subspaces of Gaussian random vectors play a key role. Subspace fitting estimates take advantage of the orthogonality between the range (or left kernel) spaces of certain matrix-valued statistics. During the last fifteen years, there has been a growing interest in these methods [12][13][14], their connection to instrumental variables [15] and maximum likelihood [16] approaches, and their invariance properties [17]. They are actually well suited for identifying the system eigenstructure.
Processing output covariance matrices is of interest for long samples of multisensor measurements, which can be mandatory for in-operation structural analysis under nonstationary natural or ambient excitation. The difference between the covariance-driven form of subspace algorithms which is described here and the usual data-driven form [12] is minor, at least for eigenstructure identification [11].

Covariance-driven subspace identification
be the output covariance and Hankel matrices, respectively; and let Direct computations of the R i 's from (1) lead to the wellknown key factorizations [18]: where are the observability and controllability matrices, respectively. In factorization (7), the left factor O only depends on the pair (H, F), and thus on the system eigenstructure of the system in (1), whereas the excitation V k only affects the right factor C through the cross-covariance matrix G.
The observation matrix H is then found in the first blockrow of the observability matrix O. The state-transition matrix F is obtained from the shift invariance property of O: Recovering F requires to assume that rank(O p ) = dim F, and thus that the number of block-rows in H p+1,q is large enough. The eigenstructure (λ, φ λ ) then results from (2).

3
The actual implementation of this subspace algorithm, known under the name of balanced realization (BR) [19] has the empirical covariances (10) substituted for R i in H p+1,q , yielding the empirical Hankel matrix Since the actual model order is generally not known, this procedure is run with increasing model orders [6,20]. The singular value decomposition (SVD) of H p+1,q and its truncation at the desired model order yield, in the left factor, an estimate O for the observability matrix O: From O, estimates ( H, F) and ( λ, φ λ ) are recovered as sketched above. The CVA algorithm basically applies the same procedure to a Hankel matrix pre-and post-multiplied by the covariance matrix of future and past data, respectively [21,22].
A minor but extremely fruitful remark is that it is possible to write the covariance-driven subspace identification algorithm under a form which involves a parameter estimating function. This is explained next.

Associated parameter estimating function
Choosing the eigenvectors of matrix F as a basis for the state space of model (1) yields the following representation of the observability matrix: where Δ Δ = diag(Λ), and Λ and Φ are as in (3). Whether a nominal parameter θ 0 fits a given output covariance sequence (R j ) j is characterized by [15,22]: H p+1,q have the same left kernel space.
This property can be checked as follows. From the nominal θ 0 , compute O p+1 (θ 0 ) using (13), and perform for example an SVD of O p+1 (θ 0 ) for extracting a matrix S such that Matrix S is not unique (two such matrices relate through a post-multiplication with an orthonormal matrix), but can be regarded as a function of θ 0 for reasons which will become clear in Section 6. Then the characterization writes For a multivariable random process Y whose distribution is parameterized by a vector θ, a parameter estimating function [23,24] is a vector function K of the parameter θ and a finite size sample of observations 1 of which the empirical counterpart defines an estimate θ as a root of the estimating equation: Since subspace algorithms exploit the orthogonality between the range (or left kernel) spaces of matrix-valued statistics, the estimating equations associated with subspace fitting have the following particular product form [14,17]: where S(θ) is a matrix-valued function of the parameter and N N is a matrix-valued statistic based on an N-size sample of data. Choosing the Hankel matrix H as the statistics N provides us with the estimating function associated with the covariance-driven subspace identification algorithm: which of course is coherent with (16). The reasoning above holds in the case of known system order. However, in most practical cases the data are generated by a system of higher order than the model. The nominal model characterization and parameter estimating function relevant for that case are investigated in [22].

Other uses of the key factorizations
Factorization (7) is the key for the characterization (16) of the canonical parameter vector θ in (3), and for deriving a residual adapted to detection purposes. This is explained in Sections 6 and 7. Factorization (6) is also the key for (i) designing various input-output covariance-driven subspace identification algorithms adapted to the presence of both known (controlled) inputs and unknown (ambient) excitations [27]; 4 EURASIP Journal on Advances in Signal Processing (ii) designing an extension of covariance-driven subspace identification algorithm adapted to the presence and fusion of nonsimultaneously recorded multiple sensors setups [20]; (iii) proving consistency and robustness results [28][29][30], including for that extension [31].
These three issues are addressed in Sections 3 to 5.

HANDLING KNOWN AND UNKNOWN INPUTS
In some applications, the key issue is to identify the eigenstructure in the presence of both a natural (unknown, unmeasured, and often nonstationary) excitation and a known (measured) input. For example, during flight tests, an aircraft is subject to both atmospheric turbulence and artificial dynamic excitations applied through control surfaces and aerodynamic vanes [32]. The corresponding model writes where (U k ) is the known input, the unknown noises (V k ) and (ε k ) are zero-mean Gaussian white noise sequences, and the three sequences (U k ), (V k ), and (ε k ) are pairwise uncorrelated. Note that the measurement noise (ε k ) does not affect the eigenstructure of the system in (21), and that a moving average sequence (ε k ) can also be encompassed [12,22]. For handling both known and unknown excitations, the use of input/output identification methods is mandatory. Within the framework of the covariance-driven subspace identification algorithm of Section 2, different types of projections can be performed for handling separately or jointly the two types of excitations. Projections are very natural tools within the subspace algorithms landscape [12]. For recovering the eigenstructure, the idea is to project system (21) onto the subspace U generated by all the known inputs, or onto its orthogonal subspace. The projections used in [27] are somewhat nonclassical. They take benefit of the factorization property (6) which holds under two different instances: as above, and where the sequence (W k ) is a measured and white input, the R i 's are the input/output cross-covariance matrices and G is the state/input cross-covariance. Five algorithms have been proposed corresponding to the following approaches.
(i) Eliminating the unknown input V by projecting (21) onto U. (ii) Eliminating the known input U by projecting (21) onto U ⊥ . (iii) Using jointly both projections of (21) onto U and U ⊥ -Variant 1.
(iv) Using jointly both projections of (21) onto U and U ⊥ -Variant 2. (v) Ignoring the presence of the known input U.
The sequence (U k ) is assumed white, except in the second approach.
The performances of these algorithms on real flight test data sets are reported in [27], together with comparisons with several frequency domain polyreference LSCF input/output and output-only eigenstructure identification methods [33,34].

HANDLING MULTIPLE MEASUREMENTS SETUPS
A classical approach in structural analysis, called polyreference modal analysis, consists in processing data measured with respect to multiple references [35]. The common practice is to collect successive data sets with sensors at different locations on the structure: Each record j contains data Y (0, j) k from a fixed reference sensor pool, and data Y ( j) k from a moving sensor pool. The number of sensors may be different in the fixed and the moving pools, and thus in each record j, the measurement vectors Y (0, j) k and Y ( j) k may have different dimensions. This setup, usually referred to as multipatch measurements setup and typically based on about 16 to 32 sensors, can mimic a situation in which hundreds of sensors are available. Processing multipatch measurements data for structural analysis is achieved today by performing eigenstructure identification for each record separately, and then merging the results obtained for records corresponding to different sensor pools. However, pole matching may be not easy in practice, and thus the result of eigenvector gluing may not be consistent.
Instead of merging the identification results, the approach in [20] achieves eigenstructure identification by merging the data of the successive records and processing them globally. The key idea is to elaborate further on the key factorizations properties (6) and (7).
To each record j (1 ≤ j ≤ J) corresponds a state-space realization in the form with a single state transition matrix F, a fixed observation matrix H 0 for the fixed sensor pool, and a specific observation matrix H j corresponding to location j of the moving sensor pool. The problem is to find how to merge the measurements in (24) and adapt the output-only covariancedriven subspace algorithm of Section 2 in order to identify the eigenstructure of F in (25). We focus on the two families Michèle Basseville et al. (26) of which empirical estimates can be computed, for lags i ≥ 0.

of covariances:
In the stationary case, the excitation covariance matrix does not depend on record j: EV , and the cross-covariance between the state and the fixed sensors output, does not depend on j either. Thus all the covariances in (26) factorize with a constant right factor: Consequently, for each lag i ≥ 0, we can stack the R j i 's into a block-column vector: which factorizes as where . The corresponding Hankel matrix factorizes as well: and the algorithm of Section 2 applies to H π . This merge fails under nonstationary excitation: if the input excitation covariance matrix depends on the record index j, the cross-covariance matrix also depends on j, and the factorizations in (28) now write with a record-dependent G: and vector R π i of stacked covariances defined in (29) no longer factorizes as in (30).
To circumvent this difficulty, the idea [20] is to rightnormalize the covariance matrices in (26), (32) to make them looking as if they were obtained with the same excitation. One interesting computational feature of the resulting algorithm is that it mainly amounts to apply the subspace identification algorithm of Section 2 to a Hankel matrix obtained by interleaving the block-columns of the "reference" Hankel matrices and the block-rows of the suitably normalized "moving" Hankel matrices. Experimental results obtained on real data recorded on a bridge have shown the relevance of this algorithm for merging multiple measurements setups and handling the nonstationarities in the data.

ROBUSTNESS TO NONSTATIONARY EXCITATION
In Sections 2 and 4, we have assumed a stationary excitation within the (each) record, with possibly record-dependent covariance matrix. A more realistic assumption is that the excitation covariance matrix is time-varying within each record. Precise mathematical results have been presented in [20,28] which justify the use of the same algorithms as introduced above, without the need for any change in the case of a nonstationary excitation. This justification can be sketched as follows.
When the excitation is nonstationary, so is the recorded signal Y k , and the empirical covariance matrices R i in (10) no longer converge to some well defined underlying R i when the sample size N grows to infinity. Instead, the matrices R i may vary in some arbitrary way. However, the following approximate factorization still holds for N large: where G = 1/N N k=1 X k Y T k and o(1) goes to zero when the sample size N grows to infinity.
Assumptions for this to hold roughly formulate as follows: the covariance matrix of the excitation has to be uniformly bounded, and the nth singular value of the empirical Hankel matrix H built using the R i 's is uniformly bounded from below, where n is the assumed model order (this is a formal version of the requirement that all the modes of the structure should be excited).
The approximate factorization in (33) is the key step in proving the consistency of the covariance-driven subspace identification method in Section 2 and its extension to multiple measurement setups in Section 4 in the case of nonstationary excitation and noises. Note that such nonstationarities may result in time-varying zeros for the underlying system. Hence, likelihood and prediction error related methods do not ensure consistency under such situation, because estimation of poles and estimation of zeros are tightly coupled (Fisher information matrix not block-diagonal). In [20,28], and using martingale techniques, it is shown that the eigenstructure estimate ( λ, ϕ λ ) provided by the subspace methods above is a consistent estimate of the true eigenstructure. Although this theoretical result holds under the assumption of known model order, experimental results suggest that it extends to the practical situation of unknown model order.
A recent generalization of this consistency result to a generic form of subspace algorithms can be found in [29,30], which separates statistical from nonstatistical arguments, therefore enlightening the role of statistical assumptions.
The main conclusion from both the theory and the practice is that the combination of the key factorization property (6) of the covariances and of the averaging operation in the 6 EURASIP Journal on Advances in Signal Processing computation (10) of their empirical estimates, allows to cancel out nonstationarities in the excitation.

BATCH-WISE CHANGE DETECTION AND MODEL VALIDATION
Change detection is a natural approach to fault/damage detection. Indeed the damage detection problem can be stated as the problem of detecting a change in the modal parameter vector θ defined in (3). It is assumed that a reference value θ 0 is available, generally identified using data recorded on the undamaged system. 2 Based on a new data sample Y 1 , . . . , Y N , the damage detection problem is to decide whether the new data are still well described by this parameter value or not. The modal diagnosis problem is to decide which components of the modal parameter vector θ have changed. The damage localization problem is to decide which parts of the structure have been damaged, or equivalently to decide which elements of the structural parameter matrices have changed.
We concentrate here on the damage detection problem for which we describe a χ 2 -test based on a residual associated with the subspace algorithm in Section 2. The modal diagnosis problem can be solved with similar χ 2 -tests focussed onto the modal subspaces of interest, using selected sensitivities of the residual with respect to the modal parameters. The damage localization problem can be solved with similar χ 2 -tests focussed onto the structural subspaces of interest [36][37][38], plugging sensitivities of the modal parameters with respect to the structural parameters of a finite elements model in the above setting. This is described in detail in [39].

Subspace-based residual
For checking whether the new data Y 1 , . . . , Y N are well described by the reference parameter vector θ 0 , the idea is to use the parameter estimating function in (20), 3 namely, to compute the empirical Hankel matrix H p+1,q in (10)- (11) and to define the vector Technical arguments for the √ N factor can be found in [43,44]. Let θ be the actual parameter value for the system which generated the new data sample, and let E θ be the expectation when the actual system parameter is θ. From (16), we know that 2 In case of nonstationary excitation, θ 0 should be identified on long data samples containing as many of these nuisance changes as possible. However, the proposed detection algorithm can be run on samples of much smaller size. 3 Building test statistics on parameter estimating functions is a widely investigated topic; see for example [40][41][42].
namely, vector ζ N (θ 0 ) in (34) has zero mean when θ does not change, and nonzero mean in the presence of a change (damage). Consequently ζ N (θ 0 ) plays the role of a residual. It turns out that this residual has highly interesting properties in practice, both for damage detection [22] and localization [39], and for flutter monitoring [45]. Even when the eigenvectors (mode-shapes) are not monitored, they are explicitly involved in the computation of the residual. It is our experience [39] that this fact may be of crucial importance in structural health monitoring, especially when detecting small deviations in the eigenstructure.

The residual is Gaussian
To decide whether θ = θ 0 holds true or not, or equivalently whether the residual ζ n is significantly different from zero, requires the knowledge of the probability distribution of ζ N (θ 0 ), which unfortunately is generally unknown. One manner to circumvent this difficulty is to assume close hypotheses: where vector δθ is unknown, but fixed. Note that for large N, hypothesis H 1 corresponds to small deviations in θ. This is known under the name of statistical local approach, of which the main result is the following [43,44,[46][47][48].
Let Σ(θ 0 ) Δ = lim N→∞ E θ0 (ζ N ζ T N ) be the residual covariance matrix (it is assumed that the limit exists). Matrix Σ captures the uncertainty in ζ N due to estimation errors: indeed the covariance matrix of the error in estimating θ 0 is this Σ(θ 0 ) as well [43,47]. It should be mentioned also that the estimation of Σ may be somewhat tricky [48,49].
Provided that Σ(θ 0 ) is positive definite, the residual ζ N in (34) is asymptotically Gaussian distributed with the same covariance matrix Σ(θ 0 ) under both H 0 and H 1 ; that is [22]: where J(θ 0 ) is the Jacobian matrix containing the sensitivities of the residual with respect to the modal parameters: As seen in (37), a deviation δθ in the system parameter θ is reflected into a change in the mean value of residual ζ N , which switches from zero (in the undamaged case) to J(θ 0 )δθ in case of small damage. Note that matrices J(θ 0 ) and Σ(θ 0 ) depend on neither the sample size N nor the fault vector δθ in hypothesis H 1 . Thus they can be estimated prior to testing, using data on the safe system (exactly as the reference parameter θ 0 ). In case of nonstationary excitation, a similar result has been proven, for scalar output signals, and with matrix Σ estimated on newly collected data [50].

test for damage detection
Let J and Σ be consistent estimates of J(θ 0 ) and Σ(θ 0 ), and assume additionally that J(θ 0 ) is full column rank (f.c.r.). Then, thanks to (37), testing between the hypotheses H 0 and H 1 in (36) can be achieved with the aid of the following χ 2test: which should be compared to a threshold. Note that the IVbased test proposed in [51] can be seen as a particular case of (39) [22]. In (39), the dependence on θ 0 has been removed for simplicity. The only term which should be computed after data collection is residual ζ N in (34). Thus the test can be computed on-board. Test statistics χ 2 N is asymptotically distributed as a χ 2 -variable, with rank(J) degrees of freedom. From this, a threshold for χ 2 N can be deduced, for a given false alarm probability. The noncentrality parameter of this χ 2 -variable under H 1 is δ θ T J T Σ −1 Jδ θ. How to select a threshold for χ 2 N from histograms of empirical values obtained on data for undamaged cases is explained in [52].
From the expressions in (34) and (39), it is easy to show that this test enjoys some invariance property: any premultiplication of the left kernel S by an invertible matrix factors out in χ 2 N [53]. This is why S defined in (15) can be considered as a function of θ 0 , as announced in Section 2.
The asymptotic properties of the test (39) have been investigated in [54] for the IV-based version, and in [55] in the case of more general (not limited to subspace) estimating functions.

criterion for model validation
In the change detection problem above, one wants to know if some fresh data {Y 0 , . . . , Y n } recorded on a structure are still coherent with a reference structural parameter vector θ 0 identified from data recorded earlier on the same structure. In that problem, a large number of independent data recording experiments is required to get information about the distribution of the damage detection test and assess whether the structural parameters have changed or not.
A different problem, known under the name of model validation, is the following. Only one experiment dataset {Y 0 , . . . , Y n } and one reference signature θ 0 are available. From this, one wants to know if the dataset and the signature match, and decide if some slight modification of the signature can better match the dataset from a statistical point of view. In that problem, a large number of signatures θ have to be tested in order to find the signature minimizing some relevant statistical criterion. This problem has received a wide attention in the identification literature [56].
The idea investigated in [57,58] consists in using the above change detection χ 2 -test as the criterion for model validation: the signatureθ and the dataset are said to match 4 where and χ 2 N is defined in (39). The end result would be either to obtain a better signature by selecting parameters which minimize the validation criterion (41), or to obtain confidence intervals depending on the variations of the validation criterion (41) around its minimum.
Experimental results, obtained on data from both simulations and a laboratory test-bed, are reported in [58] which show the relevance of the model validation criterion (41).

RECURSIVE ON-LINE CUSUM TESTS
In some applications, it is necessary to design detection algorithms working sample point-wise rather than batch-wise. For example, as explained in Section 8, the early warning of deviations in specific modal parameters is required for new aircrafts qualification and exploitation, and especially for handling the flutter monitoring problem.
A simplified although well sounded version of the flutter monitoring problem consists in monitoring a specific damping coefficient. It is known, for example from Cramer-Rao bounds, that damping factors are difficult to estimate accurately [59]. However, detection algorithms usually have a much shorter response time than identification algorithms. Thus, for improving the estimation of damping factors and achieving this in real-time, the idea is to design an on-line detection algorithm able to detect whether a specified damping coefficient ρ decreases below some critical value ρ c [45]: A good candidate for designing this test is the residual associated with subspace-based covariance-driven identification defined in (34), which can be computed recursively as follows: where

EURASIP Journal on Advances in Signal Processing
Since the hypothesis (42) regarding the damping coefficient is not local any more-compare with (36), the asymptotic local approach used in Section 6 can no longer be used for that residual, and another asymptotic should be used instead. From (37) and (43), we know that N−p k=q Z k (θ 0 )/ √ N is asymptotically Gaussian distributed, with mean zero under θ = θ 0 and J(θ 0 )δθ under θ = θ 0 + δθ/ √ N. Now, the arguments in [25,Subsection 5.4.1] lead to the following approximation: for k large enough, Z k (θ 0 ) can itself be regarded as asymptotically Gaussian distributed with zero mean under θ = θ 0 , and the Z k (θ 0 )'s are independent. Furthermore, a change in θ is reflected into a change in the mean vector ν of Z k (θ 0 ). This paves the road for the use of Cusum tests for detecting such changes, according to the type and amount of a priori information available for the parameters to be monitored [44].
For monitoring a damping coefficient (scalar parameter θ a ), the Cusum test writes  (45) and an alarm is fired when g n (θ a ) ≥ γ for some threshold γ [44,Chapter 2]. Since neither the actual hypothesis when this processing starts nor the actual sign and magnitude of the change in θ a that will occur are known, a relevant procedure consists in introducing a minimum magnitude of change ν m > 0, running two tests in parallel, for a decreasing and an increasing parameter, respectively; making a decision from the first test which fires; resetting all sums and extrema to zero and switching to the other one afterwards. This is investigated in [45].
For addressing the more realistic problem of monitoring two pairs of frequencies and damping coefficients possibly subject to specific time variations, 5 multiple Cusum tests for single parameters can be run in parallel. It turns out that the individual subspace-based tests, monitoring respectively each damping coefficient, and each frequency (or sum and difference), appear to behave in a reasonably decoupled manner, and to perform a correct isolation of the parameter which has changed [60].
The advantages and drawbacks of these recursive detection algorithms with respect to those of the recursive subspace identification algorithms described in [61] are investigated in [62].

APPLICATION EXAMPLES
The subspace-based identification and detection methods described above have proven useful in a number of simulated 5 It may be assumed indeed that two modes evolve until superimposition of each other. and real application examples [52,[63][64][65][66][67][68][69]. All these algorithms have been implemented within COSMAD, the modal module of the free INRIA software Scilab [70], and partly within LMS software environment. An overview of results obtained with the subspace-based detection algorithms on several examples is now provided.

Sports car
The proposed method has been applied [68] to detect a fatigue failure of a sports car. The method has been first applied to a reduced scale model, which consists of two vertical plates supported by a very stiff bottom plate. Between the two plates, a mass is connected by four rubber elements. The structure is vertically excited. During the endurance test, the crack initiation period is very short, the accelerometers pick up the changes very soon during the crack growth, and the resonance frequency is decreasing. The global χ 2 -test well detects this fatigue. The car endurance test has first a sports car driven on the endurance track until a fatigue problem of the gear box mounting with the car body occurs. Then a second test car is instrumented to measure the relevant strain and acceleration signals, during an endurance 4-shaker test on a body-in-white equipped with the power train. The objective of this test is to reproduce the same failure in a much shorter time and controlled conditions. The result of the test is that cracks, although less severe, are obtained in exactly the same locations as on the test track. During the test, the acceleration and strain signals are recorded every half hour in order to see whether early detection of the fatigue problem is possible. Two groups of sensors at different locations are evaluated. The first group consists of the 6 sensors on the body and the 4 sensors on the power train. The χ 2 -test value slightly increases during the crack growth, and significantly increases at the end of the crack growth.

Z24 bridge
The proposed method has been applied [52] to the Swiss Z24 bridge, a benchmark of the BRITE/EURAM project SIMCES on identification and monitoring of civil engineering structures, for which EMPA (the Swiss Federal Laboratory for Materials Testing and Research) has carried out tests and data recording. The response of the bridge to traffic excitation under the bridge has been measured over one year in 139 points, mainly in the vertical and transverse directions, and sampled at 100 Hz. The global χ 2 -test has been applied to data of the four reference stations. Thus the test has been evaluated for several data sets, for both the safe and damaged structures.
Two damage scenarios are considered: pier settlement of 20 mm and 80 mm, respectively, further referred to as DS1 and DS2. Even though the effect of the damages on the natural frequencies is really small (no more than 1% for DS1), the χ 2 -test is very sensitive: for DS1, 1000 times larger than for the safe case.

9
The implementation and tuning of an online monitoring system for automated damage detection have also been achieved. Monitoring results based on three sensors have been analyzed, from which the following conclusions have been drawn. The overall increase in the test value is slightly hidden by its daily fluctuations. These fluctuations are due to changes in the modal parameters themselves, due to variations in environmental variables such as temperature, precise hour of measurements, speed of wind, . . . and can be higher than the changes of the modal characteristics due to damage. However, modal variations due to damage imply greater variations of the test than those due to environmental changes.
Another major issue is to take care of the fluctuations of the excitation, due, for example, to changes in the traffic or neighboring activities (a new bridge was in construction a few hundred meters apart), and to avoid running the test when the excitation is clearly different from the excitation of the reference model. A good way to avoid interference between these changes and the test result is to calibrate several reference data sets corresponding to different values of the environmental variables, including excitation and temperature, and to run the test upon matching the environmental characteristics of both the reference and the fresh data sets. Another approach would be to include these variables into the model and consider them as nuisance information. This is the topic of current investigation.

Reticular structure
The method has been applied [64] on a geometrically simple test article designed, assembled and tested dynamically under impact and random shaker excitation. The test structure consists of six cylindrical bars connected in four spherical joints through screwed bolds specially designed according to the requirements of the civil building industry. In order to simulate several damage scenarios, progressive displacements are imposed on the structure by unscrewing one of the joint connectors. The most dramatic damage situations are obtained with the joint completely unscrewed first only in one location, then in two different ones. Sine sweep excitation (30-850 Hz) is applied.
New measurements are taken before and after each of the damage scenarios is applied. For each new scenario, measurements are carried out in four runs: two point locations are used as reference sensors and kept fixed while all other sensors are moved. The global χ 2 -test is applied to the data of a reduced set of sensors. This allows evaluating the test for several data sets both for the healthy and the damaged structure. The method detects damage in an early stage, and it does not require the extraction of the modal parameters from each newly collected data set. This characteristic is very well suited for monitoring purposes: it does not need continuous user interaction and it can easily be made automatic. A remarkable result is the sensitivity of the test to structural changes. The method allows detecting and separating all changes occurring on one node. Increasing stress, single and double collapses are identified by a different order of magnitude in the damage index.

Slat track
The method has also been applied during a fatigue test [69]. During experimental fatigue tests, structural health monitoring is essential to monitor the degradation of the structure with an increasing number of fatigue cycles. Moreover, especially for structures with very high fatigue We added the highlighted "period." Please check. strength, it is important that the test does not have to be interrupted. Since the above damage detection method has the advantage that it operates online, it is a good monitoring candidate for fatigue tests. It has been used in a project aimed at damage detection, life prediction and redesign of a slat track, a device which extends the surface of an airplane wing during takeoff and landing. Since the slat track has very high fatigue strength, testing times can typically take several weeks. Even though the eigenfrequencies of the test structure are not very sensitive to the fatigue crack, the global χ 2 -test above turns out to perform very well, including in comparison with other linear and nonlinear damage indicators. Moreover, the test seems to be robust against nonideal, but typical experimental and data processing issues: 50 Hz magnitude variation, violation of the white-noise assumption, and an incomplete nominal model. In addition, the approach offers the advantage that only output data are needed, and that the nominal model (in terms of modal parameters) has to be determined only once. Afterwards, fresh raw data are simply confronted with this model with these statistics.

Flutter monitoring
A crucial issue in the development of new aircrafts is to ensure the stability of the airplane throughout its operating range. For preventing from a critical instability phenomenon known under the name of aero-elastic flutter, the airplane is submitted to a flight flutter testing procedure, with incrementally increasing altitude and airspeed. The problem of predicting the speed at which flutter can occur is usually addressed with the aid of identification methods achieving modal analysis from the in-flight data recorded during these tests [71,72]. While frequencies and mode-shapes are usually the most important parameters in structural analysis, the most critical ones in flutter analysis are the damping factors, for some critical modes. Until the late nineties, most approaches to flutter clearance have led to data-based methods, processing different types of data. A combined data-based and model-based method has been introduced recently under the name of flutterometer [73].
Algorithms achieving the on-line in-flight exploitation of flight test data are expected to allow a more direct, reliable and cheaper exploration of the flight domain. One important issue is the on-line flight flutter monitoring problem, stated as the problem of monitoring some specific damping coefficients. For improving the estimation of damping factors, and moreover for achieving this in real-time during flight tests, one possible although unexpected route is to resort to detection algorithms able to decide for example whether some damping factor decreases below some critical value or not. The rationale is that detection algorithms usually have a much shorter response time than identification algorithms. This is why the on-line detection algorithms described in Section 7 have been designed. They are based on the subspace-based residual defined in (34), and on the CUSUM test [44]. The monitoring is focussed on specific parameters of interest, such as damping coefficients [45] or pairs of eigenfrequencies subject to specific time variations [60].

FURTHER COMMENTS AND CONCLUSIONS
In this paper, an overview has been provided of the design and investigation of subspace-based algorithms for solving parameter identification, change detection, model validation and data fusion problems arising in the area of modelbased structural analysis and health monitoring of structures in-operation. Some comments are in order, on open problems and ongoing research.
When it comes to vibration-based monitoring of civil engineering structures, it is well known that the dynamics of most of them is affected by the ambient temperature and other environmental effects [74]. This raises the issue of discriminating between changes in modal parameters due to damages and changes in modal parameters due to environmental effects, and in particular the effect of temperature variations. One solution to this problem that is currently investigated consists in using a model of the temperature effect on the structural dynamics, considering this effect as a nuisance parameter, and plugging in the above test a statistical nuisance rejection technique of the type discussed in [75][76][77].
As far as the flight flutter monitoring problem is concerned, the key issue is also to involve more complex models of the underlying physical phenomenon (here the flutter) within the design of the identification and monitoring algorithms. The challenge is whether the monitoring algorithms which will result from these more complex models will better solve the tradeoff efficiency/cost/robustness than the current subspace-based algorithms described in this paper.