Mario D Baquedano-Aguilar
email: m.baquedanoaguil@ufl.edu

Arturo Bretas
email: arturo.bretas@pnnl.gov

Sean Meyn
email: meyn@ece.ufl.edu

Nader Aljohani
email: nzjohani@taibahu.edu.sa

Tri-Level Linear Programming Model for Automatic Load Shedding Using Spectral Clustering

Keywords: Automatic load shedding, intentional control islanding, linear programming, phasor measurement units, spectral clustering, state estimation

Traditional load shedding schemes can be inadequate in grids with high renewable penetration, leading to unstable events and unnecessary grid islanding. Although for both manual and automatic operating modes load shedding areas have been predefined by grid operators, they have remained fixed, and may be sub-optimal due to dynamic operating conditions. In this work, a distributed tri-level linear programming model for automatic load shedding to avoid system islanding is presented. Preventing islanding is preferred because it reduces the need for additional load shedding besides the disconnection of transmission lines between islands. This is crucial as maintaining the local generation-demand balance is necessary to preserve frequency stability. Furthermore, uneven distribution of generation resources among islands can lead to increased load shedding, causing economic and reliability challenges. This issue is further compounded in modern power systems heavily dependent on nondispatchable resources like wind and solar. The upper-level model uses complex power flow measurements to determine the system areas to shed load depending on actual operating conditions using a spectral clustering approach. The mid-level model estimates the area system state, while the lower-level model determines the locations and load values to be shed. The solution is practical and promising for real-world applications.

I. INTRODUCTION

W ITH the increasing penetration levels of distributed energy resources (DERs) and the decommissioning of conventional generating units due to efforts in decarbonization policies and grid modernization, power system operations are becoming highly complex [START_REF] Rockström | A roadmap for rapid decarbonization[END_REF]. DERs are typically connected through inverters, which poses technological challenges to grid operation [START_REF] Ulbig | Impact of low rotational inertia on power system stability and operation[END_REF]. Fast and reliable protection systems are crucial for stabilizing power system dynamics and preventing blackouts. Phasor measurement units offer opportunities for advanced wide-area measurement protection and control (WAMPAC) systems, but traditional hierarchical communication architectures might not meet real-time requirements.

Load shedding is the last resort to maintain system reliability and avert blackouts caused either by power capacity shortages or stability issues [START_REF] Baquedano-Aguilar | Impact of increased penetration of large-scale pv generation on shortterm stability of power systems[END_REF]. Current load-shedding schemes comprehend manual and automatic approaches [START_REF] Omar | Under frequency load shedding (UFLS): Principles and implementation[END_REF], [START_REF] Terzija | Adaptive underfrequency load shedding based on the magnitude of the disturbance estimation[END_REF]. In manual operating mode, operators notify utilities of anticipated stressful conditions, leading to planned load shedding.

Automatic load shedding otherwise employs under-frequency relays at substations. Loads are shed in predefined stages based on peak net load estimates. Although the existing framework has ensured reliability, it may yield sub-optimal outcomes in power grids with high renewable penetration. Advanced load-shedding solutions with centralized architectures face challenges in real-life applications due to communication technology latency.

To address these issues, [START_REF] Bretas | Bi-level linear programming model for automatic load shedding: A distributed wide-area measurement system-based solution[END_REF] presented a bi-level linear programming model for automatic load shedding utilizing a distributed architecture while leveraging data from phasor measurement units (PMU). Although this model shows promising results for real-world applications, two strong assumptions were made: when determining system areas for load shedding, they were previously user-defined and after kept static. Regularly, determining areas to shed load requires the grid operator's knowledge of the system and system disturbance identification. However, recent advancements in related fields like intentional control islanding [START_REF] Quirós-Tortós | Constrained spectral clustering-based methodology for intentional controlled islanding of large-scale power systems[END_REF], have used the system topology, the admittance matrix, or transmission line ratings to pinpoint these areas [START_REF] Sánchez-García | Hierarchical spectral clustering of power grids[END_REF], [START_REF] Baquedano-Aguilar | Reduced-order models of static power grids based on spectral clustering[END_REF]. Nevertheless, these parameters remain constant in any power system and ignore the dynamic operating conditions introduced by DERs. To overcome this challenge, in this work PMU power flow measurements are considered in an upper-level model to determine system areas best to shed loads. Spectral clustering is used towards model solution. This information is then sent to a mid-level model, which employs a distributed weighted least squares constrained approach to estimate the system state. Decision variables at the mid-level model are then integrated into the lower-level model, which uses system estimated state to determine the optimal location and the load to be shed in each area. The key contributions of this work are:

• A formal linear programming tri-level model for underfrequency load shedding using a distributed architecture; • The determination of optimal areas to shed load according to the system's actual operating and disturbance conditions.

The remaining of the paper is organized as follows. Section II presents the tri-level model. A study case using the IEEE 14-bus test system is presented in Section III, while conclusions and future research topics are contained in Section IV.

II. TRI-LEVEL MODEL

The proposed tri-level linear programming model is composed of upper-level, mid-level, and a lower-level model. A schematic illustration is shown in Fig. 1. In contrast to traditional algorithm execution on a centralized computer that has communication capabilities with the load-shedding relays, we envision our application on an energy management system (EMS) that will be running in the server of the utility considering critical infrastructure protection (CIP) zones. The North American Electric Reliability Corporation (NERC) defines the CIP zones to ensure the protection of assets used to operate North America's bulk electric systems.

A. Upper-Level Model

The upper-level model estimates best system areas to shed load according to operating and disturbance conditions. These areas are determined and continuously updated in less than t < 90ms complying with the range of PMU processing times [START_REF]IEEE Standard for Synchrophasor Measurements for Power Systems[END_REF]. To avoid communication issues during fault events, the model uses the most recent set of PMU power flow measurements prior to a disturbance to construct the network's graph representation, which is essential for area determination.

1) Determining k-Shedding Areas: To determine areas where load can be shed according to actual operating conditions, we first need to build the power flow-based graph Laplacian matrix W . For this, we represent the power grid as an undirected graph, denoted as G = (V, E). The set of buses is V = {1, 2, ..., n}, where n is the total number of buses. The set of edges is E ⊂ V × V, where (i, j) ∈ E represents an edge connecting buses i and j. Edges can be transmission lines or transformers. For a network with n buses, the weighted adjacency matrix of the graph is

W =    0 ω 12 • • • ω 1n . . . . . . . . . . . . ω n1 ω n2 • • • 0    , (1)
where ω ij represents an edge weight as a non-negative function ω : V × V→ R + such that

• ω i,j = 0 if (i, j) / ∈ E, i.e.
, vertices i and j are not connected by and edge.

• ω i,j = ω j,i if i ̸ = j, i.e, edge directions are ignored since G is undirected.

Power Flow Weighted Matrix. In this work, complex power flow (in MVA) measurements are used to define edge weights.

To avoid any specific weight from leading to inaccurate clustering results due to fluctuating operating conditions reflected in the complex power flow magnitudes, we transform the resulting matrix W into a normalized Laplacian matrix.

Graph Laplacians. In this work, we use the normalized symmetric graph Laplacian matrix L n , which is commonly employed for clustering purposes (see [START_REF] Sánchez-García | Hierarchical spectral clustering of power grids[END_REF], [START_REF] Baquedano-Aguilar | Reduced-order models of static power grids based on spectral clustering[END_REF] and references therein). The entry-wise definition of both the unnormalized (2) and normalized (3) graph Laplacian matrices follows:

L ij =      d i if i = j -1 if i is adjacent to j 0 otherwise, (2)
where

d i = ω(i, i) = n j=1 ω(i, j
) is the degree of the ith vertex. This is closely related to the adjacency matrix and is sometimes written as L = D -W , where D is defined as the diagonal matrix with weighted vertex degrees d 1 , • • • , d n on the diagonal (i.e., d i ). The normalized graph Laplacian matrix L n is defined as:

L n = D -1/2 LD -1/2 = I -D -1/2 × W × D -1/2 . (3)
The eigenvalues of the normalized Laplacian matrix satisfy the inequality 0 ≤ λ Ln i ≤ 2 for all i [START_REF] Chung | Spectral graph theory[END_REF]. The interested reader is referred to [START_REF] Baquedano-Aguilar | Reduced-order models of static power grids based on spectral clustering[END_REF] for a detailed explanation of the preference of L n over L in the power systems area.

2) Choice of Dimension: In situations where information such as the number of coherent groups of generators in the system is limited, a typical approach to determining the number of areas to shed load involves calculating the ratio of the absolute value of the difference between two consecutive eigenvalues of L n and their size. This parameter serves as a similarity measure known as the relative eigengaps [START_REF] Chung | Spectral graph theory[END_REF]:

