Alexis Lebranchu
email: alexis.lebranchu@valemo.fr

Sylvie Charbonnier
email: sylvie.charbonnier@gipsa-lab.grenoble-inp.fr

Christophe Bérenguer
email: christophe.berenguer@grenoble-inp.fr

Frédéric Prevost
email: frederic.prevost@valemo.fr

A Combined Mono-and Multi-Turbine Approach for Fault Indicator Synthesis and Wind Turbine Monitoring Using SCADA Data

Keywords: Wind turbine monitoring, Wind farm monitoring, SCADA data, Fault detection, Condition monitoring, Performance evaluation

The monitoring of wind turbines using SCADA data has received lately a growing interest from the fault diagnosis community because of the very low cost of these data, which are available in number without the need for any additional sensor. Yet, these data are highly variable due to the turbine constantly changing its operating conditions and to the rapid fluctuations of the environmental conditions (wind speed and direction, air density, turbulence, ...). This makes the occurrence of a fault difficult to detect. To address this problem, we propose a multi-level (turbine and farm level) strategy combining a mono-and a multi-turbine approach to create fault indicators insensitive to both operating and environmental conditions. At the turbine level, mono-turbine residuals (i.e. a difference between an actual monitored value and the predicted one) obtained with a normal behavior model expressing the causal relations between variables from the same single turbine and learnt during a normal condition period are calculated for each turbine, so as to get rid of the influence of the operating conditions. At the farm level, the residuals are then compared to a wind farm reference in a multi-turbine approach to obtain fault indicators insensitive to environmental conditions. Indicators for the objective performance evaluation are also proposed to compare wind turbine fault detection methods, which aim at evaluating the cost/benefit of the methods from a production manager's point of view. The performance of the proposed combined mono-and multiturbine method is evaluated and compared to more classical methods proposed in the literature on a large real data set made of SCADA data recorded on a French wind farm during four years : it is shown than it can improve the fault detection performance when compared to a residual analysis limited at the turbine level only

Introduction

Renewable energy has been a growing sector for several years, because of the necessity to reduce CO2 emission in the near future. The electrical power produced by wind turbines has been multiplied by 10 in the last 10 years. A significant amount of the electricity consumed in the world now relies on the electric power produced by wind farms, which have to be operational all along the year. Failures may cause important production losses, mainly due to the damages they cause and the time it takes to repair, which are no longer acceptable. This calls for a drastic change in maintenance solutions, which must switch from periodic and corrective to condition-based. One of the motivations, and a perspective, for the work presented in this paper is to develop a fault detection procedure that can deliver information on a developing fault early enough so that it can be used for condition-based or predictive maintenance decision making : indeed an in-advance detection allows the maintenance decision-maker to better plan the maintenance operations.

From the sensor technology point of view, 3 technologies are possible to monitor a wind turbine, [START_REF] Tavner | Offshore Wind Turbines -Reliability, availability and maintenance[END_REF]:

• Using high rate data from a Condition Monitoring System (CMS) monitoring e.g. vibration: the use of these data for fault detection usually relies on signal processing methods aiming at identifying the signature of a fault, e.g. in the signal spectrum;

• Medium rate 1-second SCADA data: these data can be used for fault detection using for example model-based fault detection and isolation (FDI) methods, [START_REF] Iserman | Fault Diagnosis Systems -An Introduction from Fault Detection to Fault Tolerance[END_REF];

• Slow rate 10-minutes averaged SCADA data: in this case, the measured quantities are averaged on 10-minutes length windows.

Several condition monitoring systems (CMS) are available on the market [START_REF] Tchakoua | Wind Turbine Condition Monitoring: State-of-the-Art Review, New Trends, and Future Challenges[END_REF][START_REF] Qiao | A Survey on Wind Turbine Condition Monitoring and Fault Diagnosis -Part I: Components and Subsystems[END_REF][START_REF] Qiao | A Survey on Wind Turbine Condition Monitoring and Fault Diagnosis -Part II: Signals and Signal Processing Methods[END_REF][START_REF] Tavner | Offshore Wind Turbines -Reliability, availability and maintenance[END_REF]. They are based on vibration analysis for the most part. They require additional sensors to be placed on each nacelle to measure the vibrations of several components of the drive train. The highly sampled acceleration signals they produce must be analyzed by an expert using signal processing methods, so as to detect faults in progress. Indeed, the automatic decision systems developed still generate a high number of false alarms. These make them a costly solution to implement [START_REF] Yang | Wind turbine condition monitoring: technical and commercial challenges[END_REF][START_REF] Mcmillan | Quantification of condition monitoring benefit for offshore wind turbines[END_REF].

On the opposite, SCADA systems have been integrated in wind farms since the emergence of wind energy. They generate loads of cheap data containing useful information on the turbines state but the data low sampling rate, one average measure every 10 minutes, is not adapted to an accurate in-depth monitoring of wind turbines. However, their availability for free and the breadth of SCADA data scope have encouraged researchers to propose solutions to create SCADA based fault indicators [START_REF] Yang | Wind turbine condition monitoring by the approach of SCADA data analysis[END_REF][START_REF] Tavner | Offshore Wind Turbines -Reliability, availability and maintenance[END_REF]. The work presented in this paper is developed for 10-minutes averaged SCADA, classically recorded for production supervision purposes.

SCADA data monitoring relies on the numerous variables measured mainly for production supervision purposes. Fault monitoring approaches based on SCADA data differ in the way the data can be merged to synthesize a fault indicator residual, and can be categorized in "mono-turbine methods" or "multi-turbine methods", as proposed in [START_REF] Lebranchu | Review and analysis of SCADA data-based methods for health monitoring of wind turbines[END_REF]. Mono-turbine methods combine variables recorded on a unique wind turbine while multi-turbine methods combine variables coming from all the turbines within a wind farm or subset of a wind farm. In this work, it is proposed to use 10-minutes SCADA data in an innovative way, combining a mono-turbine (or turbine-level) and a multi-turbine (or farm-level) information processing.

