A new damage detector for bridges based on natural frequencies with missing data

Abstract

Automatic structural health monitoring can simplify the surveillance process of many structures and bridges if its underlying methods return correct interpretations of the structural state. A common method to differentiate between a damaged and undamaged state of a structure is to use its modal properties from an assumed undamaged state to build a baseline to which all new information is compared. The comparison can be performed by calculating the Mahalanobis squared distance (MSD) of natural frequencies. Considering the inherent uncertainties associated with automatic system identification, a new novelty detection algorithm is proposed in this work, intended to work with missing and randomly available natural frequency information, like the outcome of automatic operational modal analysis and mode tracking algorithms. The moments of a multivariate normal distribution used to characterize the bridge’s undamaged behavior are determined elementwise. The damage indicator measures the MSD of new data points to this distribution considering the available natural frequencies and normalizes it using the chi-squared nature of the MSD. The proposed method works as intended for two numerical cases with 25% of the natural frequency values missing at random, where all but the smallest of damages become clearly detectable. It is also tested on two real-world bridges, one of which has a small, controlled change to its structural state. The automatic operational modal analysis of the bridges’ data recordings leads to randomly missing natural frequency values. Despite this, the damage can be detected by the proposed novelty detection algorithm.

Introduction

Implementation of structural health monitoring (SHM) has commonly acknowledged positive societal and environmental impacts. ¹ Bridges are a prime target for SHM due to their strategic importance and difficult replacement in case of closure. Manual SHM techniques are labor and cost-intensive and are therefore currently not suitable to be scaled for widespread implementation. Automatic SHM methods solve this problem, but they are still an active field of research. By identifying the characteristic values of a bridge only dependent on its structural system and tracking their evolution over time, any change in these values could be highlighted as a change to the structural system and, therefore, the potential onset of damage.

Two main hurdles remain in the workflow of automatic SHM: robustly identifying characteristic values and separating normal and abnormal variations in these characteristics. The Mahalanobis squared distance (MSD) shows promise at solving the second barrier, regardless of which characteristic values are chosen. In 2011, Kiremidjian et al. ² and in 2012, Mosavi et al. ³ show two separate approaches of using autoregressive models fitted to the vibration data of bridges to identify the onset of imposed damage. Also, using autoregressive coefficients of vibration data, Svendsen et al. ⁴ show that it is possible to detect different damage types on truss railway bridges using an MSD-based novelty indicator. An advantage of auto-regressive (AR) coefficients is the ease of extracting them from recorded acceleration datasets. They are also sensitive to nonlinear behavior of the structure. More information can be found in Refs. 4–6. A second option for characteristic values of a structure is to use its modal properties, such as natural frequencies or mode shapes. Although natural frequencies are considered less sensitive to damage than AR coefficients, the advantage of this approach is that natural modes are well known and understood by researchers and partitioning engineers alike. They are also easy to establish, both numerically and experimentally, requiring few sensors. ⁷ However, correctly identifying a structure’s modal properties from its vibration recordings can be challenging without manual intervention. Automatic operational modal analysis (AOMA) algorithms are a class of algorithms capable of automatically performing system identification, Since around 2010, a growing selection of these algorithms has become available and more performant. Two notable AOMA algorithms for either automaticity and performance are those by Reynders et al. ⁸ and Kvåle and Øiseth. ⁹ Once the mode properties have been identified, a novelty detection (ND) algorithm is necessary to separate between environmental and operational variability in mode properties and damage induced changes.

Several recent ND algorithms use the natural frequencies of a structure as input for an MSD-based indicator to identify how far the latest available characteristic values are from an established baseline describing the normal behavior of the structure. Their research focuses on describing the structure’s normal behavior as precisely as possible. Multiple multivariate normal distributions, known as a Gaussian mixture model (GMM), are used to describe various normal structural states (e.g., warm and freezing weather) in Figueiredo et al. ¹⁰ and Santos et al. ¹¹ Replacing GMMs with agglomerative concentric hyperspheres is successfully demonstrated to differentiate between multiple normal operational regimes of the Z24 bridge and its damaged regime. ¹² Neural networks also show promise in learning how to identify a threshold between damaged and undamaged conditions based on an MSD value from the natural frequencies of a structure and information about the environmental conditions.^13,14 Locally robust MSD values, using concepts from local metric learning, considering the operational and environmental variability, can accurately identify damage in real-world cases.^15,16

All the methods above, except for Entezami et al., ¹⁶ have in common that every dataset contains a complete set of all characteristic values of interest. This is rarely the case when AOMA is used to extract modal parameters of a structure, even when using the best-performing algorithms. ¹⁷ Most structural modes cannot be identified from all data recordings, leading to a missing data problem, which makes the definition of a baseline normal condition of the structure ambiguous. When addressing missing data for ND, several approaches have been explored. The missing data values are in all cases imputed; missing values are substituted by estimated values. Ma et al.^18,19 leveraged two variations of principal component analysis (PCA), probabilistic PCA and mixture of probabilistic PCA, to impute missing data and to subtract the environmental and operational variability of strain measurements for ND of a sports arena with retractable roof. Natural frequencies with missing entries are also imputed and used by Xu et al.^20,21 for ND under environmental variations by the utilization of a robust PCA or a low-rank matrix approximation combined with cointegration analysis. Imputation can be rather unreliable when dealing with a significant amount of missing data, ²¹ and it has the potential to introduce bias even when the missing data rate is low.

This work proposes a ND method capable of dealing with missing and randomly available—no a priori pattern can characterize when the data will be missing—natural frequency data by defining the moments of a multivariate normal distribution in an elementwise approach. The MSD of a new measurement to this distribution is determined using the available modal parameters, and a normalized novelty probability is determined considering the chi-squared distribution that the MSD follows. The method is developed in close relationship to the AOMA and mode tracking (MT) algorithms by which all data recordings are processed before they are given to the ND algorithm. The proposed algorithm’s functionality and performance are discussed and tested through two numerical and two real-world cases, namely the Hardanger bridge and the Rognes bridge.

Theory

The ND presented in this work relies on three previous steps that pretreat the measurement data and extract the relevant information for ND. A brief overview of each of these steps, as well as the methods used in this research, is presented in this section. The first step, operational modal analysis (OMA), is shortly discussed and illustrated with the covariance-driven stochastic subspace identification (SSI) (cov-SSI) algorithm. Next, the principles for automatic interpretation of the outcome of OMA are highlighted, and the algorithm for this purpose by Dederichs and Øiseth ²² is summarized. The third step, MT, is illustrated with the MT algorithm presented in Dederichs and Øiseth. ²² Finally, the main idea for using the MSD as a damage detector is discussed before this work’s proposed ND method is presented.

Automatic operational modal analysis

Operational modal analysis

OMA determines the eigenmodes of a structure whose excitation is unknown from measured vibration accelerations at multiple points on the structure. Many methods exist to perform OMA in the time and frequency domains. OMA and its various methods are a well-known topic and will not be discussed in detail in this work; the interested reader is referred to Rainieri and Fabbrocino ²³ for more information. The time-domain cov-SSI system identification method has become the most popular method, due the quality of its outcome and its relatively low computational time, and is also the method used in this work.

Automation algorithms

The number of modes present in a data recording is, in most cases, unknown beforehand. Solving the system identification eigenproblem multiple times for different model orders helps solve this problem but introduces an element of manual decision into the process. As the mode detections appear in complex conjugate pairs, there will be at most half the model order $o$ mode detections. When increasing the model order, more mode detections will be made; however, it is unknown if a mode detection represents a physical mode or merely a mathematical artifact of the system identification procedure. To separate these two types of poles, one looks at the similarity of the poles detected at different model orders. Poles that occur persistently at many model orders are likely physical poles, and those with no similar occurrences at different model orders are mathematical poles. Visually, when creating a graph of model orders versus pole frequency, physical poles form vertical lines and mathematical ones scatter randomly throughout the plot. Such a plot is known as a stabilization diagram. The traditionally manual task of stabilization interpretation has been looking through stabilization diagrams to extract which detected modes are present and their properties while disregarding the mathematical poles. This task is time and resource-consuming and quickly becomes infeasible if many datasets need to be analyzed. For this reason, many AOMA algorithms that can perform this task automatically and without human supervision have been suggested since the early 2010s.

