Machine Learning Clustering Techniques to Support Structural Monitoring of the Valgadena Bridge Viaduct (Italy)

Abstract

The lack of precise and comprehensive information about the health of bridges, and in particular long span ones, can lead to incorrect decisions regarding maintenance, repair, modernization, and reinforcement of the structure itself. While the consequences of inadequate interventions are quite apparent, incorrect decisions can also result in unnecessary or misdirected actions. For example, an inadequate assessment of the structural health can lead to the modernization and replacement of some components that are still sound. Structural Health Monitoring (SHM) involves the use of time series derived from periodic measurements of the structure’s behavior, considered in its operational and load environment. The goal is to determine its response to various solicitations and, in particular, to highlight any critical issue in the structure’s behavior that may affect its reliability and safety due to anomalies and deterioration. This paper proposes an SHM method applied to the Valgadena bridge, one of the tallest viaducts in Italy and Europe (maximum height 160 m), located on the Altopiano dei Sette Comuni in the Province of Vicenza. Despite the fact that the viaduct itself had already been monitored during its construction using classical geometric leveling techniques, the methodology proposed here is based instead on the use of affordable dual-frequency GNSS (Global Navigation Satellite System) receivers to determine static and dynamic components of the bridge movements. Specifically, an effective combination of time series analysis methods and machine learning techniques is proposed in order to determine the vibration modes of the monitored viaduct. Monitoring is performed in regular operation conditions of the bridge (operational modal analysis (OMA)), and the use of certain machine learning methods aims at supporting the development of an effective automatic OMA procedure. To be more specific, the random decrements technique is used in order to make the vibration characteristics of the collected signals more apparent. Time-domain-based subspace identification is applied in order to determine a proper model of the collected measurements. Then, clustering methods, namely DBSCAN (Density-Based Spatial Clustering of Applications with Noise) and GMMs (Gaussian Mixture Models), are used in order to reliably estimate the system poles, and hence the corresponding vibration characteristics. The performance of the considered methods is compared on the Valgadena bridge case study, showing that the use of GMM clustering reduces, with respect to DBSCAN, the impact of the choice of certain parameter values in the considered case.

1. Introduction

Bridge monitoring is a fundamental practice for ensuring safety, functionality, and longevity of bridge structures. The importance of such practice has became even more apparent after the tragic collapse of the Polcevera viaduct in Genoa, Italy, better known as the Morandi Bridge, which occurred on 14 August 2018: inspection and monitoring of bridges, in particular long span ones, such as cable-stayed and suspension bridges, are crucial activities in order to schedule proper maintenance interventions and prevent tragic events. To such aim, geomatic monitoring techniques encompass a diverse array of methodologies aimed at comprehensively assessing various aspects of bridge behavior and condition. Among these techniques, geometric leveling stands as one of the traditional yet highly effective methods for monitoring structural deformations and movements. Geometric leveling involves the precise measurement of vertical height differentials between points using optical instruments, namely levels and levelling rods (also called “level staffs”) [1]. By comparing repeated measurements over time, engineers can detect subtle changes in the bridge’s vertical alignment, which may indicate structural deformations or settlements. While geometric leveling provides accurate results, its application can be labor-intensive and time-consuming, especially for large-scale bridge networks.

To address these challenges, Structural Health Monitoring (SHM) takes advantage of sensor data, provided by a properly designed monitoring system, in order to supplement and substitute, if possible, long human-made inspections [2]. A wide range of technologies have already been considered for such aim, ranging from more classical ones, such as strain gauges [3] and accelerometers [4], fiber optic sensors [5], to modern geomatic monitoring approaches, such as Global Navigation Satellite System (GNSS) [6], vision-based methods [7], also for degradation monitoring [8], Terrestrial Laser Scanning (TLS), and RADAR (Radio Detection and Ranging, and, in particular satellite InSAR - Interferometric Synthetic Aperture Radar [9]), where each of such technology can be appropriate for assessing different kinds of information, such as strains, displacements and deformations, and vibrations [10]. Strain gauges are traditionally used in order to assess material deformations, but vibrating wire strain gauges can be used to determine vibration characteristics as well. Since accelerometers do not make it possible to determine absolute object positions, they are typically used for analysis of local position changes, e.g., vibration-based monitoring. GNSS technology enables the continuous real-time monitoring of bridge movements and deformations by tracking the positions of GNSS receivers installed on the structure [11]. TLS facilitates high-resolution 3D mapping of bridge surfaces, allowing for detailed analysis of structural conditions and detecting localized deformations or damages [12]. Additionally, remote sensing techniques, including satellite and aerial imagery, complement geomatic monitoring efforts by providing broad-scale assessments of bridge infrastructure and surrounding environments [13]. In particular, InSAR recently proved to be an effective tool for monitoring bridge deformations [14]. Together, these monitoring methods offer a comprehensive toolkit for bridge health assessment, enabling proactive maintenance, early detection of potential hazards, and informed decision-making for infrastructure management and safety enhancement [15]. SHM is a very popular research topic, and as a comprehensive review of SHM approaches is beyond the scopes of this paper, the reader is referred to [16,17,18,19,20,21,22], and more specifically to [23,24,25,26,27,28,29] for the SHM of bridges for a more complete overview of the state-of-the-art SHM methods.

Either model-based or data-driven monitoring techniques are usually implemented [30,31]: the first aim is to properly describe the physical system with a proper finite element model, which however often may either lead to too complex and computationally inefficient models, or require the introduction of some simplifying assumptions that, on the other hand, may reduce the overall reliability of the system. Differently, data-driven methods aim at directly inferring the system characteristics from the collected data, typically leading to reliable results when the amount of available data is sufficiently large.

The quite recent success of machine (and deep) learning methods for data analysis in a number of applications motivated their usage also for SHM: for instance, machine and deep learning-based analysis of multi-sensor time series, e.g., using multiple accelerometers, has been investigated in [32,33] for damage classification, whereas vision-based methods have also been implemented for surface damage detection [34,35] and for determining the structure vibrations even at the sub-millimeter level [36,37]. It is worth noting that, despite the fact that in principle the usage of machine learning methods is not constrained to neither model-based nor data-driven-based approaches, the typical large data availability associated to data-driven methods make them quite ideal for the use of this kind of methods, as is also considered in this work. The reader is referred to [31,38,39] and references therein for more complete overviews on the use of machine learning methods in SHM applications.

In the specific case of SHM for viaduct monitoring using GNSS techniques, as considered in this work, some feasibility studies have already been conducted, such as in [40], wherein high-frequency receivers have been deployed, or integrating GNSS with other sensors, such as accelerometers [41,42], while some have developed studies on the use of mobile GNSS receivers [43]. Recently, some investigations have been dedicated to the feasibility of employing low-cost GNSS receivers [44,45], while, finally, certain studies have delved into the examination of algorithms tailored to the specific application of quality control data for bridge monitoring [46].

The movements of a bridge–viaduct are characterized by static/semi-static and dynamic components. The static component remains unchanged over time, while the semi-static component identifies the trend over a medium-length time interval. The short-term component comprises dynamic displacements and noise, wherein the former rapidly changes over time as a consequence of a physical process, and the latter is typically randomly distributed, with some possible specific characteristics due to the adopted measurement method and instruments. GNSS is frequently used to assess static and dynamic movement components, even in real time, while conventional monitoring using accelerometers and extensometers does not make it possible to determine static and semi-static displacements [47]. Estimating motion with geotechnical sensors like accelerometers presents challenges such as contamination with errors like drift and their capability to measure only relative displacements of a structure. Advances in GNSS sensors and processing algorithms mitigate positioning errors; integrating GNSS with supplementary sensors like accelerometers and extensometers improves position estimation accuracy. High-frequency GNSS sensors (>1 Hz) can be used, for instance, to exploit the RTK (Real-Time Kinematic) positioning approach; however, the obtained measurements are typically affected by errors and noise characterized by different statistical distributions (colored noise and white noise), which must be filtered out before displacement and deformation monitoring. This has been achieved in previous works by implementing time series analysis techniques, allowing for the extraction of semi-static and dynamic movement components [48].

