Domain-adaptive nonlinear probabilistic latent variable model for vibration-based structural damage detection under climate change

Abstract

A major challenge in vibration-based structural damage detection is differentiating actual damage from environmental and operational variations. Most existing studies are limited to short-term scenarios and do not consider long-term variations. With climate change driving more frequent and intense extreme weather events, long-term structural monitoring data exhibit covariate shift, for example, annual mean temperature increases over time. Such shifts pose fundamental challenges for machine learning models trained on historical datasets, as traditional algorithms often fail to extrapolate to these out-of-distribution environmental domains. This study develops a domain-adaptive nonlinear probabilistic latent variable model to advance long-term damage detection under climate change. First, Bayesian inference and kernel techniques are integrated to capture underlying nonlinear environmental effects without explicit environmental data. Second, domain adaptation is incorporated to enhance model adaptability to covariate shift scenarios due to climate change. The developed method is applied to a laboratory-tested slab and the practical Z24 bridge. Structural vibration properties are estimated under changing temperature and humidity conditions with consideration of four climate emission scenarios. The results demonstrate that the proposed method can distinguish damage-induced variations from climate-related environmental effects. It outperforms traditional methods that may result in false-positive and false-negative damage detections due to climate change and environmental nonlinear impacts.

Introduction

Structural health monitoring (SHM) is an advanced technology to continuously monitor and evaluate safety conditions of civil infrastructures. Structural damage is defined as changes in the physical parameters that adversely affect structural current or future performance. ¹ Vibration-based structural damage detection techniques have been extensively developed over the past decades,^2–4 which can be broadly classified into model-based and data-driven categories. While model-based methods rely on finite element updating, ⁵ which often suffers from substantial modeling uncertainties in boundary conditions, data-driven methods ⁶ leverage statistical pattern recognition, making them highly effective and adaptable to various types of structures for real-time monitoring without the need for physical models.

The impact of environmental and operational variations on structural vibration properties poses a great challenge to data-driven damage detection methods.^7–9 Civil structures are typically subjected to a wide range of operational and environmental conditions, such as temperature fluctuations, humidity changes, and varying loads, which may also lead to changes in structural damage-sensitive features. Failing to account for these effects can obscure actual structural damage and lead to false-positive (i.e., healthy structures are mistakenly identified as damaged) or false-negative (i.e., actual damage is misidentified as a healthy condition) detections. ¹⁰ Therefore, filtering out these environmental influences prior to damage detection is critical. The advancement of machine learning (ML) techniques, known for their strong capabilities in big data mining, has led to widespread exploration of ML-based data-driven damage detection methods. Farrar and Worden ¹¹ provided a comprehensive overview of ML techniques for SHM and highlighted their potential to enhance the accuracy and reliability of damage detection. Bao and Li ¹² summarized the ML paradigm for SHM challenges including data processing and pattern recognition, environmental uncertainties, improving damage detection accuracy and reliability.

Dimensionality reduction is a key branch of ML methods to distinguish the effects of environmental variations from actual structural damage. ¹¹ In many cases, the variations in observations are influenced by only a few underlying factors. Dimensionality reduction techniques map high-dimensional data into a lower-dimensional space, capturing essential environmental factors while filtering noise. Healthy data can be reconstructed using these reduced dimensions, and subsequent new measurements that exhibit high reconstruction errors are then flagged as damage. Representative dimensionality reduction techniques include cointegration, ¹³ independent component analysis,^14,15 linear discriminant analysis, ¹⁶ factor analysis,^7,17 non-negative matrix factorization, ¹⁸ and principal component analysis (PCA).^8,19 However, these models face two major constraints. First, they are limited to linear models, while environmental effects on structural vibration properties can be nonlinear. Second, ML models may experience performance degradation when dealing with out-of-distribution detection problems, where the model encounters data samples that follow a different distribution from the training data (i.e., distribution shift).

The first issue regarding the nonlinear environmental effects on structural vibration properties may arise due to varying temperature-sensitive mechanical properties of materials and complicated boundary conditions. ⁹ For example, Peeters and De Roeck ⁹ monitored the Z24 bridge for around 10 months and found that natural frequencies had a nonlinear relation with the temperature due to the frozen asphalt layer in cold days. Oh et al., ²⁰ Reynders et al., ²¹ and Ghoulem et al. ²² applied kernel PCA to handle nonlinear damage detection examples. However, kernel PCA necessitates the manual selection of optimal parameters for the kernel function, and the dimensionality to be retained in the feature space is typically determined based on researchers’ experiences or through a trial-and-error approach. Additionally, the non-probabilistic model does not account for uncertainties that inherently exist in SHM data. More recently, advanced nonlinear frameworks have been exploited and demonstrated powerful feature extraction capabilities, such as deep learning-based approaches. ²³ The black-box deep learning models typically require massive training data and lack interpretability. Alternatively, the Gaussian process latent variable model ²⁴ offers a rigorous probabilistic approach to nonlinear dimensionality reduction, by providing a flexible mapping from the latent space to observation space and handling uncertainty. Nevertheless, the application of these models is currently limited to short-term SHM scenarios.

The second issue is that damage is a slow, progressive, and long-term process. The environmental factors may keep changing gradually in the context of climate change. Over a structure’s lifespan, the operational data will inevitably fall outside the distribution of the historical data used to train SHM models, thus demanding more robust and adaptive methods. Climate change leads to heightened intensity of monsoons, ²⁵ increased frequency and strength of tropical cyclones, ²⁶ greater duration of heat waves, ²⁷ and elevated concentration of CO₂. ²⁸ The global average temperature has risen by more than 0.7°C from 1986 to 2016 compared to the period from 1901 to 1960. ²⁹ The extreme temperature values have increased at a much faster rate than the mean values, and some regions experienced a more increase than others. ³⁰ The impact of climate change on infrastructure has been investigated by many researchers, including the corrosion and deterioration of metals ³¹ and timber materials, ³² deterioration and creep of reinforced concrete,^33–35 and wind-induced damages to infrastructures. ³⁶ However, there is limited research related to the impact of climate change on vibration-based structural damage detection, which has recently become a critical area of investigation as the non-stationarity introduced by a changing climate can render historical training data obsolete and compromise the reliability of ML-based SHM systems. For example, Figueiredo et al. ³⁷ have demonstrated through studies on the Z24 bridge data that significant shifts in mean temperature can cause trained classifiers to produce false negatives, misidentifying damaged conditions as healthy due to outdated reference datasets. To address this challenge and develop more resilient SHM strategies, Möller et al. ³⁸ have proposed enhancing data-driven models with physical knowledge to improve robustness against unobserved environmental conditions. Furthermore, Quqa et al. ³⁹ have put forward the idea of establishing new regional-scale monitoring paradigms that integrate remote sensing and knowledge transfer to adapt to emerging climate risks. Performance degradation is a common and challenging issue when ML models encounter data that significantly differs from the training data. ⁴⁰ Many existing ML-based data-driven damage detection methods are trained using historical data. Their performance will inevitably deteriorate applied to long-term damage detection, since the test data will fall outside the range of the training data in the context of climate change. To mitigate the issue of distribution shift (or domain shift), domain adaptation techniques have been widely explored to align the statistical distributions of features across different domains to improve generalization.^40–43 However, most existing studies in SHM focus on population-based transfer (between different structures) ⁴³ or simulation-to-real ⁴² transfer. There is still a lack of specialized domain adaptation frameworks specifically designed to handle the temporal distribution shift caused by long-term non-stationary climate change within a single structure’s lifespan. Therefore, improving the robustness of ML models to long-term distribution shifts under the impact of climate change is highly necessary.

Considering the research gaps stated above, this study develops a novel dimension reduction approach for long-term damage detection in the context of climate change. The primary contributions of this study are as follows:

1. From a methodological perspective, this study establishes a rigorous nonlinear latent variable model that integrates domain adaptation directly into a probabilistic framework. Rather than simply combining techniques, the feature augmentation is embedded within the kernel function of a sparse Gaussian process model. This novel embedding allows the model to simultaneously address structural nonlinearity and retain the Bayesian probabilistic framework for efficient automatic parameter estimation without requiring explicit correlation expressions. Crucially, this formulation provides the mathematical capacity to compensate for distribution shifts, a capability lacking in standard dimensionality reduction techniques.