An AOMA algorithm generally tends to be composed of three main parts. First, there is a screening of all the poles resulting from the OMA algorithm to remove as many mathematical poles as possible. Second, similar poles representing the same structural mode are grouped and separated from those representing other modes. Third, the resulting groups are screened again to remove groups of mathematical poles that were not eliminated in the first step. Finally, representative properties from the remaining groups are selected to characterize the detected structural modes. Usually, these steps are performed using unsupervised machine learning algorithms, logical decisions, or a combination. The main steps of the Dederichs 2023 ²² AOMA algorithm are discussed to illustrate the functionality of an AOMA algorithm.

The Dederichs 2023 algorithm starts by separating the physical poles from the mathematical ones using a two-component GMM on the stability properties of the poles. A $k$ -component GMM is an unsupervised fuzzy classifier algorithm where every element is given a likelihood of belonging to each of the $k$ clusters. ²⁴ If this method is used as a discrete classifier, elements are attributed to the cluster with the highest likelihood. The pole stability properties are measured in relative difference in frequency $\mathrm{df}$ , mode phase collinearity (MPC) $d\mathrm{MPC}$ , and mode shape $d\mathrm{MAC}$ separately. The difference is measured between a pole and its nearest neighbor at the next inferior model order. The poles with low relative differences tend to be physical, as they form a vertical line in the stabilization diagram, and those with larger relative differences tend to be mathematical because they do not have similar poles at the next model order. The GMM separates the poles based on these three measures, and the poles deemed in the larger-difference cluster are disregarded from the next steps.

The next step of the algorithm combines the similar poles into groups by structural mode using single-linkage hierarchical clustering. Hierarchical clustering is an iterative process where all elements start in separate groups, and the two closest groups are combined until only a selected number of groups remain. ²⁵ The distance metric used by the AOMA algorithm between two poles $i$ and $j$ is:

$$\begin{array}{c}{d}_{i,j}=d{f}_{i,j}+1-\mathrm{MA}{\mathrm{C}}_{i,j}\hfill \end{array}$$

(1)

where $\mathrm{MA}{\mathrm{C}}_{i,j}$ is the modal assurance criterion between the mode shapes ${\varphi }_i$ and ${\varphi }_j$ of poles $i$ and $j$ , and $d{f}_{i,j}=\frac{|f_i-f_j|}{max(f_i,f_j)}$ . The iterative process is stopped once the number of clusters equals the average model order of the OMA. As the average model order of the OMA will overestimate the number of structural modes, some structural modes will be represented by multiple hierarchical clusters, which need to be recombined. Similar hierarchical clusters are merged using their average features and similarity limits extracted from Step 1.

Differentiating between the small hierarchical clusters containing leftover mathematical poles and those representing structural modes is achieved with a second two-component GMM. The inputs to the classifier are the size of each hierarchical cluster and the number of clusters of similar size. The structural modes tend to be represented by large hierarchical clusters of disparate sizes, whereas the mathematical noise clusters tend to be small and homogenous. Once the two types of hierarchical clusters are identified, the structural mode features are extracted as the mean values of the poles associated with each cluster.

Mode tracking

AOMA algorithms extract the structural modes from a single dataset, giving an insight into the structural modes at a discrete time. When multiple datasets from the same structure exist, AOMA can be run on each dataset to gain insight into the structural modes’ evolution over time. However, due to the randomness of unsupervised learning algorithms and the different excitation conditions of each dataset, the outcome from the AOMA may differ considerably. For example, the first fundamental frequency from dataset one may not be detected in the next dataset, or the labeling order of the detected structural modes may not be the same in one AOMA outcome as in the next. These difficulties are often compounded by complex structures with closely spaced modes and modes whose features vary depending on loading and environmental conditions. Furthermore, it is also important to follow the evolution of a mode through large changes, such as variations due to environmental effects or the sudden onset of damage, with either a sudden or slow onset, because this change may give valuable information about the state of the structure. A MT algorithm can identify which structural modes are worth tracking and label the AOMA detections from subsequent datasets. The functionality of the MT algorithm presented in Dederichs and Øiseth ²² is summarized here.

The MT algorithm has two stages: an initialization stage, where the structural modes worth tracking (the reference modes) are identified, and a second online stage, where the AOMA detections from newly available datasets are attributed to a reference mode. The first stage requires the user to define three parameters, which are the number of datasets attributed to identifying the reference modes $N_i$ , a limit of similarity between two mode shapes $\mathrm{ma}{\mathrm{c}}_{\lim }\in [0,1]$ , and a “keep” parameter $k\in [0,1]$ . Using hierarchical clustering to parse through all the mode detections in the first $N_i$ datasets, the algorithm groups all mode detections which are similar and likely represent a common structural mode. The similarity is measured in terms of the inverse MAC number between two modes, where an inverse MAC of 0 indicates that the modes are identical, and an inverse MAC of 1 means that two modes are as different as they can get. Any groups with more than $k·N_i$ mode detections are then considered the reference modes of the structure. Using the latest $m\in \mathbb{N}$ available datasets, a statistical representation of each reference mode is calculated, containing the average frequency, the frequency standard deviation, the average mode shape, and the average MAC number between all pairs of mode detections.

The second stage of the algorithm is the online MT phase. For each new dataset which has been processed by AOMA, every mode detection within the new dataset (with frequency $f_r$ and mode shape ${\varphi }_r$ ) is attributed to the nearest reference mode as measured by Equation (1). The attribution is only confirmed if, first (Equation (2)), the frequency of the mode detection falls inside the 98% confidence interval of the reference mode, assuming a normal distribution of the frequency, and if, second (Equation (3)), the average MAC value between the mode detection and all the mode shapes associated with the reference mode falls in the 98% confidence interval of the reference mode, assuming an exponential distribution with the average MAC number of the reference mode as parameter:

$$\begin{array}{c}{\mathrm{\Phi }}^{-1}(0.01|{\mu }_f^{(i)},{\sigma }_f^{(i)})\le f_r\le {\mathrm{\Phi }}^{-1}(0.99|{\mu }_f^{(i)},{\sigma }_f^{(i)})\hfill \end{array}$$

(2)

$$\begin{array}{c}\overline{\mathrm{mac}}({\varphi }_r)\le \mathrm{EX}{\mathrm{P}}^{-1}(0.98|{\lambda }^{(i)})\hfill \end{array}$$

(3)

${\mu }_f^{(i)}$ , ${\sigma }_f^{(i)}$ , and ${\lambda }^{(i)}$ are the representative mean natural frequency, natural frequency standard deviation, and average MAC number of reference mode $i$ , defined by $M^{(i)}$ the set composed of the latest $m$ datasets attributed to the reference mode. ${\mathrm{\Phi }}^{-1}$ and $\mathrm{EX}{\mathrm{P}}^{-1}$ are the inverse cumulative density functions of a normal distribution and an exponential distribution:

$$\begin{array}{c}{\lambda }^{(i)}=\frac{1}{m^2}\sum _{a\in M^{(i)}}\sum _{b\in M^{(i)}}(1-\mathrm{mac}({\varphi }_a,{\varphi }_b))\hfill \end{array}$$

(4)

$$\begin{array}{c}\overline{\mathrm{mac}}({\varphi }_r)=\frac{1}{m}\sum _{a\in M^{(i)}}(1-\mathrm{MA}{\mathrm{C}}_{r,a})\hfill \end{array}$$

(5)

If the attribution is confirmed, the reference mode parameters are recomputed to consider the newest $m$ available datasets. This process is repeated for all mode detections in a new dataset until no new datasets are available. The MT algorithm makes a provision in case a change of modal parameters is too abrupt to be tracked, and the mode trace is lost. At regular intervals $I_r$ , if a mode trace has been deemed to be lost—no new mode attributions happened in the last $m$ datasets—the first stage of the algorithm is run to search for reference modes again. If an updated reference mode is found to be similar to the lost one (MAC of average mode shapes within $\mathrm{ma}{\mathrm{c}}_{\lim }$ ), the lost reference mode is considered redetected, and its parameters are updated to match the updated reference mode. Mode shapes are robust to damage and do not tend to change much if it occurs. This makes the mode shapes ideal to track the evolution of the reference modes when large changes occur. A MAC value between undamaged and damaged conditions below $\mathrm{ma}{\mathrm{c}}_{\lim }$ (recommended to be between 0.75 and 0.90) would no longer be tracked, which is a downside of this tracking approach. However, a damage of that severity and urgency would undoubtedly also be reported via other channels, and likely quicker than identifiable by any modal-based surveillance algorithm. All other reference modes whose trace was not lost are not updated in this procedure. Typical values for the parameters are: $N_i=100$ , $\mathrm{ma}{\mathrm{c}}_{\lim }=0.85$ , $k=0.15$ , $m=25$ , and $I_r=200$ . The MT algorithm is not intended to define the environmental and operational variability (EOV) of the structure but merely to attribute mode detections from the AOMA stage to the correct structural reference mode. For this reason, $N_i$ can be considerably lower than the number of datasets typically required to define the EOV. The MT algorithm only needs to approximately define each mode’s properties to be able to follow their evolutions in the online stage of the algorithm. The EOV can then be defined by another algorithm (of the operator’s choice) once enough datasets have been tracked.