γ k,rel = |λ k+1 -λ k | λ k (k ≥ 2). (4)
The presence of favorable k-areas is indicated by a small value of the kth eigenvalue, denoted as λ k in the graph Laplacian. Then, a high value of γ k,rel means that a minimum of k-shedding areas can be determined. This characteristic is expected to be evident through the spectral embedding in the k-dimensional space.

3) The Algorithm: In this work, we use the normalized spectral clustering algorithm according to [START_REF] Shi | Normalized cuts and image segmentation[END_REF] with some minor modifications.

Algorithm 1 The Spectral Clustering Algorithm.

INPUT: Matrix W ∈ R n×n ; parameter k. Compute the normalized Laplacian L n . Compute the eigenvectors v 1 , • • • , v k of the generalized eigenproblem L n v = λDv
and take the first k eigenvectors Let V ∈ R n×k be the matrix containing the vectors

v 1 , • • • , v k as columns. For i = 1, • • • , n let y i ∈ R k be the vector corresponding to the ith row of V . Cluster the points (y i ) i=1,••• ,n in R k with k-means or agglomerative clustering into areas C 1 , • • • , C k . OUTPUT: Areas A 1 , • • • , A k with A i = {j|y j ∈ C i } B. Mid-Level Model
This model estimates the system state in each of the k areas determined by the upper-level model. For a power system comprised of n buses and m measurements, it is possible to represent it as a collection of linear algebraic equations, as follows (see [START_REF] Bretas | Bi-level linear programming model for automatic load shedding: A distributed wide-area measurement system-based solution[END_REF] and references therein):

z = H(x) + e, (5)
where z ∈ R m is the measurement vector, x ∈ R N is the state variables vector, H(x) : R N -→ R m , (m ≥ N ) is a linear function that relates the states to the measurements, and e is the measurement residual vector with a Gaussian probability distribution, zero mean, standard deviation σ, and N state variables such that N = 2n -1. We can now write (5) for k areas, n buses, and m measurements as follows:

z k = H k (x k ) + e. (6)
The distributed weighted least squares (WLS) state estimator aims to find optimal estimates of the system states x k by minimizing the objective function J k (x k ):

J k (x k ) = e(z k , H k ) T * R -1 k * e(z k , H k ), (7)
where in area k, the residual vector e(z k , H k ) is represented as

(z k -H k (x k ))
, and R k stands for the measurement covariance matrix of the kth area. J k (x k ) is essentially an L 2 -norm in the measurements vector space R m k of area k. Notice that R k for each area is a symmetric positive definite matrix with diagonal elements determined by the magnitudes of the corresponding measurements. The mid-level model's decision variables consist of the voltage angles of the buses within each area. The overall objective cost function for the distributed state estimation architecture is obtained by summing up the individual objectives of areas k = 1, 2, 3. min

x k ∈χ k ,x kl K k=1 e(z k , H k ) T * R -1 k * e(z k , H k ) (8)

C. Lower-Level Model

The lower-level model utilizes the mid-level estimated system states to determine both the location and the load to be shed. The primary objective is to reduce the effects of load shedding in the system. This is quantified by calculating the mismatch between the current system state obtained from the mid-level model, and the estimated state after implementing load shedding. To achieve this, the L 1 -norm of the difference between the current and the estimated system states is used. The multi-objective optimization formulation of [START_REF] Bretas | Bi-level linear programming model for automatic load shedding: A distributed wide-area measurement system-based solution[END_REF] has a secondary goal: to minimize the load flows through the tielines of the k areas. However, as a consequence of the spectral clustering approach, this additional goal is unnecessary since the algorithm already ensures that the system is divided into k areas with tie-lines carrying the minimum complex power flow, resulting in minimal energy disruption. The interested reader is referred to [START_REF] Sánchez-García | Hierarchical spectral clustering of power grids[END_REF] and references therein for a detailed explanation of determining system areas considering minimum power flow disruption. The formulation of the lower-level model corresponds to a constrained linear optimization problem, expressed as follows:

min d∈Λ ω l δ l + ω p δ p s.t. N i=1 N j=1 |∆V ij -g l | = δ l ∆V ij = | Vi -Vj | + | Vout,i -Vi | + | Vout,j -Vj | S rs -S lim rs ≤ δ p ; δ p ≥ 0 k l i=1 d i ≥ 0.9 k l i=1 d mid,i ; S = XV, (9)
where δ l ∈ R is the distance to the primary goal ( Vout,i is the mid-level model decision variable of bus i, and Vi is the estimated system state after load shedding, and g l is a previously defined goal), ω l is its normalizing weight; δ p is a user-defined parameter, which represents the complex power flow violation over the limit for each tie-line connecting the k areas modeled as a soft constraint that penalizes the objective function using a high value of ω p so that every area adjusts its tie-line complex power flows to remain within the limits;

979-8-3503-1360-4/24/$31.00 ©2024 IEEE Λ is the set of area loads; S rs is the complex power flow between buses r and s; S lim rs is the complex power flow limit between buses r and s; d ∈ R k l is the load shedding vector, and k l is the area number of loads that can be shed; d mid ∈ R k l is the load vector estimated from the mid-level model; S ∈ R N is the intra-area complex power injection bus vector; X : R N -1 → R m , (m > N ), is a linear function that relates the states to the measurements; and V ∈ R N -1 is the state variables vector.

III. CASE STUDY

In this case study, the tri-level model is applied to the IEEE 14-bus system for three operating conditions: maximum, mid, and low demand. In all cases, the loss of the generator connected at bus 8 is considered. Each area is considered to have a measurement data set comprising 12 measurements obtained from GridSandbox [START_REF] Mustafa | Cyber-power co-simulation for endto-end synchrophasor network analysis and applications[END_REF]. This results in a local global redundancy level of four.

An under-frequency load shedding (UFLS) scheme is set to operate once the frequency drops below 59.5 Hz, 59.1 Hz, and 58.7 Hz by shedding 10% of the maximum demand with a time delay of 14 cycles between frequency thresholds. The preset value for g l in ( 9) is set to zero. All relay settings are adapted from PJM's automatic under-frequency relay setting program [START_REF][END_REF]. For the study case, islanding schemes are designed to activate after all three steps of the UFLS scheme have been executed, and the frequency drops to 58.5 Hz or below. This setting is based on the Southwest Power Pool's (SPP) automatic islanding scheme [START_REF]SPP PC UFLS Plan[END_REF] for which no timedelay is considered as in [START_REF] Baquedano-Aguilar | Implementational aspects for the characterization of low-inertia in power systems[END_REF]. Obtaining the Number of Areas. Complex power flow magnitudes in the system are used to update the weighted adjacency matrix W described in Section II-A1. Then, the relative eigengaps of L n are computed as described in Section II-A2. Fig. 2 shows the relative eigengaps of L (left) and L n (right) for the power flow-based Laplacian of the IEEE 14-bus test system for different operating conditions: maximum, mid, and low demand.

The maximum values of γ k,rel for the unnormalized and normalized Laplacians among all operating conditions are approximately 1.4 and 3.1 respectively. The violation of the inequality constraint in [START_REF] Omar | Under frequency load shedding (UFLS): Principles and implementation[END_REF] indicates that the unnormalized Laplacian yields a poor similarity measure, aligning with findings in [START_REF] Sánchez-García | Hierarchical spectral clustering of power grids[END_REF], [START_REF] Baquedano-Aguilar | Reduced-order models of static power grids based on spectral clustering[END_REF], [START_REF] Luxburg | A tutorial on spectral clustering[END_REF]. Therefore, we will also use the relative eigengaps of the normalized Laplacian matrix as a reliable similarity measure. Then, γ Ln k,rel = 3.1 suggests that -at least-k = 3 shedding areas can be determined. Determining the Buses by Area. By knowing the number of areas we seek to find, the spectral clustering algorithm described in Section II-A3 can be implemented to determine and continuously update the system buses that belong to each area. The resulting areas to shed load for all operating conditions are depicted in Fig. 3 and specified in Table I. Fig. 3. k = 3 Shedding areas of the IEEE 14-bus system for three operating conditions via spectral clustering.

Due to the lack of space, the tri-level model for automatic load shedding is tested under a generator outage for the maximum demand operating condition only (see case A in Fig. 3 andTable I).

1) Scenario A -Generator outage with traditional UFLS scheme: After generator at bus 8 is lost (at t = 1 seconds), the frequency across the system rapidly drops, and complex power flows over branches 4-9, 5-6, and 13-14 increase, as shown in Fig. 4. Two of the flows exceed the threshold set for the relay, which is indicated by dashed lines of the same color. After a short delay, these two tie-lines trip, causing another significant increase in the power flow magnitude of branch 5-6, causing it to trip as well. This leaves line 10-11 to be the only tie-line connecting the system, which is rapidly overloaded and trip. This cascading failure results in a frequency instability event 979-8-3503-1360-4/24/$31.00 ©2024 IEEE as shown in Fig. 5. This situation splits the system into three unstable islands, with frequency dropping as illustrated.