Wind turbines SCADA variables are highly non-stationary because of the frequent changes in operational conditions and of the variations in external conditions. Following a classical diagnosis approach [START_REF] Iserman | Fault Diagnosis Systems -An Introduction from Fault Detection to Fault Tolerance[END_REF], mono-turbine methods aim at explaining a SCADA variable evolution as a function of other SCADA variables recorded on the same turbine and thus generate residuals (unexplained part) for diagnosis purposes. Variables may be linked by causal relations -a change in some variables induces a change in other variables-or by similarity relations -two variables evolve in the same way because they are submitted to the same excitations, [START_REF] Lebranchu | Review and analysis of SCADA data-based methods for health monitoring of wind turbines[END_REF]. The most common causal model in wind energy is the power curve [START_REF] Tavner | Offshore Wind Turbines -Reliability, availability and maintenance[END_REF], which expresses the link between the wind speed and the active power produced by the turbine. It can be used as a visualization tool by comparing the power curve built with data gathered during the current period with the reference curve provided by the constructor [START_REF] Mclaughlin | Wind farm performance assessment in the real world[END_REF]. It can also be used to create fault indicators by measuring the difference between the produced power and its value predicted by a model using the wind speed. In order to reduce the dispersion of the power data, additional input variables may be added such as the air density, the wind direction [START_REF] Cambron | Power curve monitoring using weighted moving average control charts[END_REF][START_REF] Butler | Exploiting SCADA system data for wind turbine performance monitoring[END_REF][START_REF] Yang | Wind turbine condition monitoring by the approach of SCADA data analysis[END_REF], the rotational speed or the pitch angle [START_REF] Elijorde | A Wind Turbine Fault Detection Approach Based on Cluster Analysis and Frequent Pattern Mining[END_REF]. The model may be a simple polynomial approximation, a Gaussian process model or an artificial neural network [START_REF] Li | Using neural networks to estimate wind turbine power generation[END_REF]. The fault indicators may be the difference at each sampling time between the measured power and the expected power or it can be a distance between a reference power curve and an on line curve built with the most current measures [START_REF] Yang | Wind turbine condition monitoring by the approach of SCADA data analysis[END_REF]. The limitation of this approach comes mainly from (i) the length of the training period, which has to be large in order to cover all the wind speed and temperature ranges, and (ii) the difficulty to localize the fault once it is detected with the power curve. To address the localization issue, one solution is to split the wind turbine into smaller independent systems, such as the rotor, the gearbox, the generator, the transformer and the convertor and to build models of these reduced systems, [START_REF] Schlechtingen | Wind turbine condition monitoring based on SCADA data using normal behavior models. Part 1: System description[END_REF][START_REF] Schlechtingen | Wind turbine condition monitoring based on SCADA data using normal behavior models. Part 2: Application examples[END_REF]. One popular variable to be modeled is the temperature of some components, whose variations can be explained by changes in the operational conditions or in the outside temperature. Models explaining the temperature variations use at least the produced power, the nacelle temperature and the train rotation speed as input variables. They are learnt on data measured when the turbine operates in normal conditions. The fault indicator is defined as the difference between the actual measured temperature and the expected temperature, named residual. Models may be static, i.e. they use data measured at one sampling time [START_REF] Wilkinson | Comparison of methods for wind turbine condition monitoring with SCADA data[END_REF][START_REF] Schlechtingen | Comparative analysis of neural network and regression based condition monitoring approaches for wind turbine fault detection[END_REF] or dynamic, such as ARX models [START_REF] Garlick | A model-based approach to wind turbine condition monitoring using SCADA data[END_REF][START_REF] Cross | Model-based and fuzzy logic approaches to condition monitoring of operational wind turbines[END_REF]. They can be simple linear regression models or more complex artificial neural networks [START_REF] Garcia | SIMAP: Intelligent System for Predictive Maintenance: Application to the health condition monitoring of a windturbine gearbox[END_REF][START_REF] Zaher | A multi-agent fault detection system for wind turbine defect recognition and diagnosis[END_REF][START_REF] Zaher | Online wind turbine fault detection through automated SCADA data analysis[END_REF][START_REF] Li | Normal Behavior Models for the Condition Assessment of Wind Turbine Generator Systems[END_REF][START_REF] Zhang | Comparison of data-driven and model-based methodologies of wind turbine fault detection with SCADA data[END_REF]. Let note at this point that several published works follow a classical model-based fault detection and isolation approach (model-based FDI) to build fault indicators sensitive to faults, but robust to disturbances, for wind turbine monitoring, [START_REF] Shen Yin | Data-driven design of robust fault detection system for wind turbines[END_REF][START_REF] Chen | Observer-based fdi schemes for wind turbine benchmark[END_REF][START_REF] Wang | Fault diagnosis of wind turbine blade based on robust residual error design[END_REF][START_REF] Liu | Robust fault detection for wind turbine systems[END_REF][START_REF] Borchersen | Model-based fault detection for generator cooling system in wind turbines using SCADA data[END_REF]. In these works, the residuals are built using state observers or parity equations, [START_REF] Iserman | Fault Diagnosis Systems -An Introduction from Fault Detection to Fault Tolerance[END_REF]. The proposed methods are usually implemented on wind turbine benchmarks, emulated by differential equations models where different faults can be simulated. They use high frequency SCADA data, recorded every second. Though the results obtained are very interesting, the problem addressed in these works is very different in nature from the one addressed in the present paper where 10-minutes SCADA data are used and the methods are validated on real data.

Another approach is to model the temperature of one component as a function of the temperature of some other reference components that should evolve in a similar way, such as the temperatures of two bearings but also the bearing temperature and the stator temperature [START_REF] Schlechtingen | Comparative analysis of neural network and regression based condition monitoring approaches for wind turbine fault detection[END_REF] or the temperature of the hydraulic break and the bearing temperature [START_REF] Butler | Prognostic Algorithms for Condition Monitoring and Remaining Useful Life Estimation[END_REF]. Such models using similarity may be simpler than models using causal relations but they rely on two temperatures evolving in the same way.

Mono-turbine approaches merge variables from the same wind turbine to generate fault indicator residuals that are insensitive to changes in its operational conditions. However, these residuals remain sensitive to the variations in the external environment such as the wind direction, air humidity and so on. On the opposite, multi-turbines approaches merge variables recorded from different turbines of the same wind farm in order to reduce the influence of the environmental conditions. Indeed, turbines of the same farm are submitted to the same weather conditions so variables should evolve in a similar way, somehow. [START_REF] Astolfi | Fault prevention and diagnosis through SCADA temperature data analysis of an onshore wind farm[END_REF] compares the behavior of different turbines in the same farm using curves displaying the temperature of a drive train bearing as a function of the produced power. The temperatures and powers are measured during a period of time and the curves from all the turbines are plot on the same graph. Faults can be identified visually when one of the curves deviates from the others. [START_REF] Mclaughlin | Wind farm performance assessment in the real world[END_REF] measures the difference between the power curve given by the manufacturer and an actual power curve built using data gathered during a current period of time. The difference between the two curves can exhibit a loss in performance. The differences measured for all the turbines in the farm are compared one with each other to detect a turbine with a larger loss of performance, which can be due to a fault. [START_REF] Boucher | Lowering the cost of project using simple analysis of SCADA data -a real case example[END_REF] compares the evolution in time of the temperature of a component normalized by the external temperature for turbines from the same farm. The deviation of the temperature of one turbine from the others can be the symptom of a fault. [START_REF] Wilkinson | Comparison of methods for wind turbine condition monitoring with SCADA data[END_REF] builds residuals from the differences in the normalized temperatures and concludes that changes in operational conditions can create themselves fluctuations that are too large to allow for reliable fault detection.

The literature review shows that most authors proposed methods to build fault indicators using mono-turbine approaches and so get rid of the influence of the operating conditions on the fault indicators. Few authors adopt a multi-turbine approach that allows getting rid of the influence of the external conditions, and mostly as a visualization tool. A method to synthesize a fault indicator for each turbine in a farm by comparing the temperature measured on a turbine to a farm reference (average or median of temperatures measured on all the turbines within the farm) has been proposed in [START_REF] Lebranchu | Using SCADA data for fault detection in wind turbines: local internal model vs distance to a wind farm reference[END_REF] and used in [START_REF] Cambron | Control chart monitoring of wind turbine generators using the statistical inertia of a wind farm average[END_REF] on other types of measured SCADA variables. It has been shown that such an indicator remains sensitive to the operating conditions, which can be different from one turbine to another. Hence, no solution able to get rid of the influence of both the operating conditions and external environment has been proposed thus far. Moreover, as stated by [START_REF] Tautz-Weinert | Using SCADA data for wind turbine condition monitoring -a review[END_REF], in their extensive review on wind turbines condition monitoring using SCADA data, there is a lack of published performance metrics to properly evaluate the advantage of one method from the others in terms of false alarm, true failure prediction and normal behavior prediction. To address these issues and fill this gap, we propose in this paper a hybrid multi-level synthesis method to take benefits of both approaches -mono-and multi-turbine -and to build fault indicators combining the two approaches. At the turbine level, residuals obtained with a mono-turbine model learnt during a normal condition period are first calculated for each turbine. At the farm level, these residuals are compared to a wind farm residual reference, in a multi-turbine approach. The use of mono-turbine residuals enables the influence of the operating conditions to be reduced while the use of a wind farm reference enables the changes in the environmental conditions to be accounted for. The performance of the method proposed is evaluated and compared to methods proposed in the literature on a large data set made of SCADA data recorded on a French wind farm during four years. Objective performance evaluation metrics are also proposed to compare the methods, which aim at evaluating the cost/benefit of the methods from a production manager's point of view. Fault indicators should warn of a progressing fault early enough for a maintenance operation to be scheduled. However, they should not generate false alarms with the extra cost of sending a maintenance operator team on a remote site for no reason.

The contribution of the work presented in this paper is then twofold:

• The first and main contribution of the hybrid multi-level approach proposed in this paper lies in the online real-time comparison of the prediction residual computed on one turbine (i.e. at the turbine level, using a turbine normal behavior model) to a farm-level reference prediction residual, computed using the prediction residuals from the neighbor wind turbines. This approach is said "hybrid multi-level" because it combines a mono-turbine (or turbine-level) step with a multi-turbine (or farm-level) step. To the best of the authors' knowledge, no other existing and published method proposes to compare online in real-time the prediction residual generated on a given turbine to a farm reference synthesized from the prediction residuals generated on the neighbor turbines.

• Second, the other contributions of the work presented in this paper are (i) the definition of objective performance metrics for the fault detection consistent with the applicative requirements of a wind farm manager and (ii) using these performance metrics, the evaluation of the proposed fault detection approach on a 4 years real data base from wind farm with 6 turbines, which shows that it has better performance than a residual analysis limited at the turbine level only.

