3.1. Mathematical model of structural system

The equation of motion of a linear system with a $ {n}_{\mathrm{dof}} $ degree of freedom under forced vibration can be represented by a linear second-order ordinary differential equation as

(3.1) $$ \mathbf{M}\ddot{\mathbf{u}} (t)+\mathbf{C}\dot{\mathbf{u}}(t)+\mathbf{Ku}(t)={\mathbf{S}}_p\mathbf{p}(t), $$

where $ \mathbf{u}(t)\in {\mathrm{\mathbb{R}}}^{n_{\mathrm{dof}}} $ is the vector of displacements at $ {n}_{\mathrm{dof}} $ degrees of freedom, while the matrices $ \mathbf{M} $ , $ \mathbf{C} $ , $ \mathbf{K} $ denote, respectively, the mass, damping and stiffness of the system. The external forces applied to the system are represented by $ {\mathbf{S}}_p\mathbf{p}(t) $ as a combination of loads. The Boolean input shape matrix $ {\mathbf{S}}_p\in {\mathrm{\mathbb{R}}}^{n_{\mathrm{dof}}\times {n}_{\mathrm{p}}} $ indicates the degree of freedom where the $ {p}^{th} $ load is applied, while the input vector $ \mathbf{p}(t)\in {\mathrm{\mathbb{R}}}^{n_{\mathrm{p}}} $ represents the time history of the $ {p}^{th} $ load.

When it comes to complex systems, it is a common practice to represent the equation of motion by using a reduced order model as described in (Craig Jr., Reference Craig1985; Lourens et al., Reference Lourens, Papadimitriou, Gillijns, Reynders, De Roeck and Lombaert2012a). The undamped eigenvalue problem associated with eq. (3.1) is expressed as

(3.2) $$ \mathbf{K}\boldsymbol{\Phi } =\mathbf{M}{\boldsymbol{\Phi} \boldsymbol{\Omega}}^2, $$

which yields the eigenvector $ {\Phi}_i $ as a column of the matrix $ \boldsymbol{\Phi} \in {\mathrm{\mathbb{R}}}^{n_{\mathrm{dof}}\times {n}_{\mathrm{dof}}} $ and the natural frequencies $ {\omega}_i^2 $ as diagonal of the matrix $ {\boldsymbol{\Omega}}^2\in {\mathrm{\mathbb{R}}}^{n_{\mathrm{dof}}\times {n}_{\mathrm{dof}}} $ . In the context of a reduced order model system, the dynamic behavior is captured by a reduced number of modes $ {n}_r $ and the matrices are consequently truncated into $ \overline{\boldsymbol{\Phi}}\in {\mathrm{\mathbb{R}}}^{n_{\mathrm{dof}}\times {n}_{\mathrm{r}}} $ and $ {\overline{\boldsymbol{\Omega}}}^2\in {\mathrm{\mathbb{R}}}^{n_{\mathrm{r}}\times {n}_{\mathrm{r}}} $ . By premultiplying each term by $ {\overline{\boldsymbol{\Phi}}}^{\mathrm{T}} $ and transforming $ \mathbf{u}(t)=\overline{\boldsymbol{\Phi}}\mathbf{r}(t) $ , eq. (3.1) becomes

(3.3) $$ {\overline{\boldsymbol{\Phi}}}^{\mathrm{T}}\mathbf{M}\overline{\boldsymbol{\Phi}}\mathbf{r}\ddot (t)+{\overline{\boldsymbol{\Phi}}}^{\mathrm{T}}\mathbf{C}\overline{\boldsymbol{\Phi}}\dot{\mathbf{r}}(t)+{\overline{\boldsymbol{\Phi}}}^{\mathrm{T}}\mathbf{K}\overline{\boldsymbol{\Phi}}\mathbf{r}(t)={\overline{\boldsymbol{\Phi}}}^{\mathrm{T}}{\mathbf{S}}_p\mathbf{p}(t). $$

The equation of motion can be decoupled by leveraging the orthogonality conditions corresponding to a set of mass-normalized eigenvectors, such that $ {\overline{\boldsymbol{\Phi}}}^{\mathrm{T}}\mathbf{M}\overline{\boldsymbol{\Phi}}=\mathbf{I} $ and $ {\overline{\boldsymbol{\Phi}}}^{\mathrm{T}}\mathbf{K}\overline{\boldsymbol{\Phi}}={\overline{\boldsymbol{\Omega}}}^2 $ . When assuming proportional damping, the modal damping matrix $ {\overline{\boldsymbol{\Phi}}}^{\mathrm{T}}\mathbf{C}\overline{\boldsymbol{\Phi}}=\overline{\boldsymbol{\Gamma}} $ is also diagonal. The matrix $ \overline{\boldsymbol{\Gamma}} $ includes the terms $ 2{\xi}_j{\omega}_j $ for $ j=1,\dots, {n}_r $ , where $ {\xi}_j $ is the modal damping ratio, $ {\omega}_j $ is the natural frequency, and $ {n}_r $ is the number of retained modes.

Consequently, the governing equation of motion in modal coordinates can be written as

(3.4) $$ \mathbf{r}\ddot (t)+\overline{\boldsymbol{\Gamma}}\dot{\mathbf{r}}(t)+{\overline{\boldsymbol{\Omega}}}^2\mathbf{r}(t)={\overline{\boldsymbol{\Phi}}}^{\mathrm{T}}{\mathbf{S}}_p\mathbf{p}(t). $$

Considering the forces transformed into their modal contributions by $ \mathbf{f}\left(\mathbf{t}\right)={\overline{\boldsymbol{\Phi}}}^{\mathrm{T}}{\mathbf{S}}_{\mathbf{p}}\mathbf{p}\left(\mathbf{t}\right) $ , where $ \mathbf{f}\left(\mathbf{t}\right)\in {\mathrm{\mathbb{R}}}^{n_{\mathrm{r}}} $ and defining the state vector $ \mathbf{x}\left(\mathbf{t}\right)\in {\mathrm{\mathbb{R}}}^{n_{\mathrm{s}}\times {n}_{\mathrm{s}}} $ as

(3.5) $$ \mathbf{x}(t)=\left[\begin{array}{c}\mathbf{r}(t)\\ {}\dot{\mathbf{r}}(t)\end{array}\right], $$

where $ {n}_s=2{n}_{\mathrm{r}} $ , the continuous-time second-order differential equation, detailed in eq. (3.4), can be written as a set of first-order continuous-time differential equations. These are so-called continuous-time state-space equations and are expressed as

(3.6) $$ \dot{\mathbf{x}}(t)={\mathbf{A}}_c\mathbf{x}(t)+{\mathbf{B}}_c\mathbf{f}(t), $$

where the system matrices $ {\mathbf{A}}_{\mathbf{c}}\in {\mathrm{\mathbb{R}}}^{n_s\times {n}_s} $ and $ {\mathbf{B}}_{\mathbf{c}}\in {\mathrm{\mathbb{R}}}^{n_s\times {n}_r} $ are defined as

(3.7) $$ {\mathbf{A}}_c=\left[\begin{array}{cc}\mathbf{0}& \mathbf{I}\\ {}-{\overline{\boldsymbol{\Omega}}}^2& -\overline{\boldsymbol{\Gamma}}\end{array}\right],\hskip1em {\mathbf{B}}_c=\left[\begin{array}{c}\mathbf{0}\\ {}\mathbf{I}\end{array}\right], $$

When it comes to SHM systems, it is necessary to define a measurement model to describe the relationship between the state vector and the observed variables. These observed variables can be a linear combination of displacement, velocity, and acceleration and be arranged in a vector $ \mathbf{y}(t)\in {\mathrm{\mathbb{R}}}^{n_y} $ where $ {n}_{\mathrm{y}} $ represents the number of measurements. The data vector $ \mathbf{y}(t) $ becomes:

(3.8) $$ \mathbf{y}(t)={\mathbf{S}}_a\mathbf{u}\ddot (t)+{\mathbf{S}}_v\dot{\mathbf{u}}(t)+{\mathbf{S}}_d\mathbf{u}(t), $$

where the selection matrices $ {\mathbf{S}}_a\in {\mathrm{\mathbb{R}}}^{n_a\times {n}_{\mathrm{dof}}},{\mathbf{S}}_v\in {\mathrm{\mathbb{R}}}^{n_v\times {n}_{\mathrm{dof}}},{\mathbf{S}}_d\in {\mathrm{\mathbb{R}}}^{n_d\times {n}_{\mathrm{dof}}} $ are populated according to the location where acceleration, velocity and displacement (strains) are measured, respectively.

The observation equation, expressed in eq. (3.8), can be rewritten by utilizing the definition of the state vector $ \mathbf{x}(t) $ as

(3.9) $$ \mathbf{y}(t)={\mathbf{G}}_c\mathbf{x}(t)+{\mathbf{J}}_c\mathbf{f}(t), $$

where the output influence matrix $ {\mathbf{G}}_c $ and the direct transmission matrix $ {\mathbf{J}}_c $ are expressed as

