Virtual sensing in an onshore wind turbine tower using a Gaussian process latent force model

Abstract

Wind turbine towers are subjected to highly varying internal loads, characterized by large uncertainty. The uncertainty stems from many factors, including what the actual wind fields experienced over time will be, modeling uncertainties given the various operational states of the turbine with and without controller interaction, the influence of aerodynamic damping, and so forth. To monitor the true experienced loading and assess the fatigue, strain sensors can be installed at fatigue-critical locations on the turbine structure. A more cost-effective and practical solution is to predict the strain response of the structure based only on a number of acceleration measurements. In this contribution, an approach is followed where the dynamic strains in an existing onshore wind turbine tower are predicted using a Gaussian process latent force model. By employing this model, both the applied dynamic loading and strain response are estimated based on the acceleration data. The predicted dynamic strains are validated using strain gauges installed near the bottom of the tower. Fatigue is subsequently assessed by comparing the damage equivalent loads calculated with the predicted as opposed to the measured strains. The results confirm the usefulness of the method for continuous tracking of fatigue life consumption in onshore wind turbine towers.

Impact Statement

As the first large-scale wind farms reach their end-of-lifetime, attention is drawn to the importance of having accurate knowledge on structural reserves that can enable wind farm operators to make reliable decisions about lifetime extension in the future. In this contribution, a methodology for fatigue load monitoring on the basis of a sensor network is proposed and validated. The loads exciting the tower are estimated in conjunction with its dynamic response using a Gaussian process latent force model. Estimates of the full-field strain response can be derived, and therewith fatigue consumption at critical locations in the structure. The results confirm the potential of the proposed methodology for continuous tracking of fatigue life consumption in the structures supporting wind turbine rotors.

1. Introduction

Towers of onshore wind turbines are highly susceptible to fatigue damage given the cyclic wind loading, where the cyclic nature of the experienced wind loading is a consequence of rotational sampling of the wind field by the rotor. In recent years, a lot of effort has been spent on indirect measurement of the true fatigue life consumption by means of so-called virtual sensing techniques. Virtual sensing refers to the prediction of response quantities at unmeasured, but critical locations on a structure, based on limited data from sensors installed at easily accessible locations. In the case of wind turbines, also data from the Supervisory Control and Data Acquisition (SCADA) system, present in all modern wind turbines, can be utilized for the purpose of virtual sensing.

Several different techniques for virtual sensing have been proposed. A distinction can be made between data-driven models, physics-based models, and hybrid techniques. Data-driven techniques use artificial intelligence to train black-box models from measured or simulated data. For example, Cosack (2010), Smolka and Cheng (2013), Vera-Tudela and Kühn (2017), and Schröder et al. (2018) used artificial neural networks with SCADA as input to estimate damage equivalent loads at different locations of the structure. Ziegler et al. (2017, 2019) applied a k-nearest neighbor regression algorithm to extrapolate strain measurements along a monopile substructure. Recently, d N Santos et al. (2021) estimate thrust loads from 1s SCADA data through artificial neural networks, which are then combined with high-resolution acceleration measurements to estimate tower loads. Gulgec et al. (2020) apply neural networks to estimate time-history of strain responses directly from acceleration measurements. Other types of machine learning have also been applied in the context of virtual sensing, such as a Gaussian process regression by Avendaño-Valencia et al. (2021).

In the physics-based approach, a complete numerical model of the asset is used for virtual sensing. This requires an update of the model with a limited set of measurement data. Augustyn et al. (2020) recently used in situ acceleration measurements from a jacket substructure of an offshore wind turbine to update its numerical counterpart. The updated model can then be used for aero-hydro-elastic simulations to obtain loads on unmeasured locations. Currently, there is a trend toward framing these concepts as “digital twins.” A digital twin refers to a digital replicate of the physical asset continuously updated with real-time measurements (Wagg et al., 2020). Recent works promoting the concept of digital twins for wind turbines have been published by Branlard et al. (2020), Schröder (2020), and Augustyn (2021).

In hybrid schemes, measurements are combined with simplified numerical models of the asset for virtual sensing. This includes model-based extrapolation of responses and techniques where the forces driving the response are also estimated. In the former category, one finds the modal decomposition and expansion approach (MD&E), first applied in the field of (offshore) wind turbine foundation monitoring by Iliopoulos et al. (2014, 2016) and Noppe et al. (2016). Using MD&E, a classical inverse problem is solved where modal coordinates of the measured response are estimated through a pseudo-inverse of the model-based eigenvector matrix. To improve accuracy, a system identification experiment followed by a model update can be performed to obtain the eigenvectors. Also, the inverse problems can be solved separately for different frequency bands using multi-band MD&E (Iliopoulos et al., 2017), where the set of eigenvectors used in the inversion can be more carefully tailored to the expected response from each mode within the frequency band.

The second category of virtual sensing techniques encompasses the Kalman filter based methods (Kalman, 1960), in numerous variants. In this case, the response extrapolation is performed by simultaneously estimating the driving forces and the states of the system, where the latter is often defined to be the modal coordinates. This is done in a stochastic setting, accounting for both modeling and measurement errors, and applying principles of minimum-variance unbiased estimation. In this category falls, in no particular order, the augmented Kalman filter (Lourens, Reynders, et al., 2012), the dual Kalman filter (Eftekhar Azam et al., 2015; Tatsis and Lourens, 2016), and the joint input-response estimator (Gillijns and De Moor, 2007; Lourens, Papadimitriou, et al., 2012; Maes et al., 2016). Although strongly related in terms of their basic principles, these variants differ in their modeling of the input force. In the dual and augmented filter, the force is modeled as a random walk, and its variance is tuned based on the measurement data. In the joint input-response estimator, however, no assumptions on the evolution of the force in time are made, and its time history is found solely through minimum-variance unbiased estimation in the context of a set of uncertain model equations. Although excellent results have been obtained with the Kalman filter based methods in various applications, both in the laboratory and full-scale, a major drawback remains their instability. Especially in the case of acceleration-only measurements, the more classical filters require modifications to improve their stability (Eftekhar Azam et al., 2015; Naets et al., 2015). Often these modifications come at a cost, for example, the need to include additional parameters that require a cumbersome and fragile tuning process.

In this contribution, the virtual sensing problem is solved using a Gaussian process latent force model (GPLFM). The GPLFM relieves instability problems and/or fragile tuning processes that afflict the Kalman filter based methods, while at the same time preserving the benefits of a simultaneous estimation of both the driving forces and the system states. The GPLFM was first applied to the joint input-state estimation problem by Nayek et al. (2019) and extended for the solution of joint input-state-parameter identification by Rogers et al. (2020), both drawing on work by Hartikainen and Särkkä (2010), Särkkä (2013), and Särkkä and Solin (2019). Unlike the Kalman filter based methods, it models the driving forces as Gaussian processes, driven by specified covariance functions. The parameters defining the covariance functions can be optimized for using the measured data, resulting in a flexible, yet robust framework for incorporating prior information on the time-varying behavior of the force.

In the following sections, an application of the GPLFM for virtual sensing on a full-scale wind turbine tower is presented. Section 2 provides the theoretical framework, followed by the full-scale case study in Section 3, including the validation through strain gauges installed near the bottom of the tower. For finding the parameters that define the covariance functions of the unknown forcing, a new procedure exploiting the physics underlying the GPLFM is suggested in Section 2.6 and applied in Section 3.3. It is noted here that the estimation of these unknown parameters is commonly realized by finding the parameters such that the posterior distribution is maximized. More theoretical background can be found in Särkkä (2013), where several strategies for parameter estimation are discussed, such as maximum a posteriori (applied in Nayek et al., 2019) and Markov Chain Monte Carlo (applied in Rogers et al., 2020). The downside of these approaches is the need for a robust and expensive optimization process, and furthermore, the loss of physical relevance of this process. It is shown that the suggested simpler, more intuitive procedure can be effective for finding the hyperparameters needed to characterize the covariance of the unknown forces.

