A Damage Detection Approach for Axially Loaded Beam-like Structures Based on Gaussian Mixture Model

Axially loaded beam-like structures represent a challenging case study for unsupervised learning vibration-based damage detection. Under real environmental and operational conditions, changes in axial load cause changes in the characteristics of the dynamic response that are significantly greater than those due to damage at an early stage. In previous works, the authors proposed the adoption of a multivariate damage feature composed of eigenfrequencies of multiple vibration modes. Successful results were obtained by framing the problem of damage detection as that of unsupervised outlier detection, adopting the well-known Mahalanobis squared distance (MSD) to define an effective damage index. Starting from these promising results, a novel approach based on unsupervised learning data clustering is proposed in this work, which increases the sensitivity to damage and significantly reduces the uncertainty associated with the results, allowing for earlier damage detection. The novel approach, which is based on Gaussian mixture model, is compared with the benchmark one based on the MSD, under the effects of an uncontrolled environment and, most importantly, in the presence of real damage due to corrosion.


Introduction
Structural Health Monitoring (SHM) is certainly one of the most-discussed topics in the literature of mechanical, civil and aerospace engineering, due to the central relevance of safety, maintenance and quality. By adopting the most general definition, damage is the loss of a structure to perform its intended functions [1,2] due to unwanted changes of the geometric and material properties, boundary conditions and system connectivity with respect to an initial reference condition. An early detection of these changes plays a crucial role to carry out prompt maintenance actions [3]. For this reason, structures are monitored over time by adopting sensors to measure physical quantities that are related to the structural properties [4].
According to a data-driven approach, data coming from sensors are used to extract damage-sensitive features that can be statistically analysed to determine the current state of the system health. In many real applications, damage conditions and locations greatly vary from structure to structure and databases containing monitoring data related to all the possible damage conditions are not available. As a result, damage detection must be carried out by adopting an unsupervised learning approach (i.e., an approach that does not require damage-related data). First, monitoring data are acquired during a reference period, to statistically characterize the system response under healthy condition; damage is then assessed when a meaningful deviation from what is considered the normal behaviour is observed [2]. In real applications, unsupervised learning approaches are limited by the effects of environmental and operational variations that cause changes in the structural behaviour that, in turn, reflect in changes of the damage feature. Often, these changes are greater than those caused by damage at an early stage and they must be properly accounted for [5].
In this context, this paper is devoted to the development of unsupervised learning damage-detection approaches for beam-like structures. Beams are widely adopted in engineering structures, often representing fundamental elements of complex systems. Examples can be tie-rods of arches and vaults, diagonal braces of truss girders, struts and ties of space truss structures. More specifically, the main attention is here focused on axially loaded beams, i.e., tensioned beams. A vibration-based damage detection approach is used in this work, particularly suitable for slender structures which undergo significant vibration levels under operational conditions. Vibration-based damage detection relies on a basic assumption: the presence of damage alters the structural properties of a system (mass, stiffness, damping), which in turn, changes the dynamic response of the system [6]. Axially loaded beam-like structures can serve as a challenging case study: indeed, during their normal operation, these structures undergo dramatic changes of the axial load that reflect into significant changes of the dynamic response, thus making vibration-based damage detection a hard task.
A specific case study will be presented, which is the case of tie-rods (nevertheless, the proposed damage detection strategy can be adopted on any beam-like structure subject to axial load). Tie-rods are metallic beam elements adopted to balance the lateral force transmitted to the base of arches and vaults, both in modern and historical civil structures. The integrity of these elements is important for the overall structural equilibrium, since a tie-rod failure may result in a global failure. To fulfil their function, tie-rods are subject to a tension at the time of the installation. During their life, temperature variations cause changes in the thermal coefficients of both the tie-rod and the structure, which reflects into changes of the axial load. Other tension variations are due to deformation and displacement of the connecting walls, that may be caused by terrain crawl, subsidence of foundations or seismic events [7].
Many authors focused their attention on the axial load estimation, considering only two possible failure types: a low tension, that points out a loosen tie-rod, not exerting its intended function and therefore requiring replacement, or a high tension, pointing at a risk of overloading, due to an abnormal structural deformation. In this latter case, the tie-rod axial load can provide information about changes in the overall structural health.
When operating tie-rods are considered, there is a high uncertainty related to the real axial-load, which can be directly measured only by strain gauges or load cells, installed and calibrated prior to the initial tensioning procedure. This, of course, is not a viable solution for already installed tie-rods. For this reason, many researchers focused on indirect approaches for the identification of the tensile force [8], based on static [9,10], static-dynamic [11] or fully dynamic [12,13] in situ indirect measurements. While the main shortcomings of static and static-dynamic approaches are related to the complicated experimental procedures, fully dynamic methods are based on the estimation of the tie-rod modal parameters, which can be more easily estimated due to the slenderness of these structural elements that allows, for instance, the adoption of output-only operational modal analysis [14]. However, the problem is complex because modal parameters are not only related to the axial load but also to a number of other physical variables which are affected by a high uncertainty, such as geometrical and material properties and constraint characteristics. The inverse problem is non-linear and potentially ill-conditioned [15]; thus, many approaches that can be found in the literature are based on analytical solutions [16,17], minimization between the theoretical and the measured modal parameters [18,19], model updating [20] and genetic algorithm optimization [21].
From a review of the state of the art, all the above-mentioned references consider the excessive load or deformation as the only critical problem, neglecting any possible presence of damage in the tie-rod. On the other hand, experience shows that cracks due to ageing or corrosion phenomena can lead to high-risk situations. Often, these types of damage are hard to be detected through visual inspections, which are even made more complex by the difficulties to access tie-rods and carry out in situ tests. Moreover, as already mentioned, a change of the axial load cannot be directly related to the presence of a crack in the tie-rod, since the axial load exhibits a high sensitivity to other physical variables which are not correlated to the state of health of the tie-rod (e.g., temperature).
Recently, a strategy for damage identification in tie-rods was presented in [7], where an analysis of the tie-rod dynamic behaviour is shown and the estimate of the flexural compliance is proposed as a way to highlight the differences between a healthy and a damaged tie-rod. The potential of the approach was proved with a laboratory test where the effect of a breathing crack was reproduced by a bolted joint between two beam sections. The approach can highlight the presence of a crack but it requires a comparison with a reference tie-rod under known health conditions. Moreover, environmental and operational effects were not considered, thus the validity of the approach in field applications was not verified.
To overcome the current limitations, in previous works [22,23], the authors proposed the adoption of a vibration-based data-driven approach that does not require knowledge of any physical variable of the problem, that does not require specific in situ tests and that can be adopted continuously, without human supervision. The statistical pattern recognition approach was adopted, framing the problem of tie-rod damage detection as a multivariate outlier detection one, by adopting damage feature vectors based on modal parameters. The motivations behind the potential of the strategy are described in [22]: the most important point is that the effects of damage can be distinguished from those of the environmental variations, because they cause different changes in the pattern of modal parameters of multiple vibration modes. Consequently, a damage index was defined using a well-known multivariate metrics, the Mahalanobis squared distance (MSD). The proposed damage index showed the potential to be used for automatic outlier detection.
Starting from these observations, in [23] the authors further developed the process, focusing on the automatisation of the algorithm. A multivariate damage feature based on a collection of eigenfrequencies associated with multiple vibration modes was adopted and an automatic data cleansing algorithm was proposed, to develop a robust algorithm that can be used in real applications. The proposed strategy was validated on long term monitoring data, coming from a one-of-a-kind test case in the research field of SHM: nominally identical full-scale tie-rods, monitored for several months under the effects of uncontrolled operational and environmental conditions. Moreover, while the majority of works found in the literature consider simulated damage scenarios, a peculiarity of this work lies in the validation on data acquired during an ongoing real damage, carefully generated through a corrosion process, which realistically evolved throughout different months. The results were promising and allowed for the automatic detection of damage coming from two corrosion attacks in two different sections of the tie-rod.
This work proposes a different approach that is aimed at accounting for the dispersion associated with the data acquired under real operating conditions, through the adoption of Gaussian Mixture Model (GMM). GMM is a class of finite mixture models widely used for unsupervised data clustering, successfully adopted for applications related to sound/speech recognition [24,25] and image processing [26,27]. In the context of SHM, GMMs have been adopted for crack detection in reinforced concrete [28][29][30] and laminated composites [31] based on Acoustic Emission technique. Furthermore, GMMs have been adopted to automatically cluster features extracted from Lamb wave signals in both the time and frequency domains for aircraft wing spar damage detection under time-varying boundary condition [32]. An application to bridge monitoring is that of [33], where a GMM-based algorithm for stay cables condition monitoring is presented. The authors of [33] used GMM on a damage feature which is the ratio between the tension of multiple cables, measured through load-cells. The authors of [34] proposed the adoption of GMM based on coefficients of an autoregressive model with moving average used to fit simulated vibration data of the ASCE Benchmark Structure [35]. In [36], GMM were adopted on eigenfrequencies of the Z-24 Bridge in Switzerland, observing that the approach can outperform other state-of-the-art algorithms (principal component analysis or Mahalanobis squared distance) in the presence of non-linear effects caused by the operational and environmental conditions. As mentioned, this paper presents a novel GMM-based damage index for SHM of axially loaded beam-like structures. To summarize, the main challenges are related to the development of a damage detection strategy that: • Does not require knowledge of physical variables (above all, does not require axial load estimate); • Requires a simple and cost effective set-up; • Is automatic; • Is validated under real environmental and operational conditions; • Is validated in the presence of real damage.
The innovative aspects of this research are: • The unsupervised data clustering approach to damage detection is applied for the first time to the case study of tie-rods; • A comparison between GMM-based and MSD-based approaches is shown in the presence of a real damage condition which evolves over time. It will be shown that the proposed GMM-based approach outperforms the benchmark one, based on the classical MSD, both in terms of sensitivity and uncertainty associated with the results.
The paper is organized as it follows. In Section 2, the theoretical background and the experimental set-up are described. The damage feature based on tie-rod eigenfrequencies and the benchmark damage index are recalled in Sections 2.2 and 2.3.1, respectively. The theoretical background on GMM and the novel damage index are described in Section 2.3.2. Results are presented in Section 3 and discussed in Section 4. The general conclusions are drawn in Section 5.