(3.10) $$ {\mathbf{G}}_{\mathbf{c}}=\left[{\mathbf{S}}_d\overline{\boldsymbol{\Phi}}-{\mathbf{S}}_a{\overline{\boldsymbol{\Phi}}\overline{\boldsymbol{\Omega}}}^2\hskip0.5em {\mathbf{S}}_v\overline{\boldsymbol{\Phi}}-{\mathbf{S}}_a\overline{\boldsymbol{\Phi}}\overline{\boldsymbol{\Gamma}}\right],\hskip1em {\mathbf{J}}_{\mathbf{c}}=\left[\begin{array}{c}{\mathbf{S}}_a\overline{\boldsymbol{\Phi}}\end{array}\right]. $$

Equations (3.6) and (3.8) define the reduced order model of the continuous time-space representation of the system described in eq. (3.1). However, measurement systems are not able to capture continuous-time observations $ \mathbf{y}(t) $ in reality; instead, they sample the response of the system at discrete time intervals. Consequently, the continuous-time state space model is discretized in time intervals defined at $ {t}_k=k\Delta t $ where $ k=1,\dots, N $ and expressed as

(3.11) $$ {\mathbf{x}}_{k+1}={\mathbf{Ax}}_k+{\mathbf{Bf}}_k+{\mathbf{w}}_{\mathbf{k}}, $$

(3.12) $$ {\mathbf{y}}_k={\mathbf{Gx}}_k+{\mathbf{Jf}}_k+{\eta}_{\mathbf{k}}, $$

where $ \mathbf{A}={e}^{{\mathbf{A}}_c\Delta t},\hskip1em \mathbf{B}=\left[\mathbf{A}-\mathbf{I}\right]{\mathbf{A}}_c^{-1}{\mathbf{B}}_c,\hskip1em \mathbf{G}={\mathbf{G}}_c,\hskip1em \mathbf{J}={\mathbf{J}}_c. $ The terms $ {\mathbf{w}}_k $ and $ {\boldsymbol{\eta}}_k $ account for structural process and measurement noise, respectively, and are modeled as zero-mean, white Gaussian noises. They are assumed mutually uncorrelated, i.e., $ \unicode{x1D53C}\left\{{\mathbf{w}}_k{\boldsymbol{\eta}}_k^{\top}\right\}=\mathbf{0} $ , with covariances $ {\mathbf{Q}}^x=\unicode{x1D53C}\left\{{\mathbf{w}}_k{\mathbf{w}}_k^{\top}\right\} $ and $ \mathbf{R}=\unicode{x1D53C}\left\{{\boldsymbol{\eta}}_k{\boldsymbol{\eta}}_k^{\top}\right\} $ . The structural model uncertainty, $ {\mathbf{w}}_k $ , is specified in continuous time through a white-noise driving process with spectral density $ {\mathbf{Q}}_c^x $ . The corresponding discrete-time covariance $ {\mathbf{Q}}^x $ is obtained with the latent forces through the discretisation of the augmented state-space model, as detailed in Section 3.3. It is worth noting that, since the structural modelling error is approximated as a zero-mean white process, possible systematic or coloured components of the model discrepancy are not explicitly included.

3.2. Gaussian process regression in state-space model

GP regression is a non-parametric Bayesian regression method that provides a robust and probabilistic approach to modeling data (Williams and Rasmussen, Reference Williams and Rasmussen1995; Rasmussen, Reference Rasmussen2004). A GP is defined as a collection of random variables, any subset of which follows a joint Gaussian distribution. The GP is defined as

(3.13) $$ f(t): GP\left(m(t),k\left(t,{t}^{\prime },\theta \right)\right). $$

where $ m(t) $ represents the expected value of the function, and $ k\left(t,{t}^{\prime },\theta \right) $ , also known as the kernel, defines the covariance between two points $ t $ and $ {t}^{\prime } $ of the process. The properties of the GP, including smoothness, periodicity, and the scale of variation, are determined by the kernel. The hyperparameters $ \theta $ of the kernel control these properties, enabling the GP to adapt to the data’s characteristics.

GPs with certain kernels that satisfy the Markov properties can be reformulated within a state-space model (SSM) framework (Hartikainen and Särkkä, Reference Hartikainen and Särkkä2010; Särkkä et al., Reference Särkkä, Solin and Hartikainen2013). This formulation replaces the conventional kernel-based representation with a dynamical system, where the GP is converted into a linear time-invariant stochastic differential equation (SDE) of order $ \beta $ :

(3.14) $$ \frac{d^{\beta }f(t)}{dt^{\beta }}+{a}_{\beta -1}\frac{d^{\beta -1}f(t)}{dt^{\beta -1}}+\cdots +{a}_0f(t)=w(t), $$

where $ f(t)\sim \mathcal{GP}\left(0,k\left(t,{t}^{\prime}\right)\right) $ is a GP with zero-mean and stationary covariance function $ k\left(t,{t}^{\prime}\right) $ , $ w(t) $ is zero-mean continuous-time white noise with spectral density $ {S}_w(w)={q}_w $ , and $ \mathbf{a}(t)={\left[\begin{array}{c}{a}_{\beta -1},{a}_{\beta -2},\dots, {a}_0\end{array}\right]}^{\top } $ is a set of coefficients. The function $ f(t) $ and its derivatives $ \beta -1 $ can be represented as a vector $ \mathbf{z}(t)={\left[\begin{array}{c}f(t),\frac{df(t)}{dt},\dots, \frac{d^{\beta -1}f(t)}{dt^{\beta -1}}\end{array}\right]}^{\top } $ and used to build a first-order vector Markov process as

(3.15) $$ {\displaystyle \begin{array}{l}\dot{\mathbf{z}}(t)={\mathbf{F}}_c\mathbf{z}(t)+{\mathbf{L}}_c\mathbf{w}(t),\\ {}\mathrm{f}(t)={\mathbf{H}}_c\mathbf{z}(t),\end{array}} $$

where the matrices $ {\mathbf{F}}_c\in {\mathrm{\mathbb{R}}}^{n_{\beta}\times {n}_{\beta }},{\mathbf{L}}_c\in {\mathrm{\mathbb{R}}}^{n_{\beta}\times 1},{\mathbf{H}}_c\in {\mathrm{\mathbb{R}}}^{1\times {n}_{\beta }} $ are given by

(3.16) $$ {\mathbf{F}}_c=\left[\begin{array}{cccc}0& 1& \cdots & 0\\ {}\vdots & \vdots & \ddots & \vdots \\ {}0& 0& \cdots & 1\\ {}-{a}_0& -{a}_1& \cdots & -{a}_{\beta -1}\end{array}\right],\hskip1em {\mathbf{L}}_c=\left[\begin{array}{c}0\\ {}\vdots \\ {}0\\ {}1\end{array}\right],\hskip1em {\mathbf{H}}_c=\left[\begin{array}{cccc}1& \cdots & 0& 0\end{array}\right]. $$

Solving the SDE in eq. (3.15) provides a complete probabilistic characterization of the state vector $ \mathbf{z}(t) $ , which is modeled as a Gaussian random variable with a mean $ \hat{\mathbf{z}}(t) $ and covariance $ \mathbf{P} $ , such that $ \mathbf{z}(t)\sim \mathcal{N}\left(\hat{\mathbf{z}}(t),\mathbf{P}\right) $ . For the specific case of linear time-invariant (LTI) SDEs with initial conditions $ \mathbf{z}\left({t}_0\right)\sim \mathcal{N}\left({\hat{\mathbf{z}}}_0,{\mathbf{P}}_0\right) $ , the mean and covariance of $ \mathbf{z}(t) $ are defined by the following differential equations:

(3.17a) $$ \frac{d\hat{\mathbf{z}}(t)}{dt}={\mathbf{F}}_c\hat{\mathbf{z}}(t), $$

(3.17b) $$ \frac{d\mathbf{P}(t)}{dt}={\mathbf{F}}_c\mathbf{P}(t)+\mathbf{P}(t){\mathbf{F}}_c^{\mathrm{T}}+{\mathbf{Q}}_c^f, $$

where the term $ {\mathbf{Q}}_c^f={\mathbf{L}}_c{q}_w{\mathbf{L}}_c^{\mathrm{T}} $ represents the process noise covariance, which quantifies the uncertainty introduced by the white noise $ w(t) $ in the system.

For a zero-mean stationary process, the solution of the SDE reaches a steady-state form, where the state vector follows a constant Gaussian distribution, $ {\mathbf{z}}_k\sim \mathcal{N}\left(\mathbf{0},{\mathbf{P}}_{\infty}\right) $ . Under these conditions, eq. (3.17b) takes the form of the Lyapunov equation:

(3.18) $$ 0={\mathbf{F}}_c{\mathbf{P}}_{\infty }+{\mathbf{P}}_{\infty }{\mathbf{F}}_c^{\mathrm{T}}+{\mathbf{Q}}_c^f. $$

The solution of eq. (3.18), which is a special case of Riccati equations, is usually computed using numerical algorithms to obtain the steady state covariance matrix $ {\mathbf{P}}_{\infty } $ (Grewal and Andrews, Reference Grewal and Andrews2014; Särkkä and Solin, Reference Särkkä and Solin2019).

In order to use the Bayes’ filter and perform GP regression, eq. (3.15) needs to be discretized. By leveraging the Markov property, data is processed sequentially, allowing the posterior distribution of the state vector $ {\mathbf{z}}_k $ to be estimated at each time step $ k=0,\Delta t,\dots, T $ . The discretized model is defined as

