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

Falk Hille

Laurent Mevel

Werner Rücker

Structural Health Monitoring with Statistical Methods during Progressive Damage Test of S101 Bridge

Keywords: Subspace methods, Operational modal analysis, Uncertainty bounds, Damage detection, Prestressed concrete bridge

For the last decades vibration based damage detection of engineering structures has become an important issue for maintenance operations on transport infrastructure. Research in vibration based structural damage detection has been rapidly expanding from classic modal parameter estimation to modern operational monitoring. Since structures are subject to unknown ambient excitation in operation conditions, all estimates from the finite data measurements are of statistical nature. The intrinsic uncertainty due to finite data length, colored noise, non-stationary excitations, model order reduction or other operational influences needs to be considered for robust and automated structural health monitoring methods. In this paper, two subspace-based methods are considered that take these statistical uncertainties into account, first modal parameter and their confidence interval estimation for a direct comparison of the structural states, and second a statistical null space based damage detection test that completely avoids the identification step. The performance of both methods is evaluated on a large scale progressive damage test of a prestressed concrete road bridge, the S101 Bridge in Austria. In an on-site test, ambient vi-bration data of the S101 Bridge was recorded while different damage scenarios were introduced on the bridge as a benchmark for damage identification. It is shown that the proposed damage detection methodology is able to clearly indicate the presence of structural damage, if the damage leads to a change of the structural system.

INTRODUCTION

The main objective of research work on Structural Health Monitoring (SHM) systems is to develop or improve methodologies for an effectual identification and handling of risk implying deficiencies of dynamic systems on the base of automated monitoring. Damages lead to changes in the modal parameters of a structure (natural frequencies, damping ratios, mode shapes), hence one eminent component of a decision support framework is the adequate processing of measurement data to reliably detect such changes in the systems under observation [START_REF] Worden | The fundamental axioms of structural health monitoring[END_REF]. The detection of damages is the first and essential step for automated SHM methods, before further levels in the whole damage identification process -consisting of damage detection, localization, quantification and prediction of remaining service life time [START_REF] Rytter | Vibrational based inspection of civil engineering structures[END_REF] -can be treated.

Within the last 30 years a vast literature on vibration-based damage identification methods has emerged. An introduction to the subject as well as an overview of many methods can be found in [START_REF] Farrar | Vibration-based structural damage identification[END_REF][START_REF] Carden | Vibration-based condition monitoring: a review[END_REF][START_REF] Fan | Vibration-based damage identification methods: a review and comparative study[END_REF]. Many of the currently used vibration-based methods for early damage detection performs modal system identification and compare the obtained modal parameter of the actual state with those of the undamaged state [START_REF] Kullaa | Damage detection of the Z24 bridge using control charts[END_REF][START_REF] Magalhães | Dynamic monitoring of a long span arch bridge[END_REF][START_REF] Ramos | Monitoring historical masonry structures with operational modal analysis: Two case studies[END_REF][START_REF] Rainieri | Robust output-only modal identification and monitoring of buildings in the presence of dynamic interactions for rapid post-earthquake emergency management[END_REF][START_REF] Magalhães | Vibration based structural health monitoring of an arch bridge: from automated OMA to damage detection[END_REF]. However, the sensitivity of modal damage indicators is often not satisfying for a lot of real applications in civil engineering with ambient excitation and a limited number of sensors. Varying operational or environmental influences may have a greater impact on the dynamic properties of a system than structural damage in early stages [START_REF] Magalhães | Vibration based structural health monitoring of an arch bridge: from automated OMA to damage detection[END_REF][START_REF] Balmès | Merging sensor data from multiple temperature scenarios for vibration monitoring of civil structures[END_REF]. Besides, civil structures differ from many other mechanical structures by their size and their uniqueness and therefore cannot undergo any kind of prototype testing, as it is common for products from e.g. automotive and aeronautic industry [START_REF] Brownjohn | Structural health monitoring of civil infrastructure[END_REF]. Also, they are subject to operational and environmental loading of mostly large variances and are designed for an exceptional life span compared to mechanical structures. For these reasons, civil structures are provided with a high safety margin, which implies in general a high rigidity. These facts make the development of vibration based strategies for damage detection in an early stage a complex and challenging task for civil structures. Methods that are successfully run on mechanical structures in a certain deterministic environment might not at all be applicable on civil structures.

For the detection of damages in an early stage it is important to take the uncertainty of damage related parameter estimates into account to distinguish small (physical) changes from statistical variability. Since usually only unknown ambient excitation and no artificial excitation is available and practical for continuous and automated monitoring of structures, all parameters that are estimated from the output-only vibration measurements are afflicted with statistical uncertainty. This uncertainty results from the operation conditions, such as finite number of data samples, measurement noises, non-stationary excitations, nonlinear structure, model order reduction and further sources.

Subspace methods have been proven to be very well-suited for the output-only vibration analysis of civil structures in operation due to their excellent theoretical properties [START_REF] Benveniste | Single sample modal identification of a non-stationary stochastic process[END_REF][START_REF] Van Overschee | Subspace Identification for Linear Systems: Theory, Implementation, Applications[END_REF][START_REF] Benveniste | Nonstationary consistency of subspace methods[END_REF]. In this paper two subspace-based methods for vibration monitoring are considered that take the statistical uncertainties into account. First, we use the covariance-driven subspace-based system identification [START_REF] Peeters | Reference-based stochastic subspace identification for output-only modal analysis[END_REF][START_REF] Döhler | Fast multi-order stochastic subspace-based system identification[END_REF] together with their confidence interval estimation [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] for the operational modal analysis. Thanks to this uncertainty quantification, the significance of the modal parameters and their changes during the monitoring period can be evaluated. Second, we perform the actual damage detection with a null space based statistical damage detection test [START_REF] Balmès | Merging sensor data from multiple temperature scenarios for vibration monitoring of civil structures[END_REF][START_REF] Basseville | Subspace-Based Fault Detection Algorithms for Vibration Monitoring[END_REF][START_REF] Basseville | Statistical model-based damage detection and localization: subspace-based residuals and damage-to-noise sensitivity ratios[END_REF][START_REF] Döhler | Subspace-based damage detection under changes in the ambient excitation statistics[END_REF],

where new data from a possibly damaged state is directly compared to a reference state in a test statistic. Like this, the entire structural response is compared without computing any modal parameters. Comparing the test value to a threshold indicates if the structure is damaged.

Both approaches are dealing with the SHM problem from the identification and detection points of view. They are nonetheless related by the principles of subspace methods and the same fundamental uncertainty quantification computation from the output-only measurement data.

Our contribution to these methods concerns their applicability: Being elaborate methods but lacking some feasibility in practice, both methods were recently enhanced with a strongly decreased computational burden, feasibility of high model orders and a more numerically robust computation [START_REF] Döhler | Fast multi-order stochastic subspace-based system identification[END_REF][START_REF] Döhler | Efficient multi-order uncertainty computation for stochastic subspace identification[END_REF][START_REF] Döhler | Subspace-based damage detection under changes in the ambient excitation statistics[END_REF]. The results are fast algorithms that can be applied easily to vibration data from civil, mechanical or aeronautical structures. In this paper, they are demonstrated on the large scale progressive damage test of a prestressed concrete road bridge, the S101 Bridge in Austria [START_REF]Progressive Damage Test S101 Flyover Reibersdorf[END_REF]. In an on-site test, ambient vibration data of the S101 Bridge was recorded while different damage scenarios were introduced on the bridge as a benchmark for damage identification. The verification of the effectiveness of SHM is particularly important on such structures, especially since the effort for maintenance and repair of the prestressed concrete road bridges increased dramatically in the last decades [START_REF] Zilch | Assessment of the condition of prestressed or posttensioned German concrete bridges[END_REF]. It is shown that the proposed methodology is able to clearly indicate the presence of structural damage, if the damage leads to a change of the structural system. Small damages which do not result in a system change when not activated by loading, do not lead to a modification of the dynamic response behavior and for that cannot be detected with the proposed global monitoring method. This paper is organized as follows. In Section 2, the basics of stochastic subspace identification are given. In Section 3, the confidence interval computation of modal parameters is de-scribed and in Section 4, the damage detection strategy is summarized. These algorithms are applied to S101 Bridge in Section 5, before finishing with some concluding remarks in Section 6.

STOCHASTIC SUBSPACE IDENTIFICATION (SSI)

Stochastic Subspace Identification methods are among the state of the art methods for modal parameter estimation [START_REF] Van Overschee | Subspace Identification for Linear Systems: Theory, Implementation, Applications[END_REF][START_REF] Benveniste | Nonstationary consistency of subspace methods[END_REF]. In this section, an overview of the covariance-driven identification algorithm is given, whose properties are also the foundation of the subsequent damage detection algorithm.

