On the Consistent Classification and Treatment of Uncertainties in Structural Health Monitoring Applications

Many critical engineering tasks such as model calibration, nonlinearity identification, and damage detection of complex dynamical systems are significantly affected by the presence of uncertainty. Computational tools developed within fields that are driven from dynamics, such as structural dynamics and control, system identification, and structural health monitoring (SHM), do rely on similar input–output system representations, but often employ different definitions and classifications of the various uncertainty sources, ultimately affecting how the uncertainty is quantified and communicated for supporting decision-making.

The purpose of this review work is twofold. First, to define and classify the sources of uncertainty within the context of structural dynamics, system identification, and SHM applications through the prism of the classical input–output system representation. Second, we aim to present an overview of methods that are suitable for treating the various forms of encountered uncertainties from the perspective of defining downstream tasks (e.g., prediction, inference, and damage identification) within these domains.

Figure 1 shows a high-level representation of a monitored dynamical system, as often represented in the context of dynamics, control, and SHM [1]. The box S represents a physical engineering system that is of dynamic nature and is being monitored. This system is set up to address a particular downstream task. Therefore, it is designed to capture only some aspects of the real world, often linked to the states or degrees-of-freedom (DOFs) of the dynamical system, and typically comprises an assembly of components that interact with each other. The system receives an input u(t), while concurrently being exposed to disturbances d(t) due to confounding processes, namely, environmental and operational variability (EOV) [2], as well as to deterioration processes h(t). When exposed to ${u(t),d(t),h(t)}$⁠, the system produces an output y(t), a subset of which is usually monitored via the use of sensors. Input, disturbance, and deterioration quantities can also be measurable directly or indirectly via the use of dedicated sensors. A model $M$ is employed that aims to represent the actual physical system as closely as possible. Such a model can be prescribed a priori, e.g., numerical approximation of the system on the basis of available engineering knowledge, or it may alternatively be inferred from data through the process of system identification [3]. The fidelity and complexity of this model are dependent on the task, on the availability of informative measurements, and on domain knowledge and physics insights. In its simplest representation, the system model is assumed to be linear time-invariant (LTI); an assumption which is often inadequate, due to the varying nature of the loads and the configuration of the system itself, which can be nonstationary or nonlinear [4]. It becomes evident that uncertainties are involved in every component of the essential system depicted in Fig. 1. In this work, we attempt to distinguish and classify these sources in a way that allows us to subsequently guide the selection of an appropriate method for addressing these (see Sec. 2).

From a high-end perspective, uncertainties are often classified into two broad categories, namely, aleatory and epistemic [5,6]. Aleatory uncertainty is attributed to randomness and inherent variability (spatial, temporal, part-to-part, sensor noise), or in other words to factors that cannot be described with a single value and are therefore treated as stochastic. For instance, the wind or traffic loads that a bridge experiences are subject to inherent variability both in a spatial and temporal sense, which can be modeled as an aleatoric stochastic process. Epistemic uncertainty is attributed to missing information (lack of useful data or knowledge) on model parameters or inadequacy of modeling assumptions. Aleatory uncertainty is considered irreducible, in the sense that sufficiently rich data can only refine the description of the dependency between the model parameters/inputs and the outputs, which remains stochastic in nature. Epistemic uncertainty is instead considered reducible, e.g., by gathering additional data or by refining the adopted models, it would be possible in principle to identify the underlying deterministic, but unknown, parameter value or model. The classification of an uncertain quantity as aleatory or epistemic is oftentimes a subjective and somewhat philosophical issue. Also, an uncertain quantity can change from aleatory to epistemic depending on the problem configuration. For example, prior to a structure being erected, the concrete strength can be considered as an aleatory variable because of lack of control over available fabrication and manufacturing processes and the expected differences between the properties of the designed and as-built structure. Nonetheless, the same variable can be considered as epistemic as soon as the building is built, since detailed (e.g., nondestructive) investigations can be carried out to identify the actual concrete strength [5].

This broad classification of uncertainties in terms of aleatory and epistemic classes is not comprehensive in capturing all aspects of the essential dynamical process depicted in Fig. 1, including the crucial effects linked to time-varying phenomena. In this paper, we propose the adoption of a refined classification for characterization of the various sources of uncertainty within the context of structural dynamics, system identification, and SHM, with the latter serving as the primary focus area. We first review relevant literature that discusses different sources of uncertainty in engineering contexts where some system representation and measurement data are available. We next outline the terminology adopted in this paper for addressing different sources of uncertainty. Subsequently, we outline methods that attempt to account for the individual forms of uncertainty that characterize each component of the presented flowchart. We will discuss these methods through the prism of the following defining tasks within an SHM perspective:

• System identification and model inference [7],

• Model updating [8,9],

• Capturing of environmental and operational variability [2],

• State estimation/virtual sensing [10],

• Damage identification [8,11,12],

• Remaining useful life/prognostic health management [13].

This paper is organized as follows. Section 2 discusses thoroughly the uncertainty sources that are involved in every component of the process illustrated in Fig. 1 and presents the taxonomy of uncertainties adopted in this work. Section 3 reviews methods and algorithms from existing literature that address the uncertainties encountered in fundamental SHM tasks, through the lens of the taxonomy defined in Sec. 2. Section 4 employs a model of a benchmark shear frame with hysteretic links, serving as a running example for demonstrating methods and algorithms suitable for performing these fundamental SHM tasks. Finally, Sec. 5 touches upon open challenges and opportunities in uncertainty quantification for SHM and offers some concluding remarks.

In engineering applications, various types of uncertainty sources may be overlapping, making it difficult to distinguish individual contributions. We will now attempt to break down the sources of uncertainty in relation to the essential components comprising the dynamical system depicted in Fig. 1.

Here, we discuss those sources of uncertainty that relate to the input of the monitored system, or in other words, the loads/excitation that are acting onto the system. We further summarize the methods that can be used to tackle the associated uncertainties. In terms of their temporal profile, such inputs range from stationary broadband (e.g., white noise) to nonstationary excitation sources (e.g., wind gusts or earthquake signals), which may further operate concurrently (e.g., seismic excitation under traffic loading on a bridge). In terms of their spatial profile, they can be classified to either concentrated (point) or distributed loads. Examples of the former are loads from operating devices (e.g., shakers), which are typically positioned in a fixed location, or moving loads, such as axle loads from traversing vehicles on a bridge. Large scale structures are, on the other hand, most often exposed to distributed loads (e.g., wind on an aircraft wing or a bridge deck, wave loads along wind turbine towers).

The input to the system may be measurable; however, in most typical instances of operating engineering structures, the acting loads remain unobserved. Examples of measurable inputs pertain to measured actuator forces in controlled laboratory experiments, contact forces in instrumented bearings, or the earthquake excitation reaching the foundation of seismically excited structures. It should be explained that, even when the input quantities, or a subset thereof, are measurable, uncertainty is involved in the form of noise corruption due to the employed instrumentation, which is bound by prescribed noise levels (defined in related sensors documentation) and in the form of a bias, for example, outside the linearity range of the sensor. In such a case, it is typical to simulate this noise process on the basis of an assumed, typically white noise Gaussian distribution. When records of the inputs are available, these can be used to build probabilistic models of the input processes. For instance, in Ref. [14], it is shown that earthquakes can be modeled via parameterized stochastic models, using a time-modulated filtered white noise process with the filter having time-varying characteristics. In addition, under availability of known or controlled inputs, the system identification task is facilitated. The inference of a model of the system is accomplished via techniques such as experimental modal analysis (EMA) [15] and corresponding algorithms, which rely on input–output identification methods. Since, however, most loads acting on engineered systems are of a distributed nature, they are also not possible to precisely measure in a realistic operational setting (field conditions). When the input to the system is unmeasured, various schemes exist for tackling this challenge:

• Assume that the input follows a prescribed probabilistic distribution (e.g., Gaussian), and apply corresponding identification schemes, e.g., those under the class of operational modal analysis [16].

