# Damage detection and localization method for wind turbine rotor based on Operational Modal Analysis and anisotropy tracking

Ambroise Cadoret, Enora Denimal-Goy, Jean-Marc Leroy, Jean-Lou Pfister, Laurent Mevel

► To cite this version:

Ambroise Cadoret, Enora Denimal-Goy, Jean-Marc Leroy, Jean-Lou Pfister, Laurent Mevel. Damage detection and localization method for wind turbine rotor based on Operational Modal Analysis and anisotropy tracking. Mechanical Systems and Signal Processing, 2025, 224, pp.111982. 10.1016/j.ymssp.2024.111982. hal-04902642

HAL Id: hal-04902642

https://inria.hal.science/hal-04902642v1

Submitted on 21 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Damage detection and localization method for wind turbine rotor based on Operational Modal Analysis and anisotropy tracking

Ambroise Cadoret<sup>a,b</sup>, Enora Denimal Goy<sup>b</sup>, Jean-Marc Leroy<sup>a</sup>, Jean-Lou Pfister<sup>c</sup>, Laurent Mevel<sup>b,*</sup>

<sup>a</sup>IFP Energies nouvelles, 69360 Solaize, France

<sup>b</sup>Université Gustave Eiffel, Inria, COSYS-SII, I4S, 35042 Rennes, France

<sup>c</sup>IFP Energies nouvelles, 92852 Rueil-Malmaison, France

## Abstract

Subspace-based damage detection methods are widely used for civil engineering structures modeled as linear time-invariant systems. For operating wind turbines modeled as linear time-periodic systems, these methods cannot be theoretically used, due to the inherent assumptions associated with these methods in the context of linear time-invariant systems. Based on a model approximation of time-periodic systems as time-invariant ones, these methods can still be applied and adapted to perform change detection for time-periodic systems, through a Gaussian residual built upon damage sensitive parameters coming from the identified modal parameters. The proposed method is tested and validated on data simulated with an aero-servo-elastic model of an operating wind turbine, with the detection of local stiffness reductions on a blade and the localization of the damaged blade. Furthermore, two selections of sensors are tested, to evaluate the impact of the sensor choice on the performance of the detection and localization methods.

**Keywords:** Damage detection, Wind turbine, Vibration measurements, Statistical tests, Operational Modal Analysis

## 1. Introduction

Regarding the actual installed capacity, given the prospects for growth in the number of wind farms in the coming years, and the continuing quest to reduce maintenance costs, it is important to implement reliable monitoring methods based on data collected during the operation of wind turbines. Early fault detection is a requirement in order to increase the lifetime of wind turbines and to improve the planning of maintenance operations. In [2], a study on the failure rates of each component of a wind turbine has shown that the rotor blades are one of the critical components to study. It is therefore important to implement reliable methods for damage detection for wind turbine blades.

Currently, there are several damage detection methods designed for wind turbines. The first methods are focused on detecting aerodynamic or mass imbalance of the rotor blades through tower sensors. In [26], the imbalance is detected with the PSD of the electrical power, and then the nature and the severity are determined using blade measurements. In [8], the focus is on pitch errors, based on tower measurements, where the detection localization and quantification are performed using a multi-layer perceptron. In [22], the PSD of the rotational speed is used as the damage feature and the support vector machine is chosen to perform detection and quantification of mass imbalances. All the cited works are performed on simulated data. However, in [47] the impact of

<footnote>*Corresponding author; E-mail address: laurent.mevel@inria.fr</footnote>
blade icing on the PSD of tower accelerations is shown on real measurements. In the literature, some methods defined for the detection of structural rotor faults can be found for example in [44, 43], with the detection of trailing edge opening on one blade of a small wind turbine (225 kW). In these papers, the cross-covariance of rotor sensors is used as a damage feature through a Mahalanobis distance. The issue with this method is its sensitivity to the operational conditions that must be mitigated, as it is done in [39, 32, 20]. However, this approach does not provide localization information. Finally, in all the cited papers, a mechanical actuator is installed on the monitored blade. This results in a signal with low covariance and constant external conditions, which is not relevant to the real external conditions of a wind turbine blade.

The well-established methods based on subspace identification for Linear Time Invariant (LTI) modeled structures have been widely used for civil structures. These methods must be challenged to be applied on the operating wind turbine rotors, which are modeled as Linear Time Periodic (LTP) systems.

For civil engineering structures (modeled as LTI systems), a damage detection method based on a local approach for change detection can be used as in [17], and its extension to deal with real structures has been demonstrated in [14]. With this method, based on a model parametrization, it is possible to perform change isolation, corresponding to damage localization in the civil or mechanical engineering context [5]. Moreover, the possibility to isolate any change in the parameterization can be assessed [27]. For detection and localization, the approach performs a statistical test on a data-driven residual. To be robust to environmental conditions and built upon the identified estimates of the system modes, some recent residual was directly designed in [30], using the mode shapes identified with a subspace identification. Also, the definition of the residual has been extended with a residual function of an estimated reference [46], which is important when a reference state has to be learned from the data and not from the model. In the current paper, any proposed residual should be designed with these two considerations in mind.

There are many system identification methods to deal with periodic systems, like the Periodic PO-MOESP [38], the partial Floquet analysis [40] or the periodic SSI [23], which identify exactly the periodic modes of the structures. A priori, the previously presented *local* damage detection and localization approach cannot be used for monitoring operating wind turbines (modeled as LTP systems), since this approach relies on the same theoretical assumptions than the standard subspace identification methods developed for LTI models. [10, 13] has nonetheless established that LTP systems can be approximated as LTI systems with periodic feed-through matrices. Consequently, these systems can be identified with the standard identification methods (SSI) based on Operational Modal Analysis (OMA), and the damage detection and localization method designed for LTI systems based on identification methods can be performed thanks to this model approximation. It has in particular been shown that the modes identified by the LTI subspace identification methods are the Fourier components associated with the Fourier decomposition of the periodic modes corresponding to the underlying LTP model of the wind turbine.

However, once the usage of an LTI subspace-based damage detection method for wind turbine monitoring is accepted, the exact definition of the residual should be investigated. Indeed, it has been shown that the subspace-based damage detection method defined in [15] belongs to a general class of Gaussian residual-based damage detection methods that can be redefined depending on the monitored features of interest, as in [17].

To adapt this method to the detection of faults on wind turbine blades, the damage detection methods need
to be defined as a function of the relevant damage-sensitive features of the studied system. In [18], the type of damage that occurs on a wind turbine rotor and damage detection methods have been explored. It has been shown that to monitor wind turbine rotor faults, a solution is to monitor the anisotropy or the imbalance (mass or pitch anisotropy) - see [9] for more discussion on this topic. In [42], it was stated that the phase shift and amplitudes of the edge bending modes shape are the most sensitive damage-sensitive features for rotor anisotropy detection. Therefore, these features will be used to construct a new residual for wind turbine rotor damage detection.

The main contribution of this paper is the adaptation and the application of the state-of-the-art damage detection method to operating wind turbine rotors. Furthermore, a mean residual is defined to increase the damage detectability while using fixed-length data sets.

The paper is organized as follows. Section 2 defines the OMA procedure used for the identification of the operating wind turbine rotor, which is based on the approximation of the Floquet modes as a finite sum of their non-negligible Fourier harmonics. Then, Section 3 presents the damage detection and localization methods and the adaptation for the detection and localization of wind turbine rotor faults. The damage detection approach is based on the development of a residual and the associated statistical theory, based on the aggregation of identified features averaged from multiple measurement data sets. The performance of the residual is theoretically investigated by computing the probability of detection (POD) from the Fisher information matrix and the corresponding detectability information. Finally, Section 4 applies the methods, with the detection of local stiffness losses of a blade of the DTU 10MW wind turbine model [36], and the localization of the damaged blade. Also, two selections of sensors are tested, to evaluate the impact on damage localization.

## 2. Operational Modal analysis of rotating wind turbine rotor

In this section, the same assumption as [13] are used, with more details and discussions in the aforementioned paper.

### 2.1. *Dynamic model*

The motion of a wind turbine operating with a constant rotating speed can be expressed after linearization as an LTP-system,

$$ \mathcal{M}(t)\ddot{u}(t) + \mathcal{C}(t)\dot{u}(t) + \mathcal{K}(t)u(t) = v(t), \quad (1) $$