2) Scenario B -Generator outage with the Tri-level automatic load shedding model: With the same scenario as in A, the presented tri-level model would handle the event faster and more intelligently compared to a traditional UFLS scheme.

The tri-level model was able to proactively shed load at the most impacted buses to arrest system frequency prior to occurrence of significant deviations. As shown in Fig. 5, loads at buses 6 and 11 in Area 1, as well as the load at buses 10 and 14 in Area 2 (see Table I) are shed over several short steps. As we can see in Fig. 5 (see the frequency plot for scenario B), this scheme successfully prevents the cascading failure seen in the previous case, and brings the system frequency back to stable operation. Overshooting Drawback. Given the fact that the tri-level model considers post-contingency predictive state estimates instead of frequency measurements, it continues to shed additional load so that islanding is still likely to occur due to line relaying, even though frequency appears to be recovering.

IV. CONCLUSIONS AND FUTURE WORK

In this work, a tri-level linear programming model for automatic load shedding to avoid system islanding using spectral clustering is presented. The upper-level model continuously updates the k shedding areas considering minimum power flow disruption using a spectral clustering approach. The system state was estimated in the mid-level model by minimizing the sum of the L 2 -norm of the measurement errors by area. The decision variables of this level were embedded into the lower-level model, which was formulated using a goal programming framework for which the primary goal was to minimize the impact of load shedding on the current system state. A traditional load shedding solution was built, which failed to respond quickly enough to prevent a cascading failure while the presented tri-level model was able to determine k shedding areas depending on the current operating condition via a spectral clustering approach, identify the severity of the event, and shed load to avoid a system islanding situation. Topics for future research. Areas of future work include adding a feedback control using the residual to correct the model to get a better estimate of the system states and solving the current overshooting issue as well.

Fig. 1 .

 1 Fig.1. The Tri-level model for automatic load shedding using linear programming and spectral clustering.

Fig. 2 .

 2 Fig. 2. The relative eigengap, γ k,rel , of L (left) and Ln (right).

Fig. 4 .

 4 Fig. 4. Scenario A: The complex power flow magnitudes over tie lines.

Fig. 5 .

 5 Fig. 5. Top: Frequency response for scenarios A and B. Bottom: Load shedding for scenario A (solid lines) and scenario B (dashed lines).

TABLE I SOLUTION

 I FOR k = 3 SHEDDING AREAS BASED ON THREE OPERATING CONDITIONS: MAXIMUM, MID, AND LOW DEMAND.

	Loading Level	Case A Demand = Max.	Case B Demand = 50%	Case C Demand = 33%
	Areas: Blue/Green/Red	A1	A2	A3	B1	B2	B3	C1	C2	C3
		6	7	1	4	6	1	6	1	3
		12	8	2	7	10	2	12	2	4
	Bus Number	13 11	9 10	3 4	8 9	11 12	3 5	13 -	5 10	7 8
		-	14	5	-	13	-	-	11	9
		-	-	-	-	14	-	-	-	14


Boualem Merainani
email: boualem.merainani@univ-eiffel.fr

Bian Xiong

Vincent Baltazart

Michael Döhler

Jean Dumoulin

Qinghua Zhang

Subspace-based modal identification and uncertainty quantification from video image flows

Keywords: operational modal analysis, structural health monitoring, computer vision, displacement extraction, uncertainty quantification

Vibration measurements for structural health monitoring (SHM) by operational modal analysis (OMA) are classically obtained from sensors that are embedded in or physically attached to the monitored structure, like accelerometers or strain gauges. However, the setup time of these sensors and their restricted number and space coverage limit their monitoring capabilities. Video image-based sensing methods can overcome these shortcomings. With adequate image processing methods, motion signals are extracted from video image flows, which are then processed by system identification methods to estimate modal parameters. In this way, the pixels in selected regions of interest within the images act as a dense network of contactless sensors distributed over the whole structure. In this paper, the efficiency of this video-based approach is demonstrated with laboratory experiments on a cantilever beam, in particular, by evaluating its capability for detecting weak damages mimicked by slight mass modifications. To this end, the steerable filter-based method (ST), that recovers displacements from local phase, is first extended to overcome its motion limitation of one pixel size. Then, the performance of the improved motion extraction method is compared with two other well-established methods in the context of OMA, where natural frequencies, damping ratios and mode shapes with high spatial resolution are estimated together with their uncertainty bounds using covariance-driven subspace identification. The compared methods are evaluated with the help of reference laser displacement measurements as well as a finite element model of the beam, revealing differences in the accuracy of the estimated mode shapes depending on the chosen method for motion extraction. Finally, aiming to investigate early structural damage detection, experiments are carried out under small structural changes and the results are compared to a reference state with the help of estimated uncertainties. Small but statistically significant changes in the modal parameters are detected, showing the potential of the vision based framework for SHM.

Introduction

Civil engineering infrastructures (e.g. bridges, buildings) are facing increasing challenges due to fast demographic growth, construction densification and climate change. Moreover, many structures exceed their planned lifespan. Therefore, health monitoring of these systems becomes increasingly important. Structural health monitoring aims to assess the structural behavior, to evaluate their performance and safety during the life cycle [START_REF] Farrar | An introduction to structural health monitoring[END_REF].

Operational modal analysis is a crucial part of vibration-based SHM under operating conditions by extracting modal parameters, which are essential for structural dynamic behavior analysis on the one hand, and for damage detection, localization and quantification on the other hand [START_REF] Magalhães | Vibration based structural health monitoring of an arch bridge: From automated OMA to damage detection[END_REF][START_REF] Brownjohn | Structural health monitoring of civil infrastructure[END_REF][START_REF] Avci | A review of vibration-based damage detection in civil structures: From traditional methods to machine learning and deep learning applications[END_REF][START_REF] Bhuyan | Vibration-based damage localization with load vectors under temperature changes[END_REF][START_REF] Döhler | Structural health monitoring with statistical methods during progressive damage test of S101 Bridge[END_REF]. It is therefore of great interest to provide a simple and affordable technology for reliable OMA. Traditional contact-based vibration measurement systems, e.g. by accelerometers or strain gauges, have practical shortcomings, such as the necessity of installing each of the sensors on the structure as well as the setup of a wired or wireless data transmission system. Furthermore, they can only offer a limited spatial coverage of a structure. In contrast, computer vision-based systems offer a promising alternative with a simpler setup and much higher spatial coverage and resolution. Nevertheless, the methods to process the vision data for OMA have not reached a full maturity yet and have shown some limitations, especially for small amplitude displacements under ambient excitation [START_REF] Feng | Computer vision for shm of civil infrastructure: From dynamic response measurement to damage detection-a review[END_REF]. This motivates methods for motion extraction with high accuracy on one side, and a systematic analysis of the accuracy in the OMA context on the other side. In particular, the influence of the chosen motion extraction method on the accuracy and precision of OMA is in general unclear, since a systematic evaluation is lacking so far. Besides, the performance of vision-based OMA with respect to traditional measurement techniques is rarely evaluated.

In this paper, a framework for vision-based OMA with high accuracy is proposed by combining improved image processing techniques with state-of-the-art system identification methods. Moreover, its capability to detect small structural changes is demonstrated. In addition to the proposed framework, an important contribution of the paper is a toolset that allows the analysis of the accuracy and precision of the vision-based OMA framework, and subsequently the comparison of performance for different motion extraction methods in the framework. The toolset is based on reference laser measurements, a finite element model of the investigated structure and computed uncertainty bounds that are tied to the used motion extraction and system identification methods.

This paper is organized as follows. Section 2 provides an overview of related works, followed by the objectives. Section 3 briefly presents the theory of image registration and three subpixel refinement approaches for motion extraction, as well as the subspace-based method for modal analysis. In Section 4, the experimental setup is described for vision-based measurements of a cantilever beam in the laboratory. The modal parameters evaluated based on the different approaches for motion extraction are presented in Section 5, where they are analysed with respect to their precision and accuracy, and finally their evolution under structural changes is presented for the most accurate approach in Section 6. Section 7 gives the final conclusions.

Related works and objectives

About the context to this work, Section 2.1 introduces existing vibration sensing technologies used for OMA or structural monitoring, including a brief review of vision system technology. Section 2.2 reviews methods used for motion extraction from video images and Section 2.3 details the objectives of the paper.

Vibration sensing technologies for OMA and SHM