The presented MT algorithm is very flexible and can track structural modes even if they have large and rapid changes to their properties. It can also track modes only irregularly and sporadically detected by AOMA.

Novelty detection

Mahalanobis squared distance

The MSD measure is commonly used in damage detection algorithms in civil engineering. It assumes that its features—in this case, frequencies—are described by a multivariate normal distribution. It measures how far a point is from the mean of the multivariate distribution, considering the variances and correlations within the distribution. For a multivariate normal distribution $D$ characterized by its mean vector $\overline{x}$ and its covariance matrix $\Sigma$ , the MSD, denoted $d^2$ , of a point $p$ to $D$ is given by

$$\begin{array}{c}{d}^2={(p-\overline{x})}^T{\mathrm{\Sigma }}^{-1}(p-\overline{x})\hfill \end{array}$$

(6)

A similar metric is encountered when defining the reliability of a structure using a probability of nonexceedance of a set $\beta$ threshold. The variables affecting the structural reliability are transformed to standard normal space, and the likelihood of an event is the distance from the center of the distribution to the point. Indeed, it is shown ²⁶ that if the underlying distributions in a reliability problem are normal, then the problem is equivalent to the MSD. By taking the variable substitution $y=U^{T}x$ , where $U$ is the matrix of eigenvectors from the eigenvalue decomposition of the covariance matrix $\Sigma =\mathrm{US}{U}^T$ , the mean $\overline{y}$ and covariance matrix ${\Sigma }_y$ of the substitution variable can be written as

$$\begin{array}{c}\overline{y}=U^T\overline{x}\hfill \end{array}$$

(7)

$$\begin{array}{c}{\mathrm{\Sigma }}_y=U^T\Sigma U=S\hfill \end{array}$$

(8)

Equation (8) shows that the covariance matrix of the substitution variable is diagonal, as it equals the eigenvalue matrix $S$ , which, by definition, is diagonal. Taking the inverse transform $\overline{x}=U\overline{y}$ , and $p=\mathrm{Uq}$ , and substituting them into the Mahalanobis distance (Equation (6)):

$$\begin{array}{c}{d}^2={(q-\overline{y})}^T{U}^{T}U{\mathrm{\Sigma }}_{y}U{U}^T(q-\overline{y})=\sum _{i=1}^n\frac{1}{{\sigma }_i^2}{(q_i-{\overline{y}}_i)}^2\hfill \end{array}$$

(9)

it becomes clear that the MSD reduces to the sum of $n$ independent standard normal variables, where $n$ is the number of dimensions of the feature space. The reliability problem can be seen as a generalization of the MSD, as it can have any underlying variable distributions, which are then transformed to independent standard normal space using, for example, a Rosenblatt transform.

As the MSD is a sum of independent standard normal variables, it follows a chi-squared ${\chi }^2(n)$ distribution with $n$ degrees of freedom and its expected value can be easily found:

$$\begin{array}{c}E\left[d^2\right]=E\left[\sum _{i=1}^n\frac{1}{{\sigma }_i^2}{(q_i-{\overline{y}}_i)}^2\right]=\sum _{i=1}^{n}E\left[\frac{1}{{\sigma }_i^2}{(q_i-{\overline{y}}_i)}^2\right]=n\hfill \end{array}$$

(10)

because the expected value of a standard normal variable squared is equal to its variance, which, in this case, is one. The ${\chi }^2(n)$ distributed nature of the distance measure allows a single novelty probability to be given to any point as to how much of an outlier it is from the multivariate Gaussian distribution. The novelty probability is found by computing the cumulative density function $F(d^2,n)$ of the ${\chi }^2(n)$ distribution (Equation (11)) for the MSD:

$$\begin{array}{c}{F}_{\chi }(d^2,n)=\frac{\gamma (\frac{n}{2},\frac{d^2}{2})}{\mathrm{\Gamma }(\frac{n}{2})}\hfill \end{array}$$

(11)

Here, $\gamma (a,b)$ is the lower incomplete gamma function. Combining the ${\chi }^2(n)$ cumulative density function and the expected value of the MSD, one finds that the expected probability of a point $p$ randomly drawn from the underlying multivariate normal distribution can be expressed as

$$\begin{array}{c}{F}_{\chi }(E[d^2],n)=\frac{\gamma (\frac{n}{2},\frac{n}{2})}{\mathrm{\Gamma }(\frac{n}{2})}\hfill \end{array}$$

(12)

The expected probability only depends on the number of dimensions of the multivariate distribution. It is illustrated in Figure 1 for dimensions ranging from 1 to 10.

Figure 1.

Expected ${\chi }^2(n)$ probability of a random point $p$ drawn from a multivariate normal distribution, measured by the MSD $d^2$ .

Proposed ND method

The proposed novelty method is based on the natural frequencies of a system and their evolution in time, as tracked by an MT algorithm. It is similar to many proven methods, such as ²⁷ but its main focus is dealing with incomplete and missing data from complicated structures. As discussed, AOMA algorithms are affected by the excitation of the structure, the system identification method, and the automatic stabilization diagram interpreters, resulting in less-than-perfect results. They commit errors, such as false and duplicate detections, and they do not detect all structural modes from all datasets, whereas the errors of AOMA algorithms can be ignored and compensated by an MT algorithm, and the latter cannot replace missing data. Furthermore, an MT algorithm is also prone to errors and may fail to track parts of structural modes. These errors lead to mode traces which are not complete; the structural modes do not have natural frequency information for every dataset. The proposed ND algorithm aims to valorize the natural frequency information as much as possible and does not disregard any structural modes due to poor detection rates. The initial $N_b$ -tracked datasets are used to establish the features’ baseline. A single multivariate normal distribution is used to describe the natural variations of the structural modes. Each tracked structural mode $i$ has a set $D^{(i)}$ containing the dataset numbers of all datasets containing a detection of the mode. The list of features and the natural frequencies of each structural mode $i$ is written as

$$\begin{array}{c}{f}^{(i)}=\{f_k^{(i)}\},k\in D^{(i)}\hfill \end{array}$$

(13)

a vector containing $n^{(i)}$ frequencies. The mean natural frequency of each structural mode, ${\mu }_f^{(i)}$ , is computed as the mean of all the detection instances of the frequency.

$$\begin{array}{l}{\mu }_f^{(i)}=E(f^{(i)}),k\in D^{(i)},k\le N_b\\ =\frac{1}{n_{Nb}^{(i)}}\sum _{\begin{array}{c}k\in D^{(i)}\\ k\le N_b\end{array}}{f}_k^{(i)}\end{array}$$

(14)

where $n_{Nb}^{(i)}=\mathrm{len}(\{f_k^{(i)}\},k\in D^{(i)},k\le N_b)$ is the number of datasets containing a realization of structural mode $i$ within the initial $N_b$ datasets. Although the mean natural frequencies are not necessarily calculated from the same selection of datasets, this approach maximizes the information available from each structural mode. In a similar approach, the covariance matrix of the multivariate distribution is also defined in an elementwise fashion, maximizing the information available. Each element ${\sigma }_f^{(\mathrm{ij})}$ of the covariance matrix ${\Sigma }_f$ corresponding to the covariance between structural modes $i$ and $j$ is computed using all datasets containing mode detections for both structural modes:

$${\sigma }_f^{(ij)}=\frac{1}{n_{Nb}^{(ij)}-1}\sum _{\begin{array}{c}k\in D^{(i)}\stackrel{D^{(j)}}{\cap }\\ k\le N_b\end{array}}(f_k^{(i)}-{\mu }_f^{(i)})(f_k^{(j)}-{\mu }_f^{(j)})$$

(15)