where $u(t) \in \mathbb{R}^m$ represents the displacements of the structure at the degrees of freedom (DOF) of the system, and $\mathcal{M}(t+T) = \mathcal{M}(t)$, $\mathcal{C}(t+T) = \mathcal{C}(t)$, $\mathcal{K}(t+T) = \mathcal{K}(t)$, respectively the mass, damping and stiffness matrices, with $T$ represents the rotational period. Notice that the state variable is an unmeasured unknown variable. As such, the modelling can encompass any dynamical state component that can be related to the observations. The unknown input $v(t)$ is assumed to be a Gaussian white noise. Refer to [13] for more discussions on the fact that the input should be better modeled as separated in deterministic (pitch, torque, yaw bearing angle and average wind at the hub) and stochastic (turbulence). In the following, the mechanical system is expressed in state-space form, from the definition of the state vector $x(t) \in \mathbb{R}^n$ where $n = 2m$ and the observation vector $y(t) \in \mathbb{R}^r$.

$$ x(t) = \begin{bmatrix} u(t) \\ \dot{u}(t) \end{bmatrix} \quad \text{and} \quad y(t) = C_a\ddot{u}(t) + C_v\dot{u}(t) + C_du(t), \quad (2) $$
where $C_a, C_v$ and $C_d$ are the selection matrices. A noise $w(t)$ can be added to the observation. $w(t)$ is assumed to be a Gaussian white noise, that models the small random fluctuations of the aerodynamic forces due to atmospheric turbulence. This leads to the following state-space expression:

$$ \begin{cases} \dot{x}(t) = A_c(t)x(t) + B_c(t)v(t) \\ y(t) = C(t)x(t) + D(t)v(t) + w(t) \end{cases}, \qquad (3) $$

with

$$ \begin{aligned} A_c(t) &= \begin{bmatrix} 0 & I \\ -\mathcal{M}(t)^{-1}\mathcal{K}(t) & -\mathcal{M}(t)^{-1}\mathcal{C}(t) \end{bmatrix}, \quad C(t) = \begin{bmatrix} C_d - C_a\mathcal{M}(t)^{-1}\mathcal{K}(t) & C_v - C_a\mathcal{M}(t)^{-1}\mathcal{C}(t) \end{bmatrix}, \\ B_c(t) &= \begin{bmatrix} 0 \\ -\mathcal{M}(t)^{-1} \end{bmatrix} \quad \text{and} \quad D(t) = C_a\mathcal{M}(t)^{-1}(t). \end{aligned} $$

All the matrices are periodic with period $T$, with $A_c(t) \in \mathbb{R}^{n \times n}$, $C(t) \in \mathbb{R}^{r \times n}$, $B_c(t) \in \mathbb{R}^{n \times m}$ and $D(t) \in \mathbb{R}^{r \times m}$.

Using Floquet theory [19, 41], the homogeneous part of the state vector can be expressed as a sum of $n$ eigenmodes with a periodic amplitude, namely the Floquet modes

$$ x_h(t) = \sum_{j=1}^{n} X_j(t) \exp(\mu_j t) q_j(t_0), \qquad (4) $$

with $\mu_j$ being the $j$-th characteristic exponent, $q_j(t_0)$ a scalar depending of the initial conditions and $X_j(t)$ the associated $T$-periodic amplitude.

## 2.2. Approximation of the Floquet modes

Once the Floquet mode decomposition is applied, the observation vector is expressed as a finite sum of eigenmodes to obtain the description of a time invariant system. Similarly, the Floquet mode decomposition of the observation $y_h(t)$ is

$$ y_h(t) = \sum_{j=1}^{n} Y_j(t) \exp(\mu_j t) q_j(t_0), \qquad (5) $$

where $Y_j(t) = C(t)X_j(t)$ is both the amplitude of the $j$-th Floquet mode of the observation and a periodic vector of period $T = \frac{2\pi}{\Omega}$ ($\Omega$ the rotational speed). Using Equation (5) and expanding into a Fourier series the amplitudes of the Floquet modes, it comes

$$ y_h(t) = \sum_{j=1}^{n} \sum_{l=-\infty}^{\infty} Y_{j,l} \exp((\mu_j + il\Omega)t) q_j(t_0). \qquad (6) $$

with $Y_{j,l}$ the Fourier coefficients. The contributions of the harmonics in the expansion of $y_h(t)$ are determined by the participation factor [7, 37] and its variant as described in [13]

$$ \phi_{j,l}^{y} = \frac{\|Y_{j,l}\|^2}{\sum_{l=-\infty}^{\infty} \|Y_{j,l}\|^2}. \qquad (7) $$

A minimum participation factor ($\phi_{min}^y$) is applied to select harmonics with the highest contributions. Then, an approximation of the observation is constructed as $\hat{y}(t)$ by a finite sum of eigenmodes. Notice that this participation factor is only introduced to explain how the periodic model transitions to a linear time invariant
model. It is not used in the system identification method, which is the standard SSI method for LTI systems. Refer to [13] for a complete discussion.

$$ \hat{y}_h(t) = \sum_{(j,l), \phi_{j,l}^y \ge \phi_{min}^y} Y_{j,l} \exp((\mu_j + il\Omega)t) q_j(t_0), \quad (8) $$

then, $\hat{y}_h(t)$ can then be expressed as a sum of $\tilde{n}$ eigenmodes

$$ \hat{y}_h(t) = \sum_{p=1}^{\tilde{n}} Y_p \exp(\mu_p t) q_p(t_0), \quad (9) $$

where each index $p$ corresponds to a pair $(j, l)$ and $\mu_p = \mu_j + il\Omega$. This homogeneous part of the approximation is similar to the homogeneous part of a linear time-invariant system.

Finally, the Floquet modes of an LTP system have been approximated by a finite number of terms identical to those of an LTI system, for the free response of this LTP system. Subsequently, the forced response can also be extended to such formalism as detailed in [13], where the reader is referred for further details.

## 2.3. Subspace Identification

In [13], it has been established that a LTP system can be approximated and discretized as a linear system with constant **A** and **C** matrices, and periodic **B** and **D** matrices. For conciseness, refer to [13] for the details, and assume that from the definition of an LTP system and its associated approximation, it is possible to define the following discrete state-space model

