\section{Introduction}
\label{sec:intro}

With many existing monopile-based offshore wind farms approaching the end of their service life and the frequency of extreme climate events on the rise in recent years, offshore wind turbines (OWTs) face increasingly greater risk of fatigue-induced structural failure. The fatigue life of an offshore wind turbine (OWT) is driven by its dynamic response to the rotor motion of the nacelle in addition to wind and wave loading, which fluctuates over the lifetime of the structure \cite{Maes2016}. For newly constructed offshore wind farms, where OWTs are larger and more slender in design in order to maximize power generation capacity, there is a greater risk of resonance effects from the rotor frequency coinciding with the fundamental frequency of the structure, leading to accelerated fatigue damage \cite{Devriendt2014}.
\\
\indent Structural Health Monitoring (SHM) techniques offer a means to obtain information necessary for fatigue life evaluation of existing OWTs and design validation of recently installed OWTs. In practice, the strain response at critical locations can be measured over time to predict the accumulation of fatigue damage. However, current approaches to SHM which rely on physical measurements of strains are not well adapted for OWTs, because fatigue damage in monopile-based turbines concentrates in the foundation below the mudline where accurate strain readings are impractical to obtain. The installation and maintenance of sensors at subsoil fatigue hotspots is cost-prohibitive, and for monopiles instrumented prior to installation, the sensors are often damaged during the pile-driving process \cite{Henkel2019}.
\\
\indent To overcome these challenges, \textit{virtual sensing} approaches to SHM have received considerable attention in OWT monitoring applications, since they provide a full-field estimation of an unobserved quantity of interest (e.g. strains) using measurements of a proxy response variable (e.g. accelerations) from physical sensors. These methods can be efficiently implemented using a small number of standard sensors, such as accelerometers, placed at accessible locations of an operating structure. The response at unmeasured locations of the structure is then reconstructed from the measurements through a dynamics model of the system or transfer function. This approach to SHM not only reduces sensor deployment costs significantly, but also improves the quality of data collected, as accelerometers are more reliable, less sensitive to noise, and capable of recording at higher sampling rates compared to strain gauges \cite{Noppe2016, Smyth2007}. Virtual sensing approaches may be categorized into two broad classes: 1) deterministic model-based extrapolation, where the estimation of response states beyond measured locations is obtained through a well-calibrated mechanical model, such as a finite element model, and 2) probabilistic state estimation, where strain response states are treated as random variables whose probability distributions are learned from measured data.
\\
\indent Among the deterministic techniques, the predominant method in practice is modal decomposition and expansion (MD\&E) by Iliopoulos et al. \cite{Iliopoulos2016}, which decomposes the dynamical system into mode shape components from which the response between measured points may be inferred. The accuracy of MD\&E depends on the parameters of the finite element model, which are calibrated using system identification techniques such as operational modal analysis. In applications to operational data from OWTs, MD\&E produces a good estimate of the strain response in the time domain but faces difficulty with estimating higher frequency regions of the response spectrum \cite{Maes2016}. Moreover, because MD\&E derives displacements from the double integration of measured accelerations, the strain estimate is susceptible to low-frequency noise amplification, a numerical issue termed the ``drift effect" which is detailed in \cite{Smyth2007}. To address the difficulty in the frequency domain strain estimate, Iliopoulos et al. \cite{Iliopoulos2017} introduce a multi-band method where the dynamic frequency range is separated into two sections, such that a variable number of sensors and modes may be chosen for each frequency band. Studies by Noppe et al. \cite{Noppe2016, Noppe2018} show that the multi-band method may lead to discontinuities in the estimate at the defined frequency band limits. Henkel et al. \cite{Henkel2018, Henkel2019} extends the studies in \cite{Maes2016, Noppe2016, Noppe2018} to the estimation of subsoil strains, where multi-band MD\&E is observed to result in a phase lag in the time domain and an underestimation of the first mode spectral density for periods of high wind speeds. Generally, ill performance of MD\&E can be attributed to its fully deterministic nature, such that it remains highly sensitive to discrepancies between the physical system and idealized model. In particular, MD\&E does not explicitly account for stochastic measurement noise and model error resulting from uncertainty in soil-foundation interaction properties. Moreover, since calibration of the model must be performed offline, MD\&E cannot account for shifts in the modal parameters of the system (e.g. natural frequencies, damping ratios) across operational states and loading conditions, which limits the robustness of its use. 
\\
\indent The other class of techniques takes a Bayesian inference approach to virtual sensing and has a number of notable differences from deterministic methods. First, in estimating the trajectory of the unknown state, Bayesian inference methods produce a probability distribution of values of the state over time, reflecting additional information on the level of variance in the predicted state, rather than a single estimate of the time series. Secondly, sources of uncertainty from stochastic measurement noise and model error are explicitly incorporated in the construction of the inference problem. Finally, these virtual sensing methods offer a statistical framework to account for the direct dependence of the dynamic response on qualities of the force, or ``input", acting on the system, such as its magnitude, variability, and correlation structure. As no practical means exist for direct measurement of forces, the unknown input is either implicitly modeled as a form of process noise or treated as an additional random variable whose time evolution may be derived simultaneously with that of the strain response. The latter of these approaches, in which the joint distribution of unknown inputs and states is inferred from measurement data, is referred to as ``joint input-state estimation" and constitutes the focus of this paper. In particular, this work draws a comparison of existing methods in literature regarding assumptions placed on the unknown input and proposes to implement a Bayesian inference technique with a more generic treatment of the unknown input for SHM applications. 
\\
\indent Bayesian inference approaches to virtual sensing are commonly implemented using the Kalman filter \cite{Kalman1960}, whose sequential structure to inferring states from measurements is suitable for long time series-type data. Among the Kalman filter-based techniques, one version comprises of minimum variance unbiased filters derived from optimal control theory which do not assume prior knowledge of the time-varying behavior of the input. While early methods addressed the state estimation problem exclusively \cite{Kitandis1987, Hou1998}, Hsieh et al. \cite{Hsieh2000} recognized the value of performing input estimation to improve the state estimation. The methods of \cite{Kitandis1987} and \cite{Hsieh2000} were extended in Darouach et al. \cite{Darouach2003} to include direct feedthrough between the inputs and states. These contributions lead to the development of the joint input-state estimator by Gillijns \& de Moor \cite{Gillijns2007a, Gillijns2007b}, abbreviated as GDF or JIS in literature, which formalized a filter for simultaneously estimating inputs and states in linear systems satisfying minimum variance and unbiasedness conditions. Since then, the method has been extended for cases with a time delay using Rauch-Tung-Striebel (RTS) smoothing \cite{RTS1965, Maes2018a}, for cases with a reduced order model and equivalent forces \cite{Lourens2012}, and with a method to account for estimation uncertainty and correlation between inputs and states \cite{Maes2016b}. The GDF is evaluated extensively in literature, using numerical tests in \cite{Azam2015, Tatsis2017}, experimental tests in \cite{Lourens2012, Azam2017}, and operational data in \cite{Maes2016, Lourens2012, Petersen2017}. Across the studies, it is consistently demonstrated that GDF is sensitive to the drift effect from the accumulation of integration errors when only acceleration data are used \cite{Azam2017}. The GDF is also susceptible to numerical instability due to rank deficiency of system matrices, occurring in overdetermined systems where the number of sensors exceeds the number of modes in the model \cite{Lourens2012}, and where the sensor configuration results in near linear dependence between the inputs and states \cite{Petersen2017}.
\\
\indent In contrast to the GDF, which makes no prior assumption on the input, more recent developments in Kalman filter-based estimators assume that the input takes the form of a random walk model, where the time evolution of the input is driven by zero-mean white Gaussian noise. The augmented Kalman filter, or AKF \cite{Lourens2012b},  performs simultaneous estimation of the input and state variables in augmented state-space; in comparison to Maes et al. \cite{Maes2016, Maes2016a}, which implements a classical Kalman filter for pure state estimation by assuming the process noise term encompasses contributions from the input, the AKF defines the covariance structure of the unknown input separate from that of process noise. To address issues of unobservability of the augmented state matrix when acceleration-only measurements are used, Naets et al. \cite{Naets2015} proposes the AKF-dm, supplementing the AKF with artificial "dummy displacements" as additional measurements in order to stabilize the state estimate. The dual Kalman filter, or DKF \cite{Azam2015}, has also been proposed to overcome unobservability issues by introducing a successive structure to the filter, where the input estimation is followed by the state estimation. The performance of the Kalman filter-based estimators is compared using numerical data in \cite{Tatsis2017, Nayek2019}, experimental data in \cite{Azam2017}, and operational data in \cite{Noppe2018}.
\\
\indent The primary shortcoming of AKF, AKF-dm, and DKF is that their accuracy largely depends on an appropriate choice of the covariance matrix of the unknown input, which is universally treated as uncorrelated with the dynamic states. While tuning procedures for the covariance matrix have been proposed, including heuristic rules of thumb, a maximum likelihood estimation approach \cite{Yuen2007}, and L-curve regularization \cite{Lourens2012b}, manual adjustment of the covariance matrix is generally still used to obtain a stable input estimate and evade the drift effect in the state estimate. In \cite{Azam2015, Lourens2012b}, the input covariance matrix is constrained to be diagonal with constant values in order to simplify the selection process, but places the strict assumption that all inputs have equal magnitude variance. All these methods face the additional limitation that the white noise assumption on the input may be too restrictive for real life applications. In the case of OWTs, which are assumed as linear time-periodic (LTP) systems, cyclic load sources are expected to contain harmonic components which violate the broadband assumptions of white noise input.  
\\
\indent Recently, latent force models have been proposed as a generic case of Bayesian estimators which allow for greater flexibility in the treatment of the unknown input. Latent force models mark a connection between the control theory-based perspective and machine learning-based perspective of the state estimation problem; whereas purely data-driven approaches in statistical modeling traditionally suffer in cases of extrapolation and sparse datasets, the hybrid approach of a latent force model supplements the statistical model with a mechanical model of the underlying system dynamics to guide regression in the latent variable space \cite{Alvarez2009}. In virtual sensing applications, the latent force model relates the measured accelerations, or observed variables, to the inputs and states, which comprise the unobserved or ``latent" variables of the structural system \cite{Nayek2019}. Alvarez et al. \cite{Alvarez2009} introduced the concept of the Gaussian process latent force model, or GPLFM, where a linear system is driven by an unknown Gaussian process whose covariance matrix is characterized by a kernel function. The GPLFM exploits properties of Gaussian processes to construct covariance relationships between the inputs and states, using parameters of the assumed kernel as well as correlation information supplied by the mechanical model. Posterior inference of latent variables may then be reduced to a problem of Gaussian process regression of the observed variables with joint inference of the unobserved variables. Whereas \cite{Alvarez2009} provides a convolution-based solution to Gaussian process regression, Hartikainen \& Sarkka \cite{Hartikainen2010} demonstrate that Gauss-Markhov processes may be reformulated in state-space such that regression can be more efficiently performed using Kalman filtering \cite{Kalman1960} and RTS smoothing \cite{RTS1965}. Their study is extended to the joint input-state estimation problem in \cite{Hartikainen2011}, which provides the state-space framework for the GPLFM.
\\
\indent The application of the GPLFM to structural health monitoring is presented in Nayek et al. \cite{Nayek2019}, where it is shown that the GPLFM is a generalization of the AKF which offers greater flexibility in the definition of the input in augmented state-space. Rogers et al. \cite{Rogers2020} extends the study to address system parameter estimation in addition to the joint input-state estimation problem. The study by Nayek et al. \cite{Nayek2019} serves as the first work to provide numerical validation of the GPLFM for linear structural systems and demonstration of its superior performance over GDF, AKF, AFK-dm, and DKF for joint input-state estimation. However, literature on the experimental or operational validation of the GPLFM in SHM is sparse. Petersen et al. \cite{Petersen2022} demonstrates the potential of combining analytical wind load models with the latent force model for long-span suspension bridges, but without the availability of validation data. Using data collected from an operating onshore wind turbine, Bilbao et al. \cite{Bilbao2022} evaluates the performance of the GPLFM for strain estimation across the height of the tower above the foundation. 
\\
\indent This contribution provides a critical study of the in-situ validation of the GPLFM and associated practical considerations for its use in fatigue monitoring. The fatigue assessment of OWTs serves as a suitable case study, as the GPLFM is expected to offer particular advantages in state estimation for systems subject to high model uncertainty and variable excitation sources which are not well characterized by the white noise assumption grounding other virtual sensing techniques. Moreover, this work presents the rare opportunity for validation of strain estimates in the OWT foundation below the mudline, as measurements from strain gauges embedded within the monopile are available in the considered case study. The paper is organized as follows: Section \ref{sec:theory} presents the formulation of both the mechanical model of the system and the stochastic model of the unknown input in state-space, as well as the theoretical framework of the GPLFM. Section \ref{sec:casestudy} presents the case study for in-situ validation, using subsoil data obtained in a monitoring campaign of an OWT in the Westermeerwind Park in the Netherlands. Section \ref{sec:results} presents the results of strain estimation for fatigue assessment and an evaluation of the sensitivity of the GPLFM to both measurement data fidelity and model fidelity. Section \ref{sec:conclusions} concludes and provides recommendations for future research directions.  

\subsection{Mechanical model of structural system}
\label{sec:mechmodel}

In structural dynamics, the response states of a physical system under forced vibration are described by a mechanical model which takes the form of a linear second-order differential equation:  

\begin{equation}
    \label{eq:eom}
    \mathbf{M\ddot{u}}(t)
    +\mathbf{C\dot{u}}(t)
    +\mathbf{Ku}(t)
    =\mathbf{S}_{\mathrm{p}}\mathbf{p}(t)
\end{equation}

where $\mathbf{u}(t) \in \mathbb{R}^{n_{\mathrm{u}}}$ are the displacement states at $n_{\mathrm{u}}$ degrees of freedom of a system with mass, damping, and stiffness properties defined by $\mathbf{M}, \mathbf{C}, \mathbf{K} \in \mathbb{R}^{n_{\mathrm{u}} \times n_{\mathrm{u}}}$ respectively. The system is externally excited by $n_{\mathrm{p}}$ forces $\mathbf{p}(t) \in \mathbb{R}^{n_{\mathrm{p}}}$, with spatial distribution indicated by the selection matrix $\mathbf{S}_{\mathrm{p}} \in \mathbb{R}^{n_{\mathrm{u}} \times n_{\mathrm{p}}}$ \cite{Lourens2012}. In this presentation, the system is assumed to be time-invariant, i.e. the model parameters and force locations are constant in time. 
\\ 
\indent For complex models with several degrees of freedom, the solution for the response states may be more efficiently computed using model order reduction as in \cite{Lourens2012, Maes2015}. The eigensolution to $\mathbf{K}\Phi = \mathbf{M}\Phi\bm{\Omega}^2$ returns the mass-normalized mode shapes $\phi_j$, given by the columns of matrix $\Phi \in \mathbb{R}^{n_{\mathrm{u}} \times n_{\mathrm{u}}}$, and the natural frequencies $\omega_j^2$ in rad/s, which are collected in the diagonal matrix $\bm{\Omega}^2 \in \mathbb{R}^{n_{\mathrm{u}} \times n_{\mathrm{u}}}$. In the modally reduced order formulation, the matrices are truncated into $\bar{\Phi} \in \mathbb{R}^{n_{\mathrm{u}} \times n_{\mathrm{m}}}$ and $\bar{\bm{\Omega}}^2 \in \mathbb{R}^{n_{\mathrm{m}} \times n_{\mathrm{m}}}$  to retain a subset of modes $n_{\mathrm{m}}$ which dominate the system response. By casting into modal states $\mathbf{r}(t) \in \mathbb{R}^{n_{\mathrm{m}}}$ with the transformation $\mathbf{u}(t) = \bar{\Phi} \mathbf{r}(t)$ and premultiplying by $\bar{\Phi}^\textup{T}$, Eq. \eqref{eq:eom} becomes: 

