Keywords: Model updating, Damage quantification, Uncertainty quantification, subspace methods

The purpose of model updating is to minimize the misfit between the structural response measurements and an assumed numerical model. In the context of damage quantification, this misfit is characterized by some features computed from the response data measured on the faulty structure, and its Finite Element (FE) model in the healthy condition. The FE model is parameterized so that the estimated features are related to the physical parameters of the model. Therein, the parameterization size may be large. As a consequence of low instrumentation, different parameters can have a similar effect on the estimated features, resulting in non uniqueness of the updating problem solutions, even taking into account the inherent uncertainty errors, originating both from the model and the measurements. In this paper a model updating-based damage quantification strategy is proposed. It involves the minimization of two Hankel matrices, one computed from the data and another from the optimized model. The difference between those two matrices is studied, in particular in the practical case where the ambient excitation is unknown. It yields a statistical residual, whose deviations from zero can be evaluated through a statistical test. The resulting optimization is based on the Generalized Likelihood Ratio test as an objective function and uses its 95 per cent quantile as a measure of closeness for a stopping criterion for the optimization. Moreover, the large size of the finite element model to optimize compared to the low instrumentation has to be taken into account by clustering the parameter space. This clustering is proceeded using the well known stochastic subspace-based damage localisation method. The proposed framework is validated on simulations of a simple mechanical system.

Statistical optimization for subspace-based damage quantification

Szymon Greś * Michael Döhler * * Laurent Mevel * *

INTRODUCTION

The damage diagnosis problem can be divided into damage detection, localization, quantification and lifetime prediction [START_REF] Rytter | Vibrational based inspection of civil engineering structures[END_REF]. Whereas the detection of damages from vibration measurements is well established e.g. in [START_REF] Döhler | Subspace-based damage detection under changes in the ambient excitation statistics[END_REF][START_REF] Ulriksen | Operational modal analysis and wavelet transformation for damage identification in wind turbine blades[END_REF][START_REF] Worden | Damage detection using outlier analysis[END_REF][START_REF] Bernal | Kalman filter damage detection in the presence of changing process and measurement noise[END_REF][START_REF] Yan | Null subspace-based damage detection of structures using vibration measurements[END_REF][START_REF] Gres | A comparison damage detection methods applied to civil engineering structures[END_REF], the localization and quantification of damages is more complex, and requires additional physical information on the examined structure [START_REF] Allahdadian | Towards robust statistical damage localization via model-based sensitivity clustering[END_REF][START_REF] Bernal | Subspace exclusion zones for damage localization[END_REF][START_REF] Ulriksen | Shaped input distributions for structural damage localization[END_REF][START_REF] Ulriksen | Structural damage localization by outlier analysis of signal-processed mode shapes -analytical and experimental validation[END_REF]. In this respect, few complete frameworks for damage identification exist, e.g. [START_REF] Döhler | Fault detection, isolation and quantification from Gaussian residuals with application to structural damage diagnosis[END_REF][START_REF] Hà | Localization and quantification of damage in beam-like structures using sensitivities of principal component analysis results[END_REF][START_REF] Yan | Structural damage localization by combining flexibility and stiffness methods[END_REF]. Some methods consider a specific structural type in their theory, e.g. beams [START_REF] Hà | Localization and quantification of damage in beam-like structures using sensitivities of principal component analysis results[END_REF], while other methods incorporate the required physical information from a FE model, which allows application to more complex structures. In this context, the sensitivities of a damage feature computed from the structural responses can be used to infer FE model parameter changes in [START_REF] Döhler | Fault detection, isolation and quantification from Gaussian residuals with application to structural damage diagnosis[END_REF], under the assumption of small damage extents. Contrary to these approaches, FE model optimization-based methods are usually not limited by the structural type or by the damage extent. Therein model updating is a generic term encompassing a family of many methods [START_REF] Mottershead | Model updating in structural dynamics: A survey[END_REF][START_REF] Friswell | Finite element model updating in structural dynamics[END_REF], that are often applied in the damage quantification context [START_REF] Behmanesh | Probabilistic damage identification of a designed 9-story building using modal data in the presence of modeling errors[END_REF][START_REF] Huth | Damage identification using modal data: Experiences on a prestressed concrete bridge[END_REF][START_REF] Reynders | Damage identification on the Tilff Bridge by vibration monitoring using optical fiber strain sensors[END_REF][START_REF] Teughels | Damage detection and parameter identification by finite element model updating[END_REF]. However, model updating methods are often poorly conditioned due to the possibly large FE parameter space in comparison to relatively few parameters that can be extracted from data, and are prone to statistical uncertainty errors of the estimated features. The objective of this paper is to develop a cascade framework where the damage is first localized and then quantified on the subset of the damaged parameter space, and which explicitly considers data-based uncertainties in the model optimization. An objective function is proposed based on the statistical hypothesis test of a residual composed from the Hankel matrices of output covariances of the monitored system compared to a numerical counterpart. The proposed residual unites simplicity of the damage feature, robustness towards changing ambient excitation conditions by an appropriate normalization, and a sound statistical evaluation within the local approach framework.

The paper is organized as follows. The considered model optimization problem, the background on system modeling is stated in Section 2. The objective function based on Hankel matrix difference is defined in Section 3. The damage quantification framework based on the damage localization-based parameter clustering and the stochastic optimization scheme is presented in Section 4. The application of the proposed framework on a numerical simulation of a simple mechanical system is enclosed in Section 5.

BACKGROUND

The overall goal of this paper is the estimation of the structural parameter change between a healthy reference state and the current damaged state. Let θ ∈ Θ ⊂ R p be the parameter vector that contains the damage-sensitive parameters of the structural elements of interest for the considered problem within a bounded parameter space Θ. The parametrization is user-defined and adapted to the specific monitoring problem at hand. It is assumed that the parameter vector θ 0 in the healthy state of the structure is known, and let θ * be the parameter vector in the damaged state. It is the goal to obtain an estimate θ of θ * from vibration measurements in the damaged states for estimating the damage extent δ = θ * -θ 0 . Finding the optimal θ can be formulated as an optimization problem, where an objective function F (θ), also called cost function, is minimized over the parameter space Θ.

The objective function is designed to represent the discrepancy between v, the estimate of a feature vector computed from measurement data recorded under the (unknown) system parameter θ * , and its counterpart v(θ) computed from a parametric model. The optimal solution for θ is then obtained as θ = arg min θ∈Θ F (θ).

(1)