(3.19a) $$ {\mathbf{z}}_k=\mathbf{F}\hskip0.1em {\mathbf{z}}_{k-1}+\mathbf{L}\hskip0.1em {w}_{k-1}, $$

(3.19b) $$ {f}_k=\mathbf{H}\hskip0.1em {\mathbf{z}}_k, $$

where $ \mathbf{F}=\exp \left({\mathbf{F}}_c\Delta t\right) $ , $ \mathbf{L}=\left(\mathbf{F}-\mathbf{I}\right){\mathbf{F}}_c^{-1}{\mathbf{L}}_c $ and $ \mathbf{H}={\mathbf{H}}_c $ .

The state-space representation in eq. (3.19) provides a structured framework for Gaussian Process regression within a dynamic system. Specifically, eq. (3.19a) defines the dynamics model, which governs the evolution of the latent states over time. On the other hand, eq. (3.19b) serves as the measurement model, relating samples $ y $ to the states. Equation (3.19) enables the update of the filtered posterior distribution of state vector $ {\mathbf{z}}_k $ at the time step $ k $ based on the filtered posterior distribution computed at the previous time step $ k-1 $ . This process is carried out using forward filtering methods, such as the Kalman filter (Kalman, Reference Kalman1960). However, while the Kalman filter is optimal in a forward-prediction sense, it does not account for future observations that may further refine the estimates (Särkkä, Reference Särkkä2008). To address this limitation, a backward smoother, such as the RTS smoother, can be applied (Rauch et al., Reference Rauch, Tung and Striebel1965). The RTS smoother computes the smoothing distribution of $ \mathbf{z} $ , using future measurements to enhance state estimates retrospectively.

The application of filters and smoothers requires defining the covariance $ {\mathbf{Q}}^f $ of the discretized zero-mean Gaussian white noise $ {w}_{k-1} $ , which can be expressed as

(3.20) $$ {\mathbf{Q}}^f={\int}_{t_k}^{t_{k+1}}\exp \left({\mathbf{F}}_c\left(\Delta t-\tau \right)\right){\mathbf{Q}}_c^f\exp {\left({\mathbf{F}}_c\left(\Delta t-\tau \right)\right)}^{\mathrm{T}} d\tau . $$

For a zero-mean stationary process, numerical algorithms can be used to retrieve the process noise covariance $ {\mathbf{Q}}^f $ by solving the Lyapunov equation given by

(3.21) $$ {\mathbf{Q}}^f={\mathbf{P}}_{\infty }-{\mathbf{F}\mathbf{P}}_{\infty }{\mathbf{F}}^{\mathrm{T}}. $$

3.3. Gaussian process latent force model

The GPLFM provides a probabilistic framework for estimating latent (unknown) forces acting on dynamical systems. Traditional SSMs rely on observed inputs to describe system dynamics. However, in many real-world applications, these inputs are either unknown or exhibit high complexity, making direct measurement or modeling impractical. The GPLFM addresses this issue by modeling latent forces as GPs and integrating them into a state-space formulation (Alvarez et al., Reference Alvarez, Luengo and Lawrence2009). This hybrid approach takes advantage of both the mechanical system dynamics and the stochastic nature of the latent forces to infer the evolution of both system states and unknown inputs over time.

To achieve this, the state vector of the mechanical system, $ \mathbf{x}\left(\mathbf{t}\right) $ in eq. (3.6), is augmented with the input state vector $ \mathbf{z} $ , which is derived from the GP representation in eq. (3.15). This results in an augmented state-space model, formulated as

(3.22) $$ {\displaystyle \begin{array}{c}\left[\begin{array}{c}\dot{\mathbf{x}}(t)\\ {}{\dot{\mathbf{z}}}^{\left(\mathbf{1}\right)}(t)\\ {}{\dot{\mathbf{z}}}^{\left(\mathbf{2}\right)}(t)\\ {}\vdots \\ {}{\dot{\mathbf{z}}}^{\left({\mathbf{n}}_{\mathbf{r}}\right)}(t)\end{array}\right]=\underset{{\mathbf{F}}_c^a}{\underbrace{\left[\begin{array}{ccccc}{\mathbf{A}}_c& {\mathbf{b}}_1{\mathbf{H}}_c^{(1)}& {\mathbf{b}}_2{\mathbf{H}}_c^{(2)}& \cdots & {\mathbf{b}}_{n_r}{\mathbf{H}}_c^{\left({n}_r\right)}\\ {}0& {\mathbf{F}}_c^{(1)}& 0& \cdots & 0\\ {}0& 0& {\mathbf{F}}_c^{(2)}& \cdots & 0\\ {}\vdots & \vdots & \vdots & \ddots & \vdots \\ {}0& 0& 0& \cdots & {\mathbf{F}}_c^{\left({n}_r\right)}\end{array}\right]}}\left[\begin{array}{c}\mathbf{x}(t)\\ {}{\mathbf{z}}^{\left(\mathbf{1}\right)}(t)\\ {}{\mathbf{z}}^{\left(\mathbf{2}\right)}(t)\\ {}\vdots \\ {}{\mathbf{z}}^{\left({\mathbf{n}}_{\mathbf{r}}\right)}(t)\end{array}\right]+\underset{{\mathbf{w}}^a(t)}{\underbrace{\left[\begin{array}{c}{\mathbf{w}}^{\mathbf{x}}\left(\mathbf{t}\right)\\ {}{\tilde{\mathbf{w}}}_{GP}^{(1)}(t)\\ {}{\tilde{\mathbf{w}}}_{GP}^{(2)}(t)\\ {}\vdots \\ {}{\tilde{\mathbf{w}}}_{GP}^{\left({n}_r\right)}(t)\end{array}\right]}},\\ {}\mathbf{y}(t)=\underset{{\mathbf{H}}_c^a}{\underbrace{\left[\begin{array}{ccccc}{\mathbf{G}}_c& {\mathbf{j}}_1{\mathbf{H}}_c^{(1)}& {\mathbf{j}}_2{\mathbf{H}}_c^{(2)}& \cdots & {\mathbf{j}}_{n_r}{\mathbf{H}}_c^{\left({n}_r\right)}\end{array}\right]}}\left[\begin{array}{c}\mathbf{x}(t)\\ {}{\mathbf{z}}^{\left(\mathbf{1}\right)}(t)\\ {}{\mathbf{z}}^{\left(\mathbf{2}\right)}(t)\\ {}\vdots \\ {}{\mathbf{z}}^{\left({\mathbf{n}}_{\mathbf{r}}\right)}(t)\end{array}\right]+\mathbf{v}(t),\end{array}} $$

where $ {\mathbf{b}}_1,{\mathbf{b}}_2,\dots, {\mathbf{b}}_{n_r} $ are the columns of the structural matrix $ {\mathbf{B}}_c $ , $ {\mathbf{j}}_1,{\mathbf{j}}_2,\dots, {\mathbf{j}}_{n_r} $ are the columns of the structural matrix $ {\mathbf{J}}_c $ , and $ \mathbf{v} $ represents the measurement noise modeled as white Gaussian noise. The term $ {\mathbf{w}}^a(t) $ represents the augmented continuous process noise with spectral density $ {\mathbf{Q}}_c^a $ . Assuming the structural model error and the latent force driving noise are uncorrelated in continuous time, the spectral density matrix $ {\mathbf{Q}}_c^a $ is block-diagonal:

(3.23) $$ {\mathbf{Q}}_c^a=\left[\begin{array}{cc}{\mathbf{Q}}_c^x& \mathbf{0}\\ {}\mathbf{0}& {\mathbf{Q}}_c^f\end{array}\right], $$

where $ {\mathbf{Q}}_c^x $ is the spectral density of the structural model error and $ {\mathbf{Q}}_c^f $ contains the spectral densities of the GP latent forces. For each $ j $ -th latent force, the noise vector is defined as $ {\tilde{\mathbf{w}}}_{GP}^{(j)}(t)={\mathbf{L}}_c^{(j)}{w}^{(j)}(t) $ , where $ {w}^{(j)}(t) $ is the scalar white noise driving the SDE. Consequently, $ {\mathbf{Q}}_c^f $ is constructed as a block-diagonal matrix where the $ j $ -th block is given by

(3.24) $$ {\mathbf{Q}}_c^{(j)}={\mathbf{L}}_c^{(j)}{q}_j{\left({\mathbf{L}}_c^{(j)}\right)}^{\bullet }, $$

where $ {q}_j $ is the spectral density derived from the kernel hyperparameters (e.g., variance and length-scale) of the $ j $ th latent force.

To facilitate practical implementation, the continuous-time augmented system is discretized over finite time intervals $ \Delta t $ . In shorthand, the discrete augmented SSM is expressed as

(3.25) $$ {\mathbf{s}}_k^a={\mathbf{F}}^a{\mathbf{s}}_{k-1}^a+{\mathbf{q}}_{k-1}^a, $$

(3.26) $$ {\mathbf{y}}_k={\mathbf{H}}^a{\mathbf{s}}_k^a+{\mathbf{v}}_k, $$