• Assume that it is described by a random process of known type, e.g., random walk, or a stochastic time series model, e.g., of the AutoRegressive form, which is assigned to the eXogenous input [17].

• Attempt to infer a more elaborate probabilistic model, describing the evolution of the input in the form of a stochastic differential equation, e.g., Gaussian latent force models [18]. This may call for use of sample data of the inputs.

• Attempt to reject the unknown input in the form of a disturbance or jointly estimate this along with the system's response (state), through appropriate observer designs [19,20].

We elaborate on the influence of this source of uncertainty in terms of the prescribed downstream tasks in Sec. 3. It is worth mentioning that the characterization of this input/excitation uncertainty is extremely important in the presence of nonlinearities in the system, since the level of energy injected into the system notably affects the nature of the system's response [21].

The presence of disturbances that are clustered under the class of EOV raises a significant source of uncertainty which often hinders the “health” assessment of monitored systems [22–25]. EOV, often also referred to as varying environmental and operational conditions, refers to the variability in the (i) environmental conditions within which a structure operates, such as temperature, humidity, alkalinity (pH), as well as the variance of the common operational loads, including wind speed and direction, traffic loading, and rotor revolutions per minute. The impact of operational conditions, including loading and environmental factors, on monitored structural dynamics has been extensively studied over the past years. In bridges, in particular, temperature variations affect structural stiffness, as well as support conditions (e.g., bearing performance) with notable impact on the structure's dynamic properties. Typically, this influence is in a benign and nonharmful manner.

In the domain of bridge monitoring, Farrar et al. [26] observed natural frequency variations of the order of 5% over a 24-h period for the Alamosa Canyon Bridge, attributed to changes in temperature and spatial temperature gradients. Alampalli [27] reported up to a 50% difference in natural frequencies for an abandoned bridge in Claverack, NY, due to the freezing of its supports. In the context of the SIMCES project [28], a maximum variation of 14–18% was observed in the first four natural frequencies of the Z24 Bridge over nine months. Catbas et al. [29] demonstrated that ambient temperature significantly affected the structural reliability of a long-span truss bridge in the U.S. Cross et al. [30] identified variations of up to nearly 5% in the modal properties of the Tamar suspension bridge in the UK due to temperature, traffic loading, and wind speed. Martín-Sanz et al. [31] report on the influence of operational loads in the identification of bridge structures, explaining that nonstationary behavior can arise, which calls for use of dedicated system identification tools that account for time-varying response.

Despite the EOV influence being often reported for the case of bridge structures, this variability affects further structures, including tall buildings, as reported in Ref. [32], where one-year measurements of a 22-story reinforced concrete building are analyzed, and a relationship is established between natural frequencies and environmental factors, such as temperature and humidity. Finally, perhaps the strongest case to be made for dependence of structural response on EOV is the case of wind turbine systems, where the EOV factors are logged in the form of statistical averages that are logged as supervisory control and data acquisition [33,34].

Environmental and operational variability contributes to the variability of the system's response and affects its output and thereby the acquired measurements. It can further interact with nonlinearities in the system, giving rise to complex response effects [21]. EOV impacts the capacity to solve inverse problems that require identification and tracking of the structure's response (e.g., damage detection tasks). Thus, in terms of SHM-based assessment, EOV acts as a confounding process and must therefore be properly accounted for or disentangled from other contributions to the system output. To tackle this, various model-based or data-driven methods are available, as discussed in Sec. 3.3.

Deterioration processes of various forms (e.g., corrosion, fatigue, and wear) can adversely affect the intended performance of engineering systems. Various sources of uncertainty influence the deterioration process h(t), mainly relating to environmental factors, loading conditions, or material properties. We here classify the result of these processes under the term deterioration, i.e., a category that is separate to the EOV (d(t)), since we here refer to those mechanisms that are now internal and intrinsic to the system.

Deterioration modeling is a rather complex task. With respect to the mathematical form describing the deterioration process h(t), deterioration models can be classified into two categories, empirical models or physics-based models [35]. Owing to the uncertainty of the influencing agents, a stochastic approach to deterioration modeling must be adopted. With respect to the way in which stochasticity is introduced in the mathematical form describing the deterioration process, empirical and physics-based deterioration models can be classified into random variable models and stochastic process models [36]. If a deterioration model is employed, deterioration model form uncertainty and model parameter uncertainty are introduced, as later discussed in Sec. 2.5. In certain cases, a deterioration process can be directly monitored via use of dedicated sensors—albeit in most cases indirect measurements are obtained. In either case, measurement uncertainty is naturally prevalent.

Deterioration processes directly affect the validity of the LTI assumption that is often adopted in a simplified approximation of operating systems. It becomes evident, that under influence of deterioration phenomena, model parameters can no longer be assumed as time-invariant, and dedicated uncertainty quantification strategies need to be implemented [37]. However, depending on the time-scale of the data being collected, it could be assumed that deterioration is a slow-varying phenomenon with respect to the duration of the monitoring campaign. In this case, its influence may be considered in a (long) scale that is different to the short scale of the model and measurements, as in the case of vibration-based monitoring [23].

When exposed to the set of inputs and influencing agents ${u(t),d(t),h(t)}$⁠, the dynamical system produces an output y(t), which reflects the set of quantities that describe a system's response, e.g., displacement or rotation, velocity, acceleration, reaction force, or stain. The system outputs can be fully or partially measured. In the partially measured case, we can distinguish between the measured y_s and the unmeasured subset y_u. The entire subset can form target of estimation schemes. For the monitored subset y_s, the goal of estimation may be to estimate the output quantities with a higher precision than the physical measurement, or to perform a step-ahead prediction for diagnostic and control tasks. When the main target is to estimate (predict) quantities that lie in the unmeasured subset y_u, then this is termed a virtual sensing task [38]. We regard the uncertainty surrounding the estimated quantities as an outcome of the modeling process and not a “source” of uncertainty; we thus do not classify this as a class in the taxonomy offered in Table 1.

The subset y_s reflects output quantities that are monitored via sensors, with a multiplicity of sensing systems today available for monitoring these diverse response quantities [44]. Sensor measurements are in turn associated with measurement uncertainty, which consists of random measurement noise (e.g., electronic noise, thermal noise, etc.) and/or systematic errors or biases (e.g., due to calibration issues, sensor drift, etc.) [11]. The raw measurement data are typically processed by subsequent signal processing or system identification algorithms [45], which can contribute to the presence of systematic errors or biases. The random measurement noise is an aleatory irreducible uncertainty, whereas systematic errors or biases classify as epistemic uncertainties, which can be reduced if properly accounted for. The level of random measurement noise can significantly impact SHM-based assessment, e.g., high levels of noise can make it challenging to differentiate between small changes caused by damage and the inherent variability in the measurements due to noise. Added to this is the challenge of sensor faults, which occur in a nonsystematic manner, raising false, which negatively impact diagnostic tasks [46,47].

Sensor arrangements are typically sparse, necessitating careful consideration of sensor placement in order to ensure that the selected sensor locations effectively capture relevant information on the system's behavior. Sensor placement methodologies can be classified under two main schemes, namely, information-theoretic [48,49] or decision-theoretic approaches [50,51]. The information-theoretic schemes employ rigorous criteria for optimal sensor configuration from an information-theoretic point of view (e.g., via use of the Fisher information matrix), relying on the principle that more sensors lead to higher information gained and, therefore, improved damage identification. On the other hand, decision-theoretic criteria attempt to balance the information gain with further important or pragmatic design parameters, such as the cost of the monitoring system, bringing in concepts such as the expected value of information [50]. A placement scheme, which attempts to balance information with available budget, ought to account for correlation among measurements both in the spatial (e.g., due to spatial dependencies in the system dynamics) and temporal sense (e.g., due to temporal trends) [52–54]. Furthermore, the heterogeneity of available monitoring solutions and measurements, often obtained at different sampling rates with different type of sensors, adds a further layer of complexity. Addressing these issues is essential for robust uncertainty quantification in dynamical systems.

