The Extreme Function Theory for Damage Detection: An Application to Civil and Aerospace Structures

Featured Application: This paper presents an application of Extreme Function Theory (EFT) and Gaussian Processes (GPs) to perform mode shape-based damage detection and Structural Health Monitoring (SHM). Abstract: Any damaged condition is a rare occurrence for mechanical systems, as it is very unlikely to be observed. Thus, it represents an extreme deviation from the median of its probability distribution. It is, therefore, necessary to apply proper statistical solutions, i.e., Rare Event Modelling (REM). The classic tool for this aim is the Extreme Value Theory (EVT), which deals with uni- or multivariate scalar values. The Extreme Function Theory (EFT), on the other hand, is deﬁned by enlarging the fundamental EVT concepts to whole functions. When combined with Gaussian Process Regression (GPR), the EFT is perfectly suited for mode shape-based outlier detection. In fact, it is possible to investigate the structure’s normal modes as a whole rather than focusing on their constituent data points, with quantiﬁable advantages. This provides a useful tool for Structural Health Monitoring, especially to reduce false alarms. This recently proposed methodology is here tested and validated both numerically and experimentally for different examples coming from Civil and Aerospace Engineering applications. One-dimensional beamlike elements with several boundary conditions are considered, as well as a two-dimensional plate-like spar and a frame structure.


Introduction
The identification and localisation of damage in one-, two-, or three-dimensional structures through their mode shapes are widespread in the Structural Health Monitoring (SHM) community. It is well known that by inserting a localised discontinuity, the mode shapes diverge from their usual deflection path [1]. This principle has been extensively applied to investigate damage-induced variations in the mode shapes slope or curvature for both 1-dimensional beam-like [2] and 2-dimensional plate-like [3,4] structures. In this regard, a review of classic approaches can be found in Reference [5]. These approaches consider modal curvatures [6], mode shape rotations [7], and/or several other Damage-Sensitive Features (DSFs) based on the structure's eigenvectors. In particular, mode shapesbased methods are preferred over other modal parameters for damage localisation, since they inherently have the spatial resolution needed for this specific task [8,9]. However, even having established that the mode shapes and derived quantities can be exploited as a reliable and spatially dense DSFs, it remains to define how these features can be used for anomaly detection.
The changes in the mode shapes can be detected, e.g., through a Machine Learning (ML) process, trained exclusively on the mode shapes extracted from the current state of

Gaussian Processes and Extreme Function Theory
In the specific case of interest, the problem can be stated as follows. Let i = 1, . . . , n be the index of n examples taken from a population of functional data. Each i-th example is an identified mode shape ϕ i , i.e., a vector of j = 1, . . . , n i modal coordinates ϕ i,j observed at a finite number of corresponding spatial coordinates x i,j . From a theoretical standpoint, there is no requirement that the n mode shapes have the same number of coordinates nor the same location. This is the first strong point of EFT in comparison to EVT, as the functions can be collected from different sensor arrangements-even if, for practical concerns, it is more reliable to consider a fixed number and position of the output channels. Importantly, x i,j can be a scalar for a 1-dimensional beam element, or a pair x i,j = ξ i,j , η i,j for a 2-dimensional plate-like structure (and so on). For simplicity's sake, all the rest of this discussion will focus on the 1-dimensional case but the formulation for higher dimensionalities can be straightly derived.
There are two subsequent aims. The first one is to build the "normality" model M from the pairs of training data {x i , ϕ i }, collected from a given structure excited in the current conditions (assumed as a baseline for outlier detection). Once the model is defined, the second step is to discern if a test set x * i , ϕ * i (in this case, i = 1, . . . , n * ) belongs or not to the tails of the distribution M. The first part can be achieved with the classic GP Regression, which is here briefly recalled for completeness, while the EFT comes in the second part.