The outline of this paper is as follows. Section 2 first presents different possible methods for the generation of fault indicators and focuses on the proposed combined mono-and multi-turbine approach developed in this work. The objective performance metrics considered for the evaluation of the fault detection performance are also detailed in this section. The dataset used to analyze the performances is described in Section 3 followed by a presentation and a discussion of the results obtained.

Fault detection combining turbine-level and farm-level information

This section presents the proposed two-levels methodology to synthesize fault indicators by combining both a mono-turbine and a multi-turbine approach. Each of the two levels is described : synthesis of the fleet reference at the farm-level, and elaboration of a relevant turbine-level variable for comparison to the fleet reference. Finally, the performance evaluation metrics, which are used latter in Section 3 to compare the obtained fault detection results with different fault indicators, are introduced.

Fault indicators synthesis

Farm-level : comparison of the turbine variables to a fleet reference

In wind farms, turbines are part of a fleet : they are of the same make and are subject to the same environmental conditions (wind speed, external temperatures, . . . ). Thus, the SCADA variables monitored on turbines from the same farm should evolve in a similar way when the turbines operate under normal conditions. The rationale for the method proposed in this paper relies on this assumption. SCADA variables from different turbines are compared on line. A deviation between a variable computed from a turbine and a farm reference is an indicator of an abnormal behavior.

Its general concept is presented in Figure 1.

[Figure 1 about here.] SCADA variables or model residuals computed from SCADA variables are recorded on line, averaged over a given time period and compared to a fleet reference. The distance between the variable monitored on a wind turbine and the fleet reference serves as a fault indicator for the turbine.

Let V be a variable measured on a wind turbine or a model residual synthesized on a wind turbine. V is assumed to carry information on the turbine deterioration, and several ways to choose this variable V are considered in Section 2.1.2. V j (k) is the value of the variable V from turbine j at time k. Let N T be the number of turbines in the wind farm. Let W in be an analysis window of size W . For each turbine j, with j varying from 1 to N T , a fault indicator F j (k) is built at time k, as follows :

1. For each turbine j in the farm, an averaged value of the turbine-level variable V j is calculated over the N W j samples present in the analysis window W in:

V j (k) = 1 N W j k t=k-W V j (t) (1)
For the practical implementation, it has to be noticed that, for some reasons that will be explained further, some samples may be lacking and N W j may be strictly inferior to W . If the number of samples N W j is smaller than N W L , V j (k) is not computed as it would not be representative enough.

2. The fleet reference V f leet (k) is calculated using the averaged values V j (k) at the turbine level, as follows :

V f leet (k) = median (l=1,...,N T ) (V l (k)) (2)
The fleet reference is calculated using the median and not the mean, so as to make it insensitive to abnormal values generated by faulty turbines. For the practical implementation, the fleet reference is computed if the number of averaged values V j (k) computed at time k is higher than a limit number N T L (i.e. information is available from enough turbines) and the reference should be representative of the farm normal conditions as long as more than half of the turbines used to calculate the reference operate in normal conditions.

3. A fault indicator F j (k) is then calculated for each turbine j as the distance between the monitored variable V j (k) and the fleet reference V f leet (k):

F j (k) = |V j (k) -V f leet (k)| (3)
In practice, N W L is set to W 2 + 1, which means that the variable V j (k) should be available during at least half the analysis window W . N T L is set to N T 2 + 1 , which means that V j (k) should be available for more than half the turbines of the farm for the fleet reference to be computed.

Synthesis of a turbine-level variable for comparison to the fleet reference

At the turbine level, the variable V used to create the fault indicator may be of two kinds : it can be either a residual generated as the difference between a measured value of SCADA variable and its prediction by a normal behavior model, or in a more straightforward way, it can be directly a measured SCADA variable, Figure 1.

• The variable V may be a residual, computed as the difference between a SCADA variable and its prediction by a normal behavior model. In this paper, in the same logic as using a component temperature as a fault indicator, the difference between the measured temperature of a component and its temperature predicted by a physics-related model is used.

In the generator, the active power produced by the turbine, the rotation speed of the drive shaft, and the nacelle temperature can account for a change in a component temperature. So a linear model relating these variables to the component temperature is proposed, Equation (4)

T j (k) = a.P j (k) + b.Ω j (k) + c.T nacelle j (k) + d (4)
with T j (k) the component predicted temperature, P j (k) the active power produced, Ω j (k) the shaft rotation speed, T nacelle j (k) the nacelle temperature, recorded by the SCADA system and a, b, c and d constant parameters. To built a normal behavior model, the parameters a, b, c and d of the model are regression coefficients learnt using a classical least squares algorithm [START_REF] Iserman | Fault Diagnosis Systems -An Introduction from Fault Detection to Fault Tolerance[END_REF], on data gathered during periods when the turbine is operating under normal conditions.

The variable V is then taken equal to the temperature residual Rs taken as the difference between the actual measured temperature and its predicted value, i.e. for turbine j at time k:

V j (k) := Rs j (k) = T j (k) -T j (k) (5)
The normal behavior model (4) expresses normal temperature variations due to the turbine producing electric power. Thus, the model is valid only when the turbine is in operating conditions and, for the practical implementation, all the samples gathered when the power produced by the turbine is below a production threshold T h p , ie when the turbine is stopped, are removed from the analysis window W in.

This hybrid "mono-multi-turbine" approach consisting in generating a fault indicator by comparing the residual generated for each ("mono") turbine to a ("multi-turbine") farm reference is the original fault indicator synthesis proposed and defended in this paper for its good performance and lower sensitivity to both internal (operating conditions) and external (environmental variations) influences. The main steps for the implementation of this fault indicator synthesis are sketched in Algorithm 1.

• In a more direct way, the variable V may be a SCADA variable, relevant in itself as a fault indicator, such as a component temperature. As seen in the literature -e.g. [START_REF] Butler | Prognostic Algorithms for Condition Monitoring and Remaining Useful Life Estimation[END_REF][START_REF] Schlechtingen | Comparative analysis of neural network and regression based condition monitoring approaches for wind turbine fault detection[END_REF]-, component temperatures are commonly used as fault indicators. An overheating component can be the sign of a mechanical problem or of a cooling system failure. In this case, the variable V is taken as the monitored temperature T of a considered component, i.e. for turbine j at time k:

V j (k) := T j (k) (6)
Reading or sensors errors may affect SCADA measures. Thus, in the practical implementation of this fault indicator, any temperature value outside an acceptable range is removed from the analysis window W in.

This approach based on the direct comparison of SCADA data of each individual turbine to a farm reference has already been proposed and assessed in [START_REF] Lebranchu | Using SCADA data for fault detection in wind turbines: local internal model vs distance to a wind farm reference[END_REF] for temperature measurements, and extended to other kinds of measurements in [START_REF] Cambron | Control chart monitoring of wind turbine generators using the statistical inertia of a wind farm average[END_REF] : it has been proved too sensitive to the influence of possibly different operating conditions between each individual turbine.It is considered in this paper for comparison purposes with the proposed combined mono-and multi-turbine method. For each sampling period k,

Step 1 : calculation of the averaged mono-turbine residuals

For each turbine j, with j=1 from Nt Calculate the predicted temperature T j (k) -Equation ( 5)

Calculate the corresponding residual Rs j (k) -Equation ( 6)

Calculate the averaged residual Rs j (k) on the analysis window Win -Equation (1)

EndFor

Step 2 : calculation of the farm reference for residuals comparison

Calculate the farm reference Rs f leet (k) -Equation (2)

Step The occurence of a fault on wind turbine j leads to a change in the corresponding fault indicator

F j .
Ideally, a fault indicator should be sensitive to a fault occurence, but insensitive to the variations of the farm environment and of the wind turbine operating conditions : building such a fault indicator was the motivation to propose a combined "mono-multi-turbine" approach. Under these assumptions, the statistical properties, in particular the expected value, of the fault indicator F j changes with a fault occurrence, and the fault detection can be performed using statistical hypothesis testing tools (eg testing hypothesis

H 0 : E(F j ) = µ 0j vs H 1 : E(F j ) = µ 0j
), and can be implemented by setting a threshold on F j . In theory, in this classical setting, when no information is available on the properties of the fault indicator F j in presence of a fault, the value of this threshold can be determined using the estimated or assumed statistical properties of F j when no fault is present in order to guarantee, for example, a false alarm rate. In our approach, we want to avoid any additional hypothesis on the distribution of F j under the "no fault " assumption, which could be different for each turbine and could possibly lead to a different threshold for each turbine.