Models and Parameters

The behaviour of a mechanical system is assumed to be described by a stationary linear dynamical system

( ) ( ) ( ) ( ), t t z t v z z t + + = ɺɺ ɺ M C K (1)
where t denotes continuous time, , ,

d d × ∈ ℝ M C K
are the mass, damping and stiffness matrices, high-dimensional vector d z ∈ℝ collects the displacements of the degrees of freedom of the structure and the external force v is unmeasured and considered as noise. The eigenstructure of (1) with the modes µ i and mode shapes (observed eigenvectors)

r i ∈ ℝ ψ is a solution of 2 2 det( ) 0, ( ) 0, i i i i i i i M C K M C K L + + = + + = = µ µ µ ψ µ φ φ , (2)
where matrix

r d L ×
∈ ℝ maps the r sensor locations to the d degrees of freedom of the structure.

Observing model (1) at the r sensor locations (e.g. by acceleration, velocity or displacement measurements) and sampling it at some rate 1/τ yields the discrete model in state-space form

1 k k k k k k x Ax v y Cx w + = + = + , (3)
where

n n A × ∈ ℝ is the state transition matrix, r n C × ∈ ℝ is the observation matrix, n k
x ∈ ℝ are the states of the system and r k y ∈ ℝ the output measurements at the discrete times t k = τ , where n is the system order. The vectors v k and w k are the unmeasured input and output disturbances. The eigenstructure of system (3) is given by det( ) 0, ( ) 0,

i i i i i A I A I C - = - = = λ λ φ ϕ φ . (4)
Then, the eigenstructure of the continuous system (1) is related to the eigenstructure of the discrete system (3) by e ,

i i i i = = τµ ψ λ ϕ . ( 5
)
From the eigenvalues λ i , the natural frequencies f i and damping ratios i ξ of the system are direct-

ly recovered from 2 2 2 2 100 , 2 i i i i i i i a b f a b b - = = + + ξ πτ
where

Im( | arctan |, ln | | R ) ( ) e i i i i i a b = = λ λ λ . ( 6
)

Covariance-Driven SSI

To obtain the eigenstructure of system (3) from measurements (y k ) k=1,…,N , the stochastic subspace identification algorithm is used in its covariance-driven form [START_REF] Benveniste | Single sample modal identification of a non-stationary stochastic process[END_REF][START_REF] Peeters | Reference-based stochastic subspace identification for output-only modal analysis[END_REF]. In the case of a high number of available sensors r, a subset of 0 r r ≤ sensors can be chosen as reference sensors (or socalled projection channels) in the output vector

1 N T i k i k r r k y R y N + = × ∈ = ∑ ℝ (7)
are computed for the lags i = 1, ..., p + q, where the parameters p and q are chosen such that 0 min{ , } pr qr n ≥ , with in general p + 1 = q. With these correlations, the Hankel matrix

1 2 2 3 1 1 2 q q p p p q R R R R R R R R R + + + +       =         ⋯ ⋯ ⋮ ⋮ ⋱ H (8)
T V U U         ∆ ∆ = 0 1 0 1 H , def 1/ 1 ( 1 2 ) 1 p r p n C CA U CA + ×       = ∆ =       ∈   ℝ ⋮ O . ( 9
)
From the observability matrix O the matrices C in the first block row and A from a least squares solution of

A = O O with 1 p C CA CA -       =         ⋮ O , 2 p CA CA CA       =        
⋮ O [START_REF] Magalhães | Vibration based structural health monitoring of an arch bridge: from automated OMA to damage detection[END_REF] are obtained. The eigenstructure (λ i ,φ i ) of the system (3) and the modal parameters are finally obtained from ( 4)- [START_REF] Kullaa | Damage detection of the Z24 bridge using control charts[END_REF].

In Operational Modal Analysis, the model order n in ( 9) is unknown and spurious modes appear due to the unknown noise characteristics in the ambient excitation. Based on the observation that physical modes remain quite constant when estimated at different over-specified model orders, while spurious modes vary, they can be distinguished using so-called stabilization diagrams. The system is identified truncating in (9) at multiple model orders, and frequencies from this multi-order system identification are plotted against the model order [START_REF] Peeters | Reference-based stochastic subspace identification for output-only modal analysis[END_REF][START_REF] Döhler | Fast multi-order stochastic subspace-based system identification[END_REF]. From the modes common to many models and using further stabilization criteria, such as threshold on damping values, low variation between modes and mode shapes of successive orders etc., the final estimated model is obtained. Like this, stabilization diagrams provide a GUI where the user is assisted in selecting the identified modes of an investigated structure.

CONFIDENCE INTERVALS OF MODAL PARAMETERS

The uncertainty quantification of the estimated modal parameters that are obtained from subspace identification in the previous section is necessary to assess the confidence one can have in these values, e.g. when comparing the modal parameters of different states of a structure. Modal parameters with little confidence (and hence large confidence intervals) are little useful for comparing structural states.

The uncertainties of the modal parameters at a chosen system order can be computed from the uncertainty of the Hankel matrix H by doing a sensitivity analysis. Since the modal parameters are functions of the Hankel matrix, it holds (vec ), (vec ), (vec ),

i i i i i f i f ∆ = ∆ ∆ = ∆ ∆ = ∆ ϕ ξ ξ ϕ J J J H H H (11)
for (theoretical) first-order perturbations of the frequencies f i , damping ratios ξ i and mode shapes ϕ i , where J denotes their sensitivity with respect to the stacked Hankel matrix vec 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][START_REF] Döhler | Uncertainty quantification for modal parameters from stochastic subspace identification on multi-setup measurements[END_REF]. The uncertainty of the stacked Hankel matrix is quantified by its covariance Σ H , which can be easily evaluated by splitting the available sensor data into blocks on which instances of the Hankel matrix are obtained. Then, Σ H is computed from their sample covariance. It follows cov( ) , cov( ) , cov( ) ,

i i i i i i i T T T i f f i f = Σ = Σ = Σ ξ ϕ ξ ϕ ξ ϕ H H H J J J J J J (12)
from where the standard deviations of the modal parameters are obtained. This offers a possibility to compute confidence intervals on the modal parameters directly in a fully automated way, without repeating the system identification step. In [START_REF] Döhler | Efficient multi-order uncertainty computation for stochastic subspace identification[END_REF] an efficient and fast computation scheme of the uncertainty computation laid out in ( 12) is described in detail, where confidence bounds up to high model orders can be computed within reasonable computation times. The efficient computation at multiple model orders for stabilization diagrams is also addressed in [START_REF] Döhler | Efficient multi-order uncertainty computation for stochastic subspace identification[END_REF].

STATISTICAL SUBSPACE-BASED DAMAGE DETECTION

Instead of performing damage detection by comparing the modal parameters from a new data set to a reference state, the statistical subspace-based damage detection algorithm [START_REF] Basseville | Subspace-Based Fault Detection Algorithms for Vibration Monitoring[END_REF][START_REF] Basseville | Statistical model-based damage detection and localization: subspace-based residuals and damage-to-noise sensitivity ratios[END_REF] compares different structural states by using a χ 2 test on a residual function. This residual confronts the left null space of a subspace matrix H 0 in the reference state with a subspace matrix H of a new, possibly damaged state. This algorithm has the advantage that no modal parameter extraction is necessary, while the χ 2 test checks for changes implicitly. Furthermore, the whole system response is compared implicitly and not only the first modes as it is often the case when monitoring modal parameters.

The algorithm consists of the following steps. In the reference state, a subspace matrix H 0 is computed from the output data of the system as in [START_REF] Ramos | Monitoring historical masonry structures with operational modal analysis: Two case studies[END_REF]. Note that reference scenarios with different environmental conditions can be mixed when computing H 0 , see [START_REF] Balmès | Merging sensor data from multiple temperature scenarios for vibration monitoring of civil structures[END_REF], such that the damage detection test will be robust to these environmental conditions. Then, the left null space S of H 0 is computed and the characteristic property of a system in the reference state writes S T H = 0. The null space S can be obtained with an SVD of H 0 as in ( 9) from S = U 0 . The associated residual vector [START_REF] Balmès | Merging sensor data from multiple temperature scenarios for vibration monitoring of civil structures[END_REF][START_REF] Basseville | Subspace-Based Fault Detection Algorithms for Vibration Monitoring[END_REF][START_REF] Basseville | Statistical model-based damage detection and localization: subspace-based residuals and damage-to-noise sensitivity ratios[END_REF] writes vec( )

T N N S ζ = H , ( 13
)
where N is the number of samples, on which H is computed. Furthermore, compute the residu-