The Gaussian Process (GP) Regression
Being the mode shapes ϕ a function of the space variable x, it is possible to define a GP prior over this latent variable, in the general form not dissimilarly from how a normal distribution is uniquely defined by its mean and variance, the GP of the general process f(x) can be defined uniquely by its mean function and on the covariance matrix where E[·] denotes the expectation. Throughout this whole discussion, a zero-mean and squared-exponential (SE) covariance function will be applied (the exact formulation will be discussed later; other less frequent options exist and a discussion on the topic can be found in Reference [29]).
If f = f(x) indicates the set of the function values at the training points X (note that the training vectors have been assembled into a matrix form) and f * a set of test outputs at a new set of points X * , according to the zero-mean prior it is possible to define f f * ∼ N 0, K(X, X) K(X, X * ) K(X * , X) K(X * , X * ) , where K(X * , X) ∈ R n * ×n is a matrix with the terms defined by the covariances evaluated at all pairs of test and training points (the same applies for K(X, X), K(X, X * ), and K(X * , X * )). The posterior distribution is defined by conditioning the joint prior on the observations, resulting (for an ideally noise-free scenario) in p(f * |X * , X, f) ∼ N K(X * , X)K(X, X) −1 f, K(X * , X * ) −K(X * , X)K(X, X) −1 K(X, X * ) .
By considering the measurement noise as i.i.d. Gaussian noise with variance σ 2 n , one can define ϕ = f (x) + (6) thus, Equation (5) becomes with the predictive mean and the predictive covariance matrix At this point, it is evident that the only step remaining before obtaining the normality model is to define its covariance and the added noise. It was already assumed that this latter depends solely on σ 2 n . In the simpler (i.e., 1-dimensional) case, the aforementioned SE covariance function can be formulated as where k ϕ is the covariance of the noisy target set ϕ, σ 2 f is the variance of the function, l is the length-scale parameter, δ p,q is the Kronecker delta, and x p , x q are the locations of the training data sets corresponding to the respective indices. σ 2 f , l, and σ 2 n are the free parameters who need to be set from the data and are collectively known as the GP Regression hyperparameters θ. Their setting requires an optimisation, which is generally performed by maximising the log marginal likelihood with respect to the hyperparameters. However, in general, the minimisation, rather than the maximisation, is performed; therefore, it is more practically convenient (without any theoretical drawback) to define directly the negative log marginal likelihood (NLML) as where K ϕ indicates the covariance matrix of the noisy test set. The iterative minimisation is performed through a Conjugated Gradient Optimisation (CGO) algorithm; the process stops when the difference in NLML values are smaller than a set tolerance. At that point, the hyperparameters can be considered as the best fitting for the SE covariance function on the training data considered. The GP Regression over the training data can be therefore easily achieved and the metamodel M defined.

The Extreme Function Theory and the Complete Procedure
Consider the posterior distribution as defined in Equation (7), set accordingly to the hyperparameters optimised as described in the previous Section. By definition, a GP distribution over functions is a multivariate normal (Gaussian) distribution, thus the probability distribution z = p(f * |X * , X, ϕ) will follow the form where m * = m(f * ), K * = cov(f * ), and D is the dimensionality of x * (here onwards it is assumed that a single mode shape is tested at a time; even if it is theoretically possible to test jointly a set of more mode shapes at once, this will have little practical sense). Therefore, the whole test set x * i , ϕ * i is reduced to a single scalar value corresponding to its likelihood. In turn, this can be used to fit an EV distribution over its left tail (i.e., for a very low probability of occurrence, as expected for damaged conditions). The lower 10% of the validation dataset was used to this aim. Specifically, the Cumulative Distribution Function (CDF) is here fitted on low values of the Gaussian log probability ln(z). For the GEV minima distribution, the exact formulation of the CDF is where λ, δ, and γ represent, in the same order, the location, scale, and shape parameter of the GEV distribution, which needs to be inferred from the data. Importantly, the shape parameter γ controls the specific form of the limit minima distribution, which simplifies in a Gumbel distribution for γ → 0 , in a Weibull distribution if γ < 0, or in a Fréchet distribution if γ > 0 [30]. The estimation of the GEV parameters has been achieved utilising a Differential Evolution (DE) algorithm, as suggested in Reference [31]; specifically, the Self-Adaptive Differential Evolution (described in Reference [32]) was applied, considering 10 runs, 100 generations, and population size equal to 30. The Normalised Mean Square Error (NMSE) was set as the target cost function, minimised against the empirical CDF calculated according to the Hazen plotting position; an example is depicted in Figure 1. Having constructed the CDF it is then possible to define a threshold value ln z lim in correspondence of an arbitrarily chosen quantile α. Here in this study, α = 1% was imposed. This completes the description of the whole procedure.
where , , and represent, in the same order, the location, scale, and shape parameter of the GEV distribution, which needs to be inferred from the data. Importantly, the shape parameter controls the specific form of the limit minima distribution, which simplifies in a Gumbel distribution for → 0, in a Weibull distribution if < 0, or in a Fréchet distribution if > 0 [30]. The estimation of the GEV parameters has been achieved utilising a Differential Evolution (DE) algorithm, as suggested in Reference [31]; specifically, the Self-Adaptive Differential Evolution (described in Reference [32]) was applied, considering 10 runs, 100 generations, and population size equal to 30. The Normalised Mean Square Error (NMSE) was set as the target cost function, minimised against the empirical CDF calculated according to the Hazen plotting position; an example is depicted in Figure 1. Having constructed the CDF it is then possible to define a threshold value ln lim in correspondence of an arbitrarily chosen quantile . Here in this study, = 1% was imposed. This completes the description of the whole procedure.