In this work, we adopt a modal analysis approach; that is, we focus on the study of the vibration characteristics of the bridge: the aim is that of properly formulating a methodology to assess the main vibration characteristics of the bridge, which are supposed to change in case of damage; hence, monitoring the stability of such vibration characteristics over time shall make it possible to determine when/whether maintenance or human intervention is needed. To be more specific, this work focuses on the determination of bridge vibration characteristics on a specific dataset, in order to determine the standard vibration behavior of the bridge, whereas the application of the implemented methodology for long-time monitoring will be considered in our future works.

From an operative point of view, we consider an effective combination of different processing techniques of GNSS data in order to provide a reliable modal identification to properly characterize the bridge movements, and, consequently, enable vibration-based SHM by the detection of changes in the GNSS principal modes. Despite the fact that, for several mechanical systems, proper identification of a system’s vibration characteristics can be implemented in controlled operating conditions, e.g., in laboratory tests, where both excitation signals and the system output are available, this is usually hard to replicate for periodic monitoring of several civil infrastructures, such as bridges. Hence, it is important to develop infrastructure monitoring procedures that can be effective even during their regular working conditions, by exploiting ambient vibration, e.g., excitation inputs induced by natural or human causes during the standard operations of the infrastructure. This kind of monitoring technique is usually referred to as Operational Modal Analysis (OMA) [49].

To be more specific, this paper focuses on Automatic OMA (AOMA), i.e., on procedures to automatically determine and analyze the modal characteristics of a certain infrastructure from data collected during its regular operating conditions [50,51,52,53]. Similarly to [54], this paper focuses on the use of clustering techniques (i.e., DBSCAN (Density-Based Spatial Clustering of Applications with Noise [55,56])) in order to reliably determine the most representative vibration frequencies of an infrastructure. However, differently from [54], this work investigates the use of DBSCAN on different parameters characterizing the system and evaluates different alternatives in order to properly assess the uncertainty on the estimated mode frequencies, which can be of major importance for the final aim of ensuring good monitoring performance over time. Furthermore, the use of Gaussian Mixture Models (GMMs), instead of DBSCAN, is also investigated in order to provide an assessment of the representative vibration frequencies, and of their estimation uncertainties, from the models determined with subspace identification.

This work compares the results obtained with the considered methods, applied on GNSS measurements collected with affordable receivers, on the case study of the Valdagena bridge (Section 2). The overall workflow is composed by a combination of several tools, as detailed in Section 3: (i) low-pass filtering in order to separate static/semi-static and dynamic components, (ii) random decrement as a powerful noise-reduction tool, (iii) system identification, and (iv) machine learning tools in order to support the automatic implementation of the modal analysis procedure. The obtained results are shown in Section 4, while some discussion and conclusions are drawn in Section 5.

2. Valgadena Bridge, Reference Leveling and Installed GNSS Stations

2.1. Reference Monitoring via Geometric Leveling

The Valgadena viaduct is located in the Veneto region (Italy), as seen in Figure 1a. It is a first-category bridge supported by two cement pillars, measuring 300 m in length and 160 m in height (Figure 1b). Please notice that, due to the non-orthogonal viewing angle of the satellite, the bridge direction appears slightly different in Figure 1a with respect to Figure 1b.

The Valgadena viaduct is well known to the authors, who have repeatedly been involved in its monitoring since its initial inspection, which is briefly described below. The initial inspection took place on 24 November 1990, utilizing classical geometric leveling techniques employing Wild/Leica NA2 levels with parallel plate and invar staffs spaced less than 30 m apart, thus constituting a high-precision geometric leveling (Figure 2). The inspection occurred within a single day, employing 20 vehicles with a total weight of 1000 tons, arranged in precise and progressive loading configurations to test the bridge behavior once subjected to loads exceeding those set during design. The settlements and corresponding returns to initial positions, on one hand, and local deformations in critical areas of the most stressed cross-sections of the bridge, as previously mentioned, were measured using high-precision leveling (0.1 mm on vertical displacement components compared to expected displacements on the order of centimeters) and through vibrating wire strain gauges (with a sensitivity of microns on strain gauge lengths of 500 mm), connected to two electronic data acquisition and processing units. For dynamic testing, a series of accelerometers were positioned on the bridge at predetermined points under the action of vehicular traffic, providing information on the modes and natural frequencies of vibration of the viaduct, as well as other structural parameters. The dynamic test was conducted to enable subsequent repetition, allowing for any variations in stiffness and damping parameters of the structure to be identified.

2.2. GNSS-Based Monitoring System

The Valgadena viaduct, following the just described initial inspection, has subsequently undergone multiple monitoring sessions, resulting in a well-documented record of its behavior. All monitoring activities, up to the most recent ones, have been conducted using classical, reliable but laborious, geometric leveling. Instead, the proposed methodology, which is based on a proper processing of GNSS measurements, is operationally and logistically much simpler to execute.

The movements of the structure, which is made of reinforced concrete, were assessed by using two low-cost multi-constellation (GPS—Global Positioning System, GLONASS—GLObal NAvigation Satellite System, Galileo, BeiDou, QZSS—Quasi-Zenith Satellite System) multi-band GNSS Emlid M2 receivers (each costing less than EUR 1000) as rovers, also named hereafter station 1 (St1) and 2 (St2), positioned along the viaduct bridge, as shown in Figure 1b, above the safety barriers in order to reduce the impact of the multi-path effect. St1 and St2 were positioned at approximately one-third and two-thirds of the bridge’s length, and they acquired data at 10 Hz. Such sampling frequency imposes a limit on the maximum frequency of the detectable modes with the proposed approach; nevertheless, the use of higher-rate GNSS receivers (e.g., 50–100 Hz) can overcome such limitation [6]. GNSS measurements were obtained by exploiting the corrections (in Real-Time Kinematic mode) provided by a base station, an Emlid RS2 receiver, located approximately 300 m from St1 and St2. The obtained solution takes into account of the observations provided by all the available constellations/satellites (≥14 satellites were available at each epoch). Details on the technical specifications of the GNSS receivers deployed in this work can be found on the manufacturer’s website [57].

The methodology described in the next section considers the projected coordinates determined from the data collected by St1 and St2 as inputs for the remaining of the procedure. To be more specific, the results shown in Section 4 will consider three cases: using data only from St1, only from St2, or from both of them (St1 + St2). Furthermore, it is worth to notice that, despite the fact that in this paper only the RTK case was considered, a similar procedure could be implemented with NRTK measurements as well.

3. Methodology

The proposed methodology is based on properly processing GNSS measurements, in order to significantly reduce measurement noise in the signal dynamic component and automatically determine the principal system modes. To this aim, the following workflow is proposed:

• High-frequency GNSS position measurements for St1 and St2 were obtained by exploiting corrections from a local base station very close to St1 and St2. The obtained positions were expressed in an ad hoc local (Cartesian) coordinate reference system (CRS), where the x and y axes correspond to the bridge longitudinal and transverse horizontal directions and z to the vertical one. The use of an ad hoc local reference system was previously suggested in [58,59,60] in order to ensure the best result accuracy. Alternatively, the positions can be projected in UTM (Universal Transverse Mercator) coordinates, leading however to a minor processing precision. Then, the coordinates are separately processed as described in the following steps.
• The static/semi-static component was computed by means of low-pass filtering. Hence, the dynamic component was obtained by subtracting the outcome of the low-pass filtering step from the original signal.
• Random decrement signatures were computed in order to reduce the noise impact and ease the identification of vibration modes.
• Automatic modal analysis was implemented by executing the following steps:
  • Subspace system identification of the signals collected by the GNSS receiver/s, varying the system order (i.e., the complexity of the system model), whose correct value is usually unknown.
  • Apply machine learning techniques, based on the use either of DBSCAN or GMM, in order to automatically determine the repetitive presence of certain poles in models identified at the previous step, when varying the system order, while discarding non-repetitive poles, considered as outliers, e.g., typically caused by noise. A vibration mode is associated to each cluster of poles and its statistical characteristics are determined based on the values of the poles within such cluster.