\begin{equation}
    \label{eq:eom2}
    \bar{\Phi}^\textup{T} \mathbf{M} \bar{\Phi} \mathbf{\ddot{r}}(t)
    +\bar{\Phi}^\textup{T} \mathbf{C} \bar{\Phi} \mathbf{\dot{r}}(t)
    +\bar{\Phi}^\textup{T} \mathbf{K} \bar{\Phi} \mathbf{r}(t)
    =\bar{\Phi}^\textup{T} \mathbf{S}_{\mathrm{p}}\mathbf{p}(t)
\end{equation}

With mass-normalized mode shapes, the orthogonality conditions $\bar{\Phi}^\textup{T} \mathbf{M} \bar{\Phi} = \mathbf{I}$ and $\bar{\Phi}^\textup{T} \mathbf{K} \bar{\Phi} = \bar{\bm{\Omega}}^2$ are satisfied. Assuming proportional damping, the modal damping ratios $\xi_j$ are given by $\bar{\Phi}^\textup{T} \mathbf{C} \bar{\Phi} = \bar{\bm{\Gamma}}$, where $\bar{\bm{\Gamma}} = {\rm diag}(2 \xi_j \omega_j^2)$ for $j=1,...,n_{\mathrm{m}}$. The forces are transformed into their modal contributions by $\mathbf{f}(t)=\bar{\Phi}^\textup{T} \mathbf{S}_{\mathrm{p}}\mathbf{p}(t)$, where $\mathbf{f}(t) \in \mathbb{R}^{n_{\mathrm{m}}}$ is the vector of modal force time histories. The modally reduced order model then becomes: 
\begin{equation}
    \label{eq:eom3}
    \mathbf{\ddot{r}}(t)
    +\bar{\bm{\Gamma}} \mathbf{\dot{r}}(t)
    +\bar{\bm{\Omega}}^2 \mathbf{r}(t)
    =\mathbf{f}(t)
\end{equation}

To formulate the model for filtering and inference, the continuous-time linear differential equation is converted to state-space form. By defining the state vector $ \mathbf{x}(t) = \begin{bmatrix}
\mathbf{r}(t) & 
\mathbf{\dot{r}}(t)
\end{bmatrix}^\textup{T}$, Eq. \eqref{eq:eom3} can be expressed as:

\begin{equation}
    \label{eq:dynamicsmodelfull}
    \begin{bmatrix} \mathbf{\dot{r}}(t) \\ \mathbf{\ddot{r}}(t) \end{bmatrix}
    = \underbrace{ \begin{bmatrix} \mathbf{0} & \mathbf{I} \\ -\bar{\bm{\Omega}}^2 & -\bar{\bm{\Gamma}} \end{bmatrix} }_\text{$\mathbf{A}_{\mathrm{c}}$}
    \begin{bmatrix} \mathbf{r}(t) \\ \mathbf{\dot{r}}(t) \end{bmatrix}
    + \underbrace{ \begin{bmatrix} \mathbf{0} \\ \mathbf{I} \end{bmatrix} }_\text{$\mathbf{B}_{\mathrm{c}}$}
    \mathbf{f}(t)
\end{equation}



which is succinctly written as: 
\begin{equation}
    \label{eq:dynamicsmodelshort}
    \mathbf{\dot{x}}(t) = \mathbf{A}_{\mathrm{c}}\mathbf{x}(t) + \mathbf{B}_{\mathrm{c}}\mathbf{f}(t)
\end{equation}

Eq. \eqref{eq:dynamicsmodelshort} constitutes the dynamical model of the state-space system, which describes the evolution of response states over time as a function of the system parameters, given by $\mathbf{A}_{\mathrm{c}}$, and of the input forces, whose influence is given by $\mathbf{B}_{\mathrm{c}}$. 
\\
\indent To facilitate numerical implementation, the dynamical model is converted from continuous-time to discrete-time with time step size $\Delta t$ using the zero-order hold assumption. The discretized dynamical model is as follows: 

\begin{equation}
    \label{eq:dynamicsmodeldiscr}
    \begin{gathered}
        \mathbf{x}_k = \mathbf{A}\mathbf{x}_{k-1} + \mathbf{B}\mathbf{f}_{k-1}
        \\
         \mathbf{A} = \textup{exp}(\mathbf{A}_{\mathrm{c}} \Delta t), \quad
         \mathbf{B} = [\mathbf{A} - \mathbf{I}] \mathbf{A}_{\mathrm{c}}^{-1} \mathbf{B}_{\mathrm{c}}
    \end{gathered}
\end{equation}


\indent In addition to the dynamical model, a discrete-time measurement model is defined to describe the relation between state variables and observed variables of a distinct sampling frequency. While accelerometers are predominantly used in SHM applications, axial strain readings from strain gauges may be additionally available with certain sensor configurations. Heterogeneous data fusion with displacement measurements is presented in \cite{Smyth2007, Ferrari2015}; for the case where strain measurements are available at an equal sampling rate as acceleration measurements, the measured quantities $\mathbf{y}_k \in \mathbb{R}^{n_\mathrm{y}}$ for $k=1,...,N$ are related to the original displacement, velocity, and acceleration states by: 

\begin{equation}
    \label{eq:measuremodel}
    \mathbf{y}_k = 
    \begin{bmatrix} \bm{\epsilon} \\ \mathbf{\ddot{u}} \end{bmatrix}_k
    = \begin{bmatrix} \mathbf{S}_{\mathrm{s}}\mathbf{T}_{\mathrm{s}} & 0 & 0 \\
    0 & 0 & \mathbf{S}_{\mathrm{a}} \end{bmatrix}
    \begin{bmatrix} \mathbf{u} \\ \mathbf{\dot{u}} \\ \mathbf{\ddot{u}} \end{bmatrix}_k
\end{equation}

where $\mathbf{S}_{\mathrm{a}} \in \mathbb{R}^{n_{\mathrm{acc}} \times n_{\mathrm{u}}}$ is the selection matrix of degrees of freedom where $n_{\mathrm{acc}}$ accelerations are recorded and $\mathbf{S}_{\mathrm{s}} \in \mathbb{R}^{n_{\mathrm{str}} \times n_{\mathrm{e}}}$ is the selection matrix of elements where $n_{\mathrm{str}}$ strains are recorded. $\mathbf{T}_{\mathrm{s}} \in \mathbb{R}^{n_{\mathrm{e}} \times n_{\mathrm{u}}}$ is the transformation matrix which computes axial strain deformation $\bm{\epsilon} \in \mathbb{R}^{n_{\mathrm{e}}}$ from lateral displacements of a finite element model with $n_{\mathrm{e}}$ elements based on shape functions from beam theory, provided in \ref{sec:app_disptransf}.
\\
\indent The observation equation with respect to the modal states $\mathbf{x}$ is then given by:
\begin{equation}
    \label{eq:measuremodelfull}
    \begin{bmatrix} \bm{\epsilon} \\ \mathbf{\ddot{u}} \end{bmatrix}_k
    = \underbrace{ \begin{bmatrix} \mathbf{S}_{\mathrm{s}}\mathbf{T}_{\mathrm{s}}\Phi & \mathbf{0} \\
    -\mathbf{S}_{\mathrm{a}}\Phi\bm{\Omega}^2 & -\mathbf{S}_{\mathrm{a}}\Phi\bm{\Gamma} \end{bmatrix} }_\text{$\mathbf{G}$}
    \begin{bmatrix} \mathbf{r} \\ \mathbf{\dot{r}} \end{bmatrix}_k
    + \underbrace{ \begin{bmatrix} \mathbf{0} \\ \mathbf{S}_{\mathrm{a}}\Phi \end{bmatrix} }_\text{$\mathbf{J}$}
    \mathbf{f}_k
\end{equation}

which can be reduced to the following equation, with output influence matrix $\mathbf{G}$ and direct feedthrough matrix $\mathbf{J}$: 

\begin{equation}
    \label{eq:measuremodelshort}
    \mathbf{y}_k = \mathbf{G}\mathbf{x}_k + \mathbf{J}\mathbf{f}_k
\end{equation}

Together, Equations \eqref{eq:dynamicsmodeldiscr} and \eqref{eq:measuremodelshort} give the discrete-time state-space model of the physical system: 

\begin{subequations}
    \label{eq:ssm}
    \begin{gather}
    \mathbf{x}_k = \mathbf{A}\mathbf{x}_{k-1} + \mathbf{B}\mathbf{f}_{k-1} + \mathbf{\eta}_{k-1} \\
    \mathbf{y}_k = \mathbf{G}\mathbf{x}_k + \mathbf{J}\mathbf{f}_k + \mathbf{v}_k
    \end{gather} 
\end{subequations}

The process noise term $\mathbf{\eta}_{k}$ with covariance $\mathbb{E}[\mathbf{\eta}_k\mathbf{\eta}_l^\textup{T}] = \mathbf{Q}^{\mathrm{x}}$ and measurement noise term $\mathbf{v}_k$ with covariance $\mathbb{E}[\mathbf{v}_k\mathbf{v}_l^\textup{T}] =\mathbf{R}$ are added to account for model uncertainty and measurement error, respectively. The noise terms are mutually uncorrelated white Gaussian noise processes, i.e.\ $\mathbb{E}[\mathbf{\eta}_k\mathbf{v}_l^\textup{T}] = \mathbf{0}$.
\\
\indent In its current form, the linear state-space model in Eq. \eqref{eq:ssm} cannot be used to directly solve for the response states because the forcing process $\mathbf{f}_{k}$ is unknown. Whereas Maes et al.\ \cite{Maes2016a} simplifies the problem by grouping the unknown input with the white noise term $\mathbf{w}_k$, this approach has limited robustness in capturing the frequency range of dynamic strains for events with non-white input. We are therefore motivated to develop a more flexible characterization of the force which accurately reflects in-situ conditions and improves the state estimate. 



\subsection{Gaussian process models}
\label{sec:inputmodel}

This section provides a brief overview of properties of Gaussian processes (GPs) and GP regression, which form the basis of the GPLFM methodology. In the case of a single dimensional input $t \in \mathbb{R}$, a Gaussian process defines a probability density over possible functions $h(t)$ of a random variable $t$. A GP has the property of being completely characterized by the first two moments of the process, namely its mean and covariance functions, expressed as follows: 