Test of the Performances of the Procedure
To quantitatively estimate the performances of the described approach, this was tested on numerically simulated and experimental datasets with known damage conditions. For this purpose, three datasets -training, validation, and test -are defined. The training and the validation datasets, needed to preliminary set and validate the normality model offline, are made up of functions drawn from the undamaged conditions, while the

Test of the Performances of the Procedure
To quantitatively estimate the performances of the described approach, this was tested on numerically simulated and experimental datasets with known damage conditions. For this purpose, three datasets-training, validation, and test-are defined. The training and the validation datasets, needed to preliminary set and validate the normality model offline, are made up of functions drawn from the undamaged conditions, while the test set includes both damaged and undamaged states. For clarity's sake, the several steps described in the previous Subsection are graphically depicted in Figure 2. In comparison to a preliminary version of the algorithm described in Reference [28] and tested for simple cantilever beams, the following enhancements have been included.
Firstly, a data augmentation strategy was applied to overcome the limited amount of available measured data (numerical data were generated with the same scarcity as for the real data from the experimental case study). Dealing with the feature of transverse displacements, the latent function ( ) of Equation (6) is the resulting mean of all the original mode shapes and the resulting standard deviation indicates how data are spread out from the average. Both the mean and the standard deviation were evaluated for each channel, to consider the distribution for the specific measurement point and not of the whole dataset. After generating the desired number of noise-free copies, white Gaussian noise was added by iteratively reducing the Signal-to-Noise-Ratio (SNR) until the artificially simulated data fell completely within the set limits. In this case, the range was arbitrarily chosen to be ±2σ and not to cover the whole ±3σ range, since this restriction strengthens the "normality" condition (i.e., undamaged conditions, not far from the median value). It should be noted that data coming from the augmentation process are involved in the training and validation steps, while only original data are tested to assess the efficiency of the algorithm.
As a second improvement, a step of k-fold cross-validation was added, to avoid overfitting of model parameters estimations. In contrast to the usual validation procedure, the k-th holdout fold subset was used to train the Gaussian Process and estimate its hyperparameters, while keeping the remaining − 1 folds for the CDF estimation and model In comparison to a preliminary version of the algorithm described in Reference [28] and tested for simple cantilever beams, the following enhancements have been included.
Firstly, a data augmentation strategy was applied to overcome the limited amount of available measured data (numerical data were generated with the same scarcity as for the real data from the experimental case study). Dealing with the feature of transverse displacements, the latent function f(x) of Equation (6) is the resulting mean of all the original mode shapes and the resulting standard deviation indicates how data are spread out from the average. Both the mean and the standard deviation were evaluated for each channel, to consider the distribution for the specific measurement point and not of the whole dataset. After generating the desired number of noise-free copies, white Gaussian noise was added by iteratively reducing the Signal-to-Noise-Ratio (SNR) until the artificially simulated data fell completely within the set limits. In this case, the range was arbitrarily chosen to be ±2σ and not to cover the whole ±3σ range, since this restriction strengthens the "normality" condition (i.e., undamaged conditions, not far from the median value). It should be noted that data coming from the augmentation process are involved in the training and validation steps, while only original data are tested to assess the efficiency of the algorithm.
As a second improvement, a step of k-fold cross-validation was added, to avoid overfitting of model parameters estimations. In contrast to the usual validation procedure, the k-th holdout fold subset was used to train the Gaussian Process and estimate its hyperparameters, while keeping the remaining k − 1 folds for the CDF estimation and model testing. This inversion in the common paradigm of using the larger part of data for training step is due to the need to take into consideration only the lower 10% of the available data to estimate the minima form of CDF as defined in Equation (13). Subsequently, the training (TR), validation (VA), and test (TS) subsets were assembled from the given d structural configurations considered as the "normality" model (generally d = 1, the integer measurements data set) and d * structural configurations considered as damaged, each with n actual observations and n * artificial copies. Table 1 describes the procedure settings, as applied for all the numerical and the experimental case studies for consistency. Table 1. The procedure settings of the damage detection.

Number of Observations (n) 50 for the Experimental Case Study 100 for the Numerical Case Studies
Data Augmentation (n*) 500 As said, the test set contains both undamaged (ϕ und ) and damaged (ϕ dam ) mode shapes. When a ϕ und is classified as damaged, it is recorded as a false positive; on the other hand, misclassified ϕ dam are considered as false negatives. Otherwise, as pictorially described in Figure 3, the results are labelled as true positives and true negatives.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 7 of 23 structural configurations considered as the "normality" model (generally = 1, the integer measurements data set) and * structural configurations considered as damaged, each with actual observations and * artificial copies. Table 1 describes the procedure settings, as applied for all the numerical and the experimental case studies for consistency.