where $n_{Nb}^{(ij)}=\mathrm{len}(\{f_k^{(i)}\},k\in D^{(i)}\stackrel{D^{(j)}}{\cap },k\le N_b)$ is the number of datasets containing realizations of structural modes $i$ and $j$ within the initial datasets. Establishing the covariance matrix ${\Sigma }_f$ requires that each pair of structural modes is present in at least two datasets. If this is not the case, fewer structural modes should be considered, or more datasets should be included in the baseline. There should be significantly more than two datasets available per pair of structural modes to obtain a covariance representative of the underlying multivariate distribution.

Once the properties of the assumed features’ distribution— ${\mu }_f=\left\{{\mu }_f^{(i)}\right\},\forall i$ and ${\Sigma }_f$ —have been established using the $N_b$ baseline datasets, the novelty probability of subsequent new datasets can be evaluated. Each new dataset $m>N_b$ is represented by the vector of natural frequencies it contains $f_m=\left\{f_m^{(i)}\right\},i\in M$ , where $M$ is the set of structural modes detected in the dataset. The distance that this new dataset is from the features’ distribution is found using the MSD:

$$d_m^2={(f_m-{\mu }_f^{(M)})}^T{\Sigma }_f^{(M)}{(f_m-{\mu }_f^{(M)})}^T$$

(16)

considering that ${\mu }_f^{(M)}$ is the reduced mean frequency vector containing only the natural frequencies corresponding to the structural modes $M$ present in dataset $m$ . Likewise, ${\Sigma }_f^{(M)}$ is the reduced covariance matrix containing only the rows and the columns of ${\Sigma }_f$ corresponding to the structural modes $M$ .

The novelty probability of the new dataset $P(m)$ is found from the cumulative density function of the chi-squared distribution ${\chi }^2(n_m)$ with $n_m$ degrees of freedom, where $n_m$ is the number of structural modes in $M$ :

$$\begin{array}{c}P(m)=F(d_m^2,n_m)=\frac{\gamma (\frac{n_m}{2},\frac{d_m^2}{2})}{\mathrm{\Gamma }(\frac{n_m}{2})}\hfill \end{array}$$

(17)

Figure 2 graphically illustrates the flow of the proposed ND method.

Figure 2.

Flowchart of the proposed novelty detection method.

Alternate ND method

An alternate ND method is also proposed, where the computation elements are identical to the method presented in “Proposed novelty detection method,” except that the baseline is not fixed but evolves as new datasets become available. Once the initial baseline containing the first $N_b$ datasets has been established and the novelty probability $P(m)$ of the next dataset $m=N_b+1$ has been determined, the number of baseline datasets is increased by one $N_{b,\mathrm{new}}=N_{b,\mathrm{old}}+1$ . The baseline features ${\mu }_f$ and ${\Sigma }_f$ are then recomputed considering the first $N_{b,\mathrm{new}}$ datasets. Figure 3 illustrates the flow of the alternate ND method.

Figure 3.

Flowchart of the alternate novelty detection method with an updating baseline.

The updating baseline in this method allows for more information to be considered when computing the novelty probability of a new dataset. As long as the structure behaves according to the same underlying multivariate distribution, the updating baseline will better capture the features’ distribution and return a more accurate novelty probability. At the first novelty occurrence—a change to the structure’s underlying distribution—both methods should be able to detect this, and the novelty probability of subsequent datasets should tend toward one. If the structure stays in this new configuration, the first method will continue, indicating a novelty probability tending toward one, whereas the second method will see a steady decrease in the novelty probability as the new state is slowly introduced into the baseline. The first method is advantageous if the base state of the structure is known beforehand, and the main purpose is to capture any deviation from this state. The base state need not necessarily be an undamaged state. It could also be a structure already in a damaged state, and the intention of the monitoring is to highlight any, potentially dangerous, evolution of the damage. The second method can be advantageous if one is interested in highlighting states where the structure has not yet been. A totally new state, different from all previous states, should be flagged and could be interpreted as a prompt for a more in-depth inspection of the structure. However, the more states are included in an updating baseline, the less sensitive this method will become to new conditions, as these new conditions may be too similar to other previous states to trigger a high novelty probability. If the primary concern of the structure monitoring is to highlight the first onset of damage, both methods are similar. But if the onset of damage happens very slowly, the second method may have time to adapt and learn the changes as they occur, thus never leading to a clear damage detection. For this reason, if the first onset of damage is of concern, the first method—the standard formulation of the ND algorithm—is recommended.

Ill-posed covariance matrices

Ill-defined covariance matrices can occur when creating a covariance matrix using the technique shown in this work, as each element in the covariance matrix is estimated from a different subset of the population. The covariance matrix may not be positive definite, which cannot be used to measure a distance, as needed for the MSD. A nonpositive definite covariance matrix does not necessarily mean the analysis is impossible. In real-world cases with missing data, most datasets will not contain frequency values for all the reference modes, requiring a set of rows and columns to be selected from the nonpositive definite covariance matrix to create a smaller matrix, which can still be positive definite, hence useful. For all cases where the subselection of the covariance matrix is nonpositive definite, the ND algorithm does not consider the dataset at which this occurs.

Including more datasets in the training of the ND algorithm will tend to decrease the likelihood of an ill-posed covariance matrix problem because more data will lead to better estimates of each element in the covariance matrix and decrease the size of any potential bias in their estimation. Likewise, in the presence of increased levels of missing data, more datasets will be necessary in the training of the ND algorithm to obtain a positive definite covariance matrix. Each feature pair must be present in at least two datasets for a covariance value to be established between those features.

Covariance matrix regularization, where a weighted average is performed between the estimated covariance matrix and a positive definite diagonal matrix to ensure that the outcome is also positive definite, is not applicable in this case, because it degrades and dilutes the covariance information. Another approach used to fix nonpositive definite covariance matrices, ²⁸ developed for risk management in a financial setting, only returns a semipositive definite matrix, which cannot be inverted. Finally, Pavez and Ortega ²⁹ propose an unbiased estimator of the covariance matrix for data with missing observations. For a large number of observations, the outcome of this method coincides with the outcome of the method proposed in this work.

Numerical case studies

Four scenarios for demonstrating the ND algorithm are presented in this work. The first two cases are numerical, and the last two are based on experimental data recorded from two suspension bridges, the Hardanger bridge and the Rognes bridge, both located in Norway. The first numerical case is based on example natural modes and is intended to demonstrate the features of the ND algorithm. The second numerical example is based on the natural modes of a numerical simply supported beam bridge subject to the abrupt onset of damage. The damage is simulated by a reduction in the height of the beam cross-section for a segment of the bridge. Multiple combinations of damage length and reduction in height are considered to demonstrate the sensitivity of the ND algorithm. The numerical cases are based on the simulation of the OMA outcome, not the generation of vibration time series of the structure. Simulating the outcome of OMA or automatic OMA reduces the complexity and calculation time for the simulations compared to generating vibration time series of the structures subjected to random white noise and extracting the natural modes using system identification. As 1000 datasets are considered per simulation case, employing AOMA to interpret generated vibration dataset would be prohibitive time and complexitywise; therefore, the outcome of AOMA is simulated for each dataset, greatly speeding up the process. The simulated outcome of AOMA includes the variability normally observed from AOMA on real data. Each analysis case has a list of reference structural modes described by a natural frequency and mode shape. At every simulated dataset, a selection of these reference modes is included in the dataset. The selection is decided on a random, per-mode basis; each mode has a probability of 0.75 of appearing in a dataset. This value represents the average performance of a good AOMA algorithm when dealing with real data. ¹⁷ A Bernoulli distribution describes the number of mode detections par dataset. The likelihood of the number of mode detections per dataset, assuming 9 or 10 structural modes, is shown in Table 1, when an individual mode has a probability of 0.75 of appearing.

Table 1.

Likelihood of the number of modes generated in a dataset, assuming 9 or 10 structural modes, and an individual probability of 0.75 of appearance of each mode.

Structural modes	Number of modes in dataset
	0	1	2	3	4	5	6	7	8	9	10
9	0.00	0.00	0.00	0.01	0.04	0.12	0.23	0.30	0.23	0.08	—
10	0.00	0.00	0.00	0.00	0.02	0.06	0.15	0.25	0.28	0.19	0.06

Random noise is then added to each mode to simulate the uncertainties associated with environmental and system processes. The noise is randomly chosen from two normal zero-mean distributions, one for frequency noise with standard deviation ${\sigma }_f$ and one for each element of the mode shape ${\sigma }_{\varphi }$ . The random noise is added to the natural frequency and to each term of the mode shape. After that, the mode shape is renormalized so that its maximal absolute value equals one. Specificities concerning the value of the noise is provided in each case study.