where $ {\mathbf{s}}_k^a={\left[{\mathbf{x}}_k^{\mathrm{T}},{\mathbf{z}}_k^{\mathrm{T}}\right]}^{\mathrm{T}} $ is the augmented state vector and $ {\mathbf{F}}^a=\exp \left({\mathbf{F}}_c^a\Delta t\right) $ is the discretized transition matrix. The discrete process noise vector $ {\mathbf{q}}_{k-1}^a $ aggregates the stochastic contributions from both the structure and the latent forces. Its covariance matrix $ {\mathbf{Q}}^a=\unicode{x1D53C}\left[{\mathbf{q}}_{k-1}^a{\left({\mathbf{q}}_{k-1}^a\right)}^{\mathrm{T}}\right] $ is computed by solving the continuous Lyapunov equation for the steady-state covariance $ {\mathbf{P}}_{\infty}^a $ of the full augmented system:

(3.27) $$ {\mathbf{F}}_c^a{\mathbf{P}}_{\infty}^a+{\mathbf{P}}_{\infty}^a{\left({\mathbf{F}}_c^a\right)}^{\top }+{\mathbf{Q}}_c^a=\mathbf{0}. $$

Using the computed steady-state covariance $ {\mathbf{P}}_{\infty}^a $ , the discrete process noise covariance is obtained as

(3.28) $$ {\mathbf{Q}}^a={\mathbf{P}}_{\infty}^a-{\mathbf{F}}^a{\mathbf{P}}_{\infty}^a{\left({\mathbf{F}}^a\right)}^{\mathrm{T}}. $$

Once the system is discretized, its posterior distribution can be estimated using the Kalman filter at each time step, based on the current and past measurements. To further refine the state and latent force estimates, a backward smoothing step is applied using the RTS smoother.

3.4. Choice of the covariance function and hyperparameters inference

An important step in the use of GPs involves selecting the covariance function or kernel. The kernel determines how values at different points relate to each other, influencing the smoothness, periodicity, and behavior of the GPs. To build an SSM of a GP, specific conditions regarding the kernel must be satisfied: the kernel must be positive semi-definite and should enable the GP to have a spectral density $ S(w) $ in rational form. Once these criteria are met, the state-space matrices ( $ \mathbf{F},\mathbf{L},\mathrm{and}\;{q}_w $ ) can be derived through spectral factorisation, as described in (Hartikainen and Särkkä, Reference Hartikainen and Särkkä2010).

In structural engineering applications, isotropic covariance functions, such as the periodic, squared exponential, exponential, and Matérn class, are usually employed (Vettori et al., Reference Vettori, Di Lorenzo, Peeters and Chatzi2024). The Matérn family of covariance functions is often the preferred choice in situations where there is limited prior knowledge about the loading conditions due to their versatility in modeling different types of random processes. Vettori et al. (Reference Vettori, Di Lorenzo, Peeters and Chatzi2024) explored the inner dynamic features of the most common covariance functions and validated their performance on an experimental case study involving a 3D-printed scaled titanium WT blade under several loading conditions. Vettori et al. show that when the system is subjected to an impulse, or non-smooth load condition, the exponential function tends to perform best due to its ability to capture sharp variations. However, Stein (Reference Stein2012) argues that exponential kernels may impose unrealistic smoothness assumptions when applied to measured data from operating structures, potentially leading to inaccurate estimations. As an alternative, the Matérn 3/2 or 5/2 kernels are often better suited. In this work, given the limited knowledge about the excitation, the Matérn 3/2 kernel is adopted, and its state-space representation is provided in Appendix A.

Once the kernel is chosen and its SSM representation is formulated, attention shifts to inferring its hyperparameters, which directly impact the quality of the GP model and the resulting GPLFM estimates. In the case of the Matérn 3/2 kernel, the hyperparameters govern the scale, length scale, and smoothness of the latent force estimations. A common approach is to tune them by using dynamic response measurements of the dynamical system. Ghahramani and Hinton (Reference Ghahramani and Hinton1996) proposed a methodology to select the optimal hyperparameters maximizing the marginal likelihood of the data and improving the model’s fit to the data. This method was integrated into a state-space model for dynamic systems by Nayek et al. (Reference Nayek, Chakraborty and Narasimhan2019). The method consisted of minimizing the residual norm between observed data and model predictions based on the least squares solution. However, it has been observed to be susceptible to false estimations and not able to find the global maximum (Petersen et al., Reference Petersen, Øiseth and Lourens2022). To address this, various robust methodologies have been proposed over the years. Rogers et al. (Reference Rogers, Worden and Cross2020) adopted a Bayesian framework with Monte Carlo simulations to obtain the full posterior distribution. Although the methodology guarantees convergence to the true posterior distribution of interest, it is computationally expensive (Marino and Cicirello, Reference Marino and Cicirello2023; Zou et al., Reference Zou, Lourens and Cicirello2023). To reduce the computational cost of a switching GPLFM used to identify discontinuous forces in dynamical systems, Marino and Cicirello (Reference Marino and Cicirello2023) introduced an approach that combines variational inference with active sampling Bayesian quadrature and used it in the context of mechanical systems with discontinuous linearities, the so-called Variational Bayes Monte Carlo (Acerbi, Reference Acerbi2018, Reference Acerbi2020).

Zou et al. (Reference Zou, Lourens and Cicirello2023) proposed a low-complexity and effective approach that is utilised in this study. This technique takes advantage of prior assumptions encoded in the GP kernel to align the empirical distribution of measurements, denoted as $ l=p(y)\sim \mathcal{N}\left(\mathbf{0},{\mathbf{P}}^y\right) $ , with the theoretical probability distribution of the Gaussian prior, defined as $ {l}^{\theta }=p\left(y|\theta \right)\sim \mathcal{N}\left(\mathbf{0},{\hat{\mathbf{P}}}_{\infty}^y\right) $ . Specifically, the hyperparameter inference is achieved by minimising the Hellinger distance between these two distributions, conducted as follows:

(3.29) $$ {\theta}^{\ast }=\arg \underset{\theta }{\min }D\left(l,{l}^{\theta}\right), $$

where $ D\left(l,{l}^{\theta}\right) $ represents the Hellinger distance. The empirical covariance of the measurements $ {\mathbf{P}}^y $ is computed directly from the time series of sensor observations. The theoretical covariance $ {\hat{\mathbf{P}}}_{\infty}^y $ is derived from the steady-state covariance of the augmented model, $ {\mathbf{P}}_{\infty } $ . For this inference step, the structural process noise is assumed to be negligible ( $ {\mathbf{Q}}_c^x=\mathbf{0} $ ). Therefore, $ {\mathbf{P}}_{\infty } $ is computed by solving the Lyapunov equation (eq. (3.18)) using only the spectral density of the latent forces, $ {\mathbf{Q}}_c^f $ . This ensures that the optimization of $ \theta $ is driven solely by the properties of the GP kernel. This transformation to observation space is carried out using the matrix $ {\mathbf{H}}^a $ of eq. (3.26) as

(3.30) $$ {\hat{\mathbf{P}}}_{\infty}^y={\mathbf{H}}^a{\mathbf{P}}_{\infty }{\left({\mathbf{H}}^a\right)}^{\mathrm{T}}. $$

3.5. Inference of the noise covariance matrices

The process noise covariance matrix $ {\mathbf{Q}}^x $ represents the uncertainty in the assumed dynamical model of the structure, primarily due to modeling errors. Lower (or higher) values of $ {\mathbf{Q}}^x $ indicate more (or less) confidence in the dynamic model. In contrast, the input process noise covariance matrix $ {\mathbf{Q}}^f $ represents the covariance of a hypothetical noise process that governs the evolution of the input estimates. As such, $ {\mathbf{Q}}^f $ is not independently inferred or tuned, but analytically calculated based on the chosen GP kernel and its corresponding hyperparameters.

In this study, the process noise covariance matrix $ {\mathbf{Q}}^x $ was quantified by tuning the continuous-time power spectral density matrix $ {\mathbf{Q}}_c^x $ . This matrix was treated as a tuning parameter and was optimised by minimising the Normalized Root Mean Square Error (NRMSE), whose theoretical background is explained in Appendix C, between the posterior mean structural responses and the recorded measurements. Specifically, the optimization problem is formulated as

(3.31) $$ {\mathbf{Q}}_c^{x\ast }=\arg \min {\mathbf{Q}}_c^x\mathrm{NRMSE}\left({\mathbf{y}}^m,\mathbf{y}\right), $$

where $ {\mathbf{y}}^m $ denotes the measured outputs and $ \mathbf{y} $ represents the posterior mean of the predicted measurements. The optimized $ {\mathbf{Q}}_c^{x\ast } $ is embedded in the augmented continuous-time spectral density $ {\mathbf{Q}}_c^a $ and discretized through the augmented Lyapunov-based procedure (Section 3.3), yielding the discrete augmented covariance $ {\mathbf{Q}}^a $ .

The measurement noise covariance matrix $ \mathbf{R} $ was determined directly from the experimental data collected during a resting state of the structure, prior to the impact event. Specifically, the variance of the acceleration signals was computed for each sensor during a calm period in which the structural response is assumed to be dominated by sensor noise. The matrix $ \mathbf{R} $ was subsequently constructed as a diagonal matrix, assuming that the measurement noise is uncorrelated across different sensors.