\begin{equation}
    \label{eq:gp}
    h(t) \sim GP(\mu(t), \kappa(t,t';\theta))
\end{equation}

 where $\mu(t)$ denotes the mean function and $\kappa(t,t';\theta)$ denotes the covariance function, also referred to as the kernel, representing the covariance as a function of distance between any two points $t$ and $t'$ of the process. Several options for forms of the kernel and its corresponding hyperparameters $\theta$ exist to control the smoothness and differentiability properties of $h(t)$, discussed further in Section \ref{subsec:GPkernels}. 
\\
\indent Gaussian processes are commonly employed in regression problems to define a non-parametric mapping between inputs and outputs, e. g. $h \, : X \rightarrow Y$. Consider a training dataset $D = \{ X, Y \}$ with inputs $X = \{t_1,...,t_N\}$ and corresponding outputs $Y = \{y_1,...,y_N\}$. The observed $Y$ are assumed to be noisy measurements of the underlying process, modeled as $h(X)$:

\begin{equation}
    \label{eq:gpreal}
    Y = h(X) + \bm{\nu}, \quad \bm{\nu} \sim \mathcal{N}(\mathbf{0},\sigma_\nu^2\mathbb{I})
\end{equation}

Discrepancy between measurements $Y$ and the model $h(X)$ is represented with additive white Gaussian noise $\bm{\nu}$ with variance $\sigma_\nu^2$. In regression, $h$ provides an approximation of the unobserved output corresponding to a test set of values $X^*$ which are separate (or held out) from the observed dataset $D$. In GP regression, moments of the Gaussian posterior distribution $p(h(X^*)|Y)$ may be computed with closed form expressions, as provided in \cite{Hartikainen2010, Rasmussen2006}. The main drawback of this solution to GP regression is the high order of computational complexity, scaling with respect to the number of data points $N$ by $\mathcal{O}(N^3)$, due to the inversion of a $N \times N$ matrix. Rasmussen \& Williams \cite{Rasmussen2006} provide an algorithm for the practical implementation of the closed form solution, by addressing matrix inversion with Cholesky decomposition. However, it can be shown that for stationary Gauss-Markhov processes, a recursive solution to GP regression by means of Kalman filtering and RTS smoothing may be implemented with a computational complexity scaling linearly with the dataset size, e. g. $\mathcal{O}(N)$. For particularly large datasets, as in cases with time series-type data, the recursive solution has substantial computational advantage. 

 
\subsubsection{Special case of stationary Gaussian processes}
\label{subsec:stationaryGPs} 
 
A stationary GP which satisfies Markhov properties may be formulated in state-space such that Kalman filtering and RTS smoothing can be implemented to solve for its posterior distribution \cite{Hartikainen2010}. For the zero-mean stationary case of the GP, the mean function evaluates to zero, e. g. $\mu(t) \rightarrow 0$, and the covariance function does not change over time shifts $\tau = t - t'$, becoming: 
 
 \begin{equation}
    \label{eq:gpstationary}
    h(t) \sim GP(0, \kappa(\tau; \theta))
\end{equation}

The temporal GP $h(t)$ of Eq. \eqref{eq:gpstationary} may be equivalently expressed as the output of a linear time-invariant (LTI) system driven by a white noise process $w(t)$, taking the form of a stochastic differential equation (SDE) of order $\beta$: 

\begin{equation}
    \label{eq:gpSDE}
    \frac{d^{\beta}h(t)}{dt^\beta} + a_{\beta-1}\frac{d^{\beta-1}h(t)}{dt^{\beta-1}} + \hdots + a_1\frac{dh(t)}{dt} + a_0h(t) = w(t) \ , \quad S_w(\omega) = q_w
\end{equation}

in terms of coefficients $a_0, ..., a_{\beta-1}$ and the spectral density of the white noise, $q_w$. Eq. (19) may be converted to state-space through its companion form, in which the mth order linear SDE is represented as a system of m first-order linear SDEs. Defining a vector of $h(t)$ and its derivatives up to order $\beta-1$ as $\mathbf{z}(t) = \begin{bmatrix} h(t) \quad \frac{dh(t)}{dt} \quad \hdots \quad \frac{d^{\beta-1}h(t)}{dt^{\beta-1}} \end{bmatrix}^\textup{T}$, the companion form follows: 

\begin{equation}
    \label{eq:singleforcemodel}
    \begin{bmatrix} \frac{dh(t)}{dt} \\ \frac{d^2h(t)}{dt^2} \\ \vdots \\ \frac{d^{\beta}h(t)}{dt^\beta} \end{bmatrix}
    = \underbrace{ \begin{bmatrix} 0  & 1 & 0 & \hdots & 0 \\
    0 & 0 & 1 & \hdots & 0 \\
    \vdots & \vdots & \vdots & \ddots & \vdots \\
    0 & 0 & 0 & \hdots & 1 \\
    -a_0    & -a_1   & -a_2 & \hdots & -a_{\beta-1} \\ \end{bmatrix} }_\text{$\mathbf{F}_{\mathrm{c}}$}
    \begin{bmatrix} h(t) \\ \frac{dh(t)}{dt} \\ \vdots \\ \frac{d^{\beta-1}h(t)}{dt^{\beta-1}} \end{bmatrix}
    + \underbrace{ \begin{bmatrix} 0 \\ 0 \\ \vdots \\ 0 \\ 1 \end{bmatrix} }_\text{$\mathbf{L}_{\mathrm{c}}$}
    w(t)
\end{equation}

or simply as:

\begin{equation}
    \label{eq:singleforcemodelshort}
    \mathbf{\dot{z}}(t) = \mathbf{F}_{\mathrm{c}}\mathbf{z}(t) + \mathbf{L}_{\mathrm{c}}w(t)
\end{equation}

The SDE in Eq. \eqref{eq:singleforcemodelshort} constitutes the state-space model of the GP. It must be highlighted that properties of the SDE - namely, the spectral density $q_w$, the order of the model $m$, and the value of coefficients $a_0, ..., a_{\beta-1}$ - are a function of the kernel chosen for the GP. Certain conditions must be met in order to construct a state-space model of the GP with the intended stationary kernel $\kappa(\tau)$: the kernel must satisfy positive semi-definiteness, and the kernel produces a GP $h(t)$ whose spectral density $S_h(\omega)$ is in rational form. For GPs which meet these criteria, the equivalent state-space model ($\mathbf{F}_{\mathrm{c}}, \mathbf{L}_{\mathrm{c}},$ and $q_w$) may be computed by spectral factorization, the steps of which are detailed in \cite{Nayek2019, Hartikainen2010}. 
\\
\indent From Eq. \eqref{eq:singleforcemodelshort}, our interest is to estimate the probability density of the state vector, $p(\mathbf{z}(t))$. The Fokker-Planck-Kolmogorov (FPK) equation provides the complete probability density of the solution of a generic SDE, as detailed in \cite{Sarkka2012}. However, due to properties of GPs,  computation of the FPK equation can be bypassed and the first two moments of the distribution can be directly solved to provide a complete probabilistic description of the state vector, $p(\mathbf{z}(t)) \sim \mathcal{N}(\hat{\mathbf{z}}(t), \mathbf{P}(t))$. For the case of linear time invariant (LTI) SDEs with the initial condition $\mathbf{z}(t_0) \sim \mathcal{N}(\hat{\mathbf{z}}_0, \mathbf{P}_0)$, the differential equations of the mean and covariance of $\mathbf{z}(t)$ reduce to:

\begin{subequations}
    \begin{gather}
    \label{eq:meanDE}
    \frac{d\hat{\mathbf{z}}(t)}{dt} = \mathbf{F}_{\mathrm{c}}\hat{\mathbf{z}}(t) \\
    \label{eq:covDE}
    \frac{d\mathbf{P}(t)}{dt} = \mathbf{F}_{\mathrm{c}}\mathbf{P}(t) + \mathbf{P}(t)\mathbf{F}_{\mathrm{c}}^\textup{T} + \mathbf{Q}_{\mathrm{c}}
    \end{gather}
\end{subequations}

where the spectral density of the stochastic term is transformed as $\mathbf{Q}_{\mathrm{c}} = \mathbf{L}_{\mathrm{c}}q_w\mathbf{L}_{\mathrm{c}}^\textup{T}$. The solutions are then given by:

\begin{subequations}
    \label{eq:meancovsoln}
    \begin{gather}
    \hat{\mathbf{z}}(t) = \textup{exp}(\mathbf{F}_{\mathrm{c}}(t-t_0)) \ \hat{\mathbf{z}}_0 \\
    \mathbf{P}(t) = \textup{exp}(\mathbf{F}_{\mathrm{c}}(t-t_0)) \ \mathbf{P}_0 \ \textup{exp}(\mathbf{F}_{\mathrm{c}}(t-t_0))^\textup{T} +
    \int_{t_0}^{t} \textup{exp}(\mathbf{F}_{\mathrm{c}}(t-\tau)) \ \mathbf{Q}_{\mathrm{c}} \ \textup{exp}(\mathbf{F}_{\mathrm{c}}(t-\tau))^\textup{T}\,d\tau   \label{subeq:covsoln}
    \end{gather}
\end{subequations}

Rather than computing the integral in Eq. \eqref{subeq:covsoln} directly, the covariance $\mathbf{P}(t)$ can be solved for from the matrix Riccati differential equation in Eq. \eqref{eq:covDE} using the method of matrix fraction decomposition, detailed in \cite{Sarkka2012, Grewal2001}. However, the solution is further simplified in the case of zero-mean stationary processes, where it may be deduced that the SDE has a \textit{steady-state} solution in which the state vector has a constant distribution at all time points, $\mathbf{z}_k \sim \mathcal{N}(\mathbf{0}, \mathbf{P}_\infty)$. Considering the asymptotic behavior $\frac{d\mathbf{P}(t)}{dt} \rightarrow \mathbf{0}$ and $\mathbf{P}(t) \rightarrow \mathbf{P}_\infty$ in the steady-state condition, Eq. \eqref{eq:covDE} becomes the Lyapunov equation, a special case of algebraic Riccati equations \cite{Sarkka2012}:

\begin{equation}
    \mathbf{0} = \mathbf{F}_{\mathrm{c}}\mathbf{P}_\infty + \mathbf{P}_\infty\mathbf{F}_{\mathrm{c}}^\textup{T} + \mathbf{Q}_{\mathrm{c}}
\end{equation}

The Lyapunov equation is solvable with numerical algorithms \cite{Bangert1986} to obtain the steady-state covariance matrix $\mathbf{P}_\infty$. To recover the kernel function $\kappa(\tau)$ of the process, we can define the following relationship to extract $h(t)$ from $\mathbf{z}(t)$: 

\begin{equation}
    \label{eq:singleforcemodelshort2}
    \begin{gathered}
    h(t) = \mathbf{H}_{c}\mathbf{z}(t) \\
    \mathbf{H}_{c} = \begin{bmatrix} 1 & 0 & \hdots & 0 \end{bmatrix}
    \end{gathered}
\end{equation}

Then, the kernel may be computed from the solution to the LTI SDE of Eq. \eqref{eq:singleforcemodelshort}, following the algebra in \cite{Sarkka2012}, to obtain: 

\begin{equation}
    \label{eq:covariancesoln}
    \kappa(\tau) = 
    \begin{cases}
    \mathbf{H}_{\mathrm{c}} \, \mathbf{P}_\infty \, \textup{exp}(\mathbf{F}_{\mathrm{c}}\tau)^\textup{T} \, \mathbf{H}_{\mathrm{c}}^\textup{T} \ , \quad \quad \ \ \
    \tau \geq 0 \\
    \mathbf{H}_{\mathrm{c}} \, \textup{exp}(-\mathbf{F}_{\mathrm{c}}\tau)^\textup{T} \, \mathbf{P}_\infty \,  \mathbf{H}_{\mathrm{c}}^\textup{T} \ , \quad \quad 
    \tau < 0
    \end{cases}
\end{equation}


\subsubsection{Recursive solution to Gaussian process regression}
\label{subsec:gpr_recursive}

For the special class of Gauss-Markhov processes, GP regression may be more efficiently performed using recursive Bayesian filtering. Unlike the batch solution to GP regression, which considers all data at once by constructing the joint distribution across points, the recursive solution processes the data sequentially by exploiting the Markhov property, where upon observation of the present state of the process $h(t)$, future instances of the state are independent of past observations \cite{Sarkka2012}. Let $y_{i:j} = \{y_i,...,y_j\}$ denote the subset of data $Y$ between time steps $t_i$ and $t_j$. First, the steady-state solution is placed as the prior, $p(\mathbf{z}_0) = \mathcal{N}(\mathbf{0}, \mathbf{P_\infty})$. Then, the independence condition between past and future observations allows for the construction of the recursive Bayesian filter as follows: 

\begin{equation}
    \label{eq:bayesrecursive}
    p(\mathbf{z}_k|y_{1:k}) = \frac{p(y_{k} |\mathbf{z}_k) \ p(\mathbf{z}_k | y_{1:k-1}))}
    {\int_{x} \ p(y_{k} |\mathbf{z}_k) \ p(\mathbf{z}_k | y_{1:k-1})) \ dx}
    \quad \forall \ k=1,...,N
\end{equation}

in which the filtering posterior distribution $p(\mathbf{z}_k|y_{1:k}) \ \forall \ k=1,...,N$ is obtained from a sequential "forward pass" through the data via Kalman filtering. What is worthy to note is that the prior distribution at time step $k$ is taken to be the posterior distribution from the previous time step $k-1$. Hence, Bayes' rule can be applied recursively as the probability distribution of $\mathbf{z}_k$ is updated incrementally with each data point in $y_{1:k}$.
\\
\indent Kalman filtering and RTS smoothing is a practical algorithm for performing recursive Bayesian estimation for linear Gaussian systems. First, the state-space equations of the GP are discretized to suit numerical implementation: 

\begin{subequations}
    \begin{gather}
        \mathbf{z}_k = \mathbf{F}\mathbf{z}_{k-1} + \mathbf{L}w_{k-1} \label{subeq:discreteGPstate1}
        \\
        h(t_k) = \mathbf{H}\mathbf{z}_k \label{subeq:discreteGPstate2}
    \end{gather}
\end{subequations}

with discrete matrices $\mathbf{F} = \textup{exp}(\mathbf{F}_{\mathrm{c}} \Delta t)$, $\mathbf{L} = \mathbf{L}_{\mathrm{c}}$, $\mathbf{H} = \mathbf{H}_{\mathrm{c}}$. In this sense, Eq. \eqref{subeq:discreteGPstate1} constitutes the dynamics model representing the time evolution of the states, and Eq. \eqref{subeq:discreteGPstate2} constitutes the measurement model relating samples of $h(t)$ to the states. As the process is a linear SDE with normally distributed states, its posterior distribution can be recovered through the Kalman filter, which provides the prediction $p(\mathbf{z}_k|y_{1:k})$ using past observations in sequence, followed by the RTS smoother, which provides the  prediction $p(\mathbf{z}_k|y_{1:N})$ using a backward pass through the full observation dataset. The combination of filtering and smoothing gives the exact solution to GP regression for a set of noisy samples of the process \cite{Hartikainen2010}. The filtering and smoothing algorithms are detailed in \eqref{sec:app_KFKS}. 
\\
\indent In order to implement Kalman filtering and RTS smoothing, the covariance matrix $\mathbf{Q}$ of the discrete white noise $\mathbf{L}w_k$ is required. $\mathbf{Q}$ may be calculated from the spectral density of the continuous case of $w(t)$, $\mathbf{Q_c}$, as:

\begin{equation}
    \mathbf{Q} = \int_{0}^{\Delta t} \textup{exp}(\mathbf{F}_{\mathrm{c}}(\Delta t - \tau)) \ \mathbf{Q}_{\mathrm{c}} \ \textup{exp}(\mathbf{F}_{\mathrm{c}}(\Delta t - \tau))^\textup{T} \ d\tau
\end{equation}

In the steady-state condition, a practical method of computing $\mathbf{Q}$ is by enforcing $\mathbf{P}_k \rightarrow \mathbf{P}_\infty$ in the Kalman filter prediction equation (see \ref{sec:app_KFKS}): 

\begin{equation}
    \mathbf{Q} = \mathbf{P}_\infty - \mathbf{F}\mathbf{P}_\infty\mathbf{F}^\textup{T}
\end{equation}













\subsection{Gaussian process latent force model}
\label{sec:gplfm}

In the GPLFM, the unknown input $\mathbf{p}(t)$ driving the dynamic response of a physical system is modeled probabilistically as a GP, and its posterior distribution over time may be inferred through GP regression. However, since no direct measurements of the process are available, regression of the unknown input is performed through its relation with observed dynamics of the system captured in a \textit{hybrid model} consisting of a mechanical component and a stochastic component. The hybrid model provides a complete stochastic representation of the system dynamics, therefore allowing for the joint posterior inference of all latent variables, including unobserved dynamics such as the strain states. 
\\
\indent Considering the reduced-order formulation of the mechanical model from Section \ref{sec:mechmodel}, let the modal states be denoted as $\mathbf{x}(t) = \begin{bmatrix} \mathbf{r}(t) & \mathbf{\dot{r}}(t) \end{bmatrix}^\textup{T}$, where $\mathbf{r}(t) = \begin{bmatrix} r_1(t) & ... & r_{n_{\mathrm{m}}}(t) \end{bmatrix}^\textup{T}$, and the modal forces be denoted as $\mathbf{f}(t) = \begin{bmatrix} f_1(t) & ... & f_{n_{\mathrm{m}}}(t) \end{bmatrix}^\textup{T}$. The modal forces may be expressed as mutually independent GPs. The companion form of the GP model from Section \ref{sec:inputmodel} for multiple modal forces may be combined in the following block-diagonal state-space model: 

\begin{subequations}
    \begin{gather}
    \label{eq:blockforcemodel1}
    \begin{bmatrix} \mathbf{\dot{z}}_1(t) \\ \vdots \\\mathbf{\dot{z}}_{n_{\mathrm{m}}}(t) \end{bmatrix} =
    \begin{bmatrix} \mathbf{F}_{\mathrm{c},1}   &    &  \\
                        &   \ddots & \\
                        &   &   \mathbf{F}_{\mathrm{c},n_{\mathrm{m}}} \\ \end{bmatrix}
     \begin{bmatrix} \mathbf{z}_1(t) \\ \vdots \\ \mathbf{z}_{n_{\mathrm{m}}}(t) \end{bmatrix} +
     \begin{bmatrix} \mathbf{L}_{\mathrm{c},1}   &    & \\
                        &   \ddots & \\
                        &   &   \mathbf{L}_{\mathrm{c},n_{\mathrm{m}}} \\ \end{bmatrix}
     \begin{bmatrix} w_1(t) \\ \vdots \\ w_{n_{\mathrm{m}}}(t) \end{bmatrix}
    \\
    \label{eq:blockforcemodel2}
    \begin{bmatrix} f_1(t) \\ \vdots \\ f_{n_{\mathrm{m}}}(t) \end{bmatrix} =
    \begin{bmatrix} \mathbf{H}_{\mathrm{c},1}   &    &  \\
                        &   \ddots & \\
                        &   &   \mathbf{H}_{\mathrm{c},n_{\mathrm{m}}} \\ \end{bmatrix}
     \begin{bmatrix} \mathbf{z}_1(t) \\ \vdots \\ \mathbf{z}_{n_{\mathrm{m}}}(t) \end{bmatrix}
     \end{gather}
\end{subequations}

written in shorthand using the block state vector $\mathbf{s}(t) = \begin{bmatrix} \mathbf{z}_1(t) & ... & \mathbf{z}_{n_{\mathrm{m}}}(t) \end{bmatrix}^\textup{T}$ as: 

\begin{subequations}
    \begin{gather}
    \label{eq:blockforcemodelshort1}
    \mathbf{\dot{s}}(t) = \mathbf{F}_{\mathrm{c}}\mathbf{s}(t) + \mathbf{L}_{\mathrm{c}}\mathbf{\tilde{w}}(t)
    \\
    \label{eq:blockforcemodelshort2}
    \mathbf{f}(t) = \mathbf{H}_{\mathrm{c}}\mathbf{s}(t)
    \end{gather}
\end{subequations}

where $\mathbf{z}_{j}(t), \, \mathbf{F}_{\mathrm{c},j}, \, \mathbf{L}_{\mathrm{c},j}, \, \mathbf{H}_{\mathrm{c},j}$, the white noise forcing $w_j(t)$ and its spectral density $q_{wj}$ are defined by the kernel and associated hyperparameters of the $j$th modal force. The spectral density $\mathbf{Q}_{\mathrm{c}}$ of the term $\mathbf{L}_{\mathrm{c}}\mathbf{\tilde{w}}(t)$ is a block matrix with $\mathbf{L}_{\mathrm{c},j}q_{wj}\mathbf{L}_{\mathrm{c},j}^\textup{T}$ along the diagonal.
\\
\indent Because the mechanical model and the GP model of the latent modal forces both take the form of a linear differential equation, it is possible to linearly combine the two into an augmented state-space model. Defining the augmented state vector $\mathbf{z}^{\mathrm{a}}(t) = \begin{bmatrix} \mathbf{x}(t) & \mathbf{s}(t)\end{bmatrix}^\textup{T}$, the GPLFM is given as: 

\begin{subequations}
\label{eq:gplfmfull}
    \begin{gather}
    \label{eq:gplfmfull1}
    \begin{bmatrix} \mathbf{\dot{x}}(t) \\ \mathbf{\dot{s}}(t)\end{bmatrix}
    = \underbrace{ \begin{bmatrix} \mathbf{A}_{\mathrm{c}} & \mathbf{B}_{\mathrm{c}} \\ \mathbf{0} & \mathbf{F}_{\mathrm{c}} \end{bmatrix} }_\text{$\mathbf{F}_{\mathrm{c}}^{\mathrm{a}}$}
    \begin{bmatrix} \mathbf{x}(t) \\ \mathbf{s}(t) \end{bmatrix}
    + \underbrace{ \begin{bmatrix} \mathbf{0} \\ \mathbf{L}_{\mathrm{c}}\mathbf{\tilde{w}}(t) \end{bmatrix} }_\text{$\mathbf{\tilde{w}}^{\mathrm{a}}(t)$}
    \\
    \label{eq:gplfmfull2}
    \mathbf{y}_k
    = \underbrace{ \begin{bmatrix} \mathbf{G} & \mathbf{J}\mathbf{H}_{\mathrm{c}} \end{bmatrix} }_\text{$\mathbf{H}^{\mathrm{a}}$}
    \begin{bmatrix} \mathbf{x} \\ \mathbf{s} \end{bmatrix}_k
    \end{gather}
\end{subequations}

or in the form: 

\begin{subequations}
\label{eq:gplfmshort}
    \begin{gather}
    \mathbf{\dot{z}}^{\mathrm{a}}(t) = \mathbf{F}_{\mathrm{c}}^{\mathrm{a}}\mathbf{z}^{\mathrm{a}}(t) + \mathbf{\tilde{w}}^{\mathrm{a}}(t)
    \\
    \mathbf{y}_k = \mathbf{H}^{\mathrm{a}}\mathbf{z}^{\mathrm{a}}_k
    \end{gather}
\end{subequations}

where the spectral density $\mathbf{Q}_{\mathrm{c}}^{\mathrm{a}}$ of the augmented white noise vector $\mathbf{\tilde{w}}^{\mathrm{a}}(t)$ is given as: 

\begin{equation}
    \label{eq:augspectraldensity}
    \mathbf{Q}_{\mathrm{c}}^{\mathrm{a}} = \begin{bmatrix} \mathbf{0} & \mathbf{0} \\
    \mathbf{0} & \mathbf{Q}_{\mathrm{c}} \end{bmatrix}
\end{equation}


Eq. \eqref{eq:gplfmshort} constitutes the continuous-discrete version of the GPLFM, describing the time evolution of both the latent forces and latent system states. The hybrid model adopts the form of a linear SDE driven by a steady-state Gaussian input, and because the solution to the SDE is a linear operation of the input, the resulting system states are also GPs \cite{Sarkka2019}. Therefore, joint posterior inference of the augmented states $\mathbf{z}^{\mathrm{a}}(t)$ reduces to a GP regression problem, where regression of the measured variables provides evidence to simultaneously extrapolate the latent variable response \cite{Hartikainen2011}.
\\
\indent Kalman filtering and RTS smoothing may be implemented to solve GP regression. The model is discretized and noise sources are considered for the purpose of filtering and inference: 

\begin{subequations}
    \begin{gather}
        \label{eq:gplfmshortdiscr1}
        \mathbf{z}^{\mathrm{a}}_k = \mathbf{F}^{\mathrm{a}}\mathbf{z}^{\mathrm{a}}_{k-1} + \mathbf{w}^{\mathrm{a}}_{k-1} 
        \\
        \label{eq:gplfmshortdiscr2}
        \mathbf{y}_k = \mathbf{H}^{\mathrm{a}}\mathbf{z}^{\mathrm{a}}_k + \mathbf{v}_k
    \end{gather}
\end{subequations}

where as in Eq. \eqref{eq:ssm}, $\mathbf{v}_k$ is additive white Gaussian noise with covariance $\mathbf{R}$. Process noise is considered within the term $\mathbf{w}^{\mathrm{a}}_k$, a notion which is discussed further in Section \ref{subsec:noise}. 
\\
\indent We may adopt the same steps as in Section \ref{subsec:gpr_recursive} to derive moments of the LTI SDE solution. In particular, the augmented steady-state covariance $\mathbf{P}_\infty^{\mathrm{a}}$ of the augmented states $\mathbf{z}^{\mathrm{a}}(t)$ is computed from the Lyapunov equation: 

\begin{equation}
    \mathbf{0} = \mathbf{F}^{\mathrm{a}}_{\mathrm{c}}\mathbf{P}^{\mathrm{a}}_\infty + \mathbf{P}^{\mathrm{a}}_\infty(\mathbf{F}^{\mathrm{a}}_{\mathrm{c}})^\textup{T} + \mathbf{Q}^{\mathrm{a}}_{\mathrm{c}}
\end{equation}

and the augmented covariance matrix $\mathbf{Q}^{\mathrm{a}}$ of the discrete noise $\mathbf{w}^{\mathrm{a}}_k$ follows as: 

\begin{equation}
    \label{eq:augcovmatrix}
    \mathbf{Q}^{\mathrm{a}} = \mathbf{P}^{\mathrm{a}}_\infty - \mathbf{F}^{\mathrm{a}}\mathbf{P}^{\mathrm{a}}_\infty(\mathbf{F}^{\mathrm{a}})^\textup{T}
\end{equation}


\subsubsection{Parameter inference for the GP kernel function}
\label{subsec:GPkernels}

The performance of the GPLFM in joint input-state estimation depends on an appropriate choice of the GP kernel and its hyperparameters. Generally, the kernel is chosen from a standard set of convenient kernels which meet the criteria for conversion to state-space, as discussed in Section \ref{subsec:stationaryGPs}, based on qualitative prior knowledge of the time-varying behavior (e.g. periodicity, smoothness and differentiability) of the stochastic input. One example of a suitable kernel is the Mat\'ern kernel function, which characterizes processes that are finitely differentiable \cite{Hartikainen2010}. In this work, we exclusively make use of the Mat\'ern kernel function, as it is flexible enough to accommodate a broad range of process traits and is commonly employed in engineering applications to represent physically realistic processes. In the case of the Mat\'ern kernel with smoothing parameter $p=2$, covariance relationships are governed by two hyperparameters $\theta = [\alpha, l_s]$, where the vertical length scale $\alpha$ controls the amplitude and the horizontal length scale $l_s$ controls the rate of time variability, and therefore the frequency content, of the resulting function.  For convenience, the state-space model of the Mat\'ern kernel is provided in \ref{sec:app_matern}. Other common kernel choices which are valid for state-space implementation include the white noise (or Brownian motion) kernel function, which represents each time step of the GP as independent and identically distributed; the squared exponential kernel function, which describes infinitely differentiable processes; and the periodic kernel function, which captures strong harmonic components in time series data. Both the squared exponential kernel and periodic kernel require an approximation of the spectral density to be converted into state-space form, the details of which are provided in \cite{Hartikainen2010} and \cite{Solin2014}, respectively. For certain applications, GP regression may be particularly sensitive to the choice of kernel, where the substitution of another kernel function leads to a substantially different posterior prediction. While an evaluation of GP kernel sensitivity is beyond the scope of this work, an optimization-based approach to such a sensitivity study can be found in \cite{Stephenson2021}.
\\
\indent Given a GP kernel, the corresponding hyperparameters $\theta$ which govern the length scales of the resulting process may be tuned according to the measurement dataset for a particular system. Nayek et al. \cite{Nayek2019} proposes to use maximum likelihood estimation for parameter inference, in which the optimal hyperparameters are obtained from the least squares solution which minimizes the residual norm between the observed data and model predictions. For Markhov processes, maximization of the likelihood function $p(Y|\theta)$ may be reframed as the minimization of the negative log likelihood and computed iteratively using the Kalman filter equations. However, since such linear least squares problems are frequently ill-posed \cite{Kaipio2005}, the negative log likelihood function may not result in real-valued evaluations and its optimization is known to be susceptible to false solutions at local minima \cite{Nayek2019}. A discussion of the computational issues can be found in \cite{Petersen2022}. Rogers et al. \cite{Rogers2020} suggest a MCMC Metropolis random walk algorithm to enforce a log acceptance ratio and avoid complex values in the objective function. Moreover, they note that initialization of the objective function at the prior $p(\theta)$ produces the maximum a posteriori (MAP) estimate of the hyperparameters.  
\\
\indent As the iterative computation of the objective function using the Kalman filter is computationally intensive, particularly for long time-series data, this work implements an a priori fitting procedure for obtaining optimal values of the hyperparameters $\theta^*$. Through properties of linear Gaussian systems, the GPLFM inherently places a Gaussian prior over the both latent and observed variables of the system. As prior assumptions are encoded by the kernel function assumed in the GPLFM, the hyperparameters of the Gaussian prior correspond exactly to those of the augmented model. Empirical qualities of the measurement data are then leveraged as prior information on the dynamic behavior of the latent inputs and states, used to tune the kernel hyperparameters. In this approach, model parameters are fit by minimizing the distance $d$ between the empirical distribution of the measurements, $p(Y) \sim \mathcal{N}(\mathbf{0}, \mathbf{P}^{\mathrm{y}}_0)$, and the modeled Gaussian prior on the observed states $p(Y|\theta) \sim \mathcal{N}(\mathbf{0}, \mathbf{\hat{P}}^{\mathrm{y}}_{\infty})$, where by the steady-state condition, the prior covariance is $\mathbf{\hat{P}}^{\mathrm{y}}_{0} = \mathbf{\hat{P}}^{\mathrm{y}}_{\infty}$. For shorthand notation, $p(Y|\theta)$ is referred to by $p^{\theta}(Y)$. The optimal $\theta^*$ is then solved from: 

\begin{equation}
    \label{eq:mle}
    \theta^* = \underset{\theta}{\arg\min} \ d(p, p^\theta)
\end{equation}

In this work, the distance metric $d$ between two probability measures is taken to be the Hellinger distance $H(p, p^\theta)$, which is provided in \ref{sec:app_hellinger}.
\\
\indent As the distributions are taken to be zero-mean, parameter inference reduces to a problem of matching the covariance predicted by the model to the underlying covariance in the data. The empirical measurement covariance $\mathbf{P}^{\mathrm{y}}$ is computed directly from time-series data of the sensor channels. The modeled measurement covariance, $\mathbf{\hat{P}}^{\mathrm{y}}_{\infty}$, may be derived from the augmented steady-state covariance encoded by the GP kernel function, $\mathbf{P}^{\mathrm{a}}_\infty$, by passing it through the measurement model $\mathbf{H}^{\mathrm{a}}$ of the GPLFM with the transformation: 

\begin{equation}
    \label{eq:measurementcov}
    \mathbf{\hat{P}}^{\mathrm{y}}_{\infty} = \mathbf{H}^{\mathrm{a}} \, \mathbf{P}^{\mathrm{a}}_\infty \, (\mathbf{H}^{\mathrm{a}})^\textup{T}
\end{equation}

\indent Rather than computing the likelihood function iteratively over each measurement, this approach takes advantage of the linear Gaussian nature of the GPLFM, such that the likelihood of $\theta$ may be computed efficiently using low-order moments of long time-series data. As the distance between Gaussian distributions in Eq. \eqref{eq:mle} results in a smooth, unimodal optimization surface across a defined search domain, the minimization problem may be solved efficiently using gradient descent methods.


\subsubsection{Treatment of noise in the GPLFM}
\label{subsec:noise}

A feature of the GPLFM which stands in contrast to other Kalman filter-based estimators is that the definition of process noise on the augmented states through the stochastic term $\mathbf{w}^{\mathrm{a}}_k$ is inherent to the formulation of the GP in state-space. Therefore, the noise is completely defined by the kernel and hyperparameters adopted in the augmented model. Using Eq. \eqref{eq:augcovmatrix}, an appropriate augmented process noise covariance matrix $\mathbf{Q}^{\mathrm{a}}$ may be solved for by a data-driven technique of parameter inference rather than by heuristic tuning or ad-hoc methods. Moreover, the resultant covariance matrix captures correlations between the unknown inputs and states which are generally neglected in the AKF and DKF.  
\\
\indent The underlying assumption in the GPLFM of Eq. \eqref{eq:gplfmfull} is that the mechanical model is a perfect representation of the system, and sources of noise - any dynamical behavior which is unaccounted for in the mechanical model - are treated as additional external input. In this sense, the GPLFM does not distinguish between noise contributions and input excitation, which are characterized together as a GP perturbed by additive white noise. As a consequence, the input estimation is relatively sensitive to large model uncertainty, and the strain estimation may be susceptible to overconfidence. This aspect is generally not a concern in applications where one is less interested in the input estimation for which validation is impractical. However, for applications where an accurate input estimate is necessary, or wider variance bounds on the strain estimate are expected, process noise in the system states $\mathbf{x}$ may be distinguished either by the addition of process noise in the latent system states $[\mathbf{Q}^{\mathrm{x}} \quad \mathbf{0}; \quad  \mathbf{0} \quad \mathbf{0}]$, adding $\mathbf{Q}^{\mathrm{x}} \in \mathbb{R}^{2n_{\mathrm{m}} \times 2n_{\mathrm{m}}}$ to the augmented covariance matrix in Eq. \eqref{eq:augcovmatrix}, as in \cite{Nayek2019, Petersen2022}, or by the incorporation of the spectral density of the process noise $\mathbf{\eta}(t)$ into the vector of hyperparameters $\theta$ to solve for, as in \cite{Rogers2020}. 
\\
\indent In accounting for measurement noise, rather than solving for a prior estimate of the measurement noise covariance, a method for posterior tuning of $\mathbf{R}$ using the measurement data is adopted from  \cite{Bilbao2021}. Upon obtaining a posterior estimate $p(\mathbf{z}^{\mathrm{a}}(t) | Y, \theta) = \mathcal{N}(\hat{\mathbf{z}}^{\mathrm{a}}(t), \mathbf{P}^{\mathrm{a}}(t))$ from Kalman filtering/RTS smoothing, the measurement noise is taken to be the residual $\mathbf{r}_k$ between the observed data $Y$ and prediction on the observed variables $\mathbf{H}^{\mathrm{a}}\hat{\mathbf{z}}^{\mathrm{a}}_k$ as follows: 

\begin{subequations}
    \label{eq:residual}
    \begin{gather}
        \mathbf{y}_k = \mathbf{H}^{\mathrm{a}}\hat{\mathbf{z}}^{\mathrm{a}}_k + \mathbf{r}_k
        \\
        \mathbf{r}_k = \mathbf{y}_k - \mathbf{H}^{\mathrm{a}}\hat{\mathbf{z}}^{\mathrm{a}}_k 
    \end{gather}
\end{subequations}

Measurement noise tuning then becomes a problem of minimizing the error between the assumed measurement noise covariance $\mathbf{R}$ in the model and the covariance of the residual time series $\mathbf{\hat{R}} = \mathbb{E}[\mathbf{r}_k\mathbf{r}_k^\textup{T}]$, executed through iterations of the posterior estimation. Further details are provided in \cite{Bilbao2021}. An overview of the GPLFM methodology is illustrated in Figure \ref{fig:flowchart}. 


\begin{figure}[h!]
    \centering
    \includegraphics[width=0.95\linewidth]{figures/gplfm_flowchart.png}
    \caption{Diagram of GPLFM methodology.}
    \label{fig:flowchart}
\end{figure}



\subsubsection{Practical notes for numerical implementation}
\label{subsec:practical}

This section briefly discusses numerical sensitivities particular to the presented formulation of the GPLFM and approaches to resolve such issues, which are in part sourced from \cite{Grewal2001}.
\\
\indent For physical systems with several degrees of freedom, operations on large matrices in the Kalman filter and RTS smoother is susceptible to roundoff error and challenges with the inversion of sparse matrices. With model order reduction of the mechanical model in Section \ref{sec:mechmodel}, the dimension of the system becomes equal to the $n_\mathrm{m}$ number of modes considered, thus reducing the cost of computation and controlling for stable matrix inversion. Moreover, the reduced order formulation produces an estimate of the modal contribution of any number of forces rather than the exact forces, such that knowledge of the location of the unknown inputs is not necessary \cite{Lourens2017, Lourens2019}. For instance, in the transformation of forces into their modal contributions by $\mathbf{f}(t)= \bar{\Phi}^\textup{T} \mathbf{S}_{\mathrm{p}}\mathbf{p}(t)$, $\mathbf{S}_{\mathrm{p}}$ may be taken to be the identity matrix such that the collective effect of all possible forces $\mathbf{p}(t)$ acting at any degree of freedom is considered and no error is introduced from incorrect assumptions on the force locations.
\\
\indent The Kalman filter algorithm is sensitive to ill-conditioning of the covariance matrices, which more commonly occurs in the companion state-space model where higher derivatives of a process are included in the state vector \cite{Grewal2001}. The ill-conditioning may be resolved through a) reformulating the update equation in the Kalman filter using matrix decomposition methods, such that filtering is performed using stable factors of the covariance matrix $\mathbf{P}_k$, and b) ensuring the addition of non-zero measurement noise $\mathbf{R}$ to stabilize the update equation in events of near-zero change from $\mathbf{P}_{k-1}$ to $\mathbf{P}_k$. 
\\
\indent Numerical instability may arise if there exists a large scale difference between components of the state vector $\mathbf{x}$ or measurement vector $\mathbf{y}$, as the filtering of large components may cause cancellation of the significant digits of small components. Scale differences are particularly common when heterogeneous measurement data of different units are used, as is the case with the combination of acceleration and strain measurements. Such numerical errors can be avoided through either a change of units or a rescaling of the vector. In this work, we recommend at a minimum to normalize the measured quantities $\mathbf{y}$ prior to filtering, then applying the reverse operation to the posterior estimate to recover the true signal magnitude. 