AOMA algorithms also tend to make mistakes when processing data returning duplicate detections (two or more detections of the same structural mode) and false detections (indicating a structural mode which is not one). These inconsistencies are also simulated in the generated datasets by including a random number of duplicate and false detections drawn from two exponential distributions with parameters ${\lambda }_d=0.7$ for the duplicate detections and ${\lambda }_f=1.1$ for the false detections. The numbers of duplicate and false detections for each simulated dataset are drawn from the exponential distributions, as shown in Table 2. A duplicate detection is then simulated in the same way as a correct mode detection. A false detection (which tends to be very different in terms of natural frequency and mode shape from structural modes) is simulated by adding normal, zero-mean random noise with a standard deviation 10 times as large as for the correct modes, making it likely very different from the original reference mode. The duplicate and false detections are created for each dataset from a random reference mode. Damage can be introduced at any dataset by altering the reference modes from which the datasets are simulated.

Table 2.

Occurrence probabilities of false and duplicate detections per dataset.

		Number of false or duplicate detections per dataset
Type	Parameter $\lambda$	0	1	2	3	4
False detections	1.1	0.67	0.22	0.07	0.02	0.01
Duplicate detections	0.7	0.50	0.25	0.13	0.06	0.03

Once a run of 1000 datasets is simulated, the MT algorithm described in section “Mode tracking” is used to identify the structural modes and their evolution throughout the 1000 datasets. The mode detection algorithm has been shown to perform very well in cases like these, ²² achieving above 95% of correct assimilations for each mode. Errors made by the MT algorithm are only related to not identifying a mode detection when it was present. Misidentifications of mode detections are extremely rare and are negligible. Although using an MT algorithm is not entirely necessary when working with simulated datasets because the ground truth of what has been simulated is inherently known, its use is representative of the intended use of the ND algorithm on real data, where the ground truth of the outcome AOMA is not known. The numerical cases are intended to be as close to real-world applications of the ND algorithm while remaining rapid enough to compute to provide an in-depth overview of its performance.

First numerical case study

The first numerical case study is used to show that, for simple cases, the ND algorithm performs as intended. To this end, four applications are considered: first, an example without damage or change to the natural frequencies; second, another example without damage, but where the natural frequencies have a cyclic behavior; third, an example with one significant damage; and finally, two significant and equal damages. Both formulations of the ND algorithm (the fixed training baseline datasets, and the increasing training baseline datasets) are discussed in relation to these applications.

The reference natural modes and their parameters for the first numerical case are presented in Table 3. The natural modes do not correspond, in natural frequency or mode shape, to any real structure. They are chosen to contain a selection of well-identifiable modes and a set of closely spaced modes, as is common in real structures. The functioning principles of the ND algorithm remain valid, and its expected outcomes are the same as if these reference modes were to be drawn from a numerical simulation of a structure.

Table 3.

First numerical case reference natural modes.

Mode	Naturalfrequency [Hz]	Mode shape		Naturalfrequency [Hz]	Mode shape
Mode 1	1.0		Mode 6	3.5
Mode 2	1.2		Mode 7	3.6
Mode 3	2.0		Mode 8	3.65
Mode 4	2.4		Mode 9	4.8
Mode 5	3.0		Mode 10	4.95

Four different applications are considered to illustrate the intended functions of the ND algorithm. These are first, no changes (no introduction of damage); second, sinusoidal cyclic behavior of all reference modes with ±2.5% variation in natural frequency and a period of 24 datasets, but no abrupt onset of damage on top of these variations; third, an abrupt increase in natural frequency of 2% of all modes as of dataset 750; finally, two abrupt increases in frequency of 1.5% each as of datasets 700 and 900. Twenty simulations of 1000 datasets are created for each type of novelty.

The first case (no damage) is intended to illustrate a baseline case where no changes to the modal parameters of a structure are expected. The second novelty case (cyclic behavior) simplifies a structure affected by temperature changes. A day and night cycle is represented by 24 datasets, imitating an hourly measurement system. The ND algorithm is not expected to, and should not, highlight any damage for these two cases. The third novelty case (single abrupt change) is the first case with an introduction of damage, represented by the abrupt change in the natural frequencies. The ND algorithm should react to this type of change and highlight it. With two abrupt introductions of damage offset in time, the final case is introduced to show the difference between the standard ND algorithm and its alternate formulation with a dynamic baseline. The standard ND algorithm should only be able to highlight a single change, whereas the alternate formulation should react to both changes.

All simulated mode detections, from every application case, simulation, and dataset, have random normal noise added to them to simulate the variability and uncertainty associated with system identification and AOMA. The standard deviations of the noise are ${\sigma }_f=0.02$  Hz and ${\sigma }_{\varphi }=0.06$ .

Figure 4 shows the evolution of the mode natural frequencies over the 1000 datasets (left column) and the results from the standard ND algorithm (right column) for all four applications (in descending order) for one of the 20 simulations per novelty case. For all applications, all ten reference modes are detected and tracked without problem by the MT algorithm. The first 100 datasets are used to initiate the MT. The following 500 datasets (datasets 100–600) are used to set up the initial training baseline for the ND algorithm. As intended, for the two first applications, the ND algorithm does not highlight any changes during its testing phase (datasets 600–1000), as shown by the bold blue line indicating the moving average of the novelty probability based on the latest ten datasets. The faint blue line indicating the novelty probability of each dataset is unsteady and has large variations created by the noisy underlying data. The dashed and bold green line indicated the expected novelty probability (as found by Equation (12)) for 10 natural frequencies. As predicted, the moving average of the novelty probability oscillates around the expected novelty probability, indicating that no damage is identified. As the ND algorithm used in Figure 4 works with a fixed training baseline, both abrupt damage applications look similar, with the moving average outlier probability rapidly moving above the 0.99 dashed red probability threshold after introducing the damage, as is expected. The second damage introduction in the fourth case goes unnoticed in the novelty probability, as it is indistinguishable from the first change. When looking at the MSD values, it is clear that there are two distinct changes, but both altered states result in a novelty probability so close to 1 that they become impossible to differentiate. Finally, the number of frequencies considered in each dataset $n_{\mathrm{freq}}$ can be seen to be random, as intended, but hovering around 7, which is expected with ten natural frequencies in total and when combining the original 0.75 presence rate of the reference modes and the 0.95 tracking rate of the MT algorithm. In the abrupt damage applications, $n_{\mathrm{freq}}$ can be seen to drop right after the introduction of the changes before recovering again. This is due to the MT algorithm needing to adapt to the new structural state and making fewer mode classifications during this period.

Figure 4.

Tracked natural modes (left column, (.i)) and outcome of the novelty detection algorithm (right column, (.ii)), for the cases (a) no novelty, (b) cyclic behavior, (c) one abrupt change, (d) and two abrupt changes.

Figure 5 shows the novelty probability moving average from all 20 simulations of each application, using the fixed training baseline (left column) and the increasing training baseline (the alternate formulation of the ND algorithm—right column). The no-change and cyclic behavior cases show, as expected, no onset of novelty. The novelty probability moving averages all tend to oscillate around the expected probabilities for both ND algorithm formulations. The single abrupt damage application leads to all novelty probability moving averages quickly tending toward one once the damage is introduced. In the fixed baseline formulation, it stays just below one, whereas in the alternate formulation, it slowly decreases for datasets further after introducing the change and tends back toward the expected probability. This behavior is expected because the training baseline contains more and more of the datasets with change, which normalizes the change, making it seem like normal behavior. The advantage of this alternate formulation is seen in the application with two abrupt damages. The fixed baseline ND algorithm cannot differentiate the two changes, whereas the alternate formulation tends toward one again when the second change is introduced (dataset 900). This alternate formulation does have the downside of decreased sensitivity because when the baseline is expanded, it will require larger changes to be distinct from the new baseline due to its now larger inherent variability. This sensitivity decrease is shown by the more rapid decrease in the novelty probability moving average after introducing the second change compared to the first change. The novelty probability moving average stays above 0.99 for about 40 datasets after the first change but only for 20 after the second change. These four applications illustrate that the ND algorithm can identify damages when they are present and is not affected by recurrent and natural changes in the mode frequencies.

Figure 5.

Outlier probability moving average for all 20 simulations: (a) no novelty, (b) cyclic behavior, (c) single abrupt change,(d) two abrupt changes: (i) fixed training baseline (ii) increasing training baseline (alternate formulation).

Second numerical case study

