GP-ARX-Based Structural Damage Detection and Localization under Varying Environmental Conditions

Abstract

The representation of structural dynamics in the absence of physics-based models, is often accomplished through the identification of parametric models, such as the autoregressive with exogenous inputs, e.g., ARX models. When the structure is amenable to environmental variations, parameter-varying extensions of the original ARX model can be implemented, allowing for tracking of the operational variability. Yet, the latter occurs in sufficiently longer time-scales (days, weeks, months), as compared to system dynamics. For inferring a “global”, long time-scale varying ARX model, data from a full operational cycle has to typically become available. In addition, when the sensor network comprises multiple nodes, the identification of long time-scale varying, vector ARX models grow in complexity. We address these issues by proposing a distributed framework for structural identification, damage detection and localization. Its main features are: (i) the individual estimation of local, single-input-single-output ARX models at every operational point; (ii) the long time-scale representation of each individual ARX coefficient via a Gaussian process regression, which captures dependency on varying Environmental and Operational Conditions (EOCs); (iii) the establishment of a distributed residual generation algorithm for damage detection, which produces time-series of well-defined stationary statistics, with detected discrepancies used for damage diagnosis; and, (iv) exploitation of ARX-inferred mode shape curvatures, obtained via ARX-inferred global state-space models, of the healthy and damaged states, for damage localization. The method is assessed via application on two numerical case studies of different complexity, with the results confirming its efficacy for diagnostics under varying EOCs.

1. Introduction

In the field of condition monitoring, model-based fault diagnosis (FD) [1,2,3] has been gaining ground, as a robust means for condition assessment. In contrast to other schemes, which are heavily dependent on dense instrumentation systems, the existence of a model in analytical and/or numerical form introduces inter-dependencies between the monitored quantities and, in this way, guarantees a form of redundancy [4]. When the employed model draws from first principles, e.g., a finite element representation, then this offers the added benefit of response estimation in unmeasured locations. The latter is also referred to as “virtual sensing" in recent literature [5,6]. For instance, when structural, e.g., finite element, models are coupled with observers, such as those that are based on Kalman filtering, the notion of virtual sensing can be exploited in strain estimation for fatigue assessment [7,8,9,10,11].

However, it is often non-trivial—or even practically infeasible—to assemble a precise system model that relies on first principles. Detailed representations using finite/discrete element models or multibody representations tend to suffer from the curse of uncertainty, tied to the required definition of modeling parameters, and will typically come with considerable computational toll. This is particularly true in the modeling of wind turbines components, where, due to the complexity of the involved geometry and materials, the establishment of a structural model forms an intricate task [12]. For such components in particular, FD plays an important role, as it could serve for damage prediction and preventive maintenance. For achieving diagnosis at reduced computation, it is, in this case, desirable to avoid the utilization of physics-based models. This work exploits availability of a parametric model structure, which is purely derived from availability of data, in the form of input/output, i.e., load/response information.

FD studies have been conducted for a diverse suite of systems including civil [13,14,15], mechanical [16,17], wind energy [18,19,20,21,22,23], and aerospace structures [24,25,26,27,28]. An important challenge in estimation originates from the exposure of such engineered systems, even if of diverse functionality, to varying Environmental and Operational Conditions (EOCs) [29,30,31]. Most engineered systems, and structures in particular, are subject to continually changing environments, e.g., temperature and humidity conditions. Such benign variability effects often dominate the monitored quantities, often masking the changes induced by actual damage or deterioration [32,33,34]. In order to account for this variability in a way that renders detection of actual damage robust, output-only or input-output schemes may be followed. In the first case, detection relies exclusively on measurements of structural response, whilst in the second case inference of a model is sought between the system output (response) and the measured EOCs. Following an output-only scheme, Bernal [35] discusses damage detection and localization in civil engineering structures in absence of information on the environment. Figueiredo et al. [36] examine different machine learning alternatives; Harmanci et al. [37] exploit an output-only autoregressive approach; Sohn et al. [38] and Dervilis et al. [39] successfully test the use of auto-associative neural networks; while, Lämsä & Kullaa [40,41] rely on factor analysis. Numerous works in the ouput-only front rely on principle component analysis (PCA) for the incorporation of EOCs [42]. Applications in FD for structures include damage identification in the form of both detection [43,44,45] and localization [46,47].

In this work, we propose an input–output scheme for tackling dependence of structural properties on EOC variability. Regression often forms a main tool to this end, targeting the construction of an input/output dependence function, which can be then complemented with a probabilistic treatment of the estimation error for robust outlier detection [48,49]. Alternatively, time-series representations can be adopted as models of the underlying dynamics, while explicitly incorporating dependence on influencing factors, such as EOCs. Different formulations exist to this end. Linear Parameter Varying (LPV) models utilize a two-stage approach for capturing this dependence. Firstly, local (or frozen) models are computed, which correspond to specific EOC configurations, and secondly, the parameters of the identified models are interpolated to provide a single global model [50,51], operating across EOCs. Another option lies in adoption of Functionally Pooled (FP) time-series models, which are further described in Kopsaftopoulos et al. Kopsaftopoulos et al. [52] and Sakellariou and Fassois [53]. Here, the AR coefficients are modeled as explicit functions of the EOC vector, which could for instance pertain to airspeed and angle of attack for parameters met in aerospace applications. A functional basis is used comprising bivariate polynomials, resulting as tensor products from corresponding univariate polynomials (Chebyshev, Legendre, Jacobi, and other families). Spiridonakos et al. [54] and Spiridonakos and Chatzi [55] have introduced a Polynomial Chaos basis instead in order to span the functional subspace expressing the dependence of the AR models on the EOCs [56,57] tackling variability in structural components. A further alternative pertains to use of Random Coefficient (RC) models, where the time–dependent coefficients follow a multivariate Gaussian mixture model (GMM), allowing for significant flexibility in uncertainty representation [58,59]. To further separate effects, which operate on different temporal scales, Avendaño-Valencia et al. [51] propose a hierarchical time-series model. This adopts an LPV-AR representation for capturing influences on short-term dynamics, e.g., azimuth of the rotor for operating wind turbines, while long-term influences stemming from effects such as temperature, are tackled by means of Gaussian process regression (GPR) [60]. The framework is demonstrated on the tracking of the response of the Humber bridge under long-term influence that is induced by varying acting averaged wind speeds.

Based on the above, a two-stage algorithmic process is proposed in this work, with the aim (i) to construct a data-driven model of the structure, able to reliably reconstruct the response across the full range of operation of the system in its healthy state and (ii) reliably detect and possibly localize damage, via use of a suitably defined Residual Generation Algorithm (RGA) [1,61,62,63,64] featuring appropriately selected detection thresholds. The first step requires training across a sufficient number of structural response samples that sufficiently span the "global" system response. A GP-ARX model is then constructed, which is valid across the whole range of operating conditions of the structure. Once the model is put in place, a RGA is formulated as the difference between the measured response and the noise-free part of the GP-ARX model. It is demonstrated that the RGA follows a GP-AR model: this allows for delivering the theoretical statistical thresholds and defining the range of values within which the ”healthy” state lies, as a function of the EOC variable vector. Any systematic deviation from the predefined threshold zone is then a safe alert for the presence of damage. It is shown that damage detection can be accomplished even by a single GP-ARX model, e.g., a single response measurement sensor. Localization, which inevitably depends on the available sensor network, proceeds by calculating the differences between the mode shape curvatures (MSCs) of the healthy and the damaged structure.