2. From a problem-solving perspective, this study overcomes the fundamental limitation of stationary environmental assumption in current SHM literature by applying transfer learning to the context of climate change. The model’s performance in out-of-distribution extrapolation scenarios where test data falls outside the range of the training data is enhanced. The performance is tested for climate change projections under four shared socioeconomic pathways (SSPs) scenarios 126, 245, 370, and 585, demonstrating the ability to mitigate false-positive damage detection caused by data extrapolation and false-negative detection caused by latent nonlinearities.

Domain-adaptive nonlinear probabilistic latent variable model

Denote the measured structural vibration properties (such as frequencies) as $\mathbf{D}=[{\mathbf{D}}_1,{\mathbf{D}}_{2,}...,{\mathbf{D}}_N]\in {\mathbb{R}}^{N_m\times N}$ , where $N_m$ is the dimension of each observation and $N$ is the total number of observation samples, and the environmental factors (such as temperature, humidity) that affect structural vibration properties as $\mathbf{z}$ . For a collection of observations, there are a set of variables that can describe the variability of the observations. When only the observations are available, the correlation model can be constructed by treating these variables latent, that is, $\mathbf{z}$ is unknown. The latent variable model is expressed as

$$\mathbf{\Phi }\left({\mathbf{D}}_n\right)=\mathit{f}({\mathbf{z}}_n)+\mathbf{\epsilon }$$

(1)

where ${\mathbf{z}}_n$ is an M-dimensional (M <  $N_m$ ) latent variable with prior $p\left(\mathbf{z}\right)=\mathcal{N}(\mathbf{z}|0,\mathbf{I})$ ; $\mathbf{\Phi }$ and $\mathit{f}$ denote implicit or explicit functions to reveal the correlations between $\mathbf{D}$ and $\mathbf{z}$ ; and $\mathbf{\epsilon }$ is a $N_m$ -dimensional error vector that accounts for the uncertainties and model prediction errors. $\mathbf{\epsilon }$ is generally assumed to follow a Gaussian distribution with the zero-mean and covariance matrix ${\beta }^{-1}\mathbf{I}$ , that is, $p\left(\mathbf{\epsilon }\right)=\mathcal{N}(\mathbf{\epsilon }|0,{\beta }^{-1}\mathbf{I})$ . In this model, only ${\mathbf{D}}_n$ is available, and the rest parameters $\mathbf{z}$ and $\mathbf{\epsilon }$ as well as functions $\mathbf{\Phi }$ and $\mathit{f}$ are unknown.

The formulation presented in Equation (1) has been extensively employed for structural damage detection subject to varying environmental conditions.^7,21,22 Unlike black-box models, this probabilistic latent variable formulation inherently maintains physical interpretability by modeling the underlying relationship between structural responses and unobserved environmental factors. This formulation encompasses several classical methods as special cases:

1. Linear factor analysis.^7,44 When both mapping functions $\mathbf{\Phi }$ and $\mathit{f}$ are linear, the equation reduces to ${\mathbf{D}}_n=\mathbf{W}{\mathbf{z}}_n+\mathbf{\epsilon }$ (W is a matrix for dimension reduction), which is mathematically equivalent to linear factor analysis model.^7,21 Using the Bayesian inference, the linear factor analysis can optimize the unknown variables and parameters automatically. However, the model is limited to be linear and thus cannot adequately explain the variations in observations in many practical nonlinear cases.

2. Kernel PCA.^21,22 When $\mathbf{\Phi }$ represents a nonlinear kernel mapping and $\mathit{f}$ remains linear, expressed as $\mathbf{\Phi }\left({\mathbf{D}}_n\right)={\mathbf{z}}_n+\mathbf{\epsilon }$ , the model aligns with kernel PCA. The kernel trick allows for the possibility of mapping data to a very high or even infinite-dimensional space, making the kernel PCA suitable for estimating nonlinear correlations. However, Kernel PCA requires an N × N kernel matrix calculation and eigenvalue decomposition, incurring a prohibitive O(N ³ ) computational cost that makes it impractical for long-term SHM datasets ( $N$ is significantly large).

Therefore, an ideal model should maintain probabilistic properties for automatic parameter optimization and the nonlinear flexibility of kernel PCA, but with significantly lower computational demands. Hence, this study interprets the model in Equation (1) as both probabilistic and nonlinear, to address the nonlinear effects and distribution shifts inherent in long-term damage detection. A latent variable model is defined as

$${\mathbf{D}}_n=\mathit{f}({\mathbf{z}}_n)+{\mathbf{\epsilon }}_n$$

(2)

where ${\mathbf{D}}_n$ , $\mathit{f}$ , ${\mathbf{z}}_n,$ and ${\mathbf{\epsilon }}_n$ are same as the definition in Equation (1). Since determining an explicit expression for the nonlinear function $\mathit{f}$ is often intractable, the kernel trick will be adopted that computes the inner product of two vectors in a high-dimensional feature space without explicitly computing the coordinates of the data in that space. The specific techniques are introduced as follows.

Maximum likelihood estimation of the latent variable

The estimation of unknown parameters and variables in Equation (2) begins with a tractable linear model and is then extended to the nonlinear case via kernel substitution. Consider the linear model ${\mathbf{D}}_n=\mathbf{W}{\mathbf{z}}_n+\mathbf{\epsilon }$ , where W is an $N_m\times M$ matrix whose columns span data subspaces. Unlike conventional factor analysis method, which marginalizes the latent variable $\mathbf{z}$ to optimize the parameters $\mathbf{W}$ , this study adopts a dual approach, namely, the parameters $\mathbf{W}$ is marginalized out to optimize the latent variable $\mathbf{z}$ . For simplification, a Gaussian prior conjugate to the likelihood function is adopted as $p\left(\mathbf{W}\right)=\underset{i}{\overset{D}{\Pi }}\mathcal{N}({\mathit{w}}_i|0,\mathbf{I})$ , where ${\mathit{w}}_i$ is the ith row of $\mathbf{W}$ . Then the marginalization of the likelihood with respect to W is

$$\begin{array}{c}\multicolumn{1}{c}{p({\mathbf{D}}_{d,:}|\mathbf{z},\beta )=\int p({\mathbf{D}}_{d,:}|{\mathit{w}}_i,\beta ,\mathbf{z})p\left({\mathit{w}}_i\right)d{\mathit{w}}_i}\\ \multicolumn{1}{c}{=\mathcal{N}({\mathbf{D}}_{d,:}|0,{\mathbf{z}}^{\mathrm{T}}\mathbf{z}+{\beta }^{-1}\mathbf{I})}\end{array}$$

(3)

where ${\mathbf{D}}_{d,:}$ represents the dth row of D $\in {\mathbb{R}}^{N_m\times N}$ . In corresponding, the objective function of the likelihood function in the logarithm form is obtained as

$$\begin{array}{c}\multicolumn{1}{c}{L=\ln p(\mathbf{D}|\mathbf{z},\beta )=\sum _{d=1}^{D}p({\mathbf{D}}_{d,:}|\mathbf{z},\beta )=-\frac{N_{m}N}{2}\ln 2\pi }\\ \multicolumn{1}{c}{-\frac{N_m}{2}\ln \left|{\mathbf{K}}_{\mathbf{z}}\right|-\frac{1}{2}\mathrm{tr}\left({\mathbf{K}}_{\mathbf{z}}^{-1}{\mathbf{D}}^{\mathrm{T}}\mathbf{D}\right)}\end{array}$$

(4)

where ${\mathbf{K}}_{\mathbf{z}}={\mathbf{z}}^{\mathrm{T}}\mathbf{z}+{\beta }^{-1}\mathbf{I}$ is the covariance matrix or kernel function, which is based on the linear inner product of the latent variable $\mathbf{z}$ .

By replacing the inner product ${\mathbf{K}}_{\mathbf{z}}$ with nonlinear kernel functions $\mathit{f}{\left({\mathbf{z}}_m\right)}^{\mathrm{T}}\mathit{f}\left({\mathbf{z}}_n\right)$ in the reproducing kernel Hilbert space, the nonlinear projection is achieved. In this way, the explicit mapping function $\mathit{f}\left({\mathbf{z}}_n\right)$ is not required to be known, thereby avoiding searching for the high-dimensional feature space $\mathcal{F}$ . The representative kernel functions include the polynomial, sigmoid, and radial basis function (RBF) kernel. The polynomial kernel is defined as

$${\mathit{k}}_{\mathbf{z}}({\mathbf{z}}_m,{\mathbf{z}}_n)={({\mathbf{z}}_m·{\mathbf{z}}_n+1)}^d$$

(5)

where ${\mathit{k}}_{\mathbf{z}}({\mathbf{z}}_m,{\mathbf{z}}_n)$ is the element in the mth row and the nth column of the kernel matrix ${\mathbf{K}}_{\mathbf{z}}$ , and $d$ is a positive integer. The Gaussian kernel function is defined as ²¹