The workflow is summarized in Figure 3, and, more specifically, in Figure 4 for what concerns the last part of the procedure for determining the system modes. The next subsections describe each step of the procedure in detail.

Once a model is determined based on the current data, a comparison, for instance on the values of the estimated frequencies, is made with previously identified values of the considered structure. The rationale of OMA is that structural issues usually cause changes in the bridge vibration behavior, which should hence be visible when comparing vibration modes computed in different time instants: providing a reliable assessment of the modes and of the corresponding estimation uncertainty is fundamental to such aim; hence, comparisons on the reliability of the estimated uncertainties will be considered in the following as well.

It is worth to notice that the two options that will be considered in Section 3.4 in order to implement Automatic OMA, i.e., DBSCAN and GMM, are currently considered as alternatives: GMM is proposed in this paper as a potential substitute of DBSCAN, which has already been proposed in [54]. Furthermore, in addition to such two options, two different alternatives concerning the inputs for DBSCAN and GMM will be considered and compared as well.

3.1. Separation of the Dynamic Component

This subsection aims at presenting the implemented method in order to partially reduce signal noise and extract the dynamic component from the original signal, in order to properly perform vibration analysis. The semi-static component of a structure behavior is often extracted by low-pass filtering, with several kinds of filtering approaches which have been quite effectively implemented in the literature: sufficiently long moving average (MA) filters, Chebyshev and wavelet ones [61,62], just to mention certain of the most commonly adopted. Previous studies have already shown that multipath error affecting GNSS measurements is mainly distributed on relatively low frequencies, i.e., ≤0.2 Hz [63]. Hence, a low-pass (fourth order) Butterworth filter with a cut-off frequency of $0.2$ Hz was used in order to isolate the semi-static component. Such one-dimensional Butterworth filter is separately applied to each of the projected components of the GNSS measurements. Then, the dynamic component is obtained by subtracting the semi-static part from the original signal, as shown in Figure 3.

3.2. Random Decrement

The random decrement technique for vibration signature analysis dates back to the late 1960s [64]. It has already been widely used for modal analysis of structures and infrastructures [65], even in the OMA case [66]. Furthermore, despite the fact that this technique was originally thought for linear stationary systems, it has recently been extended in order to deal with non-stationary [67,68] and non-linear [69] ones too. A brief summary of the standard version for linear stationary systems is reported in the following.

Given a certain time history of a signal, its basic idea is that of averaging a certain number of signal segments properly extracted from the available time samples. Despite the fact that several criteria can be used to choose the segments to be selected, the simplest, and most commonly used one, is that of specifying a certain signal threshold in order to trigger the selection of a segment. Let $\mathtt{x}\left(t\right)$ be the considered 1D signal at time t, with known history for instance from time $t_0$ to $t_f$. Hence, given the discrete-time nature of our measurements, the available samples are assumed to be $\{\mathtt{x}\left(t_0\right),\mathtt{x}(t_0+T_s),\mathtt{x}(t_0+2T_s),\cdots ,\mathtt{x}\left(t_f\right)\}$, where $T_s$ is the sampling period. Let $\overline{\mathtt{x}}$ be the triggering threshold, and $n_{RD}$ the length of each extracted segment, e.g., the number of samples in the segment. Then, $\mathtt{x}\left(t_1\right)$ is selected as the starting point of an extracted segment (i.e., $(\mathtt{x}\left(t_1\right),\mathtt{x}(t_1+T_s),\mathtt{x}(t_1+2T_s),\cdots ,\mathtt{x}\left(t_{1+(n_{RD}-1))}\right)$) if the signal crosses the $\overline{\mathtt{x}}$ value there. It is also worth to notice that: (i) $\overline{\mathtt{x}}$ usually is expressed as a properly scaled version of the standard deviation of $\mathtt{x}$, (ii) the triggering criteria can involve a condition for the first derivative of $\mathtt{x}\left(t\right)$ as well; hence, the triggering condition can be seen as a common initial condition (in terms of the value of $\mathtt{x}$ and of its time derivative) for all the extracted segments. The general triggering condition, introduced in [70], is as follows:

$${\overline{\mathtt{x}}}_1\le \mathtt{x}\left(t\right)\le {\overline{\mathtt{x}}}_2,{\overline{\dot{\mathtt{x}}}}_1\le \dot{\mathtt{x}}\left(t\right)\le {\overline{\dot{\mathtt{x}}}}_2$$

(1)

where ${\overline{\mathtt{x}}}_1$, ${\overline{\mathtt{x}}}_2$, ${\overline{\dot{\mathtt{x}}}}_1$, ${\overline{\dot{\mathtt{x}}}}_2$ are properly selected thresholds, e.g., ref. [69] suggests to choose ${\overline{\mathtt{x}}}_i$ and ${\overline{\dot{\mathtt{x}}}}_i$, for $i=\{1,2\}$, within the following intervals: $(0.75{\sigma }_{\mathtt{x}},2.2{\sigma }_{\mathtt{x}})$ and $(0,2{\sigma }_{\dot{\mathtt{x}}})$, where ${\sigma }_{\mathtt{x}}$ and ${\sigma }_{\dot{\mathtt{x}}}$ represent the standard deviation of $\mathtt{x}$ and $\dot{\mathtt{x}}$, respectively.

Then, the random decrement $D_{\mathtt{x}}\left(\tau \right)$ of the signal $\mathtt{x}$ at time $\tau$ is computed as follows, where $n_{avg}$ is the number of extracted segments:

$$D_{\mathtt{x}}\left(\tau \right)={\displaystyle \frac{1}{n_{avg}}}\sum _j\mathtt{x}(t_j+\tau )$$

(2)

for all the ${\tau }_j$ such that $t_j$ is in the set ${\mathcal{T}}_{RD}$ of the starting instants of an extracted signal segment, i.e., $\mathtt{x}\left(t_j\right)$ (and/or its derivative) satisfies the triggering condition. When the triggering condition is just the equality of $\mathtt{x}$ with a certain prefixed value, the random decrement $D_{\mathtt{x}}(·)$ can be interpreted as the system free vibration response [66].

However, assuming the signal to be output of a linear, time-invariant system, excited by a zero-mean Gaussian random process, then it is possible to prove that the ensemble average (of segments extracted from the same time history) corresponds to averaging over all the possible signal realizations (under an ergodicity assumption), and the above equation corresponds to computing the following conditional expected value:

$${\displaystyle \frac{1}{n_{avg}}}\sum _j\mathtt{x}(t_j+\tau )\approx \mathbb{E}\left[\mathtt{x}(t+\tau )\right|t\in {\mathcal{T}}_{RD}]$$

(3)

where $\mathbb{E}[·]$ stands for the expectation operator. This corresponds to computing the expected value of $\mathtt{x}(t+\tau )$ given a proper “initial condition” for $\mathtt{x}\left(t\right)$ (and $\dot{\mathtt{x}}\left(t\right)$, if needed).

When the triggering condition depends only on a triggering level $\overline{\mathtt{x}}$ on the $\mathtt{x}\left(t\right)$ value, the random decrement $D_{\mathtt{x}}$ results to be proportional to the autocorrelation of the signal [64]:

$$D_{\mathtt{x}}\left(\tau \right)={\displaystyle \frac{R_{\mathtt{x}}\left(\tau \right)}{R_{\mathtt{x}}\left(o\right)}}\overline{\mathtt{x}}$$

(4)

where $R_{\mathtt{x}}\left(\tau \right)$ is the autocorrelation function of the random process $\mathtt{x}$, i.e., $R_{\mathtt{x}}\left(\tau \right)=\mathbb{E}\left[\mathtt{x}(t+\tau )\mathtt{x}\left(t\right)\right]$. The reader is referred to [69,71] for the relation between random decrement and the autocorrelation function in the most general case. Cross-random decrements, involving different signals, can be considered as well [72].