The paper is organized, as follows: Section 2 defines the problem and delivers the theoretical difference equation between a single input-output pair, including the dependence of the input-output delay with respect to the structural response type monitored. Section 3 describes the structural identification stage, formulates the RGA and establishes the damage detection and localization procedure. The method is assessed in Section 4, via two numerical investigations of different complexity, namely a spring-mass-damper system and a shear frame reflecting a common representation of a structural system. Finally in Section 5, the results are summarized and directions for further research are given.

2. Problem Formulation

A linear, viscously damped, operationally varying structural system with n degrees of freedom (DOFs) may be represented by a second–order vector differential equation of the form

$$\begin{array}{c}\hfill \mathbf{M}\left(\mathbf{\xi }\right)\ddot{\mathbf{x}}\left(t\right)+\mathbf{C}\left(\mathbf{\xi }\right)\dot{\mathbf{x}}\left(t\right)+\mathbf{K}\left(\mathbf{\xi }\right)\mathbf{x}\left(t\right)=\mathbf{u}\left(t\right)\end{array}$$

(1)

where $\mathbf{x}\left(t\right)$ is the $[n\times 1]$ vibration displacement vector, $\mathbf{M}\left(\mathbf{\xi }\right)$, $\mathbf{C}\left(\mathbf{\xi }\right)$ and $\mathbf{K}\left(\mathbf{\xi }\right)$ are the real $[n\times n]$ mass, viscous damping and stiffness matrices, respectively, and $\mathbf{u}\left(t\right)$ is the $[n\times 1]$ vector of excitations.

The operational variability of the structural matrices may be attributed to diverse environmental parameters (temperature, humidity, etc.), which are herein represented by a $[L\times 1]$ vector $\mathbf{\xi }={[{\mathbf{\xi }}_1,{\mathbf{\xi }}_2,\dots ,{\mathbf{\xi }}_L]}^T$. It is assumed that their temporal evolution is sufficiently slow as compared to the variation of the structural dynamics (e.g., sufficiently lower with respect to the lowest structural vibration mode). Subsequently, Equation (1) can be considered as a local LTI model that describes the structure during a small time interval. The two different time-scales are henceforth referred to as the long time-scale, corresponding to the temporal evolution of $\mathbf{\xi }$ and notated by $t_l$, and the short time-scale, corresponding, as usual, to the temporal evolution of $\mathbf{x}\left(t\right)$ and notated by t.

Under this setting, within the short time-scale, the parameter vector $\mathbf{\xi }$ may be considered as constant, which is, as a realization of a multivariate random process $\Xi$ with joint PDF $f_{\Xi }\left(\mathbf{\xi }\right)$. Accordingly, Equation (1) may be written as

$$\begin{array}{c}\hfill \mathbf{M}(\Xi )\ddot{\mathbf{x}}\left(t\right)+\mathbf{C}(\Xi )\dot{\mathbf{x}}\left(t\right)+\mathbf{K}(\Xi )\mathbf{x}\left(t\right)=\mathbf{u}\left(t\right),\mathrm{with}\Xi \left(\omega \right)\in \mathsf{\Omega }\end{array}$$

(2)

where $\mathsf{\Omega }$ denotes the event space of the random vector $\Xi$. The state-space representation of Equation (2) reads

$$\begin{array}{cc}\hfill \dot{\mathit{\zeta }}\left(t\right)& ={\mathbf{A}}_c(\Xi )\mathit{\zeta }\left(t\right)+{\mathbf{B}}_c(\Xi )\mathbf{u}\left(t\right)\hfill \end{array}$$

(3a)

$$\begin{array}{cc}\hfill \mathbf{y}\left(t\right)& ={\mathbf{H}}_c(\Xi )\mathit{\zeta }\left(t\right)+{\mathbf{D}}_c(\Xi )\mathbf{u}\left(t\right)\hfill \end{array}$$

(3b)

in which the state and input matrices have the standard form

$${\mathbf{A}}_c(\Xi )=\left[\begin{array}{cc}{\mathbf{O}}_n& {\mathbf{I}}_n\\ -{\mathbf{M}}^{-1}(\Xi )\mathbf{K}(\Xi )& -{\mathbf{M}}^{-1}(\Xi )\mathbf{C}(\Xi )\end{array}\right],{\mathbf{B}}_c(\Xi )=\left[\begin{array}{c}\mathbf{0}\\ {\mathbf{M}}^{-1}(\Xi )\end{array}\right]$$

(4)

and the output and feedforward matrices, ${\mathbf{H}}_c(\Xi )$ and ${\mathbf{D}}_c(\Xi )$, respectively, depend on the type of the monitored structural response (vibration displacement, velocity, acceleration, strain, etc.) $y\left(t\right)$, of size $[s\times 1]$.

A representation of Equation (3) in discrete-time proceeds by adopting an appropriate sampling period $T_s$, for which the intersample behaviour of the input can be regarded as constant (e.g., zero-order hold assumption). Subsequently

$$\begin{array}{cc}\hfill \mathit{\zeta }[k+1]& =\mathbf{A}(\Xi )\mathit{\zeta }\left[k\right]+\mathbf{B}(\Xi )\mathbf{u}\left[k\right]\hfill \end{array}$$

(5a)

$$\begin{array}{cc}\hfill \mathbf{y}\left[k\right]& =\mathbf{H}(\Xi )\mathit{\zeta }\left[k\right]+\mathbf{D}(\Xi )\mathbf{u}\left[k\right]\hfill \end{array}$$

(5b)

for

$$\mathbf{A}(\Xi )=e^{{\mathbf{A}}_c(\Xi )T_s},\mathbf{B}(\Xi )=\left[\mathbf{A}(\Xi )-{\mathbf{I}}_n\right]{\mathbf{A}}_c^{-1}(\Xi ){\mathbf{B}}_c(\Xi ),\mathbf{H}(\Xi )={\mathbf{H}}_c(\Xi ),\mathbf{D}(\Xi )={\mathbf{D}}_c(\Xi )$$

(6)

Equation (5) can be transformed into a modal form as

$$\begin{array}{cc}\hfill \mathbf{\eta }[k+1]& =\mathbf{\Lambda }(\Xi )\mathbf{\eta }\left[k\right]+{\mathbf{\Pi }}^T(\Xi )\mathbf{u}\left[k\right]\hfill \end{array}$$

(7a)

$$\begin{array}{cc}\hfill \mathbf{y}\left[k\right]& =\mathbf{\Psi }(\Xi )\mathbf{\eta }\left[k\right]+\mathbf{D}(\Xi )\mathbf{u}\left[k\right]\hfill \end{array}$$

(7b)