$${\mathit{k}}_{\mathbf{z}}({\mathbf{z}}_m,{\mathbf{z}}_n)=\exp (-\frac{1}{2{\sigma }^2}{({\mathbf{z}}_m-{\mathbf{z}}_n)}^{\mathrm{T}}({\mathbf{z}}_m-{\mathbf{z}}_n))$$

(6)

where ${\sigma }^2$ is the variance parameter that controls the bandwidth of ${\mathbf{K}}_{\mathbf{z}}$ , which attains an infinite-dimensional feature space using merely a unique parameter. The RBF kernel, a more general function, is defined as

$$k({\mathbf{z}}_i,{\mathbf{z}}_j)={\theta }_{\mathrm{rbf}}\exp [-\frac{\gamma }{2}{({\mathbf{z}}_i-{\mathbf{z}}_j)}^{\mathrm{T}}({\mathbf{z}}_i-{\mathbf{z}}_j)]+{\theta }_{\mathrm{bias}}+{\theta }_{\mathrm{white}}{\delta }_{\mathrm{ij}}$$

(7)

where ${\theta }_{\mathrm{rbf}}$ is the process variance, which influences the scale of the output functions; $\gamma$ is the inverse width parameter that controls the variance; ${\theta }_{\mathrm{bias}}$ is corresponding to the prior variance; ${\theta }_{\mathrm{white}}$ is the white noise term; and ${\delta }_{\mathrm{ij}}$ is the Kronecker delta function. The RBF kernel was chosen for this study due to its universality and flexibility. As a universal approximator, the RBF kernel can model arbitrarily complex nonlinear relationships without prior assumptions about the data’s functional form. The inherent smoothness is well-suited for general-purpose nonlinear modeling.

With the kernel function selected, the next step is to optimize the latent variables ${\mathbf{z}}_n$ and parameters of the kernel function, including $\gamma$ , ${\theta }_{\mathrm{rbf}}$ , ${\theta }_{\mathrm{bias}}$ , and ${\theta }_{\mathrm{white}}$ . The likelihood of the parameters and latent variables $\mathbf{z}$ can be estimated by setting the gradient of the log-likelihood in Equation (4) to zero. In particular, the gradient with respect to $\mathbf{z}$ is computed through the chain rule as

$$\frac{\partial L}{\partial \mathbf{z}}=\frac{\partial L}{\partial {\mathbf{K}}_{\mathbf{z}}}\times \frac{\partial {\mathbf{K}}_{\mathbf{z}}}{\partial \mathbf{z}}=0$$

(8)

where $\frac{\partial L}{\partial {\mathbf{K}}_{\mathbf{z}}}={\mathbf{K}}_{\mathbf{z}}^{-1}\mathbf{D}{\mathbf{D}}^{\mathrm{T}}{\mathbf{K}}_{\mathbf{z}}^{-1}-N_m{\mathbf{K}}_{\mathbf{z}}^{-1}$ ; and $\frac{\partial {\mathbf{K}}_{\mathbf{z}}}{\partial \mathbf{z}}$ can be directly obtained according to Equation (7). However, obtaining closed-form solutions from Equation (8) is difficult as the log-likelihood L has highly nonlinear relationships with the latent variable and kernel parameters. To address this problem, the scaled conjugate gradient approach ²⁴ is utilized for the gradient-based optimization. The gradients required are computed using the chain rule, that is, the gradient of the likelihood with respect to the kernel matrix $\frac{\partial L}{\partial {\mathbf{K}}_{\mathbf{z}}}$ is derived, which is then combined with the partial derivatives of the kernel function with respect to its hyperparameters $\frac{\partial L}{\partial \gamma }$ , $\frac{\partial L}{\partial {\theta }_{\mathrm{rbf}}}$ , $\frac{\partial L}{\partial {\theta }_{\mathrm{bias}}}$ , $\frac{\partial L}{\partial {\theta }_{\mathrm{wh}\mathrm{ite}}}$ or variables $\frac{\partial {\mathbf{K}}_{\mathbf{z}}}{\partial \mathbf{z}}$ . Another challenge is the computational complexity, which renders the method inefficient when handling a large number of data samples. To improve computational efficiency, Lawrence et al.^24,45 developed the informative vector machine (IVM) as a sparsification mechanism. The active set, denoted as ${\mathbf{D}}_I$ , is sequentially selected from the overall training dataset and then used to optimize the kernel parameters and latent variables. The IVM reduces the dominant computational cost to $\mathrm{O}(d^2·N)$ , where d is the amount of data in the active set (d < N), which is much more efficient than the kernel PCA model $\mathbf{\Phi }\left({\mathbf{D}}_n\right)={\mathbf{z}}_n+\mathbf{\epsilon }$ that needs $\mathrm{O}\left(N^3\right)$ computational cost. In long-term damage detection problems, the number of data samples N can grow to be exceptionally large, rendering methods with $\mathrm{O}\left(N^3\right)$ complexity computationally impractical. The sparsification of IVM is critical for the model’s practical feasibility in the target application of long-term SHM.

Domain-adaptative sparse Gaussian process for observation reconstruction

With the latent variable $\mathbf{z}$ and model parameters estimated above, the measurement data $\mathbf{D}$ can be reconstructed. Since the explicit form of the nonlinear function $\mathit{f}$ is unknown, reconstructing the observations according to Equation (2) using the estimated latent variable is difficult. Gaussian process regression offers a non-parametric modeling tool capable of making predictions without constraining the input–output relationship to a specific form. However, standard Gaussian process regression demonstrates limited ability in out-of-distribution extrapolation. This poses a significant challenge for long-term SHM, where data is frequently subject to distribution variations caused by climate change. Such factors result in covariate shift where test data falls outside the range of the training data, thereby hindering the performance of the Gaussian process model in long-term damage detection.

Domain adaptation addresses the divergence in data distribution between source and target domains, facilitating the transfer of knowledge learned from a source domain to a related but different target domain.^43,46 To address the distribution divergence between the source domain (historical health data) and target domain (test data under climate change), the feature augmentation ⁴⁷ -based domain adaptation technique can be integrated with the Gaussian process model. The feature augmentation ⁴⁷ technique explicitly separates domain-invariant and domain-specific features by expanding the latent variable vector $\mathbf{z}$ . Suppose the source domain data ${\mathbf{D}}_s$ and target domain data ${\mathbf{D}}_t$ are described as

$${\mathbf{D}}_s=\{({\mathbf{z}}_1^s,{\mathbf{D}}_1^s),\left({\mathbf{z}}_2^s,{\mathbf{D}}_2^s\right),\dots ,({\mathbf{z}}_{N_s}^s,{\mathbf{D}}_{N_s}^s\left)\right\},$$

(9)

$${\mathbf{D}}_t=\{({\mathbf{z}}_1^t,{\mathbf{D}}_1^t),({\mathbf{z}}_2^t,{\mathbf{D}}_2^t),\dots ,({\mathbf{z}}_{N_t}^t,{\mathbf{D}}_{N_t}^t)\},$$

(10)

where $\mathbf{z}$ is the input variable vector (estimated latent variables in this study), and $\mathbf{D}$ is the output vector, and $N_s$ and $N_t$ is the amount of data in the source and target domain, respectively. The domain adaptation method implements data expansion as follows ⁴⁸ :

$${\hat{\mathbf{z}}}^{\mathrm{s}}=({\mathbf{z}}^s,{\mathbf{z}}^s,0)\mathrm{and}{\hat{\mathbf{z}}}^t=({\mathbf{z}}^t,0,{\mathbf{z}}^t),$$

(11)

where ${\hat{\mathbf{z}}}^s$ and ${\hat{\mathbf{z}}}^t$ is the expanded input vectors for the source and target domains, respectively. The expanded input and output data for Gaussian process regression are:

$$\begin{array}{l}{\mathbf{D}}_{\mathrm{all}}=\{({\widehat{\mathbf{z}}}_1^s,{\mathbf{D}}_1^s),({\widehat{\mathbf{z}}}_2^s,{\mathbf{D}}_2^s),\dots ,({\widehat{\mathbf{z}}}_{N_s}^s,{\mathbf{D}}_{N_s}^s),({\widehat{\mathbf{z}}}_1^t,{\mathbf{D}}_1^t),\\ ({\widehat{\mathbf{z}}}_2^t,{\mathbf{D}}_2^t),\dots ,({\widehat{\mathbf{z}}}_{N_t}^t,{\mathbf{D}}_{N_t}^t)\}.\end{array}$$