In this work, random decrements were computed on the dynamic component extracted from the projected coordinates, derived as described in the previous subsection. It is worth to notice that, differently from the previous subsection, random decrements have not been computed independently on each dimension of the dynamic component. To be more specific, in order to avoid losing the temporal correlation between the different dimensions of the dynamic component signal, the triggering condition was computed at each time instant on the norm of the dynamic component vector, and, consequently, the set ${\mathcal{T}}_{RD}$ is the same for all the dimensions of the considered random decrements, e.g., (2) is applied to the three dimensions of the dynamic components of the projected GNSS coordinates, but keeping for all of them the same set of starting instants ${\mathcal{T}}_{RD}$.

3.3. Subspace System Identification

This subsection describes the implementation of a time domain system identification method (first step of the procedure in Figure 4), as frequently performed in combination with the random decrement technique [66], in order to properly determine the modal characteristics of the time domain “response” computed in terms of random decrements. Nevertheless, differently from most of the previously implemented approaches, in the next subsection system identification will be combined with some automatic procedures (see also [53,54]) in order to enable automatic modal analysis. This subsection follows with a short introduction on subspace identification, and in particular on the implemented method, while the last part of the subsection is dedicated to describe the relation of subspace identification with the random decrements produced as outcomes of the procedure presented in the previous subsections. The reader is referred to [73] for a more detailed presentation of subspace identification methods.

For simplicity of notation, consider a discrete time representation of the system, where the system time period is conventionally 1, being trivial the adaptation of such system description to the general case, i.e., time period $=T_s$, as previously stated. Then, assume that the system behavior can be described (in innovation form [73,74,75,76]) as follows:

$$\begin{array}{ccc}\hfill x(t+1)& =& \hfill Ax\left(t\right)+Bu\left(t\right)+Ke\left(t\right)\hfill \end{array}$$

(5)

$$\begin{array}{ccc}\hfill y\left(t\right)& =& \hfill Cx\left(t\right)+Du\left(t\right)+e\left(t\right)\hfill \end{array}$$

(6)

where $x\left(t\right)$ is the system state at time t, $e\left(t\right)$ is assumed to be white noise, $u(·)$ is the system input, $y(·)$ the output, and A, B, C, D, K describe the system characteristics: A describes the system time dynamics, C represents the relation between state $x\left(t\right)$ and output y at any time t, in the most general case the input u can be related both to the state x by means of B and to directly to the output through D, whereas in the innovation form system representation the random noise $e\left(t\right)$ directly influences the output $y\left(t\right)$ while it is pre-multiplied by K in order to provide its contribution to $x(t+1)$. Furthermore, the size n of the state vector $x\left(t\right)$ is also called the system order. The reader is referred to [75,77] for a more in-depth description of stochastic linear system state models. System identification aims at determining the values of the parameters A, B, C, D, K given the values of system inputs and outputs in a certain time interval. Given the aims of this work, our major interest is on the eigenvalues of A, which are clearly strongly related with the system dynamic characteristics [75,77]. It is worth to notice that the state-space innovation form of a linear stochastic system, reported in (5) and (6), is strictly related to Kalman filtering theory. A detailed description of stochastic linear system characteristics and the derivation of such innovation from system representation is out of the scope of this work; nevertheless, to such aim we refer the readers to the sections dedicated to the description of the innovation representation in [73,75,78] and to [75,77,79,80] for a description of stochastic linear systems, other properties of Kalman filtering, and estimation theory.

Since in this work the random decrements, here thought as the system output, can be interpreted as the system free vibration response, inputs can be assumed to be zero hereafter, the system model can be assumed to be that in (7) and (8), and the identification problem reduces to estimating A, C, K. In accordance with such observation, this work implements the subspace system identification method described in the following. Nevertheless, other system identification approaches could be considered, in particular in order to deal with inputs, typically involving the oblique projection of the future observations onto the past (e.g., N4SID—Numerical algorithms for Subspace State Space System Identification, MOESP—Multivariable Output Error State Space, PBSID—Predictor Based Subspace IDentification [73,76]).

$$\begin{array}{ccc}\hfill x(t+1)& =& \hfill Ax\left(t\right)+Ke\left(t\right)\hfill \end{array}$$

(7)

$$\begin{array}{ccc}\hfill y\left(t\right)& =& \hfill Cx\left(t\right)+e\left(t\right)\hfill \end{array}$$

(8)

Using a notation similar to [79], let the past and future data matrices, ${\mathbf{Y}}_t^{-}$ and ${\mathbf{Y}}_t^{+}$, the state ${\mathbf{X}}_t$ and the noise ${\mathbf{E}}_t^{+}$ matrices be defined as follows:

$${\mathbf{Y}}_t^{-}=\left[\begin{array}{cccc}y\left(0\right)& y\left(1\right)& \cdots & y(N-1)\\ y\left(1\right)& y\left(2\right)& \cdots & y\left(N\right)\\ ⋮& ⋮& \cdots & ⋮\\ y(t-1)& y\left(t\right)& \cdots & y(t+N-2)\end{array}\right]$$

(9)

$${\mathbf{Y}}_t^{+}=\left[\begin{array}{cccc}y\left(t\right)& y(t+1)& \cdots & y(t+N-1)\\ y(t+1)& y(t+2)& \cdots & y(t+N)\\ ⋮& ⋮& \cdots & ⋮\\ y(2t-1)& y\left(2t\right)& \cdots & y(2t+N-2)\end{array}\right]$$

(10)

$${\mathbf{X}}_t=\left[\begin{array}{cccc}x\left(t\right)& x(t+1)& \cdots & x(t+N-1)\end{array}\right]$$

(11)

$${\mathbf{E}}_t^{+}=\left[\begin{array}{cccc}e\left(t\right)& e(t+1)& \cdots & e(t+N-1)\\ e(t+1)& e(t+2)& \cdots & e(t+N)\\ ⋮& ⋮& \cdots & ⋮\\ e(2t-1)& e\left(2t\right)& \cdots & e(2t+N-2)\end{array}\right]$$

(12)

Let $\Gamma$, the extended observability matrix, and $H_s$ be as follows:

$$\Gamma =\left[\begin{array}{c}C\\ CA\\ C{A}^2\\ ⋮\\ C{A}^{t-1}\end{array}\right]$$

(13)

$$H_s=\left[\begin{array}{ccccc}I& 0& \cdots & 0& 0\\ CK& I& \cdots & 0& 0\\ ⋮& ⋮& & \ddots \\ C{A}^{t-2}K& C{A}^{t-3}K& \cdots & CK& I\end{array}\right]$$

(14)

The observability matrix is related to the observability property of a linear system, i.e., it is related to the level at which the system state can be inferred from the outputs [75]. Instead, despite being a lower triangular matrix, $H_s$ is partially similar to a Hankel matrix [78], which is often used in stochastic realization (e.g., Ho-Kalman algorithm [81,82]) and subspace identification approaches [73,76,83].

Then, from (5) and (6), discarding the inputs:

$${\mathbf{Y}}_t^{+}=\Gamma {\mathbf{X}}_t+H_s{\mathbf{E}}_t^{+}$$

(15)

It is worth to notice that the state component $\Gamma {\mathbf{X}}_t$ in the above equation “belongs to the past”; hence, it can reasonably be computed as the part of ${\mathbf{Y}}_t^{+}$ predictable from the past.

Let ${\widehat{\Sigma }}_{+-}$, ${\widehat{\Sigma }}_{++}$, ${\widehat{\Sigma }}_{--}$ be

$$\begin{array}{ccc}\hfill {\widehat{\Sigma }}_{+-}& =& \hfill {\displaystyle \frac{1}{N}}{\mathbf{Y}}_t^{+}{{\mathbf{Y}}_t^{-}}^{\top }\hfill \end{array}$$

(16)

$$\begin{array}{ccc}\hfill {\widehat{\Sigma }}_{++}& =& \hfill {\displaystyle \frac{1}{N}}{\mathbf{Y}}_t^{+}{{\mathbf{Y}}_t^{+}}^{\top }\hfill \end{array}$$

(17)

$$\begin{array}{ccc}\hfill {\widehat{\Sigma }}_{--}& =& \hfill {\displaystyle \frac{1}{N}}{\mathbf{Y}}_t^{-}{{\mathbf{Y}}_t^{-}}^{\top }\hfill \end{array}$$