with $\mathbf{\Lambda }(\Xi )$ corresponding to the diagonal matrix of the eigenvalues of $\mathbf{A}(\Xi )$, ${\mathbf{\Pi }}^T(\Xi )$ to the matrix of modal participation vectors and ${\mathbf{\Psi }}^T(\Xi )$ to the matrix of mode shapes. It is reminded that, for an underdamped structure, the diagonal entries of $\mathbf{\Lambda }(\Xi )$ arrive in complex conjugate pairs, and so do the row vectors of ${\mathbf{\Pi }}^T(\Xi )$ and the column vectors of ${\mathbf{\Psi }}^T(\Xi )$.

For convenience, we will henceforth assume that the structure is excited by a single input, e.g., $\mathbf{u}\left[k\right]\to u\left[k\right]$. Subsequently, from Equation (7), the dynamics of a single input-output pair are described by

$$\begin{array}{cc}\hfill \mathbf{\eta }[k+1]& =\mathbf{\Lambda }(\Xi )\mathbf{\eta }\left[k\right]+{\mathbf{\Pi }}^T(\Xi )u\left[k\right]\hfill \end{array}$$

(8a)

$$\begin{array}{cc}\hfill y_{\ell }\left[k\right]& ={\mathbf{\psi }}_{\ell }(\Xi )\mathbf{\eta }\left[k\right]+d_{\ell }(\Xi )u\left[k\right]\hfill \end{array}$$

(8b)

where ${\mathbf{\psi }}_{\ell }(\Xi )$ and $d_{\ell }$ are the ℓ-th row and element of $\mathbf{\Psi }(\Xi )$ and $\mathbf{D}(\Xi )$, respectively. The application of the $\mathcal{Z}$-transform to Equation (8) yields

$$\begin{array}{cc}\hfill \mathbf{\eta }\left(z\right)& ={\left[z{\mathbf{I}}_{2n}-\mathbf{\Lambda }(\Xi )\right]}^{-1}{\mathbf{\Pi }}^T(\Xi )u\left(z\right)\hfill \end{array}$$

(9a)

$$\begin{array}{cc}\hfill y\left(z\right)& ={\mathbf{\psi }}_{\ell }(\Xi )\mathbf{\eta }\left(z\right)+d_{\ell }(\Xi )u\left(z\right)\hfill \end{array}$$

(9b)

and the transfer function between $y\left(z\right)$ and $u\left(z\right)$ is

$$G_{\ell }\left(z\right)={\mathbf{\psi }}_{\ell }(\Xi ){\left[z{\mathbf{I}}_{2n}-\mathbf{\Lambda }(\Xi )\right]}^{-1}{\mathbf{\Pi }}^T(\Xi )+d_{\ell }(\Xi )$$

(10)

or, as a modal decomposition

$$G_{\ell }\left(z\right)=\sum _{i=1}^{2n}\frac{{\mathbf{\psi }}_{\ell ,i}(\Xi ){\mathit{\pi }}_i^T(\Xi )}{z-{\lambda }_i(\Xi )}+d_{\ell }(\Xi )=\sum _{i=1}^{2n}\frac{R_{\ell ,i}(\Xi )z^{-1}}{1-{\lambda }_i(\Xi )z^{-1}}+d_{\ell }(\Xi )$$

(11)

with $R_{\ell ,i}(\Xi )$ obviously defined. In deriving a rational expression for $G\left(z\right)$, notice that the far right hand side of Equation (11) can be written as

$$G_{\ell }\left(z\right)=\sum _{i=1}^{2n}\frac{R_{\ell ,i}(\Xi )z^{-1}{\overline{A}}_i\left(z\right)}{A\left(z\right)}+d_{\ell }(\Xi )$$

(12)

for

$$\begin{array}{cc}\hfill A\left(z\right)& =(1-{\lambda }_1(\Xi )z^{-1})(1-{\lambda }_1(\Xi )z^{-1})\dots (1-{\lambda }_{2n}(\Xi )z^{-1})\hfill \\ \hfill & =1+a_1(\Xi )z^{-1}+a_2(\Xi )z^{-2}+\dots +a_{2n}(\Xi )z^{-2n}\hfill \end{array}$$

(13)

The polynomial ${\overline{A}}_i\left(z\right)$ in the numerator of Equation (12) corresponds to $A\left(z\right)$ with the i-th eigenvalue removed and it reaches up to order $-2n+1$. By combining further (i) the presence of the unit delay, $z^{-1}$; and, (ii) the term $d_{\ell }(\Xi )$, a qualitative rational expression for $G\left(z\right)$ can be derived

$$G_{\ell }\left(z\right)=\frac{B\left(z\right)}{A\left(z\right)}=\frac{b_0(\Xi )+b_1(\Xi )z^{-1}+b_2(\Xi )z^{-2}+\dots +b_{2n}(\Xi )z^{-2n}}{1+a_1(\Xi )z^{-1}+a_2(\Xi )z^{-2}+\dots +a_{2n}(\Xi )z^{-2n}}$$

(14)

with $b_0(\Xi )=d_{\ell }(\Xi )$. Based on Equations (10)–(13), it follows that the coefficients of $B\left(z\right)$ depend on both of the eigenvalues and the eigenvectors of $\mathbf{A}(\Xi )$. Also observe that, in the absence of a feedforward term in Equation (8), then $G\left(z\right)$ has a unit input-output delay, since $d_{\ell }(\Xi )=0$. In either case, the difference equation that describes the structural dynamics of the ℓ-th response in discrete-time stems from Equation (14) as

$$\begin{array}{cc}\hfill y_{\ell }\left[k\right]+a_1(\Xi )y_{\ell }[k-1]+\dots +& a_{2n}(\Xi )y_{\ell }[k-2n]\hfill \\ \hfill & =b_0(\Xi )u\left[k\right]+b_1(\Xi )u[k-1]+\dots +b_{2n}(\Xi )u[k-2n]\hfill \end{array}$$

(15)

Under this setting, the problem treated herein is the detection and the localization of damage when

-

there is no prior information about the structure,

-

input-output measurements are consistently recorded over the short scale (e.g., k), and

-

the complete set of $\Xi$ is available (measurable) over the long time-scale (e.g., $t_l$)

The proposed strategy is of a decentralized nature, in the sense that estimation is first accomplished on a local basis, by considering all individual input-output data pairs independently, as dictated by Equation (15). The adopted framework comprises two phases, e.g., the structural identification phase and the damage detection and localization phase, which are discussed in the following.

3. The Methodological Framework

3.1. The Structural Identification Phase

For simplicity, it is henceforth assumed that the operational variability is due to a single parameter, e.g., $L=1$ and $\Xi \to \Xi$. Equation (15) indicates that at the healthy state of the structure, an excitation-response pair can be effectively described in discrete-time by an ARX model of the form [65]

$$y\left[k\right]+\sum _{i=1}^p{a}_i(\Xi )y[k-i]=\sum _{i=0}^p{b}_i(\Xi )u[k-i]+e\left[k\right]$$

(16)

where $y\left[k\right]$ denotes the noise-corrupted structural vibration response of a single DOF, $u\left[k\right]$ is the excitation, $a_i(\Xi )$ and $b_i(\Xi )$ are the coefficients of the AR and exogenous polynomials, of order p, and $e\left[k\right]\sim \mathcal{N}(0,{\sigma }_e^2(\Xi ))$ is a zero-mean Gaussian white noise realization that models process and measurement uncertainty (e.g., model inconsistencies and instrumentation noise). It is noted that, in Equation (16), $b_0\ne 0$ only when vibration acceleration response measurements are available.