The Experimental Case Study
In this section, the case study is introduced. The data that will be discussed in the following come from an experimental set-up located in the Mechanical Engineering laboratory at Politecnico di Milano, in Italy. The set-up (see Figure 1) is composed by a series of nominally identical tie-rods, monitored for long time by an SHM system, for research purposes (e.g., [22,23,37,38]).
The tie-rods are made of aluminium and they are characterized by a free length W of 4 m, and a cross-section of 0.015 mm × 0.025 mm. At the two ends of the beam, there are clamps made from two steel plates, in contact with the upper and lower face of the beam, held together by bolts. During the installation, the bolted joints of one of the two clamps (clamp 1 in Figure 1) are tightened, while those of the other clamp (clamp 2 in Figure 1) are left loose: in this way, the beam is not fully constrained along the axial direction, allowing traction to be applied. After the axial load is applied through a tensioner, also the bolted joints of clamp 2 are tightened, obtaining a "clamped-clamped" configuration.
The beams are equipped with general purpose industrial accelerometers, model PCB603C01, characterized by a sensitivity of 10.2 mV/(m/s 2 ) and a full scale of ±490 m/s 2 . Vibration data, laboratory temperature and the axial load (which is measured through a calibrated Wheatstone full bridge) are acquired with NI9234 modules, with anti-aliasing filter on board, at a sampling frequency of 512 Hz. The resulting bandwidth of 200 Hz is appropriate to acquire the range of frequency significantly excited by the operational environment.
Indeed, the tie-rods are subject to environmental vibrations, in a laboratory where human activities take place and where other testing benches and machineries are working. The operating environment usually provides a broadband excitation that significantly decreases above 200 Hz. The temperature is not controlled: as an example, during the first year of monitoring, the temperature ranged from a minimum of 6°C during winter, to a maximum of 29°C during summer, with minimum and maximum daily thermal excursions approximately equal to 3 and 8°C, respectively. The operating loads and the temperature conditions were intentionally uncontrolled, to obtain a challenging environment to test the developed unsupervised learning damage detection algorithms. On some of the specimen, a chemical attack was carried out to intentionally introduce a corrosion process which progressively evolved through several months. A general corrosion caused a decrease of the nominal section of the tie-rod, due to a progressive electrochemical reaction between the metal and concentrated solutions of sodium hydroxide and hydrofluoric acids (more details on the corrosion process can be found in reference [23]). The data acquired during the deteriorative process will be used in Section 3, to compare the damage detection indexes.
As it will be explained in Section 2.2, this work focuses on a vibration based damage feature based on tie-rod eigenfrequencies. Before going into detail of the theory, some aspects are discussed about the eigenfrequency database that will be used in the following. Different operational modal analysis algorithms can be adopted to estimate the tie-rod modal parameters. Since the proposed approach is based on eigenfrequencies only, the choice was for the adoption of a simple single-output approach. This choice was made to develop a damage detection strategy that requires a very simple implementation and that could be easily translated to real applications. Eigenfrequencies were identified from the power spectrum of the response, adopting a single-degree-of-freedom (SDOF) modal identification technique [39], around every resonance of the tie-rod. The best fitting between the experimental power spectrum of the response, and the analytical power spectrum of the response of an SDOF system excited by white noise was carried out (see Appendix A). For this specific case, the tie-rods showed lightly coupled modes, not closely spaced in frequency and not heavily damped; thus, such simple and fast approach is a viable solution. However, it must be pointed out that the theory and approach developed in the following could be adopted regardless the operational modal analysis algorithm chosen.
The identification was carried out using vibration data acquired by a sensor that was placed at a distance approximately equal to 1/10 of the free length. This position was chosen because it is not a vibration node for the first six bending vibration modes. Despite being close to the constraints, where the eigenvector components are generally low, the signal-to-noise ratio allowed a sufficiently stable identification of the eigenfrequencies. The frequency averaging approach [40,41] was adopted to calculate the experimental power spectrum every hour, with a duration of the sub-records for the averaging procedure equal to 40 s (with an overlap of 50%), and adopting a Hanning window.
Being in an uncontrolled environment, it is not always that the averaging procedure allows obtaining a good reconstruction of the power spectrum, mainly due to the presence of disturbances. In such cases, the best fitting procedure may fail or converge to wrong solutions. For this reason, an automatic data cleansing approach was used to automatically discard the corrupted eigenfrequency estimates. The data cleansing algorithm is presented and discussed in detail in Ref. [23]; however, since the data adopted in the following are those obtained after the removal of wrong identifications, some details are provided also here, for the sake of completeness.
The data cleansing algorithm is performed at two stages, where two checks are carried out. The first check is carried out every time an eigenfrequency is identified, considering vibration modes one at a time. Outliers are discarded when a poor fit is observed between the analytical power spectrum of the response of a single-degree-of-freedom mechanical system excited by white noise [39] and the experimental power spectrum, estimated through Welch's approach [40,41]. This can be performed by setting an acceptance threshold on the value of the R 2 coefficient (identifications associated with an R 2 < 0.9 were discarded in this case).
At the second stage, multiple observations of the eigenfrequencies related to all the considered vibration modes are considered, assuming that damage-related phenomena are characterized by a long-term evolution, i.e., they do not significantly evolve in a short period, e.g., one or two weeks. From analytical formulations [42,43], every tie-rod squared eigenfrequency is linearly dependent on the axial load and, consequently, squared eigenfrequencies of vibration modes are linearly dependent to each other, if the axial load is the only changing variable. This linear dependency between couples of squared eigenfrequencies is exploited, removing observations that significantly deviate from the mentioned linear trend.