(18)

Then, compute the best prediction $\mathbb{E}\left[{\mathbf{Y}}_t^{+}\right|{\mathbf{Y}}_t^{-}]$ of ${\mathbf{Y}}_t^{+}$ given ${\mathbf{Y}}_t^{-}$, compute its Singular Value Decomposition (SVD) and take its best n-order approximation:

$$\Gamma {\mathbf{X}}_t\approx \mathbb{E}\left[{\mathbf{Y}}_t^{+}\right|{\mathbf{Y}}_t^{-}]={\widehat{\Sigma }}_{+-}{\widehat{\Sigma }}_{--}^{-1}{\mathbf{Y}}_t^{-}=US{V}^{\top }\approx U_n{S}_n{V}_n^{\top }$$

(19)

where $U_n$ and $V_n$ are the matrices obtained by considering only the first n columns from U and V, whereas $S_n$ is obtained by selecting the fist n rows and columns from S. Actually, as noticed in [73], the above approach can be generalized by factorizing $W_r\mathbb{E}\left[{\mathbf{Y}}_t^{+}\right|{\mathbf{Y}}_t^{-}]W_c$, where different choices for the matrices $W_r$ and $W_c$ lead to different implementations proposed in the literature.

From the equation above it immediately follows that $\Gamma {\mathbf{X}}_t\approx U_n{S}_n{V}_n^{\top }$; hence, $\Gamma$ and ${\mathbf{X}}_t$ can be set as follows:

$$\Gamma =U_n{S}_n^{1/2}$$

(20)

$${\mathbf{X}}_t=S_n^{1/2}{V}_n^{\top }$$

(21)

Finally, the A matrix, whose eigenvalues define the system modes, can be extracted from $\Gamma$ exploiting for instance the shift invariance method: let ${\Gamma }_u$ and ${\Gamma }_d$ be the portions of $\Gamma$ obtained by discarding the last and the first block, respectively. Then, it is clear from (13) that A can be obtained as $A={\Gamma }_u^{-L}{\Gamma }_d$, where ${\Gamma }_u^{-L}$ is the left inverse of ${\Gamma }_u$. The reader is referred to [73,79] for details on the computation of the other parameters of the system (7) and (8).

Subspace system identification is applied in this work to the random decrements computed as described in the previous subsection, i.e., the time series of the output signal $y(·)$ is composed by the random decrements on the different dimensions (e.g., two planimetric and the up direction), obtaining estimates of the system matrices, in order to obtain a good representation of the observations by means of the Equations (7) and (8). Since this work deals with modal analysis, focusing on determining reliable estimates of the values of the eigenvalues of A is of major importance: the next subsection will deal with this problem.

3.4. Clustering to Support Automatic OMA

This subsection describes the steps in the procedure in Figure 4 after the system identification one.

It is worth to notice that the system order n plays a fundamental role in the identified model (described by (7), (8)). The “correct” system order for an infrastructure is typically initially unknown; hence, the procedure described in the previous subsection is usually applied for different values of the system order n and, then, a suitable system order shall be selected according to a proper criterion. Despite the fact that several criteria have been proposed in the system identification literature in order to select a proper model, such criteria mostly deal with the description ability of the considered model with respect to the model complexity. Instead, as suggested in [53,54], in this work modal analysis takes into account of the most reliable vibration characteristics emerging from different identified models.

To such aim, first, the system identification procedure is applied to the random decrements for N different order values, from $n_{min}$ to $n_{max}$. Then, the list $\mathcal{L}$ of all the determined system poles $\left\{p_i\right\}$, computed as the eigenvalues of the identified A matrices, is analyzed. The rationale of considering such list of poles, is that poles associated to real oscillatory modes shall be present in a significant portion of the identified models. Differently, poles caused by noise shall be present in only a few of the computed models.

Similarly to [53,54], Hard Validation Criteria (HVC) are applied to such pole list, selecting poles according to the following:

• Only oscillatory modes, hence associated to vibrations in a free decay response, are considered, i.e., those associated to poles complex conjugated.
• Only poles associated to positive dampings can be associated to mechanical systems and infrastructures. Furthermore, infrastructures should realistically have damping ratios lower than a certain threshold, e.g., 20% in [53].

Then, in accordance with the rationale mentioned above, [54] proposed the use of DBSCAN, applied to the coordinates of poles in $\mathcal{L}$, in order to separate them in different clusters. Then, from each detected cluster is derived just one “confirmed” pole (which can be determined by, for instance, averaging those in its cluster, as will be described more precisely in the following), considered as a real one at the end of the procedure. In addition to such a possibility, this work considers other alternative options, leading to the following four cases, where the first, case 1A, corresponds to the already mentioned one:

• Case 1A: DBSCAN applied to the coordinates of poles in $\mathcal{L}$.
• Case 1B: DBSCAN applied to frequencies and damping ratios associated to the poles in $\mathcal{L}$.
• Case 2A: GMM applied to the coordinates of poles in $\mathcal{L}$.
• Case 2B: GMM applied to frequencies and damping ratios associated to the poles in $\mathcal{L}$.

Such four cases are considered in this work as alternatives, whose performance will be compared in Section 4. The four cases will be described in the following, along with a brief description of the characteristics of DBSCAN and GMM, the procedure to exploit their results in order to automate the assessment of the most representative system vibration modes, and the proposed methods to assess the uncertainties on the estimated modes.

The clustering performance of the DBSCAN algorithm relies on properly setting the values of two design parameters: $ϵ$, the radius of a point neighborhood to be considered in order to determine whether such point is an outlier (to be discarded) or part of a cluster, and $minPts$, which is the minimum amount of pole instances needed in order to support the existence of the corresponding cluster. Hence, interestingly, DBSCAN is able to quite effectively perform outlier rejection. Nevertheless, properly setting the values of $ϵ$ and $minPts$ may be difficult for the following reasons:

• Since a real mode should have been identified in a significant portion of the models, the value of $minPts$ shall be set depending on N, the number of considered models, e.g., ≈$N/3$.
• A possible choice for selecting a meaningful value for $ϵ$, see for instance [53,84], is that of using the elbow rule on the k-distance graph of the pole coordinates.

Let us consider case 1A: once clusters have been determined, each of them should be associated to a mode and hence the corresponding estimates for frequency f and damping ratio $\xi$ should be provided for each of them. To be more precise, each determined cluster ${\mathcal{L}}_j$ identifies a mode, whose characteristics, in terms of f and $\xi$, are computed from those of the poles $\left\{p_i\right\}$ in ${\mathcal{L}}_j$: let $f_i$ and ${\xi }_i$ be the frequency and damping ratio values computed for the pole $p_i$ in ${\mathcal{L}}_j$; then, the estimated frequency $f_j$ and damping ratio ${\xi }_j$ associated to such cluster are determined as the averages of the frequencies and damping ratios of all the $p_i$ in ${\mathcal{L}}_j$, i.e., $f_j={\sum }_{p_i\in {\mathcal{L}}_j}{\displaystyle \frac{f_i}{N_j}}$, ${\xi }_j={\sum }_{p_i\in {\mathcal{L}}_j}{\displaystyle \frac{{\xi }_i}{N_j}}$, where $N_j$ is the number of poles in ${\mathcal{L}}_j$. Then, uncertainties ${\widehat{\sigma }}_{f_j}$ and ${\widehat{\sigma }}_{{\xi }_j}$ on $f_j$ and ${\xi }_j$ are assessed by computing the standard deviation from the same samples.