Equation (16) indicates that the ARX($p,p$) model accounts for the variability of the response due to variations of $\Xi$ through the AR and exogenous parameters, ${\alpha }_i(\Xi )$ and $b_i(\Xi )$, respectively, as well as through the noise variance matrix ${\sigma }_e^2(\Xi )$. These are herein admitting a GPR representation as

$$\begin{array}{cc}\hfill a_i(\Xi )& =\sum _{j=1}^{p_a}{\alpha }_{ij}{S}_{a,j}(\Xi )+P_a(\Xi )+{ϵ}_a\hfill \end{array}$$

(17a)

$$\begin{array}{cc}\hfill b_i(\Xi )& =\sum _{j=1}^{p_b}{\beta }_{ij}{S}_{b,j}(\Xi )+P_b(\Xi )+{ϵ}_b\hfill \end{array}$$

(17b)

$$\begin{array}{cc}\hfill {\sigma }_e^2(\Xi )& =\sum _{j=1}^{p_{\sigma }}{\gamma }_{ij}{S}_{\sigma ,j}(\Xi )+P_{\sigma }(\Xi )+{ϵ}_{\sigma }\hfill \end{array}$$

(17c)

In light of this parametric representation, the model of Equation (16) is termed as GP-ARX($p,p$). In Equation (17), $S_{(·),(·)}(\Xi )$ denotes an appropriate set of basis functions, $P_{(·)}(\Xi )\sim \mathcal{G}\mathcal{P}(0,{\kappa }_{(·)}(\Xi ,{\Xi }^{\prime }))$, where ${\kappa }_{(·)}(\Xi ,{\Xi }^{\prime }))$ is the kernel function, being fully determined by a set of hyperparameters ${\mathbf{\vartheta }}_{(·)}$, and ${ϵ}_{(·)}\sim \mathcal{N}(0,{\sigma }_{{ϵ}_{(·)}}^2)$

Upon proper selection of the basis and kernel functions, $S_{(·),(·)}(\Xi )$ and ${\kappa }_{(·)}(\Xi ,{\Xi }^{\prime }))$, respectively, the representation of the structure through GP-ARX($p,p$) models is succeeded through the solution of an associated identification problem, in which all unknown model parameters are estimated, which is, the integers $p,p_a,p_b,\mathrm{and}{p}_{\sigma }$, the GPR coefficients ${\alpha }_{ij}$, ${\beta }_{ij}$ and ${\gamma }_{ij}$ in Equation (17), as well as the vectors of hyperparameters ${\mathbf{\vartheta }}_a$, ${\mathbf{\vartheta }}_b$ and ${\mathbf{\vartheta }}_{\sigma }$ and the variances ${\sigma }_{{ϵ}_a}^2$, ${\sigma }_{{ϵ}_b}^2$, and ${\sigma }_{{ϵ}_{\sigma }}^2$. Yet, since the operational parameter $\Xi$ varies slowly, this problem cannot be treated by pooling a large regression problem over a set of Q realizations of $\Xi$ [65,66]: in applications involving large infrastructures (bridges, buildings, wind turbines, etc.), $\Xi$ usually varies seasonally. Thus, it is rather impractical to wait until a full operational cycle becomes available before estimating a global GP-ARX model.

In contrast, the strategy chosen herein pertains to (i) the estimation of short time-scale ARX($p,p$) models, whenever data become available; (ii) the storage of the estimated AR and exogenous parameters, the estimated variance and measurement(s) of $\Xi$; and, (iii) the GP modelling of all parameters according to Equation (17), when a full operational cycle is completed. Notice that the direct modeling of the parameters $a_i(\Xi )$, $b_i(\Xi )$, and ${\sigma }_e^2(\Xi )$ in Equation (17), is a process that can be easily realized in parallel and essentially corresponds to typical data-driven uncertainty quantification problem.

The short time-scale ARX modelling step is conventionally accomplished: for any given pair of input-output data at a specific long time-scale instant $t_l$, $u_{t_l}\left[k\right]$ and $y_{t_l}\left[k\right]$, respectively, for $k=1,2,\dots ,N$, the original ARX model is transformed into its linear regression form

$$y_{t_l}\left[k\right]={\mathbf{\varphi }}_{t_l}{\left[k\right]}^T{\mathbf{\theta }}_{t_l}+e_{t_l}\left[k\right]$$

(18)

for

$${\mathbf{\varphi }}_{t_l}\left[k\right]={\left[-y_{t_l}[k-1]\dots -y_{t_l}[k-p]|u_{t_l}\left[k\right]\dots u_{t_l}[k-p]\right]}^T$$

(19)

and

$${\mathbf{\theta }}_{t_l}={[{\alpha }_{1,t_l}{\alpha }_{2,t_l}\dots {\alpha }_{p,t_l}|b_{0,t_l}{b}_{1,t_l}\dots b_{p,t_l}]}^T$$

(20)

For $k=1,\dots ,N$, Equation (18) implies

$${\mathbf{Y}}_{t_l}={\mathbf{\Phi }}_{t_l}{\mathbf{\theta }}_{t_l}+{\mathbf{E}}_{t_l}$$

(21)

where ${\mathbf{Y}}_{t_l}={\left[y_{t_l}\left[1\right]\dots y_{t_l}\left[N\right]\right]}^T$, ${\mathbf{\Phi }}_{t_l}={\left[{\mathbf{\varphi }}_{t_l}^T\left[1\right]\dots {\mathbf{\varphi }}_{t_l}^T\left[N\right]\right]}^T$ and ${\mathbf{E}}_{t_l}={\left[e_{t_l}\left[1\right]\dots e_{t_l}\left[N\right]\right]}^T$. Minimization of the quadratic least squares criterion provides the estimate

$${\widehat{\mathbf{\theta }}}_{t_l}={\left({\mathbf{\Phi }}_{t_l}^T{\mathbf{\Phi }}_{t_l}\right)}^{-1}{\mathbf{\Phi }}_{t_l}^T{\mathbf{Y}}_{t_l}$$

(22)

while the prediction errors are recovered by ${\mathbf{E}}_{t_l}={\mathbf{Y}}_{t_l}-{\mathbf{\Phi }}_{t_l}{\widehat{\mathbf{\theta }}}_{t_l}$ and their variance is calculated as ${\sigma }_{e,t_l}^2={\mathbf{E}}_{t_l}^T{\mathbf{E}}_{t_l}$.

This procedure is repeated over the whole range of available operational realizations ${\Xi }_{t_l}$, $t_l=1,\dots Q$. Under the reasonable assumption that the order of the structural system remains unaltered, the AR and exogenous polynomial orders, estimated statistically from an arbitrary single realization, can be maintained the same for every $t_l$. Subsequently, the GPR estimation of each individual parameter of Equation (17) is taking place while using the hybrid optimization method described in Dertimanis et al. [67].

3.2. The Damage Detection and Localization Phase