\section{Mathematical formulation}
\label{sec:theory}



\input{02.1_MechanisticModel}
\input{02.2_StochasticModel}
\input{02.3_GPLFM}


\section{Offshore wind turbine case study}
\label{sec:casestudy}

In this contribution, validation of the GPLFM for dynamic strain estimation is performed using measurement data from the SWT 3.0-DD offshore wind turbine on a monopile foundation, referred to as the W27 turbine hereon. The W27 turbine is a member of Westermeerwind Park, an operating near-shore wind farm situated in the IJsselmeer Lake in the Netherlands. The instrumentation scheme for the W27 turbine is rare among industry OWTs, as it contains strain sensors on the monopile foundation at varying depths below the mudline. 


\subsection{Measurement setup and data processing}
\label{sec:measurement}

The measurement set used for validation consists of the following data, which was gathered from separate data acquisition systems. All positions of sensors are given with respect to the mean water level (NAP) and illustrated in Figure \ref{fig:sensorsetup}.
\begin{itemize}
    \item Rotor-nacelle assembly (RNA) data: data obtained from the SCADA system, where associated sensors are located at the level of the rotor axle at +95.0m NAP. Accelerations in the fore-aft (FA) and side-side (SS) directions of the turbine are recorded at the RNA, where FA references the direction of the rotor axis and SS is the perpendicular direction. In addition, the RNA data consists of measurement channels for wind direction and wind speed, rotor speed, blade pitch angle, power generation level, and yaw angle (rotational position of the nacelle).
    \item Support structure (SUS) data: data from accelerometers and strain gauges installed along the height of the turbine tower and monopile foundation.  Three accelerometers are located above the waterline: Acc. 1, measuring in the east-west direction at +4.5m NAP; Acc. 2, measuring in the north-south direction at +4.5m NAP; and Acc. 3, measuring along the 39.4\textdegree - 219.4\textdegree \ line from North at +42.5m NAP. The SUS data includes nine rings of 4 strain gauges labeled A, B, C, D, measuring axial strains at four faces of the monopile; however, several of the strain gauges at positions C and D, as well as all strain gauges in the first ring level, were damaged during the pile-driving process. Therefore, strain gauge readings at positions A and B in the five rings located below the mudline (Rings 2-6) and three rings located above the mudline (Rings 7-9) are retained for validation. The inclinometer reading at +4.5m NAP is neglected in this study. 