Since DBSCAN results are influenced by the choice of $ϵ$ and $minPts$, the impact of such parameters on the estimated pole parameters is assessed by repeating DBSCAN clustering for different values of $ϵ$ and $minPts$, including in the considered ranges of all their reasonable values, determined taking into account of a visual inspection of the diagram of poles and to the number of identified systems N. Then, the average value and variability of the estimates for f and $\xi$ are assessed considering the results of DBSCAN clustering over all the considered cases. Let $f_j$ be the frequency estimate provided for the j-th cluster and let $f_{j,k}$ be the estimate of such frequency coming from the corresponding cluster on the k-th run of DBSCAN. Then, $f_j$ is obtained as the average of the $f_{j,k}$ for all the values of k corresponding to DBSCAN runs wherein such cluster has been detected (it may not have been detected in certain runs), where the detection criterion for a certain cluster is the closeness to the frequency of the corresponding mode up to the considered DBSCAN run (when a cluster detected in run k is not compatible with any of the previously detected clusters it is considered as a new one, whereas if is compatible with more than a previously considered cluster it is assigned to the closest one). The damping ratio estimate ${\xi }_j$ can be obtained similarly.

Then, variability on the determined values $f_j$ and ${\xi }_j$ is assessed in the following two alternative ways:

(i)

${\sigma }_{f_j}$ is obtained as the standard deviation of the estimates $f_{j,k}$ used to compute $f_j$:

$${\sigma }_{f_j}=\sqrt{{\displaystyle \frac{{\sum }_{k=1}^{n_j}{(f_{j,k}-f_j)}^2}{n_j-1}}}$$

(22)

where $n_j$ is the number of DBSCAN runs where a $f_{j,k}$ contributing to $f_j$ has been found. ${\sigma }_{{\xi }_j}$ is computed similarly.

(ii)

Standard deviations ${\sigma }_f^{\prime }$ and ${\sigma }_{\xi }^{\prime }$ are computed by directly taking into account of the frequencies and damping ratios of the poles in all the clusters associated to $f_j$ in the different DBSCAN runs. It is worth to notice that the same pole may be counted several times, depending on the DBSCAN iterations where it has been considered in the cluster associated to $f_j$.

In case 1B, clustering is considered on a different set of parameter values related to the poles in $\mathcal{L}$: DBSCAN is applied to the set of frequency and damping ratio parameters related to the poles in $\mathcal{L}$. Since the range of frequencies and damping ratios may be quite different, centering and normalization of such parameter values are performed before clustering. Centering is obtained just by subtracting the respective mean values, whereas normalization is obtained by dividing by the maximum absolute variation.

Similarly to case 1A, estimates are obtained for manually optimized values of the DBSCAN parameters and according to the investigation among different runs of DBSCAN, as explained above.

Since the selection of such parameter values, and in particular of $ϵ$, can have repercussions on the DBSCAN clustering results, in this work a different approach is proposed, aiming at reducing the dependence of the results on the parameter selection. To be more specific, this work investigates the use of GMMs for clustering the N identified poles (cases 2A and 2B).

In case 2A, GMM aims at representing the distribution of the pole coordinates as a mixture of Gaussian distributions:

$$p\left(\mathbf{x}\right)=\sum _{i=1}^{n_G}{w}_i\mathcal{N}(\mathbf{x};{\mu }_i,{\Sigma }_i)$$

(23)

where $p\left(\mathbf{x}\right)$ represents the probability density of the pole position $\mathbf{x}$ according to the Gaussian mixture, given by the sum of $n_G$ Gaussians, with weights $\left\{w_i\right\}$ and characterized by the means $\left\{{\mu }_i\right\}$ and covariances $\left\{{\Sigma }_i\right\}$, where the values of the mixture parameters are usually assessed by means of the Expectation Maximization (EM) algorithm [85].

Since the ideal number of Gaussians is unknown, a GMM is learned for several values of $n_G$. Then, the Akaike Information Criterion (AIC, [86]) is used for model selection among the computed GMM models.

Once a GMM has been selected among the computed ones, the Gaussian centers can be used as the identified poles and the covariances to determine the corresponding estimation uncertainties. Furthermore, the covariance matrix also enables the determination of the poles in $\mathcal{L}$ which are present in any of the clusters of interest. To be more specific, the following criterion is used to determine whether the pole with frequency $f_i$ and damping ${\xi }_i$ is within cluster k, the latter associated to the already detected mode with frequency $f_k$ and damping ratio ${\xi }_k$.

$$\left[\begin{array}{cc}{f}_i-f_k& {\xi }_i-{\xi }_k\end{array}\right]{\Sigma }_k^{-1}\left[\begin{array}{c}{f}_i-f_k\\ {\xi }_i-{\xi }_k\end{array}\right]\le thr$$

(24)

where, in our tests, $thr=3$. If a pole satisfies the above condition for more than a detected mode, then it is assigned to the one minimizing the value on the left of the above inequality. It is worth to notice that, differently from DBSCAN, GMM does not have an explicit mechanism to implement outlier rejection. To such aim, it has proven useful in our tests to discard those Gaussians in the GMM representation (23) characterized by a weight $w_i\le \overline{w}$, where the threshold value $\overline{w}$ was set to 0.05 in our tests.

Similarly to case 1B, in case 2B GMM is applied directly to the frequencies and damping ratios of the poles. It is worth to notice that the covariance matrices of the Gaussians in GMM make it possible to differentiate along different directions; hence, differently from case 1B, normalization of the variables, before applying GMM, is not considered here. With the exception of the variables on which GMM is applied, the rest of the procedure of case 2B is analogous to 2A.

Finally, it is worth to notice that, since in cases 2A and 2B the order $n_G$ selection has been obtained by means of the AIC criterion, the investigation previously conducted for DBSCAN results, varying the DBSCAN input values $ϵ$ and $minPts$, is not replicated for cases 2A and 2B.

4. Results

This section presents the results obtained on data collected over 4 h at the Valdegana bridge on 29 May 2024, in low-wind conditions. First, the semi-static component of the local CRS coordinates of St1 and St2 is extracted by applying, separately on each of the three directions, a low-pass Butterworth filter, as shown in Figure 5a for the up direction for St2, just as an example (similar figures can clearly be obtained for the other directions and for St1). Then, the dynamic component is obtained by subtracting the extracted semi-static component from the projected GNSS coordinates (see Figure 5b for checking the results on the up direction for St2). Figure 5 shows the results on the up direction being this the one of major interest in this case study, as the displacements of the dynamic component along the two planimetric coordinates are significantly smaller than those along the up direction.

The random decrement technique, described in Section 3.2, is applied to the dynamic components, leading to the results shown in Figure 6 for the up direction for St2. Similar results were obtained for St1 as well. Figure 6a shows the response decay over a 4000 s time interval, whereas certain vibration modes are more clearly distinguishable in Figure 6b, on a 200 s interval. In addition to the up direction, Figure 6b also shows the random decrements along the longitudinal and transverse directions of the bridge. Vibrations along the up direction are apparently larger than those along the other directions: in terms of signal energy, computed on the 4000 s interval, longitudinal and transverse direction ones are one order of magnitude smaller with respect to the up one. Since the maximum value of the random decrements is strictly related to the threshold used as triggering condition, in order to show results more independent with respect to such value, the vertical axes in Figure 6a,b have been normalized with respect to their maximum value.

Then, subspace system identification is used in order to properly assess the parameters of the dynamic model described by (7) and (8). In our tests, model order is varied from 2 to 30 (N = 15). Then, poles are computed, as eigenvalues of the A matrices of the identified models, and the unfeasible ones are filtered out taking into account of the HVC criteria (see Figure 7a, for the poles obtained using data from both St1 and St2). The final selection of the system poles is performed either taking into account the clusters computed with DBSCAN or the GMM results, as shown in the following, focusing first on two alternatives for DBSCAN and then on GMM.

The list $\mathcal{L}$ of pole coordinates was determined based on poles detected for both St1 and St2. Then, Figure 7b shows DBSCAN results obtained setting $ϵ$ according to the elbow rule (see Section 3.4) and $minPts$ to 5, which experimentally proved to be a reasonable choice in this case study, making it possible to determine the most representative clusters, highlighted also for most of the other choices of the $minPts$ value, but having also a good outlier rejection effect. As previously explained in Section 3.4, once clusters have been determined, estimates for frequency f and damping ratio $\xi$ associated to a cluster are obtained in terms of the average values of those of each pole instance in the cluster. Focusing, for simplicity, on the numerical results for what concerns the detected vibration mode with lowest frequency $f_{min}$,