2. Theoretical Framework

The GPLFM is built by combining a mechanical model of the wind turbine structure and a stochastic model of the unknown load. Both models are expressed in state-space form and are subsequently combined into an augmented state-space model. This augmented model is a stochastic representation of the system, including the response and the load. The following subsections describe each model and show how the resulting augmented model is handled. In particular, the computation of the prior and the posterior as well as the estimation of the hyperparameters of the model are presented.

2.1 Mechanical model

Consider the governing equation of motion of a linear system discretized in space,

$$ \mathbf{M}\ddot{\boldsymbol{u}}(t)+\mathbf{C}\dot{\boldsymbol{u}}(t)+\mathbf{K}\boldsymbol{u}(t)={\mathbf{S}}_{\boldsymbol{p}}\boldsymbol{p}(t), $$

where the mass, damping, and stiffness matrices are represented by $ \mathbf{M},\mathbf{C},\mathbf{K}\in {\mathrm{\mathbb{R}}}^{n_u\times {n}_u} $ , respectively, and the load locations are defined through $ {\mathbf{S}}_{\boldsymbol{p}}\in {\mathrm{\mathbb{R}}}^{n_u\times {n}_p} $ . The response of the system is defined as $ \boldsymbol{u}(t)\in {\mathrm{\mathbb{R}}}^{n_u} $ , and the loads as $ \boldsymbol{p}(t)\in {\mathrm{\mathbb{R}}}^{n_p} $ . The mechanical model employed in this work considers only a single load: $ {n}_p=1 $ .

The classical modal representation of this equation of motion reads

$$ \mathbf{I}{\ddot{\boldsymbol{u}}}_{\boldsymbol{m}}(t)+\boldsymbol{\Gamma} {\dot{\boldsymbol{u}}}_{\boldsymbol{m}}(t)+{\boldsymbol{\Omega}}^2{\boldsymbol{u}}_{\boldsymbol{m}}(t)={\boldsymbol{\Phi}}^{\mathrm{T}}{\boldsymbol{S}}_{\boldsymbol{p}}p(t), $$

with $ {\boldsymbol{u}}_{\boldsymbol{m}}\in {\mathrm{\mathbb{R}}}^{n_u} $ the modal coordinates, $ \boldsymbol{\Phi} \in {\mathrm{\mathbb{R}}}^{n_u\times {n}_u} $ the mass-normalized eigenvector matrix (i.e., $ {\boldsymbol{\Phi}}^T\mathbf{M}\boldsymbol{\Phi } =\mathbf{I} $ ), and $ \boldsymbol{\Gamma} \in {\mathrm{\mathbb{R}}}^{n_u\times {n}_u} $ and $ \boldsymbol{\Omega} \in {\mathrm{\mathbb{R}}}^{n_u\times {n}_u} $ diagonal matrices formed by the natural frequencies $ {\omega}_n $ and damping ratios $ {\zeta}_n $ so that $ \boldsymbol{\Omega} =\operatorname{diag}\left({\omega}_n\right) $ and $ \boldsymbol{\Gamma} =\operatorname{diag}\left(2{\zeta}_n{\omega}_n\right) $ , respectively.

A subset of modes $ {n}_m\hskip0.35em \le \hskip0.35em {n}_u $ is assumed to be sufficient to accurately approximate the response. In other words, modal reduction is employed as

$$ \boldsymbol{u}(t)\hskip0.35em \approx \hskip0.35em {\boldsymbol{\Phi}}_{\boldsymbol{r}}{\boldsymbol{u}}_{\boldsymbol{mr}}(t), $$

with the reduced mode shape matrix $ {\boldsymbol{\Phi}}_r\in {\mathrm{\mathbb{R}}}^{n_u\times {n}_m} $ and modal coordinates $ {\boldsymbol{u}}_{\boldsymbol{mr}}\in {\mathrm{\mathbb{R}}}^{n_m} $ . The reduced system matrices are written as $ {\boldsymbol{\Omega}}_r\in {\mathrm{\mathbb{R}}}^{n_m\times {n}_m} $ and $ {\boldsymbol{\Gamma}}_r\in {\mathrm{\mathbb{R}}}^{n_m\times {n}_m} $ .

Finally, the reduced-order modal model is written in state-space form

(1) $$ {\displaystyle \begin{array}{c}\underset{{\dot{\boldsymbol{x}}}_{\boldsymbol{m}}(t)}{\underbrace{\left\{\begin{array}{c}{\dot{\boldsymbol{u}}}_{\boldsymbol{m}\boldsymbol{r}}\\ {}{\ddot{\boldsymbol{u}}}_{\boldsymbol{m}\boldsymbol{r}}\end{array}\right\}}}=\underset{{\boldsymbol{A}}_{\boldsymbol{cm}}}{\underbrace{\left[\begin{array}{cc}\mathbf{0}& \mathbf{I}\\ {}-{\boldsymbol{\Omega}}_r^2& -{\boldsymbol{\Gamma}}_{\mathrm{r}}\end{array}\right]}}\underset{{\boldsymbol{x}}_{\boldsymbol{m}}(t)}{\underbrace{\left\{\begin{array}{c}{\boldsymbol{u}}_{\boldsymbol{m}\boldsymbol{r}}\\ {}{\dot{\boldsymbol{u}}}_{\boldsymbol{m}\boldsymbol{r}}\end{array}\right\}}}+\underset{{\boldsymbol{B}}_{\boldsymbol{cm}}}{\underbrace{\left[\begin{array}{c}\mathbf{0}\\ {}{\boldsymbol{\Phi}}_r^T{\mathbf{S}}_{\boldsymbol{p}}\end{array}\right]}}p(t),\end{array}} $$

with $ {\boldsymbol{x}}_{\boldsymbol{m}}\in {\mathrm{\mathbb{R}}}^{2{n}_m},\hskip0.35em {\mathbf{A}}_{\boldsymbol{cm}}\in {\mathrm{\mathbb{R}}}^{2{n}_m\times 2{n}_m},\hskip0.35em {\mathbf{B}}_{\boldsymbol{cm}}\in {\mathrm{\mathbb{R}}}^{2{n}_m\times {n}_p} $ .

2.2. Stochastic load model

The stochastic load model is defined as a stationary Gaussian process with zero mean and a Matérn covariance function with $ \nu =3/2 $ (see Rasmussen and Williams, 2006), hereafter referred to as the Matérn32 covariance function:

$$ p(t)\sim \mathcal{GP}\left(0,{\kappa}_{M32}\left(\tau; \sigma, {l}_{sc}\right)\right), $$

$$ {\kappa}_{M32}\left(\tau; \sigma, {l}_{sc}\right)={\sigma}^2\left\{1+\frac{\sqrt{3}\left|\tau \right|}{l_{sc}}\right\}{e}^{\left(-\frac{\sqrt{3}\left|\tau \right|}{l_{sc}}\right)}, $$

where $ \tau $ corresponds to the time lag $ \tau =\left|t-{t}^{\prime}\right| $ . The covariance function requires the definition of two hyperparameters: $ \sigma $ and $ {l}_{sc} $ , where $ \sigma $ controls the amplitude of the covariance function, and $ {l}_{sc} $ the extent of the correlation in time. It is highlighted that $ \sigma $ has force units (N), as it is related to the amplitude of the load, and $ {l}_{sc} $ has time units (s), since it is directly related to the expected frequency content of the load: a lower $ {l}_{sc} $ result in the presence of higher frequencies and vice versa. For further details, see Rasmussen and Williams (2006) and Hartikainen and Särkkä (2010).