4.1. Description of the Magerholm ferry quay

The Magerholm ferry quay became a benchmark facility for SHM research through a collaboration between Møre og Romsdal county and NTNU in 2023 (Sibille et al., Reference Sibille, Nord, Frøseth, Fenerci and Siedziako2025). The quay ranks as one of the busiest ferry quays in Norway, transporting approximately 950,000 vehicles across the Storfjorden fjord in 2023. Two ferries operate on this route, each measuring 111 meters in length, with a gross tonnage of 3,000, and capable of transporting up to 120 cars. The operational schedule varies throughout the week, with 55 crossings per day on weekdays and 48 crossings on weekends. The frequency of crossings changes throughout the day, with the longest break occurring at night (approximately two hours). During peak hours, departures take place every 20 minutes.

The quay, depicted in Figure 1, consists of two independent structures. The first is a concrete pier approximately 110 meters long, supported by concrete-filled steel piles. It is fitted with thirteen fenders arranged to absorb the forces generated during ferry docking. The second structure is a steel linkspan designed for loading and unloading vehicles and passengers. The linkspan measures approximately 15 meters in length and 9 meters in width, composed of five HE550B longitudinal beams and IPE240 cross beams. Due to its exposure to frequent ferry impacts and significant wear, the front beam is made from Hardox400 steel. The linkspan is vertically supported by a concrete abutment and two hydraulic lifting towers, which allow height adjustments to accommodate water level fluctuations. The longitudinal movements (surge direction) are absorbed by four conical fenders positioned between the steel structure and the concrete foundation. The quay also includes an electric charging station and a vacuum anchoring system designed to stabilize the ferry when docked, reducing its movements. For a comprehensive description of the quay, readers can refer to Sibille et al. (Reference Sibille, Nord, Frøseth, Fenerci and Siedziako2025).

Figure 1.

The Magerholm ferry quay.

4.2. Measurement campaign and data processing

This study focuses exclusively on the linkspan, while the general layout of the Magerholm ferry quay is outlined for context. The linkspan was selected due to its higher susceptibility to structural degradation, more demanding maintenance routines, and notably greater dynamic responses compared to the concrete pier.

Response data were collected on the linkspan during two distinct ferry impact events, using two different data acquisition configurations shown in Figure 2a and Figure 2b. The instrumentation consisted of G-Link-200 3-axis accelerometers, housed in waterproof enclosures provided by Parker LORD MicroStrain, and Linear Variable Differential Transformers (LVDTs) supplied by Load Indicator System AB. The accelerometers’ range was configured to $ \pm $ 2g with a noise density of $ 25\mu g/\sqrt{\mathrm{Hz}} $ . The LVDTs had a nominal stroke of 100 mm. For both configurations, all sensors transmitted data to a laptop through a WSDA-2000 gateway via a wireless connection. To reduce low-frequency noise and potential drift, all recorded acceleration signals analyzed in Section 5 underwent high-pass filtering with a cutoff frequency of 0.5 Hz to attenuate low-frequency energy that cannot be reliably captured due to sensor sensitivity limitations. The displacement data were not filtered to preserve the quasi-static component of the linkspan’s response.

Figure 2.

Sensor positions on the linkspan: (a) sensor configuration 1; (b) sensor configuration 2.

The monitoring campaign was split into two distinct acquisition setups due to the data logger’s capacity limitations, which restricted the employment of many sensors at high sampling frequencies. One setup focused on maximizing spatial information by using a larger number of accelerometers sampling at a moderate frequency (128 Hz), while the other prioritized temporal resolution by capturing acceleration and displacement responses at a higher sampling rate (1024 Hz), allowing the assessment of the GPLFM under measurements with increased high-frequency content. This approach enabled a comprehensive assessment of how sensor positions, types (accelerometers and LVDTs), and sampling frequencies influence the performance of virtual sensing. In both setups, accelerometer data were recorded only along the longitudinal axis (aligned with the direction of the expected impact) to maximize the number of operating sensors and sampling frequency within the data logger’s constraints.

The two data acquisition setups were as follows:

• • Configuration 1 (Impact event 1): this configuration involved fourteen G-Link-200 accelerometers, positioned as shown in Figure 2a. These sensors were set up to sample at 128 Hz.

• • Configuration 2 (Impact event 2): this configuration included five G-Link-200 accelerometers and two LVDTs, with their locations depicted in Figure 2b. These sensors recorded data at a sampling frequency of 1024 Hz. The LVDTs were positioned to measure the relative movement of the linkspan with respect to the concrete abutment.

The complete docking can be categorized into four distinct phases: (1) adjustment of the linkspan prior to impact, (2) ferry impact, (3) vehicle passage, and (4) undocking. The unfiltered acceleration response recorded by sensor S1 in configuration 1 during these four phases is shown in Figure 3. The adjustment phase occurs as the ferry approaches the quay. Because the clearance between the ferry and the quay varies with the tidal level, the linkspan must be manually repositioned to align with the ferry ramp. An operator onboard the ferry performs this adjustment to ensure that the impact occurs at a suitable location on the linkspan. This phase is reflected in the acceleration signals as changes in the mean value since the gravitational component is captured by the sensors. The second phase is generally short in duration, as it typically lasts only a few seconds. However, it results in the highest acceleration peaks due to the impact between the ferry and the structure. The third phase involves the unloading and loading of passengers and vehicles. This phase is significantly longer, lasting several minutes (approximately 300 seconds in this case). The acceleration response during this period is more moderate, though high accelerations can be observed when heavy vehicles, such as trucks, pass near the sensors. Finally, the ferry undocking phase marks the return of the linkspan to its initial, unloaded configuration. This phase is typically characterized by a gradual change in acceleration levels as the ferry undocks and the linkspan readjusts its position.

Figure 3.

Full time series recorded from sensor S1 in configuration 1 during a complete docking event.

The distinct operational phases, particularly the impact and vehicle passage, introduce significant complexities for system modeling. During adjustment and undocking, the changing inclination of the linkspan occurs. While variations in the stiffness of structural elements like conical fenders or lifting towers due to these inclination changes might be considered minor, the more dominant challenge arises from the interaction between the linkspan and the ferry. This interaction alters the linkspan’s support configuration. This creates time-varying boundary conditions and couples the linkspan’s dynamics with the ferry’s movements, which are influenced by unpredictable environmental factors (e.g., waves, wind) and operational variables (e.g., engine thrust; Sibille et al., Reference Sibille, Nord, Frøseth, Fenerci and Siedziako2025, Reference Sibille, Nord, Le, Siedziako and Cicirello2026). Given these complexities, and to provide a focused evaluation of the GPLFM capabilities for virtual sensing under real conditions of unknown impact loading and inherent model uncertainties, this study concentrates its analysis exclusively on the ferry impact phase. For this specific time window, the structural system of the linkspan is assumed to behave as a Linear Time-Invariant (LTI) system. This simplification is adopted because precisely modeling the rapidly evolving contact mechanics during the impact event (e.g., determining whether it behaves as a perfect impact, quantifying a coefficient of restitution, or timing the exact moment of impact) presents considerable challenges. A detailed analysis of the physics of the impact is beyond the scope of this paper, which focuses on virtual sensing. A comprehensive analysis and characterization of the ferry impacts exerted on the same case study can be found in Sibille et al. (Reference Sibille, Nord, Frøseth, Fenerci and Siedziako2025).

4.3. FE model

In this work, the system matrices $ {\mathbf{A}}_c $ and $ {\mathbf{B}}_c $ used in eq. (3.22), representative of the system’s dynamics and its relationship to external inputs, are retrieved from a FE model of the quay. The FE model is developed in Abaqus based on construction drawings and field measurements. The model includes only structural elements, such as HE550B longitudinal beams, IPE240, HE240A, and UNP220 secondary beams, all of which are modeled using S235 and S275 steel. The mesh is generated using quadrilateral finite membrane-strain shell elements. Additionally, the mass of non-structural elements, including railings and gratings, is incorporated, resulting in a total modeled structural weight of 53.64 tons. The conical fenders are modeled as longitudinal springs with a stiffness of 2000 kN/m. The lifting towers are represented as vertical springs with a stiffness set at 10000 kN/m (Siedziako et al., Reference Siedziako, Nord and Fenerci2024; Le et al., Reference Le, Siedziako, Nord and Sibille2026). The four conical fenders of the linkspan (see Figure 1) are designed to engage sequentially as a function of the imposed displacement. For the impact events analysed in this study, only two of the four fenders were activated, as the recorded displacements remained below the threshold required to engage the remaining two fenders. These events correspond to typical operational conditions of the ferry quay, characterised by relatively smooth berthing manoeuvres. Under more severe impact scenarios, larger displacements would activate the additional fenders, thereby increasing the effective stiffness of the system in the longitudinal direction and enhancing its energy dissipation capacity. Accordingly, the present study focuses on typical operational events, for which the system can be approximated as linear and supported by two fenders, which are therefore used in the construction of the FE model. The model is updated using a sensitivity-based FE model updating method with modal parameters obtained from operational modal analysis (Le et al., Reference Le, Siedziako, Nord and Sibille2026). During the optimization process, the following parameters are defined as design variables: steel density, Young’s modulus, and the stiffness of the two lifting towers.