Table 1. Estimates of $f_{min}$ obtained in case 1A, using data from either both the stations or just one of them, applying DBSCAN $minPts=6$ and $ϵ$ selected according to the elbow rule.

	$\widehat{\mathit{f}}$ [Hz]	${\widehat{\mathit{\sigma }}}_{\mathit{f}}$ [Hz]	$\widehat{\mathit{\xi }}$	${\widehat{\mathit{\sigma }}}_{\mathit{\xi }}$
St1 + St2	0.133	0.009	0.14	0.07
St1	0.132	0.002	0.18	0.02
St2	0.134	0.011	0.13	0.08

reports the obtained results, distinguishing between those obtained while using only measurements from one of the two GNSS stations and from both of them. As a side effect of clustering on the pole, real and imaginary coordinates, in certain cases damping values for pole instances forming a cluster, may vary quite significantly. This issue can clearly be reduced by reducing the $ϵ$ value, which however typically leads to a larger number of identified clusters.

Since DBSCAN results are influenced by the choice of $ϵ$ and $minPts$, as mentioned in Section 3.4, to assess the impact of such parameters on the estimated system pole parameters,

Table 2. Estimates of $f_{min}$ obtained in case 1A: analysis of the results on different runs of DBSCAN varying its input parameter values.

	$\widehat{\mathit{f}}$ [Hz]	${\mathit{\sigma }}_{\mathit{f}}$ [Hz]	${\mathit{\sigma }}_{\mathit{f}}^{\prime }$ [Hz]	$\widehat{\mathit{\xi }}$	${\mathit{\sigma }}_{\mathit{\xi }}$	${\mathit{\sigma }}_{\mathit{\xi }}^{\prime }$
St1 + St2	0.133	0.002	0.005	0.17	0.01	0.01
St1	0.132	0.001	0.001	0.17	0.01	0.01
St2	0.135	0.001	0.007	0.14	0.02	0.07

compares the results obtained for the estimate of $f_{min}$ while varying $ϵ$ and $minPts$ in the ranges of $[0.005,0.1]$ and $[3,10]$, respectively. The $minPts$ value is typically suggested to be similar to $N/3$, as previously stated in Section 3.4. Hence, the range of $minPts$ values to be investigated was chosen to be reasonably centered with respect to $N/3=5$ and $N=6$, which experimentally proved to give good results in this case study. Instead, the interval of values for $ϵ$ were chosen after a visual inspection of the diagram of poles in order to include all the range of reasonably possible values for such parameter in the considered case study. Average value and variability of the estimated values for f and $\xi$ are evaluated taking into account the estimates out-coming from DBSCAN clustering over all the considered cases, as explained in Section 3.4. Table 2 compares the results of such analysis distinguishing between those obtained while using only measurements from one of the two GNSS stations and from both of them. The results reported in Table 2 confirm those of Figure 7: variability on the estimated damping ratios is proportionally higher than that on f. Furthermore, comparing the estimated values and their corresponding variability, uncertainty (due to the unknown correct input parameter values) assessed with ${\sigma }_f$ and ${\sigma }_{\xi }$ appears to be more realistic than with ${{\sigma }^{\prime }}_f$ and ${{\sigma }^{\prime }}_{\xi }$.

Furthermore, Figure 8 shows the variability in the number of clusters with respect to the values of $minPts$ (a) and $ϵ$ (b) in case 1A, considering for simplicity only the station 1 case, whereas the variability in ${\sigma }_f$ and ${\sigma }_{\xi }$ with respect to $ϵ$ is visible in Figure 9. Since the variability in $minPts$ has a quite minor impact with respect to $ϵ$ on the values of ${\sigma }_f$ and ${\sigma }_{\xi }$, the corresponding figures are not shown here.

Then, considering case 1B, DBSCAN was applied to the pole frequency and damping values instead, obtaining the results shown in Figure 10, for data from both St1 and St2, selecting $ϵ$ according to the elbow rule and $minPts=6$. To be more precise, the point coordinates used as inputs for DBSCAN were obtained by centering and normalizing frequency and damping pole values. Normalization was obtained by dividing by the maximum absolute variation. Figure 10 shows that damping ratio and frequency variations within an identified cluster appear to be very limited in this case.

Table 3. Estimates of $f_{min}$ obtained in case 1B, using data from either both the stations or just one of them, applying DBSCAN $minPts=6$ and $ϵ$ selected according to the elbow rule.

	$\widehat{\mathit{f}}$ [Hz]	${\widehat{\mathit{\sigma }}}_{\mathit{f}}$ [Hz]	$\widehat{\mathit{\xi }}$	${\widehat{\mathit{\sigma }}}_{\mathit{\xi }}$
St1 + St2	0.130	0.013	0.175	0.011
St1	0.132	0.001	0.167	0.009
St2	0.128	0.019	0.182	0.007

shows the results for what concerns the estimate of $f_{min}$.

Then,

Table 4. Estimates of $f_{min}$ obtained in case 1B: analysis of the results on different runs of DBSCAN varying its input parameter values.

	$\widehat{\mathit{f}}$ [Hz]	${\mathit{\sigma }}_{\mathit{f}}$ [Hz]	${\mathit{\sigma }}_{\mathit{f}}^{\prime }$ [Hz]	$\widehat{\mathit{\xi }}$	${\mathit{\sigma }}_{\mathit{\xi }}$	${\mathit{\sigma }}_{\mathit{\xi }}^{\prime }$
St1 + St2	0.133	0.002	0.005	0.174	0.021	0.007
St1	0.133	0.001	0.001	0.164	0.003	0.007
St2	0.130	0.002	0.012	0.176	0.003	0.006

compares the variability in the results obtained, by applying DBSCAN to the pole frequency and damping values, for the estimate of $f_{min}$ while varying $ϵ$ and $minPts$ in the ranges of $[0.005,0.1]$ and $[3,10]$, respectively. In comparison with Table 2, the results in Table 4 show a lower difference in the estimates provided in the different cases and a lower estimated variability as well.

Figure 11 shows the variability in the number of clusters with respect to the values of $minPts$ (a) and $ϵ$ (b) in case 1B, considering for simplicity only the station 1 case, whereas the variability in ${\sigma }_f$ and ${\sigma }_{\xi }$ with respect to $ϵ$ is visible in Figure 12. Since the variability of $minPts$ has a quite minor impact with respect to $ϵ$ on the values of ${\sigma }_f$ and ${\sigma }_{\xi }$, the corresponding figures are not shown here.

For what concerns cases 2A ad 2B, several instances of GMM are run, varying the system order up to $n_G=10$, and selecting the best model according to the AIC method. The cluster detection performance obtained by using GMM, as described in Section 3.4, can be visually evaluated in Figure 13 and Figure 14 (for data from both St1 and St2), for case 2A and 2B, respectively. It is worth to notice that, in the latter case, differently from the procedure implemented for DBSCAN, frequencies and damping ratios have not been centered and normalized. Clusters were considered only if formed by at least $minPts=7$ of the poles in $\mathcal{L}$.

Focusing on $f_{min}$,

Table 5. Estimates of $f_{min}$ obtained with GMM: cases 2A and 2B, and varying the used stations.

Case	$\widehat{\mathit{f}}$ [Hz]	${\mathit{\sigma }}_{\mathit{f}}$ [Hz]	$\widehat{\mathit{\xi }}$	${\mathit{\sigma }}_{\mathit{\xi }}$
2A: St1 + St2	0.135	0.003	0.14	0.04
2A: St1	0.136	0.002	0.14	0.03
2A: St2	0.134	0.003	0.13	0.04
2B: St1 + St2	0.134	0.004	0.17	0.03
2B: St1	0.132	0.002	0.17	0.02
2B: St2	0.138	0.022	0.19	0.03