Number of Observations (n) 50 for the Experimental Case Study 100 for the Numerical Case Studies
As said, the test set contains both undamaged ( ) and damaged ( ) mode shapes. When a is classified as damaged, it is recorded as a false positive; on the other hand, misclassified are considered as false negatives. Otherwise, as pictorially described in Figure 3, the results are labelled as true positives and true negatives.

Numerical Simulations
The procedure is firstly validated on a numerically-simulated beam with different boundary conditions. Then, for more complex structures, the proposed approach is firstly validated on two other numerical datasets.

Numerical Simulations
The procedure is firstly validated on a numerically-simulated beam with different boundary conditions. Then, for more complex structures, the proposed approach is firstly validated on two other numerical datasets.
The first case study comes from an aerospace application and involves the High Aspect Ratio (HAR) aluminium spar of the XB-1 wing prototype [33]. The spar-portrayed in Figure 4a-can be seen as a thin, plate-like structural element with a peculiar shape. Due to its large flexibility, the spar may undergo large flap-wise oscillations [34,35] and it is therefore highly subject to the insurgence of fatigue damage, especially close to its clamped end and in its mid-length at the section where the width taper angle change. The simple Finite Elements (FE) model utilised here was realised in Ansys ® Mechanical APDL; 400 8-noded, 6-degrees-of-freedom-per-node quadratic shell elements were utilised. The main mechanical and geometrical properties are reported in Table 2 (the values of the parameters derive from preovius studies described in [36]). The model parameters were updated with vibrational data coming from laser Doppler vibrometer (LDV) and video acquisitions [36], extracted by means of the Fast Relaxed Vector Fitting (FRVF) approach [37,38]. Five output channels (located at the positions of the 4 Inertial Measurement Units and of the LDV target used during the experimental acquisitions, as represented in Figure 4a) were utilised to define the mode shapes. Several levels of reduced Young's modulus were applied to the portions highlighted in Figure 4b to simulate 20 scenarios [10]. In particular, the first two states (01 and 02) are intended to represent small perturbations from the undamaged baseline, where the variations are not sufficiently marked to be defined certainly as actual damage. This simulates the possible statistical variability of the identified mechanical properties; therefore, these states should preferably be labelled as false positives by a reliable (i.e., not hypersensitive) damage detection procedure. Indeed, it was found that, when trained on noiseless mode shapes only, the procedure does not discern these small perturbations from the normal (undamaged) model.
The second numerical case study comes from a civil engineering application and is intended to model a simple multi-storey frame structure, which will serve as the experimental benchmark in the next Section. The FE ( Figure 5) model, developed utilising the StaBil 2.0 MatLab Toolbox [39], was set as follows. 8 beam elements per column and per beam were used, totalling 72. The Timoshenko model was applied for all elements.
To mimic the experimental setup which will be described in the next section, the Young's modulus of the beams E beam was set as two orders of magnitude larger than its counterpart for the columns (E column ). This allows to approximate the structural system as a shear-type frame structure; therefore, a single channel per floor is enough to capture its main vibrational modes. Three output channels were considered as virtual sensors, located at nodes #4, #12, and #20 (indicated by the purple dots in Figure 5). The main technical details are reported in Table 3. The column tracts highlighted in green are the ones where the damage was inserted (as a stiffness reduction equal to 25% or 50%). ticular, the first two states (01 and 02) are intended to represent small perturbations from the undamaged baseline, where the variations are not sufficiently marked to be defined certainly as actual damage. This simulates the possible statistical variability of the identified mechanical properties; therefore, these states should preferably be labelled as false positives by a reliable (i.e., not hypersensitive) damage detection procedure. Indeed, it was found that, when trained on noiseless mode shapes only, the procedure does not discern these small perturbations from the normal (undamaged) model.      Figure 6 shows the results for two examples, the symmetrical pinned-pinned structural configuration and the asymmetrical clamped-pinned configuration. The method can be nevertheless applied to any statically determined or undetermined scenario. As it will be shown in the next Section, the same is valid for frame structures as well as 2dimensional or more complex structures. Two levels of artificially added noise have been considered-SNR = 60 dB and 40 dB, respectively. The first state depicted in Figure 6 corresponds to the undamaged structure. Consider that in this state, no damage was to be detected, so the successful damage identifications and false alarms are both zeroed. The remaining states are defined for nine crack locations, equally spaced of 10% l between 0 and l, and for a crack depth equal to 10% of the thickness (states 2 to 10) or 20% (states 11 to 19). It can be noticed how these states, being affected by larger damage, are more easily detected by the algorithm as foreseeable. Only the first mode shape is reported for conciseness, calculated at the nodes #10, #30, and #49 (numbered sequentially, left to right, with the beam discretised in 50 elements).
Appl. Sci. 2021, 11, x FOR PEER REVIEW 11 of 23 l, and for a crack depth equal to 10% of the thickness (states 2 to 10) or 20% (states 11 to 19). It can be noticed how these states, being affected by larger damage, are more easily detected by the algorithm as foreseeable. Only the first mode shape is reported for conciseness, calculated at the nodes #10, #30, and #49 (numbered sequentially, left to right, with the beam discretised in 50 elements). It can be appreciated that the rate of success strictly depends on the position of the damage. This issue is well-known for any mode shape-based approach. In general, this depends on the position of the mode shapes nodes and antinodes and can be solved by considering, e.g., different mode shapes at once [40].
For a baseline-based approach as the one proposed here, the normalisation to a unit maximum displacement also causes some regions of the damaged mode shape to be very close to the ones of the undamaged scenario, reducing the detectability of the damage. This arises from the eigenproblem being underdetermined and therefore the eigenvectors being defined only up to scale.
This difference between the normalised damaged and baseline mode shapes is here considered in terms of cumulative Euclidean distance (CED) and it is exemplified in Figure 7 (for better interpretability, the complemental similarity 1-CED is shown). The numerical example represents the 1st mode shape of a cantilevered aluminium box beam, modelled from the experimental case study analysed in Reference [9] (the concept can be extended to any structural configuration and set of boundary conditions and expanded to 2-and 3-dimensional structures). To make the damage effects more visible, an unrealistically large damage (a crack with deep equal to the 75% of the beam height) has been inserted at 20 equally spaced locations. The mode shapes are defined at any centimetre over the whole beam length (totalling 715 nodes). Figure 8 a,b represent the same concept adapted for the two structural configurations discussed previously. For indicative purposes only, damage equal to a 10% reduction in the beam stiffness was considered, roved between 0 and l at steps of 1%. Only for this graphic, the mode shapes have been defined using all the available nodes as output channels and with noise-free data, i.e., considering ideal conditions. It can be appreciated that the rate of success strictly depends on the position of the damage. This issue is well-known for any mode shape-based approach. In general, this depends on the position of the mode shapes nodes and antinodes and can be solved by considering, e.g., different mode shapes at once [40].
For a baseline-based approach as the one proposed here, the normalisation to a unit maximum displacement also causes some regions of the damaged mode shape to be very close to the ones of the undamaged scenario, reducing the detectability of the damage. This arises from the eigenproblem being underdetermined and therefore the eigenvectors being defined only up to scale.
This difference between the normalised damaged and baseline mode shapes is here considered in terms of cumulative Euclidean distance (CED) and it is exemplified in Figure 7 (for better interpretability, the complemental similarity 1-CED is shown). The numerical example represents the 1st mode shape of a cantilevered aluminium box beam, modelled from the experimental case study analysed in Reference [9] (the concept can be extended to any structural configuration and set of boundary conditions and expanded to 2and 3-dimensional structures). To make the damage effects more visible, an unrealistically large damage (a crack with deep equal to the 75% of the beam height) has been inserted at 20 equally spaced locations. The mode shapes are defined at any centimetre over the whole beam length (totalling 715 nodes). Figure 8 a,b represent the same concept adapted for the two structural configurations discussed previously. For indicative purposes only, damage equal to a 10% reduction in the beam stiffness was considered, roved between 0 and l at steps of 1%. Only for this graphic, the mode shapes have been defined using all the available nodes as output channels and with noise-free data, i.e., considering ideal conditions. (a) (b) Figure 8. The similarity between the damaged and undamaged mode shapes in terms of cumulative Euclidean distance (a) simply supported beam (b) clamped-hinged beam. Points with similarity equal or larger than 80% are highlighted in shades of red, above the dash-dotted line.