A Gaussian process with a Matérn32 covariance function can be expressed as the output of the following linear time invariant stochastic differential equation (Hartikainen and Särkkä, 2010):

(2) $$ {\displaystyle \begin{array}{c}\ddot{p}(t)+\left(2\frac{\sqrt{3}}{l_{sc}}\right)\dot{p}(t)+\left(\frac{3}{l_{sc}^2}\right)p(t)=w(t).\end{array}} $$

Note that it is a stochastic differential equation as it is driven by a white Gaussian noise $ w(t) $ with a known spectral density

$$ w(t)\sim \mathcal{N}\left(0,{\sigma}_w^2\right);{\sigma}_w^2=\frac{12\sqrt{3}{\sigma}^2}{l_{sc}^3}. $$

Equation (2) can again be written in continuous-time state-space form as

(3) $$ {\displaystyle \begin{array}{c}\underset{\dot{\boldsymbol{z}}(t)}{\underbrace{\left[\begin{array}{c}\dot{p}(t)\\ {}\ddot{p}(t)\end{array}\right]}}=\underset{{\mathbf{F}}_{\boldsymbol{c}}}{\underbrace{\left[\begin{array}{cc}\hskip1em 0& \hskip2em 1\\ {}-\frac{3}{l_{sc}^2}& -2\frac{\sqrt{3}}{l_{sc}}\end{array}\right]}}\underset{\boldsymbol{z}(t)}{\underbrace{\left[\begin{array}{c}p(t)\\ {}\dot{p}(t)\end{array}\right]}}+\underset{{\mathbf{L}}_{\boldsymbol{c}}}{\underbrace{\left[\begin{array}{c}0\\ {}1\end{array}\right]}}w(t),\end{array}} $$

with $ \boldsymbol{z}(t)\in {\mathrm{\mathbb{R}}}^2 $ , $ {\mathbf{F}}_c\in {\mathrm{\mathbb{R}}}^{2x2} $ and $ {\mathbf{L}}_c\in {\mathrm{\mathbb{R}}}^2 $ . The load $ p(t) $ is directly related to the state vector $ z(t) $ by

$$ p(t)={\mathbf{H}}_cz(t), $$

with $ {\mathbf{H}}_c=\left[\begin{array}{cc}1& 0\end{array}\right]. $

Note that the performance of the GPLFM in virtual sensing depends on an appropriate choice of the GP kernel and its hyperparameters. A large number of convenient kernels are available, where a selection is made based on prior knowledge of the time-varying behavior of the modeled process, for example, periodicity, smoothness, or differentiability. The Matérn kernel chosen here is used to characterize processes that are finitely differentiable (Hartikainen and Särkkä, 2010), and is often employed in engineering applications because of its flexibility in accommodating a broad range of process traits.

2.3. Augmented model

Both state space models (equations 1 and 3) can be combined into an augmented system of equations as

(4) $$ {\displaystyle \begin{array}{c}\underset{{\dot{\boldsymbol{z}}}^{\boldsymbol{a}}(t)}{\underbrace{\left[\begin{array}{c}{\dot{\boldsymbol{x}}}_{\boldsymbol{m}}(t)\\ {}\dot{\boldsymbol{z}}(t)\end{array}\right]}}=\underset{{\mathbf{F}}_{\boldsymbol{c}}^{\boldsymbol{a}}}{\underbrace{\left[\begin{array}{cc}{\mathbf{A}}_{\boldsymbol{c}\boldsymbol{m}}& {\mathbf{B}}_{\boldsymbol{c}\boldsymbol{m}}{\mathbf{H}}_{\boldsymbol{c}}\\ {}0& {\mathbf{F}}_{\boldsymbol{c}}\end{array}\right]}}\underset{{\boldsymbol{z}}^{\boldsymbol{a}}(t)}{\underbrace{\left[\begin{array}{c}{\boldsymbol{x}}_{\boldsymbol{m}}(t)\\ {}\boldsymbol{z}(t)\end{array}\right]}}+\underset{{\boldsymbol{w}}^{\boldsymbol{a}}(t)}{\underbrace{\left[\begin{array}{c}\mathbf{0}\\ {}{\mathbf{L}}_{\boldsymbol{c}}w(t)\end{array}\right]}}\end{array}} $$

with

$$ {\boldsymbol{w}}^{\boldsymbol{a}}(t)\sim \mathcal{N}\left(0,{\mathbf{Q}}_{\boldsymbol{c}}^{\boldsymbol{a}}\right);{\mathbf{Q}}_{\boldsymbol{c}}^{\boldsymbol{a}}=\left[\begin{array}{cc}0& 0\\ {}0& {\boldsymbol{L}}_c{\sigma}_w^2{\boldsymbol{L}}_c^T\end{array}\right], $$

Where $ {\mathbf{Q}}_{\boldsymbol{c}}^{\boldsymbol{a}}\mathbf{\in}{\mathrm{\mathbb{R}}}^{2\times 2} $ is the augmented noise covariance matrix, $ {\boldsymbol{z}}^a(t)\in {\mathrm{\mathbb{R}}}^{2{n}_m+2} $ , and $ {\mathbf{F}}_c^a\in {\mathrm{\mathbb{R}}}^{\left(2{n}_m+2\right)\times \left(2{n}_m+2\right)} $ .

Note that the augmented state-space model is a stochastic differential equation. The augmented state has mean $ {\hat{\boldsymbol{z}}}^{\boldsymbol{a}}(t)=\unicode{x1D53C}\left[{\boldsymbol{z}}^{\boldsymbol{a}}(t)\right] $ and covariance $ {\mathbf{P}}^{\boldsymbol{a}}(t)=\unicode{x1D53C}\left[{\boldsymbol{z}}^{\boldsymbol{a}}(t){\boldsymbol{z}}^{\boldsymbol{a}}{(t)}^T\right] $ . Expressions that define the derivatives of these statistics can be found through equation (4). The details can be found from Särkkä and Solin (2019):

(5) $$ {\displaystyle \begin{array}{c}{\dot{\hat{\boldsymbol{z}}}}^a(t)={\mathbf{F}}_{\boldsymbol{c}}^{\boldsymbol{a}}{\hat{\boldsymbol{z}}}^a(t),\end{array}} $$

(6) $$ {\displaystyle \begin{array}{c}{\dot{\mathbf{P}}}^a(t)={\mathbf{F}}_{\boldsymbol{c}}^{\boldsymbol{a}}{\mathbf{P}}^{\boldsymbol{a}}(t)+{\mathbf{P}}^{\boldsymbol{a}}(t){{\mathbf{F}}_{\boldsymbol{c}}^{\boldsymbol{a}}}^T+{\mathbf{Q}}_{\boldsymbol{c}}^{\boldsymbol{a}},\end{array}} $$

with $ {\hat{\boldsymbol{z}}}^a(t)\in {\mathrm{\mathbb{R}}}^{2{n}_m+2} $ , $ {\mathbf{P}}^a(t)\in {\mathrm{\mathbb{R}}}^{\left(2{n}_m+2\right)\times \left(2{n}_m+2\right)} $ .

It is highlighted that when combining the (deterministic) mechanical model (equation 1) with the stochastic model for the load (equation 3) into the augmented model (equation 4), a stochastic definition for the state of the system is obtained (augmented state $ {\boldsymbol{z}}^a $ ). In particular, $ {\boldsymbol{z}}^a $ is now a Gaussian process: the response of the structure is recognised to be stochastic since it is driven by a stochastic load.

2.4. The prior distribution