\end{itemize}

\indent The RNA and SUS data are collected continuously and organized into 10-minute increments over a 15-month period between 2016 and 2017. The turbine is monitored over various operational states of the OWT, namely the idling state (minimal rotor motion), production state (full rotor speed), and run-up state (the transitional state). Each of the 10-minute records is categorized by operational state using state curves, such as the wind speed-rotor speed curve in \cite{Sun2019}.
\\
\indent Due to the use of different data acquisition systems, signal processing steps are applied to ensure synchronization between the RNA and SUS data. First, the RNA data is downsampled to meet the 50 Hz sampling rate of the SUS data, which is determined to be sufficient for capturing the frequency resolution of all data channels, as the significant frequency content is expected to fall below 10 Hz. Resampling is followed by the time series synchronization of the data channels based on cross-correlation sequences. Full details on the time series synchronization procedure are provided in \cite{Laanes2021}. Because the primary interest in this study is to estimate dynamic strain response which drives fatigue damage accumulation, quasi-static components are removed via a bandpass Butterworth filter of order 6 with stopband frequencies of 0.12 Hz and 10 Hz. 

\begin{figure}[hbt!]
    \centering
    \includegraphics[width=0.75\linewidth]{figures/measurementsetup.jpg}
    \caption{Position of 8 strain gauge rings and 3 accelerometers along the support structure of the W27 wind turbine. Figure adopted from [Versteijlen et al.].}
    \label{fig:sensorsetup}
\end{figure}



\subsection{Modal analysis}
\label{sec:oma}

To identify in-situ modal parameters of the structural system - such as mode shapes, natural frequencies, and damping ratios - operational modal analysis (OMA) is conducted using the covariance-driven stochastic subspace identification (SSI-Cov) technique \cite{Peeters1999}. For the W27 turbine, acceleration measurements at three elevations from Acc. 1 and 2, Acc. 3, and the accelerometer in the RNA are transformed into the local RNA coordinate system to allow for separate identification of the FA and SS modes. From the stabilization diagram produced by the SSI-Cov routine, stable columns corresponding to structural modes are identified using the OPTICS clustering algorithm \cite{Ankerst1999, Boroschek2019}. 
\\
\indent To ensure accurate and consistent results, a number of considerations are made for the selection of 10-minute records to use for modal analysis. A) The dataset is filtered for records which correspond to the idling state, as loading conditions in this state most closely match the assumption of white noise excitation in the SSI-Cov method and allow for structural resonance to be discerned from noise \cite{Devriendt2014}. B) The records are further filtered based on the yaw angle of the turbine in order to allow for distinction between the FA and SS modes. The specific sensor configuration of the W27 turbine presents the limitation that Acc. 3 only provides a uniaxial measurement along the 39.4\textdegree - 219.4\textdegree \ line, meaning it cannot measure components of motion perpendicular to its axis. Therefore, the RNA must be oriented parallel to the axis of Acc. 3 for the channel to be utilized in identifying FA modes and perpendicular to its axis for SS modes, such that the measurement corresponds completely to the motion in the respective direction. Records with little to no deviation in yaw angle across the 10-minute period are chosen to reduce conflation with out-of-plane modes. C) Of the remaining records, those which correspond to low wind speeds are prioritized in the selection in order to reduce the effect of aerodynamic damping on the modal estimate. D) Records in the final selection are truncated to sections of the data with relatively consistent time-frequency content which reflects stationary conditions. Following these criteria, the FA and SS modal parameters are each derived from twelve 10-minute records taken across the 15-month duration of the monitoring campaign, with a sample record shown in Figure \ref{fig:datasample}. It should be clarified that these criteria for record selection only applies to the dataset used for modal identification by SSI-Cov; the validation records used in strain estimation are not subject to such constraints, but rather chosen broadly across operational states and wind speeds. 



\begin{figure}[hbt!]
    \centering
    \includegraphics[width=\linewidth]{figures/rawdata_2_9.png}
    \caption{Sample 10-minute record of the three acceleration channels used for modal identification, shown with the filtered time series (left column), discrete Fourier transform (center column), and spectrogram (right column) of the data. The cutoff frequencies of the bandpass filter are indicated in the Fourier spectrum.}
    \label{fig:datasample}
\end{figure}

\subsection{Finite element model}
\label{sec:fem}

A mechanical state-space model of the W27 turbine is derived from a Finite Element (FE) model which is first constructed based on known geometric and material properties, then tuned to match with the modal parameters identified through OMA. A relatively simple FE model, such as a cantilever beam model, is sufficient for capturing the dominant modal properties of the system. The FE model consists of 124 Timoshenko beam elements, each 1.0 m in length, to span across the 29.0 m monopile depth and 95.0 m tower height. The tower does not have a separate transition piece as the monopile and tower have a bolted connection. The tapered profile of the top section of the turbine is modeled using pipe cross sectional elements which diminish in diameter linearly with height. All elements of the turbine are made of structural steel with Young's modulus $E$ = 210 GPa, Poisson's ratio $\nu$ = 0.3, and mass density $\rho$ = 7850 kg/m\textsuperscript{3}. External components such as the RNA, platform, and power unit are represented with concentrated masses at their respective elevations. 
\\
\indent A nontrivial source of model uncertainty is in the subsoil properties of the foundation. Mass and stiffness contributions of the soil are modeled by adding a mass density of 1500 kg/m\textsuperscript{3} and a stiffness profile to the elements located below the mudline. A custom 1D soil stiffness profile for the W27 turbine is taken from Versteijlen et al. \cite{Versteijlen2017a, Versteijlen2017}, where detailed studies of soil characterization and modeling were conducted for the site. In their work, geophysical measurements from Cone Penetration Tests (CPT), Seismic Cone Penetration tests (SCPT), and laboratory testing of borehole samples are used to develop an initial 1D effective Winkler profile $k_{eff}(z)$ which fits elastic properties of the soil-structure system as a function of depth $z$. Then, in-situ shaker tests of the W27 monopile are performed to obtain dynamic strain measurements from Rings 2-7, which are used to solve for an optimal scale factor $\gamma$ which tunes the magnitude of the profile. It is found that while the optimal $\gamma$ varies as a function of frequency, a constant factor $\gamma=0.78$ is most appropriate for the low-frequency regime which governs the response in OWTs. The resultant $\gamma \, k_{eff}(z)$ is illustrated as a distributed springs model in Figure \ref{fig:foundation}. Note that in the absence of high-fidelity soil characterization data, the p-y curve \cite{Versteijlen2017a} is an alternative soil stiffness profile which is commonly adopted in practice.
\\
\indent Table \ref{tbl:modalparams} compares five modes of the final FE model, denoted with subscript $fem$, with those identified by OMA, denoted with subscript $id$. There is a very good match in the natural frequencies for the first four modes, with the relative error well below 1\%. The high modal assurance criterion (MAC) above 0.9 indicates a strong fit between the mode shapes, as is visible in Figure \ref{fig:modalparams}. In general, the FA and SS modes are close to symmetric; any deviations from total symmetry are likely due to differential local soil damping and stiffness. There is greater difficulty in fitting the third SS mode of the FE model to the OMA results, which is due to the fact that the dynamic response of the OWT is controlled primarily by the first two modes. Therefore, higher modes are identified with significantly greater sample variance and reduced certainty in the OMA procedure; for instance, it was not possible to identify the third FA mode from the dataset. It is expected that the first two modes in each of the FA and SS directions dominate modal contributions to the dynamic response of the OWT across operating states. 

\begin{figure}[hbt!]
    \centering
    \includegraphics[width=0.8\linewidth]{figures/foundationprofile.jpg}
    \caption{Diagram of Winkler foundation model and soil stiffness profile. Figure adopted from [Versteijlen et al.].}
    \label{fig:foundation}
\end{figure}

\begin{figure}[hbt!]
    \centering
    \includegraphics[width=0.8\linewidth]{figures/modeshape.png}
    \caption{Mode shapes identified by OMA (solid blue) and estimated by the finite element model (dashed red).}
    \label{fig:modalparams}
\end{figure}

\begin{table}[!htp]
\centering
\fontsize{8}{8}\selectfont
\setlength{\tabcolsep}{20pt}
\begin{tabular}{llllll}
\toprule
Mode & $\xi_{\mathrm{id}}$ \ ${[}\%{]}$ & $\textit{f}_{\mathrm{id}}$ \ $\textrm{{[}Hz{]}}$ & $\textit{f}_{\mathrm{fem}}$ \ $\textrm{{[}Hz{]}}$  &  $\textit{e}_{\textit{f}}$ \ ${[}\%{]}$ & MAC   \\ \midrule
FA1  & 1.932               & 0.296                  & 0.296                       & 0.01\%                & 0.999 \\
SS1  & 3.011               & 0.290                  & 0.289                       & -0.12\%               & 0.996 \\
FA2  & 2.379               & 1.869                  & 1.874                       & 0.23\%                & 0.999 \\
SS2  & 2.887               & 1.885                  & 1.883                       & -0.09\%                & 0.939 \\
SS3  & 2.604               & 4.563                  & 5.193                       & 13.79\%               & 0.743 \\ \bottomrule
\end{tabular}
\caption{Modal parameters (a) identified by OMA and (b) estimated by the finite element model.}
\label{tbl:modalparams}
\end{table}



\subsection{Hyperparameter fitting}
\label{sec:prior}