The successful completion of the previous phase provides the availability of s single-input, single-output GP-ARX${}_{\ell }(p,p)$ models, for $\ell =1,\dots ,s$, where s corresponds to the number of available DOFs’ measurements. These models are capable of tracking the evolution of the structural response of the relative DOF during the healthy structural state, for the whole expected range of variations in $\Xi$. In establishing an effective damage detection scheme, a residual generation algorithm (RGA) is formulated as shown in Figure 1, for every AFS-ARX model. The RGA is a computational machine that operates in parallel to the actual structure and issues an alert when abnormal behavior appears, i.e., different to the one estimated during the structural identification phase.

In order to formulate the RGA, the ℓ-th AFS-ARX model of Equation (16) is compactly written as

$$y_{\ell }\left[k\right]=\frac{B_{\ell }(q,\Xi )}{A(q,\Xi )}u\left[k\right]+\frac{1}{A(q,\Xi )}{e}_{\ell }\left[k\right]$$

(23)

where q is the backshift operator (i.e., $q^{-i}y\left[k\right]=y[k-i]$) and $A(q,\Xi )$, $B_p(q,\Xi )$ the AR end exogenous polynomials, respectively, of the form,

$$\begin{array}{cc}\hfill A(q,\Xi )& =1+a_1(\Xi )q^{-1}+\dots +a_p(\Xi )q^{-p}\hfill \end{array}$$

(24a)

$$\begin{array}{cc}\hfill B_{\ell }(q,\Xi )& =b_{0,\ell }(\Xi )+b_{1,\ell }(\Xi )q^{-1}+\dots +b_{p,\ell }(\Xi )q^{-p}\hfill \end{array}$$

(24b)

Notice that the first term of the right-hand side of Equation (23) corresponds to the estimated discrete-time transfer function between the applied excitation and the vibration response at the monitored DOF of the structure. To this end, the RGA can be simply formulated as the difference

$$r_{\ell }\left[k\right]=y_{\ell }\left[k\right]-\frac{B_{\ell }(q,\Xi )}{A(q,\Xi )}u\left[k\right]$$

(25)

retaining the units of the monitored quantity. In the absence of damage, the residual $r_{\ell }\left[k\right]$ is thus equal to

$$r_{\ell }\left[k\right]=\frac{1}{A(q,\Xi )}{e}_{\ell }\left[k\right]$$

(26)

implying that it follows an GP-AR$\left(p\right)$ model of the form

$$r_{\ell }\left[k\right]+\sum _{i=1}^p{a}_i(\Xi )r_{\ell }[k-i]=e_{\ell }\left[k\right]$$

(27)

The variance of $r_{\ell }\left[k\right]$ can be calculated by truncating the last expression to an infinite moving average (MA) model

$$r_{\ell }\left[k\right]=e_{\ell }\left[k\right]+\sum _{i=1}^{\infty }{\mathcal{H}}_{i,\ell }{e}_{\ell }[k-i]$$

(28)

in which ${\mathcal{H}}_{i,\ell }$ are the coefficients of the associated Green function [68]. The mean and the variance of the time-series of Equation (28) are given by

$${\mu }_{r_{\ell }}=E\left\{r_{\ell }\left[k\right]\right\}=0\mathrm{and}{\sigma }_{r_{\ell }}^2(\Xi )=E\left\{r_{\ell }^2\left[k\right]\right\}=\sum _{i=0}^{\infty }{\mathcal{H}}_{i,\ell }^2{\sigma }_{e,\ell }^2(\Xi )$$

(29)

respectively. This means that 99.7% of the damage-free RGA time-series conform to

$$-3{\sigma }_{r_{\ell }}(\Xi )\le r_{\ell }\left[k\right]\le 3{\sigma }_{r_{\ell }}(\Xi )$$

(30)

Equations (28)–(30) imply that the quantitative evaluation of the damage detection thresholds require knowledge of the infinite MA model coefficients ${\mathcal{H}}_{i,\ell }$. These are calculated in closed form via the residue expansion of the transfer function of Equation (27). The reader is referred to Box et al. [68] [Chapter 3] for further details. Here, it is simply reminded that the AR polynomial should be of minimum-phase for the whole range of operating conditions (e.g., for all values of the temperature), in order to ensure the convergence of ${\mathcal{H}}_{i,p}$, and that AR roots’ amplitudes close to unity result in a slower the decay rate of ${\mathcal{H}}_{i,p}$.

Damage localization follows the detection stage and operates on data that are safely considered as outcomes of the damaged structure, e.g., after the damage time instant detected by the RGAs. Localization requires the estimation of a new set of simple ARX models, one for each available input–output data set and it is based on the differences between the mode shape curvatures between the GP-ARX models and the estimated-based-on-damaged-data ARX ones. For both sets of models, mode shape curvature estimation proceeds by first establishing the state-space realizations of the noise-free parts of the GP-ARX and ARX models [69]

$$\begin{array}{cc}\hfill \mathit{\zeta }[t+1]& =\mathbf{A}\mathit{\zeta }\left[t\right]+\mathbf{B}u\left[t\right]\hfill \end{array}$$

(31a)

$$\begin{array}{cc}\hfill \mathbf{y}\left[t\right]& =\mathbf{H}\mathit{\zeta }\left[t\right]+\mathbf{D}u\left[t\right]\hfill \end{array}$$

(31b)

where

$$\begin{array}{cc}\hfill \mathbf{A}=& \left[\begin{array}{cccc}0& 1& \dots & 0\\ 0& 0& \dots & 0\\ ⋮& ⋮& \ddots & 0\\ 0& 0& 0& 1\\ -a_p& -a_{p-1}& \dots & -a_1\end{array}\right],\mathbf{B}=\left[\begin{array}{c}0\\ 0\\ ⋮\\ 0\\ 1\end{array}\right]\hfill \end{array}$$

(32a)

$$\begin{array}{cc}\hfill \mathbf{H}=& \left[\begin{array}{cccc}{b}_p-a_p{b}_0& b_{p-1}-a_{p-1}{b}_0& \dots & b_1-a_1{b}_0\end{array}\right],\mathbf{D}=\left[\begin{array}{c}{b}_0\end{array}\right]\hfill \end{array}$$

(32b)

and then transforming Equation (31) to the modal form of Equation (8) via the eigenvalue problem of the state matrix $\mathbf{A}$. The real part of the estimated row vector ${\mathbf{\psi }}_{\ell }$ (see Equation (8)) contains the values of the identified mode shapes at the ℓ-th sensor location. Gathering together all of the estimated ${\mathbf{\psi }}_{\ell }$ vectors, the matrix of mode shapes is constructed as

$$\mathbf{\Psi }=\left[\begin{array}{c}{\mathbf{\psi }}_1\\ {\mathbf{\psi }}_2\\ ⋮\\ {\mathbf{\psi }}_s\end{array}\right]=\left[\begin{array}{cccc}{\tilde{\mathbf{\psi }}}_1& {\tilde{\mathbf{\psi }}}_2& ⋮& {\tilde{\mathbf{\psi }}}_p\end{array}\right]\left([s\times p]\right)$$

(33)