OMA relies on vibration measurements of the monitored structure subject to natural excitations such as wind and road traffic. Conventionally, contact-based sensors are installed on the structure, e.g., linear variable-differential transformers (LVDT), accelerometers or strain gauges [START_REF] Wang | Full-scale measurements and system identification on sutong cable-stayed bridge during typhoon fung-wong[END_REF][START_REF] Park | The determination of bridge displacement using measured acceleration[END_REF][START_REF] Li | Health monitoring of joint conditions in steel truss bridges with relative displacement sensors[END_REF]. Such contact-based measurements are limited to a moderate number of points of the structure. Furthermore, the installation and maintenance of the sensors and the data acquisition system are often difficult and costly, in particular in a SHM context [START_REF] Feng | Computer vision for shm of civil infrastructure: From dynamic response measurement to damage detection-a review[END_REF]. Furthermore, the additional mass of sensors may induce changes in the dynamic behavior of some structures [START_REF] Zhang | A high-speed vision-based sensor for dynamic vibration analysis using fast motion extraction algorithms[END_REF].

Non-contact based vibration measurement techniques include Laser Doppler vibrometry (LDV) [START_REF] Nassif | Comparison of laser doppler vibrometer with contact sensors for monitoring bridge deflection and vibration[END_REF], terrestrial laser scanner [START_REF] Tzortzinis | Using 3d laser scanning for estimating the capacity of corroded steel bridge girders: Experiments, computations and analytical solutions[END_REF][START_REF] Bautista-De Castro | Integrating geomatic approaches, operational modal analysis, advanced numerical and updating methods to evaluate the current safety conditions of the historical bôco bridge[END_REF], GPS technology [START_REF] Meng | Detecting bridge dynamics with gps and triaxial accelerometers[END_REF][START_REF] Xu | Wind and structural monitoring of long span cable-supported bridges with gps[END_REF] and radar interferometry [START_REF] Negulescu | Comparison of seismometer and radar measurements for the modal identification of civil engineering structures[END_REF]. While offering alternatives to traditional measurement devices, they suffer from some drawbacks that still limit their wide application. Despite the high precision of LDV (higher than 0.05 mm) [START_REF] Uehan | Development of the u-doppler non-contact vibration measuring system for diagnosis of railway structures[END_REF][START_REF] Ri | Displacement measurement of concrete bridges by the sampling moiré method based on phase analysis of repeated pattern[END_REF], the price of such devices is very high and only displacements in the direction of the laser beam can be measured [START_REF] Gabbert | Vibration analysis of an electric wheel hub motor at stationary operating points[END_REF]. GPS sensor technology may show some limited accuracy for SHM application and/or low acquisition frequency [START_REF] Ribeiro | Non-contact measurement of the dynamic displacement of railway bridges using an advanced video-based system[END_REF][START_REF] Kumberg | A wireless sensor network using gnss receivers for a short-term assessment of the modal properties of the neckartal bridge[END_REF] .

Video cameras are an appealing technology for contactless measurements. Compared to traditional point measurement techniques, video cameras can provide full-field vibration measurements, virtually with each pixel acting as a sensor. They can range from industrial cameras with high resolution and sampling rate to consumer grade. The omnipresence of cameras, their flexible measurement distance and commonplace availability make them a promising candidate for widespread usage in SHM applications. A comprehensive review on this topic is given by Feng et al. [START_REF] Feng | Computer vision for shm of civil infrastructure: From dynamic response measurement to damage detection-a review[END_REF], and an overview on recent monitoring applications is given by Zona [START_REF] Zona | Vision-based vibration monitoring of structures and infrastructures: An overview of recent applications[END_REF].

The applied video-based technology depends on the camera system, with or without artificial targets, and on the applied image processing methods to extract the vibration data. For example, Luo et al. [START_REF] Luo | Robust vision sensor for multi-point displacement monitoring of bridges in the field[END_REF] adopted a high contrast artificial target vision-based approach for multi-point bridge displacement monitoring, and precision targets illuminated by infrared (IR) lamps or LED have been used in [START_REF] Ribeiro | Non-contact measurement of the dynamic displacement of railway bridges using an advanced video-based system[END_REF][START_REF] Tian | Remote bridge deflection measurement using an advanced video deflectometer and actively illuminated led targets[END_REF]. Other works take advantage from any natural texture in the structure as long as it has a strong contrast compared with the surrounding background. For example, Dorn et al. [START_REF] Dorn | Efficient full-field vibration measurements and operational modal analysis using neuromorphic event-based imaging[END_REF] proposed a target-free approach, which allows full-field vibration measurements and does not require any target installation on the investigated structure. The applied image processing methods for vibration data extraction are diverse and are reviewed in the following section. To validate the vision technology, comparisons to traditional displacement measurement technology are made, e.g., via laser vibrometers [START_REF] Chen | Modal identification of simple structures with high-speed video using motion magnification[END_REF] or LVDT sensors [START_REF] Ribeiro | Non-contact measurement of the dynamic displacement of railway bridges using an advanced video-based system[END_REF][START_REF] Feng | Nontarget vision sensor for remote measurement of bridge dynamic response[END_REF].

Motion extraction from video images

To extract motion signals from video images, most methods follow a two-stage approach. By maximizing the cross-correlation between two images (the reference and shifted images), the integer pixel shift is first estimated. Then, subpixel accuracy can be obtained by various refinement methods classified into three categories: gradientbased, correlation-based interpolation and local phase-based subpixel shifts.

In the first category, the optical flow method [START_REF] Lucas | An iterative image registration technique with an application to stereo vision[END_REF][START_REF] Horn | Determining optical flow[END_REF] has gained great importance for many applications. Chan et al. [START_REF] Chan | Subpixel motion estimation without interpolation[END_REF] proposed a simplified version of the optical flow. Based on this work, Javh et al. [START_REF] Javh | The subpixel resolution of optical-flow-based modal analysis[END_REF] have extracted subpixel motions of a vibrating steel beam and a cymbal, then performed modal analysis. Lu et al. [START_REF] Lu | Output-only modal parameter identification of structures by vision modal analysis[END_REF] performed structural modal parameter identification by vision modal analysis. Zhang et al. [START_REF] Zhang | A high-speed vision-based sensor for dynamic vibration analysis using fast motion extraction algorithms[END_REF] have combined the cross-correlation with modified Taylor approximation refinement to measure the displacements of a sound barrier on a bridge.

The second category is based on an interpolation in the vicinity of the cross-correlation peak to retrieve the subpixel shift [START_REF] Javh | The subpixel resolution of optical-flow-based modal analysis[END_REF][START_REF] Pan | Sub-pixel registration using quadratic surface fitting in digital image correlation[END_REF]. In [START_REF] Kim | Vision-based monitoring system for evaluating cable tensile forces on a cable-stayed bridge[END_REF], the authors have applied the digital image correlation technique to evaluate cable tensile forces on a cable-stayed bridge. In [START_REF] Thurman | Efficient subpixel image registration algorithms[END_REF], the authors have proposed to up-sample the cross correlation (UCC) between the actual image and the reference one by means of the fast Fourier transform (FFT) to locate its peak and to extract the subpixel shift. Satisfactory results were found in [START_REF] Feng | A vision-based sensor for non-contact structural displacement measurement[END_REF] when evaluating the displacements based on UCC and comparing them to high-performance laser measurements. An enhanced normalized cross-correlation method was proposed in [START_REF] Psarakis | An enhanced correlation-based method for stereo correspondence with sub-pixel accuracy[END_REF].

The last category relies on extracting the spatial local phase using the theory of steerable filters firstly introduced by [START_REF] Freeman | The design and use of steerable filters[END_REF]. The phase-based motion magnification algorithm developed by Wadhwa et al. [START_REF] Wadhwa | Phase-based video motion processing[END_REF] allows to highlight small variation in specific frequency bands. The method decomposes the video into amplitude and phase, at different scales, using the steerable pyramid filter bank. Then, the decomposed phases are band pass filtered, amplified and reconstructed back into video. Many researchers adopted the algorithm, e.g., for use in modal analysis [START_REF] Chen | Modal identification of simple structures with high-speed video using motion magnification[END_REF][START_REF] Zhu | A robust structural vibration recognition system based on computer vision[END_REF] or damage detection and localisation based on operational deflection shapes [START_REF] Civera | An experimental study of the feasibility of phase-based video magnification for damage detection and localisation in operational deflection shapes[END_REF].

Objectives

The main objective of this paper is to further complete and advance the video-based approach to SHM, by combining improved image processing techniques with state-of-the-art system identification methods, demonstrating the capability of such a framework to detect small structural changes. A particular focus is put on the evaluation of the accuracy and precision of the proposed framework in order to allow the comparison of different methods for extracting motion signals from video images in the context of OMA.

First, three methods for motion extraction -one from each of the categories above (see section 2.2) -are selected for a comparison. A constraint for motion extraction from the local phase with the steerable filter is the limitation of less than one pixel displacement in the image coordinates. Beyond this condition, the method breaks down [START_REF] Chen | Video camera-based vibration measurement of infrastructure[END_REF]. Since this condition may not always be satisfied for civil structures in operation, an algorithmic extension is proposed in this paper to extract displacements beyond the limitation of one pixel by combining the advantages of the ST with the zero-mean normalized cross-correlation (ZNCC). The approach is named EX-ST.