In this section, we refer to those uncertainties that arise on the basis of the model configuration, which can pertain to both the model form, as well as the defining parameters of a prescribed form. The seminal work of Kennedy and O'Haghan [40] presents a fundamental classification of the various sources of uncertainties that are present when using computational models to describe different processes in science and engineering. Parameter uncertainty describes the uncertainty corresponding to the model parameters. In our context, such parameters could correspond to the system's characteristics/properties, but even to parameters that are used to describe the input (e.g., spectral or amplitude properties). Model inadequacy refers to the fact that a computational model cannot perfectly capture the real process, i.e., no model is perfect. Specifically, the authors define model inadequacy as the difference between the true mean value of the real-world process and the model output evaluated at the true model input values. Residual variability refers to variations of the real process, even in cases when the conditions specified by the model input remain fixed. This may be because the real process is inherently stochastic, or because the model lacks specification of additional, unrecognized conditions. The latter case is very closely related to model inadequacy. Parametric variability relates to the fact that the model input parameters will inherently vary with time, e.g., under influence of a deterioration process. Finally, in the taxonomy of Kennedy and O'Haghan [40], code uncertainty emerges due to the fact that, typically, the mathematical model needs to be implemented in a computer code. O'Hagan and Oakley [55] further categorize the sources of uncertainty presented in Ref. [40] to either aleatory or epistemic.

In a similar vein, Simoen et al. [11] discuss the different sources of uncertainty related to the prediction model within a damage assessment setting. Specifically, the authors separate between: (i) model parameter uncertainty (same as parameter uncertainty of Kennedy and O'Haghan [40]), (ii) model structure uncertainty, and (iii) model code uncertainty, i.e., the uncertainty that results from errors in the computer implementation, which is usually ignored or assumed part of the model structure uncertainty. Model structure uncertainty arises as a result of modeling assumptions and simplifications made by the designer/analyst due to lack of knowledge or understanding of the true system (this relates both to model inadequacy and implicitly also to residual variability from the taxonomy of Kennedy and O'Haghan [40]). Examples include the uncertain extent of a numerical model, potentially invalid assumptions regarding boundary conditions, model order, and governing equations. It should be noted that models that are inferred via data-driven or system identification procedures, also suffer from the curse of inadequate model structure, and can thus be classified under this category, since the deduced form is limited to capturing the information contained in their trained datasets. Indicatively, we mention the popular use of stochastic subspace identification or other operational modal analysis (OMA) methods for inferring a linear state-space model or the modal properties of monitored dynamical systems [56], or the use of black-box methods for nonlinear model inference from monitored data (e.g., the deep Markov model [57]). Finally, Simoen et al. [11] further place finite element (FE) discretization errors or approximation errors in the model structure uncertainty category.

The aspect of model inadequacy and model structure uncertainty, which appears in the aforementioned taxonomies, has more recently been referred to as model form uncertainty; a term that is more broadly adopted in our context [58–61]. As discussed in Refs. [62] and [63], model form uncertainty can arise owing to wrong physics assumptions and physics-based model approximations, corrupted or uninformative measurements leading to an incorrect data-driven model, or a combination of both. Some parameter and model form uncertainties might be irreducible because of Ref. [62]: (i) impossibility of carrying out destructive testing, inaccessible sensors locations, lack of access to experts and/or information, unavailability of informative data, or (ii) impossibility of distinguishing confounding sources in the measurements because of ignorance. Such irreducible uncertainties might lead to poor model prediction performance (e.g., low variance and high bias), underestimation of the overall output uncertainty, and poor inference on future health conditions. This type of uncertainty hampers extrapolation potential, limiting forecasting and prognostic capabilities, thus negatively impacting decision-making under uncertainty.

We henceforth adopt the following terminology and definitions for the different sources of uncertainty within the context of dynamics and SHM applications, through the prism of the representation of our reference high-end model, depicted in Fig. 1. We recast that representation in the context of the proposed uncertainty taxonomy in Fig. 2.

• Input uncertainty; under this class, we cluster any uncertainty and associated assumption related to the dynamic inputs (load/excitation) to the system, which are assembled in the vector $u(t)$⁠. Given that the input itself can be modeled via a function, this uncertainty comprises both a statistical and an input-model form uncertainty. We represent this as a parameterized probabilistic input submodel $fθu(u|Mu)$ describing the distribution of the random vector $u$⁠, which can be dependent on a physical input submodel $Mu$⁠. The model, when supplied, describes the assumed model of the input with respect to time and space $u=Mu(t,θu)$⁠. Such a model, which is here parameterized by the parameter vector $θu$⁠, can be a mathematical model, e.g., a dynamic evolution model, a model of the spectrum of the input, or a model of its spatial distribution. Under this definition, uncertainties related, for example, to the input matrix B and feedthrough (or feedforward) matrix D of a state space model would be considered as model form uncertainty.

• Model form uncertainty arises from the unknown functional form of the underlying model of the system, which aims to reproduce the outputs $y(t)$⁠. This can be represented by a physical system model $y=M(u,x,ξ)$ describing the relationship between the input quantities $u(t)$⁠, the possibly present system's latent variables $x(t)$⁠, and the derived output quantities $y(t)$⁠. In our context, the parameter vector $ξ$ can correspond to the system's characteristics/properties, e.g., material properties, geometry, joint properties, and damping/energy dissipation mechanisms. Thus, by model form uncertainty, we here specifically refer to the uncertainty arising due to the selection of the functional form—or structure—of the system model $M$⁠, whether this is purely data-driven or physics-based. This may refer to the selected physics-based model, which may not fully capture the physics of the problem at hand (e.g., unknown term(s) in the governing equations, use of simplified model), to inadequacy of a data-driven model approximation, or to implementations which suffer errors in the computer implementation/solver. Within this definition, model form uncertainty accounts for the many plausible model choices and model parametrization. As a result, the model hyperparameters (set manually and/or to be tuned) are also included in this uncertainty. It should be noted that the effect of deterioration h(t), which is assumed to intrinsically affect the system, is also included in this class, should the modeler opt for approximating the deterioration process via a model. In a more extended context, within a Bayesian model updating (BMU) setting, it may further pertain to the choice of the likelihood function describing the discrepancy between the observations and model prediction (form of probabilistic submodel).

• Model parameter/variable uncertainty refers to the statistical uncertainty associated with the parameters and variables of the employed model of the system, once its form is selected. It can be represented through definition of the probabilistic model $fX(x)$ of the (possibly present) latent variables $x(t)$ and the probabilistic model $fΞ(ξ)$ of the parameter vector $Ξ$⁠, which parameterizes the functional form of the system model $M$⁠. The effect of deterioration can also be reflected here on the basis of its influence on the model parameters (e.g., stiffness and strength). On the other hand, the dynamic latent variables $x(t)$⁠, which depending on the model form may or may not exist, describe the system's latent dynamical state.

• Measurement uncertainty reflects the uncertainty associated with sensor noise, typically represented as a random process $v(t)$⁠, and/or systematic errors (bias). This uncertainty class affects any measurable quantity of the system described in Fig. 1, including the inputs $u(t)$⁠, outputs $y(t)$⁠, and even the disturbances $d(t)$⁠.

• Inherent variability is defined as the inherent uncertainty linked to variability in the inputs, disturbances, and deterioration processes of the system. This class pertains to irreducible aleatory uncertainties. It encompasses uncertainties that operate at a temporal scale that is slower than the scale describing the dynamics of the operating system. Similarly to the case of the input vector, we represent this uncertainty via a parameterized probabilistic input submodel $fθd(d|Md)$ describing the distribution of the random vector d, which can be dependent on a physical input submodel $Md$⁠. The model, when supplied, describes the assumed model of the disturbances/EOV parameters with respect to time and space $d=Md(t,θd)$⁠.

In Table 1, we relate the terms adopted in this work to categorize the different sources of uncertainty with the equivalent terms adopted in other original contributions from the literature.

We proceed with a review of available methods and algorithmic tools that have been put forth in existing literature for tackling the uncertainties appearing in our above-defined taxonomy, now classified under the prism of defining SHM tasks, as listed in the introductory Sec. 1.