For each record in the validation dataset, the optimal hyperparameters of the GPLFM are fit according to the minimum variance error approach of Equation \eqref{eq:mle}. Evaluation of the objective function produces a convex surface with a clear optimum $\theta^*$, illustrated in Figure \ref{fig:HPsurface}a for sample record P3 with search bounds $\alpha=[0, 50]$ and $l_s=[0, 0.8]$. To visualize the goodness of fit, the empirical distribution of the measurements (in blue) and the predicted prior distribution corresponding to the optimal model hyperparameters $\theta^*$ (in orange) are overlaid on the time series plots and histograms of the measured variables in Figure \ref{fig:HPsurface}b. The near exact match in the distributions across all three measurement channels indicates an appropriate fit of the model to prior information on the true system dynamics.
\\
\indent The optimal hyperparameters for all validation records are shown in Figure \ref{fig:HPtrends}. In the leftmost plot, the hyperparameter sets $(\alpha, l_s)$ for each of the operating states are observed to fall in distinct clusters. In particular, the idling state (records D1, ..., D8) corresponds to large values of $l_s$ and small values of $\alpha$ due to the dominance of low frequency oscillations driven by low-powered ambient excitation of the tower. The production state (records P1, ..., P8) corresponds to small values of $l_s$ and large values of $\alpha$ due to the combination of high wind speeds and higher-frequency motion resulting from interaction with the turbine rotor. From these correlations, we may draw physical interpretation of the controlling hyperparameters of the Mat\'ern kernel as they relate to qualities of the unknown input. To substantiate this reasoning, we refer to the center plot showing a positive correlation between $\alpha$ and the mean wind speed of the record, implying the vertical length scale parameter relates to the magnitude of the input as it varies with wind speed conditions. In the rightmost plot, $l_s$ is compared with a metric of relative power in the low frequency range, which is computed as the ratio of the power (integration of area under the power spectral density function) from 0.1 to 1.0 Hz with the power from 1.0 to 5.0 Hz. Larger values of this ratio indicate greater contributions of low frequencies to the response spectrum. From the positive correlation observed, it can be inferred that the horizontal length scale parameter relates to the dominance of low frequency inputs, which varies with both ambient conditions and rotor-induced loads.
\\
\indent The run-up state (records R1, ..., R8) corresponds to hyperparameters falling between the ranges of the idling and production states, as it marks the transition between these two states. Moreover, lower variance in the hyperparameters is observed compared to the other two states because run-up of the OWT always occurs in predefined wind conditions, whereas the OWT may remain stationed in the idling or production state across a wider range of wind speeds (see Figure \ref{fig:statecurve}). 


\begin{figure}[hbt!]
    \centering
    \includegraphics[width=0.9\linewidth]{figures/priorfit2.jpg}
    \caption{a) The optimization surface across the search domain $\alpha=[0, 50]$ and $l_s=[0, 0.8]$ for hyperparameters of the Mat\'ern $\nu=5/2$ kernel. b) Measurement time series and histogram of acceleration data at 3 channels for record P3, with the empirical distribution indicated in blue and modeled distribution indicated in orange.}
    \label{fig:HPsurface}
\end{figure}

\begin{figure}[hbt!]
    \centering
    \includegraphics[width=0.95\linewidth]{figures/hp_plot.png}
    \caption{Optimal $l_s$ vs. optimal $\alpha$ (left), wind speed vs. optimal $\alpha$ (center), and optimal $l_s$ vs. relative power in low-frequency range 0.1-1.0 Hz (right) for eight records in each operational state.}
    \label{fig:HPtrends}
\end{figure}
\subsection{Response estimation}
\label{sec:responseestimation}

The strain estimation results in the idling, run-up, and production states are given in Figures \ref{fig:strain_idling}, \ref{fig:strain_runup}, and \ref{fig:strain_production}, respectively. Overall, the GPLFM attains remarkable accuracy in the estimation of the frequency content of sub-soil strains; the primary discrepancy observed is in the high frequency range of the strain power spectral density (PSD) in Figure \ref{fig:strain_idling}, where low signal-to-noise ratio in the idling state results in difficulty with capturing the signal where it falls below the level of sensor noise. This discrepancy disappears for the run-up and production states, as the response is no longer obscured by the noise amplitude. Despite the good performance in capturing the dominant frequency components, an overestimation of the strain response is evident from both the time- and frequency domain plots. The error in the amplitude is most apparent in the idling state in Figure \ref{fig:strain_idling} as well as at periodic impulses caused by turbine rotor motion in the run-up and production states, shown in the time series of Figures \ref{fig:strain_runup} and \ref{fig:strain_production} respectively. 
\\
\indent Strain estimation results for all validation records is shown in Figure \ref{fig:accuracystats}, where the boxplots are computed from accuracy metrics of the strain estimate across eight records in each operating state. In general, high accuracy in the frequency domain estimate is reflected by a very high correlation coefficient (CC) near 1, while fair accuracy in the time-domain estimate is reflected by a mean relative error (MRE) generally falling in the range of 5-10\%. There is a reduction in accuracy in the run-up state due to larger variance in the wind and rotor speed (and therefore aerodynamic damping) among the records. Lower accuracy observed at strain gauge 2 may be a result of extrapolation uncertainty, as the sensor is located at the furthest depth in the foundation where uncertainty in the soilbed properties is highest. 



\begin{figure}[hbt!]
    \centering
    \includegraphics[width=\linewidth]{figures/strain_idling_4_4.png}
    \caption{Strain estimation for a record from the idling state (D7), shown with the full time series (left column), zoomed-in view of the time series (center column), and power spectral density (right column) for five strain channels.}
    \label{fig:strain_idling}
\end{figure}


\begin{figure}[hbt!]
    \centering
    \includegraphics[width=\linewidth]{figures/strain_runup_10_2.png}
    \caption{Strain estimation for a record from the run-up state (R5), shown with the full time series (left column), zoomed-in view of the time series (center column), and power spectral density (right column) for five strain channels.}
    \label{fig:strain_runup}
\end{figure}

\begin{figure}[hbt!]
    \centering
    \includegraphics[width=\linewidth]{figures/strain_production_7_2.png}
    \caption{Strain estimation for a record from the production state (P4), shown with the full time series (left column), zoomed-in view of the time series (center column), and power spectral density (right column) for five strain channels.}
    \label{fig:strain_production}
\end{figure}


\begin{figure}[hbt!]
    \centering
    \includegraphics[width=0.9\linewidth]{figures/accuracy_boxplots.png}
    \caption{Correlation coefficient (CC) and mean relative error (MRE) across the 6 sub-soil strain channels in the idling (left), run-up (center), and production (right) states.}
    \label{fig:accuracystats}
\end{figure}

\indent It is expected that discrepancy in the predicted signal amplitude is largely caused by the well-known challenge of damping identification in the given application to offshore wind turbines, as the level of contribution of aerodynamic damping, in addition to hydrodynamic damping and soil damping in offshore contexts, is expected to fluctuate widely over time \cite{Devriendt2014, vanVondelen2021}. Variation in modal damping cannot be accounted for in the present implementation of the GPLFM where SSI-Cov is employed for modal identification. As with the majority of OMA techniques, it is implicitly assumed in SSI-Cov that 1) the system is linear time invariant, meaning the modal parameters remain unchanged over time; 2) the system is driven by unobserved ambient excitation which takes the form of white Gaussian noise; and 3) the system response is not colored by nonstructural damping sources. With SSI-Cov, robust estimates of modal parameters can only be achieved by using records from the idling state where loading conditions satisfy these criteria as closely as possible. \cite{Devriendt2014}. In reality, the OWT is a linear time-periodic (LTP) system in which harmonic environmental loads and rotor interaction forces introduce variation in the modal parameters of the OWT \cite{vanVondelen2021}. Therefore, the results of OMA obtained from data from the idling state do not fully reflect modal properties of the OWT in other states and are expected to result in an underestimation of modal damping levels. In particular, underestimation of the damping ratio at the first natural frequency, which controls low frequency response, leads to overestimation in the global amplitude of the predicted signal. 
\\
\indent Erroneous damping identification is evidenced by the modal force estimate for each of the operational states of the OWT, shown in Figure \ref{fig:modalforce}. While it is not possible to obtain direct measurements of the loads for validation of the input estimate, the frequency spectrum of the prediction provides insights to the joint relationship between inputs and states defined in the GPLFM, such that the input estimate may be used as a type of heuristic indicator of model error. First, it is clear that the magnitude of the input increases from the idling state to the production state, reflected by the difference in signal amplitude in the time series plots. Higher estimation uncertainty for higher modes is shown by the larger shaded area representing the $\pm3$ standard deviation bounds on the modal force prediction. Of particular interest is the frequency content estimated at the natural frequencies of the system, noted by dashed lines in the PSD plots. For the run-up state and production state, there is a visible dip in the frequency spectrum at the first natural frequency, which is an indicator of misidentification of the damping ratio in the first mode. In particular, when an underdamped model of the structural system predicts response amplitudes which exceed that of the measured acceleration response, the GPLFM compensates for the mismatch by reducing the power at that natural frequency in the frequency spectrum of the modal force. This effects suggests that in the joint prediction of latent inputs and states, GPLFM has some capacity to accommodate for model error by attributing its source to stochastic noise, reflected by adjustments to the predicted spectrum of external excitation to the system in order to fit the data provided. However, the input estimate has limited efficacy in accounting for other sources of systematic model error, such as a mismatch in the relationship between accelerations and strains defined in the mechanical model. Since congruence in the mode shapes identified from OMA and those of the mechanical model may only be evaluated at three points where accelerations are recorded, there is a high likelihood for deviation in the section of the mode shape profile falling below the mudline, leading to error in the extrapolation of latent states from measured accelerations. Nevertheless, these results in response estimation provide strong evidence that the GPLFM overcomes several challenges with reconstructing subsoil strain response compared to a similar study conducted using dual-band MD\&E for virtual sensing in \cite{Henkel2019}, which was not able to achieve the same level of accuracy in the frequency domain estimate and phase synchronization of the time domain estimate.


\begin{figure}[hbt!]
    \centering
    \includegraphics[width=0.95\linewidth]{figures/modalforce_comp.png}
    \caption{Estimation of the first three modal forces for records in the idling, run-up, and production states, shown with the full time series on the left and power spectral density on the right.}
    \label{fig:modalforce}
\end{figure}

\subsection{Sensitivity to sensor configuration}
\label{sec:measurevariants}

In this section, the performance of the GPLFM is examined under different  assumptions with respect to the availability and configuration of measurement data. The motivation is twofold: to determine an instrumentation scheme which aids in correcting for the overestimation of strain amplitude, and to determine the minimum instrumentation scheme which provides sufficient evidence to characterize the strain response and inform future measurement campaigns for fatigue monitoring of similarly sized OWTs. Considering the practicality of in-situ implementation, we assume only sensors at accessible locations above the waterline are used for virtual sensing of subsoil response. Strain estimation with the GPLFM is conducted with three sensor configurations, labeled as: 

\begin{itemize}
    \item A3, using only the acceleration channel in the RNA (+95.0m NAP);
    \item A3+A1, using two acceleration channels with one located in the RNA (+95.0m NAP) and the other at +4.50m NAP; and
    \item A3+S8, using the acceleration channel in the RNA (+95.0m NAP) and the strain gauge ring located at +6.95m NAP.
\end{itemize}

Validation of the strain estimate with Record P3 for the three sensor configurations is given in Figure \ref{fig:sensorvar_str}. For brevity, only strain channels 2, 6, and 7 are shown. It is observed that the A3 configuration is only capable of capturing the first modal response, with the accuracy at frequencies above the first natural frequency quickly tailing off. The introduction of the second acceleration channel in A3+A1 allows for accurate estimation of higher frequency contributions to the response. However, A3 and A3+A1 lead to a comparable level of overestimation of the signal amplitude due to the inability of the model to overcome sources of error, including uncertainty in local damping conditions in the foundation system, using acceleration-only measurements. In contrast, the inclusion of a strain channel in A3+S8 leads to substantial improvements in capturing the strain amplitude across channels, in exchange for reduced accuracy in the estimation of non-governing frequencies (falling between natural frequencies of the system). Compared to A3, the single strain channel in A3+S8 adds information on not only the high frequency response (such as the second modal response) in the tower section closer to the mudline, but also the amplitude of strains, thereby overcoming local damping error.
\\
\indent The quality of the strain estimate for different sensor configurations may be compared using fatigue assessment measures. In this work, the S-N diagram (also referred to as the Woehler curve) is constructed using a rainflow counting algorithm which translates time series data into the number of cycles counted within defined stress ranges \cite{Laanes2021, ASTM2017}. The S-N diagram of Figure \ref{fig:sensorvar_SN} shows the stress cycle counts of the true signal and the predicted signals at Str. 6, where the highest bending moments in the tower-foundation system are expected to occur, across four records from the production state. From this illustration, it is clear that configurations A3 and A3+A1 lead to an overestimation of the stress magnitude for lower cycle counts, which corresponds to low frequency oscillations. For the A3 case, the contribution of low-magnitude high frequency oscillations to fatigue stress accumulation is completely miscalculated. Figure \ref{fig:sensorvar_SN} affirms the A3+S8 configuration provides the best estimate of the stress-cycle curve. This result is further reflected in the damage equivalent load (DEL) measures of fatigue stress accumulation, where the deviation of the DEL of the prediction from that of the true strain response is 106.89\% for A3,  98.43\% for A3+A1, and 16.04\% for A3+S8. In other words, the fatigue estimate from the acceleration-only configurations A3 and A3+A1 is on the order of double that of the true fatigue damage accumulation. If implemented in practice for fatigue estimation, these acceleration-only configurations produce errors which would have large implications on the performance-based design of OWTs, resulting in overly conservative design measures which reduce the economic efficiency of offshore wind energy production.  


\begin{figure}[hbt!]
    \centering
    \includegraphics[width=\linewidth]{figures/sensor_comp_strain.png}
    \caption{Strain estimation with different sensor configurations for record P3.}
    \label{fig:sensorvar_str}
\end{figure}

\begin{figure}[hbt!]
    \centering
    \includegraphics[width=0.55\linewidth]{figures/sensor_comp_SN.png}
    \caption{S-N diagram with different sensor configurations across 4 samples in production state.}
    \label{fig:sensorvar_SN}
\end{figure}


\subsection{Sensitivity to model uncertainty}
\label{sec:modelvariants}

Previously, we demonstrated that the mathematical construction of the GPLFM attributes process noise to the unknown input. This section evaluates this effect using in-situ validation data, by imposing controlled model error in the mechanical model and evaluating the consequences on the strain estimation. In this case study, as well as in general for offshore wind turbines, one significant source of uncertainty is in properties of the foundation stiffness. Whereas the stiffness profile $k_{eff}(z)$ is derived from soil characterization data which is customarily collected at wind farm sites, the subsoil strain data used to tune the magnitude of the profile $\gamma$ in \cite{Versteijlen2017a} is generally not available. Without such data, there is the potential for substantial error in the assumed soil stiffness magnitude. We examine both extremes of introducing error in the soil stiffness scaling $\gamma$ by an order of magnitude above and below the baseline stiffness magnitude, with all other FE model parameters remaining constant. The following model variants representing different levels of model fidelity are compared: 

\begin{itemize}
    \item $\gamma = 0.078 $: limiting case of underestimating the soil stiffness scaling by an order of magnitude
    \item $\gamma = 0.78 $: baseline case of soil stiffness scaling to reflect the best fit to true in-situ conditions, determined in \cite{Versteijlen2017a}
    \item $\gamma = 7.8 $: limiting case of overestimating the soil stiffness scaling by an order of magnitude
\end{itemize}