Second, the modal parameters are obtained by 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 computation of system matrices in subspace-based system identification[END_REF]. Compared to other methods for modal analysis, subspace methods have favorable properties in terms of accuracy, precision and implementation [START_REF] Döhler | Fast multi-order computation of system matrices in subspace-based system identification[END_REF][START_REF] Reynders | System identification methods for (operational) modal analysis: review and comparison[END_REF]. Moreover, they allow the computation of uncertainty bounds of the obtained modal parameter estimates [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][START_REF] Reynders | Uncertainty quantification in data-driven stochastic subspace identification[END_REF][START_REF] Greś | Uncertainty quantification of modal parameter estimates obtained from subspace identification: An experimental validation on a laboratory test of a large-scale wind turbine blade[END_REF]. These uncertainty quantification techniques are part of our proposed toolset to analyse the accuracy and precision of the vision-based OMA framework. In particular, this enables the evaluation of the precision of the modal parameter estimates obtained from the vision data, the comparison of different motion extraction methods within the framework, as well as the comparison with results from traditional measurement techniques. Up to the author's knowledge, this is the first time that uncertainty bounds of modal parameters are evaluated for vision data, showcasing the precision of the framework. The mode shapes obtained with such an approach are of particular interest due to their high spatial resolution, and thus their special usefulness, e.g. for finite element model updating or damage diagnosis. In this paper, it will be demonstrated that their accuracy depends on the accuracy of the selected method for motion extraction, which will be evaluated with the help of a finite element model. Finally, the vision-based framework is applied under small structural changes. Such changes are detected based on statistically significant changes in the obtained natural frequencies and mode shapes [START_REF] Greś | Uncertainty quantification of the modal assurance criterion in operational modal analysis[END_REF], for which the computed uncertainty bounds are again useful. All the developments of this paper are demonstrated with laboratory experiments on a cantilever beam.

Operational modal analysis using vision-based displacement measurement

This section presents the toolset which allowing to perform the vision-based modal analysis at the operational level. First and foremost, we give a presentation of three image processing methods to extract, along the structure, the motion signals from image flows. Motion signals are then processed by subspace-based identification method for modal analysis

Motion extraction from video frames

To extract motions from video images, the pixel shifts between two successive frames at time t and t + ∆t in the horizontal and/or vertical directions should be effectively revealed. Such shifts are often estimated through a twostep process, e.g. [START_REF] Zhang | A high-speed vision-based sensor for dynamic vibration analysis using fast motion extraction algorithms[END_REF], as illustrated in Fig. 1. First, the cross-correlation matrices between the reference image and every object image in the following video sequence are calculated. The coordinates of the maximum values in the cross-correlation matrices give the pixel level (integer) shifts, as shown in Section 3.1.1. Second, refinements through subpixel shifts estimation (δx, δy) are performed. In this work, three subpixel accuracy refinement methods are investigated and compared, one from each of the categories described in the introduction: the first one is based on Taylor approximation (Section 3.1.2), the second one is based on quadratic surface fitting (Section 3.1.3) and the last one is a phase-based method using steerable filters (Section 3.1.4). An algorithmic extension of the phase-based method to deal with displacements beyond the limitation of one pixel is proposed in Section 3.1.5.

Image registration using cross correlation

Commonly, the integer pixel displacement (∆x, ∆y) is estimated by locating the maximum value in the crosscorrelation matrix calculated between two image templates, of size X × Y, at two different times I(x, y, t) and I(x, y, t + ∆t) as [START_REF] Quam | Hierarchical warp stereo. in image understanding workshop[END_REF] (∆x * , ∆y * ) = arg max

-X≤∆x≤X -Y≤∆y≤Y CC(∆x, ∆y), (1)
where CC(∆x, ∆y) is the cross correlation defined as

CC(∆x, ∆y) = (x,y)∈T I(x -∆x, y -∆y, t + ∆t)I(x, y, t). ( 2
)
where T is the set of the pixel indexes of the image templates.

According to [START_REF] Sutton | Image correlation for shape, motion and deformation measurements: basic concepts, theory and applications[END_REF], the CC is sensitive to brightness and contrast variation. By normalizing Eq. ( 2) and subtracting the mean of both templates, the problem can be overcome. These operations gives rise to the zero mean normalized cross correlation given by ZNCC(∆x, ∆y) =

(x,y)∈T [I(x,y,t)-Ī(t)][I(x-∆x,y-∆y,t+∆t)-Ī(t+∆t)] (x,y)∈T [I(x,y,t)-Ī(t)] 2 (x,y)∈T [I(x-∆x,y-∆y,t+∆t)-Ī(t+∆t)] 2 (3)
where

Ī(t) = 1 XY (x,y)∈T I(x, y, t) (4)
Ī(t + ∆t) = 1 XY (x,y)∈T I(x -∆x, y -∆y, t + ∆t) (5)

Taylor approximation refinement

Once the pixel shift (∆x, ∆y) is estimated, the shift between I(x, y, t) and I(x + ∆x, y + ∆y, t + ∆t) is regarded as less than one pixel. Define x = x + ∆x and ỹ = y + ∆y, then the relationship between two consecutive templates can be written as [START_REF] Zhang | A high-speed vision-based sensor for dynamic vibration analysis using fast motion extraction algorithms[END_REF][START_REF] Chan | Subpixel motion estimation without interpolation[END_REF]

I(x, y, t + ∆t) = I( x + δx, ỹ + δy, t) ≈ I( x, ỹ, t) + ∂I( x, ỹ, t) ∂x δx + ∂I( x, ỹ, t) ∂y δy, (6)
which is the first order Taylor series approximation. The partial derivatives are estimated by finite difference between two images. In Eq. ( 6) there are two unknown parameters, δx and δy. The linear least-squares estimator minimizes the square error as ( δx, δy) = arg min (δx,δy)

Φ(δx, δy)

with

Φ(δx, δy) ≜ (x,y)∈T [I(x, y, t + ∆t) -I( x, ỹ, t) - ∂I( x, ỹ, t) ∂x δx - ∂I( x, ỹ, t) ∂y δy] 2 . ( 8
)
The solution to this standard least squares problem can be obtained by solving

∂Φ(δx, δy) ∂δx = 0, ∂Φ(δx, δy) ∂δy = 0 (9)
for δx and δy, which can be rewritten into

       x,y ( ∂I( x,ỹ,t) ∂x ) 2
x,y

∂I( x,ỹ,t) ∂x ∂I( x,ỹ,t) ∂y x,y ∂I( x,ỹ,t) ∂x ∂I( x,ỹ,t) ∂y x,y ( ∂I( x,ỹ,t) ∂y ) 2        δx δy =       x,y (I(x, y, t + ∆t) -I( x, ỹ, t)) ∂I( x,ỹ,t) ∂x x,y (I(x, y, t + ∆t) -I( x, ỹ, t)) ∂I( x,ỹ,t) ∂y       . ( 10
)
With the solution for refined subpixel level displacement, the total displacement can be written as

d x = ∆x + δx, d y = ∆y + δy. (11)

Quadratic surface fitting

Different from the Taylor approximation refinement method, the quadratic surface equation is fitted within a 3 x 3 neighborhood on the cross-correlation maximum matrix [START_REF] Pan | Sub-pixel registration using quadratic surface fitting in digital image correlation[END_REF] after estimating the pixel shift (∆x, ∆y). First, the zero-frequency component is moved to the center of the 2D DFT cross-correlation results. The quadratic surface fitted by the nine points can be written as

S (x, y) = a 0 + a 1 x + a 2 y + a 3 x 2 + a 4 xy + a 5 y 2 (12)
where x, y are relative coordinates respect to position of the maximum cross-correlation peak. The coefficients a i , i = 0 . . . 5, can be estimated through a surface model S (x, y) after nine points are substituted into the quadratic surface equation

          S (-1, 1) S (0, 1) S (1, 1) S (-1, 0) S (0, 0) S (1, 0) S (-1, -1) S (0, -1) S (1, -1)           =           CC(x i -1, y i + 1) CC(x i , y i + 1) CC(x i + 1, y i + 1) CC(x i -1, y i ) CC(x i , y i ) CC(x i + 1, y i ) CC(x i -1, y i -1) CC(x i , y i -1) CC(x i + 1, y i -1)           ( 13
)
where CC is the cross-correlation matrix between the template image and the object image, and (x i , y i ) are the coordinates of the maximum value in ZNCC. After obtaining the coefficients of the quadratic surface, the extreme point can be easily calculated by

δx = 2a 1 a 5 -a 2 a 4 a 4 2 -4a 3 a 5 , δy = 2a 2 a 3 -a 1 a 4 a 4 2 -4a 3 a 5 , (14)
and the total displacement follows as in [START_REF] Zhang | A high-speed vision-based sensor for dynamic vibration analysis using fast motion extraction algorithms[END_REF].