Model inference forms a primary task within the domain of SHM, typically described under the term system identification, which describes an inverse problem formulation that aspires to deliver a mathematical (and not necessarily physics-based) representation of the system, on the basis of the available monitoring data [1]. In this sense, this task is inherently linked to the model form uncertainty of our described taxonomy. System identification pertains to broader class of methods, attempting to infer the input–output relation, or in other words the model $M$⁠, of a considered system. These can be clustered under the broad categories of purely data-driven and hybrid, also referred to as model-based schemes.

Such methods can be classified into parametric and nonparametric schemes. Parametric methods rely on consideration of a parameterized system model (e.g., of the transfer function or state-space form). Nonparametric methods deliver defining system properties, such as modal characteristics, without imposing a specific structure, with the system regarded as a black box. The primary SHM task, for which both parametric and nonparametric methods are deployed, is operational modal analysis. A listing of typically adopted such schemes for the purpose of OMA, along with a further classification in terms of their domain of operation, i.e., the frequency versus the time domain, is offered in Table 2.

These methods are further defined by assumptions that relate to input uncertainty. As aforementioned, OMA schemes assume ambient, i.e., broadband inputs which are assumed unmeasured. On the other hand, input–output methods (such as EMA), assume availability of measurements of the input, linking to the measurement uncertainty class. Both OMA and EMA, its input-aware counterpart, impose the assumption of linear and time invariance. Uncertainty is further introduced via the assumption that the system response is a realization of a Gaussian distributed stochastic process that has zero mean, liking to the model parameter/variable uncertainty class. To what concerns linear systems in particular, a suite of methods has been put forth for tackling mode form uncertainty [58], but for further also quantifying the propagated uncertainty to the output quantities, such as modal properties [45]. The interested reader is related to the works of Döhler [77] and Greś et al. [56,78] for extensive discussions in this respect.

However, real operating structural systems often experience more complex effects, such as nonlinearity or nonstationarity. Nonlinearity is met in the case of structural response to pronounced (e.g., seismic) loads, which push structural sections beyond their yield/ultimate capacity, or can be induced by large deformations in deformable structures such as cables. Nonstationary is often linked to the nature of loads a structure (e.g., a wind turbine) is subjected to. For tackling such cases, refined alternatives of parametric and nonparametric system identification tools again exist. The parametric class includes tools such as nonlinear autoregressive (AR) models with eXogenous input [79], as well as linear parameter varying (LPV) models for time-varying structural response modeling [80]. In such a case, uncertainty is introduced by the model's residual sequence, which is modeled as a random process assumed to be described by a given probability distribution. Nonparametric techniques that are suited for nonlinear and nonstationary tasks include wavelet-based approaches and methods utilizing the short-time Fourier and Hilbert transforms, which can as well be combined with stochastic processes [81].

The former summary pertains to methods that attempt to infer a system model, $M$⁠, purely on the basis of data availability. However, within the context of SHM, hybrid constructs are also employed. Hybrid methods exploit a model of first principles in the identification loop. These methods are also referred to as model-based or, when used in conjunction with machine learning (ML)-based schemes, as physics-enhanced techniques [62,63,82,83]. Engineering applications, particularly within the domain of civil and mechanical engineering, often rely on use of numerical approximations of a system, often delivered in the form of FE representations. In these instances, system identification typically involves estimating physical model parameters to ensure that the numerical model accurately represents the monitored structural behavior. This procedure is more appropriately termed model updating or model calibration and is further elaborated in Sec. 3.2.

On the other hand, the balance sought in a hybrid construct between the selected model of the physics and a learnable representation that is captured on the basis of data-driven (often ML-based) learners is linked to model inference and the model form uncertainty class of our taxonomy. The derivation of such hybrid representations in a way that balances availability of the physics-based knowledge with the data at hand is offered in prior work of the authoring team. We refer the interested reader to Haywood-Alexander et al. [63] for a discussion on the spectrum of available methods merging physics-based models and data-driven learners, and to Cicirello [62] for a discussion as to how these models are linked to the sources of uncertainty identified in our taxonomy, including parameter and model form uncertainty. Available approaches include schemes that try to discover the underlying model form (e.g., dictionary methods [84–86]), schemes that account for model discrepancy (e.g., Gaussian processes [83,87,88]), methods that attempt to inject physics-based representations in the architecture of data-driven learners (physics-guided schemes [89]), architectures that are fit to describe the dynamics observed (neural ODEs [90] and symplectic encoders [91,92]), methods that impose the physics in a more loose way in the objective function of the problem (physics informed neural networks [93]), and more. The discussion in this chapter can be extended to include further aspects, including the aspect of model selection and/or model aggregation [94], which however, we do not address here in detail. See Ben Abdessalem et al. [95] for an elaboration on this theme.

In the context of dynamics and SHM, model updating refers to the process of calibrating the computational model of a system based on data collected from sensors [96]. The goal is to improve the accuracy of the computational model by identifying the parameters of a prescribed model form, so that the updated model would better capture (reproduce) the actual behavior of the structure. In this respect, this downstream task is more closely linked to the model parameter/variable uncertainty class of our taxonomy. The updated model can then be used to investigate the dynamic system behavior under various inputs and EOV conditions, including extreme events. Uncertainties in the model updating process are often accounted for via adoption of a Bayesian framework [97]. The goal of BMU framework is to infer latent uncertain model parameters, and quantify their posterior uncertainty [8,11,40,98]. This is achieved by specifying a “prior distribution” that encodes prior knowledge, experts' opinion, physics constraints and prior data on those parameters that cannot be directly observed or measured, and a likelihood function which models the discrepancy between some observed data and the corresponding model prediction. This discrepancy can be modeled by considering additive or multiplicative discrepancy errors and measurement noise.

Bayesian model updating frameworks typically rely on evaluating the posterior distributions via use of asymptotic methods such as the Laplace approximation [99], Markov chain Monte Carlo (MCMC) methods [9,98,100,101], or variational inference approaches [102,103]. Several works [12,104,105] make a distinction between model parameter estimation uncertainty and inherent variability (e.g., EOV). They discuss that the classical BMU framework does not properly account for inherent variability and modeling errors. Specifically, they demonstrate that the classical BMU scheme leads to underestimation of the true variability in the posterior model parameter estimates, with the uncertainty bounds becoming unrealistically narrow when multiple datasets are available. To tackle this, they propose a hierarchical BMU framework.

An important task within the context of SHM lies in quantification—or removal—of the influence of disturbances, primarily attributed to EOV, in the context of estimation and assessment tasks. This task is linked to the inherent variability uncertainty class of our taxonomy, which typically aims to effectively capture the influence of long-term trends underlying the induced dynamics [80,106]. When the connection to EOV is established via adaptation of the system model, then the relevant uncertainty taxonomies are also influenced. Different strategies have been put forth for tackling this challenge. In the case of data-driven schemes, these can be roughly classified in three main categories, namely, (i) nonstationary time series models, (ii) projection methods, and (iii) deterministic or stochastic functional dependence models.

Nonstationary time series models are constructed as output only, or input–output time series representations, which account for time variability such as smoothness priors time-varying autoregressive models or LPV models [107]. Projection-based schemes, also known as data normalization methods, aim at projecting features derived from system measurements onto a subspace where the EOV influence can be removed [108,109]. Characteristic such strategies are found in the cointegration approach proposed by Cross et al. [106], as well as the kernel principal component analysis approach suggested by Rainieri et al. [110] in the context of SHM. These methods offer the advantage of not requiring explicit measurements of the EOV parameters, as opposed to methods relying on inference of dependence models. Deterministic functional dependence models assume a deterministic functional relationship between EOV parameters and characteristic features of model $M$ describing the dynamics of the response [22,111,112]. Deterministic methods fail to account for the uncertainty associated with EOV. Random or stochastic functions are more appropriate for such a task. Instances of such models, associated with time series dynamics modeling, are found in random coefficient [113], as well as polynomial chaos expansions and Gaussian process-based time series models. Such instances include polynomial chaos expansions-AR type models [33,79] and GP-LPV-AR models [43,80] and time-varying autoregressive-GPRMs [114]. Such models are capable of accounting for the global, i.e., both short- and long-term, identification of the dynamic response of a structure by representing the short-term dynamics via a linear-in-the parameters regressive time-series model and a stochastic model, which takes on the stochastic dependence on the (measured) EOV parameters.