Results for the HAR Wing Spar (2-Dimensional) FE Model
The results for the damage scenarios depicted in Figure 4b are portrayed in Figure 9. Only the 1st mode shape is showed for brevity; the 2nd and 3rd mode shape returned similar results in terms of accuracy (with similar rates of occurrence for both false positives and false negatives). Artificial white Gaussian was added considering SNR = 60 dB and 40 dB.

Results for the HAR Wing Spar (2-Dimensional) FE Model
The results for the damage scenarios depicted in Figure 4b are portrayed in Figure 9. Only the 1st mode shape is showed for brevity; the 2nd and 3rd mode shape returned similar results in terms of accuracy (with similar rates of occurrence for both false positives and false negatives). Artificial white Gaussian was added considering SNR = 60 dB and 40 dB.
One can see that the method always successfully labels the undamaged and damaged states correctly for SNR = 60 dB, which is already a relatively large amount of noise.
Note that no damage was inserted in state 00; therefore, no damage was to be identified there, similarly to what expressed for the results of the beam structure (in Figure 6).
The states 01 and 02, which have too small variations of Young's modulus to be effectively considered as damaged, are rightly identified as false negatives. This proves the reliability of the method in discerning small perturbations, which may happen in real-life acquisitions due to the statistical variance of the acquisitions, from the actual damage. On the other hand, with even more white Gaussian noise artificially added to the input data (SNR = 40), the approach still correctly identifies the undamaged states, while it loses accuracy on the damaged mode shapes of the test set. Specifically, it remains able to recognise the most damaged situations (with a 25% or 50% reduction of Young's modulus) whenever the damage is inserted at the clamped base. As expected, also by comparison with previous works [10], the scenarios with the damage inserted solely at the mid-length cross-section (states 3 to 5) or on the sides of the base cross-section not fixed to the centre wing box (states 18 to 20). One can see that the method always successfully labels the undamaged and damaged states correctly for SNR = 60 dB, which is already a relatively large amount of noise.
Note that no damage was inserted in state 00; therefore, no damage was to be identified there, similarly to what expressed for the results of the beam structure (in Figure 6).
The states 01 and 02, which have too small variations of Young's modulus to be effectively considered as damaged, are rightly identified as false negatives. This proves the reliability of the method in discerning small perturbations, which may happen in real-life acquisitions due to the statistical variance of the acquisitions, from the actual damage. On the other hand, with even more white Gaussian noise artificially added to the input data (SNR = 40), the approach still correctly identifies the undamaged states, while it loses accuracy on the damaged mode shapes of the test set. Specifically, it remains able to recognise the most damaged situations (with a 25% or 50% reduction of Young's modulus) Regarding the computational efficiency of the algorithm, the whole simulation of the complete HAR wing dataset was performed in around 31 s. The non-optimised MatLab script elapsed on average 2.2 s to perform a single training fold, 3.9 s for CDF validated estimation and 0.02 s for the damage detection task on a single damage state, with no major differences among the 20 states, on an Intel ® Core™ i7-10750H CPU with 2.6 GHz base frequency and 15.8 GB of available RAM. Similar elapsed times were found for the other numerical and experimental cases.