Vibration-Based Damage Feature for Beam-Like Structures
In the statistical pattern recognition approach [2], each pattern is a collection of a number C of parameters, or features, extracted from the monitoring data and can be thought of as a point in a C-dimensional space. Feature selection is a very important step: the goal is to select a damage feature that allows pattern vectors associated with different structural states to occupy disjoint regions in the C-dimensional space, such that decision boundaries can be established to detect abnormal structural conditions [44]. This means that a good damage feature must be highly correlated with the severity of damage, and, ideally, lowly correlated to confounding effects, i.e., environmental and operational variations. Moreover, it is important to assess that little changes in the structural condition reflect into detectable changes in the damage feature.
Modal parameters are arguably among the most commonly adopted vibration-based damage features [45], mainly eigenfrequencies and mode shapes [46,47] (few successful examples about the use of modal damping as a damage feature can be found in the literature [48], mainly due to the complexity of the damping mechanism and the uncertainty associated with its estimate [49,50]). This paper focuses on eigenfrequencies, which can be easily estimated relying on few sensors and are usually less affected by experimental noise with respect to mode shapes [51]. The main drawback of adopting eigenfrequencies as damage features is related to the fact that they are highly sensitive to environmental variations [52]. Indeed, temperature variations can cause changes in eigenfrequency values which are greater than those caused by damage at an early stage. Moreover, usually only few low eigenfrequencies can be identified, which might not be sensitive enough to local damage (e.g., the presence of a crack), because local damage reflects on high vibration modes [2].
In the paper [22], the authors showed how a feature vector composed by the eigenfrequencies of multiple vibration modes of a tie-rod can be used to detect damage and in paper [23] the conclusions were validated in the presence of real damage. The multivariate feature vector f based on eigenfrequencies is defined as it follows: where f m (with m = 1, · · · , M) are the eigenfrequencies of M considered vibration modes.
The key point of the proposed approach is that eigenfrequencies are used to synthetically represent the state of the monitored tie-rod, since they are representative of the physical variables that mostly influence its dynamic behaviour (e.g., the axial load). Although these variables change due to environmental and operational conditions, variations associated with the confounding effects (e.g., temperature) give rise to different patterns with respect to those due to damage, in the multivariate feature space. As an example, the eigenfrequencies of the first three bending vertical modes of a healthy tie-rod are considered: a decrease of temperature would cause an increase in the values of all three the eigenfrequencies and the lower the vibration mode considered, the higher the effect [22]. If the temperature does not change but damage (e.g., a reduction of cross-section) is present at midspan, only the eigenfrequencies of the odd vibration modes would change (midspan is a vibration node for the even vibration modes) and the higher the vibration mode considered, the higher the effect [22]. The damage feature f can be adopted to exploit this observation: in the multivariate feature space defined by the eigenfrequencies of multiple vibration modes, patterns related to temperature variations differ from those related to damage. Thus, the damage feature f allows for the separation between damaged and undamaged states even in the presence of environmental and operational variations, and, for this reason, it will be adopted to define a novel damage index based on unsupervised learning data clustering, as it will be discussed in the next subsection.