State estimation is a fundamental problem in control theory, signal processing, and SHM, where the goal is to infer the internal states, i.e., the response, of a dynamical system from noisy and potentially incomplete observations y(t). This is of particular importance for SHM purposes, since it is linked to the virtual sensing task [38,115], or in other words, the inference of the response of the system in degrees-of-freedom that are not measured. Returning to our reference system in Fig. 1, this implies that we try to infer the latent dynamics contained in $x(t)$ through measurements/observations of the output, which we assume here contained in $y(t)$⁠. This task is typically accomplished by observer formulations, which are today available in the form of conventional linear or nonlinear models as well as deep learning variants. The classical formulations are typically delivered in the form of Bayesian filters, which essentially form dynamic Bayesian networks [116]. These couple a model of the system $M$⁠, which can be either physics-based or inferred via the system identification tools outlined in Sec. 3.1. Within a virtual sensing setting, these observers can be exploited to deliver a number of tasks of increasing complexity. These include (i) response (state) estimation, curbing model variable uncertainties, (ii) joint or dual state-parameter estimation, thereby curbing model parameter uncertainties [117], (iii) input-state estimation, curbing input uncertainties [118–120], and even (iv) joint state-parameter-input identification [121,122]. The delivered estimations can be further processed to achieve diagnostic tasks, such as damage detection [123,124].

in which a transition model$pθt(x1:T)$ parameterized by $θt$ is considered for describing how the dynamical system (or process equation) evolves over time. Correspondingly, an emission model$pθe(y1:T|x1:T)$ parameterized by $θe$ can be established for governing the relationship between the observed variables and the latent states. $θ=θt∪θe$ is vector of all parameters involved, formed via concatenation of the transition and emission parameters.

Under the assumption that the state vector x_k reflects a Markov process, the transition model can be prescribed according to a transition probability density $p(xk|xk−1)$⁠. This model can be expressed as $xk|xk−1=xk∼p(xk|xk−1)$⁠. The observations y_k are then described by the following conditional distribution $yk|xk=yk∼p(yk|xk)$⁠.

Bayesian filters leverage this hybrid formulation to produce a refined posterior estimate of the system's response, $xk$⁠, using a “predict” and “update” procedure. The description of the details of these implementation lie outside the scope of this work. However, what is important to stress is that the finally updated estimate is obtained as a weighted combination which accounts for the process $w(t)$ and measurement noise sources $v(t)$⁠, as well as the discrepancy (innovations) between the measurements and the model predicted estimates. In this way, the final prediction of the system's states and or parameters are affected by model, parameter, input, disturbance, and measurement uncertainties. Variants of these filters are designed to handle linear systems (Kalman filter—KF) or nonlinear systems (extended KF, unscented KF—UKF, particle filter, etc.) for various estimation tasks. Additionally, when coupled with appropriate reduced order models, Bayesian filter estimators can operate in real-time or near real-time [125,126]. More details on the implementation of these filters can be found in Refs. [117] and [127], whereas a python library is available in conjunction with a tutorial on nonlinear Bayesian filtering [128].

The conventional formulation can tackle model parameter/variable uncertainty through the joint state-parameter estimation setup, but is less robust in handling model form uncertainty. This challenge can be tackled by the deep learning variants of these filters, such as deep Kalman filters and deep Bayesian filters, where the nonlinear model $M$ equations are also inferred. This can be executed within a black box setting (no physics information) as accomplished in the deep Markov model [57], or alternatively within a gray box, physics-enhanced scheme, as proposed in the physics guided deep Markov model [89].

Moreover, it is possible to tackle input uncertainties, by modifying these filters for use in absence of input $u(t)$ measurements in what was already referred to as joint state-input estimation schemes [129]. This can be achieved on the basis of (i) a random walk evolution assumption, or—when feasible—by seeding some knowledge of the expected form of the input (e.g., periodic [130]), or via inference of a model that is suited to describe the evolution of the input. The latter has been tackled via adoption of Gaussian process latent force models [18,118,131].

Damage identification within the SHM context is typically organized across the levels defined in the hierarchy specified in the work of Rytter [132], namely, (i) detection, (ii) localization, (iii) assessment (quantification), and (iv) consequence. Most commonly, damage identification is narrowed down to the first three levels, which correspond to diagnostic tasks, whereas the consequence, which reflects a prognostic task, is separately treated. Also here, we discuss the fourth level in Sec. 3.6. Damage identification forms a complex task that can be achieved on the basis of both purely data-driven as well as hybrid schemes and typically mobilizes the whole construct depicted in Fig. 1, along with the corresponding sources of uncertainty of our defined taxonomy in Table 1.

Both purely data-driven and hybrid schemes rely on the extraction of appropriate damage sensitive features (DSFs), whose shifts with respect to a baseline serves for detecting damage. The system identification tools, which were overviewed in Sec. 2.5, serve for deriving such DSFs, as elaborated in Fig. 3. Modal properties, such as natural frequencies or associated modal shapes, serve as a primary instance of DSFs. Detection is achieved via detection of shifts in these properties when compared against the values observed under normal “healthy” conditions. Given the global nature of such properties, the sensitivity of such features against local damages is often questioned. Moreover, if a network of sensors is installed, localization is further possible and will naturally depend on the spatial resolution and type of sensors deployed. The task of quantification is often more challenging and benefits from the availability of a physics-based model (hybrid schemes) [133].

The characterization of damage forms in essence a classification task, which can be achieved using unsupervised, semisupervised, or fully supervised schemes [134]. Within the context of SHM for civil structures, the fully supervised scheme is often not exploitable, since extensive labels of damages of various types and severity are hardly available. Thus, unsupervised schemes are more valuable and often pursued with priority. These can be achieved using either more classical frameworks or methods exploiting machine learning. Classical frameworks include statistical hypothesis testing [135–138], which relies on quantification of the uncertainty of the damage sensitive features and characterization of their probabilistic distribution. Crucial to this task is the complementary task of efficiently capturing, or normalizing for, the dependence on influencing EOV parameters, as described in Sec. 3.3.

Data-driven methods are appealing for real-time damage identification, but often struggle to achieve higher levels of the damage identification hierarchy due to the absence of a physics-based model of the system. In contrast, hybrid schemes exploit the derived DSFs to estimate potential damage severity and its type/position, usually on the basis of the model updating procedure, which was described in Sec. 3.2. While these methods are less suited for real-time (or fast) damage identification tasks due to the high computational cost, the availability of a physics-based (e.g., FE) model allows for more detailed damage assessment [8,12]. It is worth mentioning that for alleviating computational costs, within an inverse setting for the purpose of damage detection, an emulator/surrogate of an expensive physics-based model can be adopted [99,139,140]. An overview of the vast literature on damage identification lies outside the scope of this paper. We refer the interested reader to Refs. [64,141], and [142] for an overview of state of the art schemes for damage detection in the context of SHM.

Prognostics focus on predicting the time instance when a dynamical system will fail to function properly. In an online monitoring setup, prognostic models are developed on the basis of SHM and related data. These models, which may be hybrid or purely data-driven, deliver predictions of the future evolution of a system's condition (e.g., predictions of the system's remaining useful life). The prognostics process is characterized by several uncertainties, necessitating the adoption of stochastic approaches [143]. However, this aspect has often been overlooked in the literature, where deterministic ML-based prognostic approaches are abundant. For a detailed discussion on the multiplicity of uncertainty sources, the reader is referred to Ref. [144]. One basic aspect to consider in prognostics relates to how the uncertainty regarding the system's current state is going to propagate in the future. This typically requires models of the future loading conditions and the deterioration processes. Health management refers to the subsequent task of decision support for maintenance planning on the basis of the prognostic predictions. Needless to say, the health management process must properly account for the uncertainty in the prognostic predictions and must be formulated as a problem of decision-making under uncertainty [145]. The overall process from prognostic predictions to health management is referred to as prognostic health management [146,147].