To facilitate computational efficiency, a reduced-order model is constructed by retaining the first seven most influential modes identified from the FE Model. These modes, whose mode shapes are presented in the Appendix B, correspond to a frequency range extending up to 10.10 Hz, as modal analysis conducted on experimental campaigns indicates that the dominant horizontal, bending, and torsional modes are contained within this frequency range. Figure 4 illustrates the modal participation at each sensor location for these seven modes, specifically for sensor configuration 1 (Figure 4a with fourteen sensors S1–S14) and sensor configuration 2 (Figure 4b with five sensors A1–A5). The color intensity in the plots represents the magnitude of the product of the sensor selection matrix ( $ {\mathbf{S}}_{\mathbf{a}} $ ) and the mode shape matrix ( $ \boldsymbol{\Phi} $ ), indicating the sensitivity of each sensor to each mode. The figure shows that the first two modes, which are predominantly horizontal, exhibit the highest influence across sensor locations in both configurations. This is consistent with the accelerometers recording primarily along the longitudinal axis. Among the other modes, which represent vertical/torsional behaviors, the third and seventh modes show the lowest influence on the sensors.

Figure 4.

Modal influence of accelerations for: (a) sensor configuration 1; (b) sensor configuration 2.

Furthermore, due to the absence of experimental modal damping estimates in the horizontal direction, a uniform modal damping ratio of 5% is assumed for all modes included in the reduced order model. This damping ratio is used to build the modal damping matrix $ \overline{\boldsymbol{\Gamma}} $ in eq. (3.4).

4.4. Challenges to data and model assumptions

The application of the GPLFM to the Magerholm ferry quay for virtual sensing faces various challenges. These arise from: (i) uncertainties in accurately representing the structural model, (ii) the unknown characteristics of the impact force, (iii) the time-variability of the system, especially during its interaction with the ferry, and (iv) spatially sparse measurements. The following points detail how each key challenge is accounted for:

• • Modeling uncertainties: The structure incorporates numerous steel components connected through joints, which are subjected to a harsh maritime environment. This exposure causes corrosion, altering material properties and joint behavior in ways that are difficult to quantify precisely. These uncertainties are partially mitigated by performing FE model updating using parameters derived from OMA data (Le et al., Reference Le, Siedziako, Nord and Sibille2026).

• • Unknown impact force characteristics: The unknown magnitude, duration, and specific time history of the impact force are directly handled by modeling the force as a GP latent force within the GPLFM. The unknown location and direction of the impact force are addressed by decomposing the latent forces $ \mathbf{f}\left(\mathbf{t}\right) $ into its modal components. This modeling approach involves two critical assumptions. First, the modal components of the forces are assumed to be uncorrelated, which simplifies the estimation process but may not always be realistic when dealing with experimental data, as structural systems often exhibit coupling between dynamic modes. Second, all modal components are modeled using the same GP prior, characterized by identical kernels and hyperparameters. Although the same state-space matrices are employed for all the latent forces and calculated using the same kernel and hyperparameters, the posterior variance of the states of the modal components of the latent forces is computed via Kalman filter and smoothing. Zou et al. (Reference Zou, Lourens and Cicirello2023) demonstrated that this assumption does not affect the estimation of the posterior estimates of the latent states and the performance of the GPLFM, as they have assigned different hyperparameters for each latent force and assessed the accuracy of the GPLFM.

• • System time-variability and boundary conditions: The system experiences considerable time-variability due to the alteration of the support configuration. As illustrated in Figure 5, the linkspan switches from a ‘disconnected state’ (Figure 5a), where it is supported by the abutment (A) and the hydraulic cylinders (C), to a ‘docked state’ (Figure 5b), where the linkspan rests on the ferry. This shift in boundary conditions not only alters the physical supports but also dynamically couples the linkspan’s behavior with the ferry’s movements, which are influenced by variable environmental conditions (e.g., wave action, wind) and operational inputs (e.g., engine thrust; Sibille et al., Reference Sibille, Nord, Frøseth, Fenerci and Siedziako2025, Reference Sibille, Nord, Le, Siedziako and Cicirello2026). These transitions alter the modal properties, such as natural frequencies, damping ratios, and mode shapes, of the linkspan within short timeframes. For the scope of this study, the system is assumed to be an LTI during the analyzed time window. While this is acknowledged as a strong assumption, it allows for avoiding introducing a switching system that would entail its own set of challenging unknowns, as the estimation of the modal properties for the different states within the impact and of the precise timing of these switches, which are beyond the primary focus of this study. Moreover, modal damping ratios are assumed quantities due to a lack of direct operational identification for all relevant modes/directions. Consequently, the sensitivity of the GPLFM’s virtual sensing performance to variations in the assumed damping ratio is evaluated to understand the potential impact of this uncertainty in Section 5.5.

• • Spatially sparse measurements and uni-axial data: The data acquisition setup comprises a limited number of sensors, leading to spatially sparse measurements across the structure. Moreover, the available sensors record acceleration data exclusively along the longitudinal direction (x-axis) of the linkspan. This uni-axial measurement constraint implies that dynamic responses in the transverse (y-axis) and vertical (z-axis) directions are not directly captured.

Figure 5.

The operational states of the linkspan: (a) disconnected state; (b) docked state.

These assumptions, along with the corresponding analyses performed in this study, are summarized in Table 1.

Table 1.

Summary of key challenges: assumptions made and performed analyses

5.1. Influence of sensor location on the accuracy of the estimated response

A critical aspect influencing the robustness of any virtual sensing methodology is the spatial relationship between the available instrumented sensors and the location where unmeasured responses are to be estimated. To explore this sensitivity within the context of the GPLFM, the accuracy of the posterior mean of the estimated acceleration response is evaluated as it varies with the target sensor location. For this investigation, sensor configuration 1, comprising fourteen accelerometers (S1–S14) as illustrated in Figure 2, was selected.

The accuracy of the posterior mean of the estimated acceleration response at each sensor location was evaluated using the root mean squared error (RMSE) metric, whose theoretical background is reported in Appendix C. The evaluation was performed using a leave-one-out cross-validation approach: for each sensor location $ {S}_i $ (where $ i $ ranges from 1 to 14) within sensor configuration 1, the measured data from that specific sensor $ {S}_i $ was temporarily excluded from the set of observations. The GPLFM was then performed using the measurement data from all other thirteen available sensors (i.e., sensors $ {S}_1,\dots, {S}_{i-1},{S}_{i+1},\dots, {S}_{14} $ ) to perform the joint input-state estimation process. Based on these measurements, the GPLFM then probabilistically estimated the acceleration response at the location of the excluded sensor $ {S}_i $ . The posterior mean of the estimated response was compared in the time domain against the measured response from $ {S}_i $ using the RMSE metric. This entire process was repeated iteratively for each of the fourteen sensor locations.

Figure 6 shows the variation in estimation accuracy across the different sensor locations. Overall, the RMSE values are generally below 0.15 m/s $ {}^2 $ for most locations. However, locations near the edges of the longitudinal beams, such as S1, S6, S9, and S14, tend to exhibit slightly higher RMSE values compared to sensors located more centrally along the beams (e.g., S2–S5, S10–S13). This may be attributed to the fact that these sensors have fewer neighboring sensors to provide spatial information. Moreover, a notable increase in RMSE is observed at locations S7 and S8, which are situated on the front beam of the linkspan.

Figure 6.

Root mean squared error (RMSE) for acceleration response estimation at each sensor location (S1–S14) using a leave-one-out approach with sensor configuration 1.

To further illustrate these performance differences, the estimated and measured acceleration responses for two representative cases are shown in the time and frequency domains. Specifically, Figure 7 shows the results for location S13 (corresponding to the lowest RMSE), while Figure 8 showcases the estimations at location S8 (corresponding to the highest RMSE).

Figure 7.

Acceleration response comparison at location S13: (a) full time history of the posterior mean estimate; (b) difference between measured response and posterior mean estimate; (c) zoomed-in time window around impact; (d) power spectral density (PSD). Black solid lines denote measurements, dashed red lines denote estimates, and the red vertical bar in the PSD plots indicates the natural frequency of the highest mode included in the reduced-order model.

Figure 8.

Acceleration response comparison at location S8: (a) full time history of the posterior mean estimate; (b) difference between measured response and posterior mean estimate; (c) zoomed-in time window around impact; (d) PSD. Black solid lines denote measurements, dashed red lines denote estimates, and the red vertical bar in the PSD plots indicates the natural frequency of the highest mode included in the reduced-order model.