Results for the Frame FE Model
The settings for the second numerical investigation are the same as described in Table 1. Figure 9 reports the results of the frame structure FE model for an increasing level of artificially added noise. Again, as for the previous case studies, consider that no damage was inserted in the baseline (state 1). For what concerns the noise level, the results in Figure 10b show how, for a realistic noisy scenario consistent with current measurement technology, a good performance is reached in the states with high damage severity (stiffness reduction~75%) or when the damage location is highly influent (states 8 and 9, i.e., with both base columns damaged). Otherwise, Figure 10a shows a high performance level obtained in correspondence of a small noise reduction, achievable with improvement in measurement accuracy or by applying data cleansing and noise reduction techniques. All these findings are consistent with what preliminary assessed in the previous works [28] for the simpler cantilevered beam model. script elapsed on average 2.2 s to perform a single training fold, 3.9 s for CDF validated estimation and 0.02 s for the damage detection task on a single damage state, with no major differences among the 20 states, on an Intel ® Core™ i7-10750H CPU with 2.6 GHz base frequency and 15.8 GB of available RAM. Similar elapsed times were found for the other numerical and experimental cases.

Results for the Frame FE Model
The settings for the second numerical investigation are the same as described in Table  1. Figure 9 reports the results of the frame structure FE model for an increasing level of artificially added noise. Again, as for the previous case studies, consider that no damage was inserted in the baseline (state 1). For what concerns the noise level, the results in Figure 10b show how, for a realistic noisy scenario consistent with current measurement technology, a good performance is reached in the states with high damage severity (stiffness reduction ~75%) or when the damage location is highly influent (states 8 and 9, i.e., with both base columns damaged). Otherwise, Figure 10a shows a high performance level obtained in correspondence of a small noise reduction, achievable with improvement in measurement accuracy or by applying data cleansing and noise reduction techniques. All these findings are consistent with what preliminary assessed in the previous works [28] for the simpler cantilevered beam model.

Experimental Validation
The experimental validation was performed on a well-known database realised at the Los Alamos National Laboratory (LANL). As for the frame numerical model described in the previous Section, the structure (represented in Figure 11) is a three-storey shear-type frame. Seventeen scenarios were considered, as described in Table 4. They include one or more stiffness reductions applied to the beams (to simulate damage) and/or one or more added masses (to simulate changing operational conditions). The dataset also included some scenarios with a bumper-column apparatus attached at the top floor. The original intent of this mechanism was to emulate the nonlinear effects induced by a breathing crack [41], which is known to act as a pointwise source of bi-linearity (a deeper discussion can be found in Reference [42]).