$$ \left\{ \begin{array}{l} z_{k+1} = \mathbf{A}z_k + \mathbf{B}_k v_k \\ y_k = \mathbf{C}z_k + \mathbf{D}_k v_k + \tilde{w}_k \end{array} \right. , \quad (10) $$

with $z_k$ the new state vector based on the approximation of the outputs, **A** and **C** two constant matrices and **B**$_k$ and **D**$_k$ two non-constant matrices, whose time evolution is periodic of period $T_d = \frac{T}{\Delta t}$ ($\Delta t$ the data time step) and with $\tilde{w}_k = w_k + \epsilon_{y,k}$, where $\epsilon_{y,k}$ is the approximation error of the observation (supposed to be with zero mean, independent with moments of same order as the noises). The complete definition of the state-space of the approximation and the demonstration of the applicability of the Stochastic Subspace Identification (SSI) method is detailed in [13]. Here, an efficient version [16] of the SSI method [45] is used. The covariance-driven SSI [45] aims to identify the eigenmodes of the system through the sample correlations and is now presented. The first step is to construct a Hankel matrix filled by correlations, directly constructed from matrices gathering the observation data

$$ \hat{H} = \mathcal{Y}^{+} (\mathcal{Y}^{-})^{T}. \quad (11) $$

where $\mathcal{Y}^{+} \in \mathbb{R}^{(p+1)r \times N}$ and $\mathcal{Y}^{-} \in \mathbb{R}^{q \times N}$ are defined in [45]. $\hat{H}$ can be seen as the Hankel matrix filled with the correlations $\hat{R}_i$, the estimate of the correlation $R_i = \mathbb{E}(y_k y_{k-i}^T) = \mathbf{C}\mathbf{A}^{i-1}\mathbf{G}$, where $\mathbf{G} = \mathbb{E}(z_{k+1} y_k^T)$. From the estimation of the Hankel matrix, by using a singular value decomposition it is possible to identify the system matrices **A** and **C**. Then, the eigenmodes can be computed with the eigenvalue decomposition of **A**

$$ \mathbf{A} = \Psi [\mu] \Psi^{-1}. \quad (12) $$

The characteristic exponents $\lambda_i$, that are identified as the continuous-time eigenvalues are found from the discrete-time eigenvalues $\mu_i$ by $\lambda_i = \frac{\log(\mu_i)}{\Delta t}$. Then the natural frequency ($f_i$) and the damping ratio ($\xi_i$) of the associated
mode are defined such that $f_i = \frac{|\lambda_i|}{2\pi}$ and $\xi_i = -\frac{\Re(\lambda_i)}{|\lambda_i|}$. Finally, the mode shape matrix is found from $\Phi = C\Psi$. It should be remembered that when studying the approximation, the identified modes are the Fourier harmonics of the Floquet modes.

The presented identification method is theoretically defined for a wind turbine with a constant rotational speed, but it can be also applied to measurement from a wind turbine with a variable rotational speed, as it has been performed in [12]. As the identified parameters are obtained by correlations, the mean structure (associated with the mean rotational speed) is identified.

In this section, it has been concluded that the identification methods defined for LTI systems can be used for the identification of the eigenmodes of the LTP systems thanks to the approximation of the Floquet modes. Consequently, the damage detection methods developed for the LTI systems can also be used for the LTP systems.

### 3. Damage detection and localization method

In this section, the damage detection based on the local approach is presented with the detection, the localization, and the clustering of the parameters. Then, the adaptation of the method for operating wind turbine rotor is explained in addition to the definition of a mean residual to increase the damage detectability.

#### 3.1. Statistical damage detection

The objective of the damage detection method is to detect a change in the system parameters (that can represent the stiffness, damping or mass of the elements model for a mechanical system), represented as the vector $\theta$, with behavior assumed as

$$ 
\begin{aligned}
H_0: \theta &= \theta_0 \quad (\text{reference system}), \\
H_1: \theta &= \theta_0 + \frac{\delta}{\sqrt{N}} \quad (\text{damaged system}),
\end{aligned} 
\qquad (13) $$

where $\delta$ is unknown but fixed and $N$ is the length of the considered signal.

The change in the parameters will be sought through a residual $\zeta \in \mathbb{R}^l$, computed from a subspace identification. The impact of a change in the parameters on the residual can be modelled as a change in the mean of the residual [17], such that

$$ \zeta \sim \begin{cases} \mathcal{N}(0, \Sigma) : H_0 \\ \mathcal{N}(\mathcal{J}\delta, \Sigma) : H_1 \end{cases}, \qquad (14) $$

where $\delta \in \mathbb{R}^h$ is related to the unknown change under $H_1$, $\mathcal{J} = \frac{\partial\zeta}{\partial\theta}$ the sensitivity matrix of the residual with respect to the parameters and $\Sigma$ the covariance matrix of the residual under both $H_0$ and $H_1$.

Considering Equations (13) and (14) together implies that, for a given statistical change $\delta$, the larger the amount of data used, the smaller the detectable variation of the parameters can be detected. Conversely, for a given change in the parameters $\Delta\theta$, the more data are used, the higher the absolute mean of the residual will be.

Among many possibilities of Gaussian residuals, the residual was previously defined in [4] as $\zeta = \sqrt{N}\text{vec}(U_2^T \hat{H})$, with the Hankel matrix $\hat{H}$ defined in Equation (11) and $U_2$ the left null part of a SVD decomposition of $\hat{H}$.

To detect a change in the residual distribution, the Generalized Likelihood Ratio (GLR) is used [4],

$$ GLR(\zeta) = -2 \log \frac{\sup_{\theta \in H_0} p(\zeta|\theta_0)}{\sup_{\theta \in H_1} p(\zeta|\theta)}. \qquad (15) $$
As the distribution of the residual is assumed to follow a statistical Gaussian law, the GLR can be expressed as

$$ GLR(\zeta) = t = \zeta^T \Sigma^{-1} \mathcal{J} (\mathcal{J}^T \Sigma^{-1} \mathcal{J})^{-1} \mathcal{J}^T \Sigma^{-1} \zeta, \qquad (16) $$

with $\delta^* = (\mathcal{J}^T \Sigma^{-1} \mathcal{J})^{-1} \mathcal{J}^T \Sigma^{-1} \zeta$ the value that maximise the GLR test [6] and as such satisfying $\partial_{\delta} GLR(\zeta) = 0$. From this test, the Fisher matrix can be defined

$$ F = \mathcal{J}^T \Sigma^{-1} \mathcal{J}. \qquad (17) $$

This matrix can be seen as a ratio between the sensitivity and the covariance of the residual, where the diagonal terms give the detectability of each parameter component and the non-diagonal terms quantify the impact of one parameter component on the others and vice versa. This means that the higher the value on the diagonal, the more detectable a change in the parameter.

The statistical test follows a $\chi^2$ distribution such as

$$ t \sim \begin{cases} \chi^2(\nu, 0) : H_0, \\ \chi^2(\nu, \lambda) : H_1, \end{cases} \qquad (18) $$

where $\nu$ is the number of degrees of freedom of the distribution equal to the dimension of the parameter space (i.e. size of $\theta$), and $\lambda$ is the non-centrality parameter such that

$$ \lambda = \delta^T F \delta. \qquad (19) $$

Based on the properties of the $\chi^2$ distribution, the mean of the test denoted as $\bar{t}$, will be $\bar{t} = \nu$ under $H_0$ and $\bar{t} = \nu + \lambda$ under $H_1$.

From the theoretical distribution of the test (under both $H_0$ and $H_1$), it is possible a priori to assess the minimal change that will be detected based on a given confidence level, 95% in this paper. Here, the method defined in [28] is used. First, the upper bound of the interval that contains 95% of the test corresponding to the reference system $[0, t_{crit}]$ (leading to a probability of false alarm of 5%) is computed, with

$$ \int_{0}^{t_{crit}} f_{\chi^2}(\nu, 0)(t)dt = 0.95. \qquad (20) $$

Then, the minimum non-centrality parameter ($\lambda_{min}$) is defined such that

$$ \int_{0}^{t_{crit}} f_{\chi^2}(\nu, \lambda_{min})(t)dt = 0.05 \qquad (21) $$

and estimated by minimizing the following function

$$ f(\lambda) = \left| \int_{0}^{t_{crit}} f_{\chi^2}(\nu, \lambda)(t)dt - 0.05 \right|. \qquad (22) $$

Then, from Equation (19) the associated minimum statistical change can be assessed and for one parameter component

$$ \delta_{h \min} = \sqrt{\lambda_{min}/F_{hh}}, \qquad (23) $$

where $F_{hh}$ corresponds to the contribution of the h-th parameter component in the Fisher matrix. From the minimum 95%-level detectable statistical change, the associated minimum 95%-level detectable change in the h-th parameter component is

$$ \Delta\theta_{h \min} = \frac{\delta_{h \min}}{\sqrt{N}}. \qquad (24) $$
Also, from the value $t_{crit}$ it is possible to compute the probability of detection of a statistical change $\delta$ and the associated change $\Delta\theta$ in the physical parameter $\theta$ ($\Delta\theta = \theta - \theta_0$). From the distribution of the test and the non-centrality parameter associated with the change (Equation (19)), the probability of detection (POD) is defined by

$$ \mathrm{POD}(\Delta\theta) = \int_{t_{\mathrm{crit}}}^{\infty} f_{\chi^2}(\nu, N\Delta\theta^T F\Delta\theta) (t)dt. \qquad (25) $$

## 3.2. Damage localization by fault isolation

Once the damage detection is performed, it is possible to determine the parameter component that has changed. For such purpose, the direct localization/isolation test defined in [4] is evaluated, where each parameter component is individually tested for change. This test assumes that the change is restricted to one statistical change component $\delta_h$ ($\mathcal{J}\delta = \mathcal{J}_h\delta_h$), where $\mathcal{J}_h$ corresponds to the sensitivity of the parameter component of index $h$. Leading to

$$ t_h = \zeta^T \Sigma^{-1} \mathcal{J}_h (\mathcal{J}_h^T \Sigma^{-1} \mathcal{J}_h)^{-1} \mathcal{J}_h^T \Sigma^{-1} \zeta, \qquad (26) $$

with the following distribution

$$ t_h \sim \begin{cases} \chi^2(1, 0) : H_0 \\ \chi^2(1, \lambda) : H_1 \end{cases}. \qquad (27) $$

For the direct test on the other parameter components (that have not changed), their respective distribution is $t_{\bar{h}} \sim \chi^2(1, \bar{\lambda})$ under $H_1$, with $\bar{\lambda}$, the non-centrality parameter defined as [17]

$$ \bar{\lambda} = \delta_h^2 F_{h\bar{h}} F_{\bar{h}\bar{h}}^{-1} F_{h\bar{h}}^T, \qquad (28) $$

where $F_{h\bar{h}} = \mathcal{J}_h^T \Sigma^{-1} \mathcal{J}_{\bar{h}}$ and $F_{\bar{h}\bar{h}} = \mathcal{J}_{\bar{h}}^T \Sigma^{-1} \mathcal{J}_{\bar{h}}$. It means that under $H_1$, when the changed parameter component is tested, the test mean (denoted as $\bar{t}_h$) will be $\bar{t}_h = 1 + \lambda$, whereas $\bar{t}_h = 1 + \bar{\lambda}$ when an unchanged parameter component is tested.

## 3.3. Damage quantification

From the direct localization test, it is possible to compute a quantification test, based on the relation between the change of the mean of the residual and the fault. So, from [17] the fault quantification test is defined as

$$ \hat{\delta}_h = (\mathcal{J}_h^T \Sigma^{-1} \mathcal{J}_h)^{-1} \mathcal{J}_h^T \Sigma^{-1} \zeta. \qquad (29) $$

Under $H_1$ and with $\mathcal{J}\delta = \mathcal{J}_h\delta_h$ the distribution of the test is

$$ \hat{\delta}_h \sim \mathcal{N}(\delta_h, F_{hh}^{-1}). \qquad (30) $$

Otherwise, if the tested parameter is not the actual changed parameter, the quantification follows

$$ \hat{\delta}_h \sim \mathcal{N}(F_{hh}^{-1} F_{h\bar{h}} \delta_{\bar{h}}, F_{h\bar{h}}^{-1}). \qquad (31) $$

So, the associated fault quantification can be used jointly with the localization test to assess the false positive. For example, let us consider a structure parametrized by its stiffness (at each element, each component, etc.) and a damage that reduces one of those stiffnesses. If a parameter (i.e. one stiffness) has a non-negligible localization
test (i.e. δ with a high value) and an associated fault with a positive value, then this test is a false positive as the parametrization and the modeling only allow for a stiffness reduction.

### 3.4. Clustering of the redundant parameters

In a considered system parametrization, some parameters can have similar sensitivity with collinear sensitivity vectors, leading to the same values in the damage localization test in Equation (26). It implies that it is not possible to make a difference between these parameters in the localization. As it is not possible to distinguish them, a solution is to group parameters with the same sensitivity into clusters. In [1, 27] a method to group parameters with similar sensitivities into clusters has been defined based on the hierarchical clustering, as it is not necessary to define the number of clusters a priori, also the results can be analyzed with a dendrogram. To perform this clustering, the following distance is used

$$d_{i,j} = 1 - \frac{\tilde{\mathcal{J}}_i^T \tilde{\mathcal{J}}_j}{\|\tilde{\mathcal{J}}_i\| \|\tilde{\mathcal{J}}_j\|}, \qquad (32)$$

where $\tilde{\mathcal{J}}_i = \Sigma^{-1/2} \mathcal{J}_i$ is the normalized sensitivity of the i-th parameter. Once the cluster have been defined, the associated sensitivities can be defined as the cluster center, with

$$\tilde{\mathcal{J}}_k^c = \frac{1}{m_k} \sum_{i \in C_k} \tilde{\mathcal{J}}_i, \qquad (33)$$

the sensitivity of the k-th cluster ($C_k$) composed of $m_k$ parameters.

As the sensitivities of the clusters are normalized, the test definitions have to be adapted. Consequently, the damage detection, localization and quantification tests function of normalized sensitivities are defined such that

$$t = \zeta^T \Sigma^{-1/2} \tilde{\mathcal{J}} \left( \tilde{\mathcal{J}}^T \tilde{\mathcal{J}} \right)^{-1} \tilde{\mathcal{J}}^T \Sigma^{-1/2} \zeta \qquad (34)$$

$$t_h = \zeta^T \Sigma^{-1/2} \tilde{\mathcal{J}}_h \left( \tilde{\mathcal{J}}_h^T \tilde{\mathcal{J}}_h \right)^{-1} \tilde{\mathcal{J}}_h^T \Sigma^{-1/2} \zeta \qquad (35)$$

$$\hat{\delta}_h = \left( \tilde{\mathcal{J}}_h^T \tilde{\mathcal{J}}_h \right)^{-1} \tilde{\mathcal{J}}_h^T \Sigma^{-1/2} \zeta. \qquad (36)$$

### 3.5. New residual for wind turbine rotor

Since the objective is to detect a structural change in a wind turbine rotor, some specific quantities can be chosen to adapt the state-of-the-art method with the definition of a new residual. Based on previous works [42] and shown in [11], the phase shift and the amplitude of the blades edgewise bending mode shapes are one of the most sensitive indicators of changes in the rotor. Those two quantities are defined as follows

$$p_j^i = \operatorname{atan2}\left(\frac{\Im(\phi_j^i)}{\Re(\phi_j^i)}\right) \quad \text{and} \quad a_j^i = \sqrt{\Re(\phi_j^i)^2 + \Im(\phi_j^i)^2}, \qquad (37)$$

with $p_j^i$ and $a_j^i$ respectively the phase shift and the amplitude of the j-th DOF of the i-th mode shape ($\phi_j^i$). $\phi_j^i$ is identified using the identification method described in Section 2.3. It has to be noted that one DOF has to be chosen to define the reference. In this paper, the first DOF is the chosen reference, with $p_1^i = 0$ and $a_1^i = 1$ for all the mode shapes. As an illustration of their relevance, those quantities are displayed in Figure 1 for the DTU 10MW wind turbine model for three different mode shapes and for different damage sizes. The damage is defined here as a global loss in stiffness (up to 5%) of the third blade. The phase shift and the modal amplitude
are displayed in the form of a polar plot, where the radial coordinate corresponds to the modal amplitude and the angular coordinate to the phase shift. The first mode displayed corresponds to the collective edge plot (see Figure 1(a)), i.e. all the blades have the same motion. This can be deduced from Figure 1, where only one point per damage can be seen, meaning that the phase shift between the blades is null and the amplitude is the same. In this case, when the fault increases, the amplitude decreases. The second considered mode is the backward edge whirl (see Figure 1(b)). The deformed shape is also given in Figure 1(d). Here, a phase shift between the blades exists (i.e. they reach their maximum vibration amplitude at different times), and it is equal to 120° for the undamaged case. When a fault appears, one blade has an increase in amplitude while another observes a decrease. The phase is also impacted as it decreases for one and increases for another. Similar analysis can be done for the forward edge mode (see Figure 1(c)). It is worth mentioning here that the multi-blade coordinate transformation has not been used in this paper. Applying it will enforce frequency separation of the backward and forward whirling modes by 2Ω, and thus makes the identification easier. However, the identification has not been a source of problems in the current paper as we proceed SSI similarly as in [13] and focus on damage detection through residuals.

In Figure 1, each subplot corresponds to the evolution of the phase shift and amplitude of the mode shape of the first edge bending modes (bending of the blade in the rotor plane), namely the Collective, Backward and Forward edge modes, which differ by the phase shift between the bending of the blades. Each color corresponds to a scenario of global loss in stiffness of the third blade, namely 1%, 2.5%, and 5% stiffness loss. As the first blade defines the reference, the associated amplitude is always equal to one and the phase to zero, whatever the scenario due to the mode shape normalization. For the collective edge, no phase shift is observed, even with an increase in fault, and a loss in amplitude is observed for the third blade only. For the case with structural isotropy (blades with the same geometrical and material properties), and for the forward and backward modes, the three blades have the same amplitude and phase shift of 120°. However, a strong change in both amplitude and phase shift is observed when a fault appears. It illustrates the high sensitivity of those indicators to faults. As an indication, the natural frequencies of these modes are only varying around 10⁻² Hz for the same faults. The choice of the reference blade will have an impact on the values of phase shift and amplitude, but the global change compared to the isotropic state will be the same.

Let us define the vector $v_i$ that gathers the phase shifts and the amplitudes of the i-th mode shape (defined with r outputs)

$$v_i = [p_2^i \ p_3^i \ \dots \ p_r^i \ a_2^i \ a_3^i \ \dots \ a_r^i]^T. \qquad (38)$$

With this vector, it is possible to express a new residual

$$\zeta_i = \sqrt{N} (\hat{v}_i - v_i) \in \mathbb{R}^{2(r-1)}, \qquad (39)$$

where $\hat{v}_i$ is an estimate of $v_i$ using a data of length N. Then, a residual that gathers the information of n mode shapes can be defined, such that

$$\zeta = [\zeta_1^T \ \dots \ \zeta_n^T]^T = \sqrt{N} (\hat{V} - V) \in \mathbb{R}^{2n(r-1)}, \qquad (40)$$

where $V = [v_1^T \ \dots \ v_n^T]^T$ is the vector regrouping the phases and amplitudes of n mode shapes.
<figure><img src="image_2.png" /><figcaption>Figure 1: Polar plot figures of the evolution of the phase shift (angular position) and modal amplitude (radial axis) of the edgewise bending mode shapes against the third blade stiffness and visualization of the backward edge mode shape</figcaption></figure>

In real operational conditions, it can be difficult to know the theoretical value of V, so an estimate ($\hat{V}^0$) of this quantity in the reference state can be used instead. In this paper, an average of ($n_f$) different estimates is used

$$ \hat{V}^0 = \frac{1}{n_f} \sum_{j=1}^{n_f} \hat{V}^{(j)}, \qquad (41) $$

where $\hat{V}^{(j)}$ is an estimate of V using the j-th data set. From this, similarly to [46], a residual ($\tilde{\zeta}$) function of the estimated reference can be defined, that considers the uncertainties in the reference state

$$ \tilde{\zeta} = \sqrt{N} (\hat{V} - \hat{V}^0) = \zeta - \sqrt{N} (\hat{V}^0 - V). \qquad (42) $$

### 3.6. Residual distributions

$\hat{V}$ is an estimate of V using data of length N. Then, the distribution of the residual $\zeta$ (Equation (40)) is defined by the following distribution

$$ \zeta = \sqrt{N} (\hat{V} - V) \sim \begin{cases} \mathcal{N}(0, \Sigma) : H_0 \\ \mathcal{N}(\mathcal{J}\delta, \Sigma) : H_1 \end{cases}, \qquad (43) $$

with the assumption that the effect of the parameter change on V is linear.

Now, let us take into account the estimated reference into the distribution of the residual ($\tilde{\zeta}$ in Equation (42)). First, as $\hat{V}^{(j)}$ is computed using a data set of length N

$$ \sqrt{N} (\hat{V}^{(j)} - V) \sim \mathcal{N}(0, \Sigma). \qquad (44) $$
As the different estimates are computed with $n_f$ disconnected data sets, they are independent, so,

$$ \sum_{j=1}^{n_f} \sqrt{N} (\hat{V}^{(j)} - V) \sim \mathcal{N}(0, n_f \Sigma). \qquad (45) $$

Consequently, using Equation (41) the distribution of the estimated reference is

$$ \sqrt{N} (\hat{V}^0 - V) \sim \mathcal{N} \left( 0, \frac{\Sigma}{n_f} \right). \qquad (46) $$

Finally, as the data used for the computation of the estimates $\hat{V}$ and $\hat{V}^0$ are different, $\hat{V}$ and $\hat{V}^0$ are independent. So, the distribution of the residual $\tilde{\zeta}$ is

$$ \tilde{\zeta} \sim \begin{cases} \mathcal{N} \left( 0, \Sigma \left( 1 + \frac{1}{n_f} \right) \right) : H_0 \\ \mathcal{N} \left( \mathcal{J}\delta, \Sigma \left( 1 + \frac{1}{n_f} \right) \right) : H_1 \end{cases} \qquad (47) $$

### 3.7. Definition of the mean residual

With the presented method, the value of the detected parameters changes depending on two quantities, namely the covariance which is related to the external noise, which is assumed to remain unchanged, and the length of the signal. The latest can be adapted according to a targeted value. Otherwise, in the wind turbine field, it is common to have a fixed measurement length to be synchronized with the SCADA. One option is to concatenate the data in a unique data set and then identify the residual. Another option is to identify the residual for each data set of fixed-length and then compute a mean residual. The latter is chosen because it is more versatile and less time consuming.

The mean residual, function of the phase shifts and amplitudes of the edge bending modes is defined as

$$ \bar{\zeta} = \sqrt{N} (\bar{V} - V) \in \mathbb{R}^{2n(r-1)}, \qquad (48) $$

with

$$ \bar{V} = \frac{1}{nb_f} \sum_{j=1}^{nb_f} (\hat{V})_j, \qquad (49) $$

the mean of the phase shifts and amplitudes, computed from the system identification results obtained from $nb_f$ measurement sessions of data length $N$.

Let us define the distribution of this new residual, as it has been done in the previous section

$$ \sqrt{N} \left( \frac{1}{nb_f} \sum_{j=1}^{nb_f} (\hat{V})_j - \frac{1}{nb_f} \sum_{j=1}^{nb_f} V \right) \sim \begin{cases} \mathcal{N} \left( 0, \frac{\Sigma}{nb_f} \right) : H_0 \\ \mathcal{N} \left( \mathcal{J}\delta, \frac{\Sigma}{nb_f} \right) : H_1 \end{cases}, \qquad (50) $$

then

$$ \bar{\zeta} \sim \begin{cases} \mathcal{N} \left( 0, \frac{\Sigma}{nb_f} \right) : H_0 \\ \mathcal{N} \left( \mathcal{J}\delta, \frac{\Sigma}{nb_f} \right) : H_1 \end{cases} \qquad (51) $$

the distribution of the mean residual, computed with $nb_f$ identified residuals. Thereafter, the minimum 95% detectable fault associated with the mean residual can be expressed as a function of the minimum detectable fault defined in Equation (23),

$$ \bar{\delta}_{h \ min} = \frac{\delta_{h \ min}}{\sqrt{nb_f}}. \qquad (52) $$
Consequently, $nb_f$ can be adjusted to detect a specific fault or parameter change.

Analogously, a mean residual computed using $nb_f$ identified residuals with a reference estimated with $n_f$ identifications is defined as

$$ \bar{\zeta}^* \sim \begin{cases} \mathcal{N}\left(0, \Sigma\left(\frac{1}{nb_f} + \frac{1}{n_f}\right)\right) : H_0 \\ \mathcal{N}\left(\mathcal{J}\delta, \Sigma\left(\frac{1}{nb_f} + \frac{1}{n_f}\right)\right) : H_1 \end{cases}. \qquad (53) $$

To conclude, from a model study leading to the sensitivities of a damage sensitive residual and a signal processing based on the GLR, it is possible to perform damage detection and localization. The whole method can be summarised with the flowchart in Figure 2. The inputs and outputs of the method are given in red. Intermediary steps are in orange. Blue steps represent the damage detection or localization tests. Green steps represent decisions.

<figure><img src="image_3.png" /><figcaption>Figure 2: Flowchart of the damage detection and localization method</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 3: CAD model of the DTU 10 MW wind turbine under study (taken from [3, 36])</figcaption></figure>

# 4. Application: local loss of stiffness of one blade

## 4.1. *Model of the wind turbine and data computation*

In this article, the damage detection method will be applied on data computed with an aero-servo-elastic wind turbine numerical model. Precisely, the DTU 10MW wind turbine is considered [3, 36], implemented within the software OpenFAST [35]. Its CAD model is represented in Figure 3. The DTU 10MW is a synthetic model of a large wind turbine with two configurations, namely onshore and offshore on a jacket. The onshore version is used in this work, some major properties are recalled in Table 1 and all details can be found in [3, 36]. The mechanical part of the wind turbine model is modelled using OpenFAST’s Elastodyn module, where the multi-body system is defined by modal decomposition of each part, with a damping defined by Rayleigh damping. The final structural system is composed of 15 dofs: four dofs represent the modal amplitude in the fore-aft and side-side directions of the two first tower modes, two dofs represent the drivetrain torsional motion and nine dofs represent the modal amplitude of the first three modes in the flap and edge directions (three dofs per blade). In the state-space form, the system is of size 30.

<table><caption>Table 1: Characteristic values of the DTU 10MW wind turbine [3, 36]</caption><tr><td>Rated power</td><td>10 MW</td></tr><tr><td>Rotational speed</td><td>4 to 9.6 rpm</td></tr><tr><td>Rated tip speed</td><td>90 m/s</td></tr><tr><td>Wind speed</td><td>cut-in: 4 m/s<br>rated: 11.4 m/s<br>cut-out: 25 m/s</td></tr><tr><td>Rotor diameter</td><td>178.3 m</td></tr><tr><td>Hub diameter</td><td>5.6 m</td></tr><tr><td>Hub height</td><td>119 m</td></tr><tr><td>Rotor mass</td><td>229 tons (~ 41 tons per blade)</td></tr><tr><td>Nacelle mass</td><td>446 tons</td></tr><tr><td>Tower mass</td><td>605 tons</td></tr></table>

The objective is to create data with an environmental variability close to what can be seen on field. So, for each data set, the external load, which is the wind, is computed based on in-field data. The field data are ten-minutes mean wind speed from a *Météo France* buoy located in the gulf of Gascogne (See Figure 4). Those data have been collected in 2021 and are free to use [34]. Then, based on the mean wind speed the full field turbulent wind is computed using the software *Turbsim* [24], where the turbulence intensity is automatically
defined based on the IEC B norm (see IEC 61400-1 norm).

<figure><img src="image_5.png" /><figcaption>Figure 4: Localization of the meteorological buoys of Météo-France with accessible data [34], the used buoy is circled</figcaption></figure>

The simulated data are 600 s long and sampled at 25 Hz. The first 100 s are not used to remove the transient response of the wind turbine, so the signal length is $N = 25000$. To ensure data representing an operating wind turbine at the same operational point (e.g. same pitch angle), only data with an average wind speed between 4 and 11 m/s are kept. This range contains almost 70% of the speed values, as shown in Figure 5, with a histogram of all the wind speed values and an evolution of the wind speed over one month. So, it can be concluded that the selected mean wind speed range represents the majority of the conditions of a wind turbine, but also ensures similar structural configurations.

<figure><img src="image_6.png" /><figcaption>Figure 5: Wind speed values for the Gascogne buoy</figcaption></figure>

Finally, to perform the identification and the different statistical tests, the edge bending modes of the rotor are identified using the edge bending blade root moments. Those outputs are used as they are closer to a real instrumentation, for which it is easier to install sensors close to the hub. Before computing the damage detection and localization tests, the covariance matrix ($\Sigma$) and the reference ($\hat{V}^0$) have to be estimated. To do so, the phase shifts ($p_i$) and the amplitudes ($a_i$) of each mode are identified $n_f = 200$ times (with 200 different samples of measured wind speed), assuming that the structure is healthy, to estimate the reference and its covariance. Then, the reference is estimated as the mean of the identified phase shifts and amplitudes over the $n_f$ data sets.
And, with the same identified values, the covariance matrix is estimated such that

$$ \hat{\Sigma} = \frac{1}{n_f - 1} \sum_{j=1}^{n_f} (\hat{V}_j - \hat{V}^0)^T (\hat{V}_j - \hat{V}^0), \qquad (54) $$

with $\hat{V}^0$ defined in Equation (41).

## 4.2. Parametrization and sensitivity analysis

For damage detection and localization, a parametrization of the system has to be defined, and the associated sensitivity matrix has to be computed. Here, each blade is divided into ten sections of the same lengths and each section is parametrized by a ratio of the stiffness in both flap and edge. This ratio is equal to one in the reference state (i.e. healthy) and lower than one in the damaged state. Thus, it represents 30 parameters in total to parametrize the full wind turbine. All those parameters are numbered section by section and in the order of the blades, so the number of the parameter of section $i$ of blade $j$ is $3i + j$ with $i \in [0, 9]$ and $j \in [0, 2]$. To compute the sensitivity matrix, the first order approximation is used, with

$$ \mathcal{J}_h = \frac{\zeta(\theta + \Delta\theta_h) - \zeta(\theta)}{\Delta\theta_h}, \quad \text{with } \Delta\theta_h = 0.05 \qquad (55) $$

Then, with the sensitivity matrix and the covariance matrix of the residual, it is possible to compute the Fisher matrix associated with the residual (defined in Equation (17)) and the minimum 95% detectable faults (Equation (24)).

<figure><img src="image_7.png" /><figcaption>Figure 6: Sensitivity analysis based on the parametrization</figcaption></figure>

In Figure 6a, the color map of the Fisher matrix is displayed. As a reminder, the higher the value of a Fisher coefficient of a parameter, the more sensitive the test will be to a variation in that parameter. Because the non-centrality of the damage detection test depends on that value. It can be seen that for parameters 21 to 29, the corresponding matrix coefficients are lower. These components correspond to the parameters at the end of the blades. However, since the detectable fault directly depends on the Fisher matrix components (see Equation (23)), the detectability will be lower, that is, the minimum detectable faults will be higher for these parameters.

This is shown in Figure 6b, where the values of the minimum detectable fault for each parameter are given (computed from Equation (24)). One can see that the values are higher for parameters 21 to 29, with values
higher than 100%. The delimitation between the elements with less than 100% and the elements over 100% is
clearly defined as the last segment. Consequently, with the selected residual, measurements, and signal length,
it is not possible to detect a loss of stiffness in the last 30% of the blades. This means that a loss of stiffness
at the tip of the blade cannot be detected from the phase shift and amplitude of the edge bending mode shapes
(nor from the frequencies of these modes after testing) without increasing significantly the total number of
data sets for the residual computation. As a consequence, the corresponding parameters are removed from the
parametrization in the following, leading to 21 parameters representing the stiffness of the first 70% of the three
blades. The incapacity of detecting a fault with this instrumentation and approach on the last 30% of the blades
is a drawback. However, as presented in [25] the end of the blade is significantly more susceptible to erosion
instead of crack appearance, and in [31] a review on the prevention of blade erosion is presented.

To avoid redundant parameters over the 21 used in the parametrization of the model, clustering is performed
to group parameters of similar sensitivity based on the method detailed in Section 3.4. Results are displayed
with a dendrogram given in Figure 7a, where the distance is obtained by ward linkage. From the clustering,
three main clusters can be distinguished, each one grouping the stiffness of a blade. Smaller clusters are not
distant enough for exploitation and propose a finer localization. The final parametrization is composed of three
parameters and is displayed in Figure 7b. It can be concluded from the clustering that it should be possible to
localize the damaged blade, however, it should not be possible to assess which part of the blade is damaged as
the sensitivity of parameters on the same blade are too similar. Those assumptions will be confirmed in the next
sections.

<figure><img src="image_8.png" /><figcaption>Figure 7: Results of the clustering of the rotor stiffness parameters</figcaption></figure>

4.3. Damage detection

To perform damage detection, the different tests are computed on data of the reference system, to assess the distributions of the tests and the value of the 95% upper bound. Then, two damaged systems are investigated. The first case corresponds to a 7.5% stiffness reduction of the second blade in the fifth section (40% to 50% of the blade from the hub). The second case corresponds to a 7.5% of stiffness reduction of the third blade in the second section (10% to 20%). The order of magnitude of the stiffness losses can be linked to a shear web
debonding, as it is presented in [33]. The eigenmodes used in OpenFast modal decomposition are changed for each damage case. For all the different damage cases, 100 test results are computed to obtain the progressive damage detection plot.

As the impact of the damage on the residual is low (theoretical POD computed using Equation (25) is lower than 50%), it is chosen to use $nb_f = 10$ to compute the mean residual and increase the damage detectability. The objective is to simulate continuous monitoring of the wind turbine rotor, so, the mean residual is computed using a simple moving average with

$$ \bar{\zeta}(k) = \frac{1}{nb_f} \sum_{p=0}^{nb_f-1} \zeta_{k+p}, k \in [1, 100 - (nb_f - 1)]. \quad (56) $$

Then, the time evolution of the damage detection test using the moving average is displayed in Figure 8. Firstly, one can see that the mean residual and the moving average cause tests that are a function of values of the reference and damaged structure, designed as transition test, with a total of $nb_f - 1$ tests. Secondly, with the first damage (Figure 8.a), it can be seen that the distribution of the test corresponding to the damaged structure has shifted, with an empirical probability of detection (POD) of 100% (i.e. all the tests are over the 95% threshold), which means that the damage is detected. For the second damage (Figure 8.b), the shift of the distribution of the test corresponding to the damaged structure is lower, but the damage is still detected with a POD of around 70%. The shift of the distribution of both simulated damages can be seen in Figure 9, where the histograms of the different tests (Figure 8) are displayed.

<figure><img src="image_9.png" /><figcaption>Figure 8: Evolution of the damage detection test of the moving average mean residual</figcaption></figure>

From the results of the damage detection, it can be concluded that the defined residual and the statistical test can be used to detect a change in a wind turbine rotor through the change in the distributions. However, for the second simulated damage, the empirical POD is lower than the theoretical one, i.e. around 70% empirically compared to more than 99% theoretically. This is due to a lower shift than expected and can be the consequence of the assumption of a linear impact of the parameter change on the mean of the residual.

Also, as the data sets are computed with different wind speeds and consequently different rotational speeds, it can be concluded that the residual and consequently the damage detection test is not sensitive to the envi-
<figure><img src="image_10.png" /><figcaption>Figure 9: Histograms of the detection test for the reference and damaged structure</figcaption></figure>

ronmental conditions, as long as the wind turbine is operating in the conditions that are described in Section 4.1.

Now the damage localization test will be exploited to confirm that the damage can be localized and confirm
the clustering choice of the parameters.

4.4. Damage localization

Along with the damage detection, the damage localization test has been computed, where each cluster (each blade) is tested. Also, the damage quantification is performed jointly with the localization (see Equation (29)). Firstly, the localization and quantification test results of the first damage case are studied. In Figure 10, the boxplot of each localization test and the mean of the quantification tests are displayed. With the localization, one can see that the distribution of the localization test corresponding to the second blade is the test with the higher shift. So the second blade can be designated as damaged. Furthermore, the quantification corresponding to the second blade is the only one with a negative value. It confirms that the parameter change is in the second blade as it is the only blade associated with a change that has a physical meaning. Also, the mean identified damage is around 6.5%, which is close to the simulated damage.

<figure><img src="image_11.png" /><figcaption>Figure 10: Results of the damage localization and quantification for a loss of 7.5% of stiffness in the section 5 of the second blade</figcaption></figure>

To confirm the previous results, localization tests are also performed for the second simulated damage, with
results displayed in Figure 11. With the box plots of the localization, it can be seen that the test corresponding
to the third blade is the test associated with the highest distribution shift. Then, the associated mean damage quantification is only negative for the third blade. So, analogously to the previous application, it can be concluded that the third blade is damaged, which corresponds to the damaged blade in the simulations. The mean value of the identified damage is 3.5%, which is the same order of magnitude as the simulated damage. However, the estimated value is further away from the simulated damage than in the previous case.

<figure><img src="image_12.png" /><figcaption>Figure 11: Results of the damage localization and quantification for a loss of 7.5% of stiffness in the section 2 of the third blade</figcaption></figure>

### 4.5. Alternative selection of sensors

The presented results have been obtained using the blade root moments as unique sensors. However, both detection and localization strongly depend on the sensor choice. The benefit of completing the instrumentation will be demonstrated by using accelerations of the blades, along the edge, at 40% and 80% of the blade span. It is expected to yield better results in terms of detectability and localization as more information is available and more distributed along the blades.

The first step is to perform the clustering once again, as the sensitivities depend on the sensor type and their location. In Figure 12a, the dendrogram of the hierarchical clustering is displayed. Firstly, one can see that there are the same sub-clusters in the three main clusters, which means that it seems that it is possible to define some clusters inside the blades, leading to future more accurate localization. Secondly, the distance between the clusters is significantly more important with this instrumentation compared to the previous one (see Figure 7a), which means that the sensitivities are less collinear than with the blade root moments. So, it is chosen to define two clusters per blade, as illustrated in Figure 12b. Also, as with the previous instrumentation, the theoretical minimum 95% detectable faults for the last 30% of the blades are too large (more than 100%), so associated parameters are removed of the parametrization.

Once the parametrization has been defined, the damage detection can be performed, once again with an estimated reference with $n_f = 200$ and a mean residual with $nb_f = 10$. Once the test is performed, an empirical POD of 100% is obtained for the two simulated damages. Histograms are given in Figure 13, where the shift between the distributions is obvious. Compared to the previous instrumentation (see Figure 9), the shift of the test corresponding to the damaged structure is more important and always over the 95% threshold. This means, that using two accelerations per blade instead of one root moment gives more information and gives better damage detectability.
<figure><img src="image_13.png" /><figcaption>Figure 12: Results of the clustering of the rotor stiffness parameters, use of accelerations</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 13: Histograms of the detection test for the reference and damaged structure, use of the blade accelerations</figcaption></figure>

The final step is to confirm that the localization has been improved. From the clustering, it has been chosen that the parametrization will be composed of two clusters per blade. So, for the damage localization, each cluster is tested using Equation (26), with conjointly the estimation of the associated parameter change. In Figure 14, the results of the damage localization and quantification of one damage case are displayed. Thus, it can be seen that the damage has been well localized, where the cluster which includes the position of the simulated damage is the cluster with the maximum shift. Also, the damage quantification gives an equivalent estimation of the damage compared to the use of the blade root moments, which is close to the simulated damage. Similar results are obtained for the second simulated damage, with the localization of the cluster which includes the position of the simulated damage and an associated quantification close to the one obtained with the moments.

It can be concluded that, with a larger number of sensors and different measured quantities, it is possible to improve the damage detectability and to localize the damage with a precision of almost half of a blade. Furthermore, it shows that despite a limited number of sensors, the blade root moment gives an equivalent
<figure><img src="image_15.png" /><figcaption>Figure 14: Results of the damage localization and quantification for a loss of 7.5% of stiffness in the section 5 of the second blade, use of the accelerations</figcaption></figure>

damage detectability and enables the localization of the damaged blade. It also has been tested to couple blade root moments and accelerations. If one acceleration per blade is used with the moments (located at 40% or 80% of the blade length) there is no major improvement in the damage detection and the damage localization is equivalent to the one obtained using only the blade root moments. Then it has been tested to couple the blade root moments with both accelerations located at 40% and 80% of the blade length. The damage detection and localization are equivalent to those obtained using only both accelerations on the blades, with just a higher shift of the damage detection test. So, using both accelerations on the blades (with or without the moments) is the only configuration that enables localizing the damage with the best resolution.

## 5. Conclusion

With the identification of the LTP system based on the approximation of the Floquet modes, it is concluded that it is possible to adapt the method based on the local approach to perform damage detection and localization on an operating wind turbine rotor. For this purpose, a new residual based on the most promising damage-sensitive features of a wind turbine rotor, i.e. the phase shift and amplitudes of the edge bending modes has been defined. To increase the damage detectability, a mean residual is defined. Thus, it has been possible to detect faults in the wind turbine rotor and to identify the damaged blade. The damage quantification associated with the damage localization does not allow us to have an exact estimation of the simulated damage, but gives a good estimation of the order of magnitude of this damage. This is the consequence of the assumption of the linear impact of the parameter on the residual. Nevertheless, the results of the quantification can be used as an initialization point for model updating [21]. In the present work, the wind turbine modal properties remain rather limited. As a perspective, it would be relevant to investigate the robustness of the proposed methodology when more complex wind turbine models are considered (using BeamDyn or HAWC2 for example). This would be a preliminary step before considering real-world measurements. Also, a challenge for future research would be to assess the validity of the method when the wind turbine experiences several damages.

The performance of the presented method relies on the sensors used in the identification to construct the residual. This explains why better damage detectability and localization resolution have been obtained using two accelerations per blade instead of the blade root moments. This means that the rotor instrumentation can
be chosen to maximize the performance of the damage detection and localization [29].

Damage detection and localization can also be performed for other types of defaults, like an addition of mass or a pitch misalignment on a blade, with the parametric or the non-parametric test (defined in [17]) depending on the availability of the sensitivity regarding such damages. A challenge for future research would be to combine several types of damage, and to assess the validity of the approach on real-world measurements.

## References

[1] S. Allahdadian, M. Döhler, C. Ventura, and L. Mevel. Towards robust statistical damage localization via model-based sensitivity clustering. *Mechanical Systems and Signal Processing*, 134:106341, 2019.

[2] E. Artigao, S. Martín-Martínez, A. Honrubia-Escribano, and E. Gómez-Lázaro. Wind turbine reliability: A comprehensive review towards effective condition monitoring development. *Applied energy*, 228:1569–1583, 2018.

[3] C. Bak, F. Zahle, R. Bitsche, T. Kim, A. Yde, L. C. Henriksen, M. H. Hansen, J. P. A. A. Blasques, M. Gaunaa, and A. Natarajan. The dtu 10-mw reference wind turbine. In *Danish wind power research 2013*, 2013.

[4] M. Basseville, M. Abdelghani, and A. Benveniste. Subspace-based fault detection algorithms for vibration monitoring. *Automatica*, 36(1):101–109, 2000.

[5] M. Basseville, L. Mevel, and M. Goursat. Statistical model-based damage detection and localization: subspace-based residuals and damage-to-noise sensitivity ratios. *Journal of sound and vibration*, 275(3-5):769–794, 2004.

[6] A. Benveniste, M. Basseville, and G. Moustakides. The asymptotic local approach to change detection and model validation. *IEEE Transactions on Automatic Control*, 32(7):583–592, 1987.

[7] C. L. Bottasso and S. Cacciola. Model-independent periodic stability analysis of wind turbines. *Wind Energy*, 18(5):865–887, 2015.

[8] S. Cacciola, I. M. Agud, and C. L. Bottasso. Detection of rotor imbalance, including root cause, severity and location. In *Journal of Physics: Conference Series*, page 072003. IOP Publishing, 2016.

[9] S. Cacciola, C. Riboldi, and A. Croce. Monitoring rotor aerodynamic and mass imbalances through a self-balancing control. In *Journal of Physics: Conference Series*, volume 1037, page 032041. IOP Publishing, 2018.

[10] A. Cadoret, E. Denimal, J.-M. Leroy, J.-L. Pfister, and L. Mevel. Linear time invariant approximation for subspace identification of linear periodic systems applied to wind turbines. In *11th IFAC Symposium on Fault Detection, Supervision and Safety for Technical Processes - SAFEPROCESS*, Pafos, Cyprus, 2022.

[11] A. Cadoret, E. Denimal, J.-M. Leroy, J.-L. Pfister, and L. Mevel. Mode shape phase change detection in wind turbine under anisotropy variation. In *ISMA 2022-International Conference on Noise and Vibration Engineering*, Leuven, Belgium, 2022.

[12] A. Cadoret, E. Denimal, J.-M. Leroy, J.-L. Pfister, and L. Mevel. Periodic system approximation for operational modal analysis of operating wind turbine. In *EWSHM 2022 - 10th Workshop on Structural Health Monitoring*, Palermo, Italy, 2022.

[13] A. Cadoret, E. D. Goy, J.-M. Leroy, J.-L. Pfister, and L. Mevel. Linear time periodic system approximation
based on floquet and fourier transformations for operational modal analysis and damage detection of wind turbine. *Mechanical Systems and Signal Processing*, 212:111157, 2024.

[14] M. Döhler, F. Hille, L. Mevel, and W. Rücker. Structural health monitoring with statistical methods during progressive damage test of s101 bridge. *Engineering Structures*, 69:183-193, 2014.

[15] M. Döhler, F. Hille, L. Mevel, and W. Rücker. Structural health monitoring with statistical methods during progressive damage test of s101 bridge. *Engineering Structures*, 69:183-193, 2014.

[16] M. Döhler and L. Mevel. Fast multi-order computation of system matrices in subspace-based system identification. *Control Engineering Practice*, 20(9):882-894, 2012.

[17] M. Döhler, L. Mevel, and Q. Zhang. Fault detection, isolation and quantification from gaussian residuals with application to structural damage diagnosis. *Annual Reviews in Control*, 42:244-256, 2016.

[18] Y. Du, S. Zhou, X. Jing, Y. Peng, H. Wu, and N. Kwok. Damage detection techniques for wind turbine blades: A review. *Mechanical Systems and Signal Processing*, 141:106445, 2020.

[19] G. Floquet. Sur la théorie des équations différentielles linéaires. In *Annales Scientifiques de L'École Normale Supérieure*, volume 8, pages 3-132, 1879.

[20] D. García Cava, L. D. Avendaño-Valencia, A. Movsessian, C. Roberts, and D. Tcherniak. On explicit and implicit procedures to mitigate environmental and operational variabilities in data-driven structural health monitoring. In *Structural Health Monitoring Based on Data Science Techniques*, pages 309-330. Springer, 2022.

[21] S. Greś, M. Döhler, and L. Mevel. Statistical model-based optimization for damage extent quantification. *Mechanical Systems and Signal Processing*, 160:107894, 2021.

[22] G. R. Hübner, H. Pinheiro, C. E. de Souza, C. M. Franchi, L. D. Da Rosa, and J. P. Dias. Detection of mass imbalance in the rotor of wind turbines using support vector machine. *Renewable Energy*, 170:49-59, 2021.

[23] A. Jhinaoui, L. Mevel, and J. Morlier. A new ssi algorithm for lptv systems: Application to a hinged-bladed helicopter. *Mechanical Systems and Signal Processing*, 42(1-2):152-166, 2014.

[24] B. Jonkman. Turbsim user's guide v2. 00.00. *Natl. Renew. Energy Lab*, 2014.

[25] P. Kaewniam, M. Cao, N. F. Alkayem, D. Li, and E. Manoach. Recent advances in damage detection of wind turbine blades: A state-of-the-art review. *Renewable and Sustainable Energy Reviews*, 167:112723, 2022.

[26] J. Kusnick, D. E. Adams, and D. T. Griffith. Wind turbine rotor imbalance detection using nacelle and blade measurements. *Wind Energy*, 18(2):267-276, 2015.

[27] A. Mendler, M. Döhler, C. Ventura, and L. Mevel. Clustering of Redundant Parameters for Fault Isolation with Gaussian Residuals. In *IFAC 2020 - 21st International Federation of Automatic Control World Congress*, Berlin, Germany, July 2020.

[28] A. Mendler, M. Döhler, and C. E. Ventura. A reliability-based approach to determine the minimum detectable damage for statistical damage detection. *Mechanical Systems and Signal Processing*, 154:107561, 2021.

[29] A. Mendler, M. Döhler, and C. E. Ventura. Sensor placement with optimal damage detectability for statistical damage detection. *Mechanical Systems and Signal Processing*, 170:108767, 2022.

[30] A. Mendler, S. Greś, M. Döhler, and S. Keßler. On the probability of localizing damages based on mode shape changes. In *European Workshop on Structural Health Monitoring*, pages 233-243. Springer, 2023.
[31] L. Mishnaevsky Jr, C. B. Hasager, C. Bak, A.-M. Tilg, J. I. Bech, S. D. Rad, and S. Fæster. Leading edge erosion of wind turbine blades: Understanding, prevention and protection. *Renewable Energy*, 169:953-969, 2021.

[32] A. Movsessian, D. G. Cava, and D. Tcherniak. An artificial neural network methodology for damage detection: Demonstration on an operating wind turbine blade. *Mechanical Systems and Signal Processing*, 159:107766, 2021.

[33] N. Myrent, D. E. Adams, and D. T. Griffith. Wind turbine blade shear web disbond detection using rotor blade operational sensing and data analysis. *Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences*, 373(2035):20140345, 2015.

[34] Météo-France. Observations en mer. https://donneespubliques.meteofrance.fr/?fond=produit&id produit=95&id_rubrique=32. Accessed: 2022-04-27.

[35] National Renewable Energy Laboratory - NREL. Openfast v2.5.0. https://github.com/OpenFAST/openfast/releases/tag/v2.5.0. Accessed: 2021-06-03.

[36] T. U. of Denmark DTU. dtu-10mw-rwt. https://rwt.windenergy.dtu.dk/dtu10mw/dtu-10mw-rwt. Accessed: 2021-06-03.

[37] R. Riva, S. Cacciola, and C. Bottasso. Periodic stability analysis of wind turbines operating in turbulent wind conditions. *Wind Energy Science*, 1:177-203, 2016.

[38] R. Riva, S. Cacciola, and A. Croce. High-resolution periodic mode shapes identification for wind turbines. In *Journal of Physics: Conference Series*, volume 1037, page 062002. IOP Publishing, 2018.

[39] C. Roberts, D. G. Cava, and L. D. Avendaño-Valencia. Understanding the influence of environmental and operational variability on wind turbine blade monitoring. In *European Workshop on Structural Health Monitoring: Special Collection of 2020 Papers-Volume 1*, pages 109-118. Springer, 2020.

[40] P. F. Skjoldan and O. A. Bauchau. Determination of modal parameters in complex nonlinear systems. 2011.

[41] P. F. Skjoldan and M. H. Hansen. On the similarity of the coleman and lyapunov-floquet transformations for modal analysis of bladed rotor structures. *Journal of Sound and Vibration*, 327(3):424-439, 2009.

[42] D. Tcherniak. Rotor anisotropy as a blade damage indicator for wind turbine structural health monitoring systems. *Mechanical Systems and Signal Processing*, 74:183-198, 2016.

[43] D. Tcherniak and L. L. Mølgaard. Active vibration-based structural health monitoring system for wind turbine blade: Demonstration on an operating vestas v27 wind turbine. *Structural Health Monitoring*, 16(5):536-550, 2017.

[44] M. D. Ulriksen, D. Tcherniak, and L. Damkilde. Damage detection in an operating vestas v27 wind turbine blade by use of outlier analysis. In *2015 IEEE Workshop on Environmental, Energy, and Structural Monitoring Systems (EESMS) Proceedings*, pages 50-55. IEEE, 2015.

[45] P. van Overschee and B. de Moor. Subspace algorithms for the stochastic identification problem. *Automatica*, 29(3):649-660, 1993.

[46] E. Viefhues, M. Döhler, F. Hille, and L. Mevel. Statistical subspace-based damage detection with estimated reference. *Mechanical Systems and Signal Processing*, 164:108241, 2022.

[47] B. Wondra, J. Rupfle, A. Emiroglu, and C. Grosse. Analysis of icing on wind turbines by combined wireless and wired acceleration sensor monitoring. In *European Workshop on Structural Health Monitoring: EWSHM*
2022-Volume 1, pages 143–155. Springer, 2022.