Consequently, in the following, our approach is to set the detection threshold for the whole wind farm, on the basis of the whole available historical 4-years database for the farm in order to get a given false alarm level.

Performance evaluation indicators

Fault indicators are to be used in wind farms to assist maintenance operators. They should allow the detection of faults occurring on a wind turbine without generating too many false alarms.

However, in the literature on wind turbine fault indicators, no paper addresses the issue of the performance evaluation of the indicators proposed. Most of them merely report the ability of their indicator to detect a particular fault. In this paper, we chose to evaluate the fault indicators proposed from a wind farm production manager's point of view. Fault indicators should definitely be able to detect faults on wind turbines, to avoid a major degradation and then a costly repair.

However, they should be able to do it early enough for a maintenance operation to be scheduled.

Another major issue is the cost of false alarms. A false alarm results in a useless maintenance action: a maintenance team has to be sent to the production site to inspect the supposedly faulty turbine. Because wind turbines are usually located in remote sites, the inspection takes time and consequently costs a lot of money. So, in this paper, we propose 3 performance evaluation indicators, able to make a fair and objective assessment of the indicator usefulness: the detection time before failure, the useless maintenance action number and the indicator persistence, which gives to the operator some indications about the relevance and performance of the fault detection process based on the proposed indicator.

Performance evaluation indicators are computed from a data base gathered on a wind farm during a period of length D, where faults occurred. The data base is previously split into normal conditions periods and faulty periods. A normal conditions period is a period of time where no fault occurred on a wind turbine. A faulty period is a period of time where a fault was continuously present on a turbine. It ends when the failure occurs, ie when the turbine stops functioning.

Let T h be a detection threshold, set on the fault indicators F j (k). When the indicator exceeds T h, an alarm is raised. T h is set for the whole wind farm, ie the same value is set for all the turbines fault indicators, F j (k) .

Detection time before failure

The detection time before failure measures the time separating the first time the fault indicator exceeded the detection threshold T h during a faulty period from the failure time, ie the end of the faulty period. It gives an estimation of the time left to the maintenance team to repair the fault before the failure.

Number of useless maintenance actions

A false alarm occurs when the fault indicator exceeds the detection threshold T h during a normal conditions period. The false alarm rate is converted into an equivalent "number of useless maintenance actions" so as to better consider its cost. A "false alarm day" is a day when at least one false alarm occurred on at least one of the turbines. A "false alarm period" is a period made of consecutive false alarm days. The shorter possible length for a false alarm period is one day.

During a false alarm period, a maintenance team has to be sent once to the wind farm location.

So each false alarm period results in a useless maintenance action. The useless maintenance action number is the number of times a maintenance team has to be sent to the wind farm during the period D ; it is equal to the number of false alarm periods.

Indicator persistence

The persistence measures the percentage of time during which the fault indicator remains above the detection threshold T h during a faulty period. A persisting indicator makes the maintenance operator more confident in the occurrence of the fault than an indicator constantly being set on and off.

Results and discussion

Data base

The methods presented in Section 2 are applied on a set of real data, gathered on a French wind farm located in the south of France. The 6 wind turbines forming the farm are identical. They are of the same make, conceived to produce 2 MW, with a horizontal axis. SCADA data were recorded • Fault on two bearings : 2 bearings broke down on 2 generators, because of a lack of lubrification. One of the generators had to be replaced. These 2 episodes are named "faulty bearing WT6 (FB_WT6)" and "faulty bearing WT9 (FB_WT9)", for faulty bearing on wind turbines 6 and 9.

• Faults on two stator windings : two generators were stopped because of a short-cut on the stator winding, induced by an over-heating period which damaged the winding insulating material. The generators had to be replaced. These episodes are named "faulty winding WT9 (FW_WT9)" and "faulty winding WT11 (FW_WT11)".

• Faults on two cooling systems : a fault occurred on the cooling system of two machines, due to a fixation screw getting loose. These episodes are named "faulty cooling system WT7 (FCS_WT7)" and "faulty cooling system WT10 (FCS_WT10)".

[Figure 2 about here.]

The SCADA data used to implement and compute the fault indicators are the bearing temperatures, the stator temperatures, the nacelle temperatures, the active power produced and the shaft rotational speed on the 6 turbines.

Implementation of fault indicators

Four fault indicators are implemented and compared using the performance evaluation indicators presented in Section 2.2 calculated on the 4 years of data : the first considered fault indicator corresponds to the proposed combined mono-and multi-turbine approach, the three others are more classical and are implemented for comparison purposes. Each fault indicator is calculated every 10 minutes, using a sliding window of size W = 144 samples, ie 24 hours. N W L is thus set to 72 and N T L to 4. When relevant, the limit production threshold T h p , above which the temperature can be predicted, is set to 50 kW, mainly because production levels below 50 kW correspond to starting or shutting down transient behaviors. In the following, depending on the considered fault to monitor, T j can be either the bearing temperature or the stator temperature.

Implementation of the proposed combined mono-and multi-turbine approach

The indicator RES multi j (k) = Rs j (k) -Rs f leet (k) is computed for each turbine j, following the proposed combined mono-and multi-turbine approach, see Algorithm 1 :

• The residual Rs j (k) is generated using Equation ( 5), with T j the predicted bearing or stator temperature using the normal behavior model in Equation ( 4). The temperature variations due to changes in the power production are accounted for by the model. The model is learnt for each turbine, during a normal conditions period of two months at the beginning of the 4 years recording. Whenever a fault occurs on a turbine, the turbine model is learnt just after the necessary repair is made, on a fixed period of 2 months.

• This model residual Rs j (k) is then averaged to obtained Rs j (k), using Equations ( 1) with

V j = Rs j ;
• The average residual Rs j (k) is finally compared to the fleet reference Rs f leet (k) computed using Equation ( 2) with the averaged residuals from all the turbines, to obtain the fault indicator RES multi j (k).

Implementation of classical fault detection indicators for comparison

For each turbine j, three other fault indicators are considered : two of them are generated using a mono-turbine approach, the third one is generated by a simple multi-turbine comparison of the raw measurements (no residual generation at the turbine level).

• The two indicators T mono j (directly based on raw measurements) and RES mono j (based on residuals) are implemented following a classical mono-turbine approach, ie using SCADA data recorded from the turbine j alone :

-T mono j (k) = T j (k) = 1 N W j k t=k-W T j (t)
where T j can be either the bearing temperature or the stator temperature.

In this case, the SCADA temperature is averaged over 24 hours and thresholded, in a way very similar to the traditional SCADA detection system.

-RES mono j (k) = Rs j (k) = 1 N W j k t=k-W Rs j (t) with Rs j (k) = T j (k) -T j (k)
where T j can be either the bearing temperature or the stator temperature.

In this case, the residual generated using Equation ( 5) and the model in Equation ( 4) is averaged and thresholded. The temperature variations due to changes in the power production are accounted for by the model. The model is learnt for each turbine, during a normal conditions period of two months at the beginning of the 4 years recording.

Whenever a fault occurs on a turbine, the turbine model is learnt just after the necessary repair is made, on a fixed period of 2 months.

• An indicator T multi j (directly based on measurements) is implemented using a multi-turbine approach, using a fleet reference 1) and ( 2) with V j = T j , where T j is the bearing temperature or the stator temperature. The component temperature is averaged and compared to the fleet averaged temperature reference.