Damage Indexes Based on Unsupervised Learning Approach
Once the damage feature is selected, a training phase is needed to learn the relationship between a pattern and the health condition of the monitored structure. Labelled data must be used for this scope, i.e., data for which the health state of the structure is known. If training data are available for all structural conditions (healthy as well as damaged), the problem is in the area of supervised learning. In the context of real operating structures, training data related to damage conditions are often not available. In this case, the problem is of unsupervised learning: the objective is learning the intrinsic relationship between the pattern and the normal structural condition only, defining a single class (undamaged state) and then test whether new data are still consistent with that class. According to this approach, damage detection can be framed as a problem of outlier detection (or novelty detection) or, alternatively, of unsupervised data clustering. These two approaches are compared in this work, where a novel damage index based on GMM is introduced. In more detail, the benchmark model will be the one based on the classical MSD, which is often adopted for unsupervised learning outlier detection, and it will be introduced in Section 2.3.1. The novel damage index will be introduced in Section 2.3.2, where the basics of GMM are also recalled.

The Benchmark Approach Based on Outlier Detection
The outlier analysis is based on computing a discordancy measure between the damage feature calculated on new data and a reference training set of damage features calculated when the structure is in what is assumed to be a healthy condition. If the discordancy measure exceeds a threshold, an abnormal condition is detected. The case of interest for this work is the multivariate outlier detection, since f is a feature vector. In case the damage feature is multivariate, the discordancy measure is given by the MSD. For a generic multivariate feature vector f new , its MSD from a generic matrix [ f ] where every row contains an observation of the feature vector, can be calculated according to the following formula [53]: where µ f and Σ f are the multivariate mean vector and the covariance matrix of [ f ] respectively, and the suffix "−1" means the inverse.
In case the MSD is adopted, a damage index can be obtained calculating the MSD between the candidate feature vector f new and the baseline matrix [ f ] base containing the normal condition features, i.e.,: Using the Gaussian assumption for the normal condition, the threshold can be calculated in terms of chi-squared-statistic [54] or adopting a numerical method. The latter approach is here adopted, based on the Monte Carlo method described in [53], and consisting in the following steps:

1.
Populate a (R × C) matrix with randomly generated numbers from a zero-mean and unit-standard-deviation normal distribution, where R is the number of samples in the baseline matrix and C is the number of variables of the damage feature.