(12)

Such a feature augmentation process explicitly separates domain-invariant and domain-specific features. The ML model thus learns domain-invariant patterns and adapts to the target domain through the target-specific features. ⁴⁷ This approach was specifically chosen over other common domain adaptation techniques, such as those based on geometric alignment.^40–43 Unlike geometric alignment methods that try to force distributions to overlap, which can distort the physical meaning of environmental variables in regression tasks, feature augmentation preserves the physical structure of the data while learning a mapping that accommodates the shift. This preservation is particularly crucial for our goal of observation reconstruction, which is a regression task rather than simple classification.

A critical challenge in this study is that the target domain ${\mathbf{D}}_t$ (test data under climate change to be classified) is unlabeled, whereas standard feature augmentation requires labeled data in both domains to learn the domain-specific offsets. Under covariate shift, although the marginal distribution of inputs changes (i.e., the latent variables of the target domain may fall outside the range of the source domain), the underlying physical input–output mechanism remains invariant across environments. In other words, the shift mainly manifests as $p(\mathbf{z})$ changing with environmental conditions, while the conditional relationship between latent variables and measurements is governed by consistent structural behavior. Based on this assumption, we improve the supervised feature augmentation strategy to a semi-supervised calibration scheme. Specifically, we assume access to a massive source dataset ${\mathbf{D}}_s$ and a small set of labeled data from the out-of-distribution environment, referred to as the calibration set ${\mathbf{D}}_c$ . The key role of ${\mathbf{D}}_c$ is to provide minimal but essential supervision in the shifted environment, enabling the Gaussian process to preserve the domain-invariant mapping learned from $D_s$ , and simultaneously learn how this mapping should be adjusted in the new environment via the domain-specific subspace introduced by feature augmentation. Accordingly, the Gaussian process is trained on the union of the source and calibration data:

$$\begin{array}{c}\multicolumn{1}{c}{{\mathbf{D}}_{\mathrm{train}}=\{({\hat{\mathbf{z}}}_1^s,{\mathbf{D}}_1^s),({\hat{\mathbf{z}}}_2^s,{\mathbf{D}}_2^s),\dots ,({\hat{\mathbf{z}}}_{N_s}^s,{\mathbf{D}}_{N_s}^s),({\hat{\mathbf{z}}}_1^c,{\mathbf{D}}_1^c),}\\ \multicolumn{1}{c}{({\hat{\mathbf{z}}}_2^c,{\mathbf{D}}_2^c),\dots ,({\hat{\mathbf{z}}}_{N_c}^c,{\mathbf{D}}_{N_c}^c)\},}\end{array}$$

(13)

and the expanded latent variables of the training data are

$${\hat{\mathbf{z}}}^{\mathrm{s}}=({\mathbf{z}}^s,{\mathbf{z}}^s,0)\mathrm{and}{\hat{\mathbf{z}}}^c=({\mathbf{z}}^c,0,{\mathbf{z}}^c).$$

(14)

This construction explicitly decomposes the representation into a shared (domain-invariant) component and a domain-specific component. The abundant source data ${\mathbf{D}}_s$ stabilize the learning of the shared structure, while the labeled calibration data ${\mathbf{D}}_c$ anchors the model in the shifted region of the input space and identifies the environment-specific deviation. The Gaussian process trained on $({\mathbf{D}}_s\cup {\mathbf{D}}_c)$ can leverage the learned domain-invariant patterns together with the calibration-informed domain-specific adjustment to produce predictions for the unlabeled target inputs. Finally, during the testing phase, the unlabeled target data ${\mathbf{D}}_{\mathrm{t}}$ is mapped using the same transformation as the calibration set to generate predictions:

$${\mathbf{D}}_{\mathrm{t}}=({\mathbf{z}}_1^t,{\mathbf{z}}_2^t,\dots ,{\mathbf{z}}_{N_t}^t)\mathrm{where}{\widehat{\mathbf{z}}}^t=\left({\mathbf{z}}^t,0,{\mathbf{z}}^t\right).$$

(15)

With the expanded inputs, the Gaussian process model expresses the observation ${\mathbf{D}}_n$ as

$${\mathbf{D}}_n=\mathit{\psi }({\mathbf{z}}_n)+{\mathbf{\xi }}_n$$

(16)

where the distributions of $\mathit{\psi }(\mathbf{z})$ can be specified by the Gaussian distribution. Given a new latent variable ${\mathbf{z}}_{\mathrm{new}}$ , the corresponding prediction in the measurement space is denoted as ${\mathbf{D}}_{\mathrm{new}}^{*}$ . With the optimized parameters of the kernel function and the latent variables of the active data set, the latent variable ${\mathbf{z}}_{\mathrm{new}}$ and ${\mathbf{D}}_{\mathrm{new}}^{*}$ , including the inactive subset of the training data and the test data, are obtained from the scaled conjugate gradient ²⁴ optimization. The combination of observations $\mathbf{D}$ and the prediction ${\mathbf{D}}_{\mathrm{new}}^{*}$ follows a joint Gaussian distribution expressed as

$$\left\{\begin{array}{c}\mathbf{D}\hfill \\ \multicolumn{1}{c}{{\mathbf{D}}_{\mathrm{new}}^{*}}\end{array}\right\}~\mathcal{N}[\left(\begin{array}{c}\mathit{m}\left(\mathbf{z}\right)\hfill \\ \multicolumn{1}{c}{\mathit{m}\left({\mathbf{z}}_{\mathrm{new}}\right)}\end{array}\right),\left(\begin{array}{cc}\mathbf{K}\hfill & {\mathbf{K}}_{*}\hfill \\ \multicolumn{1}{c}{{\mathbf{K}}_{*}^{\mathbf{T}}}& {\mathbf{K}}_{**}\hfill \end{array}\right)]$$

(17)

where $\mathbf{K}=\left[\begin{array}{ccc}k({\mathbf{z}}_1,{\mathbf{z}}_1)\hfill & \cdots \hfill & k({\mathbf{z}}_1,{\mathbf{z}}_N)\hfill \\ \multicolumn{1}{c}{⋮}& ··\hfill & ⋮\hfill \\ \multicolumn{1}{c}{k({\mathbf{z}}_N,{\mathbf{z}}_1)}& \cdots \hfill & k({\mathbf{z}}_N,{\mathbf{z}}_N)\hfill \end{array}\right]$ , ${\mathbf{K}}_{*}=\left[\begin{array}{c}k({\mathbf{z}}_1,{\mathbf{z}}_{\mathrm{new}})\hfill \\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{k({\mathbf{z}}_N,{\mathbf{z}}_{\mathrm{new}})}\end{array}\right]$ , and ${\mathbf{K}}_{**}=k({\mathbf{z}}_{\mathrm{new}},{\mathbf{z}}_{\mathrm{new}})$ with the kernel function defined in Equation (7). In practice, a long-term damage detection problem involves a vast amount of SHM data, leading to high computational cost of Gaussian process. To improve computational efficiency, the sparse Gaussian process ⁴⁸ is adopted which contains two general strategies, namely, the subset of data points or regressors. The former activates a subset of the dataset for model parameter estimation. The latter replaces the kernel function with an approximated one based on the active set, which is expressed as follows:

$$\hat{k}(\hat{\mathbf{z}},{\hat{\mathbf{z}}}_{\mathrm{new}})=\sum _{\mathrm{\rho }=1}^{\omega }{\alpha }_{\mathrm{\rho }}k(\hat{\mathbf{z}},{\hat{\mathbf{z}}}_{\mathrm{new}})$$

(18)

where ${\alpha }_{\mathrm{\rho }}$ is the coefficient for the linear combination, which approximates the element in kernel function $k({\mathbf{z}}_i,{\mathbf{z}}_j)$ of the full Gaussian process by a linear combination of multiple kernel functions based on active sets.

Finally, the reconstruction of the target data is performed by calculating the posterior likelihood. The sparse Gaussian process model is trained based on the expanded training data in Equation (13), and then directly applied to the latent variables ${\hat{\mathbf{z}}}^t$ of the target domain for data reconstruction. The likelihood of ${\mathbf{D}}_{\mathrm{new}}^{*}$ generated from ${\mathbf{z}}_{\mathrm{new}}$ can be calculated using the Gaussian process:

$$p({\mathbf{D}}_{\mathrm{new}}^{*}|{\mathbf{z}}_{\mathrm{new}},\mathbf{z},\mathbf{D})=\mathcal{N}({\mathbf{D}}_{\mathrm{new}}^{*}|{\mathit{\mu }}_{\mathrm{new}},{\mathbf{\sigma }}_{\mathrm{new}}^2\mathbf{I})$$

(19)

where ${\mathit{\mu }}_{\mathrm{new}}={\mathbf{D}}_I{\mathbf{K}}_I^{-1}\mathit{k}({\mathbf{z}}_I,{\mathbf{z}}_{\mathrm{new}})$ , in which ${\mathbf{K}}_I^{-1}$ is the kernel matrix calculated from the latent variables of the active set; $\mathit{k}({\mathbf{z}}_I,{\mathbf{z}}_{\mathrm{new}})$ is a column vector with elements of kernels between each sample in the active set and the new data sample; and ${\mathbf{\sigma }}_{\mathrm{new}}^2=\mathit{k}({\mathbf{z}}_{\mathrm{new}},{\mathbf{z}}_{\mathrm{new}})-{\mathit{k}}^{\mathrm{T}}({\mathbf{z}}_I,{\mathbf{z}}_{\mathrm{new}}){\mathbf{K}}_I^{-1}\mathit{k}({\mathbf{z}}_I,{\mathbf{z}}_{\mathrm{new}})$ in which $\mathit{k}({\mathbf{z}}_{\mathrm{new}},{\mathbf{z}}_{\mathrm{new}})$ is the value of the kernel function between ${\mathbf{z}}_{\mathrm{new}}$ and itself. The likelihood of ${\mathbf{D}}_{\mathrm{new}}^{*}$ follows the Gaussian distribution, whose mean corresponds to the most probable value. Therefore, ${\mathit{\mu }}_{\mathrm{new}}$ in Equation (19) can be used as the prediction result, that is, ${\mathbf{D}}_{\mathrm{new}}^{*}={\mathit{\mu }}_{\mathrm{new}}$ . This framework allows the model to leverage the abundant historical data and the small calibration set to accurately reconstruct observations in the target domain under climate change effects.

Damage index based on reconstruction errors

Through the proposed domain-adaptive sparse Gaussian process, the model is trained on the union of historical source data and the calibration set from the target domain. It is noteworthy that both source dataset ${\mathbf{D}}_s$ and calibration set ${\mathbf{D}}_c$ are in structural health condition. This allows the model to learn the underlying physical mechanism of healthy structure while explicitly accounting for the distribution shift caused by climate change. If the structure remains healthy, the test data even under climate change conforms to the domain-invariant physical laws and the domain-specific environmental patterns learned during the calibration phase. The model can thus accurately reconstruct these observations, resulting in minimal reconstruction errors. In contrast, structural damage will alter the input–output relationship and lead to a significant increase in the reconstruction error. Accordingly, the reconstruction error serves as the damage index and is calculated as:

$${\mathit{e}}_n={\mathbf{D}}_{\mathrm{new}}-{\mathbf{D}}_{\mathrm{new}}^{*}$$

(20)

where ${\mathbf{D}}_{\mathrm{new}}$ is the measurement data, and ${\mathbf{D}}_{\mathrm{new}}^{*}$ is the prediction result of the nonlinear probabilistic latent variable model in Equation (19). The squared prediction error (SPE) statistic is defined as

$${\mathrm{SPE}}_n={\left({\mathit{e}}_n\right)}^{\mathrm{T}}{\mathit{e}}_n={({\mathbf{D}}_{\mathrm{new}}-{\mathit{\mu }}_{\mathrm{new}})}^{\mathrm{T}}({\mathbf{D}}_{\mathrm{new}}-{\mathit{\mu }}_{\mathrm{new}})+{\mathbf{\sigma }}_{\mathrm{new}}^2$$

(21)

To determine the decision threshold for the SPE statistic, the kernel density estimation (KDE) ⁴⁹ is employed. Assume that the SPE statistic follows the unknown probability density function $p(x)$ :

$$p(x)=\frac{1}{\mathrm{hL}}\sum _{n=1}^N\phi \left(\frac{x-{\mathrm{SPE}}_n}{h}\right)$$

(22)

where $h$ is the kernel width and $\phi (·)$ denotes the selected kernel function. The Gaussian kernel function $\phi \left(y\right)=\frac{e^{-y^2/2}}{\sqrt{2\pi }}$ is generally used. The control limit ${\mathrm{SPE}}_{\alpha }$ is then computed by

$${\int }_{-\infty }^{{\mathrm{SPE}}_{\alpha }}p\left(x\right)\mathrm{dx}=1-\alpha$$

(23)

where $\alpha$ is the significance level.

Comparison with baseline models

The domain-adaptative nonlinear probabilistic latent variable model and basic models discussed above are compared in Table 1. The model in this study offers a probabilistic and nonlinear framework that effectively manages uncertainty and variability in the data. Nonlinearity enables the model to capture complex relationships within the data, which is particularly important in fields where interactions between variables are intricate or influenced by multiple factors simultaneously. Additionally, the model demonstrates high efficiency and is well suitable for real-time analysis or involving large datasets, such as long-term damage detection. As highlighted in Table 1, the key advantage of the proposed method is its $\mathrm{O}(d^2·N)$ complexity. Since the computational cost scales linearly with the total number of samples N, the method remains efficient even when processing the extensive data volumes generated during long-term monitoring. This stands in stark contrast to kernel PCA, whose $\mathrm{O}\left(N^3\right)$ cost becomes prohibitive as N increases, making it unsuitable for the big data challenge of long-term SHM.

Table 1.

Comparison of model characteristics and efficiency.

Models	Probabilistic	Nonlinear	Dominant computational cost	$\mathrm{SPE}$ threshold
Linear factor analysis	Yes	No	$\mathrm{O}(N_m·M·N)$	$g·{\chi }_{h,\alpha }^2$
Kernel PCA	No	Yes	$\mathrm{O}\left(N^3\right)$	$\mathrm{KDE}$
This study	Yes	Yes	$\mathrm{O}(d^2·N)$	KDE

Note. KDE: kernel density estimation; PCA: principal component analysis. $N_m$ is original dimension of the measurement data, $M$ is the dimensions to be retained for data compression ( $M\le N_m)$ , $N$ is the number of training data samples, and $d$ is the amount of data in the active set (d < N).

To further illustrate this practical advantage in scalability, we ran a comparative test on a representative scenario with N = 3000 training samples and an active set of d = 400. The analysis was conducted on a desktop computer equipped with an Intel Core CPU, 48 GB of RAM, and an NVIDIA GeForce RTX 2080Ti GPU for acceleration. Using kernel PCA, the entire process, which includes the time-consuming trial-and-error method for selecting the optimal kernel parameter, took approximately 135 min. In contrast, the method developed in this study, which automates dimensionality determination and parameter optimization, completed the same task in approximately 9 min. This reduction to just 7% of the time required by kernel PCA underscores the proposed model’s superior efficiency and its advantage in automating parameter estimation.

Long-term damage detection method under climate change

Based on the developed domain-adaptive nonlinear probabilistic latent variable model, this study aims to propose a long-term damage detection method considering the impact of climate change, with a particular focus on whether the model can distinguish between the effects of climate change and structural damage. Therefore, the model’s training data will consist of structural vibration properties in a healthy state under varying environmental conditions. The test data, on the other hand, will include structural vibration properties under environmental conditions that fall outside the range of the training data.

To generate the test data for model performance validation, we need to predict structural vibration properties based on future climate data. The NASA NEX-GDDP-CMIP6 dataset ⁵⁰ provides global, high-resolution, bias-corrected climate projections suitable for regional impact studies. Complementing this effort, the Coupled Model Intercomparison Project Phase Six (CMIP6), initiated in 2014, provides climate projections to understand past, present, and future climate changes. This study utilizes the NEX-GDDP-CMIP6 dataset, which comprises global downscaled climate scenarios derived from the CMIP6 framework, produced by 35 different global climate models from different Institutions or Countries, as listed in Table 2.

Table 2.

Different global climate models in NEX-GDDP-CMIP6 dataset.