As a demonstration example, we consider a FE model of a three-dimensional two-story shear frame with nonlinear joints, each exhibiting a Bouc–Wen type hysteretic behavior [148,149]. This nonlinear multidegree-of-freedom simulator has been published as an open-access benchmark in Refs. [150] and [151] and is available through the following GitHub repository. The simulator is flexibly reconfigurable allowing for assumption of a different number of storeys and bays. The main features of this benchmark system are first described. Subsequently, it is employed to explore uncertainty quantification challenges in five downstream tasks.

Each node of the frame comprises 6DOFs, corresponding to three displacements and three rotations. A graphical illustration of the shear frame is given in Fig. 4, where the beam local axes are noted and where the hysteretic joints are modeled as virtual elements of zero length, while a negligible mass instead of zero mass is also assumed to avoid ill-conditioned system matrices. The elements are illustrated at an eccentric distance from the reference node in Fig. 4 solely for demonstration purposes. Each story of the deployed frame has a height of $h=3.2 m$⁠. The configuration which was adopted has two frames along the x axis, each of $l=7.5 m$ length and one of $w=5 m$ along the width. All beam and column elements comprise a rectangular steel cross section (⁠$40 cm×40 cm$⁠). The Young modulus for the examples demonstrated herein is set equal to $E=210 GPa$⁠, the Poisson ratio is v = 0.30, and the density $ρ=8000kg/m3$⁠. Furthermore, 4% mass proportional damping is assumed, and the structure is assumed fully constrained at the bottom nodes that represent the ground.

where the shape, smoothness, and amplitude of the hysteretic curve that characterize the dynamic behavior of each joint is determined by $β,γ,w$⁠, and A. The terms $ν(t),η(t)$ are introduced to additionally capture strength deterioration and stiffness degradation effects via the corresponding coefficients $δν$ and δ_eta. In turn, their evolution in time depends on the absorbed hysteretic energy, $ϵ(t)$⁠. This formulation allows for a structural dynamics simulator, which can be parametrized with respect to system properties and traits of the joints' behavior. For a more detailed elaboration on the physical connotations of the Bouc–Wen model parameters in terms of yielding, softening, and hysteretic behavior effects, the reader is referred to Chatzi and Smyth [117] and Vlachas et al. [150]. Regarding the forward-in-time integration of the governing equations, the simulator employs a Newmark scheme, while a forward Euler integrator is utilized to evaluate the hysteretic response of the joints. The respective time discretization is chosen to ensure stability and convergence.

For the purpose of system identification, the frame is excited on its linear regime, i.e., it is ensured that the aforementioned nonlinearities are not activated during simulation. Within its linear subdomain, the structure is characterized by 72 normal modes, ranging from 3.85 to 433.70 Hz. This band contains several pairs of closely spaced modes, e.g., three separated less than 0.05 Hz, which renders the inference process a quite challenging task, in terms of EMA. Energy dissipation is modeled via the Rayleigh damping scheme, by adopting 2% modal damping ratio for the first two modes.

In demonstrating the propagation of uncertainty during the system identification procedure, we conduct a Monte Carlo campaign, consisting of 100 independent experiments. Each single experiment is realized using a zero-mean Gaussian process as the base excitation (with $θ=π/4$ and standard deviation equal to 0.01 g), and the absolute vibration accelerations along x and y directions at point O2 (see Fig. 4), corrupted at 5% noise-to-signal (N/S) ratio, as the available structural responses. To avoid aliasing, we set the integration step as $Ts=0.0078$ s. The simulations are conducted for 255 s, and the resulted input–output data records consist of 32,768 samples.

We limit our analysis to the $[0,50]$ Hz band, which already contains 32 translational and rotational normal modes. To this end, all input–output data pairs are first lowpass filtered (via a 17th order Chebychev Type II digital filter, with less than 3 dB of ripple in the passband, 80 dB attenuation in the stopband, and 50 Hz cutoff frequency), and accordingly resampled at F_s = 128 Hz. The final datasets have length equal to 4096 samples.

Figure 6 displays the Monte Carlo coherence functions, estimated via Welch's spectral method (window size 1024, zero-padding at 2048 samples, 50% overlap, Hanning window). High coherence is observed within the frequency band of interest, whereas the dispersion of the Monte Carlo estimates remains rather low. Yet, there exist two low-coherence areas in both directions (around 5.44 Hz and 4.18 Hz, for the x and y axes, respectively), associated with rotational DOFs, pointing out a potential inefficacy of the data in modeling these modes.

The Monte Carlo frequency response function (FRF) estimates, again estimated via Welch's spectral method (window size 1024, zero-padding at 2048 samples, 50% overlap, Hanning window), are depicted in Fig. 7. As in the coherence estimate, the nonparametric FRFs show small dispersion in both amplitude and phase within the frequency band of interest, while there exist several peaks, with the magnitude in y axis appearing “richer.” However, not all 32 normal modes are visible in either axes; this is attributed to (i) the “local” behavior of some modes, which can be detected only by placing sensors at certain nodes of the frame, and (ii) to the adoption of Rayleigh damping, which causes the modal damping ratio being increased with increasing frequency.

It is noted that the above spectral estimates are characterized by inherent uncertainty, since they are based on realizations of stochastic processes and therefore they are also random variables. The uncertainty associated with these estimates has been quantified as bias and random error (i.e., statistical sampling errors of two types). Formulas for evaluating these errors can be found, for example, in Ref. [153] (Table 19.4, Chap. 19).

and N = 3000 is the length of the estimation dataset. The latter is formulated after dropping the initial 512 data (2 s), in order to avoid transient dynamics. ARX model estimation is carried out using linear least squares.

The linear least squares estimation of ARX$x(30,30,0)$ and ARX$y(42,42,0)$ models for the x and y axes, respectively, returns the parameter vectors visualized in Fig. 9, in terms of mean values ± standard deviations for every AR and exogenous polynomial coefficient. Recall that (i) the linear least squares problem, which is treated by formulating the normal equations from the associated loss function (the variance of the residuals), is a convex optimization one and, as such, is characterized by a global minimum; and (ii) the statistical conditions required for the successful ARX model estimation (stationary and persistent excitation) are fulfilled. Therefore, one would ideally expect that the estimated ARX models from the Monte Carlo simulation data would return “identical.” Still, the standard deviation of many parameters is non-negligible and demonstrates how uncertainty propagates from (noise-corrupted) data to models.

The uncertainty in the parameter vector is further inherited to the estimated vibration modes, shown in Fig. 10. We observe a certain consistency in the estimate of natural frequencies, although there exist several erroneous ones above 50 Hz, which appear as structural. Damping ratio returns quite scattered; an issue that is well reported in the literature.

under the assumption of stationary excitation. The estimation procedure is identical to the ARX one and conducted also via linear least squares. For illustrative purposes, we herein report the resulted vibration modes along the x direction, displayed in Fig. 11. The natural frequencies appear more dispersed, compared to the ARX-based ones in Fig. 10(a), whereas the modal damping ratios are completely scattered.

in which $Tspan=30 °C$⁠, f_l is the long time-scale frequency, and $wT[tl]∈U(0,2)$⁠. Figure 12(a) plots a realization of $T[tl]$ for $tl=1,2,…,365$ and $fl=1/365$⁠, while its sample and kernel-based probability density functions (PDFs) are displayed in Fig. 12(b). The induced modulus is shown in Fig. 12(c), and the variability of the first three vibration modes is depicted in Fig. 13. It is noted that Eq. (9) is a fictitious equation, which is nonetheless chosen to represent the decrease in stiffness with increasing temperature, a feature that dominates the mechanics of most structural systems.