The second numerical case is developed to show the ND algorithm’s capabilities and limitations in detecting a physically representative damage imposed on a beam bridge. Twelve different damage sizes, represented by concrete cracking on a varying section length and depth on the underside of the bridge deck, are simulated to test the algorithm. Only the standard formulation of the ND algorithm is considered in this case study, as only a single damage is imposed.

This case considers natural mode properties extracted from a numerical model of a simply supported beam bridge. This model allows for physically representative changes to be introduced to the mode properties, from which a sensitivity analysis of the ND algorithm can be performed. The 20 m long bridge has a massive section with a 4 m width and 0.5 m height. It is assumed to be constructed from concrete, and due to the small displacements of any vibrations, only the concrete contributes to the stiffness of the bridge. The concrete is assumed to be lightly cracked and therefore $E=1.5·{10}^{10}$  N/m² and $\nu =0.3$ . The bridge is discretized into 40 equal elements of 0.5 m in length. The elements are Bernoulli beams, and the torsion constant is found using the Saint–Venant approach for massive sections. Although this approach yields 41 nodes with 6 degrees of freedom each, the mode shapes extracted from the bridge are reduced to 9 nodes with two degrees of freedom (displacement in lateral and vertical directions) to better correspond to the outcome observed from a monitoring system. Five of the nine nodes are located on one side of the bridge and the other four on the other. The bridge’s setup and the nodes’ location from which the mode shapes are extracted are shown in Appendix A (Figure A1). The nine lowest-frequency global modes are extracted from the model to form the reference modes of the simulations. Six modes have a predominantly vertical behavior, two have a predominant lateral movement, and one follows a torsional deflection. Their properties and mode shapes are reported in Appendix A (Table A1).

Small alterations to the same model are made to extract reference damaged states. Twelve damage cases of cracked concrete on a varying section of the underside of the bridge are considered. The cracked concrete section is assumed to no longer contribute to the stiffness of the structure. The 12 cases result from the combinations of four different affected section heights and three different affected area lengths located close to the midspan of the bridge. The damaged section heights are 0.005, 0.010, 0.025, 0.050 m, representing 1%–10% of the initial section height. The affected areas of the bridge are 9.5–10.0 m, 9.5–10.5 m, 9.0–11.0 m, representing 1.25%–5% of the length of the bridge.

As for the first numerical case, the system identification step is bypassed to avoid the associated computation time and complexity. The natural frequency and mode shape of each generated mode are corrupted with random normal noise. Three noise corruption cases are considered to simulate high, medium, and low-quality data, and the corruption representative size is shown in Table 4. The valuation of data does not refer to measurement quality alone but encompasses all sources of error, from excitation of the structure to interpretation by an AOMA algorithm. These noise levels correspond to experimentally observed variations in the performance of AOMA algorithms.

Table 4.

Representative noise size for the high, medium, and low-quality simulated datasets.

Noise	High quality	Med. quality	Low quality
${\sigma }_f$ [Hz]	0.007	0.014	0.03
${\sigma }_{\varphi }$ [—]	0.02	0.04	0.08

All normal and expected variations encompassed in the EOV, for example, temperature, humidity, and live load-driven variations, of the structure are not described by this added noise. It is challenging to realistically describe and simulate the scale, time evolution, and mode correlation impacts of EOV in a numerical example. Each of these elements has a large possible value range and can be the subject of discussion. As the focus of this work is not to propose a ND method rivaling the state of the art in terms of differentiation between EOV and damage, but to showcase a simple damage detection method capable of dealing with missing natural frequency information, the simulations are simplified by not including EOV other than the normal random noise.

An overview of the changes in the frequency of all structural mode for each damage case and each data quality is presented in Appendix A (Table A2). The changes are also described as a factor of the standard deviation of the added noise to emphasize their relative importance, even though relative frequency changes are small. The changes resulting in a damage superior to four times the standard deviation of the noise are highlighted in gray in the table. This size of damage is clearly differentiable (e.g., zoomed graph of Figure 6(b), where the difference is 5.5 times ${\sigma }_f$ for the light green mode) for an observer looking at a MT plot, given the correct zoom level.

Figure 6.

Mode tracking results for three selected cases of damage and noise level: (a) largest damage L090-110 H0050 with low noise, (b) medium damage L095-105 H0025 with medium noise, and (c) smallest damage L095-100 H0005 with high noise. The second row shows a zoom around dataset 700 for the two highest frequency modes.

By design, the natural mode frequencies in this case study will tend to be uncorrelated; therefore, even small changes in natural frequencies contribute individually to increasing the MSD and the novelty probability. In a case where EOV leads to some degree of correlation between natural frequencies, a large change in one frequency may not contribute much to increasing the novelty probability if the other frequencies change according to their correlations identified in the training baseline. The mode shapes are only slightly affected by the damage cases, and because the ND algorithm does not use mode shapes to determine outliers, their change is not reported.

The persistence of all modes is set constant at 0.75 for all simulations. This means that all modes have a uniform random probability of presence in a dataset of 0.75. On average, over the 1000 datasets in each simulation, each structural mode should be present 750 times and the distribution of the number of mode detections per dataset is reported in Table 1.

As for the first numerical case, 20 simulations of 1000 datasets are created for each damage state. The first 100 datasets are always used to initiate the MT, and the 500 subsequent datasets are used as the baseline to train the ND algorithm. Only the fixed baseline variant of the ND algorithm is considered in this analysis because both variants are equally good at identifying the first onset of novelty, which is the task at hand for this second numerical analysis. The remaining 400 datasets are used for testing the ND algorithm. A total of 100 of these testing datasets are undamaged, and 300 contain damage, as all changes are introduced at dataset 700. The MT is performed as described in section “Mode tracking.” In total, 36 cases (12 damage states and 3 noise parameters) are created with 20 runs of 1000 datasets.

The outcome of the MT for three separate cases is shown in Figure 6. First (a) high-quality data and the maximal damage level L090-110 H0050; second, (b) medium-quality data and medium damage L095-105 H0025; and finally, (c) low-quality data and the smallest damage L095-100 H0005. A zoom of the MT illustrates the differences between the three cases. These illustrate what should be the hardest case for the ND algorithm, the medium case, and the easiest case. Overall, the MT can be seen to perform very well; in all cases, it tracks the nine modes from dataset 100 until the last dataset. Only when a large novelty onset is introduced (zoom of case (a)) does the MT algorithm temporarily lose track of some of the modes until a redetection of the lost modes is performed at dataset 800. In other cases (light green mode of case (b)), the MT algorithm temporarily loses track of the mode after the change is introduced but does not require a redetection to find the mode again, resulting in fewer datasets lost. As the MT algorithm performs as intended in all cases, it will never permanently limit the information available to the ND algorithm. The most difficult case for the MT algorithm, should become the easiest case for the ND algorithm, because the damage is most pronounced in that case. Similarly, the easiest case for the MT algorithm, with the smallest damage, should become the most challenging case for the ND algorithm, as the changes are very poorly differentiable. Indeed, visually, it is not possible to identify any change to the natural frequencies for that case (Figure 6(c)) due to the small scale of damage impact.

The ND algorithm process is shown in Figure 7 for the same three cases as the MT. The first case leads to a distinct increase in the moving average novelty probability and stagnates just below 1 as of dataset 750. The second case also tends toward a novelty probability of one, but not as rapidly and distinctly as the first case. The change is more clearly visible in the moving average than in the novelty probability for each dataset. The final case does not tend toward one, looking like a case with no change, as seen in the first numerical case study. Furthermore, the novelty probability is very jagged, as is expected with poor-quality data. The number of frequencies present per dataset shows when the damage introduced was significant enough to require the MT algorithm to redetect certain modes. Before any change is introduced, the number of frequencies per dataset hovers around 6, as expected, but in case (a), $n_{\mathrm{freq}}$ decreases to around three between datasets 700 and 800, indicating that many modes are no longer tracked. After the redetection is performed at dataset 800, it returns to around 6. A very slight decrease can be observed for case (b), but it does not require waiting until dataset 800 to return to normal. For case (c), the damage is so small that a decrease in detected frequencies cannot be remarked.

Figure 7.

Novelty detection process for cases: (a) largest damage L090-110 H0050 with low noise, (b) medium damage L095-105 H0025 with medium noise, and (c) smallest damage L095-100 H0005 with high noise. The second row shows a zoom around dataset 700 for the two highest frequency modes.