Model	Institution/country	Model	Institution/country
ACCESS-CM2	CSIRO-ARCCSS/Australia	HadGEM3-GC31-MM	MOHC/UK
ACCESS-ESM1-5	CSIRO/Australia	IITM-ESM	CCCR-IITM/India
BCC-CSM2-MR	BCC/China	INM-CM4-8	INM/Russia
CanESM5	CCCma/Canada	INM-CM5-0	INM/Russia
CESM2	NCAR/USA	IPSL-CM6A-LR	IPSL/France
CESM2-WACCM	NCAR/USA	KACE-1-0-G	NIMS-KMA/Republic of Korea
CMCC-CM2-SR5	CMCC/Italy	KIOST-ESM	KIOST/Republic of Korea
CMCC-ESM2	CMCC/Italy	MIROC6	MIROC/Japan
CNRM-CM6-1	CNRM-CERFACS/France	MIROC-ES2L	MIROC/Japan
CNRM-ESM2-1	CNRM-CERFACS/France	MPI-ESM1-2-HR	MPI-M, DWD, DKRZ/Germany
EC-Earth3	EC-Earth-Consortium/European	MPI-ESM1-2-LR	MPI-M, AWI, DKRZ, DWD/Germany
EC-Earth3-Veg-LR	EC-Earth-Consortium/European	MRI-ESM2-0	MRI/Japan
FGOALS-g3	CAS/China	NESM3	NUIST/China
GFDL-CM4	NOAA-GFDL/USA	NorESM2-LM	NCC/Norway
GFDL-CM4_gr2	NOAA-GFDL/USA	NorESM2-MM	NCC/Norway
GFDL-ESM4	NOAA-GFDL/USA	TaiESM1	AS-RCEC/Taiwan, China
GISS-E2-1-G	NASA-GISS/USA	UKESM1-0-LL	MOHC/UK
HadGEM3-GC31-LL	MOHC, NERC/UK

The NEX-GDDP-CMIP6 dataset encompasses historical data from 1995 to 2014 and future projections from 2015 to 2100, with a spatial resolution of 0.25°  × 0.25° (approximately 27.75 × 27.75 km at the equator, decreasing poleward). The dataset includes projections across four “Tier 1” greenhouse gas emissions scenarios, namely, SSP126, SSP245, SSP370, and SSP585, each representing a different trajectory of greenhouse gas emissions and socioeconomic development. Specifically, these represent a trajectory scale ranging from SSP126 (stringent emission mitigation and green technologies) through intermediate scenarios (SSP245, SSP370) up to SSP585 (the highest emission scenario with unabated fossil fuel use).

Utilizing the NEX-GDDP-CMIP6 dataset, structural vibration properties under anticipated future climate conditions can be forecasted. These forecasts serve as test data to evaluate the model’s performance. Crucially, climate change is treated here as a driver of long-term distribution shift, not as a coupled physical deterioration model. Figure 1 illustrates the steps of the developed damage detection method, with the highlighted pink block specifically emphasizing the generation of test data through these predictions. It is noteworthy that in future real applications to structural damage detection under climate change, these test data are directly available, and the steps outlined in the pink block can be bypassed.

Figure 1.

Flowchart of the developed damage detection method.

The detailed procedures in Figure 1 are as follows:

1.  Extract structural dynamic properties collected in the undamaged state from the entire dataset as the training data.2.  Use the training data to estimate unknown variables and parameters in the nonlinear probabilistic latent variable model. Select the active set, estimate the parameters of the kernel function automatically. Apply the feature augmentation process for domain adaptation.3.  Reconstruct or regenerate the training data based on the estimated model parameters. Calculate the SPE monitoring statistic and the corresponding threshold given the significance level.4.  Estimate structural dynamic properties in the future warming climate based on the NASA NEX-GDDP-CMIP6 dataset. Each model in Table 2 is used individually.5.  Apply the estimated model to the test dataset. Calculate the SPE statistics and determine the damage condition by comparing these values against the pre-established threshold.

Case study 1: A reinforced concrete slab

Structural description

The reinforced concrete (RC) slab, as detailed in Refs. 51, has dimensions of 6400 × 800 × 100 mm and two equal spans of 3000 mm with an additional 200 mm at each end, as depicted in Figure 2. The structure was monitored from June 2003 to March 2005, during which 136 sets of modal properties and environmental factors including temperature and humidity were recorded. The previous study showed that the RC slab’s frequencies had a linear relation with temperature and humidity, while the mode shapes were insensitive to temperature and humidity variations. Figure 2(b) plots the first four frequencies plotted against the sample numbers.

Figure 2.

Near 2-year measured frequencies of the RC slab. (a) RC slab and (b) Frequency variations from June 2003 to March 2005.

Figure 3 shows the linear relation of frequencies to temperature. Although Figure 3 exhibits a certain degree of scattering, this variability is attributed to measurement noise and the coupled influence of humidity, as detailed in the prior long-term monitoring study of this slab. ⁵¹ Based on the estimated linear relations between RC slab’s frequencies and temperature and humidity, ⁵¹ structural frequencies under the impact of future climate change can be predicted using NASA’s climate change data including humidity and temperature.

Figure 3.

Relation of frequencies to air temperature. (a) Mode 1, (b) Mode 2, (c) Mode 3, and (d) Mode 4.

NASA temperature and humidity data and structural frequency projections

The NEX-GDDP-CMIP6 dataset (0.25° × 0.25° resolution) was utilized. The nearest downscaling point (31.8750° S, 115.8750° E) was selected to represent the RC slab’s location at The University of Western Australia (31.9789° S, 115.8179° E). The daily near-surface temperature and relative humidity data of the point were then downloaded, involving 35 different global climate models listed in Table 2. Table 3 shows the predicted near-surface temperature and relative humidity data of ACCESS-ESM1-5 and EC-Earth3 models under four SSP scenarios in Year 2040, 2070, and 2100.

Table 3.

Predicted temperature (°C) and humidity (%) of different climate models for RC slab.

Climate models	Scenario	Factors	Year 2040		Year 2070		Year 2100
			Lowest	Highest	Lowest	Highest	Lowest	Highest
ACCESS-ESM1-5	SSP126	Temperature	9.61	37.12	11.36	36.34	9.43	33.98
		Humidity	26.95	93.43	25.36	91.73	30.81	88.53
	SSP245	Temperature	9.81	36.33	9.89	35.28	12.11	35.81
		Humidity	29.67	89.96	25.40	90.41	25.51	96.97
	SSP370	Temperature	10.51	33.24	10.85	36.36	13.64	39.36
		Humidity	30.69	90.60	28.25	91.55	27.70	86.70
	SSP585	Temperature	9.70	32.87	11.32	34.81	12.83	39.245
		Humidity	29.61	91.15	22.88	85.92	21.57	90.17
EC-Earth3	SSP126	Temperature	7.60	36.93	9.42	36.94	9.59	33.89
		Humidity	26.08	91.51	26.12	93.80	20.52	94.73
	SSP245	Temperature	8.35	34.44	11.26	33.53	11.09	35.14
		Humidity	24.81	91.08	23.84	91.38	27.69	90.13
	SSP370	Temperature	10.36	33.69	10.45	37.23	11.93	36.25
		Humidity	26.61	92.48	27.06	92.20	22.35	94.03
	SSP585	Temperature	9.97	34.06	10.32	33.84	12.28	40.68
		Humidity	30.89	93.40	30.40	93.48	29.98	93.20

The NEX-GDDP-CMIP6 dataset is utilized to forecast future observations of the RC slab’s frequencies, based on the fitted linear relationship between environmental factors and frequencies measured from June 2003 to March 2005. Notably, as this relationship is built on the air temperature near the RC slab, whereas the NEX-GDDP-CMIP6 dataset only provides a coarse spatial resolution of 0.25° × 0.25°, an additional projection is necessary to map the air temperature data provided by the NEX-GDDP-CMIP6 dataset to the specific air temperature near the local RC slab. Quantile mapping bias correction ⁵² is adopted to improve the alignment of climate model outputs with observations in Year 2004. As the measured dataset and climate dataset in 2014 have different sizes, interpolation methods (such as linear interpolation) are used to estimate the quantiles of the reference dataset at the cumulative probabilities corresponding to the quantiles of the target dataset. A mapping function will be developed to relate each quantile of the target dataset to the corresponding quantile of the reference dataset. Subsequently, for each value in the target dataset, the mapping function will be used to find the corresponding value in the reference dataset. Based on the bias-corrected climate data, corresponding structural frequencies can be estimated. Figure 4 illustrates the predicted frequencies in different years under different SSPs of the ACCESS-ESM1-5 climate model.

Figure 4.

Predicted frequencies of the RC slab of ACCESS-ESM1-5 model. (a) Predicted frequencies of SSP245, (b) Predicted frequencies of SSP585, (c) SSP245 in Year 2040, (d) SSP245 in Year 2100, (e) SSP585 in Year 2040, and (f) SSP585 in Year 2100.