al's covariance Σ = E[ζ N ζ N T
] in the reference state, where E denotes the expectation operator. This is done by cutting the data into blocks and computing the residual on each of the data blocks, analogous to Σ H in the previous section, and it holds Σ = N (I ⊗ S T ) Σ H (I ⊗ S), where I is the identity matrix and ⊗ denotes the Kronecker product. The covariance of the null space S is not computed, since it was neglected in the derivation of the original detection methods [START_REF] Basseville | Subspace-Based Fault Detection Algorithms for Vibration Monitoring[END_REF][START_REF] Basseville | Statistical model-based damage detection and localization: subspace-based residuals and damage-to-noise sensitivity ratios[END_REF],

where S is assumed to be model-based and not data-based. Note that while case studies (e.g. in [START_REF] Balmès | Merging sensor data from multiple temperature scenarios for vibration monitoring of civil structures[END_REF][START_REF] Döhler | Subspace-based damage detection under changes in the ambient excitation statistics[END_REF]) performed well without using the covariance of S, it should be taken into account in future research to increase the robustness of the detection method, requiring a non-trivial extension of the approach (see e.g. [START_REF] Dong | Robust fault detection with statistical uncertainty in identified parameters[END_REF]).

The residual vector ζ N is asymptotically Gaussian with mean zero in the reference state. It manifests itself to damage by a change in its mean value, corresponding to an increase of the

χ 2 test 2 1 T N N χ ζ ζ - = Σ . ( 14
)
The monitoring of the system is achieved by calculating the residual vector ζ N in (13) on the Hankel matrices estimated from newly recorded output data and subsequently the χ 2 test value in [START_REF] Van Overschee | Subspace Identification for Linear Systems: Theory, Implementation, Applications[END_REF], which is compared it to a threshold. Such a threshold is obtained empirically from the χ 2 test values of several data sets in the reference state of the structure for a desired type I error. A significant increase in the χ 2 value indicates that the system is no more in the reference state.

The residual covariance matrix Σ is usually high-dimensional and its inversion Σ -1 costly and numerically fragile. In [START_REF] Döhler | Subspace-based damage detection under changes in the ambient excitation statistics[END_REF] the authors proposed an efficient and numerically robust computation. To further increase the robustness of the computation, while decreasing the computation time, so-called reference sensors or projection channels [START_REF] Peeters | Reference-based stochastic subspace identification for output-only modal analysis[END_REF] are used in the computation of the subspace matrix as indicated in Section 2.2, reducing the size of the involved matrices significantly. With these improvements a fast and robust computation of the damage indicator (14) was made.

MONITORING OF S101 BRIDGE DURING PROGRESSIVE DAMAGE TEST

Progressive Damage Test of S101 Bridge

Within the European research project "Integrated European Industrial Risk Reduction System (IRIS)" a prestressed concrete bridge was artificially damaged [START_REF]Progressive Damage Test S101 Flyover Reibersdorf[END_REF]. The intention was to provide a complete set of monitoring data during a defined loss of structural integrity for testing and evaluation of various SHM methods and applications. Therefore, the static and dynamic behavior of the structure was measured permanently during the 3-day damaging process. The progressive damage campaign was planned and organized by the Austrian company VCE. The characteristics of the structure under observation, the measurement campaign and the introduced damages are presented briefly.

The S101 Bridge

The S101 Bridge was a prestressed concrete bridge from the early 1960s and therefore a characteristic representative of the partly invalid highway infrastructure asset in Europe. Despite a general lack of experience in the time of their design, prestressed concrete was a very popular construction type in those days. Retrospectively, some of the major design assumptions proofed to be erroneous and after short periods of operation significant and characteristic damage patterns occurred at the structure [START_REF] Zilch | Assessment of the condition of prestressed or posttensioned German concrete bridges[END_REF][START_REF] Polder | Non-traditional assessment and maintenance methods for aging concrete structures -technical and non-technical issues[END_REF]. The load bearing capacity and especially the durability of the bridges remained mostly limited despite costly retrofitting activities [START_REF] Grünberg | Fatigue life of concrete bridges under consideration of increasing traffic loads (in German)[END_REF][START_REF] Khan | Bridge and Highway -Structure, Rehabilitation and Repair[END_REF]. In addition, an increase of heavy load vehicle traffic has been ongoing since the start of operation. Thus, prestressed bridge structures in central Europe and other countries have been exposed to degradation processes initiated by poor design, while the operational loading and the associated dynamic stressing increases steadily in the same time [START_REF] Brady | COST 345 -Procedures Required for the Assessment of Highway Structures[END_REF][START_REF] Colditz | Bridge strengthening -a necessary prerequisite for a reliable highway network (in German)[END_REF].

Figure 1. Bridge S101 during damage test [START_REF]Progressive Damage Test S101 Flyover Reibersdorf[END_REF].

Damage Description

The progressive damage test took place from 10-13 December 2008. During the test the highway beneath the bridge was open in one direction. The second direction was closed for traffic because of construction works which in addition took place near the bridge.

Two major damage scenarios were artificially induced. First, a damage of one of the four columns was inserted by cutting through the column on its lower end. With this action a change in the global structural system was implemented since an extra hinge was formed just above the foundation. During the cutting process no signs of extra cracking (noise emission) were observed. After a second cut a 5 cm thick slice of the column was removed and the column was lowered for altogether 3 cm, while moderate noising was noticed. At the neighboring column a horizontal crack occurred. Additional cracking in the bridge deck was not observed. The column was settled until reaching elastic limits of the bridge deck. Because of the hydraulic jack there was no complete loss of support. The column was secured in its location by inserted steel plates.

Afterwards the column was uplifted again to its original position and secured there again by steel plates.

Second, prestressing tendons of one of the beams were cut successively for a further damage scenario. All in all three and a quarter of a wire bundle were cut through. Between each intersection pauses of several hours were kept to let the structural system change into a new state of equilibrium, while no indication of changing structural conditions like cracking was observed.

Since the bridge was not loaded by traffic, it is assumed that cutting the limited number of tendons reduced only the extra prestressing margin.

Pictures of both scenarios of the progressive damaging are shown in Figure 2. In Table 1 all damage actions are sorted in chronological order.

Measurement Description

For the vibration measurements a BRIMOS ® measurement system with a permanent sensor grid was used. The grid consisted of 15 sensor locations on the bridge deck (see Figure 3), in each location three sensors for measurements in the bridge deck's vertical, longitudinal and transversal direction. Altogether, 45 acceleration sensors were applied. The sensor layout was motivated by the symmetry of the bridge, so that 14 locations were set only on one side of the bridge deck, while the 15 th location on the other side has the purpose to distinguish between bending and torsional modes.

The data were recorded permanently with a sampling frequency of 500 Hz. During the three days measurement campaign 714 data files each containing 45 channels with 165 000 data points were produced, corresponding to 5.5 min of measurements each. Note that there were hardly any temperature changes during the measurement period, where misty winter weather just below freezing was dominant, so nuisance from temperature changes can be excluded in the measurements.

Figure 3. Location of the sensors on the bridge deck.

Modal Parameters and their Confidence Intervals

Being interested in the identification of the first modes in a frequency range [0-18 Hz], the data was downasampled and decimated from sampling rate 500 Hz by factor 5. System identification and the covariance computation was done with the covariance-driven subspace identification detailed in Section 2 with parameters p + 1 = q = 35 at model orders n = 1,…,100. All r = 45 sensors were used and r 0 = 3 reference sensors were chosen. The covariance computation on all identified modes at the different model orders was accomplished with the strategy explained in Section 3. For the covariance estimate on the Hankel matrix, the data was cut into 100 blocks.

System Identification Results in Reference State

First, a modal analysis of the structure in the reference state with the computation of the uncertainty bounds of the modal parameters as in Section 3 was performed. In Figure 4(a), the stabilization diagram of the natural frequencies vs. the model order is presented, where a confidence interval (± one standard deviation σ f ) of each frequency is plotted as a horizontal bar.

In Figure 4 In Figure 5, the frequency and damping ratio of first mode identified at the different model orders is magnified for a better visibility of the confidence bounds. Using the information from these confidence bounds, values for the frequency and the damping ratio of a mode can be chosen that are optimal for different model orders, as the true model order is unknown.

f f f = ⋅ ɶ σ σ
, the damping ratios ξ and their variation coefficient From the stabilization diagrams, the modal parameters are chosen. In Table 2, the system identification results with their variation coefficients are summarized. As expected, the variation coefficients of the frequency estimates are very low (lower than 0.5%), while the estimates of the damping ratios show much higher variation coefficients (up to 20% in this case).