The prior distribution of the augmented state $ \mathcal{p}\left({\boldsymbol{z}}^a(t)\right) $ can be obtained from the stationary solutions of the expressions for the statistics (equations 5 and 6, first derivatives equal to zero):

(7) $$ {\displaystyle \begin{array}{c}{\mathbf{F}}_{\boldsymbol{c}}^{\boldsymbol{a}}{\hat{\boldsymbol{z}}}_{\infty}^a=\mathbf{0},\end{array}} $$

(8) $$ {\displaystyle \begin{array}{c}{\mathbf{F}}_{\boldsymbol{c}}^a{\mathbf{P}}_{\infty}^a+{\mathbf{P}}_{\infty}^a{{\mathbf{F}}_{\boldsymbol{c}}^a}^T+{\mathbf{Q}}_{\boldsymbol{c}}^a=\mathbf{0},\end{array}} $$

Thus, as a result of equation (7), the stationary expected value is zero $ {\hat{\boldsymbol{z}}}_{\infty}^a=\boldsymbol{0} $ . The stationary covariance $ {\mathbf{P}}_{\infty}^a $ can be found numerically by solving the continuous-time Lyapunov equation, equation (8). The prior distribution is then defined as

$$ \mathcal{p}\left({\boldsymbol{z}}^a(t)\right)=\mathcal{N}\left(\boldsymbol{0},{\mathbf{P}}_{\infty}^a\right). $$

This is a prior for the augmented state which contains information regarding the response of the system as well as the load. Therefore, when defining a prior for the load, a prior for the response is also defined.

2.5. The posterior distribution

The posterior distribution of the augmented state refers to the conditional distribution upon some observations $ {\boldsymbol{y}}_{\left[1:N\right]} $ of the response

$$ \mathcal{p}\left({\boldsymbol{z}}_{\left[k\right]}^a|{\boldsymbol{y}}_{\left[1:N\right]}\right)=\mathcal{N}\left({\hat{\boldsymbol{z}}}_{\left[k\right]}^a,{\mathbf{P}}_{\left[k\right]}^a\right), $$

where discrete time notation is used. Kalman filtering and smoothing can be applied to efficiently compute the posterior distribution, as suggested by Hartikainen and Särkkä (2010). This requires the definition of the discrete version of the augmented model and the observation equation.

2.5.1. Discrete augmented model

The expressions found for the statistics (equations 5 and 6) can be solved, as these are ordinary differential equations. The solution can be expressed using discrete notation, evaluating between time steps, as

$$ {\hat{\boldsymbol{z}}}_{\left[k+1\right]}^a=\underset{{\mathbf{F}}^a}{\underbrace{{\boldsymbol{e}}^{{\boldsymbol{F}}_{\boldsymbol{c}}^a\Delta t}}}{\hat{\boldsymbol{z}}}_{\left[k\right]}^a, $$

$$ {\mathbf{P}}_{\left[k+1\right]}^a={\mathbf{F}}^a{\mathbf{P}}_{\left[k\right]}^a{{\mathbf{F}}^a}^T+\underset{{\mathbf{Q}}_{\boldsymbol{d}}^a}{\underbrace{\int_{k\Delta t}^{\left(k+1\right)\Delta t}{\boldsymbol{e}}^{{\mathbf{F}}_{\boldsymbol{c}}^{\boldsymbol{a}}\left(\left(k+1\right)\Delta t-\tau \right)}{\mathbf{Q}}_{\boldsymbol{c}}^{\boldsymbol{a}}{{\boldsymbol{e}}^{{\mathbf{F}}_{\boldsymbol{c}}^{\boldsymbol{a}}\left(\left(k+1\right)\Delta t-\tau \right)}}^T d\tau}}. $$

The expressions obtained for the statistics define an equivalent discrete version of the system

$$ {\boldsymbol{z}}_{\left[k+1\right]}^a={\mathbf{F}}^a{\boldsymbol{z}}_{\left[k\right]}^a+{\boldsymbol{w}}_{\left[k\right]}^a, $$

with $ {\boldsymbol{w}}_{\left[k\right]}^a\sim \mathcal{N}\left(0,{\mathbf{Q}}_{\boldsymbol{d}}^{\boldsymbol{a}}\right) $ . Defining the discrete version of the noise covariance $ {\mathbf{Q}}_{\boldsymbol{d}}^a $ through the integral is difficult and impractical. An alternative method of computing it is by realizing that the stationary solution is also valid for the discrete formulation.

$$ {\mathbf{P}}_{\infty}^a={\mathbf{F}}^a{\mathbf{P}}_{\infty}^a{{\mathbf{F}}^a}^T+{\mathbf{Q}}_{\boldsymbol{d}}^a $$

Since the stationary covariance can be computed by solving the continuous-time Lyapunov equation (equation 8), the noise covariance can then be obtained as

2.5.2. Observation equation

Let $ \boldsymbol{y}(t)\boldsymbol{\in}{\mathrm{\mathbb{R}}}^{n_y} $ be a subset of the system response,

$$ {\displaystyle \begin{array}{c}\boldsymbol{y}(t)={\mathbf{S}}_{\boldsymbol{a}}\ddot{\boldsymbol{u}}(t)+{\mathbf{S}}_{\boldsymbol{v}}\dot{\boldsymbol{u}}(t)+{\mathbf{S}}_{\boldsymbol{d}}\boldsymbol{u}(t).\end{array}} $$

The matrices $ {\mathbf{S}}_{\boldsymbol{a}},{\mathbf{S}}_{\boldsymbol{v}},{\mathbf{S}}_{\boldsymbol{d}}\in {\mathrm{\mathbb{R}}}^{n_y\times {n}_u} $ are the selection matrices for accelerations, velocities, and displacements, respectively. These relate the observed (measured) response $ \boldsymbol{y}(t) $ of the structure to the mechanical model.

Considering the modal formulation, and the equations of motion, the observations $ \boldsymbol{y}(t) $ can be defined in terms of the state and load applied to the system, see for example (Lourens, Papadimitriou, et al., 2012),

with dimensions: $ {\mathbf{G}}_{\boldsymbol{cm}}\in {\mathrm{\mathbb{R}}}^{n_y\times 2{n}_m};{\mathbf{J}}_{\boldsymbol{cm}}\in {\mathrm{\mathbb{R}}}^{n_y\times {n}_p} $ . Using the state space representation of the load, this can be written in terms of the augmented state

$$ \boldsymbol{y}(t)=\underset{{\mathbf{H}}_{\boldsymbol{c}}^{\boldsymbol{a}}}{\underbrace{\left[\begin{array}{cc}{\mathbf{G}}_{\boldsymbol{c}\boldsymbol{m}}& {\mathbf{J}}_{\boldsymbol{c}\boldsymbol{m}}{\mathbf{H}}_{\boldsymbol{c}}\end{array}\;\right]}}{\boldsymbol{z}}^{\boldsymbol{a}}(t). $$

Rewriting to discrete-time, and accounting for measurement noise, the following form of the observation equation is employed:

(9) $$ {\displaystyle \begin{array}{c}{\boldsymbol{y}}_{\left[k\right]}={\mathbf{H}}_{\boldsymbol{c}}^{\boldsymbol{a}}{\boldsymbol{z}}_{\left[k\right]}^a+{\boldsymbol{r}}_{\left[k\right]},\end{array}} $$

with $ {\boldsymbol{r}}_{\left[k\right]} $ a stochastic variable representing the measurement error $ {\boldsymbol{r}}_{\left[k\right]}\sim \mathcal{N}\left(0,\mathbf{R}\right) $ , where the noise matrix $ \mathbf{R}\in {\mathrm{\mathbb{R}}}^{n_y\times {n}_y} $ is assumed a diagonal matrix populated with the noise levels for each sensor $ {\sigma}_{R,s}^2 $ where “ $ s $ ” refers to sensor “s”.