Experimental Validation
The experimental validation was performed on a well-known database realised at the Los Alamos National Laboratory (LANL). As for the frame numerical model described in the previous Section, the structure (represented in Figure 11) is a three-storey sheartype frame. Seventeen scenarios were considered, as described in Table 4. They include one or more stiffness reductions applied to the beams (to simulate damage) and/or one or more added masses (to simulate changing operational conditions). The dataset also included some scenarios with a bumper-column apparatus attached at the top floor. The original intent of this mechanism was to emulate the nonlinear effects induced by a breathing crack [41], which is known to act as a pointwise source of bi-linearity (a deeper discussion can be found in Reference [42]).

Case Description
Linear States 1 Undamaged Baseline 2 Added mass of 1.2 kg at the base 3 Added mass of 1.2 kg at the 1st floor 4 87.5% stiffness reduction in one column of the 1st inter-storey 5 87.5% stiffness reduction in two columns of the 1st inter-storey 6 87.5% stiffness reduction in one column of the 2nd inter-storey 7 87.5% stiffness reduction in two columns of the 2nd inter-storey 8 87.5% stiffness reduction in one column of the 3rd inter-storey 9 87.5% stiffness reduction in two columns of the 3rd inter-storey Nonlinear States 10 Distance between bumper and column tip 0.20 mm 11 Distance between bumper and column tip 0.15 mm 12 Distance between bumper and column tip 0.13 mm 13 Distance between bumper and column tip 0.10 mm 14 Distance between bumper and column tip 0.05 mm 15 Bumper 0.20 mm from column tip, 1.2 kg added at the base 16 Bumper 0.20 mm from column tip, 1.2 kg added on the 1st floor 17 Bumper 0.10 mm from column tip, 1.2 kg added on the 1st floor Generally, the presence of such nonlinearities is used to detect and localise surface cracks in otherwise linear systems [43]. However, these data are here used with another intent. As it can be noticed from Figure 12c,d, the insertion of the mechanism causes an  87.5% stiffness reduction in one column of the 1st inter-storey 5 87.5% stiffness reduction in two columns of the 1st inter-storey 6 87.5% stiffness reduction in one column of the 2nd inter-storey 7 87.5% stiffness reduction in two columns of the 2nd inter-storey 8 87.5% stiffness reduction in one column of the 3rd inter-storey 9 87.5% stiffness reduction in two columns of the 3rd inter-storey Nonlinear States 10 Distance between bumper and column tip 0.20 mm 11 Distance between bumper and column tip 0.15 mm 12 Distance between bumper and column tip 0.13 mm 13 Distance between bumper and column tip 0.10 mm 14 Distance between bumper and column tip 0.05 mm 15 Bumper 0.20 mm from column tip, 1.2 kg added at the base 16 Bumper 0.20 mm from column tip, 1.2 kg added on the 1st floor 17 Bumper 0.10 mm from column tip, 1.2 kg added on the 1st floor Generally, the presence of such nonlinearities is used to detect and localise surface cracks in otherwise linear systems [43]. However, these data are here used with another intent. As it can be noticed from Figure 12c,d, the insertion of the mechanism causes an increase in the system stiffness. Therefore, the natural frequencies grow as well, while a breathing crack would decrease them. However, the bumper-column mechanism affects the mode shapes making them detectable as a deviation from the undamaged baseline. Moreover, since the structural nonlinearities generate noise-like distortions in the frequency response, the damage scenarios 10-17 can be utilised as well to prove the robustness of the procedure when dealing with distorted signals. increase in the system stiffness. Therefore, the natural frequencies grow as well, while a breathing crack would decrease them. However, the bumper-column mechanism affects the mode shapes making them detectable as a deviation from the undamaged baseline. Moreover, since the structural nonlinearities generate noise-like distortions in the frequency response, the damage scenarios 10-17 can be utilised as well to prove the robustness of the procedure when dealing with distorted signals. The acquisition procedure can be summed up as follows (more details can be found in the original source [44]). Fifty band-limited  white Gaussian noise realisations were applied as input for 25.6 s at the structure base. The system response was recorded at four points corresponding to the three levels plus the base, with a sampling frequency = 320 Hz. The mode shapes were then extracted from the frequency response functions (FRFs) defined between the acceleration output time histories and the force input utilising the FRVF procedure described previously.
Importantly, for states 10-17, the concept of "mode shapes" in a nonlinear context could be not totally accurate. It is necessary to remark that the specific experimental setup (single-input multi-output acquisitions with a random noise as input) only allows linear system identification [45]. The experimental identification of nonlinear normal modes would require more complex procedures which are beyond the objective of this study. The acquisition procedure can be summed up as follows (more details can be found in the original source [44]). Fifty band-limited  white Gaussian noise realisations were applied as input for 25.6 s at the structure base. The system response was recorded at four points corresponding to the three levels plus the base, with a sampling frequency f s = 320 Hz. The mode shapes were then extracted from the frequency response functions (FRFs) defined between the acceleration output time histories and the force input utilising the FRVF procedure described previously.
Importantly, for states 10-17, the concept of "mode shapes" in a nonlinear context could be not totally accurate. It is necessary to remark that the specific experimental setup (single-input multi-output acquisitions with a random noise as input) only allows linear system identification [45]. The experimental identification of nonlinear normal modes would require more complex procedures which are beyond the objective of this study. Therefore, the extracted mode shapes should be considered as the ones of the underlying linear system for these eight nonlinear states [45].
The Data Augmentation procedure described in Section 2 was applied to increase the number of training and validation data. For illustrative purposes, Figure 13 shows the resulting model obtained by applying the GP Regression over the three mode shapes identified from the experimental data. The results of the experimental validation are shown in Figure 14. As for all the previous examples, no damage was to be identified in the first state, therefore there are no false alarms nor successful damage identifications. The EFT-based procedure correctly identifies all the structural changes except for the configurations with the bumper-column gap larger than 10 mm and no added mass (i.e., states 10 to 12). However, as it can be seen from Figure 12c, the nonlinear distortions are minimal for this larger gap and the response of the structure is almost indistinguishable from the pristine baseline. The algorithm struggled to recognise the last two states (16 and 17) as well when fed with the first mode shape; this issue did not arise either with the second or the third eigenmode. On the other hand, due to the larger confidence intervals (Figure 13c), the fitting of the third mode shape returned a relatively larger number of false positives (almost constantly 17% for any damage scenario).