/ 100 = ⋅ ɶ ξ ξ ξ σ σ mode f (Hz) f ɶ σ ξ (%) ɶ ξ σ 1
In Figure 6, the real parts of the obtained mode shapes in the vertical direction at the 14 sensors of one side of the bridge deck are displayed with their uncertainty bounds (±2 standard deviations). From the 15 th sensor on the other side of the bridge deck, whose contribution is marked as the red point in Figure 6, information about the kind of the mode is obtained: Modes 1, 3 and 5 are vertical bending modes and modes 2 and 4 are torsional modes. Since the 15 th sensor was very close to the middle of the bridge, the distinction of the mode shapes 3 and 4 was not easy due to symmetry. However there is a small but significant shape difference of more than the two standard deviations at the 15 th sensor between modes 3 and 4, where the former is negative and the latter is positive, while the mode shape at the other sensor in the middle of the bridge deck is negative in both cases in Figure 6.

System Identification Results in Maximal Damaged State

In the next step, the modal analysis with the confidence interval computation was performed in the maximal damaged state on an example record after the settling of the column and the insertion of the steel plates, corresponding to damage event F in Table 1. The analysis is performed with the confidence interval computation analogously to the previous section and the results are presented in Table 3 and Figure 7. In comparison to the modes in the reference state in the previous section, it can be noted that the mode shapes have changed significantly especially at the location of the bridge that was settled (see Figure 3). These changes are higher than the associated confidence bounds. They are also confirmed by a comparison of the modal assurance criterion (MAC) between the modes in Figure 8. Mode 3, which was a vertical bending mode in the reference state, is not present anymore.

Instead, two torsional modes appear that are denoted as 3a and 3b, where mode 3b has no similarity to any of the modes from the reference state as can be seen in Figure 8. Note that modes 3 and 4 from the reference state are hardly distinguishable in the comparison of MAC values. This is because they are bending and torsion modes that only show a difference at the one sensor on the other side of the bridge deck, where however the deflection is close to zero for these modes.

Torsional mode 4 and bending mode 5 in the damaged state show similar behavior as in the reference state with changes around the settled column. Concerning the natural frequencies in comparison to the reference state, it can be observed that modes 1 and 2 show a frequency drop (0.5% and 6%, respectively), while the frequencies of modes 4 and 5 show an increase (7% and 3%, respectively). The damping ratios of the modes 1, 2, 4 and 5 have increased.

Results on system identification during all damage scenarios are presented in the next section.

Frequency Monitoring During Progressive Damage Test

During the progressive damage test of the S101 Bridge, more than 700 datasets were available.

Some of them contained erroneous data due to destruction work on the bridge or other influences that were excluded from the analysis. On the left 680 datasets, an automated monitoring procedure was applied, which did the system identification and confidence interval computation automatically for each dataset. This means that a stabilization diagram was built for each dataset. In the diagrams, only modes in accordance with reasonable thresholds for the damping estimates and confidence interval bounds were preselected. The automatic mode selection was then per-formed based on a simple statistic clustering algorithm, where each cluster contains modes with small deviations of the frequency, damping and MAC values between successive model orders.

Clusters with a minimal number of modes were chosen and the cluster centers were used as the final results. Note that it was not our purpose to develop refined automated Operational Modal Analysis algorithms (as e.g. in [START_REF] Magalhães | Vibration based structural health monitoring of an arch bridge: from automated OMA to damage detection[END_REF][START_REF] Reynders | Fully automated (operational) modal analysis[END_REF]), but to demonstrate the confidence interval computation of the modal parameters during the progressive damage test. In this context, no tracking of the individual modes was realized, which would also prove to be difficult given the strong changes in the modes and even in the number of modes, as described in the previous section.

The results of the frequency monitoring of all datasets during the progressive damage test are shown in Figure 9, where the respective damage scenarios were explained in Table 1. Especially the frequency change can be clearly seen when one column of the bridge was lowered before it was lifted up again (between A and G). This affected mainly the second, third and fourth mode, while the frequency changes in the first mode were less significant. An evaluation of the changes in the fifth mode is difficult due to its high uncertainty in the automated evaluation, but its change was shown on a dataset in the damaged condition in the previous section. The change in the third mode is remarkable, as it was replaced by two modes during the lowering of the column, with one lower and one higher frequency than before.

The change in the frequencies when cutting the tendons (between G and K) is not significant.

Only after the uplifting of the column and before cutting the first tendon (between G and H), some of the frequencies are dropping, probably due to the settling of the structural system after the uplifting. However, no significant frequency change can be observed afterwards.

Damage Detection

Data Analysis

The vibration measurements were recorded throughout the whole three day long damage test, including the nights. The measurement campaign started approximately 12 hours before the damage of the column was executed to monitor the undamaged state for an adequate time period.

The reference state of the undamaged structure was set up by computing and averaging the reference matrices from several datasets. With this approach [START_REF] Balmès | Merging sensor data from multiple temperature scenarios for vibration monitoring of civil structures[END_REF], a possible disturbance by single excitation events or different environmental conditions is minimized. The size of the Hankel matrix depends on the number of chosen output channels as well as on the time lags, to be accounted for computing the output covariance matrices.

The size of the analyzed null space S depends on the order of the system and can be generally determined by the rank of the matrix. Because of noise, the singular values will not drop to ex-actly zero and the system order must be chosen. The measured response signals of S101 Bridge were dominated by around ten modes, which led us to a choice of n = 20 as the system order.

For calculation of the residual covariance matrix Σ it is important to use as much data sets as are required for the matrix to have full rank and a precise estimation. It is convenient to divide available datasets into subsets to provide the necessary amount of data sets, whereas the number of data points for one residual vector should not be too small to be significant.

In the test stage, the χ 2 test statistics in ( 14) is computed for every data set, which in real time means an indicator of damage for every 5.5 min. Note that no modal parameter estimation or tracking is necessary for this damage detection approach.

Results of Damage Detection

Figure 10 shows a bar plot of values as damage indicators of all consecutive tests within the three days campaign, where all 45 sensors (r = r 0 = 45) were used. For the computation of the residual covariance matrix 100 datasets of the undamaged state were used. The abscissa of the plot describes the chronological sequence of the damage activities as noted in Table 1 as well as the 6am and 6pm points of time for orientation. The reaction of the damage detection test at the main damage events can be easily observed.

From the measurements in the reference state (before damage event A) a threshold of the damage indicator can be set easily, such that all the subsequent scenarios after damage event B would be correctly classified as damaged. This qualifies the proposed method for damage detection and in the following the obtained results are analyzed in more detail.

It can be seen that the damage indication is interfered by noise in the ambient excitation of the bridge. The χ 2 values periodically swell up in the morning and ebb away in the evening. It is assumed that the traffic going underneath the bridge and/or the construction work nearby are the source of the noise. Since the measurements in the reference state were only taken during 12 hours in the night, an extension of the measurement time period of the undamaged state to a whole day circle might reduce these disturbances significantly when taken into account in the residual covariance matrix. Influences by solar radiation and/or temperature alternation can be excluded, since during the 3-days campaign misty winter weather with only moderate temperature changes just below freezing was dominant. As can be seen in Figure 11, the three steps of the column settlement action are very distinctive in their influence on the computed damage indicator. Obviously, the dynamic system has changed to quite some extent and the elastic settlement of altogether 27mm can be clearly detected. Although not that strongly visible, the cutting of the column (A+B) also caused an in-crease of the indicator of approximately 75%. However it has to be mentioned that the absolute effect is superimposed by noise effects. The inlet in Figure 11 shows a detail of that time period. The column remained in the settled condition for approximately one day and was then uplifted again in its former position (event G). The effect of the uplifting is again clearly visible in Figure 10 by a drop of the damage indicator. However, the indicator did not drop completely to its original value, which is certainly due to the fact that the lowering of one column has led to cracking within the concrete structure to some extent and hence to a change of the dynamic signature of the system.

Detection of Column Damage

Detection of Tendon Damage

As one can see in Figure 10 the cutting of the prestressing tendons did not lead to a significant change in the damage indicator after the single cutting steps. Nevertheless, a distinctive increase of the indicator could be observed at the end of the measurement. Figure 12 shows the last time period in detail. It is possible that a change of the bridge's structural system took place with a time delay after cutting partly the fourth tendon. For a prestressed concrete structure the loss of prestressing is a major damage which comes along with a significant loss of its load bearing capacity.

One reason that the cutting of the tendons does not affect the proposed damage indicator might be that the overall prestressing is designed for combinations of dead and traffic loads.

Since the bridge deck of S101 has a quite slender cross section, the dead load is not that high in comparison with the maximum design traffic loading. An additional dead load reduction comes from the removal of the asphalt surface before the damage test. Furthermore it has to be recalled that the prestressing was designed with adequate safety margins. Also, the cutting results only in local loss of prestressing because of the bound between tendons and concrete. For these specified reasons it is assumed that the cutting of the tendons during the damage test did not lead to a significant change of the global structural system, since the loading of the structure was not high enough to activate the damage right after its insertion.