2.5.3. Kalman filtering and smoothing

With the discrete version of the augmented model and the observation equation well defined, Kalman filtering and smoothing can be applied to recursively compute the posterior. The relevant expressions are summarized in Table 1.

Table 1. Kalman filter and smoother equations for the augmented model.

2.6. Hyperparameter estimation

So far, it has been assumed that the hyperparameters $ \sigma, \hskip0.35em {l}_{sc} $ of the Matérn32 covariance function that represents the load are well defined. However, given that the load is unknown, so is the covariance function that defines it stochastically. Furthermore, the noise levels $ {\sigma}_{R,s}^2 $ in $ \mathbf{R} $ are also unknown.

As presented in Särkkä (2013), the Bayesian manner to obtain these unknown parameters is either through maximum a posteriori or maximum likelihood estimation (a method used by Nayek et al., 2019), or by employing Markov Chain Monte Carlo methods (a method used by Rogers et al., 2020). An important disadvantage of these methods lies, however, in their computational cost. In what follows, an alternative, more practical method to compute the optimal values of the hyperparameters is suggested, where empirical qualities of the data are used to directly link the prior model to the measurements.

The suggested procedure for the hyperparameter estimation is divided into two parts: estimating the hyperparameters of the load ( $ \sigma, \hskip0.35em {l}_{sc} $ ), and estimating the noise levels of the sensors $ {\sigma}_{R,s}^2 $ .

2.6.1. Estimating the load hyperparameters (σ, l_sc)

In Section 2.4, it was shown that when a prior for the load is defined—through the hyperparameters $ \left(\sigma, {l}_{sc}\right) $ —then a prior for the response of the whole system is also defined. This includes the acceleration at the measured locations. The proposed method involves defining $ \left(\sigma, {l}_{sc}\right) $ so as to have a good match between the prior of the accelerations and the measurements.

The prior distribution of the accelerations can be easily computed from the prior distribution of the augmented state through the observation equation, disregarding the noise, so that

$$ {\boldsymbol{y}}_{\left(\mathrm{noisefree}\right)}(t)={\mathbf{H}}_{\boldsymbol{c}}^{\boldsymbol{a}}{\boldsymbol{z}}^{\boldsymbol{a}}(t), $$

and

$$ \mathcal{p}\left({\boldsymbol{y}}_{\left(\mathrm{noise}\hbox{-} \mathrm{free}\right)}(t)\right)=\mathcal{N}\left({\hat{\boldsymbol{y}}}_{\infty },{\mathbf{P}}_{y,\infty}\right), $$

$$ {\hat{\boldsymbol{y}}}_{\infty }={\mathbf{H}}_c^a{\hat{\boldsymbol{z}}}_{\infty}^a=\boldsymbol{0}, $$

$$ {\mathbf{P}}_{y,\infty }=\unicode{x1D53C}\left[{\boldsymbol{y}}_{\left(\mathrm{noise}\hbox{-} \mathrm{free}\right)}(t){\boldsymbol{y}}_{\left(\mathrm{noise}\hbox{-} \mathrm{free}\right)}{(t)}^T\right]={\mathbf{H}}_c^a{\mathbf{P}}_{\infty}^a{{\mathbf{H}}_c^a}^T, $$

where $ {\hat{\boldsymbol{y}}}_{\infty }=0 $ is the expected value and $ {\mathbf{P}}_{y,\infty } $ the associated covariance of the accelerations. In particular, the variance for the acceleration at each measured location “ $ s $ ”, $ {\sigma}_{y,s}^2 $ , will be defined by the diagonal elements of this covariance matrix $ {\mathbf{P}}_{y,\infty } $ . The correlation between the observed points is defined through the off-diagonal elements.

It is highlighted that, when calculating $ {\boldsymbol{y}}_{\left(\mathrm{noise}\hbox{-} \mathrm{free}\right)}(t)\sim \mathcal{N}\left({\hat{\boldsymbol{y}}}_{\infty },{\mathbf{P}}_{y,\infty}\right) $ through the augmented model, no information about the measurements has been used other than their location within the mechanical model. This is done through $ {\mathbf{H}}_c^a $ , which is constructed using the selection matrices for the measurement locations as well as the dynamic characteristics of the system (cf., Section 2.5.2.).

Since the load is modeled as a Gaussian process, the response obtained from the model is also a Gaussian process as shown above. Therefore, it is possible to compare the variance obtained from the model $ {\sigma}_{y,s}^2 $ with the variance estimated from the signal $ {\sigma_{y,s}^{\ast}}^2 $ . The best hyperparameters for the load $ \left(\sigma, {l}_{sc}\right) $ are such that $ {\sigma}_y^2\hskip0.35em \approx \hskip0.35em {\sigma_y^{\ast}}^2 $ . An exact match may, however, not be possible. This may for instance be due to (a) the model not being a perfect representation of the real structure, (b) an inaccurate estimation of the variance from the records, or (c) due to a high amount of noise in the signals.

To match the variances, an individual objective function is defined for each sensor “ $ s $ ”: $ {\mathcal{q}}_s=\mathcal{q}\left({\sigma}_{y,s}^2\right) $ . The individual objective function is defined as a lognormal distribution function. The peak will correspond to the computed variance of the measured records $ {\sigma_{y,s}^{\ast}}^2 $ , where $ {\mathcal{q}}_s^{\ast }=\mathcal{q}\left({\sigma_{y,s}^{\ast}}^2\right) $ , and the spread of the distribution is defined so as to have the 95% interval of confidence between half and twice the measured variance. The full definition and illustration can be found in Figure 1.

Figure 1. Definition of the objective function as a log-normal distribution.

When the hyperparameters $ \left(\sigma, {l}_{sc}\right) $ are defined, the variance of the accelerations at each sensor location $ {\sigma}_{y,s}^2 $ are defined as a consequence, from which each individual objective function can be evaluated: $ {\mathcal{q}}_s=\mathcal{q}\left({\sigma}_{y,s}^2\right) $ . The optimum hyperparameters $ \left(\sigma, {l}_{sc}\right) $ will be such that the individual objective function for each sensor is maximized: $ {\mathcal{q}}_s={\mathcal{q}}_s^{\ast } $ . Instead of evaluating each individual objective function, a combined objective function is defined by their multiplication: $ {\mathcal{q}}_{s1,\hskip0.35em s2,\hskip0.35em \dots, \hskip0.35em sn}={\mathcal{q}}_{s1}\cdotp {\mathcal{q}}_{s2}\cdotp \dots \cdotp {\mathcal{q}}_{sn} $ . Note that this combined objective function is the same as defining it as the joint (independent) lognormal distribution.

The best hyperparameters will therefore maximize this objective function $ {\mathcal{q}}_{s1,\hskip0.35em s2,\hskip0.35em \dots, \hskip0.35em sn} $ . The benefit of using this function is that it has some properties that are convenient for the analysis: it will take high values for variances relatively close to the estimated variance from the signal, and low values for variances that are far from either half or twice the estimated variance. Furthermore, each estimation can be assessed individually through each individual objective function. Alternative ways of defining the objective function or error measure between $ {\sigma}_{y,s}^2 $ and $ {\sigma_{y,s}^{\ast}}^2 $ may be adopted and or explored in future research.

2.6.2. Estimating the noise levels $ {\sigma}_{R,s}^2 $

The noise levels for each signal, $ {\sigma}_{R,s}^2 $ , contained in $ \mathbf{R} $ , also needs to be defined. This will be realized by studying the posterior distributions of the measured accelerations. The acceleration measurements are related to the model through the observation equation which includes noise (equation 9). An estimation for the noise process can therefore be obtained by using the expected value for the augmented state, obtained from the posterior