The modal parameters of the three model variants are shown in Figure \ref{fig:modelvariants} and reported in Table \ref{tbl:modalparams_modelvar}, where the relative frequency error $e_f$ and MAC values are with respect to the modal parameters identified through OMA. It is observed that the lowest two mode shapes in each of the FA and SS directions, which dominate the OWT response, are relatively unaffected by scale differences in the soil stiffness profile. A larger effect is observed in the third SS mode. In both $e_f$ and MAC, the error appears greater in the case of underestimation of $\gamma$ than in the case of overestimation. In certain instances, the scale difference may result in a better fit with the modal parameters from OMA; for example, the case of overestimation gives a better fit of the SS2 mode shape, while the case of underestimation does so for the SS3 mode shape. 
\\
\indent Figure \ref{fig:modelvar_str} shows the strain estimate obtained with each model variant using the acceleration-only measurement setup. The difference in the accuracy of the strain estimate varies with the depth of the strain location: there is comparatively lower change in the strain estimate at higher elevations (Str. 6 and 7) than at the lowest elevation (Str. 2), where underestimation of stiffness leads to underdamped response and the overestimation of stiffness leads to overdamped response. The difference is non-uniform across channels due to the fact that there is significantly lower stiffness at the upper section of the profile, where $k_{eff}(z)$ is on the order of $10^7 \ N/m^3$, than at lower sections of the profile, where $k_{eff}(z)$ is on the order of $10^{11} \ N/m^3$ (see Figure \ref{fig:foundation}). The scale difference therefore leads to a greater change in subsoil properties at the greatest depths. While the scaling of the predicted signal is affected by the model error, the distribution of power over the frequencies remains relatively unchanged up to the second natural frequency across all channels. 
\\
\indent The difference in prediction by the model variants is more apparent in the modal force estimate, shown in Figure \ref{fig:modelvar_force}. Here, the discrepancy in the natural frequencies between the model variants is illustrated with dashed lines in the frequency spectrum. What is particularly notable is the dip in the frequency spectrum at the $j$th natural frequency for the $j$th modal force (first natural frequency in the first modal force, second natural frequency in the second modal force, etc.), where the natural frequencies differ for each model variant. These dips reflect adjustments to the frequency contributions of the input by the GPLFM in order to compensate for error between the dominant frequencies of the measured response with that of the predicted response by the perturbed model. For instance, where the GPLFM expects greater presence of the perturbed second natural frequency of the model than what exists in the measurement data, the contribution of that frequency in the corresponding input estimate is reduced accordingly. In this manner, the GPLFM accounts for model error through the stochastic input estimate, thus retaining the relatively robust accuracy observed in the strain estimate. 
\\
\indent 
To evaluate the effect of model fidelity as a function of the instrumentation scheme, Figure \ref{fig:modelvar_SN} shows the S-N diagram corresponding to four records from the production state, produced from the estimate of Str. 6, for two sensor configurations: the original scheme of using all three acceleration measurement channels (A1+A2+A3) and the best performing scheme from the investigation in Section \ref{sec:measurevariants} (A3+S8). For both instrumentation schemes, the three model variants result in relatively comparable stress-cycle curves, indicating relative invariance to perturbation in the model. As before, A3+S8 provides the most robust reproduction of the true stress-cycle curve while the acceleration-only scheme in A1+A2+A3 universally overestimates the stress magnitude of low-frequency oscillations. The improvement gained from the use of A3+S8 over the acceleration-only scheme is reflected in the DEL values reported in Table \ref{tbl:del}, where A3+S8 results in both lower error and lower variance in the DEL estimate across levels of model perturbation.  
    


\begin{figure}[hbt!]
    \centering
    \includegraphics[width=0.85\linewidth]{figures/modeshape_modelvar.png}
    \caption{Modal parameters of model variants.}
    \label{fig:modelvariants}
\end{figure}


\begin{table}[!htp]
    \centering
    \fontsize{8}{8}\selectfont
    \setlength{\tabcolsep}{8pt}
    \begin{tabular}{lllllllllllll}
        ~ & ~ & ~ & $\gamma$ = 0.078 & ~ & ~ & ~ & $\gamma$ = 0.78 & ~ & ~ & ~ & $\gamma$ = 7.8 & ~ \\ \hline
        mode & ~ & $\textit{f}_{fem}$ \ [Hz] & $\textit{e}_f$ [\%] & MAC & ~ & $\textit{f}_{fem}$ \ [Hz] & $\textit{e}_f$ [\%] & MAC & ~ & $\textit{f}_{fem}$ & $\textit{e}_f$ [\%] & MAC \\ \cline{1-1} \cline{3-5} \cline{7-9} \cline{11-13}
        FA1 & ~ & 0.279 & -5.82 & 0.999 & ~ & 0.296 & 0.01 & 0.999 & ~ & 0.307 & 3.55 & 0.999 \\ 
        SS1 & ~ & 0.269 & 6.70 & 0.990 & ~ & 0.289 & -0.12 & 0.996 & ~ & 0.300 & 3.74 & 0.998 \\ 
        FA2 & ~ & 1.636 & -12.47 & 0.991 & ~ & 1.874 & 0.23 & 0.999 & ~ & 2.041 & 9.18 & 0.994 \\ 
        SS2 & ~ & 1.665 & -11.66 & 0.886 & ~ & 1.883 & -0.09 & 0.939 & ~ & 2.025 & 7.44 & 0.964 \\ 
        SS3 & ~ & 4.486 & -1.69 & 0.987 & ~ & 5.193 & 13.79 & 0.743 & ~ & 5.745 & 25.88 & 0.310 \\ \hline
    \end{tabular}
\label{tbl:modalparams_modelvar}
\caption{Modal parameters from each model variant, where accuracy metrics $e_f$ and MAC are computed with respect to the modal parameters identified by OMA.}
\end{table}


\begin{figure}[hbt!]
    \centering
    \includegraphics[width=\linewidth]{figures/model_comp_strain.png}
    \caption{Strain estimation for varying levels of model fidelity using all acceleration measurement channels for record P2.}
    \label{fig:modelvar_str}
\end{figure}

\begin{figure}[hbt!]
    \centering
    \includegraphics[width=0.9\linewidth]{figures/model_comp_force.png}
    \caption{Modal force estimation for varying levels of model fidelity using all acceleration measurement channels for record P2. The natural frequencies of each model variant are noted by dashed lines in the frequency spectrum.}
    \label{fig:modelvar_force}
\end{figure}

\begin{figure}[hbt!]
    \centering
    \includegraphics[width=\linewidth]{figures/model_comp_SN.png}
    \caption{S-N diagram for varying levels of model fidelity for records from the production state, for sensor configuration A1+A2+A3 on the left and sensor configuration A3+S8 on the right.}
    \label{fig:modelvar_SN}
\end{figure}


\begin{table}[!htp]
\centering
\fontsize{8}{8}\selectfont
\setlength{\tabcolsep}{10pt}
    \begin{tabular}{llllll}
        ~ & \multicolumn{2}{c}{A1+A2+A3} & ~ & \multicolumn{2}{c}{A3+S8} \\ \hline
        ~ & DEL & error [\%] & ~ & DEL & error [\%] \\ \cline{2-3} \cline{5-6}
        true & 2.946 & -- & ~ & 2.946 & ~ \\ 
        gamma = 0.078 & 7.131 & 142.1\% & ~ & 4.295 & 45.8\% \\ 
        gamma = 0.78 & 6.652 & 125.8\% & ~ & 3.418 & 16.0\% \\ 
        gamma = 7.8 & 4.994 & 69.5\% & ~ & 2.227 & -24.4\% \\ \hline
    \end{tabular}
\caption{Damage equivalent load (DEL) values corresponding to Str. 6 of each model variant for two sensor configurations, A1+A2+A3 and A3+S8. The relative error is computed with respect to the true value of DEL.}
\label{tbl:del}
\end{table}


\section{Application to virtual sensing}
\label{sec:results}

This section discusses validation of the GPLFM for virtual sensing using operational data from the W27 turbine. First, the fitting of kernel hyperparameters using the approach of Equation \eqref{eq:mle} is presented. The results of dynamic strain and modal force estimation using acceleration-only measurement data are shown for sample records from each of the turbine's operational states. Next, the sensitivity of the GPLFM to the sensor configuration is examined to determine the minimum sufficient measurement setup for characterizing the dynamic strain response of the OWT. Finally, perturbation to the foundation stiffness profile is introduced to evaluate the ability of the GPLFM to accommodate model uncertainty. The performance of the GPLFM is calculated using the correlation coefficient (CC) and mean relative error (MRE) of the predicted strain time series in relation to the validation measurements, as well as by converting the strain estimate into a Damage Equivalent Load (DEL) \cite{Laanes2021, ASTM2017} to measure the impact of estimation error on fatigue assessment. 
\\
\indent In this application of the GPLFM, the finite element model of the W27 turbine is translated into a reduced-order model retaining the first $n_{\mathrm{m}}=3$ modes in each of the FA and SS directions of the structure. The stochastic model of the input assumes the three modal forces are independent and identically characterized by a Mat\'ern kernel with $\nu=5/2$, such that $\kappa_1(\tau) = \kappa_2(\tau) = \kappa_3(\tau)$. The corresponding augmented model is therefore driven by two hyperparameters, $\theta=[\alpha, l_s]$. 
\\
\indent The acceleration-only measurement dataset consists of acceleration time series transformed into FA and SS coordinates at three elevations of the structure; these elevations are referred to as A1 (combination of orthogonal Acc. 1 and Acc. 2 at +4.50m NAP), A2 (uniaxial Acc. 3 at +42.50m NAP), and A3 (accelerometer in the RNA at +95.0m NAP). The validation dataset consists of strain time series recorded at position A for FA motion and position B for SS motion within the rings at six depths below the mudline, referred to by their corresponding ring numbers as Str. 2-7 (refer to Figure \ref{fig:sensorsetup} for an illustration). Validation is shown for eight records from each of the operating states of the turbine, with each record identified by a numeral preceded by the letter 'D' for the idling state, 'R' for the run-up state, and 'P' for the production state. The records are selected with the intention to encompass a wide range of wind and rotor conditions in each operating state, as shown on the wind speed-rotor speed state curve in Figure \ref{fig:statecurve}.


\begin{figure}[hbt!]
    \centering
    \includegraphics[width=0.6\linewidth]{figures/statecurve.png}
    \caption{Wind speed-rotor speed state curve of 10-minute records taken over the period of March 2016 to June 2017 for the W27 turbine. The records selected for validation in each of the operating states are labeled. }
    \label{fig:statecurve}
\end{figure}


\input{04.1_PriorTuning}
\input{04.2_ResponseEstimation}
\input{04.3_SensorVariants}
\input{04.4_ModelVariants}
\section{Conclusions}
\label{sec:conclusions}

This work is among the first studies to provide in-situ validation of the GPLFM as a virtual sensing technique. In contrast to other Kalman filter-based techniques to virtual sensing, the GPLFM characterizes the unknown input as a Gaussian process which can describe a greater breadth of time-varying behavior than the standard assumption of white noise input. The joint posterior distribution of latent inputs and states is determined by Gaussian process regression of acceleration measurements, which is performed sequentially in augmented state-space using Kalman filtering/RTS smoothing. Rather than using heuristic or ad hoc tuning procedures, the process noise covariance in the dynamics model is solved for in a data-driven fashion, in which the structure of the GPLFM is utilized to define the joint relationship between input and state variables and the kernel hyperparameters governing the process noise covariance matrix are determined based on a priori fitting to acceleration data. These properties of the GPLFM lend it notable advantages over other Kalman filter-based techniques such as the AKF, AKF-dm, DKF, and GDF in numerical experiments for joint input-state estimation, as shown in \cite{Nayek2019}.
\\
\indent In this contribution, the GPLFM is applied to the fatigue monitoring of offshore wind turbines. The validation data is taken from a monitoring campaign of an operating nearshore wind turbine from the Westermeerwind Park in the Netherlands, where six strain gauge rings embedded along the profile of the monopile foundation provide subsoil strain measurements for validation. This dataset presents the opportunity to evaluate the GPLFM for its ability to address the challenge in fatigue monitoring of OWTs with high model uncertainty in the soil-structure system of the monopile foundation. First using acceleration-only data from an offshore wind turbine, the GPLFM implemented with a Mat\'ern covariance kernel is able to reconstruct the frequency spectrum of the subsoil strain response with excellent accuracy across operating conditions of the OWT. In the time domain, overestimation of the strain amplitude is attributed to the sensitivity of subsoil strains to local soil damping conditions, as well as the inability to identify and account for dynamic contributions to modal damping. In order to test the sensitivity of the GPLFM to the level of model fidelity, perturbation to the scale of the foundation stiffness model is introduced. It is demonstrated that with induced errors in the natural frequencies and mode shapes of the mechanical model, the strain estimate by the GPLFM remains relatively insensitive, particularly at depths of interest near the mudline where the highest bending moments are expected to occur. This result illustrates that the GPLFM has the ability to compensate for such model error by attributing their effects to qualities of the stochastic input estimate, thereby achieving robustness in the strain estimate. Through a study of the influence of the sensor configuration on the performance of the GPLFM, it is evident that estimation bias resulting from model error may be systematically accounted for through the inclusion of a strain channel in the measurement data. Not only does the strain data aid in correcting the predicted signal amplitude along all depths of the foundation, but it may be used to supply important information on the quasi-static behavior of the strain response which cannot be captured by high-frequency acceleration measurements alone \cite{Noppe2016}. The best performing instrumentation scheme minimally consists of a single acceleration channel and single strain channel, where the GPLFM achieves a relative error of 16\% in the DEL, mean correlation coefficient of 0.968, and mean relative error of 2.28\% in the estimated subsoil strain time series across records from the production state of the OWT. Overall, the GPLFM attains strong robustness and transferrability properties in estimating the subsoil strain response to variable operational and environmental loads.
\\
\indent It is expected that the performance of the GPLFM for subsoil strain estimation would be improved with the explicit consideration of dynamic evolution of modal properties of the structure, such as modal damping ratios. This work demonstrates the limitation of employing traditional OMA techniques, whose formulation relies on the assumption of white noise input, to the modal identification of offshore wind turbines in the run-up and production states, where wind and wave loads and the rotor motion of the turbine introduce harmonic components to the excitation. Due to the treatment of modal parameters of the mechanical model as time-invariant, strain predictions by the GPLFM exhibit an overestimation of the signal amplitude as a result of underestimation of damping in the first structural mode. While discrepancy in the time series amplitude is mitigated in part by real-time compensation for model error in the input estimate, the stochastic component of the GPLFM which models unknown inputs is not fully sufficient to correct for process error in the mechanical component of the GPLFM. Therefore, the application of virtual sensing to offshore wind turbines either requires a heterogeneous dataset of acceleration and strain measurements to curtail model uncertainty, or warrants additional developments in online modal identification techniques which account for non-white excitation in the system. For the latter, the authors recommend further evaluation of proposed adjustments and alternatives to damping identification techniques for LTP systems which are reviewed in \cite{vanVondelen2021}. In \cite{Rogers2020}, the concept of modal identification is treated with a Bayesian approach, where inference of system parameters such as mass, stiffness, and damping is performed jointly with estimation of the unknown ambient excitation by means of the GPLFM. This approach enables posterior state inference and modal identification to be performed jointly with a fully cohesive model of the input and may be applied across a greater range of contexts compared to traditional OMA methods. In general, the representation of time-varying properties of the system in the mechanical model is expected to produce among the greatest gains in the performance of the GPFLM for the virtual sensing of structures in offshore environments. 
\\
\indent There is potential for further improvements to the strain estimation through model choices in the GPLFM. In particular, methods of addressing model uncertainty through the definition of process noise in the latent force model would be expected to reduce error in the time series prediction. As discussed previously, process noise covariance could be defined in the spectral density of the augmented white noise term of the GPLFM, and the scale parameter of the corresponding covariance matrix could be solved for jointly in the optimization routine for parameter inference. Secondly, there are opportunities to investigate the role of the choice of kernel function in the GPLFM for response data which is driven by inputs of varying degrees of periodicity and temporal correlation. Beyond model choices, the results of this work promote further studies on on efficient and practical sensor placement schemes in monopile-based offshore wind turbines, to select the most informative data for the purpose of virtual sensing-based fatigue monitoring. 