A classic feature vector v for the design of the objective function F (θ) in ( 1) is based on the modal parameters of the monitored system [START_REF] Simoen | Dealing with uncertainty in model updating for damage assessment: A review[END_REF], and their statistical uncertainties obtained using e.g. [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] Mellinger | Variance estimation of modal parameters from output-only and input/output subspace-based system identification[END_REF]. In this work an objective function is proposed based on the statistical hypothesis test of a residual composed from the Hankel matrices of the monitored system compared to a numerical counterpart. An immediate advantage is to avoid the computation and the correspondence between the numerical modes and the estimated ones.

A decision whether the current parameter of the numerical model θ corresponds to the unknown parameter θ * needs to take into account the uncertainty information of the estimated feature and its sensitivity towards the considered parameterization. An approach to assess this equivalence is to define a residual whose properties can be evaluated properly as the number of samples goes to infinity. For this, two hypotheses are defined

H 0 : θ = θ * (model matched), (2)
H 1 : θ = θ * + κ/ √ N (model mismatched),
where κ = √ N (θθ * ) is a change vector between the assumed system parameter θ and the the sought parameter θ * . Notice that the hypotheses H 0 and H 1 are reversed from the classical damage detection setup in [START_REF] Allahdadian | Towards robust statistical damage localization via model-based sensitivity clustering[END_REF], where H 0 represented the healthy reference state, whereas here it represents the damaged state, which is now a reference for the model updating.

System and models

Assume that the dynamics of the monitored system can be modeled as linear time-invariant (LTI) with d degrees of freedom (DOF), which are described by the differential equation of motion M θ q(t) + D θ q(t) + K θ q(t) = u(t)

(3) where t denotes the continuous time, and the matrices M θ , D θ and K θ ∈ R d×d denote the mass, damping and stiffness matrices, respectively, which depend on the system parameter θ. The vectors q(t) ∈ R d and u(t) ∈ R d denote the continuous-time displacements and the unknown external forces, respectively. Observed at r sensor positions, e.g. by acceleration, velocity or displacement sensors at discrete time instants t = kτ (with sampling rate 1/τ ), system (3) can be transformed into the discrete-time state-space model [START_REF] Juang | Applied System Identification[END_REF] x

θ k+1 = A θ x θ k + w k y θ k = C θ x θ k + v k (4)
where x θ k ∈ R n are the states, y θ k ∈ R r are the outputs, vectors w k and v k denote the white process and output noise respectively, A θ ∈ R n×n , C θ ∈ R r×n are the state and observation matrices, and n = 2m is the model order. The process noise v k is assumed to be a stationary process with zero mean and covariance matrix Q = E(v k v T k ), w k denotes the zero-mean output noise with covariance matrix R = E(w k w T k ), and the covariance between v k and

w k is S = E(v k w T k )
, where E(•) denotes the expectation operator.

For simplicity, the (•) θ notation is dropped in the remainder of this paragraph. Let R i = E(y k y T k-i ) = CA i-1 G be the theoretical output covariances of the measurements, where

G = E(x k+1 y T k ) = AΣ s C T +S and Σ s = E(x k x T k ) = AΣ s A T + Q. The collection of R i can be stacked to form a block Hankel matrix H =     R 1 R 2 . . . R q R 2 R 3 . . . R q+1 . . . . . . . . . . . . R p+1 R p+2 . . . R p+q     ∈ R (p+1)r×qr , (5)
where p and q are chosen such that min(pr, qr) ≥ n with often p+1 = q. Matrix H enjoys the factorization property

H = O(C, A) C(A, G), (6)
where the observability and controllability matrices O(C, A) and C(A, G) are defined as

O(C, A) =     C CA . . . CA p     , C(A, G) = G AG . . . A q-1 G , (7)
where (A, C) can be easily computed directly from M, D, K [START_REF] Juang | Applied System Identification[END_REF], and G is obtained based on the chosen noise properties after [START_REF] Van Overschee | Subspace Identification for Linear Systems[END_REF].

Consistent estimates Ĥ can be obtained from the output covariances of the measurements {y k } k=1,...,N +p+q e.g. by

Ĥ = Y + Y -T , (8)
where the data Hankel matrices Y + and Y -contain future and past time horizons

Y + = 1 √ N     y q+1 y q+2 . . . y N +q y q+2 y q+3 . . . y N +q+1 . . . . . . . . . . . . y p+q+1 y p+q+2 . . . y p+q+N     , Y -= 1 √ N     y q y q+1 . . . y N +q-1 y q-1 y q . . . y N +q-2 . . . . . . . . . . . . y 1 y 2 . . . y N     .
The estimates of O(C, A) and C(A, G) are classically obtained from a singular value decomposition (SVD) of Ĥ thanks to the factorization property (6).

Hankel matrix normalization scheme

Model-based Hankel matrix parametrized with θ, and Hankel matrix computed from data generated under parameter θ * can be compared by a simple objective function for a model optimization-based damage quantification. Hankel matrices, however, are not only dependent on the dynamics of the system, but also on the unknown covariance of the noise processes Q. Then, since the modelbased Hankel matrix can not be defined under the same Q as its data-based equivalent, their comparison in search of θ is meaningless. To overcome this problem, a proper normalization of the model-based Hankel matrix is recalled after [START_REF] Gres | Hankel matrix normalization for robust damage detection[END_REF], where this normalization was introduced to compensate change in the noise properties between different data sets. Here, since there is only one data set, the matrix Q is not changing but is still unknown, and the same normalisation is useful.

Let H θ * data and H θ * model be two exact Hankel matrices of rank n for a system in the state θ * , subjected to process noise with different covariances Q data and Q model . A SVD of the juxtaposed matrices H data and H model writes

H θ * data H θ * model = [U s U ker ] D s 0 0 0 V T s V T ker , (9)
where rank

H θ * data H θ * model = n, U s ∈ R (p+1
)r×n contains the left singular vectors, D s ∈ R n×n contains the non-zero singular values and V s ∈ R 2qr×n contains the right singular vectors, which are split into

V T s = V T s,data V T s,model corresponding to H θ * data and H θ * model re- spectively. Now define Z data = D s V T s,data , Z model = D s V T s,model , (10)
where both Z data and Z model are full row rank. The exact Hankel matrices share the same image in the reference state

H θ * data H θ * model = U s [Z data Z model ] . (11)
To compare H θ * data with H θ * model an appropriate normalization is given by

H model = H θ * model Z † model Z data , (12)
where H model now shares the same C(A, G) as H θ * data .

Hankel matrix difference-based residual for damage diagnosis

Let H θ model be the model-based Hankel matrix generated under some chosen process noise covariance Q model . Next, let Ĥθ * data be obtained from a data set of length N generated under a process noise covariance Q data , with Q data different from Q model . Then, after [START_REF] Ulriksen | Structural damage localization by outlier analysis of signal-processed mode shapes -analytical and experimental validation[END_REF] it yields