2.
Calculate the MSD between every row of the matrix and the matrix itself and store the maximum distance.

3.
Repeat the first two steps for a high number of trials (e.g., 1000) and store the resulting maxima in an array. The critical values for 1% or 5% tests of discordancy are given by the MSDs in the array above which 1% or 5% of the trials occurs, obtaining the socalled "inclusive threshold" t inc (i.e., the presence of data coming from the damaged structure in the baseline set is admitted).

4.
If the baseline set only includes data coming from the undamaged structure (as in the case considered in this work), the so-called exclusive threshold can be calculated according to the following expression [2]: MSD-based damage detection is one of the most popular methods in SHM due to its ease of use and computational efficiency [55,56]. One of the major advantages of the MSD is that it can be used to filter out the environmental effects on damage features, provided that the baseline set contains a wide range of environmental conditions [57].

The New Approach Based on Unsupervised Data Clustering
A new approach to damage detection in tie-rods is proposed in this section, by framing the problem as unsupervised data clustering. As mentioned above, a damage feature can be thought of as a point in a C-dimensional space. If an appropriate damage feature is selected, in a way that feature vectors associated with different structural states occupy disjoint regions in the C-dimensional domain, a data clustering approach can be adopted to carry out damage detection. The main premise is that there is a migration of groups (or clusters) of feature vectors as damage occurs to the structure. In this work, for the first time in the literature of tie-rod damage detection, the deterministic physics-based behaviour of tie-rod natural frequencies is exploited: a new multivariate domain is defined where the existence of clusters of damage-related data can be automatically spotted, setting the stage for unsupervised data clustering algorithms to be successfully adopted. Before going into details of the specific case-study, a brief recall of the theory related to the adopted algorithm (i.e., GMM) is presented.
Unsupervised data clustering can be carried out according to different approaches, e.g., statistical, neural and syntactic methods [58]. As all SHM problems are subject to various degrees of uncertainty, the statistical approaches appear to stand out as a natural choice for damage detection in real operating structures. Finite mixtures allow addressing the problem in a probabilistic way, since they inherently account for the uncertain nature of the problem. According to these approaches, every observation of the damage feature is considered as generated by one among a set of alternative random sources. Inferring the parameters of these sources allows one associating a given observation to the probability to have been generated by one specific source. In this way, data can be automatically clustered [59].
A GMM is a parametric probability density function made by a weighted sum of a number K of Gaussian densities (as it will be shown in the following, this probability density function is particularly suitable to fit the probability density function of the adopted damage feature). In the following, the generic matrix [x] is considered, which is a random sample set composed by a number R of independent random samples. A generic C-dimensional sample in the sample set is indicated by x r , with r = 1, . . . , R; x r is a row of the matrix [x] and it is one observation of the feature vector. By assuming that x r follows a GMM, the probability density function can be written as a mixture of Gaussian distributions according to the following expression: where θ i are the weights and each N i x r |µ i , [Σ] i is a C-variate Gaussian density defined as it follows: fully characterized by the i-th mean vector µ i and i-th covariance matrix [Σ] i . The weights θ i must satisfy the conditions θ i ≥ 0 for i = 1, . . . , K and ∑ K i=1 θ i = 1. For notation convenience, all the parameters that specify the mixture are collected in the set of parameters λ, defined as it follows: In order to adopt GMM for unsupervised data clustering, a set of parameters λ must be estimated so that the probability density function of the GMM best fits the distribution of the samples in the sample set [x]. The maximum likelihood estimate of the parameters can be obtained by maximizing the following expression: p(x r |λ) (8) which is the likelihood of the GMM, given a population containing a number R of feature vectors (a (R × C) matrix [x]), assuming independence among different observations. The maximization of the term in Equation (8) cannot be carried out analytically and, for this reason, an iterative procedure called expectation maximization (EM) algorithm [60] (see Appendix B) is adopted.
Now that the basics of GMM have been recalled for a general case, the proposed damage index is discussed for the specific application. First, the assumed hypothesis is that squared eigenfrequencies of different vibration modes are related to each other with a linear relationship, if the axial load is the only changing variable. If a number M of vibration modes is considered, the baseline matrix set where 1 is an all-ones column vector of the same size as s base m or s base 1 . Once the coefficients a 1m and b 1m are known, the residuals of the linear regression can be calculated according to next equation: The same procedure can be applied to the data contained in a matrix [ f ] new , containing the eigenfrequencies of the M considered vibration modes, in a period different from that of the baseline. In this case, the coefficients a 1m and b 1m , which are those estimated by considering the baseline set, will be used to calculate the residuals, i.e.,: finally obtaining a matrix [ε] new . The residuals calculated on the baseline set and those calculated on the new data can be put together in the same matrix [ε], according to the next expression: which is a population containing a number R of observations of the C-dimensional feature vector ε r , where C = M − 1.
The problem of damage detection can now be framed under a probabilistic point of view, studying the probability density function of the generic C-dimensional sample ε r . When no damage is present, the hypothesis is that the probability density function of ε r follows a single C-variate Gaussian density, i.e., p(ε r |γ) = N ε r |µ, [Σ] (13) where: are the parameters that maximize the likelihood of the model: When the set [ε] new contains data coming from the damaged structure, the hypothesis is that the probability density function of ε r follows a GMM composed by a mixture of two C-variate Gaussian densities: one that best fits the baseline residuals [ε] base and the other that best fits the damaged set [ε] new . According to Equation (5), ε r has a probability density function that can be written as: fully characterized by the parameters which can be estimated adopting the EM algorithm. The parameters contained in λ maximize the likelihood of the GMM made by two Gaussian distributions given the population of features [ε], which can be calculated as: The validity of the two hypotheses (single versus double C-variate Gaussian distribution) can be assessed by comparing L( γ) and L( λ). To do so, the parameters γ and λ, estimated adopting the EM algorithm, are used to calculate L( γ) and L( λ) given the data [ε] through Equations (15) and (18), respectively. For numerical convenience, the logarithm is used to transform products of potentially small likelihoods into a sum of logarithmic values, which are easier to distinguish from 0, in computation. Furthermore, optimization algorithms often search for the minima of functions; thus, the negative log-likelihood values will be adopted hereinafter, defined by the following Equations (19) and (20).
The general idea is that if L( λ) is significantly higher than L( γ) or, equivalently, L 2 is significantly lower than L 1 , there is a high likelihood of damage. To show the validity of the approach, data coming from the experimental set-up are now discussed. In the following, a real damage condition is considered, which is that associated with the effects of general corrosion close to the constraints.
In this case, three eigenfrequencies are considered, which are those associated with the bending vertical modes number 4, 5 and 6. The use of three eigenfrequencies (M = 3) returns a matrix of residuals with two columns (C = M − 1 = 2), which allows for a convenient representation adopting two-dimensional scatterplots.
First, the set of [ε] base is considered, with data referring to 1646 hours collected in the periods 19 December 2019 to 3 January 2020, and 22 April 2020 to 14 August 2020, after the adoption of the data cleansing algorithm described in [23]. A scatterplot of ε base 12 versus ε base 13 is reported in Figure 2. Below the x-axis and to the left side of the y-axis, the histograms of the data of ε base 12 and ε base 13 are reported, respectively. The histograms of Figure 2 suggest that the residuals in ε 12 and ε 13 are Gaussian distributed, as confirmed by a visual comparison between the empirical cumulative distribution function (CDF) and the Gaussian CDF for ε 12 and ε 13 , reported in Figure 3a,b, respectively. Thus, the data set contained in [ε] base is characterized by a bi-variate Gaussian probability density function.  The proposed damage detection strategy is explained through three figures (Figures 4-6), which are representative of three different tie-rod conditions. Each of the three figures is composed of three sub-plots, "a", "b" and "c". The sub-plot "a" shows the data in the matrix [ε], adopting black crosses and red triangles to indicate data belonging to [ε] base and [ε] new , respectively. The sub-plot "b" shows the probability density function of a single bivariate Gaussian distribution, adopting contour lines to help visualizing the single centroid of the distribution; the value of L 1 is reported above the figure (in the following, L 1 and L 2 will be used to define a damage index). The sub-plot "c" shows the probability density function resulted from the adoption of the EM algorithm, again adopting contour lines to help visualizing the centroids of the two bi-variate Gaussian distributions composing the mixture. Data are divided into two clusters: data belonging to "cluster 1" are represented with black crosses, while data belonging to "cluster 2" are represented with red triangles. The value for L 2 is reported above the figure.
When data contained in [ε] new are associated with a period when no damage is present (see Figure 4), the hypothesis of a single bi-variate Gaussian distribution that was observed for [ε] base is still valid. Indeed, data of [ε] new seems to be generated by the same probability density function of [ε] base , as confirmed by the fact that red triangles and black crosses are overlapped and mixed in Figure 4a. By comparing the probability density functions of the single bi-variate Gaussian distribution (Figure 4b) with that of the mixture (Figure 4c), the clustering solution proposed in Figure 4c seems to be due to over fitting. Moreover, it is worth noting that the value of L 2 is almost the same of that of L 1 (L 2 is 0.16% smaller than L 1 ).
Data presented in Figure 5 are associated with the presence of damage at an early stage. By a visual check of data in Figure 5a, the hypothesis that data of [ε] new belong to a different cluster with respect to that of the baseline data seems more likely. The EM algorithm proposes the clustering reported in Figure 5c. The value of L 2 is a 2% smaller than that of L 1 .
When data in [ε] new are associated with a period when a severe state of damage is present, points of [ε] new and [ε] base are disjointed (see Figure 6a). The EM converges to the solution presented in Figure 6c, showing that data in [ε] new and [ε] base belong to two different clusters. In this case, the value of L 2 is a 9% smaller than that of L 1 .
The three cases discussed here suggest that the strategy can be automatized by using a measure of the likelihood increment obtained by adopting a mixture of bi-variate Gaussian distributions in place of a single bi-variate Gaussian distribution to model the probability density function of [ε].  According to this idea, in a more general framework, the evaluation of the state of the structure can be carried out by assembling a matrix [ε], running the EM algorithm to fit a GMM with two components and calculating the values of L 1 and L 2 . A novel damage index is defined, according to the following Equation (21): LI can be adopted to quantify the likelihood increment obtained by passing from a single C-variate Gaussian distribution to a mixture of two C-variate Gaussian distributions (the procedure to calculate the index LI is described in the flowchart of Figure 7). The difference between L 1 and L 2 is normalized by dividing its value by L 2 , which is that associated with the model that is always able to best fit the underlying probability density function and, thus, set as reference. Indeed, the mixture can also approximate the single bi-variate Gaussian distribution (see the contour lines of Figure 4c as an example). On the contrary, the model based on a single C-variate Gaussian distribution fails to describe the probability density function of [ε], when more than one cluster exists in the data (see the contour lines of Figure 6b as an example). In the following section, DI (see Equation (3) in Section 2.3.1) and LI will be compared, to assess which approach can allow for an earlier damage detection, considering long term monitoring data and the effects of real damage on a tie-rod.