-T multi j (k) = T j (k) -T f leet (k), using Equations (

Performance results

[Figure 3 the energy of the indicator during a faulty period on the energy of the indicator during a normal condition period) is rather low, which makes the detection of the fault with a minimum amount of false alarms using T mono j quite difficult. The evolution of RES mono j is rather similar. Since it is a residual, its evolution varies around 0. Its standard deviation is reduced compared to T mono j since the variations in temperature due to the change in power production are accounted for by the model and thus removed but its mean remains affected by seasonal temperature variations.

Abnormal changes due to faults are more visible on the indicators but the fault detection using a fixed detection threshold seems still difficult to do.

The use of a multi-turbine approach makes it possible to remove the external temperature variations (seasonal effects) affecting the mono turbine indicators since the temperature variations are also present in the fleet reference and consequently substracted from the residuals. The faulty indicator to normal indicator contrast ratio seems to be increased and it is all the more so when RES multi j is used, since the residual is rid of variations due to external temperature changes and variations due to changes in power production. Whenever a fault occurs on a turbine, it is preceded by an increase in the corresponding turbine residual, which clearly stands out from the other residuals. A fault detection strategy using a fixed detection threshold seems much more feasible.

To compare the four indicators performances, two performance curves plotting the persistence or the detection time before failure as a function of the useless maintenance action numbers are built for each of the 6 faults. The curves are similar to traditional ROC curves, [START_REF] Van Trees | Detection, Estimation, and Modulation Theory, Part I[END_REF]. For a given value of the detection threshold, the detection time before failure (advance detection time), the persistence and the useless maintenance action number are calculated and form a point in the performance curves.

[Figure 4 about here.]

[Figure 5 about here.]

The two performances curves are displayed in Figures [START_REF] Qiao | A Survey on Wind Turbine Condition Monitoring and Fault Diagnosis -Part I: Components and Subsystems[END_REF] and [START_REF] Qiao | A Survey on Wind Turbine Condition Monitoring and Fault Diagnosis -Part II: Signals and Signal Processing Methods[END_REF]. The performance curves of multi-turbine indicators based detection are displayed in plain lines and mono-turbine ones in dotted lines. One can see that for the faults FB_WT6, FB_WT9, FW_WT11 and FW_WT9, the detection implemented with multi-turbine indicators outperforms the detection based on the mono-turbine ones, as could be expected from Figure 3. Using these multi-turbine indicators allows detecting the faults earlier with a reduced number of useless maintenance actions. If the detection threshold is set at its highest value, ie set to obtain zero useless maintenance actions, with RES multi j it is still possible to detect all the faults at least 500 minutes ahead of time ie about 4 days. If the threshold is decreased to allow 3 useless maintenance actions, the faults can be detected as early as 1200 minutes ahead of time ie about 8 days, which is sufficient to plan a maintenance intervention.

When mono-turbine indicators are used, FB_WT6 and FW_WT11, two major faults, cannot be detected. The detection time before failure remains close to zero whatever the value of the detection threshold. FB_WT9 and FW_WT9 can be detected two weeks ahead of time, with RES mono or T mono , but at the very high cost of more than 15 useless maintenance actions. The two faults where the mono-turbine indicators performances equal the multi-turbine indicators are the faults on the cooling systems, FCS_WT7 and FCS_WT10, which occurred suddenly and generated a sudden and large increase in the recorded temperature.

The detection time before failure as a function of the useless maintenance action number curves show that RES multi j performs slightly better than T multi j . The performances are comparable on FB_WT6, FW_WT9, FCS_WT11, FCS_WT7 and FCS_WT10 but improved with RES multi j on FB_WT9. However, the persistence as a function of the useless maintenance action numbers curves shows the superiority of RES multi j . Indeed, even if the detection time before failure is equivalent, the persistence is globally higher when RES multi j is used. This is an important criterion because the higher the persistence is, the higher the confidence the operator may have in the indicator. Based on this criterion, RES multi j clearly outperforms all the other fault indicators on the considered real data set, especially the mono-turbine ones. This is illustrated in Figure [START_REF] Yang | Wind turbine condition monitoring: technical and commercial challenges[END_REF] where a zoom in time is made on the period preceding FB_WT9. Indicators are presented in the same order as in Figure [START_REF] Tchakoua | Wind Turbine Condition Monitoring: State-of-the-Art Review, New Trends, and Future Challenges[END_REF]. The colored lines are the detection thresholds which set the number of useless maintenance actions at 10, 20 or 30. Their value is adapted to each indicator.. One can see that, though a slight increase in T mono j can be visually observed during the period preceding the fault, this increase cannot be detected with a fixed threshold. RES mono j can detect the fault but only because it exceeds the detection threshold for a very short period of time. This could be figured out from the persistence curve whose value remained close to zero. On the contrary, the multi-turbines indicators show an obvious increase, which can be easily detected. The fault symptom appears earlier in RES multi j and the indicator remains more often above the detection threshold, which results in a higher persistence.

[Figure 6 about here.]

Conclusion

In this paper, we propose a new hybrid multi-level strategy to build fault indicators to monitor turbines within a wind farm, following a combined mono-and multi-turbine approach. The approach proposed reduces the effect of changes in the operational and environmental conditions on the fault indicators. Residuals computed from a mono-turbine model are compared to a fleet reference residual to form a multi-turbine fault indicator. The mono-turbine residuals are less sensitive to the internal variations because the SCADA temperatures are first processed using a model explaining the variations of the temperature in function of the operational conditions. Rises in the temperature due to increases in the power produced or in the rotational speed are accounted for by the model and removed from the residuals. The mono-turbine residuals are further post processed in a multiturbine approach, which enables the impact of the variations in the environmental conditions to be reduced. Indeed, turbines in a same wind farm are part of a fleet. They are of the same make and model, and their variables evolve in a similar way since they are submitted to the same environmental conditions. New performance evaluation criteria were proposed to analyze the indicator performances for fault detection purpose, keeping in mind the cost of any useless maintenance action. The performance of the combined mono-and multi-turbine fault indicators were compared to more traditional mono-turbine indicators on a data base made of SCADA data recorded every ten minutes from a real wind farm on a four year period. As far as we know, this is the first time a paper presents such an extensive analysis of fault indicators performance for fault detection, on real data with a wind farm production manager point of view. The results clearly showed that the effects of changes in both environmental and operational conditions on the combined mono-and multi-turbine residuals were reduced during the fault-free periods. This results in the combined mono-and multi-turbine fault indicators being able to detect all the faults occurring on the farm with a limited number of useless maintenance interventions and with an increased persistence.

List of Figures

Multi-turbine approach : generic principle for the synthesis of fault indicator F j for wind turbine j. The choice of the variable V as a mono-turbine residual corresponds to the combined mono-and multi-turbine method proposed in this paper. . . . . . .

Algorihm 1 EndForPhase 2 -

 12 General algorithm for the combined mono-and multi-turbine fault indicator generation Phase 1 -Offline learning phase For each turbine j, with j=1 from Nt Learn (by a least squares algorithm) the model parameters a, b, c d of the normal behavior model in Equation(5) Online implementation phase

3 : 4 :EndFor 2 . 1 . 3 .

 34213 calculation of the fault indicator For each turbine j, with j=1 from Nt Calculate the multi-turbine fault indicator -Equation (3) EndFor Step Alarm set off For each turbine j, with j=1 from Nt If the fault indicator is higher than the detection threshold Set an alarm off for the turbine i EndIf EndFor Fault detection based on the generated fault indicators

  every 10 minutes during 4 years, from November 2009 till December 2013. During this 4 years, 6 single faults affecting the generator occurred on different turbines at different times, some of them generating a major failure and, consequently, the machine shutdown for several weeks, see Figure(2) :

  Figure (3) shows the evolution of the 4 fault indicators during the four years when the SCADA data were recorded on the 6 turbines of the wind farm. The upper figures display the evolution of

24 2Figure 2 :Figure 3 :Figure 4 :

 24234 Figure 2: Occurrence of failures and "failure periods" on the 6 wind turbines of the considered farm during the 4 years recorded database -Bearing failures in red, cooling system failures in green and windings failures in blue

Figure 5 :

 5 Figure 5: Performance characteristics curves : fault indicator persistence as a function of the number of useless maintenance interventions for 6 different faults : FB_WT6, FB_WT9, FCS_WT7, FW_WT11, FCS_WT10, FW_WT9

Figure 6 :

 6 Figure 6: Evolution of the 4 indicators over the period preceding FB_WT9

Figure 1: Multi-turbine approach : generic principle for the synthesis of fault indicator Fj for wind turbine j. The choice of the variable V as a mono-turbine residual corresponds to the combined mono-and multi-turbine method proposed in this paper.


Abderrahim Halimi

Paul Honeine

Jose Bioucas-Dias

ROBUST HYPERSPECTRAL UNMIXING ACCOUNTING FOR RESIDUAL COMPONENTS

Keywords: Hyperspectral imagery, robust unmixing, Bayesian estimation, coordinate descent algorithm, Gaussian process, gamma Markov random field

This paper presents a new hyperspectral mixture model jointly with a Bayesian algorithm for supervised hyperspectral unmixing. Based on the residual component analysis model, the proposed formulation assumes the linear model to be corrupted by an additive term that accounts for mismodelling effects (ME). The ME formulation takes into account the effect of outliers, the propagated errors in the signal processing chain and copes with some types of endmember variability (EV) or nonlinearity (NL). The known constraints on the model parameters are modeled via suitable priors. The resulting posterior distribution is optimized using a coordinate descent algorithm which allows us to compute the maximum a posteriori estimator of the unknown model parameters. The proposed model and estimation algorithm are validated on both synthetic and real images showing competitive results regarding the quality of the inferences and the computational complexity when compared to the state-of-the-art algorithms.

INTRODUCTION

Spectral unmixing (SU) is a source separation problem consisting of recovering the spectral signatures (endmembers) of the materials present in the scene, and quantifying their proportions (abundances) within each hyperspectral image pixel. The linear mixture model (LMM) is the widely used model for SU mainly because of its simplicity. However, this model can be inappropriate for some hyperspectral scenarios, namely in presence of nonlinearity (such as multiple scattering or intimate mixtures), endmember variability, or outliers. In addition, the estimated LMM abundances can be badly affected by the propagated errors in the signal processing chain such as a bad estimation of the endmembers or their number (especially under a supervised unmixing, i.e., estimating only the abundances). These effects emphasize the need for a robust hyperspectral unmixing strategy to deal with these mismodelling effects [START_REF] Altmann | Robust linear spectral unmixing using anomaly detection[END_REF][START_REF] Févotte | Nonlinear hyperspectral unmixing with robust nonnegative matrix factorization[END_REF].

This paper introduces a residual component (RC) mixture model for robust hyperspectral unmixing, the model generalizes the wellknown LMM by accounting for the illumination variation and the presence of a residual term. This term aims at reducing the effect of smooth spatial/spectral outliers that can be due to the presence of NL [START_REF] Heylen | A review of nonlinear hyperspectral unmixing methods[END_REF][START_REF] Dobigeon | Nonlinear unmixing of hyperspectral images: Models and algorithms[END_REF], EV [START_REF] Zare | Endmember variability in hyperspectral analysis: Addressing spectral variability during spectral unmixing[END_REF] or to signal processing chain errors. Estimating the abundances associated with this mixture model is a challenging problem. We propose here a hierarchical Bayesian model to estimate the parameters and hyperparameters associated with the RC model. This hierarchical model introduces prior distributions that enforce known physical constraints on the estimated parameters such as the sum-to-one and positivity of the abundances, and the smooth spectral behavior of the ME. Moreover, the spatial correlation of the residual term has been introduced by considering Markov random fields (MRF) [START_REF] Dikmen | Gamma markov random fields for audio source modeling[END_REF]. Using the likelihood and the considered prior distributions, the joint posterior distribution of the unknown parameter vector is then derived. The minimum mean square error (MMSE) and maximum a posteriori (MAP) estimators of these parameters cannot be easily computed from the obtained joint posteriors. In this paper, the MAP estimator is evaluated by considering a coordinate descent algorithm (CDA) [START_REF] Bertsekas | Nonlinear programming[END_REF][START_REF] Sigurdsson | Hyperspectral unmixing with lq regularization[END_REF][START_REF] Halimi | Bayesian estimation of smooth altimetric parameters: Application to conventional and delay/Doppler altimetry[END_REF] that sequentially updates the abundances, the noise variances and the residual term. This provides a reduced computational cost when compared to the Markov chain Monte-Carlo based algorithms [START_REF] Altmann | Bayesian nonlinear hyperspectral unmixing with spatial residual component analysis[END_REF][START_REF] Altmann | Residual component analysis of hyperspectral images: Application to joint nonlinear unmixing and nonlinearity detection[END_REF][START_REF] Halimi | Unsupervised unmixing of hyperspectral images accounting for endmember variability[END_REF]. The proposed Bayesian model and estimation algorithm are validated using synthetic and real hyperspectral images. The obtained results are very promising and show the potential of the proposed mixture and Bayesian models and their associated inference algorithm.

The paper is structured as follows. Section 2 introduces the proposed mixture model to deal with the ME. The proposed hierarchical Bayesian model and its estimation algorithm are introduced in Sections 3 and 4. Section 5 is devoted to testing and validating the proposed techniques using synthetic images with known ground truth. Section 6 shows results obtained using a real hyperspectral image. Conclusions and future work are finally reported in Section 7.

MIXTURE MODEL

The LMM is widely used because of its simplicity. However, there are a lot of situations where the linear model is not valid because of the presence of EV, NL or other ME (due to physical outliers or signal processing errors). This paper proposes a new formulation based on a residual component model [START_REF] Kalaitzis | Residual components analysis[END_REF] that is expressed as the sum of a linear model and a residual term. The general observation model for the (L × 1) pixel spectrum y i,j is given by y i,j = Si,jai,j + φ i,j (Si,j, ai,j) + ei,j,

where ai,j = (a1,i,j, • • • , aR,i,j) T is an (R × 1) vector of abundances associated with the pixel (i, j) and satisfying the positivity and sum-to-one constraints (ar,i,j ≥ 0, ∀r ∈ {1, . . . , R} and R r=1 ar,i,j = 1), R is the number of endmembers, ei,j ∼ N (0L, Σ) is an additive centered Gaussian noise with a diagonal covariance matrix

Σ = diag σ 2 , σ 2 = σ 2 1 , • • • , σ 2 L
T is an (L × 1) vector containing the noise variances, L is the number of spectral bands, Si,j(M ) = Si,j is the endmember matrix that depends on each pixel to introduce EV effect, M is a fixed endmember matrix that is assumed known (extracted using an endmember extraction algorithm) and φ i,j (Si,j, ai,j) is a residual term that might depend on the endmembers or the abundances to model NL or EV effects. This paper accounts for mismodelling effects or the presence of outliers by considering a residual term φ M E i,j = di,j that shows spatial and spectral correlations. Moreover, it accounts for the illumination variation effect by considering Si,j = ci,jM . Under these considerations, the observation model (1) reduces to

y i,j = ci,jM ai,j + di,j + ei,j, (2)
where di,j is an L × 1 smooth spectral vector. This term gathers the shape EV effects that can be approximated by a smooth function as already highlighted in [START_REF] Halimi | Hyperspectral unmixing in presence of endmember variability, nonlinearity or mismodelling effects[END_REF]. It also gathers the multiple scattering NL effects when the endmembers are themselves smooth. We refer the reader to [START_REF] Halimi | Hyperspectral unmixing in presence of endmember variability, nonlinearity or mismodelling effects[END_REF] for more details about the relation between ( 2) and other NL and EV models, and the possible interpretations of di,j.

Note that model ( 2) reduces to the LMM when di,j = 0L, and ci,j = 1, ∀i, j. Note also that other models have been introduced in the literature to account for the effect of outliers such as [START_REF] Altmann | Robust linear spectral unmixing using anomaly detection[END_REF] which proposed spatial/spectral correlated outliers by considering discrete MRF and [START_REF] Févotte | Nonlinear hyperspectral unmixing with robust nonnegative matrix factorization[END_REF] which proposed positive-sparse outliers that have no spatial or spectral correlation. In contrast to [START_REF] Altmann | Bayesian nonlinear hyperspectral unmixing with spatial residual component analysis[END_REF][START_REF] Altmann | Residual component analysis of hyperspectral images: Application to joint nonlinear unmixing and nonlinearity detection[END_REF] that deal only with NL effect, and [START_REF] Halimi | Unsupervised unmixing of hyperspectral images accounting for endmember variability[END_REF] that deals with EV, model (2) considers both effects at the same time.

HIERARCHICAL BAYESIAN MODEL

This section introduces a hierarchical Bayesian model to estimate the unknown parameter vector Θ = A, c, D, σ 2 , , where A (resp. D) gathers all the abundance vectors (resp. residual terms), and is a hyperparameter described later in the text.

Likelihood

Using the observation model ( 1), the Gaussian properties of the noise sequence ei,j, and exploiting independence between the observations in different spectral bands, yield the following Gaussian distribution for the likelihood

y i,j |ai,j, ci,j, di,j, Σ ∼ N µ i,j , Σ . (3)
where µ i,j = ci,jM ai,j + di,j and ∼ means "is distributed according to". The independence between the observed pixels leads to

f (Y |Θ) = i j f (y i,j |ai,j, ci,j, di,j, Σ)
, where Y gather all the pixels.

Parameter priors

This section introduces the prior distributions that we have chosen for the parameters of interest A, c, D and Σ.

Abundance matrix A

Due to the positivity and sum-to-one (PSTO) constraints, the abundance vector should live in the following simplex S = ai,j ar,i,j ≥ 0, ∀r and

R r=1 ar,i,j = 1 . ( 4
)
Since there is no additional information about this parameter vector, we propose to assign a uniform prior in the simplex S to the abundances [START_REF] Dobigeon | Semisupervised linear spectral unmixing using a hierarchical Bayesian model for hyperspectral imagery[END_REF][START_REF] Halimi | Nonlinear unmixing of hyperspectral images using a generalized bilinear model[END_REF].

Prior for c

The variation in illumination is introduced by the parameter ci,j that is pixel-dependent. In many works, this parameter has been fixed to the value #1, which represents the absence of illumination variability [START_REF] Halimi | Nonlinear unmixing of hyperspectral images using a generalized bilinear model[END_REF][START_REF] Heinz | Fully constrained least-squares linear spectral mixture analysis method for material quantification in hyperspectral imagery[END_REF][START_REF] Bioucas-Dias | Alternating direction algorithms for constrained sparse regression: Application to hyperspectral unmixing[END_REF]. In others the sum-to-one (STO) was relaxed leading to an unconstrained parameter ci,j (free in R). In this paper, we relax the STO constraint but include some prior knowledge on ci,j. Indeed, we allow this parameter to fluctuate around #1 by considering a conjugate Gaussian prior as follows

ci,j ∼ N 1, η 2 , ( 5
)
where η 2 is a small fixed variance (η 2 = 0.01 in the rest of the paper). For simplicity, we denote "x|θ ∼ ...", by "x ∼ ..." when the parameter θ is a user fixed parameter. Note finally that the joint prior of c is obtained by assuming a priori independence between the coefficients ci,j, as follows f (c) = i,j f (ci,j).

Mismodelling coefficients di,j

The smooth property of the vector di,j is introduced by considering a conjugate Gaussian prior as follows

di,j| 2 i,j ∼ N 0L, 2 i,j H , ( 6
)
where H is the squared-exponential covariance function given by

H( , ) = exp - ( -) 2 (L/2) 2
, which introduces the spectral smoothness on di,j. A spatial correlation is also introduced by enforcing a smooth variation of the residual energies d T i,j H -1 di,j . This is achieved by considering a specific prior for 2 i,j , as explained in Section 3.3. Finally the joint prior of D is obtained by assuming a priori independence between the mismodelling coefficients, i.e., f (D| ) = i,j f di,j| 2 i,j .

Noise variances

The noise variances are assigned a conjugate inverse gamma distribution given by

f σ 2 = L =1 f σ 2 , with σ 2 ∼ IG (ϕ , ψ ) , (7)
where σ 2 are assumed a priori independent. The hyperparameters ϕ and ψ are fixed to approximate the HySime estimated variances [START_REF] Bioucas-Dias | Hyperspectral subspace identification[END_REF]. Note finally that the parameters can also be set to ϕ = ψ = 0 in absence of prior knowledge about σ 2 , leading to a noninformative Jeffreys' prior.

Hyperparameter priors

Due to the spatial organization of hyperspectral images, we expect the energies of the mismodelling coefficients di,j to vary smoothly from one pixel to another. This behavior is obtained by introducing an auxiliary variable w (of size Nrow × Ncol) and assigning a gamma Markov random field prior for ( , w) (see [START_REF] Dikmen | Gamma markov random fields for audio source modeling[END_REF][START_REF] Altmann | Bayesian nonlinear hyperspectral unmixing with spatial residual component analysis[END_REF] for more details regarding this prior). This prior ensures that each 2 i,j is connected to four neighbor elements of w and vice-versa (see Fig. 1). The energies 2 i,j are conditionally independent and the 1st order neighbors (i.e., the spatial correlation) is only introduced via the auxiliary variables w. An interesting property of this joint prior is that the i+1,j

2 i+1,j+1 2 i+1,j-1 2 i,j 2 i-1,j w 2 i+1,j+1 w 2 i+1,j 2 i,j+1 2 i,j-1 w 2 i,j w 2 i,j+1 2 i-1,j+1 2 i-1,j-1 ζ ζ ζ ζ ζ ζ ζ ζ ζ ζ ζ ζ ζ ζ ζ ζ
Fig. 1. Proposed 1st order gamma-MRF neighborhood structure. conditional prior distributions of and w reduce to conjugate inverse gamma (IG) and gamma (G) distributions as follows

2 i,j |w, ζ ∼ IG (4ζ, 4ζρ1,i,j(w)) , w 2 i,j | , ζ ∼ G (4ζ, 1/(4ζρ2,i,j( ))) , (8)
where

ρ1,i,j(w) = (w 2 i,j + w 2 i+1,j + w 2 i,j+1 + w 2 i+1,j+1 )/4, ρ2,i,j( ) = ( -2 i,j + -2 i-1,j + -2 i,j-1 + -2 i-1,j-1 )/4, (9)
and ζ is a fixed coupling parameter that controls the amount of spatial smoothness enforced by the gamma-MRF.

Posterior distributions

The joint posterior distribution of the hierarchical Bayesian model can be computed from the following Bayes' rule

f (Θ|Y ) ∝ f (Y |Θ)f (A) f σ 2 f (c) f (D| ) f ( , w) , (10)
where ∝ means "proportional to" and we have assumed a priori independence between the parameters of each model. The MMSE and MAP estimators associated with the posterior [START_REF] Altmann | Bayesian nonlinear hyperspectral unmixing with spatial residual component analysis[END_REF] are not easy to determine. In this paper, and akin to [START_REF] Halimi | Bayesian estimation of smooth altimetric parameters: Application to conventional and delay/Doppler altimetry[END_REF], we propose to evaluate the MAP estimator by using an optimization technique maximizing the posterior (10) w.r.t. the parameters of interest.

COORDINATE DESCENT ALGORITHM

Because of the large number of parameters to estimate, we propose a coordinate descent algorithm [START_REF] Bertsekas | Nonlinear programming[END_REF][START_REF] Sigurdsson | Hyperspectral unmixing with lq regularization[END_REF][START_REF] Halimi | Bayesian estimation of smooth altimetric parameters: Application to conventional and delay/Doppler altimetry[END_REF] that sequentially updates the different parameters. In each step of the proposed CDA-ME algorithm, the posterior distribution is maximized w.r.t. one parameter, the other being fixed. Thus, the algorithm iteratively updates each parameter by maximizing its conditional distribution as follows (see [START_REF] Halimi | Hyperspectral unmixing in presence of endmember variability, nonlinearity or mismodelling effects[END_REF] for more details regarding the conditional distributions):

• Conditional of A: truncated Gaussian distribution (whose one maximum is obtained with SUNSAL-FCLS 1 [START_REF] Bioucas-Dias | Alternating direction algorithms for constrained sparse regression: Application to hyperspectral unmixing[END_REF]) 1 SUNSAL-FCLS satisfies the PSTO constraints while SUNSAL-CLS only ensure the positivity constraint.

• Conditional of D: Gaussian distribution (analytical expression of the mean)

• Conditional of : inverse gamma distribution (analytical expression of the maximum)

• Conditional of w : gamma distribution (analytical expression of the maximum)

• Conditional of σ 2 : inverse gamma distribution (analytical expression of the maximum)

• Conditional of c: Gaussian distribution (analytical expression of the mean)

Regarding the sequence generated by the coordinate descent algorithm, the proposition 2.7.1 in [START_REF] Bertsekas | Nonlinear programming[END_REF] asserts that its limit points are stationary points of [START_REF] Altmann | Bayesian nonlinear hyperspectral unmixing with spatial residual component analysis[END_REF] provided that the maximum of that function w.r.t. Θ along each coordinate is unique. This is easily checked for all the parameters since they have unimode conditional distributions (Gaussian, gamma and inverse gamma distributions). Note however that the cost function is not convex, thus, the solution obtained might depend on the initial values that need to be chosen carefully. Therefore, the abundances A are initialized with SUNSAL-FCLS [START_REF] Bioucas-Dias | Alternating direction algorithms for constrained sparse regression: Application to hyperspectral unmixing[END_REF], the residual terms are initialized by 0, the noise variance is initialized by HySime [START_REF] Bioucas-Dias | Hyperspectral subspace identification[END_REF], the illumination coefficient c is initialized by considering the sum of the abundances that are estimated using only the positivity constraint with SUNSAL-CLS [START_REF] Bioucas-Dias | Alternating direction algorithms for constrained sparse regression: Application to hyperspectral unmixing[END_REF]. With these initializations, the proposed algorithm reached minima of "good quality" in the considered simulations (see Sections 5 and 6). Thus, the CDA algorithm constitutes a good balance between computational efficiency (obtained by solving the simple problems associated with each descent step) and the quality of the solution (experimentally observed when considering a good initialization).

SIMULATION RESULTS ON SYNTHETIC DATA

Four synthetic images of size 100 × 100 pixels and L = 207 spectral bands were generated using R = 3 endmembers extracted from the ENVI software library [START_REF]ENVI User's guide Version 4.0[END_REF]. All images were corrupted by i.i.d. Gaussian noise for a fair comparison with SU algorithms using this assumption. The images were generated using different mixture models. The first image I1 is composed of LMM pixels. The second image I2 is composed of a mixture of the LMM and 3 nonlinear models (K = 4 spatial classes were generated using a Potts-MRF with granularity parameter β = 0.8). The third image I3 is partitioned into 4 spatial classes (the same partition map as for I2). In each class, the pixels were generated according to the LMM while using a different set of spectra associated with the same 3 materials, to study the EV effect. The fourth image I4 is partitioned into K = 2 classes using a Potts-MRF. The pixels of the first class were generated according to the RCA-ME model ( 2) with 2 = 0.002. The pixels of the second class were generated with 4 endmembers to simulate the effect of a bad estimation of the number of endmembers. The abundances are uniformly distributed in the simplex S for I1. However, to simulate a highly mixed scenario, the abundances of I2, I3 and I4 are clustered inside the simplex using a Dirichlet distribution whose parameters are selected randomly in the interval [START_REF] Altmann | Robust linear spectral unmixing using anomaly detection[END_REF][START_REF]ENVI User's guide Version 4.0[END_REF]. The illumination coefficient c varies in [0.9, 1.15] for all images. Table 1 compares the proposed CDA-ME algorithm with the FCLS [START_REF] Heinz | Fully constrained least-squares linear spectral mixture analysis method for material quantification in hyperspectral imagery[END_REF], SUNSAL-CLS [START_REF] Bioucas-Dias | Alternating direction algorithms for constrained sparse regression: Application to hyperspectral unmixing[END_REF], SKHype [START_REF] Chen | Nonlinear unmixing of hyperspectral data based on a linear-mixture/nonlinearfluctuation model[END_REF], and RCA-MCMC [START_REF] Altmann | Residual component analysis of hyperspectral images: Application to joint nonlinear unmixing and nonlinearity detection[END_REF] algorithms. Note that SUNSAL-CLS approximates the variation of illumination by relaxing the sum-to-one constraint. The criteria used are the distance between the actual and estimated abundances RMSE (A) = This section illustrates the performance of CDA-ME when applied to a real hyperspectral image. The considered image, denoted as Madonna, was acquired in 2010 by the Hyspex hyperspectral scanner over Villelongue, France (00 03'W and 4257'N). The dataset contains L = 160 spectral bands recorded from the visible to near infrared (400 -1000nm) with a spatial resolution of 0.5m [START_REF] Sheeren | Mapping ash tree colonization in an agricultural mountain landscape: Investigating the potential of hyperspectral imagery[END_REF]. It has already been studied in [START_REF] Halimi | Unsupervised unmixing of hyperspectral images accounting for endmember variability[END_REF][START_REF] Altmann | Unsupervised post-nonlinear unmixing of hyperspectral images using a Hamiltonian Monte Carlo algorithm[END_REF] and is mainly composed of forested areas (see Fig. 2 (top)). This image contains 100 × 100 pixels and is composed of R = 4 components: tree, grass, soil and shadow. The Bayesian UsLMM algorithm [START_REF] Dobigeon | Joint Bayesian endmember extraction and linear unmixing for hyperspectral imagery[END_REF] was used to estimate R = 4 endmembers. The estimated CDA-ME abundance maps are in good agreement with the state-of-the-art algorithms and are not shown here for brevity. Fig. 2 (middle-right) shows the energies of the difference between the reconstructed signal and the linear model (i.e., ||ŷ i,j -M âi,j||) . It is clear from Fig. 2 (middle-right) that the MEs are mainly located in regions with multiple physical components, and in presence of relief (near trees) where multiple scattering effect might occurs. These energies are in good agreement with the nonlinear coefficients estimated using the RCA-MCMC algorithm as shown in Fig. 2 (middle-left). However, in contrast to RCA-MCMC that interprets these energies as nonlinearities, the proposed CDA-ME highlights the presence of an illumination variation (see Fig. 2 (bottom-left)) and an additional residual term (see Fig. 2 (bottom-right)). Note that the residual term gathers NL, EV and outlier effects and is mainly located near the trees. To summarize, this section shows that the CDA-ME algorithm captures similar spatial residual effects than RCA-MCMC but gives a different interpretation. Most of these residuals appear in region of intersection between the physical elements and in presence of vegetation such as trees. CDA-ME is very flexible and can capture different physical effects such as NL, EV and outliers. Note finally that more details regarding the CDA-ME results, including the processing of another real image (the Moffett scene), are provided in [START_REF] Halimi | Hyperspectral unmixing in presence of endmember variability, nonlinearity or mismodelling effects[END_REF] and not shown here for brevity.

CONCLUSIONS

This paper introduced a hyperspectral mixture model and its associated Bayesian algorithm for robust hyperspectral unmixing. The proposed model was introduced under a general formulation that can be adapted to account for other physical effects. A hierarchical Bayesian model was proposed to introduce the known constraints on their parameters. Those parameters were estimated using a coordinate descent algorithm that showed a reduced computational cost when compared to state-of-the-art algorithms. The proposed algorithm showed good performance when processing synthetic data generated with the linear model or other more sophisticated models. Results on real data confirmed the good performance of the proposed algorithm and showed its ability to extract different features in the observed scenes. Future work includes the detection of the presence of outliers using hypothesis tests.

  arccos(•) is the inverse cosine operator and y n , ŷn denotes the #nth measured and estimated pixel spectra. CDA-ME provides the best performance for the LMM, EV and ME based images while it shows good performance in presence of NL effects. These results highlight the robustness of CDA-ME to the different effects that can affect hyperspectral images. Finally, apart from the LMM based algorithms (FCLS and SUNSAL-CLS), the proposed CDA-ME always provides the fastest results.

Fig. 2 .

 2 Fig. 2. Top: Real Madonna image. Middle-left: nonlinearities coefficients estimated using the RCA-MCMC algorithm [11]. Middleright: Square root of the energies of the difference between the reconstructed signal and the linear model for the proposed CDA-ME algorithm (||ŷ i,j -M âi,j||). Bottom-left: estimated illumination variation |1-ci,j| for the CDA-ME algorithm. Bottom-right: square root of the energies of the residual terms di,j for the CDA-ME algorithm.

Table 1 .

 1 Results on synthetic images.

				(×10 -2 )		Time
			RMSE	RE	SAM	(s)
		FCLS	7.78	3.58	6.24	1.2
	I1	SUNSAL-CLS	3.42	2.27	5.62	0.1
	(LMM,	SKHype RCA-MCMC	1.41 4.31	--	--	541 6737
	K = 1)	CDA-ME	1.35	2.27	5.62	88
		FCLS	24.76	15.74 10.64	1.7
	I2	SUNSAL-CLS	16.55	4.17	7.57	0.07
	(LMM, NL	SKHype RCA-MCMC	5.87 5.66	--	--	547 9009
	K = 4)	CDA-ME	6.61	2.89	6.17	37
		FCLS	10.22	2.91	6.0	1.3
	I3	SUNSAL-CLS	10.16	2.53	5.93	0.14
	(EV,	SKHype RCA-MCMC	7.69 13.33	--	--	623 11025
	K = 4)	CDA-ME	4.29	2.35	5.51	316
		FCLS	12.59	5.01	7.19	1.21
	I4	SUNSAL-CLS	12.65	2.91	6.94	0.12
	(RCA-ME,	SKHype RCA-MCMC	11.68 18.20	--	--	374 6029
	K = 2)	CDA-ME	6.07	2.29	5.56	315
		6. RESULTS ON REAL DATA

This work was supported in part by the HYPANEMA ANR Project under Grant ANR-12-BS03-003, and in part by the Portuguese Fundacão para a Ciência e Tecnologia (FCT), grant UID/EEA/5008/2013