Phase-based displacement extraction using steerable filters

The steerable filter method have been proposed in [START_REF] Chen | Modal identification of simple structures with high-speed video using motion magnification[END_REF] to extract the displacement between two images from the phase change. In its original version, the method assumes limited displacement magnitude of the order of one pixel or smaller [START_REF] Chen | Video camera-based vibration measurement of infrastructure[END_REF][START_REF] Collier | Accuracy of phase-based optical flow for vibration extraction[END_REF]. As for the Taylor-based method, the steerable filter is based on the assumption of constant brightness of the scene and considers the change of pixel intensity related to image translation [START_REF] Chen | Modal identification of simple structures with high-speed video using motion magnification[END_REF][START_REF] Chen | Video camera-based vibration measurement for civil infrastructure applications[END_REF]. The video frames are convoluted with steerable quadrature filters (G θ + iH θ ) to determine, at each pixel, the local phase and amplitude. In this paper, the real part is selected to be the frequency response of the second derivative of a Gaussian (G 2 ) and the imaginary part is its Hilbert transform (H 2 ) [START_REF] Freeman | The design and use of steerable filters[END_REF]. In Fig. 2, the 9-tap steerable filter kernels at orientations 0 and π/2 are shown.

For an image sequence I(x, y, t), the local amplitude A and the local phase ϕ at the orientation θ are computed as

A θ (x, y, t)e iϕ θ (x,y,t) = (G θ 2 + iH θ 2 ) ⊗ I(x, y, t) (15)
where ⊗ is the convolution operator.

Under some assumptions, the local phases can be used to calculate the velocity in units of pixels, [START_REF] Chen | Video camera-based vibration measurement of infrastructure[END_REF][START_REF] Fleet | Computation of component image velocity from local phase information[END_REF][START_REF] Gautama | A phase-based approach to the estimation of the optical flow field using spatial filtering[END_REF], by

u = -( ∂ϕ 0 (x,y,t) ∂x ) -1 ∂ϕ 0 (x,y,t) ∂t , v = -( ∂ϕ π/2 (x,y,t) ∂y ) -1 ∂ϕ π/2 (x,y,t) ∂t , (16)
where u and v are the velocity in the x and y directions, respectively. Based on these expressions, the horizontal and vertical displacements δx and δy are directly obtained by [START_REF] Chen | Modal identification of simple structures with high-speed video using motion magnification[END_REF][START_REF] Chen | Video camera-based vibration measurement of infrastructure[END_REF] δx = -( ∂ϕ 0 (x,y,t) ∂x ) -1 (ϕ 0 (x, y, t) -ϕ 0 (x, y, t 0 )), δy = -(

∂ϕ π/2 (x,y,t) ∂y ) -1 (ϕ π/2 (x, y, t) -ϕ π/2 (x, y, t 0 )). ( 17
)
In fact, the obtained displacements are valid in the image regions with strong enough texture. Thus, instead of processing the entire images, the motion computation may be limited to some pixels, labelled as active pixels in the following. These pixels are usually selected on the first frame, from the extremum of the steerable filter output magnitude. These pixels form a mask which is then applied on the whole image frames for motion extraction.

Extension of the local phase-based method to multi-pixel displacement

When the motion magnitude of the structure is beyond one pixel, the phase signal goes beyond a change of π in phase, and the method breaks down. One common adopted solution is to down-sample the video, before being filtered, in order to recover the subpixel condition limit, which however reduces the image contrast and smooths the image texture. As a consequence, the number of active pixels is dramatically reduced. Furthermore, high frequency information may be lost when down-sampling [START_REF] Xiao | Invertible image rescaling[END_REF]. To address this limitation, we propose to relate the steerable approach to the general two-steps process, presented in Section 3.1.2, in order to determine multi-pixel displacement without down sampling, as follows. After defining the mask, the ZNCC given in Section 3.1.1 is employed to estimate the integer pixel shift (∆x, ∆y) between the reference image I(x, y, t 0 ) and the actual image I(x, y, t). If a pixel shift exists, the following steps are taken: 1 Shift the mask by (∆x, ∆y). 2 To compute the phase change on the current frame, Replace temporarily the reference phases ϕ 0 (x, y, t 0 ) and ϕ π/2 (x, y, t 0 ), computed on the first frame and given in Eq. ( 17) by the phases of the previous frame at time t -∆t, i.e I(x, y, t -∆t), which are ϕ 0 (x, y, t -∆t) and ϕ π/2 (x, y, t -∆t).

3 Accumulate the estimated subpixel motion as

d x (t) = δx(t) + d x (t -∆t), d y (t) = δy + d y (t -∆t)
Finally, when the pixel shift is zero or comes back to zero, the reference phases are reconsidered as ϕ 0 (x, y, t 0 ) and ϕ π/2 (x, y, t 0 ). The flowchart of the proposed approach is shown in Fig. 3.

Operational modal analysis with subspace-based system identification

Assume that the vibration behavior of the investigated structure can be modelled by a linear time-invariant system, then its dynamics can be described by the discrete-time state space model [START_REF] Reynders | System identification methods for (operational) modal analysis: review and comparison[END_REF]

x k+1 = Ax k + w k y k = Cx k + v k , ( 18
)
where k is the integer time step corresponding to the system at time t = k∆t, y k ∈ R r is the output vector containing the computed displacements d x (k∆t) and d y (k∆t) for the r selected templates or active pixels of the video, and x k ∈ R n is the state vector. The state noise w k ∈ R n is related to the unknown ambient excitation, and vector v k ∈ R r describes the measurement noise. The modal parameters are related to the eigenvalues and eigenvectors (λ i , ϕ i ), i = 1, . . . , n, of A and to C by where µ i is an eigenvalue of the corresponding continuous-time system, f i is the natural frequency, ζ i is the damping ratio and φ i is the mode shape at the output coordinates. To estimate the system matrices A and C from the output data y k of length N, k = 1, . . . , N, and consequently the modal parameters in [START_REF] Ri | Displacement measurement of concrete bridges by the sampling moiré method based on phase analysis of repeated pattern[END_REF], the reference-based covariance-driven subspace algorithm [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] is used, and the standard deviations of the modal parameter estimates are computed as in [START_REF] Döhler | Efficient multi-order uncertainty computation for stochastic subspace identification[END_REF].

µ i = log(λ i ) ∆t , f i = |µ i | 2π , ζ i = -Re(µ i ) |µ i | , φ i = Cϕ i , (19)

Experimental setup for validations

The performance of the different methods for motion extraction and the efficiency of the video-based approach to SHM were verified by a set of experiments conducted on a cantilever beam structure of size 30 mm × 2000 mm × 5 mm mounted on a shake table in the Structures and Integrated Instrumentation Laboratory at Gustave Eiffel University.

Beam model

The finite element (FE) model of the beam will mainly be used for the comparison of the analytic mode shapes with the estimated ones from data. The steel cantilever beam is modelled using the Timoshenko beam theory. The model is adopted from [START_REF] Przemieniecki | Theory of matrix structural analysis[END_REF] and illustrated in Fig. 4 (left). The Young's modulus, mass density and the Poisson's ratio of the beam material is 2.1 • 10 11 Pa, 7850 kg/m 3 and 0.21, respectively.

The model is composed of 160 beam elements, which will conveniently correspond to the number of regions of interests (ROIs) into which the video frames will be divided. Every node has three rotational and three translational degrees of freedom.

The beam model is also used to define positions of additional masses that will be used to emulate different small damages in the experiments. With a preliminary analysis, eight positions and an additional mass of 0.044 kg (corresponding to 1.86% of the beam mass) were chosen, such that changes in the modal parameters are moderate. The chosen positions are designated by their distance in mm from the support as as P2000, P1803, P1507, P1277, P949, P686, P456 and P260, and are illustrated in Fig. 4 (right).

Experimental setup

The camera used for recording the video flows was a PHANTOM MIRO 320S, equipped with a CMOS sensor with resolution 1920 × 1200 and mounted with a NIKKOR lens with a fixed 50 mm focal length.

The camera was set on a tripod. The effective distance from the camera to the beam was fixed such that the full beam length will be captured. The camera was used in portrait mode to take advantage of the full 1920 pixels along the beam in vertical direction.

The number of pixels for the horizontal direction was adjusted to 136 pixels to gain space in the camera's dynamic random access memory, while ensuring that the beam remains in the scene when vibrating. The frame rate was chosen as 512 frames per second, ensuring a good identification of the first six modes. The exposure time was selected relatively short, so that the camera can collect enough light. Three additional lights (Halogen Spotlight, MAZDA CORMORAN type of 1500 W) and a white background were used to provide optimal brightness and contrast conditions. Three black and white artificial targets are mounted at three particular positions on the beam (top (1/1), middle (1/2) and 50 cm from the base (1/4)) at the positions of three high-accuracy laser displacement sensors (LDS) for comparison between displacements extracted from vision measurements and displacements measured by LDS. Two Wenglor (top and middle) sensors and one Keyence (bottom) sensor were used, whose specifications are given in Tab. 1 based on the data sheets, except for the values for the repeatability that are experimentally obtained when the beam was at rest. The experimental installation is shown in Fig. 5.