As shown in Figure 7 (location S13), the GPLFM demonstrates high accuracy in reconstructing the posterior mean of the acceleration response. It is worth noticing that the associated posterior uncertainty bounds are visually negligible in the time-domain plots, as the measurement noise covariance R is small relative to the response amplitudes during the impact event. The influence of increased measurement noise on the visibility of uncertainty bounds is investigated in Section 5.3. The posterior mean of the estimated time history closely matches the measured signal and both the amplitude and timing of the peak acceleration during the ferry impact are accurately captured. Furthermore, the power spectral density (PSD) of the estimated response shows high agreement with that of the measured data across the entire analyzed frequency spectrum. In contrast, the estimation for location S8, shown in Figure 8, presents inaccuracies, particularly in correspondence with the high acceleration peaks of the impact. However, despite the high RMSE, the timing of the acceleration peak of the ferry impact is reasonably well-identified, although its amplitude is not well reconstructed. In the frequency domain, the estimated PSD for S8 aligns reasonably with the measured PSD, below approximately 10 Hz. Beyond this frequency, the estimated PSD deviates, failing to capture the energy content at higher frequencies. While an initial inaccuracy is noticeable around 10 Hz, which roughly corresponds to the upper frequency limit of the reduced-order model used in this study (retaining modes up to approximately 10.10 Hz), the discrepancy becomes particularly pronounced after approximately 20 Hz.

The observed lower accuracy in reconstructing the response at locations S7 and S8, both situated on the front beam of the linkspan, may be attributed to: (i) ill-conditioning with distance from impact/observation locations; (ii) sensor distribution across different longitudinal beams. As established in previous research on Kalman filter-based methods (Lourens et al., Reference Lourens, Reynders, De Roeck, Degrande and Lombaert2012b; Aucejo et al., Reference Aucejo, De Smet and Deü2019), the accuracy of state and input estimation can degrade when measurements are taken from sensors that are relatively far from the location of the input force, or when trying to estimate the response in a location that is poorly observed by the available sensor network. In this specific scenario, when virtual sensing is performed for S7 or S8 (located on the front beam, near the presumed impact zone), these sensors are “left out.” If the majority of the remaining active sensors are situated at a relative distance from the front beam and the area of impact, the system becomes more susceptible to ill-conditioning (Lourens et al., Reference Lourens, Reynders, De Roeck, Degrande and Lombaert2012b). Beyond this ill-conditioning, the accuracy of virtual sensing can also be influenced by the distribution of the observing sensors across different structural elements. During the leave-one-out process for S7 or S8, no other sensors remain on that same longitudinal beam. If the system connectivity between the longitudinal beams is low, the model may struggle to reconstruct the response accurately (Erkaya, Reference Erkaya2018). This implies that local beam dynamics may become challenging to estimate from distant sensors on parallel members.

5.2. Optimal real-world sensor placement for virtual sensing

A common question when instrumenting a structure for virtual sensing involves determining the optimal placement of a limited number of sensors to maximize the accuracy of response estimation in the impact direction at a specific, uninstrumented target location. This section outlines a methodology to leverage the developed DT framework to guide an Optimal Sensor Placement (OSP) strategy. The objective is to identify a subset of sensor locations that provides the most informative data for reconstructing the response at a desired target point of interest.