The increase of the indicator at the end of the test series might be the result of a delayed rearrangement of the structural system by a reduction of stresses under generation of cracks in the concrete of the bridge deck. Though, an evidence of that assumption, for instance by an increase of the measured bridge deck deflection, could not be found.

Projection Channels

With the application of so-called reference sensors or projection channels (see Section 2.2) the computation costs of the analysis of the damage indicators was reduced massively. The full analysis of all 45 channels required an extent of memory space which commonly is not provided on desktop computers.

On S101 several numbers and constellations of projection channels were analyzed. Figure 13 shows a χ 2 value plot for only r 0 = 4 well-chosen projection channels (the fifth sensor from the right with all three directions and the 15 th sensor on the other side of the bridge deck in the vertical direction, see also Figure 3). As one can see, almost equal information about the damage indication could be achieved compared to Figure 10, and at the same calculation the computing time was cut to a fifth compared to a complete sensor analysis.

CONCLUSION

In this work, two approaches have been presented for the monitoring of dynamic systems. Both are based on Stochastic Subspace Identification and take the uncertainties in the output-only vibration data into account. With the first approach, modal parameter estimates and their confi-dence intervals are computed, which are essential for knowing the significance of the obtained estimates. With the second approach, a statistical damage indicator is computed to automatically evaluate changes in the structural response. Both approaches represent in fact the same paradigm, with the difference that the first approach evaluates the uncertainty in each modal parameter individually to be able to compare them directly, and the comparison is still to be done by the engineer, while in the second approach the change in the modal parameters is automatically assessed and quantified by the implicit use of these uncertainties in a damage detection index.

Newly developed improvements and extensions of the basic methodologies are integrated to increase numerical robustness, decrease computation efforts as well as to refine results. With these improvements, both theoretically profound methods can be easily used in practical applications. They were successfully applied to output-only vibration data from the progressive damage test of the S101 Bridge, where it was shown how confidence intervals of the modal parameters improve the assessment of system identification results. With the recent improvements made on the method, confidence intervals can now be efficiently computed up to high model orders. Together with an automated modal analysis (e.g. as in [START_REF] Magalhães | Vibration based structural health monitoring of an arch bridge: from automated OMA to damage detection[END_REF]) the confidence intervals provide a supplementary tool to analyze the structural condition individually. In the progressive damage test, a clear link between the change in the modal parameters and the destruction states can be made.

With the statistical damage detection test, the change in the system response of the entire structure is evaluated without the need for an individual modal analysis. The link between the structural changes due to the artificially introduced damage cutting/lowering of one column and the behavior of the damage indicator at each test stage was clearly shown. For a second damage scenario, the cutting of single prestressing tendons, an early stage change of the dynamic response behavior could not be indicated. It is presumed that the locally acting damage did not activate a significant change of the overall structural system due to the absence of operational loading during the test. Under the assumption of an adequate preliminary monitoring of the un-damaged structure and taking into consideration the preceding assertion about damage activation, the statistical damage detection method proves to be feasible for structural health monitoring of civil engineering structures.

In our case study no relevant temperature or other environmental changes were present, so the robustness of the presented methods to such changes could not be shown. Note that the rejection of environmental influences in the modal parameters was e.g. done in [START_REF] Magalhães | Vibration based structural health monitoring of an arch bridge: from automated OMA to damage detection[END_REF], and the confidence interval computation in the context of such a rejection approach would be an interesting extension of the presented approach for health monitoring through operational modal analysis. The rejection of environmental influences for the statistical damage detection approach was already described and successfully performed in [START_REF] Balmès | Merging sensor data from multiple temperature scenarios for vibration monitoring of civil structures[END_REF]. Thus, our presented statistical methods provide a powerful basis for Structural Health Monitoring tools that can be efficiently and reliably used on civil engineering structures.

  ∈ ℝ . Selecting reference sensors increases the computation time and improves the quality of the results[START_REF] Peeters | Reference-based stochastic subspace identification for output-only modal analysis[END_REF]. Then, the output correlations

  is filled, which possesses the factorization property H = O C into observability matrix O and the stochastic controllability matrix C. The observability matrix O is obtained from H by a singular value decomposition (SVD) and truncation at the desired model order n ( )

Figure 2 .

 2 Figure 2. Progressive damage of column and tendons on bridge S101 [25].

  (a) it can be seen that the (true) structural modes seem to have much lower uncertainty bounds than spurious modes. With this observation, a threshold of 1% was put on the variation coefficient of the frequencies (standard deviation divided by frequency), which leads to a much clearer diagram in Figure 4(b). Such a cleaned diagram is also easier to evaluate for an automated modal parameter extraction.

Figure 4 .

 4 Figure 4. Stabilization diagrams containing the natural frequencies of the first five modes with

Figure 5 .

 5 Figure 5. Zoom on first mode: estimated natural frequencies (left) and corresponding damping

Figure 6 .

 6 Figure 6. First five mode shapes in vertical direction on one side of the span (line) and the sensor

Figure 7 .

 7 Figure 7. Mode shapes in range [0-18 Hz] in vertical direction on one side of the span (line) and

Figure 8 .

 8 Figure 8. MAC values between the mode shapes from the reference and damaged states.

Figure 9 .

 9 Figure 9. Natural frequencies with confidence bounds during progressive damage test of S101

Figure 10 .

 10 Figure 10. Course of the damage indicator over the 3-day damage test.

Figure 11

 11 Figure11shows the damage indicator during the several steps of the first damage scenario, the

Figure 11 .

 11 Figure 11. Damage indicator for cutting and settling of one bridge column, detail of damage in-

Figure 12 .

 12 Figure 12. Detail of damage indicator at the end of the tendon cutting process.

Figure 13 .

 13 Figure 13. Damage indicator on S101 with four projection channels.

Table 1 .

 1 Notation of consecutive damage actions and their effects with time index after start of the measurements.

	Time Index Damage Action	Damage Effect
	A 14h 52m	Begin of cutting through	• no signs of extra cracking (noise) or in-
		column	crease of existing cracks
	B 16h 15m		• formation of an extra hinge just above the
		End of second cut through
			foundation, which itself is a constructive
		column
			fixed support
	C 17h 43m	Lowering of column -1 st	• moderate noise
		step (10mm)

Table 2 :

 2 Overview of the estimated first 5 modes with natural frequencies f , their variation coef-

	ficient	/ 100

Table 3 :

 3 Overview of the first modes in the damaged state F (notation see Table2).

	mode	f (Hz)	ɶ σ	f	ξ (%)	ɶ ξ σ
	1	4.017	0.29	1.0	36
	2	5.891	0.27	1.2	25
	3a	9.347	0.28	0.97	29
	3b	10.89	0.33	2.1	29

mode 3b, 10.89 Hz, 2.1% mode 4, 14.25 Hz, 1.9% mode 5, 16.25 Hz, 2.7%

ACKNOWLEDGEMENTS

This work was partially supported by the European project FP7-PEOPLE-2009-IAPP 251515 ISMS. We also thank VCE for providing the data from S101 Bridge.


Michael Döhler
email: michael.doehler@bam.de

Laurent Mevel
email: laurent.mevel@inria.fr

Falk Hille
email: falk.hille@bam.de

Efficient Computation of Minmax Tests for Fault Isolation and Their Application to Structural Damage Localization

Keywords: Fault isolation, residual evaluation, statistical tests, numerical computation, mechanical systems

Fault detection and isolation can be handled by many different approaches. This paper builds upon a hypothesis test that checks whether the mean of a Gaussian random vector has become non-zero in the faulty state, based on a chi2 test. For fault isolation, it has to be decided which components in the parameter set of the Gaussian vector have changed, which is done by variants of the chi2 hypothesis test using the so-called sensitivity and minmax approaches. While only the sensitivity of the tested parameter component is taken into account in the sensitivity approach, the sensitivities of all parameters are used in the minmax approach, leading to better statistical properties at the expense of an increased computational burden. The computation of the respective test variable in the minmax test is cumbersome and may be illconditioned especially for large parameter sets, asking hence for a careful numerical evaluation. Furthermore, the fault isolation procedure requires the repetitive calculation of the test variable for each of the parameter components that are tested for a change, which may be a significant computational burden. In this paper, dealing with the minmax problem, we propose a new efficient computation for the test variables, which is based on a simultaneous QR decomposition for all parameters. Based on this scheme, we propose an efficient test computation for a large parameter set, leading to a decrease in the numerical complexity by one order of magnitude in the total number of parameters. Finally, we show how the minmax test is useful for structural damage localization, where an asymptotically Gaussian residual vector is computed from outputonly vibration data of a mechanical or a civil structure.

INTRODUCTION