We conduct 365 independent experiments, corresponding to the available temperature values of Fig. 12(a). Simulations are performed using the same parameters as in the LTI case and the same lowpass filtering and resampling procedures, for limiting the analysis in the $[0,50]$ Hz band.

In the above set of equations, $S(·),(·)(T)$ is a set of basis functions, $P(·)(T)∼GP(0,κ(·)(T,T′))$⁠, where $κ(·)(T,T′))$ is a kernel function determined by a set of hyperparameters $ϑ(·)$⁠, and $ϵ(·)∼N(0,σϵ(·)2)$⁠. The associated identification problem involves the estimation of all unknown model parameters of Eqs. (11) and (12), upon proper selection of the basis and kernel functions. For details on the procedure, the reader is referred to Tatsis et al. [154]. Here, we outline the core steps:

• Short time-scale ARX$(n,n,0)$ models are estimated whenever data become available.

• The estimated AR and exogenous parameters, as well as the estimated noise variance and measurement(s) of T are stored.

• The stored parameters are modeled via GPR in accordance to Eq. (12), when a full EOV cycle is completed.

where $ŷ[t]$ denotes the model output, and μ_y is the mean value of the measured output y[t]. The BIC returns n = 42 as the most appropriate order, while the fitness is best for n = 50. However, the corresponding fitness value at the BIC order is quite close to the best one, and thus we select n = 42. Figure 15 illustrates the behavior of the ARX(42, 42, 0) prediction errors, from where the hypothesis of Gaussian white noise process can be safely adopted.

We proceed by estimating ARX(42, 42, 0) models for $tl=1,…,365$ and accordingly storing their parameter vectors and prediction error variances for the GPR modeling stage. As Fig. 16 indicates, the parametric identification process is characterized by high consistency, since the percentage fitness is concentrated around 98.95%.

Following the structural identification stage, we fit GPR models to the available AR, exogenous, and variance data. The selected basis functions and GPR orders are listed in Table 4 of Tatsis et al. [154]. Indicative results are illustrated in Fig. 17 and demonstrate the high fitting to the parameters. It is noted that not all AR/X parameters return structured from the estimation stage, which implies that their GPR modeling may fail. This is an issue already reported by the authors [154]. In any case, the GPR–ARX model may be interpreted as a global description of the system in closed-form, when EOV is present, and may be further applied for SHM purposes [154].

The goal of BMU is to infer uncertain model parameters based on data and to quantify the uncertainty in the estimation. In this section, we demonstrate the BMU process with the aid of the benchmark shear frame structure.

For the sake of simplicity, we assume only one uncertain parameter in the model, specifically the parameter a of the Bouc–Wen model of each link. For this parameter, we assume an uninformative prior uniform distribution, $πpr(a)∼U(0.05,1.0)$⁠. Accelerometers are assumed available at three locations on the frame structure. The input loading on the structure is an assumed Gaussian white noise excitation applied at the base. We generate synthetic acceleration time series data via use of the full-order FE model $MFOM$ for an assumed underlying “true” realization of the parameter a = 0.2 and for slightly perturbed values of the Young's modulus E along the elements of the structure. The acceleration time series dataset is then contaminated with Gaussian white noise of 10% root-mean-square (RMS) noise-to-signal ratio, simulating sensor measurement error. This results in the monitoring dataset $d={z˜k∈RNt,k=1,…,Ndof}$⁠, where $Ndof$ is the number of monitored DOFs, and N_t is the dimension of the measured acceleration time series vector, which depends on the monitoring signal sampling rate. The synthetic dataset d is used in the BMU process.

where $N(· ;0,Σ)$ denotes the value of the multivariate normal probability density function with zero mean vector and covariance matrix $Σ$ evaluated at a specific location. It should be noted that in many applications, particularly when using data obtained from spatially dense sensor grids sampling at high frequencies, spatial and temporal correlation may be present in the discrepancies between measurements and model predictions which can have a significant influence on the estimated posterior distribution. However, accounting for correlation can be challenging due to the difficulty of specifying an appropriate description in the likelihood function, and the additional computational effort required to infer the corresponding parameters [53,54].

The improved transitional MCMC method [157] is used to perform the BMU process. The result is shown in Fig. 18. The posterior distribution has significantly reduced uncertainty compared to the prior distribution and has shifted toward the underlying true value of the parameter α. The posterior uncertainty is non-negligible.

In this section, we utilize the $MFOM$ in order to generate synthetic data of the two-storey frame structure, under earthquake excitation. We use a scaled version of the (1986) Kalamata earthquake record (Kalamata OTE—Building N80E Station) sampled at 100 Hz, applied at a 45 deg angle. The record is scaled at an unrealistically large amplitude in order to excite nonlinear response for this toy problem, which assumes properties that are not necessarily up to the real scale (the assumed nonlinear link stiffness is here assumed $k=1.75×108$ N/m). The numerical analysis was carried out using a time-step $dt=1/100=0.01$ s. The reference simulation hysteretic parameters were set as $α=0.25$⁠, Alpha = 1.0, N = 1, β = 3, and γ = 2, with zero stiffness and strength deterioration effects assumed for all nonlinear joints (links).

In demonstrating the task of virtual sensing, we will here apply a joint state-parameter estimation strategy and attempt to infer both unmeasured signals, such as displacements and velocities, as well as a characteristic system parameter; namely, the hysteretic variable α. We assume that the accelerations in both horizontal directions are measured only for the top storey of this frame. The measured accelerations are produced by contaminating the reference $MFOM$ simulation with a 2% RMS (RMS) noise-to-signal ratio. We additionally, assume that the nonlinear hysteretic parameter α is unmeasured and attempt to jointly identify this on the fly as data is attained (online). In order to achieve this task, as in the former example, we will utilize the $MROM$ as our system model. For this simple example, the computational gains of the ROM are not severe; however, the reduced dimension allows us to deploy a (nonlinear) state space $x(t)$ of reduced dimension (⁠$2×9=18$ displacement and velocity states, versus 348 of the original model), which significantly facilitates the identification task. The UKF is here deployed for the purpose of joint state-parameter estimation. To this end, the uncertain hysteretic parameter is augmented to the original state vector $x(z)$ as an additional state, whose evolution is governed by a random walk assumption. The filter initiates from an off initial guess $α̂0=08$⁠, with assumed process and measurement noise covariance sources defined as $w∼N(0,Q), Q=1×102I18×18$ and $v∼N(0,R), R=1.5×107I12×12$⁠.

Figure 19 illustrates the UKF-estimated response, in the form of displacements and velocities, at characteristic DOFs of the first and second storey of the benchmark frame, namely, nodes o₁ and O₂ as annotated in Fig. 1. The filter is shown to closely approximate these quantities, which are not directly measured, succeeding in the state estimation/virtual sensing task. In addition, Fig. 20 displays the convergence in the estimation of the uncertain hysteretic parameter α. It is worth noting that nonlinearity is not activated in this system until after the third second of the analysis, when the ground motion amplitude is significant enough to induce nonlinear behavior. Due to the nature of its definition (3), the hysteretic parameter α remains unobservable until nonlinearity is excited. This is also manifested in the convergence of the estimation, which is activated only after this parameter becomes observable, i.e., one nonlinearity is induced. Indeed, the tasks of state estimation, parameter estimation, and input estimation require fulfillment of specific formal conditions, namely, observability, identifiability, and invertibility, respectively. These can be defined for linear [158], as well as nonlinear systems [159,160]. It is possible to incorporate the principle of observability within the Bayesian filter construct, as suggested in the discontinuous unscented Kalman filter [161].

This section demonstrates the application of a probabilistic approach for damage identification for the benchmark shear structure. In this approach, determining the existence and severity of damage in the structure is cast as an inference problem where the damage is described by a set of uncertain parameters whose posterior distribution is estimated from data using Bayesian inference. The Young's modulus of all elements in the healthy state before the onset of damage, denoted as $Eh$⁠, is considered a priori known and equal to 210.0 GPa. Damage in the structure is simulated as a uniform reduction of the stiffness of the first and second storey columns, applied at t = 2.5 s. Denoting the Young's modulus of the damaged columns at the ith storey as $Ed,i$⁠, the structure is parameterized by the ratio of the healthy and damaged Young's modulus per storey, yielding the parameters $θ1=Ed,1/Eh$ and $θ2=Ed,2/Eh$ for the first and second storey, respectively.