The training data, as illustrated in Figure 2(b), encompass a range of frequencies, such as the first frequency in the range of (17.33, 18.26) Hz. However, the predicted frequencies in Figure 4 obviously extend beyond the range of these measured frequencies, such as the frequencies in Figure 4(b) with the range of (17.07, 18.02) Hz. This discrepancy highlights an out-of-distribution problem, where the data goes beyond the range of training data. To conduct the domain adaptation introduced in Equations (13) and (14), the predicted frequencies of the first ten years (i.e., from Year 2015 to 2024) are set to be the calibration set ${\mathbf{D}}_c$ , and the rest datasets (i.e., from Year 2025 to 2100) are used as the test data ${\mathbf{D}}_t$ . The model developed in this study is trained on both training and calibration data to learn domain-invariant patterns and make calibration-informed domain-specific adjustments. Subsequently, the model’s performance is validated on the test data.

Long-term structural damage detection

The developed damage detection method is applied to the RC slab dataset for damage detection. In this case study, the near 2-year monitored structural frequencies will be used as the training data. The frequencies predicted under the impact of climate change are then used as the test data to evaluate the model performance.

In the model training stage, parameters of the kernel functions need to be initialized first. The RBF kernel in Equation (7) is employed, and the initial parameters are ${\theta }_{\mathrm{rbf}}=\gamma =1$ and ${\theta }_{\mathrm{bias}}={\theta }_{\mathrm{white}}=\exp (-1)$ . As the total number of data samples ( $N=136)$ is not significant in this example, the active set d in the IVM algorithm is directly set to the number of data samples. M is set to 2. Kernel parameters are optimized in an iterative manner based on the scaled conjugate gradient approach. ⁴² When all data in the healthy state are used as training data, the optimized parameters are ${\theta }_{\mathrm{rbf}}=8.6735$ , $\gamma =0.5136$ , ${\theta }_{\mathrm{bias}}=0.1349$ and ${\theta }_{\mathrm{white}}=1.5\times {10}^{-3}$ . With optimized parameters, the measurement data are regenerated using the Gaussian process. The SPE statistics are then calculated according to Equation (21). The SPEs of the training and test data in different years under different SSPs are plotted in Figure 5 for ACCESS-ESM1-5 model and Figure 6 for EC-Earth3 model. The threshold is computed with the confidence limit of 99%, that is, $\alpha =0.01$ . Due to the space limitation, only SSP245 and SSP585 scenarios are illustrated here.

Figure 5.

Reconstruction error of ACCESS-ESM1-5 model. (a) SSP245 and (b) SSP585.

Figure 6.

Reconstruction error of EC-Earth3 model. (a) SSP245 and (b) SSP585.

As Figures 5 and 6 show, the SPEs of the test data remain comparable with those of the training data (i.e., healthy data), indicating that the test data reflect a healthy condition, and the variations in structural frequencies are due to changes in environmental conditions rather than alterations in the structural conditions. The results accord with the reality that the test data are predicted in structural health condition under climate change. It is noteworthy that under SSP585, where the environmental conditions change most significantly and fall outside the range of the training data, the developed method still demonstrates good extrapolation performance.

The comparison study

To evaluate the performance and highlight the advantages of the proposed domain-adaptive nonlinear model, we conduct a direct comparison with the traditional linear factor analysis method,^7,44 a widely used baseline for data-driven damage detection. The linear factor analysis model is trained based on the measured 136 data samples and then applied to the test data. The SPEs of the training and test datasets for the ACCESS-ESM1-5 model are plotted in Figure 7. The dash line corresponds to the threshold at the significance level $\alpha =0.01.$ The results show that the test data with the frequency variations uncovered in the training data exceed the threshold and are falsely identified as damaged state. This false-positive damage identification becomes more pronounced for SSP 585, the more severe climate change scenario. Similarly, Year 2100 has more false detection cases than Year 2040 as the climate data in the former deviates more than the latter. when the climate data deviates further from the training data. This false-positive damage identification is avoided in Figure 5 by using the method developed in this study. The results demonstrate the superior performance of the developed method in accurately distinguishing between environmental and structural changes under the impact of climate change.

Figure 7.

Reconstruction error of ACCESS-ESM1-5 model using the linear factor analysis. (a) Year 2040 of SSP245, (b) Year 2100 of SSP245, (c) Year 2040 of SSP585, (d) Year 2100 of SSP585.

Case study 2: Z24 bridge

Bridge description

Z24 bridge is a post-tensioned concrete box-girder bridge consisting of a 30 m long main span and two 14 m side spans, as depicted in Figure 8. The bridge, located in Switzerland connecting Utzenstorf and Koppigen, was monitored from 11 November 1997 to 11 September 1998. ⁹ A total of 49 sensors were installed on the bridge to record the variations of environmental conditions, including the temperature, wind, humidity, etc. Another 16 accelerometers were installed to record acceleration responses. In the later monitoring period, progressive damages were artificially introduced to the bridge in a controlled manner, beginning with pier settlement on 10 August 1998, ⁵³ as detailed in Table 4.

Figure 8.

Side view of the Z24 bridge.

Table 4.

Progressive damage test to Z24 bridge. ⁵³

Date (1998)	Scenario
4 August	Undamaged condition
9 August	Installation of pier settlement system
10 August	Lowering of pier, 20 mm
12 August	Lowering of pier, 40 mm
17 August	Lowering of pier, 80 mm
18 August	Lowering of pier, 95 mm
19 August	Lifting of pier, tilt of foundation
20 August	New reference condition
25 August	Spalling of concrete at soffit, 12 m2
26 August	Spalling of concrete at soffit, 24 m2
27 August	Landslide of 1 m at abutment
31 August	Failure of concrete hinge
2 September	Failure of 2 anchor heads
3 September	Failure of 4 anchor heads
7 September	Rupture of 2 out of 16 tendons
8 September	Rupture of 4 out of 16 tendons
9 September	Rupture of 6 out of 16 tendons

The bridge’s responses under the ambient excitations were recorded, and the stochastic subspace identification method ⁵⁴ was used to extract the modal parameters. A total of 5624 sets of the first four frequencies were obtained during the entire 304 monitoring days. The variations of the first four frequencies and the temperature of the bridge over time are plotted in Figure 9. The blue data represent the healthy condition, where the purple ones correspond to the damaged state. The first damage was artificially introduced around day 266, corresponding to data No. 4789. The relation between the air temperature and natural frequencies is also plotted in Figure 10, which reflects that the natural frequencies have a bilinear (piecewise linear) relation with the temperature. The bilinear effect is attributed to the asphalt layer on the bridge’s surface, which froze on cold days (lower than 0°C) and significantly increased the stiffness of the structure. ⁵⁵ This well-documented bilinear effect, physically attributed to the asphalt layer freezing, provides a critical test for any damage detection model. A model that cannot capture this nonlinear physical phenomenon is likely to fail to detect the damage. The proposed method’s ability to model this specific nonlinearity will be demonstrated in the following sections.

Figure 9.

Frequency and temperature variations over time. (a) Variation of frequencies and (b) Variation of the air temperature.

Figure 10.

The relation between the frequency and temperature. (a) Second frequency with temperature and (b) Fourth frequency with temperature.

In this case study, the frequency data in the healthy state are divided into training and calibration sets for domain adaptation in Equations (13) and (14). As shown in Figure 9(b), temperatures recorded between days 166 and 265 exceed the range of the initial 165-day period. Consequently, data from the first 165 days are used for source training dataset ${\mathbf{D}}_s$ , while data from days 166 to 265 serve as the calibration set ${\mathbf{D}}_c$ . The bridge’s performance in the damaged state serves as the test data ${\mathbf{D}}_t$ , including the measured damaged data and the predicted damaged data under a future warming climate. Since damaged data are only available for temperatures above 0°C, and the measurement data for the final damage scenario are limited, predicting the damaged bridge’s frequencies in the future is challenging due to insufficient data. To address this, a Gaussian process is employed to model the relationship (green line) between frequency and temperature in the healthy state, as shown in Figure 11. This model is then fine-tuned to align with the damaged measurement data. The method proposed by Figueiredo et al. ³⁷ is adopted to generate the test data. Specifically, using the derived relationship between structural frequencies and temperature in the final damaged state, the bridge’s frequencies under future climate change impacts are predicted and used as test data.

Figure 11.

The fitted relation between the frequency and air temperature in the healthy state. (a) First frequency and (b) Second frequency.

NASA temperature data and bridge frequency projections