Fault detection and isolation (FDI) are of practical importance in many branches of automatic control and have found a strong interest in connected areas [START_REF] Hwang | A survey of fault detection, isolation, and reconfiguration methods[END_REF]. In structural engineering, they have found an application for vibration monitoring of civil, aeronautical or mechanical structures in order to detect damages (fault detection) and to locate damages (fault isolation) [START_REF] Basseville | Subspace-based fault detection algorithms for vibration monitoring[END_REF][START_REF] Basseville | Statistical model-based damage detection and localization: subspace-based residuals and damage-to-noise sensitivity ratios[END_REF][START_REF] Döhler | Subspace-based fault detection robust to changes in the noise covariances[END_REF]. For such applications, FDI is based on a Gaussian residual vector with zero mean in the reference state and non-zero mean in the faulty state. Based on a parameterization of the investigated system, hypothesis tests are performed on the Gaussian vector for FDI [START_REF] Basseville | Information criteria for residual generation and fault detection and isolation[END_REF]Nikiforov, 1993, Basseville, 1997].

For fault isolation, it has to be decided which components in the parameter set of the Gaussian vector have changed, which can be done by variants of the hypothesis test using the so-called sensitivity or minmax approaches. While only the sensitivity of the tested parameter component is taken into account in the sensitivity approach, the sensitivities of all parameters are used in the minmax approach, leading to better statistical properties at the expense of an increased computational burden. Moreover, the estimated covariance of the Gaussian vector may be badly conditioned due to noise and insufficient data, which leads to an ill-conditioned computation of the respective test variable in the FDI tests. In [START_REF] Zhang | Advanced numerical computation of χ 2 -tests for fault detection and isolation[END_REF], a numerically robust computation of both statistical tests was presented. The robust computation is based on the pertinent use of the QR decomposition applied to the sensitivity and covariance parameters of the Gaussian vector. While significantly improving the computation of these tests and making them possible for practical applications, the computations can further be optimized. Furthermore, fault isolation often requires the repetitive calculation of the test variable for the different parameters that are tested for a change, in order to identify which of these parameters actually have changed, which may represent a significant computational burden. In this paper, we improve the numerical efficiency of the computation of the test variable in the minmax test from [START_REF] Zhang | Advanced numerical computation of χ 2 -tests for fault detection and isolation[END_REF] by reducing the number of necessary QR decompositions from three to one. Based on this development, we propose an efficient test computation for large parameter sets, leading to a decrease in the numerical complexity by one order of magnitude in the total number of parameters.

For structural vibration monitoring, output-only measurements are taken from a linear system subject to state and measurement noise [START_REF] Juang | Applied system identification[END_REF]. Subspace methods have shown excellent properties for the treatment of such data [START_REF] Van Overschee | Subspace Identification for Linear Systems: Theory, Implementation[END_REF][START_REF] Peeters | Reference-based stochastic subspace identification for output-only modal analysis[END_REF][START_REF] Basseville | Subspace-based algorithms for structural identification, damage detection, and sensor data fusion[END_REF][START_REF] Döhler | Fast multi-order computation of system matrices in subspace-based system identification[END_REF]]. An asymptotically Gaussian subspace-based residual vector for FDI has been proposed in [START_REF] Basseville | Subspace-based fault detection algorithms for vibration monitoring[END_REF], which is used for damage localization in [START_REF] Basseville | Statistical model-based damage detection and localization: subspace-based residuals and damage-to-noise sensitivity ratios[END_REF][START_REF] Balmès | Statistical model-based damage localization: a combined subspace-based and substructuring approach[END_REF], corresponding to fault isolation using a parameterization that is linked to local structural properties. However, in these works only the sensitivity test was used, whereas the minmax test has more appropriate properties. In this paper, we extend the previous results on structural damage localization with the minmax test. In this context, both the numerical robust computation and the efficient test computation for many parameters are of importance, the former due to noisy data and the latter since parameter sets can be huge in this application. This paper is organized as follows. In Section 2, the sensitivity and the minmax tests are recalled. In Section 3, the numerically robust computation of these tests is recalled from [START_REF] Zhang | Advanced numerical computation of χ 2 -tests for fault detection and isolation[END_REF]] and a more efficient and fast approach for multiple parameter tests is derived. These tests are shown in the context of structural damage localization in Section 4. In Section 5, an application on a simulated flexural beam is shown, where the minmax test outperforms the sensitivity test, before stating some concluding remarks in Section 6.

BASIC FAULT ISOLATION TESTS

We consider the problem of isolating changes in a parameter θ ∈ R l of a Gaussian residual vector ζ = ζ(θ) ∈ R h in a faulty state, following the distribution

ζ ∼ N (0, Σ) under H 0 N (J δθ, Σ) under H 1 (1)
where the null hypothesis H 0 corresponds to a reference state and H 1 to a faulty state. The vector δθ ∈ R l denotes the parameter change, the matrix J ∈ R h×l with full column rank denotes the sensitivity of the residual vector, and the positive definite matrix Σ ∈ R h×h is its covariance. Many methods have been proposed in the literature for fault detection (decide if there is a change in θ) and isolation (decide which elements of vector θ have changed) by the means of such a Gaussian residual vector ζ, see e.g. [START_REF] Basseville | Detection of Abrupt Changes -Theory and Applications[END_REF], [START_REF] Basseville | Information criteria for residual generation and fault detection and isolation[END_REF].

For fault isolation, different partitions of the vector δθ into two subvectors are considered, where each time a decision is made if the first subvector is zero or not. Without loss of generality, let this partition be

δθ = δθ a δθ b . (2)
Let F = J T Σ -1 J be the Fisher information matrix of the parameter θ contained in vector ζ, and let J and F be partitioned accordingly as

J = [J a J b ], F = F aa F ab F ba F bb = J T a Σ -1 J a J T a Σ -1 J b J T b Σ -1 J a J T b Σ -1 J b . (3)
In the following, two isolation tests are recalled as stated in [START_REF] Zhang | Advanced numerical computation of χ 2 -tests for fault detection and isolation[END_REF], where a decision is made between δθ a = 0 and δθ a = 0 in partition (2). Both tests will be applied to a subspace-based residual vector with property (1) to perform fault isolation in the context of vibration-based structural damage localization [START_REF] Basseville | Statistical model-based damage detection and localization: subspace-based residuals and damage-to-noise sensitivity ratios[END_REF][START_REF] Balmès | Statistical model-based damage localization: a combined subspace-based and substructuring approach[END_REF] in Section 4.

Sensitivity Test

The simplest possibility for testing changes δθ a is to assume δθ b = 0 and thus ζ ∼ N (J a δθ a , Σ). The corresponding Generalized Likelihood Ratio (GLR) test writes as [START_REF] Basseville | Information criteria for residual generation and fault detection and isolation[END_REF]

t sens = ζ T Σ -1 J a J T a Σ -1 J a -1 J T a Σ -1 ζ, (4
) which is called sensitivity test. The test variable t sens is χ 2 -distributed with dim(θ a ) degrees of freedom and noncentrality parameter δθ T a F aa δθ a , if δθ b = 0 is actually true. For a decision, the test variable is compared to a threshold.

Minmax Test

Instead of assuming δθ b = 0, the variable δθ b is replaced by its least favorable value for a decision about δθ a , as follows. Define the partial residuals as

ζ a def = J T a Σ -1 ζ, (5a)
ζ b def = J T b Σ -1 ζ, (5b)
the robust residual as [START_REF] Basseville | Information criteria for residual generation and fault detection and isolation[END_REF]

ζ * a def = ζ a -F ab F -1 bb ζ b (6) and F * a def = F aa -F ab F -1 bb F ba .
(7) Then, the mean of the robust residual ζ * a is sensitive to changes δθ a , but blind to δθ b , and it holds

ζ * a ∼ N (F * a δθ a , F * a )
. The corresponding GLR test writes as

t mm = ζ * T a F * -1 a ζ * a , (8)
which is called minmax test. The test variable t mm is χ 2 -distributed with dim(θ a ) degrees of freedom and noncentrality parameter δθ T a F * a δθ a , independently of δθ b . For a decision, the test variable is compared to a threshold. Note that the invertibility of all matrices in the computation is guaranteed, since J is assumed to have full column rank and Σ is positive definite.

NUMERICALLY EFFICIENT COMPUTATION

Both the sensitivity and the minmax test require a number of matrix inversions, which may be numerically critical due to possible ill-conditioning of the covariance matrix Σ. In [START_REF] Zhang | Advanced numerical computation of χ 2 -tests for fault detection and isolation[END_REF]] numerically more robust computations of both tests were suggested, making use of the decomposition

Σ -1 = (Σ -1/2 ) T Σ -1/2
and a number of QR decompositions. Note that an efficient computation of Σ -1/2 is detailed in [Döhler andMevel, 2011, Döhler et al., 2014], and let n r × h be the dimensions of Σ -1/2 that may be different from each other.