and its columns correspond to the shapes of each vibration mode. In establishing $\mathbf{\Psi }$, attention should be paid to the proper distinction of structural from “erroneous” modes, as well as to the matching of eigenvectors associated to the same mode. This is accomplished by applying dispersion analysis [70,71,72] and formulating the modal assurance criterion (MAC).

When $\mathbf{\Psi }$ is set up to contain the structural vibration modes from all of the identified ARX models, the entries of the ℓ-th column of the mode shape curvature matrix are calculated by [73]

$${\left[{\tilde{\mathbf{\psi }}}_{\ell }^{″}\right]}_i=\frac{{\left[{\tilde{\mathbf{\psi }}}_{\ell }\right]}_{i+1}-2{\left[{\tilde{\mathbf{\psi }}}_{\ell }\right]}_i+{\left[{\tilde{\mathbf{\psi }}}_{\ell }\right]}_{i-1}}{h^2}$$

(34)

with h denoting the distance between two successive measurement nodes. Finally, damage localization is examined by the difference

$$\Delta {\tilde{\mathbf{\psi }}}_{\ell }^{″}={\tilde{\mathbf{\psi }}}_{\ell ,H}^{″}-{\tilde{\mathbf{\psi }}}_{\ell ,D}^{″}$$

(35)

where ${\tilde{\mathbf{\psi }}}_{\ell ,H}^{″}$ and ${\tilde{\mathbf{\psi }}}_{\ell ,D}^{″}$ correspond to the ℓ-th MSC matrix columns established by the GP-ARX models and the estimated-based-on-damaged-data ARX ones, respectively.

4. Case Studies

4.1. Damage Detection on a Spring-Mass-Damper System

A simple structural system with six DOFs, as illustrated in Figure 2, is chosen as the first case study. It is assumed that the stiffness is exponentially dependent on temperature as $K_i\left[t_l\right]={10}^{5-3K_T(T\left[t_l\right]-T_{ref})}$, where $t_l$ corresponds to long (e.g., sufficiently lower with respect to the lowest structural vibration mode) time-scale, $K_T=5\times {10}^{-3}$ is the thermal coefficient, $T_{ref}$ = 15 °C designates a reference temperature, and $T\left[t_l\right]$ a discrete-time stochastic process described as

$$T\left[t_l\right]=T_{ref}+T_{span}sin\left(2\pi f_l{t}_l\right)+w_T\left[t_l\right]$$

(36)

in which $T_{span}$ = 20 °C, $f_l$ is the long time-scale frequency and $w_T\left[t_l\right]\in \mathcal{U}(0,2)$. A stiffness-proportional damping is adopted as $c_i=2\times {10}^{-3}{K}_i$. It should be noted that this is a fictitious equation, which is nonetheless chosen to represent the decrease in stiffness with increasing temperature, met in most civil engineering works. Figure 3 displays a realization of $T\left[t_l\right]$ for $t_l=1,2,\dots ,365$ and $f_l=1/365$, along with its sample and kernel-based probability density function, respectively. For this realization, the induced stiffness and structural vibration modes (natural frequencies and damping ratios) are shown in Figure 4 and Figure 5, respectively.

Under a single excitation acting at DOF #3, all instances of the “healthy” structure are discretized at $T_s=0.02$ s ($F_s=50$ Hz) via the zero-order hold principle. Accordingly, the structural vibration acceleration responses of the discretized instances ${\mathcal{S}}_d\left[t_l\right]$ under independent realizations of zero-mean, Gaussian white noise excitation ($u\left[k\right]\in \mathcal{N}(0,1.87\times {10}^4)$, and data length $N=2^{14}$) are obtained, zero-mean subtracted and noise-corrupted at 1% noise-to-signal-ratio. For the subsequent analysis, it is assumed that only measurements from DOFs #1, #3, and #5 are available. For these, Figure 6 displays the amplitudes of the estimated frequency response functions, along with their theoretical counterparts, while Figure 7 shows the driving point FRF at reference temperature.

The establishment of GP-ARX models for describing the input–output dynamics at measured DOFs proceeds by (i) determining an appropriate order for the AR and exogenous polynomials; (ii) estimating ARX models for the whole long time-scale estimation set, consisting of the first 300 realizations (e.g., $t_l=1,\dots ,300$); (iii) performing GPR modelling for all the associated AR and exogenous polynomials’ coefficients, as well as the variance of the prediction errors; and, (iv) assessing the resulted AFS-ARX models by applying the RGA to the validation set.

Model order selection is accomplished by considering the driving point measurement at a randomly chosen long time instant $t_r\in [1,300]$, e.g., ${\mathcal{S}}_d\left[t_r\right]$. To this end, ARX$(p,p)$ models with zero delay are estimated for $p=10,12,\dots ,50$ and compared while using the Bayesian Information Criterion (BIC) and the Mean Squared Error (RSS), defined as

$$\begin{array}{cc}\hfill BIC& =ln\left({\sigma }_e^2\right)+(2p+1)\frac{ln\left(N_e\right)}{N_e}\hfill \end{array}$$

(37a)

$$\begin{array}{cc}\hfill MSE& =\frac{1}{N_e}\sum _{k=1}^{N_e}{e}^2\left[k\right]\hfill \end{array}$$

(37b)

respectively. In Equation (37), $N_e$ corresponds to the short time-scale estimation set, herein defined as the subset of $N=1,\dots ,\mathrm{16,384}$ that lies in $[1010,\mathrm{15,565}]$.

The performance of the model order selection criteria is depicted in Figure 8a,b. The BIC and the MSE both decay slowly and stabilize at high orders ($p>30$). To elaborate further on the model order selection, Figure 8c shows the frequency stabilization diagrams (FSDs), corrected using the dispersion analysis method described in Dertimanis and Chatzi [71]. The correction is based on the energy associated to each estimated structural vibration mode and on a low threshold (herein set as greater than 1% of the mode with the highest energy content) that distinguishes structural from erroneous, or “artificial” modes [70]. The dispersion-corrected FSD indicates that mode and dispersion stabilization occurs at $p=38$, which is the finally selected order for all three ARX models. Figure 9 illustrates the behaviour of the ARX${}_3(38,38)$ prediction errors, from where the hypothesis of Gaussian white noise process can be safely adopted.

Based on these results, ARX${}_{1,t_l}(38,38)$ (unit delay), ARX${}_{3,t_l}(38,38)$ (zero delay), and ARX${}_{5,t_l}(38,38)$ (unit delay) models are subsequently estimated for $t_l=1,\dots ,300$ and their parameter vectors and prediction error variances are stored for the GPR modelling stage. As

Table 1. Spring-mass-damper: statistics of percentage fitness to the short time estimation data set during ARX modelling.

	ARX${}_{1,{\mathit{t}}_{\mathit{l}}}(38,38)$	ARX${}_{3,{\mathit{t}}_{\mathit{l}}}(38,38)$	ARX${}_{5,{\mathit{t}}_{\mathit{l}}}(38,38)$
${\mu }_{fit}$	98.91%	98.93%	98.91%
${\sigma }_{fit}$	$4.2\times {10}^{-2}$	$1.9\times {10}^{-2}$	$5.0\times {10}^{-2}$

indicates, a quite satisfying degree of consistency may characterize the parametric identification process, since the percentage fitness