For data collection with the camera and the LDS, the shake table (Model: NEBS GR-63, by team vibration testing equipment) was driven in random mode horizontally, perpendicular to the camera axis. The excitation noise was band limited between 1 Hz and 200 Hz. Datasets of 33 seconds length were recorded without and with the added mass at the different positions shown in Fig. 4 (right). For processing with the different motion extraction methods, the video frames were separated to 160 ROIs, in total, of size 12 × 37 pixels. The horizontal displacements from each of the ROIs were estimated by the three introduced methods, namely Taylor approximation (Section 3.1.2), quadratic surface fitting (QSF, Section 3.1.3) and the proposed EX-ST (Section 3.1.5) that all share the ZNCC technique (Section 3.1.1) as a prior step. The computation time per frame of each algorithm is shown in Table 2, indicating that EX-ST takes

Evaluation of motion extraction methods in OMA context

In this section, the different motion extraction algorithms are evaluated in the context of modal analysis using data from the beam in the reference condition (without added mass), and compared to results from the laser data. First, a direct comparison of the extracted displacement data to the laser data is made in terms of the root mean squared error. Second, the estimated natural frequencies and damping ratios and their uncertainty bounds obtained with the covariance-driven subspace algorithm are evaluated. Finally, the quality of the mode shape estimates is evaluated for the different methods based on the FE mode shapes.

Direct comparison to laser data

Since time synchronization between the laser sensors and the camera was not possible, the laser data was synchronized a posteriori with the image flows by using cross correlation between signals. Furthermore, the extracted displacements from the video were converted from pixels into millimeters by using the pinhole camera model through scaling the distance between two targets in the world coordinate with its pixel distance in the image coordinate. The scaling factor was, thus, estimated to be 1.007 mm/pixel.

In Fig. 6 the motion signals extracted from the vision and laser systems are shown at the upper laser sensor (position 1/1), corresponding to ROI 2 in the motion extraction. The signals are represented on top of each other for comparison. It is observed that the displacements obtained from any of the three motion extraction algorithms seem to match very well with the reference laser data for the whole dataset of 33 s length. To quantify the differences between the vision and the laser data, the root mean squared error 3, where it is indicated that the proposed EX-ST method provides a lower RMSE compared to ones obtained with the QSF and Taylor methods.

RMSE =        1 N N k=1 y (vision) k -y (LDS) k 2        1/2 (20)

Operational modal analysis

The extracted displacements from the vision data, using the EX-ST, QSF and Taylor algorithms, as well as the LDS data, have been used to obtain the respective sets of modal parameters for each of the algorithms. The referencebased covariance-driven subspace method [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] has been applied for this task, together with the related uncertainty quantification to obtain the confidence intervals of the estimates [START_REF] Döhler | Efficient multi-order uncertainty computation for stochastic subspace identification[END_REF]. With the considerable number of r = 160 available virtual sensors from the vision data, three reference sensors were selected for the analysis, namely the ones coinciding with the laser sensor positions. Since the first mode of the structure is expected around 1 Hz and the sampling frequency is at 512 Hz, the number of time lags for the subspace algorithm needed to be set relatively high, and altogether 300 lags have been used for the computation of the covariance-driven Hankel matrix. The identification has been carried out at model orders 1, . . . , 100, and the resulting stablization diagram of the natural frequencies is shown together with uncertainty bounds and the power spectral density in Fig. 7 for the laser data and the vision data from the proposed EX-ST algorithm. Since physical modes remain stable and have much lower uncertainty bounds than spurious modes due to noise, a threshold of 1.5% was put on the coefficient of variation of the frequencies (standard deviation divided by frequency), leading to a clear and easily interpretable diagram. Obviously, six stable modes are clearly distinguished. Their quality is analysed in the following subsections with respect to the underlying motion extraction methods.

Precision of the natural frequencies and damping ratio estimates

The estimated frequencies and damping ratios and their ±2σ confidence intervals, obtained from the vision data with each of the three motion extraction methods as well as the LDS, are stated in Tab. 4 and visualized in Fig. 8. The confidence intervals are in absolute terms and have the same units as the frequencies and damping ratios, respectively. Figure 8: Natural frequencies and damping ratios with their ±2σ confidence intervals obtained from the LDS data and the from the vision data using the EX-ST, QSF and Taylor algorithms.

It can be noted that the first mode has the highest uncertainties, especially its damping ratio. This is due to higher uncertainties in damping estimation in general [START_REF] Gersch | On the achievable accuracy of structural system parameter estimates[END_REF][START_REF] Bernal | First mode damping ratios for buildings[END_REF], and due to the difficult estimation of this mode in particular, which is at a considerably lower frequency than all other modes while the dataset is relatively short. For all six modes the results show close correspondence although the modal parameter estimates vary slightly. The confidence intervals obtained from the different datasets overlap, and each of the frequencies and damping ratios is within all the confidence intervals of the other three estimates, showing that the frequency and damping ratio estimates obtained from the three motion extraction methods as well as the LDS are statistically coherent. Moreover, the size of the confidence intervals is similar for the most part for the different methods, with a bit more variation for some of the damping estimates which is naturally to be expected. This shows that the precision of the modal parameter estimates is also similar for the different methods; hence the choice of the motion extraction method does not seem to have a significant influence on the natural frequency and damping ratio estimates.

Accuracy of the vision-based mode shapes

Fig. 9 compares the experimentally obtained mode shapes determined by the LDS and vision data using the three motion extraction methods with those obtained from the FE model. At first glance, the vision-based mode shapes follow a similar pattern as the numerical ones and have close correspondence. Those obtained from the LDS suffer from low spatial resolutions with only three measurement points, while the vision-based mode shapes have a high Table 4: Natural frequency and damping ratio estimates with their ±2σ confidence intervals.

Mode

Laser sensor EX-ST Taylor QSF A close examination, however, reveals that errors in motion estimation have a direct impact on the identified mode shapes. For instance, the mode shapes computed from the motion data from the Taylor approximation method present more fluctuations that are visible on the beam top, as can be seen in the zoom in Fig. 9. Those computed based on QSF are smoother but the gap is larger near the anti-node positions. The mode shapes obtained based on the proposed EX-ST method seem to match the best with the FE predicted ones.