After recalling the robust computations from [START_REF] Zhang | Advanced numerical computation of χ 2 -tests for fault detection and isolation[END_REF] for completeness, we improve the computation of the minmax test in this section. While the computation in [START_REF] Zhang | Advanced numerical computation of χ 2 -tests for fault detection and isolation[END_REF] needs three QR decompositions, we suggest a more efficient computation that just needs one simultaneous QR decomposition. Based on this development, we furthermore propose a fast computation for multiple minmax tests, where different subvectors of δθ are tested.

3.1 Results from [START_REF] Zhang | Advanced numerical computation of χ 2 -tests for fault detection and isolation[END_REF] Sensitivity Test

Using the thin QR decomposition [START_REF] Golub | Matrix computations[END_REF] Σ -1/2 J a = QR, the sensitivity test in (4) writes as

t sens = α T α, where α = Q T Σ -1/2 ζ.
Here, the number of numerically critical matrix inversions and multiplications is limited to a minimum compared to a direct computation in (4).

Minmax Test Using the thin QR decompositions 8) writes as t mm = βT β, where β = QT c Σ -1/2 ζ. Again, the number of numerically critical operations is strongly reduced.

Σ -1/2 J a = Qa Ra , Σ -1/2 J b = Qb Rb , (I -Qb QT b ) Qa = Qc Rc , the minmax test in (

Minmax Test Revisited

In the following theorem, a more direct and more efficient computation than in the previous section is shown. Theorem 1. Let the thin QR decomposition of

Σ -1/2 [J b J a ] = [Q b Q a ] R bb R ba 0 R aa (9)
be given and partitioned accordingly. Then, the minmax test in (8) writes as

t mm = β T β, where β = Q T a Σ -1/2 ζ. Proof. From the QR decomposition (9) it follows Σ -1/2 J a = Q b R ba + Q a R aa , Σ -1/2 J b = Q b R bb , where Q T a Q a = I, Q T b Q b = I and Q T b Q a = 0, and thus in (3) F = F aa F ab F ba F bb = R T ba R ba + R T aa R aa R T ba R bb R T bb R ba R T bb R bb . (10)
The partial residuals in (5a) and (5b) write as

ζ a = R T ba Q T b Σ -1/2 ζ + R T aa Q T a Σ -1/2 ζ, ζ b = R T bb Q T b Σ -1/2
ζ, and the robust residual in (6) yields

ζ * a = ζ a -R T ba R bb (R T bb R bb ) -1 ζ b = R T ba Q T b Σ -1/2 ζ + R T aa Q T a Σ -1/2 ζ -R T ba R bb (R T bb R bb ) -1 R T bb Q T b Σ -1/2 ζ = R T aa Q T a Σ -1/2 ζ.
Similarly, substituting the elements of F * a in ( 7) with ( 10) yields

F * a = (R T ba R ba + R T aa R aa ) -R T ba R bb (R T bb R bb ) -1 R T bb R ba = R T aa R aa . Finally, replacing ζ *
a and F * a in the minmax test (8) yields

t mm = ζ T (Σ -1/2 ) T Q a R aa (R T aa R aa ) -1 R T aa Q T a Σ -1/2 ζ = ζ T (Σ -1/2 ) T Q a Q T a Σ -1/2 ζ, leading to the assertion.
This theorem makes the operation of the minmax test directly visible, when writing it equivalently as

t mm = βT β with β = Q a Q T a Σ -1/2 ζ. Note that Q a Q T a
defines an orthogonal projection. So in fact the squared norm of a vector β is computed, which is the orthogonal projection of the covariance-normalized residual into the subspace of the normalized "change directions" Σ -1/2 J a that are of interest and that are orthogonal to the normalized change directions Σ -1/2 J b . Like this the projection is blind to changes in the parameters θ b , while only possible changes in parameters θ a remain in the projected vector. Note that this means in practice that the sensitivities of all parameters in J need to be sufficiently pairwise distinct, otherwise the effect of a tested parameter in δθ a might be removed by another parameter in the projection. This also corresponds to the necessity of (numerical) invertibility of the R matrix in (9) and thus of R aa and R bb in the proof of Theorem 1, which is in theory given through the full column rank of J and the positive definiteness of Σ. However, if these matrices are too badly conditioned because some parameters in θ are too close, the minmax test might not perform as expected, which is even more the case in its original formulation in Section 2.2.

Compared to the computation proposed in [START_REF] Zhang | Advanced numerical computation of χ 2 -tests for fault detection and isolation[END_REF], Theorem 1 provides a more direct approach, where only one instead of three QR decompositions and only about half the numerical operations are necessary. Moreover, a fast iterative update of the QR decomposition for testing multiple subsets becomes possible thanks to the direct computation in the new approach, as described in the following section.

Minmax Tests on Multiple Parameter Subsets

When not only one subset θ a of parameter θ is tested for a change, but many different subsets, the QR decomposition in ( 9) is necessary for each tested parameter subset, which is a significant computational burden if the number of parameters l is high.

Without loss of generality, we assume θ a to be onedimensional in the following for simplicity of notation. Denote the elements of parameter θ = θ 1 θ 2 . . . θ l T . Thus, when testing parameter θ i for a change, denote θ a and θ b from the previous section as θ i a = θ i , (11a)

θ i b = θ i+1 • • • θ l θ 1 • • • θ i-1 T , (11b)
and the corresponding sensitivity matrices as J i a and J i b , respectively.

Instead of computing the QR decomposition in (9) for each parameter i, the fact can be used that the columns of the sensitivity matrix J are the same all the time, but only differently ordered in J i b J i a . In the following, an efficient iterative computation of the QR decomposition for the different parameters is derived.

Assume that the thin QR decomposition is given for parameter i in (9) as

K i def = Σ -1/2 J i b J i a = Q i R i , (12)
where

K i ∈ R nr×l . At this stage, the vector Q i a in Q i = Q i b Q i a
is used for the minmax test of parameter θ i a in Theorem 1 in the previous section.

In the next step, it is convenient to compute the QR decomposition for parameter i+1 to perform the respective minmax test, where the decomposition

K i+1 = Σ -1/2 J i+1 b J i+1 a = Q i+1 R i+1 .
(13) is required. It can be obtained efficiently as follows. The matrices K i and K i+1 are related by a permutation of one column due to the definition in (11a)-(11b): partitioning 12) and thus

K i = [k 1 K 2 ] with k 1 being the first column leads to K i+1 = [K 2 k 1 ]. Partitioning R i = [r 1 R 2 ] analogously, it follows K i = Q i r 1 Q i R 2 from (
K i+1 = Q i [R 2 r 1 ] .
The matrix [R 2 r 1 ] is upper Hessenberg due to the upper triangular structure of R i = [r 1 R 2 ]. Hence, there are Givens rotations G 1 , . . . , G l-1 , such that the product [START_REF] Golub | Matrix computations[END_REF], and the QR decomposition of the next iteration in (13) can be obtained efficiently from 14) This update for Q i+1 and R i+1 takes only 6n r l flops and 3l 2 flops, respectively, compared to 4n r l 2 -4 3 l 3 flops for the direct computation of the thin QR decomposition [START_REF] Golub | Matrix computations[END_REF].