$$ {\hat{\boldsymbol{r}}}_{\left[k\right]}={\boldsymbol{y}}_{\left[k\right]}-{\mathbf{H}}_{\boldsymbol{c}}^{\boldsymbol{a}}{\hat{\boldsymbol{z}}}_{\left[k\right]}^a. $$

The measurements $ {\boldsymbol{y}}_{\left[k\right]} $ are deterministic, as the values are directly obtained from the accelerometers. An estimation of the variance of the noise can now be obtained as

$$ {\mathbf{R}}^{\ast }=\operatorname{var}\left({\hat{\boldsymbol{r}}}_{\left[k\right]}\right). $$

Assuming that the mechanical model (equation 1) and load definition (equation 3) are representative of the true structure and loading, and furthermore assuming that the sensor disturbances can be modeled as white Gaussian noise (equation 9), then the noise levels $ {\sigma}_{R,s}^2 $ defined for the system in $ \mathbf{R} $ should match with the noise levels observed from the posterior $ {\sigma_{R,s}^{\ast}}^2 $ in $ {\mathbf{R}}^{\ast }. $ Note that in order to estimate the noise levels $ {\mathbf{R}}^{\ast } $ a noise level must first be defined in $ \mathbf{R} $ .

The methodology employed to match the noise levels is iterative, as illustrated in Figure 2, and described as follows:

1. 1. An arbitrary amount of noise is defined for the signals. A starting point is taken equal to the variance of the signal $ \mathbf{R}=\operatorname{var}\left({\boldsymbol{y}}_{\left[k\right]}\right) $ .

2. 2. The posterior distribution is computed, from which an estimation for the noise is obtained: $ {\mathbf{R}}^{\ast } $ .

3. 3. The noise is defined as the estimated noise $ \mathbf{R}={\mathbf{R}}^{\ast} $ .

4. 4. Steps 2 and 3 are repeated until the maximum relative difference between the variances is less than a user-defined tolerance: $ \max \left(\operatorname{diag}\left\{\frac{\left|{\mathbf{R}}^{\ast }-\mathbf{R}\right|}{\mathbf{R}}\right\}\right)<\mathrm{tol} $ . If the relative difference between the variances starts to increase, the process is stopped and the achieved tolerance is reported.

Figure 2. Iterative process employed to match the noise levels $ {\sigma}_{R,s}^2 $ defined for the system in $ \boldsymbol{R} $ with the noise levels $ {\sigma}_{R,s}^2 $ observed in the system through $ {\boldsymbol{R}}^{\ast } $ .

3. Case Study: Onshore Wind Turbine

In what follows, the GPLFM is used for virtual sensing in an onshore wind turbine tower, with hyperparameter optimization performed as described in the previous section. First, specifics of the measurement campaign are introduced, followed by a presentation of the mechanical model employed, and its validation through system identification. Results are presented in Section 3.3.4 and onward.

3.1. Measurement campaign

Measurement data from one of EnBW’s modern onshore wind turbines were analyzed in this case study. It is a variable-speed, pitch-controlled, and three-bladed horizontal axis wind turbine located in Germany. The wind turbine is equipped with a dedicated structural health monitoring system. The sensor locations and local axes are illustrated in Figure 3. All sensors are mounted on the inner surface of the tower. The acceleration sensors are bidirectional, as shown in Figure 3. Four unidirectional strain gauges are mounted at the bottom of the tower; these are combined to obtain bending strains. All sensors measure at a sampling frequency of 20 Hz. The records considered are from December 9, 2020 to January 30, 2021. Both acceleration and strain measurements are divided into 10 min records for analysis. All records are filtered using a fourth-order Butterworth high-pass filter with a cut-off frequency of 0.1 Hz as the quasi-static information cannot be extracted from the acceleration measurements (quasi-static amplitudes are lower than the noise threshold).

Figure 3. Sensor locations and local axes.

Ten-minute statistics from the SCADA system are not used in the methodology itself, but for the interpretation of results obtained in different operational conditions in the following subsections.

Records are selected so as to cover the design situations for fatigue as defined in standard design codes (International Electrotechnical Commission, 2019). The design situations depend upon the cut-in and cut-out wind speed. For the wind turbine under consideration, the cut-in wind speed is $ {V}_{in}=3\;\mathrm{m}/\mathrm{s} $ and the cut-out wind speed $ {V}_{out}=25\;\mathrm{m}/\mathrm{s} $ . Table 2 summarizes the specific records considered in the upcoming analyses.

Table 2. Selected records.

3.2. Mechanical model and system identification

The mechanical model is based on a two-dimensional finite element model of the turbine, consisting of a set of beam elements for the tower and a point mass for the Rotor Nacelle Assembly. The foundation is assumed to clamp the tower at the bottom. The geometry of the tower varies in height. To capture this geometry, $ n=300 $ uniform beam elements are used. Each element differs in average diameter and thickness. Figure 4 summarizes the properties of the mechanical model. The exact same model is considered for both the FA and SS directions, implying that model accuracy will deteriorate for higher-order modes that are not dominated by tower deformation. The loading applied to the tower is modeled as a concentrated point load acting at the top of the tower. Note that the point load is defined stochastically as per equation (3), with its hyperparameters yet to be defined. The damping matrix is defined by assuming a damping ratio of 0.75% for all modes. This damping value is chosen on the basis of the system identification results discussed below.

Figure 4. Details of the mechanical model.

To obtain an estimate of the accuracy of the mechanical model, system identification is employed. Covariance-driven stochastic subspace identification (SSI-Cov, Peeters and De Roeck, 1999) is applied, in combination with the OPTICS clustering algorithm (OPTICS, Ankerst et al., 1999). The approach employed is the one presented in Boroschek and Bilbao (2019). The system identification is applied to distinguish between the tower Fore–Aft (FA) and tower Side–Side (SS) direction of the wind turbine, which are obtained by transforming the local coordinates using the yaw angle provided by the SCADA data. The identified modal properties and their comparison to the mechanical model are shown in Table 3. The three tower modes are well observed by the data. Given the good correspondence with the mechanical model, changes in the mechanical model to better fit the identification results were deemed unnecessary.

Table 3. Frequency ( $ f $ ), damping ratio ( $ \zeta $ ), and Modal Assurance Criterion (MAC) values between mechanical model and identification results based on record PK1.

3.3. Strain and fatigue estimation on the onshore turbine tower

The GPLFM is now applied on the mechanical model, informed by each set of acceleration records. The estimated states and forces are used to obtain the strain at the location where the strain gauges are located, to compare the estimations with the measurements.

In what follows, the selected records are first illustrated and described, and subsequently, the hyperparameter estimation procedure is shown for both the prior distribution (to obtain the hyperparameters $ \sigma $ and $ {l}_{sc} $ ) and the posterior distribution (to obtain the noise levels $ {\sigma}_{R,s}^2 $ in $ \mathbf{R} $ ).

3.3.1. Description of selected records

The power production records (OP1 to OP7) are characterized by relatively stationary behavior of the response, illustrated in Figure 5a,b. As a general trend, operational records associated with lower wind speeds will have lower amplitudes of vibrations, and vice versa.

Figure 5. Acceleration records, sensor s1. (a) OP1-x (mean wind speed 4 m/s). (b) OP7-x (mean wind speed 16 m/s). (c) PK2-x (quantization issues due to lack of resolution in the time signal). (d) SU-3 start-up condition. (e) SD2-x (shut-down condition).

Idling records (PK1–PK3) are also stationary but are characterized by very low amplitudes of vibrations given the low wind speeds that the structure is subject to. These records were observed to have quantization issues, illustrated in Figure 5c.