shows the estimates and the corresponding uncertainties obtained in case 2A and 2B, when considering any of the two GNSS stations or both of them. It is worth to notice that, in the first case, the uncertainty was evaluated as the sample standard deviation of the frequencies (or damping ratios) of the poles in $\mathcal{L}$ which are in a cluster (similarly to ${\sigma }_f^{\prime }$ and ${\sigma }_{\xi }^{\prime }$ for DBSCAN). Instead, in case 2B, uncertainties were directly derived from the covariances $\left\{{\Sigma }_i\right\}$ of the estimated Gaussians (see (23)), as the square roots of the corresponding diagonal elements. It is worth to notice that the results in terms of estimates and uncertainties in both case 2A and 2B appear to be quite similar, those of 2B being slightly more similar to those obtained with DBSCAN, and the use of two stations instead of just one is clearly more reliable.

Finally, Figure 15 and Figure 16 show the variability in ${\sigma }_f$ (a) and ${\sigma }_{\xi }$ (b) with respect to $n_G$, in case 2A and 2B, respectively. For simplicity, these figures are shown only for the case involving the use of both the stations.

To conclude,

Table 6. Estimates of the three modes with largest energy obtained with DBSCAN and GMM from measurements of both St1 and St2.

Case	$\widehat{\mathit{f}}$ [Hz]	$\widehat{\mathit{\xi }}$	$\widehat{\mathit{f}}$ [Hz]	$\widehat{\mathit{\xi }}$	$\widehat{\mathit{f}}$ [Hz]	$\widehat{\mathit{\xi }}$
1A	0.133	0.17	0.92	0.005	1.897	0.0018
1B	0.133	0.174	0.92	0.005	1.897	0.0018
2A	0.135	0.14	0.93	0.008	1.900	0.0011
2B	0.134	0.17	0.92	0.004	1.901	0.0009

summarizes the characteristics of the three estimated modes with higher energies, comparing the results obtained in the four considered cases. The frequencies of the estimated modes are apparently quite similar, whereas the variability on the estimated damping ratios is higher, as visible in certain of the previous results as well.

5. Discussion and Conclusions

Despite the fact that the combination of random decrements and time-domain-based identification methods have already been considered in previous works, in this paper its application to GNSS measurements acquired with affordable receivers, with corrections provided with a very short baseline, confirmed the possibility of obtaining quite reasonable results for what concerns determining the main vibration frequencies of the considered bridge; for instance, in all the considered cases, the mode with minimum frequency $f_{min}$ was detected, despite the fact that the estimation uncertainty was quite different depending on the considered case.

The acquisition rate of the GNSS receiver clearly introduces a limit on the maximum frequency of the detectable modes. However, on the one hand the impact of such limitation can be reduced by using higher-rate GNSS receivers [6], and, on the other hand, since a significant portion of the first bridge’s natural frequencies usually are ≤5 Hz [63], then even the use of the considered affordable GNSS receivers should not represent a significant limitation from this point of view.

A detailed analysis on the application of clustering in order to support automatic operational modal analysis showed that the results of the previously proposed use of DBSCAN [54], applied on the detected pole coordinates, are influenced by the selection of the $ϵ$ and $minPts$ parameter values used by DBSCAN. In particular, Figure 8 and Figure 11 show a significant impact of $minPts$ on the numbers of detected clusters, whereas $ϵ$ had a more significant impact on the variability in the estimated parameter values (frequencies and damping ratios), as expected. In particular, variations on the estimated damping ratio values may be quite significant among the poles within some of the identified clusters, as can be noticed for instance by visual inspection of Figure 7b and Figure 9b, or by checking the numerical results for what concerns the uncertainty ${\widehat{\sigma }}_{\xi }$ in Table 1.

In order to reduce the latter issue, in case 1B, DBSCAN was directly applied to the frequencies and damping ratio values of the poles in $\mathcal{L}$. Figure 10 shows a reduction, with respect to Figure 7b, in the variability in the damping and frequency values for what concerns the poles within a cluster. Such reduction is confirmed by the overall decrease in the uncertainty values ${\sigma }_f$ and ${\sigma }_{\xi }$ in Table 4 with respect to Table 2 (and in Figure 12b), showing a probable convenience in applying DBSCAN to the pole frequencies and damping ratio values, as in case 1B, which will be better investigated in our future works in different case studies.

In both case 1A and 1B, but more remarkably in 1A, the assessed parameter value variations determined as ${\sigma }_f$ and ${\sigma }_{\xi }$ appeared to be more realistic than those estimated as ${{\sigma }^{\prime }}_f$ and ${{\sigma }^{\prime }}_{\xi }$ (see Table 2 and Table 4, where the uncertainties should be such that estimates obtained with any of the two stations or by both of the two of them should be compatible), hence proposing ${\sigma }_f$ and ${\sigma }_{\xi }$ as more reliable estimates of such parameter uncertainties. Despite the fact that this work currently only focuses on the extraction of vibration mode characteristics from GNSS data, the comparison of such characteristics is of fundamental importance in OMA applications; hence, also properly assessing the uncertainties in the estimated parameters is remarkably relevant.

Furthermore, the use of GMM has been investigated in order to reduce the dependency of the clustering results on the choice of the $ϵ$ and $minPts$ parameter values. In particular, only $\overline{w}$ (which plays a role similar to $minPts$) was used in this case in order to perform cluster outlier rejection, i.e., to discard unreliable clusters.

Similarly to DBSCAN, GMM was applied both on the pole coordinate domain (case 2A) and directly on pole frequencies and damping ratios (case 2B). Differently from case 1B, since the values in the covariances of the Gaussians in GMM are automatically adapted to the values of the variables input to the GMM estimator, in case 2B, GMM was applied to frequencies and damping ratios without centering and normalizing them. Furthermore, in accordance with the results reported in Table 5, and with Figure 15 and Figure 16, the direct application of GMM to frequencies and damping ratios (case 2B) seems to have a minor impact with respect to the 1B case. Indeed, the uncertainties in the estimated values appear to be quite similar in the two cases in Table 5. The motivation for such slight different behavior with respect to case 1A–1B is probably the already mentioned greater flexibility of GMM in terms of the values of its covariances along the two different directions. It is also worth to notice that, thanks to the estimation of such covariances, the estimated parameter uncertainties can be conveniently determined directly from them in case 2B. Overall, comparing 2A and 2B, case 2B provides results more coherent with those of 1A and 1B, as shown in Table 6. Furthermore, Table 6 confirms that, in all the considered cases, the estimated modes with highest energies have comparable frequencies, and instead some more variability is shown in the detected damping ratios, in particular for smaller values. This aspect will be the subject of our future investigations.

To conclude, the investigation presented in this paper confirmed that clustering could be effectively used to determine vibration modes from the poles coming from a sequence of identified models, as previously proposed in [54]. However, this work investigated more in depth the impact of the input parameter values of DBSCAN on the obtained results, and two alternative options in order to assess the uncertainties of the estimated parameters. Furthermore, the use of GMM was proposed to address the impact of input parameters of the clustering algorithm on the obtained results.

Overall, the obtained results (in particular in terms of estimated value of frequency and damping ratio of $f_{min}$, the system vibration mode with lowest frequency) were quite similar in all the considered cases. In our tests, the use of DBSCAN clustering directly on the frequency and damping ratio values was more convenient (case 1B) than applying it to the original pole coordinates, whereas such convenience was only partially confirmed when using GMM. The use of GMM instead of DBSCAN reduced the dependence of the clustering results from the input parameters of the clustering procedure, and in particular from $ϵ$, whose value showed to have an impact on the DBSCAN results, e.g., on the uncertainties on the estimated parameters. While stating the overall superiority of GMM from just a case study is unreliable, given the above-mentioned convenience and the similar obtained results, GMM proved to be a valid alternative to DBSCAN in this case study, and, potentially, more in general in AOMA applications. Overall, the provided analysis indicated both certain critical aspects and some advantages in using clustering methods in AOMA applications, including some considerations on the evaluation of the uncertainty level of the estimated parameters.

Future investigations will involve the extension of the use of the developed procedure on other case studies in order to better highlight the pros and cons of the different considered options, and to confirm or disprove the conclusions derived by the results obtained in this work. Furthermore, the integration of the proposed methodology within a fully operative monitoring application, hence including the analysis and comparison of the system characteristics with data acquired at different epochs, will be considered as well.