EFT vs. EVT
To conclude this work, a comparison between the extreme value and the extreme function theory has been performed on all the numerical and experimental data. This comparison can be done directly since the same spatial coordinates were utilised to define all

EFT vs. EVT
To conclude this work, a comparison between the extreme value and the extreme function theory has been performed on all the numerical and experimental data. This comparison can be done directly since the same spatial coordinates were utilised to define all the damaged and undamaged mode shapes. That means that any i-th mode shape is defined on the same n i points. For EFT, this results in n vectors, which return a single distribution over functions, i.e., 1-dimensional data; while in EVT, there are n i distributions over scalar values, each one with n samples.
The results are reported in Appendix A (Tables A1-A5) for all the numerical and the experimental case studies presented in the previous Sections (in the same order). Everywhere, the increment is calculated as FALSE POSITIVE EVT-FALSE POSITIVE EFT, expressed in percentage.
As expected, the incidence of false positives decreased significantly, with an improvement everywhere larger than 6% for the HAR spar Finite Element Model (even with a double-digit increment in the false alarm reduction for the case with SNR = 40 dB) and still very relevant for the second numerical case study. The EVT outperformed the EFT in terms of fewer false positives only for the simply supported beam and for the third mode shape with the lowest SNR. This is most probably due to the combination of the larger confidence intervals of the GP regression over the 3rd mode shape, the larger variability induced by the noise, and the structural symmetry. The experimental data confirmed the key findings encountered on the simulated dataset. Specifically, the percentage of false positives decreased for all the three mode shapes.

Conclusions
This study investigated the validity of EFT as a framework for mode shape-based SHM in 2-dimensional plate-like and frame structures, which are representative of common applications in Aerospace and Civil Engineering, respectively. The rationale is that damage is a rare occurrence and therefore the data-driven models defined over the undamaged conditions must take this consideration into account. Experimental and numerical data were utilised to this aim for different damage locations and severity levels. The robustness of the procedure to artificially added white Gaussian noise has been methodically studied and the results have been benchmarked against the results of the same algorithm trained with pointwise EVT values, showing a statistically relevant decrease of the number of false positives. Moreover, the use of EFT allows to compare points where the data were not directly collected; furthermore, it assigns only one possible outcome ("normal" or "abnormal") to the whole function. On the other hand, the single components of the mode shape, if taken one by one as in the EVT, could be under the threshold at some modal coordinates and over it at other locations for the same identification.
This study provides a strong foundation for future works in the field of EFT and EVT for structural health monitoring purposes. It will be important to validate the proposed approach in situ on real-life civil structures. Another related field of research, which the Authors are intended to investigate deeper in the next future, regards the scarcity of data from damaged structures and the need to compensate with numerically simulated data.