$$fit=100\left(1-\frac{\parallel y\left[t\right]-\widehat{y}\left[t\right]\parallel }{\parallel y\left[t\right]-{\mu }_y\parallel }\right)(\%)$$

(38)

to the short time estimation data set is sufficiently high.

Following the structural identification stage, GPR models in the form of Equation (17) are fitted to a total number of 228 parameters. It is noted that no modelling is performed for the $b_0$ parameters: this is due to the fact that the AFS-ARX${}_1(38,38)$ and AFS-ARX${}_5(38,38)$ models are characterized by unit delay, while the $b_0$ parameter of the AFS-ARX${}_3(38,38,0)$ model is consistently estimated at 0.01 during the structural identification stage, with negligible fluctuations around this value.

The selected basis functions and GPR orders, determined via an initial trial procedure, are listed in

Table 2. Spring-mass-damper system: GPR basis functions, kernels and orders.

Quantity	AR Parameters	Exogenous Parameters	Noise Variance
Basis Function Type	Discrete Fourier		Polynomial
Basis Function Formula	$S_{a,j}(\Xi )=cos\left(\frac{2\pi j}{300}\Xi \right)$, $j=0,2,4,\dots$		$S_{b,j}(\Xi )={\Xi }_j$, $j=0,1,2,\dots$
	$S_{a,j}(\Xi )=sin\left(\frac{2\pi j}{300}\Xi \right)$, $j=1,3,5,\dots$
Kernel Type	Matérn 3/2
Kernel Formula	$\kappa ({\Xi }_i,{\Xi }_j)={\sigma }^2\left(1+\frac{\sqrt{3}({\Xi }_i-{\Xi }_j)}{\rho }\right)e^{-\frac{\sqrt{3}({\Xi }_i,{\Xi }_j)}{\rho }}$
Order	$p_a=3$	$p_b=17$	$p_{\sigma }=3$

, while the percentage fitness results are shown in Figure 10 and

Table 3. Spring-mass-damper: GPR modelling results for the noise variance of the associated ARX models.

	GP-ARX${}_1(38,38)$	GP-ARX${}_3(38,38)$	GP-ARX${}_5(38,38)$
$fit$	91.56%	83.75%	85.37%

. In general, sufficiently high fitting to the parameters has occurred, reaching, in many cases, more than 99%. There exist some AR and eXogenous parameters with low modelling quality, as for example the $a_3$ of GP-ARX${}_3(38,38)$ (Figure 10a, second row) and the $b_{34}$ and $b_{36}$ of the same model (Figure 10b, second row); these exhibit wider scattering as compared to others, therefore rendering the fitting process more demanding. Still, an effective GPR representation for these is succeeded. This can be visualized in Figure 11, where indicative behaviour of the fitted GPR models is shown. In specific, the associated plots correspond to AR, eXogenous and variance parameter models with the lowest and highest percentage fitness.

The behaviour of the established GP-ARX models is now assessed using available information that was not used during the estimation process. In specific, the RGA stage is implemented to two long time-scale instants, which are chosen from the set $t_l=301,\dots ,365$. Notice that this corresponds to a form of extrapolation, since training was not applied to the whole long time-scale cycle. Figure 12 displays the RGA time-series of the available GP-ARX models for $t_l=320$ and $t_l=355$. A quite stable and stationary behaviour is observed, indicating that no damage has occurred within this long time-scale. Yet, in some cases the associated thresholds might result underestimated, as, for example, the ones of Figure 12d. This issue is rather attributed to the fact that the inherited uncertainties associated with the GPR models (e.g., the GPR prediction intervals) are not taken into account. Despite this inconsistency, the GP-ARX models are still capable of tracking the current state of the structure, as indicated in Figure 13, where the predicted time-series are plotted over the actual ones for the RGA of Figure 12d (e.g., the one with the “worst” behaviour).

Finally, Figure 14 illustrates the behaviour of the RGA time-series under damage, during the long time-scale $t_l=320$. The latter is introduced as a local stiffness change, which occurs abruptly at one-third of the total simulation time. The first four rows of the figure display the $r_i\left[k\right]$′s for 1% damage in $K_1$, $K_3$, $K_5$ and $K_7$, respectively, from where the detection effectiveness of the algorithm even at such low levels is demonstrated. Notice that there’s no correlation of the RGA time-series amplitudes to the location of damage, since all remain around the same levels. This can be also argued for the damage severity: the RGA time-series of the last row, which corresponds to 50% damage in $K_4$, result with increased amplitudes (around double the ones of the first four rows), yet, again this change is quite disproportional to the size of the damage (50% over 1%).

4.2. Damage Detection and Localization on a Shear Frame

The second case study pertains to the six-storey shear frame of Figure 15 under base (e.g., earthquake) excitation. It is again assumed that the stiffness is exponentially depended on temperature via the same relationship as in the spring-mass-damper case and that the temperature varies as in Equation (36). This implies that Figure 3 and Figure 4 are also valid here, with the only difference being that $K_i$ is expressed now on a MN/m scale and that modal damping is kept constant at 2%. The natural frequencies of the frame are shown in Figure 16.

In deriving the theoretical form of the GP-ARX models to be estimated, the equation of motion of the frame obeys

$$\begin{array}{c}\hfill \mathbf{M}\ddot{\mathbf{z}}\left(t\right)+\mathbf{D}\dot{\mathbf{z}}\left(t\right)+\mathbf{K}\left(T\right)\mathbf{z}\left(t\right)=-\mathbf{M}\mathbf{1}{\ddot{x}}_g\left(t\right)\end{array}$$

(39)