f (Hz) ζ (%) f (Hz) ζ (%) f (Hz) ζ (%) f (Hz) ζ (
These results are confirmed by a comparison of the modal assurance criterion (MAC) between the 6 numerical and the vision-based mode shapes obtained from the three different motion extraction methods in Fig. 10, where the MAC between two complex-valued mode shape vectors φ and ψ is defined by [START_REF] Greś | Uncertainty quantification of the modal assurance criterion in operational modal analysis[END_REF]

MAC(φ, ψ) = |φ H ψ| 2 φ H φψ H ψ . (21)
The MAC value is bounded between 0 and 1, and values near 1 indicate a strong similarity between both mode shapes.

In Fig. 10 it can be seen that all MAC values are very high (> 0.99); however the mode shapes based on the proposed EX-ST algorithm stand out among the others and offer the highest MAC values among the different motion extraction methods for all modes. Consequently, the proposed EX-ST algorithm was adopted for further analysis in the next section.

Application to detect structural changes through added mass

Finally, the vision-based technology is evaluated in a monitoring context, where a mass is attached at different positions to the beam (see Fig. 4, right) to emulate different damages. Based on the chosen EX-ST algorithm, the modal parameters are estimated for each structural state and compared to the reference state.

For the modal parameter comparison in the different structural states, a simple approach is adopted from [START_REF] Greś | Uncertainty quantification of the modal assurance criterion in operational modal analysis[END_REF] based on the natural frequencies and mode shapes together with their uncertainty bounds to infer if changes with respect to the reference are statistically significant. To this end, let the reference modal parameters ( f ref i , φi ref ), i = 1, . . . , 6, with their uncertainties (frequency standard deviation σ ref f i , mode shape covariance Σ ref φ i ) be estimated in the mass-free condition, and let the modal parameters ( fi , φi ) with their uncertainties (σ f i , Σ φ i ) be computed from the vision test data under a structural change. Then, a statistically significant change in frequencies or mode shapes can be inferred for a confidence level of 95% as follows:

• Frequency change: The confidence bound for each frequency in the reference state is

[ fi ref -2σ ref f i , fi ref + 2σ ref f i ],
since frequency estimates are approximately Gaussian distributed. If fi is not in this interval, then label it as changed.

• Mode shape change: The MAC equals one for identical mode shapes. Due to statistical uncertainties it is only close to one when considering estimates of identical mode shapes. To check the hypothesis that φref i and φi correspond to the same, unchanged mode shape, MAC( φref i , φi ) and its uncertainty is computed. As shown in [START_REF] Greś | Uncertainty quantification of the modal assurance criterion in operational modal analysis[END_REF], the MAC distribution is not Gaussian under the hypothesis of identical mode shapes, but approximated by a scaled and shifted χ 2 . Then, the quantile t MAC of the distribution can be computed for the chosen 95% confidence level based on the mode shape covariances Σ ref φ i and Σ φ i , leading to the MAC confidence bound [t MAC , 1] as detailed in [START_REF] Greś | Uncertainty quantification of the modal assurance criterion in operational modal analysis[END_REF]. If the computed MAC is outside this interval, i.e. MAC( φref i , φi ) < t MAC , then label φi as changed.

The results of the monitoring using vision data sets of all beam conditions (mass-free reference condition (Ref), and with the added mass at the positions indicated in Fig. 4) are shown for each of the modes in Fig. 11. For each mode, it is evaluated if the frequency estimate lies within the confidence bound of the reference, as outlined above, and a change in the corresponding mode shape as indicated by the confidence bound of the MAC is marked by an empty circle, while no significant mode shape change is marked by a full circle.

First of all, it can be observed that the structural changes induce changes in some of the frequencies and in some of the mode shapes, but not in all of them. For example, there is no significant change in the first mode shape of mode 1, however the frequency changes are significant for some of the mass positions. On the other side, the frequencies of mode 3 show no significant change for the added mass at P1803, P949 and P260 whereas the corresponding mode shapes do change. This also shows the importance of evaluating changes both in the frequencies and in the mode shapes, since some damages do not induce change in the frequency but in the mode shape, or vice versa.

From Fig. 11 it can also be observed that for most of the modes the mode shapes change for most of the damage scenarios. It should be noted that only the statistical evaluation of the MAC revealed that these changes are actually significant, since the MAC values are relatively high, being above 0.97. For the cases where no mode shape change was indicated, a numerical study actually confirmed that the change in the numerical mode shapes is insignificant with mostly MAC > 0.999. To illustrate the actual mode shape changes, Fig. 12 shows the mode shape estimate of mode 4 together with its confidence bound in the reference condition, as well as the mode shape estimates after adding the mass at P2000 and at P1803. It can be seen that the estimate corresponding to P1803 remains in the confidence bound of the reference, and indeed no mode shape change is indicated in Fig. 11. In contrast, the mass attached in P2000 has significantly modified the mode shape estimate that is leaving the confidence bound of the reference, and it is correctly classified as changed. In summary, the structural changes through the added mass could be detected for all of the scenarios in at least some of the modes with the proposed vision technology, either by frequency or mode shape changes, or both. Finally, it should be noted that the change of the frequency and/or MAC of a single mode beyond their reference confidence interval is sufficient to make a decision about the beam health.

Conclusion

In this paper a video-based approach for operational modal analysis was presented and evaluated on experimental data of a cantilever beam. The approach combines an improved motion extraction method with state-of-the-art subspace identification and uncertainty quantification. For motion extraction, the proposed extension of the local phase method that relies on steerable filters (EX-ST) is able to consider both small and large motions in the source videos. Furthermore, two widespread reference motion extraction algorithms namely Taylor and QSF were used for comparison. The observations of this study can be summarized as follows:

• The displacements extracted with the proposed EX-ST method achieve a better accuracy than Taylor approximation and QSF, at the expense of a moderately increased computation time.

• The application of the covariance-driven subspace method with uncertainty quantification allows to evaluate the performance of the vision-based modal parameters estimates. It was concluded that estimates obtained from the displacements extracted from 160 ROIs with motion extraction algorithms showed no statistically significant differences compared to those identified from the reference LDS, and that confidence bounds were similar.

• The identified mode shapes from the vision data have a very high spatial resolution and closely match the numerical ones obtained from the FE model. The mode shape estimates obtained from EX-ST motions were found smoother and closest to the numerical ones, which was shown by the MAC estimates, while the estimates based on the Taylor approximation and QSF methods were less accurate.

• The proposed framework was able to detect small structural changes by an added mass at different positions of the beam. Statistically significant changes in the frequencies as well as in the high resolution mode shapes were found, considering the estimation uncertainty related to the vision data, showing the potential of the framework for monitoring.

The developments in this paper show that a thorough evaluation of the motion extraction method as well as the consideration of statistical uncertainties is important to ensure good accuracy and precision in video-based operational modal analysis. Such an evaluation has been made in a controlled laboratory environment as a necessary validation step before actual application to real-world civil structures in operation.

After this successful evaluation with laboratory experiments, the presented approach will be tested on a bridge and other outdoor structures in the next step, so that its performance will be evaluated under real world operation conditions. Our future researches will address additional challenges in outdoor environments, including ambient light variations over time and atmospheric distortion, to ensure robust motion extraction for vision-based operational modal analysis.

Figure 1 :

 1 Figure 1: Illustration of overall-shift template image.

Figure 2 :

 2 Figure 2: Representation of the Filters used to compute local phase and local amplitude, (a) real horizontal(G 0 2 ), (b) imaginary horizontal(H 0 2 ), (c) real vertical (G π/22 ), (d) imaginary vertical(H π/2 2 ).

Figure 4 :

 4 Figure 4: Illustration of cantilever beam FE model (left); choice of positions for adding mass for damage detection experiments (right).

Figure 5 :

 5 Figure 5: Experimental setup: (a) video cameras in the foreground, the cantilever beam in the background, (b) closer view of the cantilever beam set, (c) raw frame from the video on which the regions of interest are indicated. Three of them are represented on the right hand side.

Figure 6 :

 6 Figure 6: Comparison of the beam displacements at the top (position 1/1) measured by the LDS and extracted from the vision data, with zoom on a time interval.

Figure 7 :

 7 Figure 7: Stabilization diagrams with uncertainty bounds (horizontal bars) of the natural frequencies, together with power spectral density. (a): Obtained from laser data, (b): Obtained from vision-data with EX-ST algorithm.

Figure 9 :

 9 Figure 9: Comparison between the identified mode shapes with those from FE model

Figure 10 :

 10 Figure 10: MAC values between vision-based mode shapes based on the three motion extraction methods and numerical mode shapes from FE model.

Figure 11 :

 11 Figure 11: Changes of the natural frequencies related to the first six modes. The green dashed horizontal lines represent the 95% confidence interval of the reference natural frequency

Figure 12 :

 12 Figure 12: Changes in the fourth mode shape due to added mass on two different positions.

Table 1 :

 1 Specifications of the LDS.

			Model
	Features	Keyence IL-030	Wenglor YP11MGVL80
	Measuring range	20 to 45 mm	50 to 100 mm
	Reference distance	30 mm	50 mm
	Wavelength of light source	655 nm	655 nm
	Resolution	n/a	25 µm
	Repeatability	1 µm	6 µm
	Camera
	Light
	Beam
	Shake table

Table 2 :

 2 Computation time to extract displacements with the EX-ST, QSF and Taylor algorithms.

	EX-ST	QSF	Taylor
	Time per frame (ms) 1.01	0.53	0.48
	about twice the time than Taylor or QSF.

Table 3 :

 3 RMSE of the extracted displacements with the EX-ST, QSF and Taylor algorithms with respect to the reference laser data.

					RMSE (mm)			Motions extraction algorithms
					Position			EX-ST	QSF	Taylor
					(1/1)				0.0404	0.0430	0.0461
					(1/2)				0.0342	0.0348	0.0422
					(1/4)				0.0438	0.0451	0.0481
	0	10	20	30	40	50	60	70	80	90	100
					Frequency (Hz)

  Reference f with the corersponding 95% CI f, where mode shapes are classified as reference based on MAC uncertainty f, where mode shapes are not classified as reference based on MAC uncertainty

		1.02	Mode 1	6.8	Mode 2
				6.75
		1		6.7
	(Hz)	0.98	(Hz)	6.65
	f		f
				6.6
		0.96
				6.55
		0.94		6.5
		Ref	P2000 P1803 P1507 P1277 P949 P686 P456 P260	Ref	P2000 P1803 P1507 P1277 P949 P686 P456 P260
			Mode 3
		19
			37.2
	f (Hz)	18.8	36.8 37
		18.6	36.6
			36.4
		18.4	36.2
		Ref	P2000 P1803 P1507 P1277 P949 P686 P456 P260

Acknowledgements

Part of this work was funded by the Research Fund for Coal and Steel under grant agreement No 800687 in the framework of the H2020 DESDEMONA project. The authors are grateful to Ivan Guéguen, Jean Philippe Gourdon and Fabrice Blaineau for the assistance provided for the experimental set up.