The Z24 bridge connects Koppigen (47.1340° N, 7.6002° E) and Utzenstorf (47.1277° N, 7.5596° E). The nearest NEX-GDDP-CMIP6 downscaling point is located at 47.1250° N, 7.6250° E. The daily near-surface air temperature data of the point are downloaded from 35 different global climate models. Table 5 presents the predicted and downscaled near-surface temperatures of the bridge from the ACCESS-ESM1-5 and EC-Earth3 models. Figure 12 further illustrates the specific temperature variations in certain years, highlighting the differences in temperature trends under different scenarios. The temperature projections under the SSP 585 scenario exhibit a more pronounced increasing trend compared to those under SSP 126.

Table 5.

Predicted temperature (°C) data of the Z24 bridge using different models.

Models	Scenario	Year 2040		Year 2070		Year 2100
		Lowest	Highest	Lowest	Highest	Lowest	Highest
ACCESS-ESM1-5	SSP126	− 5.17	26.10	− 4.90	24.78	− 7.55	26.53
	SSP245	− 5.53	27.20	− 5.16	30.06	− 3.65	28.53
	SSP370	− 5.98	27.47	− 4.95	30.63	− 5.64	31.27
	SSP585	− 9.21	31.11	− 5.39	29.35	− 4.63	31.62
EC-Earth3	SSP126	− 5.98	27.08	− 2.42	26.96	− 5.39	25.63
	SSP245	− 5.73	27.52	− 6.22	29.73	− 0.13	26.39
	SSP370	1.08	27.04	− 6.30	31.00	− 8.03	29.73
	SSP585	− 8.50	26.02	1.79	27.32	1.20	38.58

Figure 12.

Downscaled temperature data from ACCESS-ESM1-5 model. (a) Downscaled temperature of SSP126, (b) Downscaled temperature of SSP585, (c) Year 2040 of SSP126, (d) Year 2100 of SSP126, (e) Year 2040 of SSP585, and (f) Year 2100 of SSP585.

The temperature dataset is then used to project future observations of natural frequencies of Z24 bridge, based on the fitted relationship between temperature and frequencies. The quantile mapping bias correction method mentioned in Case study 1 is employed to correct the bias of datasets. Figure 13 illustrates the predicted second frequency of the bridge under SSP 126 and SSP 585, based on the ACCESS-ESM1-5 climate model. Compared with Figure 10, the predicted frequencies do not obviously exceed the established range of the training data, because the upper frequency limit was established by training data under temperature below 0°C, whereas the climate change leads to warmer test data and yields lower frequencies. Nevertheless, the test data exhibit the distributional shift (e.g., mean shift) that clusters around new warmer conditions. Notably, the sudden spikes in the predicted frequencies shown in Figure 13(c)–(f), particularly prominent during the winter months (the beginning and end of the 1-year period), correspond to days when the temperature drops below 0°C. This reflects the physical bilinear effect caused by the freezing of the asphalt layer, significantly increasing the structural stiffness.

Figure 13.

Predicted bridge frequency (test data) using ACCESS-ESM1-5 model. (a) Frequencies prediction of SSP126, (b) Frequencies prediction of SSP585, (c) Year 2040 of SSP126, (d) Year 2100 of SSP126, (e) Year 2040 of SSP585, and (f) Year 2100 of SSP585.

Long-term bridge damage detection

All SHM data under the healthy condition, specifically the first 4788 data points, are used as the training data. The RBF kernel in Equation (7) is employed, and the initialized kernel parameters remain the same as Case study 1. The active set d in the IVM algorithm is set to $d=$ 400 in this example. The latent dimension M is set to $2$ . After training, parameters of the kernel function are optimized to ${\theta }_{\mathrm{rbf}}=1.1046$ , $\gamma =1.3556$ , ${\theta }_{\mathrm{bias}}=0.0301$ , and ${\theta }_{\mathrm{white}}=2.4832\times {10}^{-4}$ when all healthy data are used as training data. The trained model is then applied to the test data. It is noteworthy that the predicted test data for the future years (2040, 2070, and 2100) represent the bridge in its damaged condition under varying climate scenarios. With the estimated parameters, the measurement data are regenerated using the sparse Gaussian process. The SPE statistics under four SSPs of the ACCESS-ESM1-5 model are then calculated using the developed method and plotted in Figure 14. The dash line corresponds to the threshold at the significance level $\alpha =0.01.$ All calculations in this study are carried out on a desktop with a CPU of Intel Core and 48 GB RAM. The computation of SPE statistics for each figure takes around 5 min.

Figure 14.

Reconstruction error of ACCESS-ESM1-5 model using the proposed method. (a) SSP126, (b) SSP245, (c) SSP370, and (d) SSP585.

As illustrated in Figure 14, the damage index in the damaged state (i.e., red dots) can be clearly distinguished from that in the healthy state, even for the slight damage at the initial stage. As damage progresses across more locations, the damage index rises and significantly surpasses the threshold. The magnitude of the damage index also indicates the severity of the damage. Regarding data under future climate conditions, the damaged state is accurately classified as such and can be distinctly separated from the healthy state. This example underscores the robustness and reliability of the developed method in identifying structural damage at early stages with latent nonlinear relationships. The efficiency and effectiveness of the developed method in addressing damage detection under environmental nonlinear impacts and climate change is also demonstrated.

The comparison study

For comparison, the classical linear factor analysis method is also applied. As illustrated in Figure 15, linear factor analysis fails to clearly separate most damaged states from the healthy condition. This leads to dangerous false-negative detections, particularly during the early stages of progressive damage ⁵³ or under climate-induced distribution shifts, allowing structural issues to worsen undetected. In contrast, the method developed in this study effectively addresses this challenge and clearly identifies the damaged condition in Figure 14. The developed method is capable of accurately detecting even minor damage, reducing the risk of false negatives.

Figure 15.

Reconstruction error of ACCESS-ESM1-5 model using the linear factor analysis. (a) SSP126, (b) SSP245, (c) SSP370, and (d) SSP585.

Conclusions

This study develops a novel domain-adaptive nonlinear probabilistic latent variable framework for long-term damage detection, specifically designed to address distribution shifts caused by climate change. By optimizing existing ML models to handle out-of-distribution problems, the research significantly enhances the robustness and accuracy of damage detection results. The introduction of a nonlinear latent variable model, combined with the Bayesian technique and kernel trick, allows for the modeling of complex physical relationships, such as the nonlinear temperature effects observed in the Z24 bridge. By building upon a transparent mathematical framework rather than a black-box approach, the model successfully distinguishes damage from these intricate environmental patterns. Furthermore, the feature augmentation method employed in this study expands the feature space rather than forcing a complex distribution projection. Hence, the model avoids the sensitivity to initial parameters often associated with traditional transfer learning, ensuring consistent performance across different scenarios. The model’s ability to operate without explicit environmental data further extends its utility to scenarios where sensor coverage is limited. The underlying principle that separating statistical environmental anomalies from structural changes is applicable to a wide range of civil infrastructure subject to thermal expansion and climate variability.

The approach has been tested on two case studies under various climate change projections. In the RC slab example, where test data exceeded the training temperature range (e.g., SSP585), the traditional linear factor analysis resulted in a high false-positive detection rate, with damage indices consistently exceeding the threshold. In contrast, the proposed method effectively eliminated these false positives, maintaining the damage index of healthy condition entirely within the 99% confidence threshold. In the Z24 bridge case study involving latent nonlinearities, the linear baseline suffered from significant false-negative detections that masked damage signals. The proposed framework achieved a clear separation between healthy and damaged states, successfully identifying damage scenarios even under extreme climate projections. The advancement of this method represents a significant contribution to the field of SHM, offering a robust damage detection approach to address the evolving challenges posed by climate change.

Nevertheless, climate change in this study is treated primarily as a driver of long-term distribution shift (i.e., covariate shift) rather than through a coupled physical deterioration model. It is important to acknowledge that while the latent variables provide a statistical convenience for capturing environmental variability, they represent abstract statistical dimensions and do not strictly correspond to isolated physical parameters. The proposed framework focuses on mitigating the statistical divergence caused by evolving environmental conditions to prevent false-positive and false-negative damage detection, rather than explicitly modeling material degradation phenomena such as creep, shrinkage, or corrosion. While the current algorithm can flag deviations in the structure-environment correlation, whether caused by acute damage or gradual aging, the performance relies on the premise that solving the distribution shift is a prerequisite for robust monitoring. Future research will aim to bridge this gap by integrating physics-based aging models into this domain-adaptive framework, thereby enhancing physical interpretability and enabling a more detailed classification that distinguishes between reversible environmental effects, gradual structural aging, and acute damage.