\section*{CRediT authorship contribution statement}

\textbf{J. Zou:} Methodology, Data curation, Software, Formal analysis, Investigation, Visualization, Writing - original draft. \textbf{E. Lourens:} Conceptualization, Supervision, Writing - review \& editing. \textbf{A. Cicirello:} Supervision, Writing - review \& editing.

\section*{Acknowledgements}

This research was funded by a Fulbright Student Research Grant awarded by the Stichting Fulbright Commission of the Netherlands. The data used in this paper is obtained from the DISSTINCT project, a collaboration between Siemens Wind Power, Fugro, Delft University of Technology, and DNV-GL which was in part funded by the Dutch government (project number TKIW02001). The authors gratefully acknowledge the people at Siemens-Gamesa Renewable Energy for providing the data, Joaquin Bilbao-Neiva for sharing code for the implementation of the SSI-Cov/OPTICS algorithm, and Mikk Laanes for sharing data preprocessing code.
\section{}


\subsection{Transformation of lateral displacements to axial strain}
\label{sec:app_disptransf}

From elastic beam theory, the lateral deformation of a 2D beam element $v$ is represented with a third order polynomial of the form: 

\begin{equation}
    v(\xi) = a_0 + a_1\xi + a_2\xi^2 + a_3\xi^3
\end{equation}

$\xi = \frac{x}{L}$ is the normalized axial distance along length $L$ of the element. The unknown coefficients $a_0, a_1, a_2, a_3$ can be solved for by enforcing boundary conditions which relate the deformation to the vector of degrees of freedom of the element $\mathbf{d}_{\mathrm{e}}$, one translational and one rotational degree of freedom at each end:

\begin{equation}
    \label{eq:strain0}
    v(\xi) = 
    \begin{bmatrix}
    N_1(\xi) & N_2(\xi) & N_3(\xi) & N_4(\xi)
    \end{bmatrix}
    \begin{bmatrix}
    v_1 \\ \theta_1 \\ v_2 \\ \theta_2
    \end{bmatrix} = 
    \mathbf{N}_{\mathrm{e}} \mathbf{d}_{\mathrm{e}}
\end{equation}

The solution provides the element shape functions $\mathbf{N}_{\mathrm{e}}$, which approximate the continuous spatial distribution of deformation between finite element nodes. The derivations of shape functions for an Euler-Bernoulli beam element and Timoshenko beam element are detailed in \cite{Liu2003} and \cite{Bazoune2003}, respectively. 

For small deformations, the axial displacement $u$ may be approximated by the rotational deformation $\theta_z$ and distance from the element centroidal axis $y$ as: 

\begin{equation}
    \label{eq:strain1}
    u = -y\theta_z
\end{equation}

The axial strain $\epsilon_{xx}$ is related to the rate of change of the axial displacement with respect to $x$, which is given as: 

\begin{equation}
    \label{eq:strain2}
    \epsilon_{xx} = \diffp{u}{x} = -y\diffp{\theta_z}{x} = -y\diffp{v}{{x^2}}
\end{equation}

Substituting the expression for $v$ from Equation \eqref{eq:strain0}, Equation \eqref{eq:strain2} becomes:

\begin{equation}
    \label{eq:strain3}
    \epsilon_{xx} = -\frac{y}{L^2}\mathbf{N}_{\mathrm{e}}''(\xi)\mathbf{d}_{\mathrm{e}} = 
    \mathbf{B}_{\mathrm{e}}(\xi)\mathbf{d}_{\mathrm{e}}
\end{equation}

where $\mathbf{B}_{\mathrm{e}}(\xi) = \begin{bmatrix} B_1(\xi) \quad B_2(\xi) \quad B_3(\xi) \quad B_4(\xi) \end{bmatrix}$ is the element strain matrix.
\\
\indent For a finite element model with a chain of $n_{\mathrm{e}}$ identical elements and free end conditions, the global strain vector can be computed as the  matrix product $\mathbf{\epsilon} = \mathbf{Td}$, which is written in expanded form as follows:

\begin{equation}
    \label{eq:strain4}
    \begin{bmatrix}
        \epsilon_{xx}(\xi_1) \\
        \epsilon_{xx}(\xi_2) \\
        \vdots \\
        \epsilon_{xx}(\xi_{n_{\mathrm{e}}})
    \end{bmatrix} = \begin{bmatrix}
        B_1(\xi_1) & B_2(\xi_1) & B_3(\xi_1) & B_4(\xi_1) &  &  &  &  &   \\
        &  & B_1(\xi_2) & B_2(\xi_2) & B_3(\xi_2) & B_4(\xi_2) &  &  &   \\
        &  &  &  &  \ddots &  &  &  & \\
        &  &  &  &  & B_1(\xi_{n_{\mathrm{e}}}) & B_2(\xi_{n_{\mathrm{e}}}) & B_3(\xi_{n_{\mathrm{e}}}) & B_4(\xi_{n_{\mathrm{e}}}) 
    \end{bmatrix}
    \begin{bmatrix}
        v_1 \\ \theta_1 \\
        v_2 \\ \theta_2 \\
        \vdots \\
        v_{n_{\mathrm{e}}} \\ \theta_{n_{\mathrm{e}}} \\
    \end{bmatrix}
\end{equation}

Each entry in the global strain vector provides the strain estimate at a distance $\xi_i$ along the $i$th element, for $i = 1,...,n_{\mathrm{e}}$. If the elements are non-identical, then each row of the global strain matrix $\mathbf{T}$ corresponds to the strain functions of the $i$th element.


\subsection{Kalman filtering and RTS smoothing equations}
\label{sec:app_KFKS}

Given a discrete linear time-invariant state-space model of the form

\begin{subequations}
    \begin{align}
        \mathbf{z}_k = \mathbf{F}\mathbf{z}_{k-1} + \mathbf{w}_{k-1} 
        \\
        \mathbf{y}_k = \mathbf{H}\mathbf{z}_k + \mathbf{v}_k
    \end{align}
    

where the state vector has a prior distribution $\mathbf{z}_0 \sim N(\hat{\mathbf{z}}_{0|0}, \mathbf{P}_{0|0})$ and the additive white Gaussian noise terms have covariances $\mathbb{E}[\mathbf{w}_k\mathbf{w}_l^\textup{T}] = \mathbf{Q}$, $\mathbb{E}[\mathbf{v}_k\mathbf{v}_l^\textup{T}] =\mathbf{R}$, and $\mathbb{E}[\mathbf{w}_k\mathbf{v}_l^\textup{T}] = \mathbf{0}$. 
\\
\indent Let $\mathbf{z}_{p|q}$ denote the state predicted at time step $p$ given all time steps $1:q$. The Kalman filter of Algorithm \ref{alg:KF} solves for the posterior filtering distribution, where the mean and covariance in $\mathbf{z}_{k|k} \sim N(\hat{\mathbf{z}}_{k|k}, \mathbf{P}_{k|k})$ is predicted from a "forward pass" through time steps $k=1,...,N$ and updated at time steps where measurements are available \cite{Kalman1960}. The RTS smoother of Algorithm \ref{alg:KS} solves for the posterior smoothing distribution, where the predicted mean and covariance from the Kalman filter (denoted with superscript "-" for clarity) is updated with a fixed-interval "backward pass" through time steps $k=N-1,...,1$ \cite{RTS1965}. The final distribution $\mathbf{z}_{k|N} \sim N(\hat{\mathbf{z}}_{k|N}, \mathbf{P}_{k|N})$ constitutes the complete posterior distribution of the state given all measurements in a sequence of $N$ steps to estimate. From this distribution, a posterior estimate of the measurement states may be recovered as $\mathbf{y}_{k|N} \sim N(\hat{\mathbf{y}}_{k|N}, \mathbf{P}_{k|N}^y)$.

\begin{algorithm}[H]
\caption{Kalman Filter}\label{alg:KF}
\begin{algorithmic}
\For{$k = 1,...,N$}
    \State Filter prediction equations:
    \State $\hat{\mathbf{z}}^-_{k|k-1} = \mathbf{F}\hat{\mathbf{z}}_{k-1|k-1}$ 
    \State $\mathbf{P}^-_{k|k-1} = \mathbf{F}\mathbf{P}_{k-1|k-1}\mathbf{F}^\textup{T} + \mathbf{Q}$  
\If{$\mathbf{y}_k$ is observed}
    \State Filter update equations:
    \State $\mathbf{e}_k = \mathbf{y}_k - \mathbf{H}\hat{\mathbf{z}}^-_{k|k-1}$ \Comment prediction error
    \State $\mathbf{S}_k = \mathbf{H}\mathbf{P}^-_{k|k-1}\mathbf{H}^\textup{T} + \mathbf{R}$ 
    \State $\mathbf{K}_k = \mathbf{P}^-_{k|k-1}\mathbf{H}^\textup{T}(\mathbf{S}_k)^{-1}$ \Comment Kalman gain
    \State $\hat{\mathbf{z}}_{k|k} = \hat{\mathbf{z}}^-_{k|k-1} - \mathbf{K}_k\mathbf{e}_k$
    \State $\mathbf{P}_{k|k} = \mathbf{P}^-_{k|k-1} - \mathbf{K}_k\mathbf{H}\mathbf{P}^-_{k|k-1}$
\EndIf
\EndFor
\end{algorithmic}
\end{algorithm}

\begin{algorithm}[H]
\caption{Rausch-Tung-Striebel (RTS) Smoother}\label{alg:KS}
\begin{algorithmic}
\For{$k = N-1,...,1$}
    \State Smoothing update equations:
    \State $\mathbf{G}_k = \mathbf{P}_{k|k}\mathbf{F}^\textup{T} \left( \mathbf{P}^-_{k+1|k} \right) ^{-1}$
    \State $\hat{\mathbf{z}}_{k|N} = \hat{\mathbf{z}}_{k|k} + \mathbf{G}_k \left( \hat{\mathbf{z}}_{k+1|N}-\hat{\mathbf{z}}^-_{k+1|k} \right)$ 
    \State $\mathbf{P}_{k|N} = \mathbf{P}_{k|k} + \mathbf{G}_k \left( \mathbf{P}_{k+1|N} - \mathbf{P}^-_{k+1|k} \right) \mathbf{G}_k^\textup{T}$  
    
    \State Moments of measurement states $\mathbf{y}$:
    \State $\hat{\mathbf{y}}_{k|N} = \mathbf{H}\hat{\mathbf{z}}_{k|N}$ 
    \State $\mathbf{P}_{k|N}^y = \mathbf{H}\mathbf{P}_{k|N}\mathbf{H}^\textup{T} + \mathbf{R}$ 
\EndFor
\end{algorithmic}
\end{algorithm}

\end{subequations}





\subsection{SSM for Mat\'ern kernel}
\label{sec:app_matern}

The Mat\'ern family is a class of covariance functions of the form

\begin{equation}
    \kappa(\tau; \nu, \alpha, l_s) = \alpha^2 \frac{2^{(1-\nu)}}{\Gamma(\nu)} \left( \frac{\sqrt{2\nu}}{l_s} \tau \right)^\nu K_\nu \left( \frac{\sqrt{2\nu}}{l_s} \tau \right)
\end{equation}


where $\Gamma(\nu)$ is the Gamma function and $K_\nu$ is the modified Bessel function of the second kind. The Mat\'ern kernel defines the covariance between two points $t, t'$ of a process to be a function of only the distance $\tau = t' - t$ between them, resulting in stationary stochastic processes. The hyperparameters influence qualities of the covariance relationship derived from the distance: the vertical length scale $\alpha$ controls the magnitude of variance at $\tau=0$, the horizontal length scale $l_s$ controls the distance $\tau$ over which covariance diminishes to zero, and the parameter $\nu$ controls the smoothness of the resulting process. Note that in its limits, a Mat\'ern kernel with $l_s \rightarrow 0$ obtains the Dirac delta function for the covariance, which corresponds to a white noise kernel, and $\nu \rightarrow \infty$ obtains the squared exponential covariance function corresponding to infinitely differentiable processes \cite{Rasmussen2006, Hartikainen2010}. 
\\
\indent This section presents the form of the Mat\'ern kernel with smoothing parameter $\nu=p+1/2=5/2$, where $p=2$, leading to processes which are twice differentiable. The associated hyperparameters are then $\theta=[\alpha, l_s]$. The covariance function becomes:

\begin{equation}
    \kappa_{5/2}(\tau; \alpha, l_s) = \alpha^2 \left( 1 + \frac{\sqrt{5}}{l_s}\tau + \frac{5\tau^2}{3l_s^2} \right) \textup{exp} \left(- \frac{\sqrt{5}}{l_s} \tau \right)
\end{equation}

A Gaussian process $h(t) \sim GP(0, \kappa_{5/2}(\tau))$ may be expressed as a third-order linear stochastic differential equation of the form: 

\begin{equation}
    \label{eq:maternSDE}
    \frac{d^{3}h(t)}{dt^3} +
    3 \lambda \frac{d^{2}h(t)}{dt^{2}} +
    3 \lambda^2 \frac{dh(t)}{dt} +
    \lambda^3 h(t) = w(t)
\end{equation}

where $\lambda = \sqrt{5}/l_s$. The system is driven by white noise with spectral density:

\begin{equation}
     q_w = S_w(\omega) = \frac{400 \sqrt{5} \alpha^2}{3l_s^5}
\end{equation}

Consequently, we can use the form of the Mat\'ern kernel to derive the spectral density of $h(t)$ as $S_h(\omega) = q_w (\lambda^2 + \omega^2)^{-(p+1)}$.
\\
\indent The companion form of \ref{eq:maternSDE} gives rise to the state-space model of process $h(t)$, with the continuous-time matrices \cite{Nayek2019}: 


\begin{equation}
    \label{eq:maternforcemodel}
    \mathbf{F}_{c} = \begin{bmatrix} 0  & 1 & 0 \\
    0 & 0 & 1 \\
    -\lambda^3    & -3\lambda^2   & -3\lambda \\ \end{bmatrix}, \quad 
    \mathbf{L}_{c} = \begin{bmatrix} 0 \\ 0 \\ 1 \end{bmatrix}, \quad 
    \mathbf{H}_{c} = \begin{bmatrix} 1 & 0 & 0 \end{bmatrix}
\end{equation}



\subsection{Hellinger distance}
\label{sec:app_hellinger}

The Hellinger distance is a type of f-divergence which integrates the error between two probability densities $p(x)$ and $q(x)$ as follows: 

\begin{equation}
    \label{eq:hellingercont}
    H^2(p, q) = \frac{1}{2} \int \left( \sqrt{p(x)} - \sqrt{q(x)} \right)^2 dx = 1 - \int \sqrt{p(x) \, q(x)} \ dx
\end{equation}

In the discrete case, the Hellinger distance between the probability distributions $p = \{ p_1, ..., p_m \}$ and $q = \{q_1, ..., q_m \}$ may be computed as a Riemann sum over $m$ subdivisions: 

\begin{equation}
    \label{eq:hellingerdiscr}
    H(p, q) = 1 - \sum_{i=1}^{m} \sqrt{p_i \, q_i}
\end{equation}