To address this, a numerical simulation was carried out using a sensor selection algorithm. A triangular force, shown in Figure 9, was applied at an arbitrary location on the front beam (Figure 9c of the FE model to simulate a representative ferry impact. Simulated acceleration responses were deterministically generated at fourteen predefined sensor locations, corresponding to those used in Configuration 1 of the experimental campaign (see Figure 2a). The motivation was to ensure comparability with experimental results, allowing validation of the numerical framework against real data. These sensor locations were selected prior to this study based on physical accessibility and to spatially cover the structure as uniformly as possible. No prior OSP strategies had been applied.

Figure 9.

Simulated triangular load used for the OSP numerical study: (a) full time series of the force; (b) power spectral density of the force; (c) sensor and force locations.

With these simulated responses, a Backward Sequential Sensor Placement (BSSP) strategy (Papadimitriou, Reference Papadimitriou2004) was implemented to identify the optimal sensor subsets for probabilistically estimating the response at a chosen target location. In this case, sensor location S8 was selected as the target for virtual sensing. The BSSP algorithm works iteratively:

1. 1. Initially, all potential sensor locations (excluding the target S8, so thirteen sensors) are considered. The GPLFM is used to estimate the posterior mean response at S8 using these thirteen simulated sensor inputs, and the RMSE between this estimate and the “true” simulated response at S8 is calculated.

2. 2. One sensor is removed at a time from the current set. For each possible removal, the GPLFM is re-run with the reduced sensor set to estimate the posterior mean response at S8, and the corresponding RMSE against the true simulated response is calculated.

3. 3. The sensor whose removal results in the highest RMSE value is permanently removed from the set. This means the removed sensor was the least critical for estimating the response at S8 given the remaining sensors.

4. 4. This process is repeated, reducing the number of active sensors by one at each iteration, until a predefined minimum number of sensors is reached, in this case, three sensors.

Figure 10a showcases the resulting RMSE values for virtual sensing at location S8 as the number of optimally selected sensors varies from 13 down to 3. As expected, Figure 10a illustrates that decreasing the number of sensors leads to an increase in the RMSE, reflecting a loss of spatial information available to the GPLFM for reconstructing the response. Notably, the RMSE remains relatively low and similar for sensor set sizes from 13 down to approximately five–six sensors. This suggests that beyond five–six optimally placed sensors, adding more sensors yields small returns in terms of improving the estimation accuracy at S8 for this specific simulated load case. However, a more significant increase in RMSE is observed when the number of sensors drops below five, particularly with only three or four sensors, indicating that these smaller sets lack sufficient spatial information to accurately estimate the posterior mean response at S8.

Figure 10.

(a) Root mean squared error (RMSE) for estimated acceleration response at location S8 as a function of the number of optimally placed sensors; (b) five optimal sensor locations (in red) used to estimate the response at location S8 (in blue).

Further, the OSP approach is validated with real-world data of configuration 1. The optimal sensor configuration identified in the numerical study for five sensors was selected. These five optimal sensor locations were identified as S1, S2, S6, S9, and S14 (see Figure 10b). The GPLFM was then applied using the measured experimental data from these five selected sensors to estimate the posterior mean of the acceleration response at location S8. Figure 11 shows the estimated posterior mean of the response at position S8 against the measured one in the time and frequency domains. The Sprague and Geers (S&G) metric, whose theoretical background is explained in Appendix C, was used as a metric to evaluate the accuracy of the estimated response in the time domain. This metric jointly accounts for amplitude and phase discrepancies between the estimated and measured signals and provides a normalized, physically interpretable measure of reconstruction quality. For the virtual sensing of location S8 using five optimally placed sensors, the combined S&G metric yields a value of 0.55. When using thirteen sensors, the S&G metric computed between the posterior mean response and the measured signal yields a value of 0.58, indicating that comparable reconstruction accuracy can be achieved with a smaller number of optimally selected sensors.

Figure 11.

Acceleration response comparison at location S8 using five optimally positioned sensors: (a) full time history of the posterior mean estimate; (b) difference between measured response and posterior mean estimate; (c) zoomed-in time window around impact; (d) PSD. Black solid lines denote measurements, dashed red lines denote estimates, and the red vertical bar in the PSD plots indicates the natural frequency of the highest mode included in the reduced-order model.

Observing the time histories in Figure 11, the posterior mean of the estimated acceleration response obtained using the five optimal sensors closely replicates the timing of the peaks seen in the measured response, although the magnitudes differ, similarly to what was observed with thirteen sensors (Figure 8). However, there are clear differences in the frequency domain. The PSD of the response estimated with five sensors begins to deviate significantly from the measured PSD after approximately 6-7 Hz. In contrast, when 13 sensors were used, this deviation occurred at around 10-20 Hz. This suggests that while the five optimally chosen sensors capture lower frequency content reasonably well, they may lack the spatial resolution to accurately reconstruct some of the higher frequency modes.

5.3. Sensitivity to measurement noise matrix R

The posterior covariance of the estimated states and responses is directly influenced by the assumed measurement noise covariance matrix $ \mathbf{R} $ within the GPLFM framework. In particular, lower values of $ \mathbf{R} $ lead to increased confidence in the measurements and, consequently, to narrower posterior uncertainty bounds, whereas higher values of $ \mathbf{R} $ result in larger posterior variances and wider confidence intervals.

In this study, the measurement noise covariance matrix $ \mathbf{R} $ was estimated directly from experimental data recorded during the resting (non-impact) condition of the structure. For configuration 1, the variance of each sensor signal was computed individually. These sensor-specific variances were then used to populate the diagonal entries of $ \mathbf{R} $ under the assumption of uncorrelated measurement noise. The resulting average acceleration variance of the diagonal entries of $ \mathbf{R} $ was approximately $ {10}^{-6}\hskip0.1em {\left(\mathrm{m}/{\mathrm{s}}^2\right)}^2 $ .

During ferry impact events, however, the acceleration response amplitudes reach values on the order of $ \mathrm{m}/{\mathrm{s}}^2 $ , making the estimated measurement noise variance several orders of magnitude smaller than the structural response. As a consequence, the associated posterior uncertainty bounds are not visually distinguishable in the time-domain plots presented in previous sections.

To illustrate the sensitivity of the posterior uncertainty to the assumed measurement noise level, Figure 12 shows the estimated acceleration response at location S13 when the measurement noise variance is increased to $ {10}^{-3}\hskip0.1em {\left(\mathrm{m}/{\mathrm{s}}^2\right)}^2 $ . In this case, the uncertainty bounds become clearly visible. The shaded region corresponds to $ \pm 3\sigma $ around the posterior mean, representing a $ 99.7\% $ confidence interval. While the uncertainty bounds are now evident, the reconstruction accuracy is reduced, with the combined S&G metric decreasing from approximately $ 0.97 $ (using the noise variance inferred from resting conditions) to $ 0.95 $ .

Figure 12.

Acceleration response at location S13 with increased measurement noise covariance matrix $ \mathbf{R} $ : (a) full time history of the posterior mean estimate; (b) difference between measured response and posterior mean estimate; (c) zoomed-in time window around impact; (d) PSD. Black solid lines denote measurements, dashed red lines denote estimates, the red shaded area corresponds to the $ \pm 3\sigma $ confidence interval (99.7% confidence level), and the red vertical bar in the PSD plots indicates the natural frequency of the highest mode included in the reduced-order model.

5.4. Influence of combining acceleration with displacement data

Traditional Kalman filter-based state estimation methods are well known to be susceptible to drift when relying solely on acceleration data (Azam et al., Reference Azam, Chatzi and Papadimitriou2015). Accelerometers are insensitive to very low frequencies and quasi-static components. Incorporating displacement data can mitigate drift and improve overall estimation accuracy in such Kalman filter frameworks (Azam et al., Reference Azam, Chatzi and Papadimitriou2015). However, it has also been highlighted that GPLFMs, due to their probabilistic modeling of latent forces and inherent Bayesian inference structure, can often circumvent these drift problems when using only acceleration data (Nayek et al., Reference Nayek, Chakraborty and Narasimhan2019; Petersen et al., Reference Petersen, Øiseth and Lourens2022).

For this analysis, sensor configuration 2 was employed, as depicted in Figure 2. This setup comprises five accelerometers (A1–A5) and two LVDTs. The data were originally sampled at 1024 Hz. To facilitate a consistent comparison with analyses performed using sensor configuration 1 (which was sampled at 128 Hz) and to isolate the effect of sensor type rather than sampling rate, the data from configuration 2 were decimated to match the 128 Hz sampling frequency of configuration 1. A detailed investigation into the influence of the sampling frequency is presented later in Section 5.6.

Virtual sensing was performed using a leave-one-out approach for two target accelerometer locations: A4 and A3 (Figure 2). When estimating the response, the measured data from that specific accelerometer were excluded from the GPLFM and used for validation. Two scenarios were considered for each target: (i) using only the accelerometer data and (ii) using both the accelerometer data and the data from the two LVDTs.

Figure 13 shows the results for location A4, while Figure 14 displays the results for location A3.

Figure 13.

Acceleration response comparison at location A4 using acceleration data only and combined acceleration–displacement data: (a) full time history of the posterior mean estimates; (b) difference between measured response and posterior mean estimates; (c) zoomed-in time window around impact; (d) PSD. Black solid lines denote measurements, dashed red lines denote responses estimated using only acceleration data, dashed blue lines denote responses estimated using combined acceleration and displacement data, and the red vertical bar in the PSD plots indicates the natural frequency of the highest mode included in the reduced-order model.

Figure 14.

Acceleration response comparison at location A3 using acceleration data only and combined acceleration–displacement data: (a) full time history of the posterior mean estimates; (b) difference between measured response and posterior mean estimates; (c) zoomed-in time window around impact; (d) PSD. Black solid lines denote measurements, dashed red lines denote responses estimated using only acceleration data, dashed blue lines denote responses estimated using combined acceleration and displacement data, and the red vertical bar in the PSD plots indicates the natural frequency of the highest mode included in the reduced-order model.

The results indicate that adding LVDT data yields similar posterior mean of the acceleration response estimation compared to using acceleration data only, with S&G values also showing minimal differences between the two scenarios for each location. For location A4, the S&G value was 0.80 when using only acceleration data and 0.82 when including LVDT data. For location A3, the S&G was 0.61 for acceleration data only and 0.66 with the addition of LVDTs. This general similarity further supports the finding that the GPLFM is robust against drift effects even when only accelerometer measurements are available. For the estimation at location A4 (Figure 13), both scenarios produce high-quality reconstructions in the time domain. In the frequency domain, slight deviations from the measured PSD are observed for locations A3. The responses estimated using combined acceleration and displacement data show higher agreement with the measured PSD at higher frequencies compared to the acceleration-only case. It is important to note that the LVDTs would be much more influential if the objective of this study was the estimation of the latent forces and their static component. For completeness, the posterior mean of the estimated modal forces as well as the associated uncertainty bounds have been added to Appendix D. Consistent with the findings in Section 5.2, the accuracy of the estimation is generally lower for the location situated on the front beam where no other sensors are placed on the same primary longitudinal beam.

5.5. Sensitivity to an assumed/identified modal damping ratio

The analysis employed the GPLFM framework, which, in this implementation, relies on system matrices $ {\mathbf{A}}_c $ and $ {\mathbf{B}}_c $ derived from a time-invariant representation of the structure. This assumption implies that the underlying physical properties, such as stiffness, damping, and mass, can be assumed to be time-invariant throughout the duration of the dynamic event being modeled. Nonetheless, the identification of damping in complex real-world structures is typically very challenging, and much less accurate than the natural frequency and mode shapes, even when the LTI assumptions hold (Farrar and Worden, Reference Farrar, Worden, Deraemaeker and Worden2010). For this case study, a modal damping ratio of 5% was assumed for all modes included in the model.

To evaluate the sensitivity of the GPLFM to potential inaccuracies in this assumed damping ratio, two perturbation scenarios were investigated using sensor configuration 2 and measurement data decimated to 128 Hz: one underestimating the damping ratio at 2% and another overestimating it at 8%. The resulting estimated acceleration responses at position A4 for the baseline (5%) and perturbed damping models (2% and 8%) are compared against the measured data in Figure 15.

Figure 15.

Acceleration response comparison at location A4 with varying modal damping ratio: (a) full time history of the posterior mean estimates; (b) difference between measured response and posterior mean estimates; (c) zoomed-in time window around impact; (d) PSD. Black solid lines denote measurements, dashed red lines denote responses estimated assuming a modal damping ratio of 5%, dashed blue lines denote responses estimated assuming a modal damping ratio of 2%, dashed green lines denote responses estimated assuming a modal damping ratio of 8%, and the red vertical bar in the PSD plots indicates the natural frequency of the highest mode included in the reduced-order model.

Observation of the time domain responses in Figure 15 shows that the estimated accelerations are similar across all three assumed damping ratios (2%, 5%, and 8%). All models accurately reconstruct the overall shape of the measured response, and specifically, the timing and approximate amplitude of the main acceleration peak corresponding to the ferry impact are well captured. This suggests a high degree of robustness in the time-domain estimation to the assumed level of modal damping. Notably, even for the case with a 2% damping ratio, corresponding to a markedly underdamped system, the reconstructed response remains very similar to those obtained with higher damping assumptions. The corresponding S&G values (0.82 for 5% damping, 0.80 for 8% and 2% damping) further confirm the high level of agreement between estimated and measured signals. In the frequency domain, the PSD of the estimated responses for all three damping scenarios closely matches the measured PSD up to approximately 20 Hz. This suggests that the GPLFM can effectively compensate for damping inaccuracies.

5.6. Sensitivity to sampling frequency

The choice of sampling frequency in data acquisition directly impacts the ability to capture the full spectrum of a dynamic event. This is particularly important when dealing with impact events of short durations, which would excite high frequencies. Higher sampling rates can capture transients and higher-frequency components that might be aliased or lost at lower rates. This section investigates the performance of the GPLFM when utilizing high sampling frequency data and compares estimations made with and without LVDT data. For this analysis, sensor configuration 2 was employed, with data recorded at its original sampling frequency of 1024 Hz. The virtual sensing estimation was performed for accelerometer location A4, with its own data excluded from the GPLFM input and used for validation. The estimated and measured responses are shown in Figure 16. Although the reduced-order model retains modes only up to approximately 10 Hz, the higher sampling rate is intentionally adopted to assess how the accuracy of the GPLFM response reconstruction deteriorates at higher frequencies. While these higher-frequency components are not explicitly represented by the structural modes, the estimated response may still be informative for specific structural health monitoring purposes, such as impact characterisation or event detection.

Figure 16.

Acceleration response comparison at location A4 sampled at 1024 Hz: (a) full time history of the posterior mean estimates; (b) difference between measured response and posterior mean estimates; (c) zoomed-in time window around impact; (d) PSD. Black solid lines denote measurements, dashed red lines denote responses estimated using only acceleration data, and dashed blue lines denote responses estimated using combined acceleration and displacement data.

In the time domain, the posterior mean of the estimated responses using acceleration data only and with LVDTs provide good and similar reconstructions of the measured acceleration at A4. The S&G value was 0.77 when using only acceleration data, and slightly higher at 0.78 when LVDT data was also included. These values are slightly lower than those obtained for the same location A4 when using data decimated to 128 Hz (which were 0.80 for accelerometers only and 0.82 with LVDTs). This performance is likely attributable to the fact that most of the response energy of this ferry impact is concentrated below 50 Hz. At higher frequencies, the GPLFM reconstructs portions of the measured spectrum through the inferred latent force. This explains the increasing deviation between estimated and measured spectra beyond approximately 80 Hz.