A synthetic dataset is obtained by subjecting $MFOM$⁠, with ground truth values for θ₁ and θ₂ equal to 0.8, to a Gaussian white noise base excitation. The parameters of the Bouc–Wen hysteretic model are taken as $k=8.5×107$ N/m, $α=0.15$⁠, N = 0.8, β = 3, and γ = 4. The acceleration response is recorded at two DOFs, and the measurements are subsequently contaminated with Gaussian white noise corresponding to a $10%$ RMS noise-to-signal ratio. The recorded acceleration signals are truncated to the interval between t = 2.5 and t = 5.0 s, corresponding to 250 measurements per sensor, yielding the dataset $d={z˜k∈RNt,k=1,…,NDOF}$ with $NDOF=2$ and $Nt=250$ (Fig. 21).

Both θ₁ and θ₂ are assigned a Gaussian prior distribution $N(1.0,0.2)$⁠. The lower tail of the Gaussian prior is truncated at a small positive value in order to restrict the prior and posterior distributions of θ₁ and θ₂ to a positive support. Additionally, c will be considered as an uncertain parameter and inferred alongside the stiffness parameters. The prior distribution $Γ(1.0,0.05)$ is placed over c, where $Γ(kp,sp)$ denotes the Gamma distribution with shape and scale parameters $kp$ and $sp$⁠, respectively. It is noted that in practice, specifying prior distributions over uncertain parameters of the probabilistic model is often challenging or even infeasible due to the limited available information or due to the difficulty of expressing this information in terms of a prior distribution.

The VBMC approach is sample-efficient, making it suitable for computationally expensive likelihood functions. However, the reduction in the number of model evaluations (compared to typical MCMC methods) comes at the cost of additional uncertainty on the posterior distribution of the parameters due to the approximation of the likelihood function and the posterior. The approximate posterior distribution $qϕ(θ1,θ2,c)$ is shown in Fig. 22. It can be seen that additional information in the measurements results in a significant reduction in uncertainty over the distribution of the parameters.

The existence of damage in the structure can be determined directly from samples of the uncertain parameters from the posterior distribution, or from samples of the response from the posterior predictive distribution [12,163]. For this application, the damage factor for the ith storey is defined as $DFi=Eh,i−Ed,i/Eh,i$⁠. Samples from the approximate posterior are used to compute the probability of damage in a given storey exceeding a given value of the damage factor df, i.e., $P[DFi>df]$⁠. The results are shown in Fig. 23. The overestimation of the first storey damage factor compared to the known ground truth can be attributed to the limited amount of data used to perform the damage identification, as well as the combined influence of measurement uncertainty, model form uncertainty, and the uncertainty over the posterior distribution due to the variational approximation.

In this paper, the definition and classification of various sources of uncertainty in structural dynamics, system identification, and SHM domains was provided through the prism of the classical input–output system representation. A taxonomy is proposed for consistent classification of the different sources of uncertainties into: (i) input uncertainty; (ii) model form uncertainty; (iii) model parameter/variable uncertainty; (iv) measurements uncertainty; and (v) inherent variability. We further offer an overview of available strategies for tackling associated uncertainties, with a focus on key downstream tasks within these domains. We demonstrate a select subset of such schemes by means of their implementation on a benchmark shear frame with nonlinear dynamic behavior, induced via hysteretic joints. The running example is deployed to investigate the uncertainty quantification challenges in five key downstream tasks in SHM: (i) model inference; (ii) environmental and operational variability characterization; (iii) model updating; (iv) state estimation/virtual sensing; and (v) inverse problem formulation for damage identification. Through this review, a set of open challenges emerges, including:

• The definition of prior assumptions affecting the uncertainty quantification. Including: (i) inaccurate choice of constitutive behavior (model form and parameter uncertainty), such as the dissipation model form and/or parameterization, or dissipation identification strategy, leading to inaccurate quantification of driving material behavior or energy dissipation mechanisms and their uncertainty. (ii) Invalid simplifying modeling assumptions, such as that of LTI behavior (model form uncertainty), often leading to unidentifiable nonlinearities, damage or degradation, and/or wrong quantification of the uncertainty in the latent model parameters. (iii) Handling of the statistical nature of the measurement uncertainty: e.g., measurements/data values are often assumed as independent and identically distributed, free from unknown external influences and/or biases; inappropriate Gaussian-like errors assumption. (iv) Choice of a data-driven architecture (model form uncertainty), which may not be appropriate for the task to be addressed, leading to underfitting/overfitting and poor generalization to unseen conditions.

• Metrics to assess uncertainty quantification accuracy, benchmarks, and validation strategies. Different SHM algorithms might return similar uncertainty bounds on training data; however, their performance on unseen data (test data) might be very different, as well as their ability to quantify the contribution of each source of uncertainty, to identity reducible and irreducible uncertainty, and to properly quantify uncertainties of slowly/fast time-evolving parameters. There is a need to develop metrics, complex benchmark case studies, and validation strategies that can be used to rate the performance of different uncertainty quantification algorithms on key tasks for SHM applications. We suggest the adoption of the proposed taxonomy for establishing such uncertainty class-conscious metrics.

• Scalable uncertainty quantification solutions for complex problems. Including: (i) mapping the expensive-to-evaluate system model in a low-dimensional space to speed-up computation while retaining accuracy (see suggested use of reduced order modeling schemes); (ii) high-precision learning from small informative datasets or large, heterogeneous, and (spatially or temporally) correlated dataset; (iii) dealing with large number of uncertainties; (iv) automatically reducing the number of uncertainties based on the task at hand; and (v) efficiently quantifying uncertainties in the quantity of interest and their impact on decision-making.

• Strategies for reducing uncertainties. Including: (i) distinguishing and handling confounding influences (including EOV, sensor failure, quality of data, and noise in the data) to quantify and potentially remove otherwise irreducible uncertainties; (ii) scalable solutions for monitoring deployment to reduce uncertainty in the latent parameters by accounting for different measurements type, resolutions, and accuracy; and (iii) development of value of information (VoI) strategies for quantifying the cost/benefit of carrying out additional measurements, information, inspection, maintenance, and so on and improving remaining useful life quantification by combining physics and domain knowledge with data-driven strategies.

By putting forth the suggested taxonomy for classification of uncertainties in the realm of SHM, we aspire to put in place a tool that is useful for organizing the treatment of these sources in a consistent manner, while further serving as a basis for systematically tackling the above-listed open challenges.

Professor Chatzi wishes to acknowledge the contribution of the French-Swiss project MISTERY funded by the French National Research Agency (ANR PRCI Grant No. 266157) and the Swiss National Science Foundation (Grant No. 200021 L 21271). The authors further wish to acknowledge the MSCA Staff Exchanges 2021 project ReCharged (Grant ID: 101086413), and the project LiveQuay: Live Insights for Bridges and Quay walls (Project No. NWA.1431.20.002) of the research programme NWA UrbiQuay which is (partly) funded by the Dutch Research Council (NWO). This project has received funding from the Horizon Europe Programme under the Marie Skłodowska-Curie Staff Exchanges Action (GA No. 101086413). Co-funded by the UK Research and Innovation, and the Swiss State Secretariat for Education, Research and Innovation.

• French National Research Agency (ANR PRCI Grant No. 266157).

• Swiss National Science Foundation (Grant No. 200021 L 21271).

• Dutch Research Council (NWO).

• Horizon Europe Programme under the Marie Skłodowska-Curie Staff Exchanges Action (GA No. 101086413).

• UK Research and Innovation (Funder ID: 10.13039/100014013).

• Swiss State Secretariat for Education, Research and Innovation.

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.