Tables 5–7 show, for each of the 36 analysis cases, the share of simulations which led to the detection of novelty and the average dataset number when this detection was made. ND is achieved when the moving average cumulative probability exceeds 0.99, identified by the dashed red line in Figure 7. The dataset number corresponding to the ND is the earliest dataset at which this happens. Only the first exceedance of a cumulative probability of 0.99 is considered.

Table 5.

Novelty detection rate (share of the 20 simulations to return a novelty detection) and average onset for the low measurement noise case. (High quality data).

	L095_100		L095_105		L090_110
	Detection rate	Averagedetectiondataset	Detectionrate	Averagedetectiondataset	Detectionrate	Averagedetectiondataset
H0005	0	—	0	—	1.0	745.9
H0010	0	—	1.0	754.2	1.0	719.8
H0025	1.0	730	1.0	723.7	1.0	723.2
H0050	1.0	770.8	1.0	722	1.0	711.5

Table 6.

Novelty detection rate (share of the 20 simulations to return a novelty detection) and average onset for the medium measurement noise case. (Medium quality data).

	L095_100		L095_105		L090_110
	Detectionrate	Averagedetectiondataset	Detection rate	Averagedetectiondataset	Detectionrate	Averagedetectiondataset
H0005	0	—	0	—	0	—
H0010	0	—	0	—	1.0	767.0
H0025	0.85	737.8	1.0	728.3	1.0	716.8
H0050	1.0	734.2	1.0	735.1	1.0	724.3

Table 7.

Novelty detection rate (share of the 20 simulations to return a novelty detection) and average onset for the high measurement noise case. (Low quality data).

	L095_100		L095_105		L090_110
	Detectionrate	Averagedetectiondataset	Detectionrate	Averagedetectiondataset	Detection rate	Averagedetectiondataset
H0005	0	—	0	—	0	—
H0010	0	—	0	—	0	—
H0025	0	—	0.2	754.3	0.9	742.1
H0050	0.65	742.5	1	722.8	1	716.5

The three smallest damages (L095-100 H0005, L095-105 H0005, and L095-100 H0010) never return a ND with a probability exceeding 0.99, regardless of the data quality. As the damage amplitude increases, it becomes gradually detectable for decreasing data quality. For example, the H0005 height reduction is only detectable with the highest quality data. For the lowest quality data, both the H0005 and H0010 height reductions are undetectable. Increasing damage amplitude usually leads to an earlier detection point because more mode combinations clearly indicate novelty. However, some of the largest damages take slightly longer to detect than smaller counterparts. This is due to the MT algorithm not immediately being able to track the structural modes with the largest changes after the onset. The ND is then left with only the modes with the smallest changes from which to identify novelty, leading to a less rapid detection. Larger damages requiring longer to be detected is a drawback in a damage detection method. This behavior is entirely due to the MT algorithm—the ND algorithm also introduces a detection lag while the moving average of the novelty probability increases, but this lag is independent of damage size and defined by the length of the moving average period (in this case, 10 datasets). For small changes in modal properties, the MT algorithm will continuously detect the reference modes (Figure 6(c.ii)), introducing no lag. For larger changes, the MT algorithm will lag before realizing a change has taken place (Figure 6(b.ii) light green mode); the larger the change, the longer that lag tends to be. However, once the change’s amplitude requires a redetection of the reference modes (Figure 6(a.ii)), the maximal detection lag is bound by the parameters selected for the MT algorithm (assuming a successful redetection). The maximal lag will not exceed $I_r+m$ , where the selected refresh interval $I_r$ is combined with $m$ the memory. No detections must occur for $m$ datasets for the trace to be deemed lost, and then, in the worst-case scenario when the latest refresh interval has just occurred, a further $I_r$ will be required before the next redetection and the discovery of the damage. Assuming a constant probability of damage occurring at a given dataset, the expected leg in the detection of large damage, requiring a reduction of the mode traces, due to the MT algorithm is $P_{\mathrm{LAG}}(t)=\frac{I_r}{2}+m$ . Although this delay in detection is not ideal, the existence of a maximal bound, controllable by the user, partially mitigates this downside.

Further insight into the ND process can be obtained by looking at the receiver operating characteristics (ROC) curves from each simulated case, shown in Figure 8. ROC curves are established to evaluate the performance of the ND algorithm. They represent the relative trade-offs between true positives, or the probability of detection, and false positives, or the probability of a false alarm. Considering the moving average novelty probability as the damage detector, the ND algorithm’s performance at distinguishing between undamaged and damaged states for each damage case and data quality can be assessed. The faint lines are the ROC curves of each simulation, and the bold lines are the average ROC curves for each case. The larger the damage and the better the data quality, the more precisely the ND algorithm can distinguish between cases, as shown by the increasing areas under the ROC curves. Data quality significantly impacts medium damage states because high- and medium-quality data lead to clear distinguishability (L090-110 H0005, L095-105 H0010, and L095-100 H 0025), whereas low-quality data does not. Only the smallest damage amplitude does not lead to a ROC curve indicating good distinguishability between undamaged and damaged states for the highest-quality data. When comparing to the ND in Table 5, one realizes that the damage cases L095-100 H0010 and L095-105 H0005 do not have any damage detections but have a ROC curve indicating good distinguishability between damage states, meaning that the probability threshold of 0.99 is too high to identify damage for those cases, but an unambiguous detection should be possible with a lower threshold. The medium- and low-quality data also have damage states where no damage is detected with a probability above 0.99, but where the ROC curves indicate that a clear detection should be possible for a lower threshold.

Figure 8.

ROC curves. Faint lines represent every randomly generated case. Bold lines represent the average curve for all cases with identical parameters.

The second numerical case provides an insight into which size of damage is identifiable by the proposed ND algorithm in the presence of missing and randomly available natural frequency data. The impact of the AOMA outcome is highlighted, as it affects the minimum size of damage detectable. The MT algorithm is intrinsically linked to the performance of the ND algorithm, as the outcome of the former is the sole input into the latter. This dependence is especially visible in cases where the damage is large, and the MT algorithm requires time to regain the tracking of structural modes after its onset, which leads to added time before damage is detected. In some cases, this temporary hindrance of the ND algorithm, and associated delay in damage detection may be critical. To mitigate this risk, the parameters of the MT algorithm should be chosen according to the design of the measurement intervals to avoid excessive detection delays.

The case is limited by the difficulty in adequately and representatively being able to simulate environmental and operational variations in the natural frequencies. Only the expected variability of AOMA algorithms is included in the simulations. The damages presented lead to small relative differences in frequency, but for some modes, especially under larger damages, this difference is multiple times the standard deviation of the variations introduced to imitate the AOMA outcome. The natural frequencies are, however, uncorrelated; therefore, any change in frequency, even small, to any mode makes the dataset an outlier from the assumed multivariate normal distribution created by the training baseline, and a sign of damage. In real-world cases, where EOV can be larger than the expected damage, some correlation between the natural frequencies is expected. This correlation can allows a mode to have a relatively large change in natural frequency without it raising the suspicion of damage if (1) enough datasets have been included in the training baseline for the behavior and correlation of the natural frequencies to be quantified for the full range of the EOV; (2) the other modes also undergo a change, as expected by the correlation of the natural frequencies. On the other hand, assuming that point (1) is met, but point (2) is not, even a small change in natural frequency of a mode can lead to a damage detection, because the dataset becomes an outlier. The ND algorithm is tested on two experimental cases in the following section to show that, even in the presence EOV, it reaches the correct decision on the state of the structure.

Experimental example

Hardanger bridge

The Hardanger bridge (Figure 9) is a 1310 m long suspension bridge crossing the Hardanger fjord in South-Western Norway. Since 2013, it has been equipped with 20 triaxial CUSP-3D series strong motion accelerometers with a ±4 g measurement range and nine sonic WindMaster Pro 3D anemometers. Data are recorded for 30 min if any of the anemometers measure a wind speed above a set threshold. The data are recorded locally before being uploaded to a server at the Norwegian University of Science and Technology (NTNU). The sensor layout can be found in Appendix B (Figure B1 and Table B1) and in Fenerci et al., ³⁰ and the recorded datasets are openly available. ³¹ There are 13 well-documented^32–34 (illustrated in Appendix B (Table B2)) natural modes below 0.425 Hz, which has been tracked using the MT algorithm discussed in section “Mode tracking.” A total of 1100 datasets covering 2013 and 2014 are considered in this work. During that period, according to the bridge owners who have carried out regular manual inspection of the bridge, there are no known changes to the structural system of the bridge. Therefore, no detection of novelty is expected. Nonetheless, using experimental data on the ND algorithm provides valuable insight into its functionality and output. The data recordings from the Hardanger bridge are processed with the fully automatic OMA algorithm presented in the Theory section. The MT results are shown on the left side of Figure 10. The first 100 datasets are used to set up the reference modes, and the following 500 datasets are used as the training baseline for the ND algorithm. Twelve of the 13 modes are tracked, making the assumed multivariate normal distribution approximated by the training phase of the ND algorithm twelve-dimensional. The number of frequencies $n_{\mathrm{freq}}$ present over all the datasets is quite random, but averages around 10, and only rarely dips below 4, showing the good performance of the AOMA and MT algorithms. The distribution of the number of detected natural frequencies over all the datasets is shown in Figure 11. Only 13.6% of the training and testing datasets has detections for all structural modes; most often, a dataset will contain 10 or 11 natural frequencies.