Start-up and shut-down records (SU1–SU3 and SD1–SD3, respectively) are characterized by transient behavior due to the operational change—compare Figure 5d and 5e.

3.3.2. Estimation of the load hyperparameters $ \left(\boldsymbol{\sigma}, {\boldsymbol{l}}_{\boldsymbol{sc}}\right) $

The hyperparameters $ \sigma $ and $ {l}_{sc} $ are defined for each analyzed record. This is done through the objective function described in Section 2.6.1. The objective function for record PK1-x is shown in Figure 6, where a clear maximum is observed, which corresponds to the values used for the hyperparameters. Figure 7 compares the prior defined through the chosen hyperparameters $ \left(\sigma, {l}_{sc}\right) $ with the distribution observed from the acceleration measurements. Because of the good fit observed, the assumption for the load location is validated.

Figure 6. Record PK1-x. $ \sigma $ and $ {l}_{sc} $ estimation. Objective function plot: $ {\mathcal{q}}_{s1,s2,s3} $ . Optimum found for maximum value: $ {l}_{sc}=0.318\;\mathrm{s} $ , $ \sigma =2556\;\mathrm{N} $ .

Figure 7. Record PK1-x. Prior distribution visual check using optimum hyperparameters: $ {l}_{sc}=0.318\;\mathrm{s} $ , $ \sigma =2556\;\mathrm{N} $ . (a) Sensor s1. (b) Sensor s2. (c) Sensor s3.

The comparison of the prior with the distribution observed from the acceleration measurements for record OP6-x is shown in Figure 8. Note that in this case, the fit is not as good as in the idling case. This is justified given the change in the operational condition of the turbine, which implies that the mechanical model—compared only in idling conditions—might not be as representative.

Figure 8. Record OP6-x. Prior distribution visual check using optimum hyperparameters: $ {l}_{sc}=0.05\;\mathrm{s} $ , $ \sigma =20845\;\mathrm{N} $ . (a) Sensor s1. (b) Sensor s2. (c) Sensor s3.

For transient conditions (shut-down and start-up records) the acceleration signal is not stationary. Therefore, the prior distribution will over- and/or underestimate the acceleration response, as illustrated in Figure 9 for record SU3-x.

Figure 9. Record SU3-x. Prior distribution visual check using optimum hyperparameters: $ {l}_{sc}=0.034\;\mathrm{s} $ , $ \sigma =866\;\mathrm{N} $ . (a) Sensor s1. (b) Sensor s2. (c) Sensor s3.

Tables 4 and 5 summarizes the objective function values obtained for each record in local x and y directions, respectively. They include, besides the values defined for $ \sigma $ and $ {l}_{sc} $ , the likelihood measured as a fraction of the maximum likelihood defined for each lognormal distribution and for the joint distribution. Therefore, a value close to 100% indicates a very good fit for the prior distribution.

Table 4. Estimation of load hyperparameters $ \left(\sigma, {l}_{sc}\right) $ .

Note. Prior distribution match measured through objective function $ \mathcal{q} $ . Local direction x.

Table 5. Estimation of load hyperparameters $ \left(\sigma, {l}_{sc}\right) $ .

Note. Prior distribution match measured through objective function $ \mathcal{q} $ . Local direction y.

It is observed that for parking conditions (PK1–PK3), the relative likelihood is high, indicating that the mechanical model provides a good representation of the structure in idling conditions. For all power production conditions (OP1–OP7) the relative likelihood is lower, which indicates that the issue previously discussed for record OP3-x is valid for all other power production records: the mechanical model considered is less representative under operating conditions (possibly due to aerodynamic damping, increased relevance of higher modes, etc.). Finally, for start-up and shut-down records (SU1–SU3 and SD1–SD3, respectively) the relative likelihood varies. For these cases, and as previously mentioned for record SU3-x, the response is not stationary. In the future, one could investigate redefining the load assumption for these particular cases.

3.3.3. Estimating the noise levels $ {\boldsymbol{\sigma}}_{\boldsymbol{R},\boldsymbol{s}}^{\boldsymbol{2}} $

The assumed noise levels for each signal are defined iteratively, as described in Section 2.6.2, with a relative tolerance of 1%. This iterative process was observed to diverge in approximately 30% of the records. The process is illustrated using records PK1-x, OP3-x, and SU3-x.

Figure 10a shows the noise levels obtained in each iteration of record PK1-x. Note that the starting value is the variance of the signal, and it quickly decays into a constant value that corresponds to the estimated variance of the noise. This is further illustrated by the relative error obtained in each iteration as shown in Figure 10b. Figure 11 shows a visual check of the resulting noise variance compared with the estimated noise $ \hat{\boldsymbol{r}} $ .

Figure 10. Record PK1-x. Noise variance estimation $ {\sigma}_{R,s}^2 $ . (a) Noise variance estimation per iteration. (b) Error measure: $ \mathit{\operatorname{diag}}\left(\frac{\hat{\boldsymbol{R}}-\boldsymbol{R}}{\boldsymbol{R}}\right) $ .

Figure 11. Record PK1-x. Noise variance visual check using optimized noise variances.

Figure 12. Record OP3-x. Noise variance estimation $ {\sigma}_{R,s}^2 $ . (a) Noise variance estimation per iteration. (b) Error measure: $ \mathit{\operatorname{diag}}\left(\frac{\hat{\boldsymbol{R}}-\boldsymbol{R}}{\boldsymbol{R}}\right) $ .

Figure 13a shows the noise levels defined for each iteration of record OP3-x. Similar to record PK1-x, the noise levels are observed to decrease to a stable value. Figure 13b shows the relative error found for each iteration. Note that in this case, the relative error for accelerometer s3 increases in the third iteration. The iteration is not stopped, however, as the maximum error at each iteration is still decreasing. It should be noted that in some cases the part of the signal considered as noise may in fact contain a part of the response, complicating the search for an optimal value.

Figure 13. Record SU3-x. Noise variance estimation $ {\sigma}_{R,s}^2 $ . (a) Noise variance estimation per iteration. (b) Error measure: $ \mathit{\operatorname{diag}}\left(\frac{\hat{\boldsymbol{R}}-\boldsymbol{R}}{\boldsymbol{R}}\right) $ .

Figure 14a shows the noise levels per iteration of record SU3-x. Note that this record was diverging given the increase in the relative error, as illustrated in Figure 14b. Given the transient nature of the response, the predicted error is not Gaussian, and as discussed in Section 3.3.2, the model employed is not representative for the particular record. The resulting high error variance and lack of convergence attest to inaccurate modeling of the measurements in this case.

Figure 14 Record PK1-x. Dynamic strain estimation.

A summary is presented in Tables 6 and 7. The noise level obtained for each signal is presented in terms of the noise-to-signal ratio (NSR), defined as follows:

$$ NSR(s)=\frac{\sigma_{R,s}^2}{\sigma_{y,s}^2-{\sigma}_{R,s}^2}, $$

with $ {\sigma}_{y,s}^2 $ referring to the measured acceleration variance, and the difference $ \left({\sigma}_{y,s}^2-{\sigma}_{R,s}^2\right) $ corresponding to the noiseless signal variance. The $ NSR $ thus indicates the estimated relative amount of noise with respect to the noiseless response. If the $ NSR $ is larger than 100%, then the variance of the noise is higher than the variance of the noiseless signal, and vice versa. For idling conditions (PK1 to PK3), the process never failed to converge. This can be understood given a representative mechanical model for the structure. The found noise levels are quite high for these records, however, which can be justified by the considerable amount of quantization issues. For power production records (OP1–OP7), it is observed that for most cases the process converges. It is debatable, however, if the defined noise levels are representative of the true measurement error as the mechanical model is less accurate. The nonconvergence of some of these records can be justified by the inaccuracy of the mechanical model employed, as previously discussed. Finally, for start-up and shut-down records (SU1–SU3 and SD1–SD3, respectively) the response is transient, and therefore the estimation of the error is questionable, as the model employed was shown to be less representative. Nevertheless, as will be shown in the following sections, the procedure provides a computationally inexpensive way to obtain hyperparameter values that allow for accurate estimation results.