Results
The two damage indexes are compared considering two case studies. The first is the most challenging one since damage occurs close to the constraints, at a distance equal to 1 10 W (W is the beam free length) from one of the fixed ends of the tie-rod. Indeed, when damage is close to the constraints, the eigenfrequencies are less sensitive, as already observed in [22]. The damage scenario is associated with the real effects of the corrosive attack induced on the tie-rod surface, that affected a length of the tie-rod of approximately 5 cm. Some pictures of different stages of the corrosion process are reported in Table 1, and they are labelled with letters A1, B1, C1 and D1. To provide not just a qualitative picture, but also a quantitative one, the different cross-section reductions ∆h are reported (e.g., ∆h = 6% means a reduction of 6% of the cross-section height, measured at the centre of the corroded area).
Data are divided into "Baseline", "Validation" (i.e., data not included in the baseline set, but still related to a condition when no damage is present) and "Corrosion" sets, according to the timeline of Figure 8. For both damage indexes, the results presented are obtained considering the eigenfrequencies associated with three vibration modes, more specifically those of the fourth, fifth and sixth bending vibration modes in the vertical plane (see Figure 9). It is recalled that the eigenfrequencies are obtained by adopting operational modal analysis, thus exploiting the excitation coming from the operational environment (see Section 2.1).    (c)    Figure 10. Damage index LI for tie-rod 1.
In Figure 11, a comparison with the benchmark approach (DI, based on the MSD) is shown. To allow for a direct comparison of the two trends, both DI and LI are normalized to the respective damage detection threshold and they are shown in Figure 9a and Figure 9b, respectively. In more detail, the damage threshold for DI is calculated according to the procedure explained in 2.3.1, adopting Equation (4), where t inc is obtained considering the critical value for 1% test of discordancy. For both LI and DI, black asterisks are used to represent the damage index; due to the scatter of DI, significantly higher than that of LI, also the trend obtained with a moving average (window of one day, overlap 1 hour) is reported with green points in Figure 9a, to allow for a clearer visualization and a direct comparison of the two trends. A zoom is made on the data sets "Validation" and "Corrosion": a black-dashed vertical line indicates the start of the corrosion attack and, again, vertical lines labelled according to Table 1  The results related to another case study are reported in the following. In this second example, a different specimen is considered, where the corrosive attack was induced at a distance equal to 5 8 W from one of the fixed ends. As for the previous case, different stages of the deteriorative process are reported in Table 2, together with the labels and cross-section reductions. Data are divided into "Baseline", "Validation" and "Corrosion" sets, according to the timeline of Figure 12. For both damage indexes, the results presented are obtained considering the eigenfrequencies associated with the fourth, fifth and sixth bending vibration modes in the vertical plane (see Figure 13). The comparison between DI and LI, both normalized to the respective damage detection threshold, is shown in Figure 14. Data are divided into "Baseline", "Validation" and "Corrosion" sets, according to the timeline of figure 12. For both damage indexes, the results presented are obtained considering the eigenfrequencies associated to the fourth, fifth and sixth bending vibration modes in the vertical plane (see figure 13). The comparison between DI and LI, both normalized to the respective damage detection threshold, is shown in figure 14.