H θ model Z † model Z data -H θ * data = 0 iff θ = θ * , (13)
and

H θ model Z † model Z data -H θ * data = 0 iff θ = θ * , (14)
where Z data and Z model are obtained from an SVD of [H θ * data H θ model ] as in ( 9) but truncated at order n. Based on both ( 13) and ( 14), the change detection residual is defined as

ζθ def = √ N vec(H θ model Ẑ † model Ẑdata -Ĥθ * data ), (15)
where Ẑdata and Ẑmodel are defined from the following SVD and partitioned at order n 15) is zero if and only if the parameters θ and θ * correspond, up to some statistical considerations that will be discussed in the next section.

Ĥθ * data H θ model = Û s Û ker Ds 0 0 Dker V T s,data V T s,model V T ker,data V T ker,model , as Ẑdata = Ds V T s,data , Ẑmodel = Ds V T s,model . Then, (
Alternatively or concurrently to the difference of Hankel matrices, other metrics can be used, based on modal parameters or the kernel of the Hankel matrix used for damage localization in [START_REF] Allahdadian | Towards robust statistical damage localization via model-based sensitivity clustering[END_REF]. The present metric has the merits to relate to the classical Manhalobis distance [START_REF] Gres | Hankel matrix normalization for robust damage detection[END_REF].

OBJECTIVE FUNCTION BASED ON HANKEL MATRIX DIFFERENCE

In this section an objective function for model optimization is expressed by a difference between two Hankel matrices. The statistical properties of the residual defined in Section 2.2 are derived and a hypothesis test for datamodel correspondence is established. Then, an objective function for optimization of the numerical model is constructed. Finally, a criterion for assessing the minima of the objective function is derived based on the quantile of the developed statistical distribution.

Residual distribution

To characterize the distribution of the residual (15) the asymptotic local approach for change detection [START_REF] Benveniste | The asymptotic local approach to change detection and model validation[END_REF] is used. Firstly, thanks to the Central Limit Theorem (CLT), the local approach ensures that Ĥθ data is asymptotically Gaussian, and it holds

√ N vec( Ĥθ * data -H θ * data ) L -→ N (0, Σ data ), (16)
where Σ data is the asymptotic reference Hankel matrix covariance. Note that an estimate of Σ data can be easily obtained from the sample covariance as in e.g. [START_REF] Döhler | Subspace-based fault detection robust to changes in the noise covariances[END_REF].

Then, it can be proved that the residual in ( 15) is asymptotically Gaussian with the following properties

H 0 : ζθ * L -→ N (0, Σ ζ ), (17)
H 1 : ζθ L -→ N (J ζ θ * κ, Σ ζ ) , (18) where
Σ ζ = J ζ H data Σ data (J ζ H data ) T , (19)
and

J ζ H data = I (p+1)r ⊗ U ker U T ker .
The sensitivity of the residual with respect to the chosen parameterization J ζ θ * is defined as

J ζ θ * = (Z † data Z data ) T ⊗ U ker U T ker J H data θ * , (20)
where

J H data θ * = ∂vec (H data )/∂θ (θ * ).
Notice that the sensitivity matrix J H data θ *

is not computed from the unknown model θ * , but from the data alone. In that sense, the sensitivities are expressed with respect to a data-driven parameterization, which consists in the modal parameters. As such, a consistent estimate of J H data θ * can be computed using the modal parameter estimates, as in [START_REF] Balmès | Statistical model-based damage localization: A combined subspace-based and substructuring approach[END_REF].

Objective function for model optimization

Let J and Σ respectively be consistent estimates of J ζ θ * and Σ ζ . Then a Generalized Likelihood Ratio (GLR) test to decide between H 0 and H 1 writes as

GLR θ = ( ζθ ) T Σ -1 J J T Σ -1 J -1 J T Σ -1 ζθ . (21)
Assuming Σ ζ to be invertible, under H 0 , GLR θ follows a χ 2 distribution with rank(J ζ θ * ) degrees of freedom. The test value from ( 21) can be directly used in the objective function for the model optimization as

F (θ) = GLR θ . ( 22
)
The quantile q χ 2 of the underlying χ 2 distribution, satisfying

q χ 2 0 f χ 2 (x)dx = γ,
where γ is a desired confidence level, can be used to define an acceptance region

Θ = {θ : GLR θ ≤ q χ 2 }, (23)
which comprises all parameter vectors θ that yield H θ model to comply with the estimated reference Ĥθ * data under the statistical distribution of the GLR. Thus, Θ comprises the statistically acceptable solutions for θ * with regards to the considered models, and a stopping criterion of the optimization search can be formulated for θ ∈ Θ.

Notice that the stopping criterion and the objective function do not involve computation or matching of the numerical modes for any model. Only the modes identified from the data set are needed in the computation of J . It yields a very direct and simple computational scheme for optimization.

DAMAGE QUANTIFICATION STRATEGY

In the previous section, an objective function has been elaborated whose set of minima corresponds to models that are statistically compliant with the Hankel matrix computed from measurements. As such, looking for the minima of this function will yield the sought models. This, however, is not trivial.

For example, some optimization strategies might be stuck in local minima due to the rugged nature of the proposed objective function. Moreover, the lack of identifiability for parameters of large FE models requires clustering in order to estimate the change over subsets of similar parameters, and consequently to quantify the possible damage. In this paper, it is proposed to chain a damage localization approach with a model optimization procedure to minimize the designed objective function [START_REF] Reynders | Uncertainty bounds on modal parameters obtained from stochastic subspace identification[END_REF] over the subset of parameters that have been recognized as damaged by the localization approach. In this way, the search domain for the optimization approach can be reduced efficiently. The resultant damage quantification framework is summarized in Figure 1.

For damage localization, the subspace-based damage localization approach [START_REF] Allahdadian | Towards robust statistical damage localization via model-based sensitivity clustering[END_REF] is chosen, and the Covariance Matrix Adaptation Evolution Strategy (CMA-ES) [START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF] is adapted as the optimization approach. Both methods are outlined in the following subsections. Similarly, other localization methods and other optimization methods can be chosen for this framework. The damage quantification procedure comprising the aforementioned strategies is summarized in Algorithm 1. Due to a possibly large FE model-based parametrization θ on one side and limited measurements regarding the number of sensors on the other side, the sensitivity of the optimized feature vector with respect to some components of θ may be equal or be very close. Thus, such parameter components are indistinguishable, and clustering of the parameters is performed to estimate their total change over some localized subset.

For this purpose the statistical sensitivity clustering, which is a part of the subspace-based damage localization [START_REF] Balmès | Statistical model-based damage localization: A combined subspace-based and substructuring approach[END_REF][START_REF] Allahdadian | Towards robust statistical damage localization via model-based sensitivity clustering[END_REF][START_REF] Mendler | Clustering of redundant parameters for fault isolation with Gaussian residuals[END_REF] is used. The output of such a procedure is a subset of the parameter vector θ, denoted by θ loc , which contains all parameters of the clusters that are recognized as damaged. Hence, it can be considered that the remaining parameters in θ are unchanged. Subsequently, the damage extent can be quantified by updating only the parameter subset θ loc . Thus, linking the damage localization with an adequate model optimization method should lead to a damage quantification scheme, where the dimension of the search domain is reduced by the localization approach.

Clustering the parameter space is essential for large parameterization. However, this will not be detailed in the considered numerical applications for proof of concept where all elements are well separated statistically.

Model optimization with modified CMA-ES

Starting with the initial value θ = θ init corresponding to the model in the reference state, the CMA algorithm consists in generating λ model candidates θ g j , j = 1, . . . , λ, in each population g, by sampling a multivariate Gaussian distribution. The model candidates are only modified on the parameter subset θ loc that is given by the localization approach. The sampling is carried out on the considered parameter subset for the subsequent population g + 1 as (θ loc j ) g+1 = m g + ε j , ε j ∼ σ g N (0, C g ) , (24) where j = 1, . . . , λ and m g is a weighted mean of the model candidates (θ loc j ) g in the parent generation. Then, the parameter subset of the full parameter vector θ g+1 j is updated with (θ loc j ) g+1 . Equation ( 24) represents a mutation and recombination into offsprings, for which the CMA-ES algorithm adapts the parameters C g and σ g in each generation. The covariance matrix C g of the added Gaussian noise represents the amplitude for the sampling to occur, and the scaling factor σ g determines the range of the considered mutation. Consequently, the optimization continues and the best parent solutions replace the offspring until it converges to a solution. For CMA-ES, the covariance matrix C g is incrementally updated with rank-one matrices representing the direction between the best parent solutions at two consecutive generations, such that the likelihood of previously successful search steps is increased [START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF].

For the convergence to a solution, a stopping criterion is included in the algorithm that is adapted to the acceptance region [START_REF] Döhler | Efficient multi-order uncertainty computation for stochastic subspace identification[END_REF]. Since the acceptance region consists of all models that are equally optimal in the statistical sense, the algorithm can stop once a number t opt of population model candidates are inside the region, where t opt ≥ 1. A higher value leads to more confidence in the set of retained solutions. Once inside the acceptance region, there is no need to further minimize the objective function, avoiding unnecessary additional computations.

Then the final solution is computed as the mean of the selected model candidates θ g j k ∈ Θ in the last population g that are lying in the acceptance region, with

θ sol = 1 topt topt k=1 θ g j k (25)
Finally, the change in the parameter vector is evaluated for damage quantification. Since the changes in the parameter components of each cluster are indistinguishable as described in the previous section, only the global change for each cluster can be evaluated, while the separate evaluation of changes of different parameters within the same cluster is impossible. Then, the change in any damaged cluster c can be evaluated as

δc = pc k=1 θ ic(k) sol -θ ic(k) init (26)
where i c (1), . . . , i c (p c ) are the indices of the components of parameter vector θ that correspond to cluster c, and p c is the number of elements in c.

NUMERICAL APPLICATION

In this section, the proposed damage quantification framework is applied in a numerical experiment on a 6 DOF mechanical chain-like system that for any consistent set of units is modeled with spring stiffness k 1 = k 3 = k 5 = 100 and k 2 = k 4 = k 6 = 200, mass m i = 1/20 and a proportional damping matrix such that all modes have a damping ratio of 3%. The chain is illustrated in Figure 2. The system is excited by a white noise acting at all DOFs whose covariance yields Q = bb T where b ∈ R 6×6 is a matrix whose entries are drawn from a standard normal distribution. The structural accelerations are simulated at DOFs 1, 3 and 5 at a sampling frequency of 50 Hz, and white measurement noise with 5% of the standard deviation of the output is added to each response measurement. Each measurement comprises N = 200,000 samples. In each data set the damage is modeled as a stiffness reduction of the 4-th spring by 5%. Hereafter let θ be the collection of the element stiffness k 1-6 . Albeit the damages considered herein relate to the stiffness changes, the proposed approach is general and can be formulated to quantify changes in any parameter of the numerical model for which the modal parameter sensitivities are nonzero. For damage quantification the procedure outlined in Algorithm 1 is used.

Due to the simplicity of the examined system, each parameter is reported in a separate cluster, indicating that the sensitivities of the estimated residual w.r.t. the considered system parameters are distinguishable. Moreover, a large probability for damage in the 4-th spring is reported by the subspace-based damage localization method, correctly mapping the location of damage. Since the parameterization is small enough, the full parameter set can be used for the damage quantification with Algorithm 1. In practice, however, the optimization would be only run on the detected clusters. The parameters used to initialize the optimization algorithm are depicted in Table 1.

Firstly, the proposed objective function is examined. The function F (θ) in [START_REF] Reynders | Uncertainty bounds on modal parameters obtained from stochastic subspace identification[END_REF] It can be viewed that F (θ) is convex in the considered parameter region. The global minimum of F (θ) is wellaligned with the value F (θ * ) at the exact parameter pair (100, 190) for the parameters k 3 and k 4 and (190, 100) for the parameters k 4 and k 5 . A group of parameters around the true value satisfies the uncertainty-based stopping criterion. This is a consequence of the statistical uncertainty introduced by estimating H θ * data in F (θ). Secondly, the optimization results based on a single measurement set are analyzed. Each optimization population contains λ parameters, from which t opt parameters are required to be within the acceptance region for the algorithm to stop. As such, the estimated change can be examined for every parameter set in each iteration of the algorithm, illustrating its performance. This is depicted in Figure 5.

Figure 5 shows the iterations of the estimated change δk3 and δk4 until the last iteration. Notice that within the last iteration not all the parameters are compliant with the stopping criterion. Those parameters are characterized with a higher objective function value, which indicates that at least one considered model is not corresponding to measurements. The model optimization takes 12 iterations to converge to a feasible solution, which together with the acceptance region is visualized in Figure 6. A zoom on δk3 and δk4 illustrated in Figure 6, corresponds to the parameters complying with the stopping criterion. The estimates are spread around the exact value, whereas the mean solution is close to the exact value.

Next, consider estimating δk for each of the 1000 measurement realizations. Each data set yields different estimates Ĥθ * data governed by the same statistical distribution. Consequently, the variability of each δk stems from the statistical uncertainty of Ĥθ * data . Thus computing δk from independent data sets should yield a histogram whose variance corresponds to the variance of the estimate of δ k . The histograms of δk for each system parameter are plotted in Figures 789.

It can be viewed that for each parameter case the mean value of δk from all histograms aligns with the true change value, and the 5% damage in k 4 is correctly quantified.

Figure 1 .

 1 Figure 1. Damage quantification flow chart

Figure 2 . 6

 26 Figure 2. 6 DOF chain system sketch In total 1000 data sets for each case of Q are realized.In each data set the damage is modeled as a stiffness reduction of the 4-th spring by 5%. Hereafter let θ be the collection of the element stiffness k 1-6 . Albeit the damages considered herein relate to the stiffness changes, the proposed approach is general and can be formulated to quantify changes in any parameter of the numerical model for which the modal parameter sensitivities are nonzero. For damage quantification the procedure outlined in Algorithm 1 is used.

Figure 3 and (k 4 , k 5 ) in Figure 4 .

 3454 Note that the remaining spring stiffnesses are (k 1 = 100, k 2 = 200, k 6 = 200).

Figure 3 .

 3 Figure 3. Objective function F (θ) for the parameter pair k 3 and k 4 .

Figure 4 .

 4 Figure 4. Objective function F (θ) for the parameter pair k 4 and k 5 .

Figure 5 .

 5 Figure 5. CMA-ES optimization of the change in parameters k 3 and k 4 .

Figure 6 .

 6 Figure 6. Zoom on the changes in k 3 and k 4 corresponding to the acceptance region in the last iteration.

Figure 7 .

 7 Figure 7. CMA-ES optimization of the change in parameters k 1 and k 2 .

Figure 8 .

 8 Figure 8. CMA-ES optimization of the change in parameters k 3 and k 4 .

Figure 9 .

 9 Figure 9. CMA-ES optimization of the change in parameters k 5 and k 6 . 6 CONCLUSIONS In this paper, a damage quantification strategy has been developed that uses the difference between a data-driven and a model-based Hankel matrices. A stopping criterion based on the quantile of the statistical law of a Mahalanobis metric based on these two matrices has been used in the CMA-ES-based optimization. Preliminary localization results have clustered the parameter space and shown what clusters were potentially damaged. It has then been shown that the damage extent is clearly identified for damaged elements in different clusters. To show the impact of clustering and how damage can not be identified and

Algorithm 1 :

 1 Model optimization Input : output measurements {y k } k=1,...,N +p+q ; physical parameter θ init ; parameterized model K θ , M θ , D θ ; optimization parameters from Table 1 ; Output: estimated parameter value θ sol and estimated total change δ for each cluster ; 1 compute Ĥθ * data from (8) and the estimate of Σ data ; 2 create K θinit , M θinit , C θinit and compute H θinit get λ model candidates θ 1 , . . . , θ λ in (24) ; 7 repeat steps 2, 3 and 4 with θ 1 , . . . , θ λ ; count the number t Θ of model candidates θ j with θ j ∈ Θ in the acceptance region, see (23) ; 10 until t Θ ≥ t opt ; 11 compute θ sol in (25) as the mean of all model candidates in Θ ;

		(6) and (7) ;						model from
	3 compute the residual ζθ (15), the estimate of Σ ζ (19)
		and J ζ θ * (20) ;
	4 evaluate F (θ init ) in (22) ;
	5 repeat
	6
	8	update CMA-ES parameters after [31] ;
	12 compute the change in the parameter value δc for
		each damaged cluster c in (26)
		k1	m1	k2	m2	k3	m3	k4	m4	k5	m5	k6	m6
		sensor 1			sensor 2			sensor 3

9

Table 1 .

 1 CMA-ES optimization parameters

	σ init	θ i init	θ i min-max	λ	topt
	60	150	75-225	1000 60% λ
	parametrized with spring stiffnesses k j . The corresponding
	function F (θ) is displayed for the parameter pair (k 3 , k 4 ) in

maps the p-dimensional parameter space to a multidimensional hyperplane whose shape can easily be illustrated for p = 2. Consider the estimate of H θ * data computed from the available data, and H θ model


Boualem Merainani
email: boualem.merainani@inria.fr

Bian Xiong
email: bian.xiong@inria.fr

Vincent Baltazart

Jean Dumoulin

Michael Döhler

Qinghua Zhang

Experimental investigation of structural modal identification using pixels intensity and motion signals from video-based imaging devices: performance, comparison and analysis

Keywords: Vision modal decomposition, Pixels intensity, Motion extraction, Subspace identification, Uncertainty bounds

This paper aims to experimentally evaluate a simplified vision-based method for structural health monitoring (SHM). Contrary to conventional solutions that rely on extracting motions through image processing, this paper proposes to conduct the SHM analysis by the direct processing of pixel brightness without extracting the motion signals beforehand. After some pre-processing steps, it is shown that the brightness data reveal essential information about the dynamic characteristics of the monitored vibrating structure. Furthermore, the low-level information of the pixel is compensated by an efficient selection of the so called "active pixels" throughout the image time flow. Finally, a subspace system identification-based method is applied to the brightness data, so that the modal parameters with uncertainty bounds are estimated within a large range of model orders displayed in a stabilization diagram. The experiment database consists of image time flows of a cantilever beam excited by a shake table driven by a band limited random noise. Modal vibrations range from 1 to 173 Hz.

INTRODUCTION

Many infrastructures of our society, including bridges and buildings, are facing increasing challenges due to climate change. To this end, extensive research works have been carried out over the past decades for structural health monitoring (SHM). [START_REF] Burgos | Damage identification in structural health monitoring: A brief review from its implementation to the use of data-driven applications[END_REF][START_REF] Zona | Vision-based vibration monitoring of structures and infrastructures: An overview of recent applications[END_REF][START_REF] Shokravi | Application of the subspace-based methods in health monitoring of civil structures: A systematic review and meta-analysis[END_REF] Operational modal analysis (OMA) is a mature field in SHM. Its goal is to identify inherent structural information such as natural frequencies, damping ratios, and mode shapes from output-only vibration measurements on a structure under unknown ambient excitation. OMA is of high interest for many applications, e.g. damage diagnosis, [START_REF] Döhler | Structural health monitoring with statistical methods during progressive damage test of s101 bridge[END_REF][START_REF] Guo | Damage identification for frame structures using vision-based measurement[END_REF][START_REF] Allahdadian | Towards robust statistical damage localization via model-based sensitivity clustering[END_REF] structural model updating, [START_REF] Doebling | A summary review of vibration-based damage identification methods[END_REF] identification of excitation forces. [START_REF] Feng | Identification of structural stiffness and excitation forces in time domain using noncontact vision-based displacement measurement[END_REF] The process of OMA starts with data collection. Conventional sensing technology was based on setting accelerometers on the structure. Then, methods like frequency domain decomposition, 9, 10 auto-regressive moving average (ARMA), [START_REF] Carden | Arma modelled time-series classification for structural health monitoring of civil infrastructure[END_REF] stochastic subspace identification [START_REF] Peeters | Reference-based stochastic subspace identification for output-only modal analysis[END_REF][START_REF] Döhler | Fast multi-order computation of system matrices in subspace-based system identification[END_REF] and the natural excitation technique [START_REF] James | The natural excitation technique (next) for modal parameter extraction from operating wind turbines[END_REF][START_REF] Siringoringo | System identification of suspension bridge from ambient vibration response[END_REF] can be used to obtain the modal parameters with high accuracy. However, with traditional sensors, vibrations are measured only at (relatively few) discrete points in the structure, providing limited information for damage diagnosis or model updating. [START_REF] Fan | Vibration-based damage identification methods: a review and comparative study[END_REF][START_REF] Ewins | Modal testing: theory, practice and application[END_REF] Non-contact sensing technology, such as laser sensor [START_REF] Nassif | Comparison of laser doppler vibrometer with contact sensors for monitoring bridge deflection and vibration[END_REF] and ground-based synthetic aperture radar, [START_REF] Negulescu | Comparison of seismometer and radar measurements for the modal identification of civil engineering structures[END_REF] have been successfully applied as an alternative to traditional devices. Whereas numerous applications benefit from these measurement technologies, they are expensive and not adapted to all types of structures.

With both technological and processing improvements, computer vision-based techniques have emerged in the past decade as an alternative, cheaper and easier solution to perform SHM. Commonly, this is done in two steps. First, the vibrating motion is estimated at long-range from video image flows. Various image processing techniques have been developed to extract the local motion signals from subtle changes in the pixel brightness over time. One can classify the approaches to video-based motion extraction into optical flow , [START_REF] Guo | Vision-based measurement for rotational speed by improving lucas-kanade template tracking algorithm[END_REF][START_REF] Javh | The subpixel resolution of optical-flow-based modal analysis[END_REF] including phase based optical flow, [START_REF] Wadhwa | Phase-based video motion processing[END_REF] and the image template matching, including the widespread digital image correlation (DIC). [START_REF] Wang | Frequency response functions of shape features from full-field vibration measurements using digital image correlation[END_REF] Second, the most dominant vibrating frequencies of the structure are determined from the power spectral density of the motion time history. More advanced SHM tools rely on using the stochastic subspacebased identification (SSI) method. [START_REF] Döhler | Efficient multi-order uncertainty computation for stochastic subspace identification[END_REF] By contrast, this paper proposes to conduct SHM by the direct processing of pixel intensity time variations without extracting the motion signals beforehand. The proposed approach is believed to be less computationally expensive. To circumvent the low-level information of the pixel intensity, the proposed approach is able to represent the brightness time variations within each image blocks (ROI : region of interest) by a single time signal obtained after selecting the so called 'active pixels' and averaging their time variations. Then, the stochastic subspace-based identification method [START_REF] Döhler | Efficient multi-order uncertainty computation for stochastic subspace identification[END_REF] that provides more accurate modal parameters along with uncertainty bounds was adopted for refined SHM analysis.

Experiments in a laboratory setting on a vertical cantilever beam are performed to verify the proposed intensity-based processing approach. A comparison has been made with conventional motion-based processing in terms of identified spectral lines and on the stability diagrams of modal identification.

EXPERIMENTAL SETUP

A vertical cantilever beam of 2 m long, 30 mm wide, 5 mm thick, and of density 7800kg/m 3 was tested to evaluate the performance of the proposed method for vision modal identification. Details of the experimental setup are given bellow and shown in Fig. 2: • The vertical beam was fixed at its bottom to a shake table (model: NEBS GR-63), which was driven in random mode horizontally in the direction perpendicular to the axis of the camera.

• Video images were recorded at the frame rate 512 f ps and at the frame resolution 1920 × 136 with a PHANTOM camera (model :MIRO 320S), mounted with a NIKKOR lens with a fixed 50 mm focal length.

The camera was set on a tripod and installed at a distance of about 5.6 m from the beam, so that the full beam length was in the scene. The scaling factor was 1.02 mm/pixel.

• Three additional lighting (Halogen Spotlight, MAZDA CORMORAN type of 1500 W) were used to improve the brightness conditions; in addition, a white curtain was used as a backdrop to enhance the contrast with the structure edges.

• For motion tracking, three paper markers were mounted at some particular positions along the beam, namely, at the top (1/1), middle (1/2) and 0.5 m (1/4) from the cantilever beam bottom. Please see Fig. 4 (top left) for illustration of the marker.

• Three laser displacement sensors were located at the same heights of the markers to serve as the motion benchmark.

IMAGE PROCESSING FOR SHM

In this section, the images of the vibrating cantilever are decomposed into image blocks (or ROIs). Two processing are performed on each ROI assuming a solid state motion without deformation. The first processing aims at computing a specific averaged of the brightness time variation. The second one estimates motion signal by the conventional video-based technique. For both cases, The time signals over the set of ROIs will be processed by SSI method for SHM purposes in Section 4.

Proposed Image Brightness Processing

The proposed image processing includes the selection of active pixels, some phase correction of the brightness time variations, before computing the final averaged time signal.

Selection of active pixels

Active pixels are supposed to bring the most information about the vibrating motion. Active pixels usually correspond to the edge between adjacent contrasted areas (The amount of active pixels increase when the image is highly textured, e.g., speckle pattern or random dot pattern). The brightness time variation of such pixels is then closely correlated to the structure motion at that location. According to Fig. 2, this is easily seen when the motion range is within a pixel. When the motion covers a few pixels range, the brightness time variations of each pixel is truncated in amplitude. The correlation with the structure motion is then recovered on the averaged brightness over some pixel neighborhood.

Many algorithms were developed in the literature to locate the edges from the detection of significant local spatial changes in images by either the gradient operators, [START_REF] Hartley | A. zisserman multiple view geometry in computer vision[END_REF] Gaussian derivatives, [START_REF] Ji | Nontarget stereo vision technique for spatiotemporal response measurement of line-like structures[END_REF] or the steerable filters [START_REF] Freeman | The design and use of steerable filters[END_REF] as illustrated in Fig. 4 (top row).

Within the scope of application, the edge computation is usually performed once from the first image frame. These methods are then better tailored to small motion range, namely within one pixel. However, in practical condition, the motion amplitude range is not known beforehand and can be beyond one pixel.

By contrast, a spectral approach is proposed in this section to select the active pixels without any prior knowledge. Indeed, it can be used to determine the pixel-wise motion range. Active pixels are selected within each ROI on the basis of the spectral magnitude of the pixel brightness time history. The proposed selection of pixels is based on the following three-steps processing. • First, the Beam Resonance Frequencies (BRF) of the mechanical structure are determined. If the finite element modelling of the structure is available, BRF can be theoretically determined. Otherwise, the main BRFs are determined from the peak magnitudes of the pixel spectrogram. The histogram of the frequency peak location is illustrated in Fig. 3. At this step, the spectrogram is approximated because the pixel brightness may bring truncated information (see Fig. 5). Nevertheless, it is sufficient to determine the set of BRFs.

× 3 Pixels intensity

Line progression

• Second, the set of BRFs is used to compute the energy of ROI pixels according to Algorithm 1. The maximum energy value serves as the reference level to select active pixels.

• Finally, the pixels having energy above 20% of the latter reference energy level, are selected as active pixels.

Algorithm 1: Active pixel selection within N × P ROI 1 Inputs: -(x, y) pixel intensity time history, I I I(x, y, t t t), where, t t t is a L × 1 time vector; The results of active pixels selection procedure are illustrated in Fig. 4 (bottom row) and by the red pixels in Fig. 5. Compared with the conventional edge pixels selection at the top row, the proposed spectral approach achieves more accurate detection. On the one hand, the target edges which are perpendicular to the motion were accurately identified. On the other hand, the width of the detected edge is closely related to the pixel-wise motion amplitude, measured by the laser sensor.

Averaged brightness ROI signal

For convenience, the brightness time variations of the selected active pixels (as shown by red pixels in Fig. 5) are arranged in a matrix form I Pa where each column I I I(x i , y i , t t t) represents the centered brightness time vector of the i th selected active pixel at the coordinates (x i , y i ) in the ROI image.

At this stage, the pixel brightness time variations may bring truncated information in case of motion magnitude larger than one pixel, as shown on Fig. 5 (right column). In addition, the brightness time variations have the same polarity than the brightness slope at the edges. [START_REF] Javh | The subpixel resolution of optical-flow-based modal analysis[END_REF] When the ROI includes different dark-to-light and light-to-dark transitions, the pixel brightness time variations show various polarity, which tend to cancel when they are summed up. As shown in, [START_REF] Javh | The subpixel resolution of optical-flow-based modal analysis[END_REF] phase correction should be performed to achieve a reliable averaged signal with the full vibrating content.

To this aim, the following polarity correction processing is proposed :

∀i ∈ [2, P a ], I I I(x i , y i , t t t) = I I I(x i , y i , t t t) × sign(I I I(x i-1 , y i-1 , t t t) T I I I(x i , y i , t t t)).

(

) 1
where P a is the amount of selected active pixels. The phase correction ensures that brightness of all selected active pixels are in phase before averaging. For illustration, Fig. 6 represents the processing result on the active pixels in Fig. 5.

Conventional video-based motion extraction

Generally, the motion extraction from video images is performed in two steps: pixel shift (∆n, ∆p) estimation using the cross-correlation (CC) and subpixel shift refinement (δx, δy). Let I 0 = I(n, p, t0) and I 1 = I(n, p, t) be the reference and a template images, respectively. The pixel shifts can be obtained by maximizing the CC as :

(∆n * , ∆p * ) = arg max -N ≤∆n≤N -P ≤∆p≤P u,v F (u, v)G * (u, v)e i2π( u∆n N + v∆p P ) , (2)
where F (u, v) , and G(u, v), are the Fourier transforms of I 0 and I 1 respectively.

Chan et al. [START_REF] Chan | Subpixel motion estimation without interpolation[END_REF] proposed a simplified version of classical optical flow estimation [START_REF] Irani | Improving resolution by image registration[END_REF] to determine the 2D subpixel displacements by solving the system of equations given in 3. It is noteworthy that, in case of large displacements, the reference image must be shifted as I(n + ∆n, p + ∆p, t0) to stay in the range of subpixel displacements.

δx δy = x,y ( ∂I0 ∂x ) 2
x,y ∂I0 ∂x ∂I0 ∂y

x,y ∂I0 ∂x ∂I0 ∂y

x,y ( ∂I0 ∂y ) 2 -1

x,y (I

1 -I 0 ) ∂I0 ∂x x,y (I 1 -I 0 ) ∂I0 ∂y , (3)
Comparison between the two processings Figure 6, shows and compares the time variations of the average ROI brightness (within the second ROI) and the displacements measured by both the vision and the reference laser sensors. The results shows a good consistency between the three signals.

On Fig. 7, the correlation coefficients (CorCoef) between the displacements and the brightness time variations are computed. It is seen that, the CorCoef varies from 0.996 to 0.9999 and shows the strong relationship between the intensities and the displacements as expected. It shall be noted, however, that the correlation is slightly lower for the upper ROIs, in which the structure motion is larger than one pixel. This is likely due to either insufficient active pixels in the ROI, non-uniform lighting conditions in the upper ROIs or inappropriate threshold. Figure 8, shows the frequency spectrogram associated to the time signals in Fig. 6. For comparison, the spectrum computed from the average of 3 adjacent active pixels, represented in Fig. 5, was also superimposed on that of the average ROI brightness over all active pixels. The peak prominent indicator 30 was adopted to quantify the local signal to noise ratio of each resonance frequency. Figure 9 compares the results. On the one hand, the peaks, related to the first 4 resonant frequencies, in the laser sensor-derived spectrum, have a larger prominence compared to those derived from the spectra of the vision sensor and the proposed LI-derived vision data approach. On the other hand, our approach shows better peaks prominence, related to the fifth and sixth It is worth to mention that, higher modes of the beam were detected in the spectra. However, The laser sensor (Model : YP11MGVl80) cut-off frequency was 100Hz. Beyond this value, the frequency content of the measured displacements may be affected. This is why, the comparison was limited over the [0-100]Hz frequency range.

SYSTEM IDENTIFICATION RESULTS

In this section, the computation of the modal parameters of the beam and their uncertainty bounds is performed using the covariance-driven subspace identification algorithm. [START_REF] Döhler | Efficient multi-order uncertainty computation for stochastic subspace identification[END_REF] With this method, the system matrices of a linear time-invariant state-space model are identified from each type of the measurements (LIs, displacement extracted from video images and laser sensor signal). From their eigenvalues and eigenvectors, the modal parameters (natural frequencies, damping ratios, mode shapes) are retrieved. Since the model order of the system (corresponding to twice the number of modes that are present in the data) is unknown in practice, the system is identified for many successive model orders. The resulting frequencies are plotted against the model order in a stabilization diagram, which allows to distinguish "true" structural modes (which are stable for successive orders) from spurious noise modes.

In Fig. 10, the stabilization diagrams are shown together with the spectra, where the ±σ confidence interval of each frequency is plotted as a horizontal bar. It can be seen that, for the three analysis, real physical poles clearly manifest as stable poles from low model order (< 40). However, a high order implies the appearance of spurious modes that may be associated to the noise content of the signals. In addition, the stabilization diagram of the LI data, in Fig. 10.a, shows a stable mode at 100Hz but it has no physical meaning. In fact, this frequency is caused by the source light flicker.

From the stabilization diagrams, the modal parameters are chosen. In Table 1, the system identification results with their standard deviations are summarized. One can observe that between the identification results from the laser sensor displacements, vision based measurements and the LIs, the maximum error is around 0.6% of in terms of natural frequency estimation.

CONCLUSION

This paper presents a framework for structural modal identification using vision measurements. The key idea lies in the fact that the intensities at some pixels (titled active pixels) carry the useful information of the structural vibration. Therefore, an algorithm has been proposed to select such pixels. Modal parameters are then identified from the collected vision data using the output-only subspace identification method. Experiments concerning the video of a vertical cantilever beam excited with a shake table is studied and results show that • Only pixels that exhibit intensity changes relating to the beam vibration are selected. Therefore, the number of active pixels change with respect to motions amplitude (see Fig. 4).

• Strong relationship exists between the time variation of the intensities and displacements extracted from Vision sensor with a minimum CorCoef of 99, 6%.

• The beam modal parameters including the natural frequencies and the damping ratios with their uncertainties were identified. A comparison between the SSI results when implemented on the intensity data, laser and vision data were also studied. The results showed around 2% of maximum error in terms of natural frequency estimation. The comparison indicates that the proposed approach has the potential to identify the structure modal parameters with reasonable levels of accuracy.

The proposed approach can estimate the modal properties. However, it is critical to efficiently select the active pixels. If they are insufficient or over-considered, the intensity time variations will be incomplete or corrupted by noise. Furthermore, the proposed approach depends highly on the illumination conditions. A constant illumination during tests is required to obtain reasonable results.

Authors are grateful to Research Fund for Coal and Steel for funding part of this work under grant agreement No 800687 in the framework of DESDEMONA project, and to Ivan Guéguen, Jean-Philippe Gourdon and Fabrice Blaineau for the essential assistance provided for the experimental set up.

Figure 1 :

 1 Figure 1: Experimental setup: (left) video cameras in the foreground, the cantilever beam in the background, (right) a close-up view of the lighting conditions, the cantilever beam and the shaking table.

Figure 2 :

 2 Figure 2: Schematic representation of the pixel brightness variation w.r.t. the motion of some step edge in front of the camera. Top row figures: the right edge of a black object moving from the left to the right along 3 pixels. Bottom row figures: the corresponding pixel brightness variations of the camera.

Figure 3 :

 3 Figure 3: Frequency peaks histogram computed from the brightness time variations of the pixels within the selected ROI.

2 -K

 2 Beam resonance frequencies BRF ; 3 Output: Pixels energy; 4 for x = 1, . . . , N do 5 for y = 1, . . . , P do 6 for k = 1, . . . , K do 7 Band-pass filtering around BRF (k) : I I I k,f iltered ; 8 Mean intensity computation I I I k,f iltered ; 9 end 10 Energy computation E x,y = sum(|I k,f iltered | 2 ); 11 end

Figure 4 :

 4 Figure 4: Pixels energy calculation based proposed approach prior active pixels selection in case of 4 different motions amplitude (Bottom), target (Top left), gradient (G x 2 ) (Top middle) and steerable filter (Top right).

Figure 5 :

 5 Figure 5: Illustration of the pixel selection: (left) Selected active pixels in red circles are superimposed on the target pattern; (right) Brightness time variations of pixels 1-3 within the boxcar.

Figure 6 :

 6 Figure 6: Averaged brightness time variations (right y-axis) superimposed on the displacements measured by both the laser and the vision sensors (left y-axis); full time record on the left, zoom on two short time intervals on the right.

Figure 7 :

 7 Figure 7: Correlation coefficients between the brightness time variations and vision-based displacements

Figure 8 :

 8 Figure 8: FFT spectra of the average LI time history of all active pixels and 3 adjacent active pixels, showing them on top of each other (right y-axis), FFT spectra of the displacements measured by the laser sensor superimposed on that obtained by the vision sensor (left y-axis).

Figure 9 :

 9 Figure 9: Comparison between the peaks prominence, of the first 6 resonant frequency lines obtained from the spectra of the laser sensor, the vision sensor and the proposed LI-derived vision data

Figure 10 :

 10 Figure 10: Stabilization diagrams containing the natural frequencies of the first eight modes with their confidence intervals.

Table 1 :

 1 Overview of the first eight modes

	Mode		Intensity		Displacement measured by Laser sensor	Vision-based displacement measurements
		f (Hz)	σ f (Hz)	ζ(%)	σ ζ (%)	f (Hz)	σ f (Hz)	ζ(%)	σ ζ (%)	f (Hz)	σ f (Hz)	ζ(%)	σ ζ (%)
	1	0.9944	0.0149	0.4018	1.6236	1.0006	0.0071	0.0145	0.8407	0.9982	0.0092	0.0476	0.9035
	2	6.7410	0.0110	0.2561	0.2118	6.7421	0.0099	0.1881	0.1346	6.7419	0.0124	0.2412	0.1644
	3	19.0076	0.0359	0.2196	0.1487	19.0173	0.0127	0.1044	0.0537	18.9960	0.0216	0.2620	0.1547
	4	37.3424	0.0195	0.1331	0.0482	37.3443	0.0182	0.1537	0.0465	37.3406	0.0185	0.1319	0.0443
	5	62.0247	0.0220	0.1379	0.0363	62.0263	0.0272	0.1263	0.0497	62.0245	0.0221	0.1335	0.0391
	6	92.5723	0.0224	0.0823	0.0206	92.5681	0.0168	0.0522	0.0159	92.5742	0.0288	0.0858	0.0456
	7	129.6552	0.0312	0.0791	0.0255	129.6541	0.0288	0.0737	0.0244	129.6795	0.0220	0.0749	0.0217
	8	172.1089	0.0263	0.0465	0.0146	172.1081	0.0255	0.0570	0.0148	172.1133	0.0322	0.0502	0.0164