G T l-1 • • • G T 1 [R 2 r 1 ] is upper triangular
Q i+1 = Q i G 1 • • • G l-1 , R i+1 = G T l-1 • • • G T 1 [R 2 r 1 ]. (
The QR decomposition is the dominant operation in the computation of the minmax test, and it holds n r ≥ l. Testing all parameters i = 1, . . . , l with the minmax test requires thus O(n r l 2 ) flops with the iterative scheme in ( 14), compared to O(n r l 3 ) flops for the direct computation.

VIBRATION-BASED STRUCTURAL DAMAGE LOCALIZATION

Based on the asymptotic local approach for change detection [START_REF] Benveniste | The asymptotic local approach to change detection and model validation[END_REF], a subspace-based residual vector has been proposed in [START_REF] Basseville | Subspace-based fault detection algorithms for vibration monitoring[END_REF] that asymptotically satisfies relation (1). It found an interesting application in vibration monitoring of civil, aeronautical or mechanical structures for damage detection and localization [START_REF] Basseville | Statistical model-based damage detection and localization: subspace-based residuals and damage-to-noise sensitivity ratios[END_REF][START_REF] Balmès | Statistical model-based damage localization: a combined subspace-based and substructuring approach[END_REF]. For the localization problem, a structural parameterization vector θ in connection with a finite element model (FEM) of the structure is used, whose components refer to local structural properties. Damage localization thus corresponds to isolating the changed components of θ by means of the residual vector.

In this section, the subspace-based damage localization method is recalled from [START_REF] Basseville | Statistical model-based damage detection and localization: subspace-based residuals and damage-to-noise sensitivity ratios[END_REF][START_REF] Balmès | Statistical model-based damage localization: a combined subspace-based and substructuring approach[END_REF], where, however, only the sensitivity test was used for fault isolation. We extend these results with the minmax test in this section.

The Residual Vector

The behavior of linear time-invariant dynamical structures subject to unknown ambient excitation can be described by system

M Ẍ (t) + C Ẋ (t) + KX (t) = υ(t) ( 15
)
where t denotes continuous time; M, C, K ∈ R m×m are mass, damping, and stiffness matrices, respectively; the (high dimensional) state vector X (t) ∈ R m is the displacement vector of the m degrees of freedom of the structure; and υ(t) is the external unmeasured force (noise). Observing system (15) at r sensor positions (e.g. acceleration or displacement sensors) at discrete time instants t = kτ (with sampling rate 1/τ ), it can be transformed to an equivalent discrete-time state space system [START_REF] Juang | Applied system identification[END_REF]

x k+1 = Ax k + v k y k = Cx k + w k (16)
with the states x k ∈ R n , the measured outputs y k ∈ R r , the state transition matrix A ∈ R n×n and the observation matrix C ∈ R r×n , where n = 2m is the model order and r the number of outputs. The state noise v k and output noise w k are unmeasured and assumed to be Gaussian, zero-mean, white. The eigenvalues and eigenvectors at the observed degrees of freedom of model ( 15) constitute a canonical parameterization that can be equivalently obtained from the eigenstructure of model ( 16) [START_REF] Juang | Applied system identification[END_REF].

The modal parameters (frequencies, damping ratios and mode shapes) are deduced from the eigenstructure, and their changes indicate damages (faults). While fault detection can thus be built on changes in the eigenstructure, fault isolation for damage localization is built on a parameterization containing local structural parameters later on.

First, a subspace-based residual vector is defined as follows [START_REF] Basseville | Subspace-based fault detection algorithms for vibration monitoring[END_REF]. Let G = E(x k+1 y T k ) be the cross-covariance between the states and the outputs, let

R i = E(y k y T k-i ) = CA i-1 G (i ≥ 1
) be the theoretic output covariances and let

H p+1,q def =     R 1 R 2 . . . R q R 2 R 3 . . . R q+1 . . . . . . . . . . . . R p+1 R p+2 . . . R p+q     def = Hank (R i ) (17)
be the theoretic block Hankel matrix with the well-known factorization property

H p+1,q = O p+1 C q (18)
into observability and controllability matrices

O p+1 =     C CA . . . CA p     , C q = G AG . . . A q-1 G .
Parameters p and q are chosen such that rank(O p ) = rank(C q ) = n. In the reference state, the left null space S of observability matrix O p+1 is obtained, such that S T O p+1 = 0. Then, S defines also a basis of the left null space of H p+1,q in the reference state because of factorization property (18), and the characteristic property of a system in the reference state writes [START_REF] Basseville | Subspace-based fault detection algorithms for vibration monitoring[END_REF] S T H p+1,q = 0.

(19) From the outputs {y k : k = 1, . . . , N }, a consistent estimate H p+1,q = Hank R i is obtained from the estimated output covariances R i = 1/N N k=1 y k y T k-i . Based on property (19), the residual vector

ζ N = √ N vec S T H p+1,q (20)
is defined, which has the property E(ζ N ) = 0 iff the system is in the reference state.

Fault Isolation Tests for Damage Localization

The evaluation of the residual is possible based on the asymptotic local approach for change detection [START_REF] Benveniste | The asymptotic local approach to change detection and model validation[END_REF], assuming the close hypotheses H 0 : θ = θ 0 (reference system),

H 1 : θ = θ 0 + δθ/ √ N (faulty system), (21)
where vector θ 0 denotes the system parameter in the reference state and vector δθ is unknown but fixed. With this statistical framework, very small changes in the system parameter θ can be detected if N is large enough. The subspace-based residual ζ N yields [START_REF] Basseville | Subspace-based fault detection algorithms for vibration monitoring[END_REF]]

ζ N d -→ N (0, Σ) under H 0 N (J δθ, Σ) under H 1 (22)
for N → ∞, where J and Σ are the asymptotic sensitivity and covariance, respectively. The estimation of J is described in detail in [START_REF] Basseville | Statistical model-based damage detection and localization: subspace-based residuals and damage-to-noise sensitivity ratios[END_REF][START_REF] Balmès | Statistical model-based damage localization: a combined subspace-based and substructuring approach[END_REF] for a chosen structural parameterization θ, and the estimation of Σ is described in [START_REF] Döhler | Robust subspace based fault detection[END_REF], both in the reference state.

Thanks to property ( 22), both the sensitivity and the minmax test can be applied to residual ζ N (asymptotically, N large enough) for each of the structural parameters in θ = θ 1 . . . θ l . A test value for a parameter above a certain threshold corresponds to a damage location in the structure associated with the parameter. In this respect, the preferred test should be the minmax test instead of the sensitivity test used in [START_REF] Basseville | Statistical model-based damage detection and localization: subspace-based residuals and damage-to-noise sensitivity ratios[END_REF][START_REF] Balmès | Statistical model-based damage localization: a combined subspace-based and substructuring approach[END_REF], since the condition δθ b = 0 can never be guaranteed when testing δθ a = δθ i for some parameter i. As all parameters of θ are tested successively, the efficient test computation of Section 3.3 is beneficial especially for large parameterizations.

As noted in Section 3, care must be taken that J is wellconditioned for the minmax test, i.e. that the structural parameters in θ are sufficiently distinct. This might require a prior clustering step as already performed in [START_REF] Basseville | Statistical model-based damage detection and localization: subspace-based residuals and damage-to-noise sensitivity ratios[END_REF][START_REF] Balmès | Statistical model-based damage localization: a combined subspace-based and substructuring approach[END_REF].

APPLICATION

The model of a flexural beam structure is considered for the application of the damage localization tests. The model comprises 40 equally sized finite 2D elastic beam elements (Figure 1) allowing two translational and one rotational displacements at each node. The beam has a fixed support on one end and a movable support on the other end, so the system is statically indeterminate and four of the altogether 123 degrees of freedom of the model are restricted. The beam has invariable cross-sectional and material properties in the undamaged state. The mass is distributed and defined such that the first frequency of the system is at 5 Hz and damping is defined such that all modes have a damping ratio of 2%.

Damage is modeled by a 20% reduction of Young's modulus for two damage cases, first in element 10, and second both in elements 10 and 18, having a linear impact on the stiffness at the affected elements. The parameter vector θ is of dimension 40 and consists of Young's modulus at each of the 40 elements.

For both the reference and the damaged states, a data set with N = 50 000 acceleration samples at ten sensor positions was generated at a sampling frequency of 500 Hz from white noise excitation. 5% white noise were added on the generated outputs as measurement noise. All parameters of the localization tests (S, J , Σ) were estimated based on the data set in the reference state and information from the FEM using only the first eight modes. In the damaged state, both the sensitivity and the minmax tests with their efficient numerical computation from Section 3 were applied to the data sets of the two fault scenarios. No comparison was made to the naive computation of the test variables (Sections 2.1 and 2.2), since the ineffectiveness of this computation was already shown in [START_REF] Zhang | Advanced numerical computation of χ 2 -tests for fault detection and isolation[END_REF].

The respective test values were computed for each of the 40 structural elements and are shown in Figures 2 and3. In both fault scenarios, the damaged elements have the highest test values, thus correctly indicating the damage locations. However, the test values corresponding to the undamaged elements in the sensitivity test are much higher (in relative terms) than in the minmax test. This behavior was expected since the sensitivity test value for some parameter may also be influenced by a change in another parameter. Thus, the minmax test indeed offers a better "contrast" for fault isolation in our application on damage localization and should be the statistical test of choice for this application.

CONCLUSIONS

In this paper, two statistical tests for fault isolation based on a Gaussian residual vector were recalled. A new, numerically more efficient computation for the minmax test was derived in Theorem 1 together with an associated fast computation when testing the different parameters of a large parameter set in Section 3.3, which are the main theoretical results of this paper. An application of this test was shown for vibration-based structural damage localization, where changes in local parameters of a structure indicate damage. In this setting, the minmax test proved to be more effective than the previously used sensitivity test and seems to be promising for applications in structural health monitoring.

  This work was supported by the European project FP7-PEOPLE-2009-IAPP 251515 ISMS.

  Fig. 1. Flexural beam model with element numbers and ten sensors.

Fig. 2 .Fig. 3 .

 23 Fig. 2. Sensitivity tests (left) and minmax tests (right) for damage at element 10.