Discussion
A first consideration that is worth mentioning is related to the eigenfrequency trends reported in Figures 9 and 13: a high variability can be observed, characterized by both daily and long-term trends related to temperature conditions. The magnitude of these fluctuations does not allow spotting the existence of damage, even when the eigenfrequency associated with the highest vibration mode (i.e., the sixth vibration mode) is considered. Indeed, even if the final stage of the corrosion process is considered (condition D1 in Figure 13c and condition D2 in Figure 9c), the eigenfrequency values are compatible with those observed in the previous year and they seem to follow the seasonal trend.
The effectiveness of the proposed approach is shown by the trend of LI in Figure 10: as it is possible to see, LI is exceeding the damage threshold when damage is present (red triangles in Figure 10). Damage is detected way before when the tie-rod is in condition B1 (see Table 1), which is remarkable considering that the corroded portion of the tie-rod is close to the fixed end. Furthermore, the blue circles associated with the validation set remain below the threshold, and this is consistent with the fact that the tie-rod is still in the same condition of the baseline set.
Coming to the comparison with the benchmark model, the damage index LI, based on the GMM, outperforms DI, based on the MSD, as it is shown in Figure 11. Indeed, DI is characterized by a higher variability (compare the black asterisks in Figure 11a with those of Figure 11b) which causes the damage index to be scattered around the damage detection threshold for a long time before it remains steadily above the threshold (approximately, between conditions C1 and D1). If the averaged trend of DI (green points in Figure 11b) is compared with LI (black asterisks in Figure 11b), it is clear how the novel GMM-based damage index has a higher sensitivity to damage, potentially allowing for a prompter damage detection. LI shows a clear trend that points out the presence of damage and its evolution over time, thanks to a higher robustness to the effects of environmental and operational variations.
The same remarks can be made regarding the second example here discussed, i.e., damage farther from the constraints. The adoption of a multivariate damage feature based on tie-rod eigenfrequencies allows spotting the presence of damage, with both DI and LI that cross the respective damage detection threshold, in Figure 14a and Figure 14b, respectively. LI is clearly detecting damage when the tie-rod is not yet in condition B2. It is worth noting that, at this stage, the damage effects are barely noticeable through a visual inspection (see Table 2). Once again, DI detects damage later than LI, staying constantly above the damage detection threshold approximately between conditions C2 and D2 (see Figure 14a).
In this second case study, the different performances between the two damage indexes can be clearly highlighted. First, if the averaged time-trend of DI is considered (green points in Figure 14b), fluctuations can be detected, while the time-trend of LI (black asterisks in Figure 14b) shows a more stable trend, again proving that the GMM-based damage index is more robust to environmental variations than the MSD-based one. Second, regarding the comparison between the trend of DI with no moving average and LI, DI is more sensitive to noise in the data than LI. If the black asterisks in Figure 14a are considered, the scatter of the results is partially associated with the uncertainty related to the modal identification (e.g., errors due to a poor signal-to-noise ratio of the original vibration data), which plays a significant role when the MSD is calculated between a single observation of the damage feature and the baseline set. On the contrary, when LI is estimated, the uncertainty is actually modelled by the GMM: the uncertainty in the data causes the spread of the damage features in the multivariate feature space (e.g., how wide the clouds of Figure 6 are), but does not reflect on the scatter of LI. Indeed, LI relies on the validity of the hypothesis related to the underlying probability density function, which significantly changes only when the deterministic effect of damage causes a migration of the centroid of the population of new data from the region occupied by the baseline data. So, even if a lower noise in the data would make it easier to distinguish different clusters, thus increasing the sensitivity of the method, the index LI is much more robust to noise in the data than DI.
To resume, the main advantages of the adoption of the GMM-based strategy are related to a higher sensitivity to damage and a lower uncertainty associated with the results with respect to the MSD-based strategy, as emerged from the experimental campaign. These characteristics are crucial when the damage is characterized by a progressive evolution in time, as in the case of corrosion. On the other hand, the main advantage of the MSDbased strategy with respect to the the GMM-based one is that it is potentially prompter in detecting fast events (e.g., a sudden failure of the constraint system). Indeed, if a single observation of the multivariate damage feature related to a fast damage event is considered, its effect can cause a high MSD from the baseline set which can immediately make DI cross the detection threshold. However, this could not be enough to cause a significant migration of the centroid of the population of new data with respect to the baseline cluster, until a sufficient number of damage-related observations is available, causing a delay between the occurrence of the failure and the moment when LI exceeds the threshold.