where $\mathbf{z}\left(t\right)$ is the vector DOFs’ displacements relative to the ground motion, e.g., $z_i=x_i-x_g$ (see Figure 15 and $\mathbf{1}$ is a column vector filled with ones. The state-space representation of Equation (39) for absolute acceleration output is

$$\begin{array}{cc}\hfill \dot{\mathit{\zeta }}\left(t\right)& ={\mathbf{A}}_c\left(T\right)\mathit{\zeta }\left(t\right)+{\mathbf{B}}_c{\ddot{x}}_g\left(t\right)\hfill \end{array}$$

(40a)

$$\begin{array}{cc}\hfill \mathbf{y}\left(t\right)& ={\mathbf{H}}_c\left(T\right)\mathit{\zeta }\left(t\right)\hfill \end{array}$$

(40b)

indicating that the equivalent digital transfer function between the acceleration of the storeys and the ground acceleration is characterized by unit delay (refer to the analysis of Section 2).

Following the steps of Section 4.1, the “healthy” structure is discretized at $T_s=0.02$ s ($F_s=50$ Hz) via the zero-order hold principle. Accordingly, the structural absolute vibration acceleration responses under independent realizations of zero-mean, Gaussian white noise base excitation (${\ddot{x}}_g\left[k\right]\in \mathcal{N}(0,1.87\times {10}^4)$, data length $N=2^{14}$) are obtained and zero-mean subtracted, while it is assumed that measurements from all DOFs are available. Figure 17 shows the absolute vibration acceleration FRF of the first story at reference temperature.

Model order selection is based on the estimation of unit delay ARX($p,p$) models for $p=10,12,\dots ,50$, while using the excitation and absolute vibration acceleration response of the first storey at a randomly chosen long time instant. The performance of the model order selection criteria and the associated dispersion-corrected FSD is depicted in Figure 18, which qualify $p=12$ as the best order. This complies with the theoretical order of Equation (15), a result that is expected due to the absence of noise-corrupted data. The behaviour of the unit delay ARX(12,12) prediction errors returns quite similar to the one of Figure 9. The subsequent identification of unit delay ARX${}_{i,t_l}(12,12)$ models for $i=1,2,\dots ,6$ and $t_l=1,\dots ,300$ returns very high mean fitness with low dispersion, as shown in

Table 4. Shear frame: statistics of the percentage fitness to the short time estimation data set during ARX modelling.

	ARX${}_{1,{\mathit{t}}_{\mathit{l}}}(12,12)$	ARX${}_{2,{\mathit{t}}_{\mathit{l}}}(12,12)$	ARX${}_{3,{\mathit{t}}_{\mathit{l}}}(12,12)$	ARX${}_{4,{\mathit{t}}_{\mathit{l}}}(12,12)$	ARX${}_{5,{\mathit{t}}_{\mathit{l}}}(12,12)$	ARX${}_{6,{\mathit{t}}_{\mathit{l}}}(12,12)$
${\mu }_{fit}$	99.97%	99.98%	99.99%	99.98%	99.98%	98.98%
${\sigma }_{fit}$	$5.7\times {10}^{-3}$	$4.6\times {10}^{-3}$	$4.2\times {10}^{-3}$	$3.9\times {10}^{-3}$	$3.6\times {10}^{-3}$	$3.0\times {10}^{-3}$

.

The GPR modelling phase is carried out using the basis functions, kernels, and orders of

Table 5. Shear frame: GPR basis functions, kernels, and orders.

Quantity	AR Parameters	Exogenous Parameters	Noise Variance
Basis function type	polynomial		discrete Fourier
Basis function formula	as in Table 2		as in Table 2
Kernel type	Matérn 3/2
Kernel formula	as in Table 2
Order	$p_a=3$	$p_b=3$	$p_{\sigma }=16$

.

Table 6. Shear frame: GPR modelling results for the noise variance of the associated ARX models.

	ARX${}_{1,{\mathit{t}}_{\mathit{l}}}(12,12)$	ARX${}_{2,{\mathit{t}}_{\mathit{l}}}(12,12)$	ARX${}_{3,{\mathit{t}}_{\mathit{l}}}(12,12)$	ARX${}_{4,{\mathit{t}}_{\mathit{l}}}(12,12)$	ARX${}_{5,{\mathit{t}}_{\mathit{l}}}(12,12)$	ARX${}_{6,{\mathit{t}}_{\mathit{l}}}(12,12)$
$fit$	16.42%	16.49%	16.70%	16.79%	16.83%	16.74%

and Figure 19 show the percentage fitness, confirming the efficacy of the GPR models in fitting of the AR and eXogenous parameters. The poor results of the ARX model variances are again attributed to the noise-free data used during the identification case. The typical behaviours of the GPR models are illustrated in Figure 20.

Figure 21 displays the short time-scale temporal evolution of the RGA time-series for $t_l=305$ under no damage. All of the residuals are stationary and within the specified thresholds, which herein result a bit overestimated, due to the low modelling quality of the respective GPR models for the ARX variances. In contrast, for the same short time-scale instant, Figure 22 illustrates the same quantities under damage, introduced as a 1% stiffness change in $K_1$, occurring abruptly at one-third of the total simulation time. One can here observe that there might be an indication about the location of the damage, by observing the descending dispersion from $r_1\left[k\right]$ to $r_6\left[k\right]$.

However, a better interpretation is offered by examining the MSCs differences, as described in Section 3.2 and Equation (35). These are plotted in Figure 23, from where it follows that the damage is indeed successfully localized in the differences of the first three modes, whereas the three trailing ones appear more “noisy” failing to localize the damage in $K_1$. To investigate this result, Figure 24 presents the “theoretical” MSCs differences, calculated via the eigenvalue problem of the exact structural matrices in healthy and damaged states. The close resemblance of the estimated MSCs differences to their theoretical counterparts is apparent, implying that (i) the reproduction of the modal space via the distributed GP-ARX models is to a great extend accurate; and, (ii) a damage may be localized only in lower modes.

5. Conclusions

This study tackles the damage detection and localization problem in structural systems that are amenable to environmental variations, and for which no prior information is available. The solution involves the estimation of distributed GP-ARX models for the representation of individual single input-single output pairs and the subsequent establishment of corresponding RGAs. These operate in parallel to the monitored system, are characterized by well-defined stationary statistics in the healthy state over the whole operational regime, and deviate from these in the presence of damage. In the latter case, localization is achieved by truncating all local GP-ARX models into a global state-space realization, deriving the modal space and calculating the MSC discrepancies between the healthy and damaged state-space representations.

The assessment of the method through the case studies indicates that the proposed framework is capable of detecting and localizing damage at quite low levels (“early” enough). The detection succeeds, even for a single input-output pair, while localization presumes a proper sensor network. It would be interesting to relax this requirement and attempt to integrate the damage localization task into the RGA, by configuring the residuals in some directional, or structured form; this forms current work of the authors. Nevertheless, it seems that the establishment of RGAs with localization attributes requires qualitative knowledge of structural dynamics, so that the GP-ARX parameter vector can be linked to the structural parameters.

Another issue that is observed from the case studies pertains to numerical stability. Indicatively, in the shear frame paradigm, it is found that some operating points from the validation set result in unstable GP-ARX models, even if the GPR modelling stage produces very accurate parametric representations. An alternative approach might lie in investigation of vector GP-ARX models, in order to bypass the local-to-global state-space realization, which could induce errors. Because vector models retain a considerably richer structure, an open research question would be the potential of performing GPR modelling on a reduced parameter space, similarly to what is discussed in [74].

It must be noted that structural damage may appear either as gradual degradation at slow rates (e.g., crack initiation/propagation), or as faults after “events” (e.g., an earthquake, strong wind/wave gusts, etc.). Our work is mostly set to treat this latter class: we assume distinct operational points along the long time-scale, we estimate the current state of the structure and subsequently forming/updating the RGA time-series. However, this does not imply that our method is incapable of treating the first class as well, since it has been demonstrated that the RGA can detect damage early enough, while it can easily operate in continuous mode, as long as measurements are available. However, the “footprint” of any gradual degradation on the RGA requires additional research, as mentioned above. Currently, our RGA does not correlate with damage evolution: it is just a yes/no answer to detection, which considers the long time-scale structural variability.

The encouraging results on the efficacy of the method suggest further research steps towards a more realistic SHM context. These include the addition of further free parameters into the GPR metamodeling stage, related to further environmental and operational parameters (operating points, boundary conditions, etc.), aside from temperature. The relaxation on the availability and the stationarity assumption of the input(s) forms a further critical aspect that is frequently encountered in practice: here, the investigation of scalar/vector GP-AR/ARMA models and their ARIMA counterparts form a straightforward task, which is not yet adequately treated in existing literature. Finally, an experimental validation is currently underway by the authors.