Table 6. Noise variance estimation.

Note. Noise to signal ratios. Local direction x.

Figure 15. Record OP3-x. Dynamic strain estimation.

Table 7. Noise variance estimation.

Note. Noise-to-signal ratios. Local direction y.

3.3.4. Response estimation

Having analyzed and defined the hyperparameters for each record, the estimated response in terms of the dynamic strains is now compared with the measured dynamics strains.

Figure 16a shows the dynamic strain estimation for record PK1-x, along with the measured and filtered strain for comparison. The mean refers to the expected value obtained for each point in time through the posterior distribution for the strain, and $ 3\sigma $ is the associated 99.7% interval of confidence. Figure 16b shows the associated frequency contents of the measurement and estimation. The unfiltered measurements are also included in the spectra, to show the influence band of the filtering. The dynamic strain estimation is observed to match the measured strain, with the mean prediction closely following the measured strain. The differences in the frequency domain are mainly observed in the low amplitude and high-frequency regions; the estimation for the dominant frequencies is observed to be accurate.

Figure 16. Record SU3-x. Dynamic strain estimation.

Figure 11a,b presents the estimation of record OP3-x in the time and frequency domain, respectively. The estimation is again observed to be a close match, though to a lesser degree when compared with the estimation obtained for record PK1-x. This is justified by the analysis exposed in Sections 3.3.2 and 3.3.3: the model employed is less accurate for this operational condition.

Figure 12a,b shows the estimation of record SU3-x in the time and frequency domain, respectively. As expected, larger errors are observed given the transient nature of the response, but the estimations are still reasonably accurate outside these transient parts. It is noted that even though the transient nature is not part of the load assumptions, it is still well estimated in the load, as shown in Figure 17.

Figure 17. Record SU3-x. Load estimation.

To quantify the estimation errors, two error metrics are considered: The Time Response Assurance Criterion (TRAC) and the Mean Absolute Error (MAE). These are defined as follows:

$$ TRAC\left(\boldsymbol{\varepsilon}, \hat{\boldsymbol{\varepsilon}}\right)=\frac{{\left(\hat{\boldsymbol{\varepsilon}}{\boldsymbol{\varepsilon}}^T\right)}^2}{\left(\hat{\boldsymbol{\varepsilon}}{\hat{\boldsymbol{\varepsilon}}}^T\right)\left(\boldsymbol{\varepsilon} {\boldsymbol{\varepsilon}}^T\right)}, $$

$$ MAE\left(\boldsymbol{\varepsilon}, \hat{\boldsymbol{\varepsilon}}\right)=\frac{1}{N}{\sum}_k\left|{\hat{\varepsilon}}_{\left[k\right]}-{\varepsilon}_{\left[k\right]}\right|, $$

with $ \boldsymbol{\varepsilon} \in {\mathrm{\mathbb{R}}}^N $ the measured and filtered strain in vector form representing the values in time, and $ \hat{\boldsymbol{\varepsilon}}\in {\mathrm{\mathbb{R}}}^N $ the estimation. The TRAC represents a measure of correlation between the estimation and the measurement. In simple words, a value close to 100% indicates a very strong similarity of the shapes, and a value close to 0% indicates the estimations and measurements are very different. The MAE indicates the average error of the estimation, in units of the measurement. The results for each analyzed record are summarized in Table 8.

Table 8. Error metrics.

3.3.5. Fatigue load estimation

The underlying aim of this work is to be able to predict the damage that the tower has suffered in terms of fatigue. In this section, it is shown how the accuracy obtained for the strain estimates translates to accuracy in fatigue damage predictions. Damage equivalent loads are used to compare design fatigue load levels with fatigue loads the asset experienced on site. The damage equivalent load, in terms of stresses, is defined by, see for example (Cosack, 2010):

$$ \Delta {\sigma}_{eq}=\sqrt[m]{\frac{\sum_i\Delta {\sigma}_i^m{N}_i}{N_{ref}}}, $$

with $ \Delta {\sigma}_i $ stress ranges and $ {N}_i $ the corresponding number of cycles. $ m $ is a material constant referring to the idealized straight slope observed on a double-logarithmic scaled SN-curve. $ \Delta {\sigma}_{eq} $ is the damage equivalent load, and $ {N}_{ref} $ the number of cycles of the equivalent load.

The stress ranges and the corresponding number of cycles are determined using the rainflow-counting method. The rainflow-counting method is applied to each record, and subsequently combined to generate the stress ranges and the corresponding number of cycles for all the records analyzed. This was done by counting the number of cycles for each record, and then adding them together. The result is shown in Figure 18. This figure also includes the rainflow-counting results for the unfiltered strain signals showcasing the general lack of the dynamic strains to capture the quasi-static high-amplitude stress ranges. The damage equivalent load and relative errors are presented in Table 9. A more detailed analysis of the damage equivalent load derived for each record is provided in Table 10. Note that stresses are computed by multiplying the strains with the assumed modulus of elasticity $ E=\mathrm{200,000}\;\mathrm{MPa} $ , the material constant is set as $ m=4 $ and the reference number of cycles is set as $ {N}_{ref}={10}^7 $ . The Wöhler slope m = 4 is commonly used for welded steel as simplification of its SN-curve with two slopes of respectively 3 (applies to high stress / low number of cycles) and 5 (applies to low stress/high number of cycle).

Figure 18. Stress ranges and number of cycles. Rainflow-counting method.

Table 9. Damage equivalent load.

Note. All records summary.

Table 10. Damage equivalent load.

Note. All records.

Table 11. Damage equivalent load.

Note. All records. y-direction.

From Table 10 and 11, it is clear that the estimation error in terms of damage equivalent loads varies greatly across the different records, from less than 1% to almost 30%. Averaging over the records, errors of less than 10% are found—compare Table 9.

4. Conclusions

The GPLFM was used to estimate the dynamic strains and subsequently predict damage equivalent loads in an onshore wind turbine tower, based only on a set of acceleration measurements. Comparisons were made with measured and filtered strains, and the damage equivalent loads were computed from them. To find the hyperparameters for the GPLFM, a new, intuitive procedure was suggested. It is highlighted that quasi-static loads and strains are not estimated, since quasi-static information cannot be extracted from the acceleration measurements. Since these loads should, however, be considered when assessing the life of a turbine, it is suggested to include quasi-static information through inclinometer measurements in the future.

The results show that the GPLFM in combination with the suggested hyperparameter optimization procedure provides a robust Bayesian framework for response estimation. The dynamic strains were estimated with, on average, high accuracy, with slight differences in accuracy depending on the operating condition. The highest accuracy is achieved for the parked plant due to a good agreement between model assumptions and the operating condition. The shut-down and start-up operating states, on the other hand, are captured with the lowest precision owing to the transient nature of these operating states. The average TRAC value per operational condition was found to be 90.8, 92.0, 98.0, and 92.3 for the power production, start-up, shut-down, and idling conditions, respectively.

The dynamic strains estimated through GPLFM were further utilized to assess the damage suffered by the tower in terms of damage equivalent loads. Contrary to the results of the dynamic strain analysis, the closest match between measured and GPLFM-based damage equivalent loads is observed for the turbine during standard operation.

Factors that were shown to influence the accuracy of the estimation include the accuracy of the mechanical model employed, and the validity of the assumptions made when formulating the stochastic model representing the load. Future research should aim to improve these factors.