Figure 9.

Hardanger suspension bridge.

Figure 10.

Mode tracking results (left) and novelty detection process (right) for the Hardanger bridge.

Figure 11.

Distribution of the number of detected natural frequencies per dataset using AOMA for the Hardanger bridge. Right of the dashed gray line are the datasets where all modes are detected.

The outcome of the ND algorithm does not, as expected, detect any novelty, despite the changes in natural frequencies being superior to those tested in the second numerical case study. The EOV of the bridge is properly characterized by the 500 training datasets, and all the variations in a mode’s natural frequency in the testing phase are met by the expected change in the natural frequencies of the other modes. Only very few datasets lead to a novelty probability above 0.99, and the 10-dataset moving average probability oscillates between 0.25 and 0.75, which is the range of values that is expected. Overall, the ND algorithm performs as expected for a real-world case, where the recorded data are pretreated using AOMA and MT algorithms.

Rognes suspension bridge

The Rognes suspension bridge (Figure 12) is a small 80 m suspension bridge carrying a single traffic lane located in Rognes in central Norway. The bridge was equipped from December 2022 to May 2023 with nine S-LYNKS wireless three-dimensional MEMS accelerometers that regularly recorded 20-min datasets. The sensor layout can be found in Appendix B (Figure B2, Table B3, and Table B4). These were then uploaded to the S-LYNKS cloud platform via 4G network. The raw acceleration data were then retrieved and stored on a server at NTNU. Each data recording was split into ten 4-min recordings overlapping 50% of the previous 4-min segment. In total, 569 datasets were processed using the AOMA algorithm of section “Automatic operational modal analysis.” The nine structural modes with lowest natural frequency are of interest for this study; they are illustrated in Appendix B (Table B4). For a series of eight 20-recordings in April 2023, a 3350 kg vehicle was parked at midspan on the bridge to increase its weight by approximately 2% to create a controlled change in the response of the bridge. During the measurement campaign, the weather conditions varied widely, ranging from −25°C with up to 75 cm of snow and +10°C without snow. The bridge had to be always kept open to vehicles, implying snow clearing and the associated piling of compact snow along both sides of the bridge, increasing the loading on the bridge and the transient excitation generated by vehicles crossing the bridge. Additionally, Figure 13 shows the distribution of the number of detected modes per dataset. Nearly no datasets contain a detection of all nine modes. For these reasons, the mode traces shown in Figure 14 are extremely irregular and undergo large variations. Nonetheless, given the circumstances, the MT algorithm performs a decent job of tracking the modes.

Figure 12.

Rognes suspension bridge with the test vehicle parked at midspan.

Figure 13.

Distribution of the number of detected natural frequencies per dataset using AOMA for the Rognes bridge. Right of the dashed gray line are the datasets where all modes are detected.

Figure 14.

Mode tracking and novelty detection process of the Rognes bridge. The gray shaded area represents the datasets where the vehicle was parked on the bridge.

Fifty datasets were used to detect the reference modes, and the following 300 datasets form the training baseline of the ND algorithm. However, the relatively low number of training datasets combined with the large variations in mode natural frequencies leads to 48% of the 219 remaining datasets’ covariance matrices being nonpositive definite and hence unusable. The moving average novelty probability of the valid datasets, as seen in Figure 14, crucially only reaches above 0.99 probability threshold when the vehicle is parked on the bridge (gray-shaded area). Before the vehicle is moved onto the bridge and after it is removed, the moving average novelty probability varies between 0.35 and 0.9. The novelty probability of the unaffected datasets may not tend toward the expected probability as closely as for the Hardanger bridge, but given the variability of the natural frequencies beforehand, this cannot be expected. Notably, the distinction between the two structural states can only be achieved when both cold and snowy and mild, and no snow conditions are considered in the training baseline of the ND algorithm. When this is the case, the ND algorithm can highlight the small change in structural behavior, even in these challenging conditions with a large EOV and many missing frequency datapoints.

Summary

This work proposes a new ND algorithm based on natural frequencies and the MSD intended to work with missing and randomly available natural frequencies. AOMA algorithms are not perfect and do not always identify every structural mode from every dataset. This ND algorithm uses this partial data as part of its assessment of a structure. The main findings in this work are as follows:

• Two formulations of the ND algorithm. The standard ND algorithm has a fixed training baseline, and the alternate ND algorithm considers all datasets, except the latest available one, as part of the baseline. The standard formulation is best suited when the identification of the first onset of damage is critical. The alternate formulation can highlight multiple damages occurrences but is not sensitive to slow changes and its sensitivity decreases after each new structural change. Both ND algorithm formulations are shown to work as intended.

• The first numerical study uses simplistic cases to show that the ND algorithms work as intended for controlled cases. When no damage, or only cyclic environmental and operational variations are simulated, they do not detect damage. When damage is introduced abruptly, both formulations detect damage.

• In the second numerical study, using a simulation of a simply supported bridge subjected to concrete cracking damage on its underside, and varying the size of this damage, the standard formulation of the ND algorithm is shown to be sensitive to all but the smallest damage, when, on average, only six of the nine structural modes are present in each dataset. The second numerical study is limited by not including the effects EOV other than that linked to AOMA algorithms. Its main conclusions remain valid if the assumption is made that EOV leads to a degree of correlation between natural frequencies and that damage to a structure induces changes which are not expected by this correlation.

• The Hardanger bridge provides a real-world case without any known damage. The standard formulation of the ND algorithm does not highlight any damage on this bridge, showing that it does not misinterpret the EOV in the natural frequencies.

• The Rognes bridge, to which mass was added to simulate damage, provides a challenging setting for the ND algorithm because of the large EOV, relatively small damage, and many missing natural frequency datapoints. The ND algorithm manages to detect the imposed damage, albeit in a less clear fashion than in the numerical studies.

• A large and diverse training baseline is necessary to ensure the best results. First, including many datasets in the baseline improves the estimation of the covariance matrix for the MSD, resulting in fewer problems with ill-posed selections of the covariance matrix. Second, the full range of EOV of a structure needs to be included in the baseline to avoid misclassifying normal variations as damage.

The proposed method has two major limitations:

• Defining the covariance matrix in an elementwise fashion from different subpopulations can, in some cases, lead to a nonpositive definite covariance matrix which cannot be used to measure a distance using the MSD.

• The MSD is measured to a single distribution of features, meaning that the natural variations in features need to be described by a multivariate Gaussian distribution, which may not always hold true.

Conclusion

A new ND algorithm is proposed based on the MSD and uses the natural frequencies of a structure as features. The ND algorithm is intended to work with missing and randomly available natural frequencies to maximize the available information, because AOMA algorithms do not always provide frequency estimates for every structural mode of interest. The performance of the ND algorithm is not expected to match state-of-the-art damage detectors, as its focus does not lie therein. Nonetheless, its performance is shown to be more than satisfactory on two numerical cases, and two experimental bridges. Due to the complexity in adequately and representatively imitating EOV, the numerical case used to test the ND algorithm foregoes simulating EOV. The precise performance quantification of the algorithm from this study are to be considered with this in mind. However, it remains that damage is detectable even for very small changes in natural frequency, if these exceed the intrinsic variations of the natural frequency. Future work should focus on developing a framework in which to simulate the outcome of AOMA algorithms combined with the EOV of a structure, in order to provide a more realistic numerical study of ND algorithms’ performances. When EOV is introduced in the experimental cases, damage remains detectable even if the EOV is large, there are many missing frequency datapoints, and the damage is small. Good estimations of the correlation of natural frequencies due to the full range of EOV are key for damage much smaller than the EOV to be detected.