Conclusions
This paper presented a novel damage index for automatic damage detection in operating tie-rods, based on unsupervised learning data clustering through GMM. The damage index was compared with a benchmark one, estimated through one of the most adopted multivariate metrics for outlier detection: the MSD. Both damage indexes are calculated on the same database containing multiple observations of the feature vector composed by tie-rod eigenfrequencies. Although the MSD is effective in detecting damage, as already shown by the authors in previous works, the novel damage index based on GMM proved to be affected by a significantly lower scatter and to provide a higher sensitivity, which allows for an earlier damage detection. Indeed, the MSD considers each new observation of the damage feature against the entire reference population; thus, the results are characterized by a scatter which is mainly related to the variability of each identification of the damage feature. Consequently, the scatter is higher than that obtained with the GMMbased approach, where a population of data is tested against the baseline population. The GMM-based strategy, which can spot the onset of new clusters in the data in a completely unsupervised way, proved to be suitable to detect evolutive deteriorative phenomena, such as corrosion, outperforming the MSD-based approach. On the other hand, for other types of failure associated with faster phenomena (e.g., a sudden failure of the constraint system), the MSD-based approach should be more suitable to detect abrupt changes in the structural condition. A further development of this work will be a comparison of the two approaches on other types of damage. This will require carrying out new experiments, where other different damage mechanisms are introduced on the tie-rods, e.g., the failure of the constraints. Funding: The research presented in this paper has been funded by the Italian National Research Program, in a project named "Life-long optimized structural assessment and proactive maintenance with pervasive sensing techniques" PRIN 2017.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.
In the expectation step, the model parameters λ (z) are considered as fixed and the posterior probability is calculated for every feature vector x r and for every cluster i: which is the a posteriori probability that the feature vector x r belongs to the cluster i. The posterior probability of all the samples in each Gaussian component is calculated according to the following expression: In the maximization step, the model parameters are updated, obtaining λ (z+1) according to the following expressions: A commonly adopted stopping criterion is defined on the log-likelihood increment between two subsequent iterations, according to the next expression: where is a user defined tolerance (a tolerance of 1e-6 is adopted in this work); otherwise the steps are recursively performed until the maximum number of iterations is reached (a maximum number of iterations equal to 500 is used in this work). The EM algorithm is strongly dependent on the initialization, i.e., on the choice of model parameters λ for the first iteration. Common solutions are the use of multiple random starting points and then choosing the final estimate with the highest likelihood, or initializing by other clustering algorithms, e.g., the K-means algorithm [63] (the latter approach is that adopted in this work). Due to the high sensitivity to the initialization, the convergence of the EM algorithm is guaranteed to a local maximum but not to a global maximum [60].