Student Member, IEEE, J. Chanussot A Villa

Senior Member, IEEE J A Benediktsson

Fellow, IEEE C Jutten

Spectral Unmixing for the Classification of Hyperspectral Images at a Finer Spatial Resolution

Keywords: Spatial resolution improvement, Hyperspectral data, Source separation, Simulated annealing, Spatial regularization

The problem of classification of hyperspectral images containing mixed pixels is addressed. Hyperspectral imaging is a continuously growing area of remote sensing applications. The wide spectral range of such imagery, providing a very high spectral resolution, allows to detect and classify surfaces and chemical elements of the observed image. The main problem of hyperspectral data is the (relatively) low spatial resolution, which can vary from a few to tens of meters. Many factors make the spatial resolution one of the most expensive and hardest to improve in imaging systems. For classification, the major problem caused by low spatial resolution are the mixed pixels, i.e., parts of the image where more than one land cover map lie in the same pixel. In this work we propose a method to address the problem of mixed pixels and to obtain a finer spatial resolution of the land cover classification maps. The method exploits the advantages of both soft classification techniques and spectral unmixing algorithms, in order to determine the fractional abundances of the classes at a sub-pixel scale. Spatial regularization by Simulated Annealing is finally performed to spatially locate the obtained classes. Experiments carried out on synthetic real data sets show excellent results both from a qualitative and quantitative point of view.

I. INTRODUCTION

Land cover classification of remote sensing data is an important application of image analysis. It is used in many practical applications, such as precision agriculture, monitoring and management of natural disasters, issues related to security and defense. The continuously growing availability of hyperspectral imagery, which records hundreds of images corresponding to different wavelength channels, has opened new possibilities in the field of image analysis and classification [START_REF] Landgrebe | Signal Theory Methods in Multispectral Remote Sensing[END_REF]. Hyperspectral sensors are characterized by a very high spectral resolution and a spatial resolution which can vary from few to tens of meters. One of the major issues of hyperspectral images is that as the spectral resolution increases, the spatial resolution usually gets worse, and low spatial resolution is quite common, especially in case of high altitude sensors or instruments which cover wide areas [START_REF] Vane | The airborne visible/infrared imaging spectrometer (aviris)[END_REF]. These sensor limitations can affect the performances of algorithms used to process hyperspectral data. In classification, the relatively low spatial resolution can lead to the challenging problem of mixed pixels (Fig. 1), i.e., pixels containing more than one land cover type [START_REF] Cracknell | Synergy in remote sensing: What's in a pixel?[END_REF]. Also in case of high spatial resolution, a hyperspectral image is often a combination of pure and mixed pixels. A number of so-called full pixel techniques, based on the assumption that each pixel corresponds to the spectral signature of one predominant land cover type, have been proposed during the last decades for the classification of hyperspectral images. [START_REF] Landgrebe | Signal Theory Methods in Multispectral Remote Sensing[END_REF], [START_REF]Kernel methods for Remote Sensing Data Analysis[END_REF]. These techniques are not suitable for the analysis of mixed pixels and will inevitably lead to a high error rate when used for scenarios with a high number of sites with mixtures of land cover classes. The issue of mixed pixels has been considered in several works. A widely investigated approach is the use of soft classification techniques [START_REF] Nachtegael | Soft Computing in Image Processing[END_REF]. These classifiers do not assign a pixel to only one class, but they produce a set of images (one per class) that express for each pixel the degree of membership in the class in question [START_REF] Nachtegael | Soft Computing in Image Processing[END_REF]. However, the membership degree does not necessarily reflect the fractional abundance of a class within a mixed pixel, and the probability of a pixel to belong to one class does not necessarily correspond to the fractional part of the pixel covered by the considered class. Linear spectral mixture analysis (SMA) [START_REF] Keshava | Spectral unmixing[END_REF] is a soft classification technique explicitely designed to address this problem. Following the spectral mixing model, the spectral signature of a mixed pixel is assumed to be the weighted sum of some constituent spectra, also called endmembers. Spectral unmixing is the procedure by which the measured spectrum of a pixel is decomposed into a collection of endmembers, and a set of corresponding fractions, or abundances, that indicate the proportion of each endmember within the pixel. A number of techniques, exploiting both statistical and geometrical properties of the data, was proposed over the last few years [START_REF] Plaza | A quantitative and comparative analysis of endmember extraction algorithms from hyperspectral data[END_REF], [START_REF] Plaza | Optical Remote Sensing -Advances in Signal Processing and Exploitation Techniques[END_REF]. These techniques can partially overcome the weakness of full pixel methods when analysing mixed pixels. However, when used to obtain crisp classification maps, the endmembers selection and the abundances determination are negatively affected by spectral variability [START_REF] Li | Wavelet-based feature extraction for improved endmember abundance estimation in linear unmixing of hyperspectral signals[END_REF], and common hard classification methods are more suitable in such a case [START_REF] Wang | Integration of soft and hard classifications using extended support vector machines[END_REF]. Illustration of the problem of mixed pixels: When the spatial resolution is not fine enough, several land cover classes lie in the same pixel. In this case, a hard classification process cannot give an accurate information about the pixel coverage, leading inevitably to a loss of information be seen as complementary methods, since the former are more suitable for the classification of pixels dominated by a single land cover class, while the latter are devoted to the mixed pixels analysis. Since a hyperspectral image usually contains many areas with pure pixels and some others with mixed, the combination of these two techniques can be seen as an interesting approach for the analysis of hyperspectral data. In spite of the fact that full pixel techniques and spectral unmixing methods could be combined in order to obtain improved classification maps in terms of accuracy, only few attempts have been made so far in order to jointly use these techniques, especially in [START_REF] Wang | Integration of soft and hard classifications using extended support vector machines[END_REF], [START_REF] Plaza | A new approach to mixed pixel classification of hyperspectral imagery based on extended morphological profiles[END_REF]. Plaza et al. [START_REF] Plaza | A new approach to mixed pixel classification of hyperspectral imagery based on extended morphological profiles[END_REF] proposed an unsupervised approach for the classification of mixed pixels in hyperspectral imagery, based on the generalization of the concept of extended morphological profiles to case of multi bands images. The comparison of the method with commonly used techniques when applied to images with a low number of land cover classes or in case of high classes spectral separability showed interesting results. However, due to its unsupervised nature, the algorithm is expected to encounter difficulties when applied to more challenging data sets. An extension of the SVM classification technique to address the problem of mixed pixels was recently proposed in [START_REF] Wang | Integration of soft and hard classifications using extended support vector machines[END_REF]. This extension provides interesting results when applied to synthetic data.

Hard classification and spectral unmixing techniques can

In this paper, a new supervised technique, which takes advantage of both probabilistic classification and spectral unmixing mapping techniques, is proposed in order to handle the issue of mixed pixels. The concept of sub-pixel mixing is also considered, in the attempt of obtaining land cover maps with an improved spatial resolution. The idea of subpixel mapping was first presented by Atkinson in [START_REF] Atkinson | Innovations in GIS[END_REF]. He proposed to use the output of a soft classification technique in order to obtain a super-resolution mapping, trying to maximize the spatial correlation of the land cover classes to determine subpixels spatial locations. Since then, a number of techniques focused on better estimating sub-pixel fractional abundances determination and obtaining land cover maps with higher spatial resolution have been proposed [START_REF] Mertens | Sub-pixel mapping: A comparison of techniques[END_REF]. In this paper, we propose the use of Simulated Annealing (SA) for this purpose, due to its simplicity and ease of use. This method has shown good results in a number of optimization and real problems, and its wide range of parameters grants a high flexibility with respect to the analysed problem. In multi-hyperspectral remote sensing, it has successfully been used for classification [START_REF] Bardossy | Fuzzy-rule based classification of remotely sensed imagery[END_REF], [START_REF] Robin | Unsupervised subpixelic classification using coarse-resolution time series and structural information[END_REF], and abundances estimation [START_REF] Debba | Abundances estimation of spectrally similar minerals by using derivative spectra in simulated annealing[END_REF].

The method proposed in this paper is in three steps. In a first step, a coarse classification is performed, based on the probabilistic output of an SVM. Every pixel can be assigned to a class, if the probability value obtained in the classification process is greater than a chosen threshold, or be unclassified. Pixels with a low probabilistic output are either mixed pixels or pixels hard to classify due to spectral variability, and their classification is addressed in a second step. In the second step, spectral unmixing is performed on the unclassified pixels by considering the preliminary results of the coarse classification step and by applying a Fully Constrained Least Squares (FCLS) method to every unlabeled pixel, in order to obtain the abundances fractions of each land cover type. Finally, in a third step, spatial regularization by SA is performed to obtain the resolution improvement. Experiments are carried out on synthetic and real hyperspectral data sets. The results are excellent both numerically and visually and show that the proposed method clearly outperforms traditional hard classification methods when the data contain mixed pixels.

The remainder of the paper is organized as follows. Section II presents in greater details the proposed approach. Section III shows the experiment on a synthetic data set, while Section IV illustrates the experimental results on real hyperspectral data. Section V finally draws the conclusions.

II. METHODOLOGY

The flow chart scheme of the proposed approach is presented in Fig. 2. The hyperspectral data are used as input for the hard classification method, in order to obtain a preliminary classification of all the pixels considered as "pure". The results of this step are the input (along with the original hyperspectral image) for the spectral unmixing, so that an appropriate set of endmembers can be found and the negative impact of spectral variabiality on the classification map minimized. In the last step, the results obtained are processed with a Simulated Annealing algorithm. Based on the assumption of spatial correlation of the land cover classes, SA is used to optimize a function where spatial proximity of pixels belonging to the same land cover class are preferred to the opposite case.

A. Pixel-wise classification

The first step of the proposed method consists in performing a pixelwise classification of the hyperspectral image, in order to obtain, for every pixel, a probability value for it to belong to one of the land cover classes. The pixels with a probability higher than a chosen treshold are considered as pixels where a single class is represented, and thus assigned to the considered class. These pixels are going to provide a preliminary classification map, where only the pixels containing a predominant 2. Flow chart scheme of the proposed approach land cover class are labeled. All the other pixels are not labeled, and their classification will be addressed in a second step. As long as it can provide a probabilistic output, every classifier can be used for the data analysis. In this work, we propose to use a probabilistic Support Vector Machine (SVM) classifier [START_REF] Vapnik | Statistical Learning Theory[END_REF], due to the good performances shown in the classification of hyperspectral data, also in case of limited training sets [START_REF] Melgani | Classification of hyperspectral remote sensing images with support vector machine[END_REF]- [START_REF] Fauvel | Spectral and spatial classification of hyperspectral data using svms and morphological profiles[END_REF]. In the following, we give a short description of the principles on which SVM is based. Due to the lack of space, we refer the reader to [START_REF] Vapnik | Statistical Learning Theory[END_REF] and [START_REF] Melgani | Classification of hyperspectral remote sensing images with support vector machine[END_REF] for further details on the theory of SVM and its application for hyperspectral images classification.

1) Probabilistic Support Vector Machine: The SVM is surely one of the most commonly used kernel learning algorithm. It performs robust non-linear classification of samples using the kernel trick [START_REF] Aronszajn | Theory of reprodusing kernel[END_REF]. The idea is to find a separating hyperplane in some feature space induced by the kernel function while all the computations are done in the original space [START_REF] Vapnik | Statistical Learning Theory[END_REF]. Given a training set S = {(x 1 , y 1 ), . . . , (x , y )} ∈ R n × {-1; 1}, the SVM computes a decision function f(x) such that sign(f(x)) can be used to predict the label of any test sample x. The decision function is found by solving the convex optimization problem:

max α g(α) = i=1 α i - 1 2 i,j=1 α i α j y i y j k(x i , x j ) subject to 0 ≤ α i ≤ C and i=1 α i y i = 0 (1)
where α are the Lagrange coefficients, k(x i , x j ) = k(x i , x j ) + δ ij /C, k the kernel function, C a constant that is used to penalize the training errors, δ ij a function such that δ = 1 if i=j, δ = 0 otherwise. To be an acceptable kernel, k should be a positive semi-definite function [START_REF] Aronszajn | Theory of reprodusing kernel[END_REF].

In [START_REF] Platt | Advances in Large Margin Classifiers[END_REF] Platt proposes approximating the posterior class probabilities P (y = 1|x) by a sigmoid function:

P (y = 1|x) ≈ P A,B (f ) = 1 1 + exp(A f + B) (2)
where f is an estimation of the decision function f(x) computed by the SVM, A and B two parameters that need to be optimized. The best parameter setting (A , B ) is determined by solving the following regularized maximum likelihood problem (with N + of the y i 's positive, and N -negative):

min z=A,B F (z) = l i=1 (t i log(p i ) + (1 -t i )(1 -log(p i )), (3)
where

p i = P A,B (f i ), t i = N++1 N++2 if y i = +1 1 N-+2 if y i = -1 , i = 1, . . . , l.
A detailed description of the method can be found in [START_REF] Platt | Advances in Large Margin Classifiers[END_REF]. In this work, we have used an improved implementation of the above algorithm [START_REF] Lin | A note on platt's probabilistic outputs for support vector machines[END_REF], which is included in the LIBSVM library [START_REF] Chang | Libsvm: A library for support vector machines[END_REF].

In this first classification step, we consider two outputs:

1) A complete probability map, containing the probability estimates for each pixel to belong to the assigned class. 2) A coarse classification map of the pixels considered as 'not mixed', containing class labels for the samples with a probability belonging to the class higher than a chosen treshold.

It is not a straightforward task to choose the treshold to determine if a pixel should be considered as pure or mixed. When labelling the pure pixels, we are interested in correctly classifying most of the pure pixels, because of two already mentioned reasons: i) When dealing with pure pixels, the full pixel methods work better than spectral unmixing, ii) this preliminary classification will be the input of the second step, and a large number of correctly classified pure pixel helps to provide suitable endmember candidates for the mixed pixels. For this same reason, the misclassification of pure pixels could lead to critical issues and cause a large error in the spectral unmixing step, thus a tradeoff is observed. The experiments carried out to investigate this issue have shown that in general a high treshold (close to 80%) allows to obtain a higher classification accuracy, since only the pixels which are reasonably sure to belong to a class are labeled. In Section III, it will be shown that the choice of this parameter is not crucial for the classification accuracy of the proposed method.

B. Spectral Unmixing

After obtaining a coarse classification map, where some pixels considered as "pure" (due to the high probability to belong to the assigned class) were classified, the labelling of the other pixels is addressed in the second step. Spectral mixture analysis (SMA) techniques have overcome some of the weaknesses of full pixel approaches by using linear statistical modeling and signal processing techniques [START_REF] Keshava | Spectral unmixing[END_REF], [START_REF] Plaza | Optical Remote Sensing -Advances in Signal Processing and Exploitation Techniques[END_REF]. They are inherently either nonlinear techniques or linear techniques. Nonlinear mixed pixel analysis estimate multiple scattering effects that may arise when the different materials form intimate association at microscopic level [START_REF] Keshava | Spectral unmixing[END_REF]. Although they can be useful for some types of analysis, in the majority of applications a linear mixing model can be considered without significant loss of information [START_REF] Keshava | Spectral unmixing[END_REF]. The key task in linear SMA is to find an appropriate set of pure spectral constituents -called "endmembers" in hyperspectral analysis terminology-, which are then used to estimate the fractional abundances of each mixed pixel from its spectrum and the endmember spectra by using a linear mixture model. In the Linear Mixture Model (LMM), the spectrum of a mixed pixel is represented as a linear combination of component spectra (endmembers). The weight of each endmember spectrum (abundance) is proportional to the fraction of the pixel area covered by the endmember. If there are M spectral bands, the spectrum of the pixel and the spectra of the endmembers can be represented by M -dimensional vectors. Therefore, the general equation for LMM is described as a linear regression form

z = L i=1 a i s i + e = As + e ( 4
)
where z is an M × 1 column pixel vector which describes the spectrum of the mixed pixel, s = [s 1 s 2 . . . s L ] is an M ×L endmember matrix of material signature, s i (i = 1, 2, . . . , L) are the M -dimensional spectra of the endmembers, a is an L × 1 column vector and is composed of abundance coefficients a i (i = 1, 2, . . . , L), e is an M -dimensional error vector accounting for lack-fit and noise effects, and L is the number of the endmembers. Due to physical reasons, (4) has to respect the following constraints of non-negativity (abundance fractions within a pixel cannot be negative) and sum to one (the sum of all the abundances fraction within a pixel must have 1 as a result):

a i ≥ 0 (5) L i=1 a i = 1. (6)
In recent years, several algorithms have been developed for automatic or semi-automatic extraction of spectral endmembers directly from the image data and to determine their fractional abundances within each pixel [START_REF] Plaza | A quantitative and comparative analysis of endmember extraction algorithms from hyperspectral data[END_REF]. Assuming that ground truth is available, we do not need to determine the endmembers composing the data, but simply the abundance of each land cover type within the pixels. In this case, a major issue is how to handle the spectral variability which affects the data. As shown in [START_REF] Foody | Variability in soft classification prediction and its implications for sub-pixel scale change detection and super-resolution mapping[END_REF], soft classification of hyperspectral images covering wide areas is negatively related to the intraclass spectral-variability, and the assumption that a single endmember could extensively represent a class is generally far from reality. The choice of appropriate endmembers is very important in order to correctly estimate the fractional abundances. If the endmembers do not represent the land cover classes well, the estimates of the sub-pixel coverage can be highly biased and lead to misclassification errors. In order to overcome this problem, we propose an adaptive approach to select the best endmember candidates for each pixel. This approach is based on two main assumptions:

1) The spatial correlation of the classes, i.e., for each pixel, it is probable that the best endmember candidates lie in the spatial proximity of the considered pixel.

2) The probabilistic output provided by the SVM, i.e., if a candidate is not spatially close to the selected pixel, but the probabilistic value of the class to which it belongs is high, it is presumably a good candidate. For each mixed pixel which has to be classified, we consider a set of 10 different spectra, that represent the endmember candidates. These candidates are chosen from the labeled samples of the training data and the results of the preliminary classification of step one, considered as a set of pure pixels correctly classified. If one of the land cover classes has a high probabilistic output (we consider a probabilistic output as high if its difference from the treshold chosen at step 1 is smaller than 5%), at least five spectra of this class are considered, otherwise all the 10 candidates are selected from the spectral signatures spatially closest to the considered pixel, after the coarse classification step. Once the spectral signatures representative of each class are extracted from the image, the abundance fraction of the elements within each pixel should be determined. Several algorithms have been developed for the linear mixing model according to the required constraints of abundances fractions. The fully constrained least squared unmixing algorithm is a widely adopted practical solution to avoid the appearance of theoretical problems, such as negative fractional abundances or abundances that sum up to more than one. Due to its computational efficiency, we have chosen this algorithm, which satisfies both abundance constraints and is optimal in terms of least squares error [START_REF] Heinz | Fully constrained least squares linear spectral mixture analysis method for material quantification in hyperspectral imagery[END_REF].

After applying FCLS, we obtain the fractional abundances of each endmember. Due to the fact that in many cases several endmembers represent the same class, by summing the fractional abundances of all the endmembers belonging to the same land cover class, we obtain the cover percentage of a class within a mixed pixel. It should be noted that not all the pixels to be classified in this second step are mixed, but there are many 'pure' pixels not labeled because of the low probability output provided by the SVM. However, the proposed method allows to label them as 'pure' pixels also in this second step.

C. Improving Spatial Resolution

Spectral unmixing is useful to describe the scene at a subpixel level, but can only provide information about proportions of the endmembers within each pixel. Since the spatial location remains unknown, spectral unmixing does not perform any resolution enhancement. In this paper, we propose a sub-pixel mapping technique, which takes advantage of the information given by the spectral mixing analysis and uses it to enhance the spatial resolution of thematic maps. Our proposed approach is as follows: In a first step, each pixel is divided in a fixed number of sub-pixels, according to the desired resolution enhancement. Every sub-pixel is assigned to an endmember, in conformity with its fractional abundance within the pixel. For example, if we want to have a zoom factor of N, we have to divide each pixel into N×N sub-pixels. For each pixel, the number of subpixels n to assign to the class i is computed according to the equation:

n i = round abd i 1/N , (7)
where abd i is the fractional abundance of the class i within the considered pixel estimated with the FCLS and round(x) returns the value of the closest integer to x.

A Simulated Annealing (SA) mapping function is then used, to create random permutation of these sub-pixels, in order to minimize a chosen cost function. Relying on the spatial correlation tendency of landcovers, we assume that each endmember within a pixel should be spatially close to the same endmembers in the surrounding pixels. Therefore, the cost function C to be minimized is chosen as the perimeter of the areas belonging to the same class:

C = I i=1 Ni j=1 P j , ( 8
)
where I is the number of the classes, N i is the number of connected components of the class i, and P j is the perimeter of the connected component j, computed according to the 8connected border pixels model [START_REF] Kulpa | More about areas and perimeters of quantized objects[END_REF]. SA is a well established stochastic technique originally developed to model the natural process of crystalization [START_REF] Metropolis | Equations of state calculations by fast computing machines[END_REF]. This process is based on an analogy from thermodynamics where a system is slowly cooled in order to reach its lowest energy state. More recently, SA has been proposed to solve global optimization problems [START_REF] Kirkpatrick | Optimization by simulated annealing[END_REF], and it has been used in various fields. The basic idea of the method is that, in order to avoid to be trapped in local minima, uphill movements, i.e., points corresponding to worse values of the objective function could, sometimes, be accepted for the following iteration. As with a greedy search, it accepts all the changes that improve the solution. Changes degrading the solution can be accepted, but with a probability that is inversely proportional to the size of the degradation (small degradations are accepted with a higher probability). This probability also decreases as the search continues, or as the system cools down, allowing eventual convergence to the optimal solution. An example of how SA spatial regularization works can be seen by looking at Fig. 3, where 3 (e) represents the initial sub-pixel distribution and 3 (f) the optimal one. First, a mixed pixed is selected, according to the information provided by the spectral unmixing step. Then, a random permutation of the subpixels within the chosen pixel is performed by SA. If this permutation leads to a decrease of the cost function (which is in our case the perimeter of connected components, i.e., components belonging to the same class), the change is accepted. Otherwise, as decribed above, the change will be probably rejected. The algorithm stops when minimum cost is reached, that is when a previously fixed number of consecutively rejected changes is reached (in our case, we set the number of consecutive operations to 100.000, since it is large enough to avoid sub-optimal solutions and Matlab takes only a few seconds to perfom this computation).

III. EXPERIMENTS ON SIMULATED DATA

The first experiment was carried out on a synthetic data set. The advantage of using a synthetic data set is the perfect knowledge that we have about the analysed image. This experiment has two main aims: the first one is to verify the assumption that the pixels classified by the SVM with a high probability value effectively correspond to 'pure' pixels. The second one is to validate the proposed method with a data set known in details, in order to evaluate the performances under different conditions of mixtures. The creation of a synthetic data set for classification purposes is not a trivial task, especially in this case where the spatial information is very important to fully exploit the potentiality of the proposed method. In order to build a data set as realistic as possible, we have considered a thematic map of an AVIRIS image, taken over the area of San Diego, with 9 land cover classes. We have created a hyperspectral data cube by substituting every class with a spectrum taken from the USGS spectral library, available on-line [START_REF] Clark | Usgs digital spectral library splib06a[END_REF]. The classes were chosen mainly from the vegetation library, in order to make more difficult the discrmination. The chosen classe were: shadow, asphalt, green grass, dry grass, maple, pine lodgepole, pine white, pinon and rub. The original image created was 400 × 400 pixels. Gaussian noise was added in order to reach a SNR of 30 dB. To have the possibility to analyse data sets where the ground truth cover is known in details, and to evaluate the obtained results from a quantitative point of view, we decided to use the original ground truth data only to compare the obtained results, and to decrease the spatial resolution of the image by applying an 4×4 low pass filter, so that we obtained an image of the same area with a resolution degraded of a factor 4. This way, we have the possibility to test the proposed method on a data set know in details, where there are pure, close to pure, and mixed pixels. The complete description of the data set is given in Table I. In order to compare the proposed approach with a common hard classification method, the same data were also classified with an SVM wih Gaussian Kernel, One vs One multiclass strategy and 10 fold cross-validation. Among the several multi-class strategies available for the SVM, we have chosen the One vs One because of its good performances in terms of robustness and computational burden [START_REF] Hsu | A comparison of methods for multi-class support vector machines[END_REF]. The comparison of the low resolution map obtained with SVM to the high resolution ground truth was not possible due to different number of pixels of the two images. However, we know that every pixel of the low resolution image corresponds to n × n pixels of the high resolution image. By comparing a pixel of the low resolution classification map with the n × n corresponding in the high resolution ground truth map, we can compute per-pixel classification accuracy. By doing this, we have to keep in mind that in case of a mixed pixel the hard classification method will inevitably lead to an error, because it will assign the corresponding high resolution n × n pixels to just one class, considered as predominant within the mixed pixel. However, this is exactly the issue that the proposed method is expected to address.

The performances of the two methods were compared in terms of overall accuracy (OA), that is the number of correctly classified test samples with respect to the total number of test samples, average accuracy (AA), which represents the average of the classification accuracies for the individual classes, and the Kappa coefficient of agreement (κ), that is a parameter that estimates the correct percentage classification without the amount that could be expected due to chance alone [START_REF] Cohen | A coefficient of agreement for nominal scales[END_REF]. In addition to this, we have computed the number of mixed pixels correctly classified, in order to show the performances of the two methods when dealing with mixtures of classes, and the number of pixels correctly labeled after the spectral unmixing step but uncorrectly positioned after Simulated Annealing.

For each class, 2% of the labeled samples were selected for training the algorithm. The treshold between pure and mixed pixels has been set to 0.7. The results of the experiment are shown in Table II. The first three columns represent the results obtained on the entire data set in terms of OA, κ and AA, and show that the proposed method provide an improvement in the overall accuracy classification. The second part of the Table II show a comparison of the performance of the two methods over groups of pixels with different degrees of purity, varying from pure pixels (where the predominant class has an abundance larger than 95%) to highly mixed pixels (where the predominant class has an abundance smaller than 55%). It can be noticed that while in case of pure and close to pure pixels the results of the two methods are quite similar, the proposed approach provided a dramatic increase of classification accuracy for mixed and highly mixed pixels, where a traditional classifier completely fails while the proposed method improves the accuracy of up to 35 percentage points. In order to evaluate the correctness of the assumption that pure pixels are classified with a high probability value, we have computed the mean probability of the maximum value of each group of pixels and the number of pixels of each group which are classified with a probability value higher than 70% (two leftmost columns of Table II). It can be easily noticed that the larger is the maximum value of abundance within the pixel, the higher is the probabilistic output provided by the SVM. The 80% of pure pixels were classified with a probability larger than 0.7, while in case of mixed and highly mixed pixels this quantity The experiments on real data were carried out considering three different data sets from two hyperspectral images. The first two data sets are from an Airborne Visible/Infrared Imaging Spectrometer (AVIRIS) image taken over NW Indiana's Indian Pine test site in June 1992. This image has been widely used in the remote sensing community for both classification and spectral unmixing purposes, and thus represents an interesting benchmark for the proposed method. According to [START_REF] Gualtieri | Support vector machine classifiers as applied to aviris data[END_REF], we considered for the first experiment a part of the scene, consisting of pixels × for a size of 86 × 68, which contains four labeled classes (the background pixels were not considered for classification purposes). We will refer to this data set as the "Subset scene". The second experiment was carried out on the whole AVIRIS data set. Sixteen land cover classes were considered. The original image is composed by 145 × 145 pixels. The calibrated data are available online 1 with detailed ground-truth information.

Finally, the third study site is the region surrounding the central-volcano Hekla in Iceland, one of the most active volcanoes in the country. Since 1970, Hekla has erupted quite regularly every 10 years, in 1970, 1980-81, 1991 and in 2000. The volcano is located on the South-Western margin of the Eastern volcanic zone in South Iceland. Hekla's products are mainly andesitic and basaltic lavas and tephra. AVIRIS data that were collected on a cloud-free day, June 17 1991, were used for the classification. The AVIRIS sensor operates in the visible, near-and mid-infrared portions of the electromagnetic spectrum, its sensitivity range spanning wavelengths from 0.4 µm to 2.4 µm. As on the previous case, the sensor system has 224 data channels, utilizing four spectrometers, whereas each spectral band is approximately 10nm in width. During the image acquisition, spectrometer 4 was not working properly. This particular spectrometer operates in the wavelength range from 1.84 µm to 2.4 µm (64 bands). These 64 bands were deleted from the imagery along with the first channels for all the other spectrometers, and the remaining 157 data channels were left. A subset of 180 × 180 pixels has been used for this experiment. In order to address the issue of the random choice of the training samples, for each data set we have repeated the experiment with ten different training sets. As in the previous experiment, due to the difficulty to have a perfect knowledge of the fractional abundances of each land cover type, we decided to use the original ground truth data only to compare the obtained results, and to decrease the spatial resolution of the image by applying an n × n low pass filter, where n varies according to the considered data set. This way, we know exactly the quantity of each class within a pixel, and we can use the low resolution image obtained after filtering as input for the proposed method. The information about the classes, the training and the test sets can be found in Table III.

A. AVIRIS subset

The first experiment was carried out on the AVIRIS subset image. The goal of this experiment is to illustrate the effectiveness of the method when used for the analysis of a simple hyperspectral data set. In this subimage, composed by 86 × 68 pixels, there are four classes with uneven number of labeled samples, namely, Corn-notill, Grass/Trees, Soybeans-no till, and Soybeans-min. The complete description of the training and test sets can be found in Table III. A 3 × 3 low-pass filter was applied to the original image, so that a new image with lower spatial resolution was obtained. The new image was composed by 28 × 23 pixels, and it was used as input for the proposed method. The low resolution image obtained after filtering and the ground truth can be seen in Fig. 5 (ab). Twenty pixels per class, considered as "pure" in the low resolution image, were randomly chosen and used for training the SVM classifier. The result of the classification with the SVM is presented in Fig. 5 (c). As it can be seen from the classification map, the two main problems are represented by the mixed pixels, which make hard to distinguish the border between different land cover areas, and the high spectral variability, which results in a noisy classification map. The proposed method provides an overall improvement for both issues (classification map shown in Fig. 5 (f)). It can be seen at the top of the image that the border of the Corn-no till field (represented in light blue in the map) is estimated with improved accuracy when compared to the traditional SVM classification, thus assessing the effectiveness of the proposed approach also to provide classification maps with a better spatial resolution. The problem of spectral variability is also solved, since the classification map is much less noisy and, therefore, a lower number of sparse pixels is observed.

The quantitative results obtained with the two methods are presented in Table IV, along with the results of the other data sets. The proposed method provides an improvement of the overall accuracy which is greater than 10%. As it will be in the other experiments, the improvement of the classification accuracy of the mixed pixel with respect the SVM is larger than for the whole data set. This demonstrates that the proposed method is effective in improving the results of data sets with mixtures of land cover classes. To assess the effectiveness of SA to locate sub-pixels in the classification map, we have also computed the number of sub-pixels correctly classified after the spectral unmixing step, but uncorrectly located after the spatial regularization. In this case, we can see that the error due to bad positioning of sub-pixels is extremely low.

B. AVIRIS complete

The second experiment was carried out on the whole AVIRIS data set. Sixteen land cover classes were considered for classification. The original image is composed by 145 × 145 pixels, and it was used as reference data. After applying a 2×2 low pass filter, an image composed by 72×72 pixels was obtained. The land cover ground truth can be seen in Fig. 6 (a). For training set, we have randomly selected, for each class, 15% of all the samples which were considered as "pure" in the low resolution image (that would correspond to about 10% of pixels of each class in the high resolution image). To have the possibility to compare the results of the proposed method with the available ground truth, we chose a zoom factor equal to 2, lower than in the previous case. However, the higher number of classes and their spectral similarity make this data set more challenging than the first one.

Figure 6 (b) and (e) shows the classification maps obtained with a conventional SVM and the proposed method. Also in this case, an improvement can be clearly seen in the classification maps, resulting in a less noisy map and an improved detection of the borders of spatial structures (in this case, agricultural fields). To have a quantitative comparison of the results obtained with the two methods, the overall accuracy of pixels correctly classified has been compared. The mean overall accuracy obtained in the five experiments with the SVM is 72.31%. As in the previous case, the low value of accuracy is due to two main factors, which are the impossibility of a common hard classification technique to distinguish different land cover classes at a sub-pixel level, and the difficulty to handle the high spectral variability. The proposed method obtained an average overall accuracy of 91.10%, showing the capability of the proposed approach to better deal with the aforementioned two main issues. By comparing Figure 6 (d) and (e), it can be noticed the effectiveness of the proposed spatial regularization with Simulated Annealing.

C. Hekla data set

For the last experiment, we consider a subset of the Hekla data set, located in the top-left corner of the scene. This subset is composed by 180 × 180 pixels, and it contains nine classes of interest. Also in this case a 2×2 low pass filter was applied to the original image, leading to a low resolution image of 90 × 90 pixels. Due to the insufficient availability of ground truth to quantify the results provided by the proposed method, we have considered as ground truth the classification map obtained by a spectral-spatial method, proposed in [START_REF] Tarabalka | Segmentation and classification of hyperspectral images using minimum spanning forest grown from automatically selected markers[END_REF], where the overall accuracy computed on the reference test set was close to 100%. Thus, also if we have to keep in mind that the results are estimated by comparison with a classification map and not with a selected land cover ground truth map, this classification map seems to be a reliable source of knowledge about the land class coverage of the area. The original image, the classification map obtained in [START_REF] Tarabalka | Segmentation and classification of hyperspectral images using minimum spanning forest grown from automatically selected markers[END_REF] and the obtained classification maps are shown in Fig. 7. 15% of the labeled pixels of the low resolution data were randomly selected from each class and used to train the classifier, and the experiment repeated ten times with different training sets.

The quantitative results in this experiment confirm those obtained in the previous ones: the proposed method provides not only a better classification map from a qualitative point of view, but also a large improvement of the overall accuracy of correctly classified pixels. Due to the irregular spatial structures in which the land cover classes are grouped, the spatial regularization method proposed in this paper was expected to be less effective than in the previous cases. The quantitative results agree with this supposition, also if the overall accuracy is penalized by only 2 percentage points. More advanced techniques could be investigated in our future works, in order to have an improvement of the classification accuracy of the mixed pixels.

D. Discussion about the choice of parameters and computational burden

In the proposed method, the parameters having an influence on the overall classification accuracy obtained (apart from the parameters of the SVM, which are automatically selected through cross-validation) are the treshold to distinguish between pure and mixed pixels and the number of 'endmember candidates' to consider for the spectral unmixing, in the second step. How the classification accuracy changes by changing the value of the parameters can be seen in Fig. IV-B. It can be noticed that the proposed method outperforms the traditional SVM in terms of accuracy over the whole range tested, being the choice of the parameters not crucial for the classification. As could be expected, a high value of the treshold to determine if a pixel can be considered as 'pure' provides a higher accuracy, since only the most reliable pixels are labeled for the preliminary classification. By setting a low value of the treshold parameter, the preliminary classification map will tend to be like the hard classification map obtained with a traditional SVM, thus decreasing the interest of the proposed method. The number of 'endmember candidates' considered in the spectral unmixing step points out the importance of the spatial information. The best results are in general obtained by considering a low number of candidates (which are the spatially closest to the considered pixel). When setting a larger value of the parameter, endmember candidates spatially far from the analysed pixel can be selected, introducing useless information and thus leading to a slight decrease in the classification accuracy.

The computational burden of the proposed method can be seen in Table V. The training of the SVM, which quadratically depends on the size of the training set, is the most computationally expensive step of the proposed approach. The spectral unmixing step depends on the number of the mixed pixel to unmix, while the Simulated Annealing regularization depends on the number of mixed pixels and on the zoom factor desired. When requiring a larger zoom, the number of possible sub-pixel combinations grows exponentially, thus requiring a heavier computational burden to reach the optimal configuration. In case of desired zoom factor equal or higher than 4, the computational burden of the spatial regularization is expected the be the most important.

V. CONCLUSIONS

Classification of hyperspectral images in presence of mixed pixels was addressed in this paper. A new method for the improvement of the spatial resolution of the classification maps was proposed. The method exploits the advantages of both soft classification techniques and spectral unmixing algorithms, in order to determine the fractional abundances of the classes at a sub-pixel scale. After the fractional abundances have been determined, spatial regularization by Simulated Annealing is finally performed to spatially locate the land cover classes within each pixel. Experiments were carried out on three different data sets and show that the proposed method clearly outperforms classical classification techniques when areas with mixtures of materials are located in the scene, providing excellent results both from a visually and quantitative point of view. Further research will be devoted to the investigation of advanced methods to better discriminate pure and mixed pixels, and of the possibility of alternative techniques of spatial regularization.

  Fig. 1.Illustration of the problem of mixed pixels: When the spatial resolution is not fine enough, several land cover classes lie in the same pixel. In this case, a hard classification process cannot give an accurate information about the pixel coverage, leading inevitably to a loss of information

Fig. 3 .

 3 Fig. 3. Basic steps of the proposed method: (a) A probabilistic classification map is computed for each class. (b)The pixels with highest probability greater than a chosen treshold are considered as pure and classified (in the figure, we set the treshold to 70%). The other pixels are considered as mixed (MIX in the figure). (c) For each mixed pixel, a set of possible endmembers is selected, considering the results of the preliminary classification. The other pixels, pure or mixed, are just ignored. (d) Spectral unmixing provides information about the abundance fraction of a class within each pixel. (e) Pixels are split into n sub-pixels, according to the desired zoom factor, assigned to an endmember and randomly positioned within the pixel. The number of sub-pixels assigned to each class reflects the fractional value estimated in the previous step. (f) Simulated annealing performs random permutations of the sub-pixels position until minimum cost is reached.

Fig. 4 .

 4 Fig. 4. Flow chart scheme of the proposed spatial regularization approach based on Simulated Annealing.

Fig. 5 .

 5 Fig. 5. AVIRIS subset data: (a) Low resolution image (band 30) obtained after applying a 3×3 filter to the high resolution image. (b) Ground truth of the high resolution image. Unknown pixels are represented in black. (c) Classification map obtained with traditional SVM. (d) Results of the preliminary classification. (e) Classification map obtained before applying the spatial regularization. (f) Final classification map obtained with the proposed method.

Fig. 6 .

 6 Fig. 6. (a) Ground truth data (b) Classification map obtained with an one versus one SVM (training set 1) (c) Classification of the ground truth pixels (d) Classification map obtained before applying the spatial regularization. (e) Final classification map obtained with the proposed approach. (f) Classification of the ground truth pixels.

Fig. 8 .

 8 Fig. 8. Variation of the overall classification accuracy versus the value of the parameter treshold to determine if a pixel can be considered as 'pure' for (a) AVIRIS Indian Pine subset (b) AVIRIS Indian Pine complete (c) AVIRIS Hekla data sets

TABLE II RESULTS

 II OF THE EXPERIMENT ON THE SYNTHETIC DATA SET. ABD REPRESENT THE maximum VALUE OF ABUNDANCE WITHIN A PIXEL. THE THIRD COLUMN DESCRIBES THE AVERAGE OF THE MAXIMUM PROBABILITY VALUE WITHIN A GROUP OF PIXELS WITH THE SAME PURITY DEGREE. THE LAST COLUMN REPRESENTS THE NUMBER OF PIXELS WITH A MAXIMUM PROBABILITY VALUE GREATER THAN 0.7, ACCORDING TO THEIR PURITY DEGREE.

		SVM	Proposed Method	Average Probabilistic SVM Output	Max(prob)>0.7
	OA	88.21%	90.84%	-	-
	κ	87.34%	89.52%	-	-
	AA	88.32%	90.25%	-	-
	abd>95%	99.50%	99.50%	0.74	5317
	95%>abd>85%	92.25%	91.44%	0.68	328
	85%>abd>75%	79.57%	77.57%	0.54	27
	75%>abd>65%	54.12%	61.17%	0.39	5
	65%>abd>55%	34.82%	50.33%	0.30	0
	abd<55%	15.36%	40.22%	0.26	1
	drops to 0.4%.
	IV. EXPERIMENTS ON REAL DATA

TABLE III INFORMATION

 III ABOUT THE TRAINING AND THE TESTING SET OF THE THREE CONSIDERED DATA SETS. IT HAS TO BE NOTICED THAT THE TRAINING SET IS SELECTED FROM THE LOW RESOLUTION IMAGE USED AS INPUT OF THE METHOD, WHILE THE TEST SET IS SELECTED FROM THE HIGH RESOLUTION REFERENCE DATA USED FOR COMPARISON.

		AVIRIS Indian Pine Subset	AVIRIS Indian Pine Complete	AVIRIS Hekla
	No.	Name	Train	Test	Name	Train	Test	Name	Train	Test
	1	Corn-no till	20	1434	Alfa Alfa	4	54	Andesite lava 1970	24	672
	2	Grass-Trees	20	747	Corn-no till	44	1434	And. lava 1980 I	126	3350
	3	Soybean-no till	20	727	Corn-min till	25	834	And. lava 1980 II	523	11916
	4	Soybean-min till	20	1926	Corn	7	234	And. lava 1991 I	220	4709
	5	-	-	-	Grass-Pasture	15	497	And. lava 1991 II	279	6918
	6	-	-	-	Grass-Trees	23	747	Lava tephra covered	103	2310
	7	-	-	-	Grass-Mowed	4	26	Rhyolite	6	181
	8	-	-	-	Hay-windrowed	17	489	Scoria	51	1286
	9	-	-	-	Oats	4	20	Firn-glacier ice	42	1058
	10	-	-	-	Soybean-no till	32	968	-	-	-
	11	-	-	-	Soybean-min till	83	2468	-	-	-
	12	-	-	-	Soybean-clean t	19	614	-	-	-
	13	-	-	-	Wheat	5	212	-	-	-
	14	-	-	-	Woods	42	1294	-	-	-
	15	-	-	-	Bldg-Trees-Drive	12	380	-	-	-
	16	-	-	-	Stone-Steel Tower	4	95	-	-	-

TABLE V COMPUTATIONAL

 V BURDEN FOR THE THREE DATA SETS.

		AVIRIS Indian Pine subset	AVIRIS Indian Pine complete	AVIRIS Hekla
	Approach	SVM	Proposed Method	SVM	Proposed Method	SVM	Proposed Method
	Classification	30 s	43 s	14 min	15 min	35 min	37 min
	Spectral Unmixing	-	15 s	-	1.5 min	-	2 min
	Simulated Annealing	-	50 s	-	3 min	-	4 min
	Total elapsed time	30 s	88 s	18 min	23 min	35 min	43 min

http://dynamo.ecn.purdue.edu/∼biehl/

Acknowledgments: This work has been supported by the European Community's Marie Curie Research Training Networks Programme under contract MRTN-CT-2006-035927, Hyperspectral Imaging Network (HYPER-I-NET).


Yacine Chitour
email: yacine.chitour@lss.supelec.fr

Sami Tliba
email: sami.tliba@lss.supelec.fr

H ∞ -feedback design for linear systems subject to input saturation

Introduction

In this paper, we study linear systems subject to input saturation, i.e. systems of the type (Σ) sat ẋ = Ax + bσ(u),

where x ∈ R n , u ∈ R and σ : R → R is a saturation function whose prototype is the standard saturation function σ 0 (s) = s max (1,|s|) . In this paper, we deliberately focus on single input systems, eventhough several results quoted below hold true for muti-inputs as well. The saturation-free linear system associated to (Σ) sat is naturally defined as

(Σ) L ẋ = Ax + bu, (2)
Our goal regards robustness issues associated to the (global asymptotic) stabilization to the origin of systems [START_REF] Angeli | Robust Stabilization via Saturated Feedback[END_REF]. More precisely, assume that there exists an non empty set K of static feedback laws (i.e. u = k(x)) so that (a) for each feedback law k ∈ K, the closed loop system is GAS with respect to the origin;

(b) the L 2 -gain associated to each k ∈ K is finite, i.e.

γ sat k := sup d∈L 2 (R + ,R) x d 2 d 2 < ∞,
where x d is the trajectory of ẋ = Ax + bσ(k(x) + d) with initial condition x(0) = 0.

Note immediately that if k is a linear feedback, then γ sat k ≥ Sγ L k , where γ L k is the L 2 -gain associated to k for the linear system (Σ) L and S is a positive constant only depending on σ.

Then, the issue we address is the following Estimate inf

k∈K γ sat k , (3)
eventually in terms of L 2 -gains of (Σ) L . In particular, if (A, b) is controllable, recall that inf k linear γ L k = 0. Our ultimate concern will be to achieve the same performance in the presence of saturation.

At the starting point of our study, (A, b) must be at least stabilizable, and by restricting the state space to the controllable space, we can (and will) assume, with no loss of generality, that the pair (A, b) is controllable. Moreover, it is well-known that (1) is stabilizable if and only if the eigenvalues of A must have non positive real part. Most delicate issues arise when the spectrum of A lies on the imaginary axis. The stabilization issue is already non trivial except for two dimensional systems which can be stabilized by linear feedbacks u = k T x. However, it was proved by Sussman and Yang ( [11] and also Fuller []) that the nth-integrator, n ≥ 3 cannot be stabilized by linear feedbacks u = k T x. Thanks to Teel [START_REF] Teel | Global stabilization and restricted tracking for multiple integrators with bounded controls[END_REF] and Sussmann, Yang and Sontag [START_REF] Sussmann | A general result on the stabilization of linear systems using bounded controls[END_REF], general and explicit static feedbacks were first constructed using nested saturations. However, as soon as the feedback requires at least two nested saturations, the corresponding L 2 -gain is infinite, eventhough one can obtain L 2 robustness results for these feedbacks at the price of a non linear concept of gain (cf. [START_REF] Teel | On L 2 performance induced by feedbacks with multiple saturations[END_REF]). One should also mention the construction of another type of feedbacks based on "minimal" ellipsoids and due to Megretsky (cf. [START_REF] Megretski | L 2 output feedback stabilization with saturated control[END_REF]). The evaluation of the L 2 -performances of these feedbacks is still an open problem. Finally, let us mention that, in the context of semi-global stabilization, the L 2 -performances of linear systems with saturations were resolved by Z. Lin and his coworkers, using a low-and-high gain design technique (cf. [START_REF] Hu | Control Systems with Actuator Saturation: Analysis and Design[END_REF] and references therein).

At the light of the above discussion, a reasonable chance to address the issue of estimating inf k∈K γ sat k only occurs when the family of stabilizing keedbacks K is large enough. To the best of the authors knowledge, this is the case for particular cases only, besides the use of nested saturations. For instance, when A is skew-symmetric (and more generally marginally stable), it was proved that the feedbacks u = -rb T x, r > 0 verify both Items (a) and (b) (cf. [7]) but inf r>0 γ sat r > 0, where γ sat r is the associated L 2 -gain. This is obtained by simply comparing γ sat r with the corresponding γ L r . However, for 2D and 3D systems and increasing saturation functions, in [START_REF] Angeli | Robust Stabilization via Saturated Feedback[END_REF], the family of stabilizing feedbacks K was enlarged to at least all the linear feedbacks for (Σ) L and it was also proved that the L ∞ -gains for these feedbacks were finite. It is not difficult to see that the same holds true for the L 2 -gains but no gain attenuation result seems straitgthforward from the results of [START_REF] Angeli | Robust Stabilization via Saturated Feedback[END_REF].

In the present paper, we return to the case of the 2D oscillator subject to input saturation and we determine an appropriate set K of linear feedbacks so that inf k∈K γ sat k = 0. We next deal with the case of the double integrator subject to input saturation i.e. ẍ = -σ(u), where x ∈ R. Eventhough the latter system is stabilizable by means of linear feedbacks, it was shown in [7] that the corresponding L 2 -gains are infinite for the main output map y = x. The same conclusion was reached for the output map y = ẋ in [START_REF] Kao | On the L 2 -gain of a double integrator with feedback loop saturation nonlinearity[END_REF] but it was proved in [START_REF] Chitour | esaim cocv[END_REF] that the L 2 -gain for the output map y = ẍ was finite, partially solving a Problem 36 as posted in [START_REF] Blondel | Open problems in Mathematical Systems and Control Theory[END_REF]. Still in [START_REF] Chitour | esaim cocv[END_REF], was provided a class of nonlinear feedbacks having a finite L 2 -gain for the output map y = x. In this paper, we show how to select a subset K of this class of feedbacks so that to get that inf k∈K γ sat k = 0. As future works, several directions seem interesting. The first one consists in extending first the present results to the general case of A skew-symmetric. We conjecture that linear feedbacks should be enough to get inf k∈K γ sat k = 0. It is definitely more difficult to reach the same results for marginally instable A and we suspect that new (non linear) feedbacks must be deviced for that purpose. Another line of research would be to investigate attenuation of incremental gains since partial results on that issue are available for A skew-symmetric ([8]).

Notations

We use L + 2 to denote the Hilbert space L 2 (R + , R) and, for d ∈ L + 2 , we use d 2 to denote the corresponding L 2 -norm of d.

Definition 1 We call σ : R → R a (normalized) saturation function (or an S-function) if for all t, t ′ ∈ R (i) |σ(t) -σ(t ′ )| ≤ inf(1, |t -t ′ |); (ii) |σ(t) -t| ≤ tσ(t).
Note that (i) is equivalent to the fact that σ is bounded and globally Lipschitz. On the other hand, (ii) is equivalent to

σ ′ (0) = 1, tσ(t) > 0 for t = 0, lim inf |t|→∞ |σ(t)| > 0, lim sup t→0 σ(t) -t t 2 < ∞.
For instance, arctan, tanh and σ 0 are examples of saturation functions. We use Σ(s) to denote, for s ∈ R the function defined by s 0 σ(t)dt. Note that Σ(s) > 0 for s = 0 and if σ is odd then Σ is an even function.

We also recall that for a linear system ẋ = Ax + bu, with (A, b) controllable, one has inf k linear γ L k = 0 where γ L k is the L 2 -gain associated to the output map u → x and a stabilizing linear feedback k. Indeed, to see that, it is enough to prove it for the Brunovsky form associated to ẋ = J n x + be n where J n is the n-th Jordan block and e n = (0 • • • 0 1) T . For the latter system, a direct computation using the family of dilations D µ := diag(µ n-i ) 1≤i≤n , µ > 0 yields the conclusion (see also Remark 4.3).

Oscillator in dimension two

We consider the linear control system subject to input saturation given by

(Σ) sat ẋ = A 0 x -e 2 σ(u), (4)
where

A 0 = 0 -1 1 0 , e 2 = 0 1
and σ is a real valued-function of "saturation" type.

The corresponding saturation free system is given by

(Σ) L ẋ = A 0 x -e 2 u. (5)
According to the Routh-Hurwitz criterion, the linear stabilizing feedbacks for (Σ) L are row vectors (k 1 k 2 ) with k 1 < 1 and k 2 > 0. It was proved in [START_REF] Angeli | Robust Stabilization via Saturated Feedback[END_REF] that these feedbacks also render (Σ) sat globally asymptotically stable with respect to the origin.

Proposition 3.1 (cf. [START_REF] Angeli | Robust Stabilization via Saturated Feedback[END_REF]) Let S be the subset of R 2 whose elements are the row vectors (k 1 k 2 ) with k 1 < 1 and k 2 > 0. Then, for every K ∈ S, the saturated system

ẋ = A 0 x -e 2 σ(k 1 x 1 + k 2 x 2 ),
is globally asymptotically stable with respect to the origin.

Proof. Consider the Lyapunov function given by

V (x 1 , x 2 ) = x 2 + 2k 1 k 2 1 + k 2 2 Σ(k 1 x 1 + k 2 x 2 ).
Then, the derivative of V along trajectories of ẋ = A 0 x -e 2 σ(k

1 x 1 + k 2 x 2 ) is equal to V = - 2k 2 k 2 1 + k 2 2 ξσ(ξ) 1 -k 1 σ(ξ) ξ ,
where

ξ := k 1 x 1 + k 2 x 2 . Since k 1 < 1 and 0 < σ(ξ) ξ ≤ 1 for ξ ∈ R, one gets that V ≤ -2k 2 (1-k 1 ) k 2 1 +k 2 2
ξσ(ξ) and then concludes by applying Lasalle's principle.

For every linear stabilizing feedback K ∈ S (Σ) L , the linear gain γ L K associated to the output map y = x is defined as

γ L K := sup d =0 x d 2 d 2 , (6)
where

d ∈ L 2 (R + , R) and x d is the unique solution of ẋ = (A 0 -e 2 K T )x -e 2 d with x(0) = 0.
It is well-known that γ L K is finite and verifies the classical Riccati criterium given next (cf. [START_REF] Green | Linear Robust Control[END_REF]): for every γ > γ L K , there exists a unique symmetric definite positive matrix P solution of the following Ricatti equation,

P (A 0 -e 2 K T ) + (A 0 -e 2 K T ) T P + 1 γ 2 P e 2 e T 2 P + I 2 = 0, (7)
and such that A + 1 γ 2 e 2 e T 2 P is Hurwitz. Moreover, the infimum of γ L K over all the linear stabilizing feedback K is equal to zero.

The goal of this note is to first determine a large set S ′ ⊂ S of linear stabilizing feedbacks K for (Σ) sat having a finite (saturated) gain γ sat K (associated to the output map y = x) i.e.,

γ sat K := sup d =0 x d 2 d 2 < ∞, (8)
where d ∈ L 2 (R + , R) and x d is the unique solution of ẋ = A 0 -e 2 σ(K T x + d) with x(0) = 0. Then, we will prove that the infimum of γ sat K over all the linear stabilizing feedback K of S ′ is equal to zero.

To proceed, we operate a linear change of variable for each feedback K ∈ S. We consider now the following saturated system

(Σ sat ) K ẋ = A 0 -bσ(ke T 2 x + d), (9)
where b := c θ s θ and k > 0. Here we use c θ and s θ to denote respectively cos θ and sin θ. The corresponding linear system (i.e., saturation free) is equal to

(Σ L ) K ẋ = (A 0 -kbe T 2 )x -bd, (10)
and the stability condition reads s θ > 0 and 1 + kc θ > 0.

Moreover, the Ricatti equation ( 7) becomes now

P (A 0 -kbe T 2 ) + (A 0 -kbe T 2 ) T P + 1 γ 2 P bb T P + I 2 = 0, (11)
In the sequel, we will use repeatedly the notation R α , α ∈ R, to denote the rotation of angle α, i.e.,

R α := c α s α -s α c α .
In particular, b = R α e 1 , where e 1 = 1 0 .

We next sketch the rest of the argument. We will first compute the linear gain γ L defined in (6) for θ < π/2 and close to π/2 and k large enough. We then compute the symmetric positive definite matrix P defined in (7) as well as the vector P b corresponding to the above mentionned of k and θ. In a last step, we provide an upper bound of γ sat in terms of γ L and conclude. Lemma 3.2 Assume that k > 0, θ ∈ (0, π/2) and 4 + 4kc θ + k 2 (c 2 θ -s 2 θ ) < 0. Then the linear gain γ L K of ẋ = (A 0 -kbe T 2 )x -bd, associated to the output map y = x, is equal to

γ L K = 1 1 + kc θ . ( 12
)
In particular, the conditions on k and θ in the above statement imply that θ ∈ (π/4, π/2). Moreover, the subscript K stands for a choice of (k, θ) in the range defined previously.

Proof. We proceed by a standard computation of the H ∞ -gain of the linear system ẋ = (A 0 -kbe T 2 )x -bd associated to the output map y = x. The closed-loop transfer mapping G(s) is given by

G(s) = 1 s 2 + ks θ s + 1 + kc θ -c θ s + s θ -s θ s -c θ .
Then one has the standard formula (cf. for instance [START_REF] Green | Linear Robust Control[END_REF]) asserts that

γ L K = G(s) ∞ = sup ω∈R σ(G(jω)),
where σ(G(jω)) is used to denote the maximum singular value of G(jω) for ω ∈ R. One easily gets that

σ(G(jω)) = 1 + ω 2 (ω 2 -1 -kc θ ) 2 + (ks θ ω) 2 .
One has therefore to determine the maximum of the function f R + → R + defined by

f (X) = 1 + X (X -1 -kc θ ) 2 + (ks θ X) 2 .
A straitghforward computation shows that, if 4 + 4kc θ + k 2 (c 2 θ -s 2 θ ) < 0, then the maximum is reached at X = 0 and is equal to

1 (1+kc θ ) 2 . Lemma 3.3 Assume that θ ∈ (0, π/2) and 4 + 4kc θ + k 2 (c 2 θ -s 2 θ ) < 0.
Then, the unique symmetric positive definite matrix P , solution of the Riccati equation, defined in (11) and such that A + 1 γ 2 e 2 e T 2 P is Hurwitz, is equal to

P = 1 ∆ R θ+ϕ p 1 p 2 p 2 p 3 + O( 1 k 2 ) R -(θ+ϕ) , (13)
In the above equation, O(•) stands the standard "Big O" notation as k or kc θ tends to infinity. All constants involved in this notation do not depend on k or θ ∈ (π/4, π/2).

Proof. . For γ > γ L K , the symmetric solution of the Ricatti equation P (γ) is positive definite. Define Q(γ) := P -1 (γ). Then Q(γ) is the symmetric positive definite solution of the equation

(Q + A)(Q + A) T = AA T - 1 γ 2 bb T . (17)
Assume that the above equation admits a symmetric positive definite solution Q for γ = γ L K . By a simple continuity argument, one gets that P (γ) tends to Q -1 as γ tends to γ L K , for fixed k and θ in the admissible range. Thus P (γ L K ) is well defined (as symmetric positive definite solution of (11) and equal to Q -1 ).

We can therefore solve (17) for γ = γ L K and, as we will notice ultimately, the symmetric solution to be found will be positive definite.

We first find that

AA T - 1 (γ L K ) 2 bb T = X 2 R θ+ϕ e 1 e T 1 R -(θ+ϕ) ,
where X = √ 1 + ks θ and ϕ is the angle defined by tan ϕ = -1 ks θ (here and after k is supposed to be large with respect to 1). Set S := R θ+ϕ e 1 e T 1 R -(θ+ϕ) . Then, there exists a rotation R η such that

Q + A = XSR η = XR θ+ϕ e 1 e T 1 R -(θ+ϕ)+η . Since Q is symmetric, one has A -A T = X(SR η -R -η S) = XR θ+ϕ (e 1 e T 1 R η -R -η e 1 e T 1 )R -(θ+ϕ) .
The matrix in parentheses in the last expression is skew-symmetric and therefore commutes with the rotation R θ+ϕ . That implies that the last equation reduces to

A -A T = X(e 1 e T 1 R η -R -η e 1 e T 1 )
.

A simple computation then yields

s η = 2 + kc θ X . (18)
On the other hand, one has

Q = - A + A T 2 + X SR η + R -η S 2 .
After some computations, one gets

Q = R θ+ϕ k(s θ R -η e 1 e T 1 R η + c θ 2 R -η (e 1 e T 2 + e 2 e T 1 )R η + X 2 (e 1 e T 1 R η + R -η e 1 e T 1 ) R -(θ+ϕ) ,
which in turn yieds that

Q = R θ+ϕ q 1 q 2 q 2 q 3 R -(θ+ϕ) ,
where

q 1 = ks θ c 2 ϕ + kc θ c ϕ s ϕ + Xc η , (19)
q 2 = -ks θ c ϕ s ϕ + kc θ 2 (c 2 ϕ -s 2 ϕ ) - Xs η 2 , (20)
q 3 = -kc θ c ϕ s ϕ -ks θ s 2 = c θ s θ + 1 ks θ + O( 1 k 2 ). ( 24
)
Setting ∆ := det Q, one gets at once (14) and then ( 13) and (15) for an analytical expression of P .

Recalling that b = R θ e 1 , one gets that

P b = c ϕ ∆ R θ+ϕ q 3 -q 2 -q 2 q 1 1 1 ks θ . Using further that c ϕ R ϕ = k 2 s 2 θ 1+k 2 s 2 θ (I 2 -1 ks θ A 0 ), one derives (16).
We can now address the question of H ∞ -feedback design for (Σ) sat . In order to access to upper bounds of γ sat K , we resort to a dissipative inequality involving an appropriate storage function. More precisely, assume that θ ∈ (0, π/2), 4 + 4kc θ + k 2 (c 2 θ -s 2 θ ) < 0 and k is large. We prove the following theorem. Theorem 3.4 There exists θ 0 ∈ (π/4, π/2) and k 0 > 0 such that, for every k ≥ k 0 , one has

γ sat K ≤ 30γ L K . (25)
In particular, for that choice of K (i.e., choice of (k, θ 0 ) with k ≥ k 0 ), (Σ) sat has a finite L 2 -induced gain for the output y = x and the infimum of these gains over all the stabilizing feedbacks is equal to zero.

Proof. . Along trajectories of ((Σ sat ) K ) given in [START_REF] Megretski | L 2 output feedback stabilization with saturated control[END_REF], consider the Lyapunov function W defined by

W (x) = x T P x + C 0 3 x 3 , ( 26
)
where P is symmetric positive definite matrix provided by Lemma 3.3 and C 0 > 0 is a constant to be determined. Set ξ := ke T 2 x + d. Then one has

(x T P x) • = -x 2 - 1 γ 2 L (x T P b) 2 -2d(x T P b) + 2(x T P b)(ξ -σ(ξ)). (27)
It is the last term of the above equality which cannot be controlled, neither by the term "-x 2 " nor by a completion of squares involving -1

γ 2 L (x T P b) 2
. This is why one must add another term to x T P x in order to achieve to a dissipation inequality.

Consider the following

( x 3 /3) • = x x T ẋ = -x (c θ e T 1 xσ(ξ) + ξ -d k s θ σ(ξ)). (28)
We next take C 0 = 2k P b s θ . The above equation yields

( C 0 3 x 3 ) • ≤ C 0 c θ x 2 + 2 P b s θ x |d| -2 x P b ξσ(ξ). (29)
Using Item (ii) in Definition 1, we can get rid of the last term of (27) with the last term of (29). One then gets

Ẇ ≤ -(1 -C 0 c θ ) x 2 + 2 P b s θ x |d| - 1 γ 2 L (x T P b) 2 + 2|d||x T P b|. (30)
By choosing k large enough, one has

C 0 c θ ≤ 3c θ s 2 θ and 2 P b s θ ≤ 6 ks 2 θ .
By choosing θ = θ 0 close enough to π/2 and independently of k, one also has

C 0 c θ 0 ≤ 1/2 and 2 P b s θ 0 ≤ 7γ L K .
For the choice θ = θ 0 and for k large enough, (30) becomes

Ẇ ≤ - 1 2 x 2 + 7γ L K x |d| - 1 (γ L K ) 2 (x T P b) 2 + 2|d||x T P b|.
After completion of squares, the above inequality can be written as

Ẇ ≤ - 1 4 x 2 + (1 + 14 2 )(γ L K ) 2 |d| 2 . ( 31
)
It is then immediate to derive (25) from the above inequality.

Remark 3.5 The choice of the Lyapunov function W as defined in (26) was already performed in [7]. However, in the abovementionned reference, the symmetric positive definite matrix P was taken as solution of the Lyapunov equation A T P + P A = -I 2 . Moreover, the feedback was taken with the choice k > 0 and θ = 0. The corresponding linear gain was always equal to one. , which then allows to use the angle θ as an extra degree of freedom in order to get C 0 c θ < 1. It was therefore important to have a "good" upper bound of P b . One can derive a first estimate of P b without an exact computation of P itself. Indeed, by taking the trace in (11), one deduces that

P b ≤ k(γ L K ) 2 + γ L K (kγ L K ) 2 -2.
Then the smallest value for C 0 c θ is obtained as θ tends to π/4 and is equal to 4, which is not "good" enough for our purposes.

Double integrator

In this section, we consider the double integrator system subject to input satutation, i.e. the linear control system defined by (DI) sat ẋ = J 2 x -e 2 σ(u),

where x = (x 1 , x 2 ) T , J 2 is the 2D Jordan block, i.e. J 2 = 0 1 0 0 and σ is a real valuedfunction of "saturation" type.

In [START_REF] Chitour | esaim cocv[END_REF], non linear feedbacks were introduced to obtain finite L 2 -gains. For that purpose, we define the class of F-functions as follows:

Definition 2 A function F : R → R is an F-function if F is C 1 , odd, F ′ (0) = 0 and there exist r ≥ 1 such that, for |x 2 | ≥ 1, we have

F ′ (x 2 ) ≥ sup 3|x 2 |, r F (x 2 ) x 2 . ( 33
)

Remark 3 . 6

 36 One of the critical inequalities obtained previously was C 0 c θ := 2kc θ P b s θ ≤ 3c θ s 2 θ

We now recall the result of [START_REF] Chitour | esaim cocv[END_REF] to be used later. Proposition 4.1 Let us consider the control system (DI) sat with σ an increasing saturation function, and the feedback k F (x) given by

where F an F-function. Then k F is a feedback stabilizer for (DI) sat with finite L 2 -gain for the output map y = x.

For instance

2 are examples of k F -feedbacks. Then, the main result of this section is the following. Theorem 4.2 Consider the double integrator system subject to input satutation as defined by in (32) and a feedback k F defined by some F-function. For µ > 0, consider the feedback

If γ(µ) is used to denote the L 2 -gain associated to k µ F for the output map y = x, then, for µ ≥ 1, one has

where γ( 1) is L 2 -gain associated to k F for the output map y = x. As a consequence,

Proof. For µ > 0, consider the feedback k µ F given in (35). Make the following change of variable and time

The state

According to Proposition 4.1, one has

Since one has

one concludes readily.

Remark 4.3 Actually, the above technique of dilations is standard and can be generalized as well for the n-th integrator subject to input saturation, in the case where there exists a stabilizing (static) feedback law with finite L 2 -gain. It works as follows.

Assume that k is a stabilizing (static) feedback law with finite L 2 -gain i.e.,

where x d is the solution of the perturbed system ẋ = J n x + e n σ(k(x) + d) with x(0) = 0. For µ > 0, set D µ = diag(µ n-i ) 1≤i≤n . Then, one has

Also set k µ (•) := k(µD µ •) for µ > 0. Then

where z d is the solution of the perturbed system ẋ = J n x + e n σ(k µ (x) + d) with x(0) = 0. Indeed, x µ d (•) := µD µ z d (•/µ) is the solution of ẋ = J n x + e n σ(k(x) + d(•/µ)), with x µ (0) = 0. For µ ≥ 1, one immediately gets,

and then deduces that γ sat µ ≤ γ sat /µ and concludes readily.


Bin Luo

Xavier Ceamanos

Sylvain Douté

Jocelyn Chanussot

MARTIAN AEROSOL ABUNDANCE ESTIMATION BASED ON UNMIXING OF HYPERSPECTRAL IMAGERY

Classical linear unmixing approaches are not valid if the atmosphere and the aerosol are present in hyperspectral data, since the mixture model is no longer linear. In this paper, we present an iterative approach for estimating the abundance of aerosol based on unmixing of Martian hyperspectral data. On one hand, the results can provide the information on the aerosol of the Mars, which is very difficult to obtain. On the other hand, we can use the result to remove the effect of the aerosol on the original image and obtain more accurate linear unmixing results on the surface reflectance.

INTRODUCTION

Visible and near infrared imaging spectroscopy is a key remote sensing technique to study the Earth and other planets. As solar light is partially transmitted, reflected and diffused back by interaction with the different constituents of the atmosphere and the surface, the analysis of reflectance spectra may allow the identification, quantification, and characterization of the chemical species. In [START_REF] Luo | Unsupervised endmember extraction: Application to hyperspectral images from mars[END_REF], the authors proposed a linear unmixing scheme for unmixing the spectra of the pixels of Martian hyperspectral images in order to identify and quantify the chemical species on the surface of the Mars. However, the mixture of the spectra of the pixels are supposed to be linear. This linear mixture model is valid for the simple surface reflectance. For the intimate mixture of the surface reflectance, the model is no longer linear. More complicated approaches are needed. In [START_REF] Broadwater | A comparison of kernel functions for intimate mixture models[END_REF], it is proposed to use the kernel models in order to quantify the chemical species which are intimately mixed. In addition, the surface reflectance is transmitted by the atmosphere and the aerosol. The transmission is also non linear. In this paper, the spectra received by the sensor are assumed to be the simple surface reflectance (which is usually a linear mixture) transmitted by the atmosphere and the aerosol. The transmission of the gases on the Martian atmosphere is known, since the effect of the atmosphere mainly depends on the topography of the planet. However, the transmission of the aerosol depends on the thickness of the aerosol, which varies according to the captured scene. The objective of this paper is to estimate the abundance of the aerosol as the This work is funded by french ANR project VAHINEES. optical depth by using linear unmixing approach. The scheme of this paper is presented in Figure 1. Given a hyperspectral image, where the spectra are the surface reflectance transmitted by the atmosphere and the aerosol. We at first correct the transmission effect of the atmosphere and the aerosol by assuming that the thickness of the aerosol is constant in the image. Afterwards, we use the linear unmixing approach to separate the surface reflectance. We then use the results of the linear unmixing for estimating the thickness of the aerosol on each pixel of the image. The above steps make up the first iteration of the algorithm. We then repeat the iteration to improve the accuracy of the estimation of the thickness of the aerosol.

The outline of the paper is as follows: in Section 2, we introduce the mixture model of the received spectra which takes into account the transmission of the atmosphere and the aerosol. In Section 3, the algorithm for estimating the abundance of the aerosol on each pixel is presented. In Section 4, the results of the estimation on synthetic data are presented. Conclusions are prensented in Section 5. 978-1-4244-8907-7/10/$26.00 ©2010 IEEE

NONLINEAR MIXTURE MODEL OF THE ATMOSPHERE AND THE AEROSOL TRANSMISSION ON THE MARS

If we take into consideration the atmosphere and the aerosol transmission, the received reflectance on a pixel x without noise is:

I(x; λ) = m AT M (x; λ) (x) R(x; λ), (1)
where λ is the wavelength. R(x; λ) is the surface reflectance, which is supposed to be a linear mixture of N c endmembers, i.e.

R(x

; λ) = Nc i=1 m(x; i)s(i; λ) (2)
where m(x; i) is the abundance of the ith endmember at pixel x and s(i; .) is its spectrum. m AT M (x; .) is the transmission of the atmosphere at pixel x which mainly depends on the topography of this pixel. (x) is the parameter related to the effect of the aerosol and gas at pixel x to be estimated, which will be detailed in Section 2.1.

Parameter (x)

can be further decomposed into two terms :

(θ 0 , θ, φ, τ aer , H scale ) = ψ(airmass). β(θ 0 , θ, φ, τ aer , H scale ) (3)
where θ 0 , θ and φ are the solar incidence angle, the spacecraft emergence angle and the azimuth angle, respectively. In our case, the geometry can only be explained by θ 0 since OMEGA (carried by Mars Express, 256 spectral bands, 300m resolution) mostly works at nadir. As a consequence, θ and φ are negligible. τ k0 aero is the column integrated opacity of the aerosols and H scale is the scale height of the distribution. In the experiment, H scale was set equal to 11.5 km, which is a typical value of Martian polar regions.

The airmass represents the distance that the photons travel until reaching the spacecraft and is given by

airmass = 1 cos(θ o ) + 1 cos(θ) (4)
Since we know perfectly φ(airmass), the estimation results of allow us to deduce the effect of the aerosol β(θ 0 , θ, φ, τ aer , H scale ), which corresponds directly to the abundance of the aerosol.

ALGORITHM FOR THE ESTIMATION OF (X)

We use an iterative scheme for estimating the parameter (x).

In a first step, the effect of the atmosphere and the aerosol is corrected by assuming that (x) is a constant everywhere in order to obtain the surface reflectance. In a second step, the linear unmixing approach is used to separate the surface reflectance. In a third step, the results of the linear unmixing are used for estimating the abundance of the aerosol. The result of the estimation is again used to correct the effect of the atmosphere and the aerosol for the next iteration. After a given number of iterations, the algorithm stops. The scheme of the algorithm is shown in Figure 1.

More concretely, we suppose that j (x) is the (x) estimated after the jth iteration. We use this value to correct the received spectra in order to obtain the surface reflectance R j (x; λ) according to Equation (1):

R j (x; λ) = I(x; λ) m AT M (x; λ) j .
(

) 5
Since we suppose that the surface reflectance is a linear mixture of the spectra of several endmembers (see Equation (2), we can separate the surface reflectance R j by using the linear unmixing method -the Minimum Volume Simplex Analysis (MVSA) [START_REF] Li | Minimum volume simplex analysis: A fast algorithm to unmix hyperspectral data[END_REF]. The sum of the abundances of all the endmembers for each pixel is equal to one, i.e. ∀k, Nc n=1 s n,k = 1, the data vectors R j are always inside a simplex of which the vertex are the spectra of the endmembers if there is no noise. The Vertex Component Analysis (VCA) which extracts the extrema of the simplex as the endmembers is proposed in [START_REF] Nascimento | Vertex component analysis: A fast algorithm to unmix hyperspectral data[END_REF] as an effecient algorithm for extracting the endmembers. However, the VCA requires at least one pure pixel (which contains only one endmember) for each endmember. MVSA can find the optimal simplex contains all the data vectors of R j and extract its extrema as endmembers, even though there is no pure pixel.

By using MVSA, we can obtain the spectra of endmemebers mj (.; i), as well as their spatial abundances ŝj (i; .). We then reconstruct the surface reflectance by using Equation ( 2)

Rj (x; λ) = Nc i=1 mj (x; i)ŝ j (i; λ). (6)
Finally we use this reconstruction to re-estimate the parameter j+1 (x). We suppose the optimal j+1 (x) minimise the reconstruction error between I and m AT M (x; λ) j+1 (x) R j (x; λ), i.e.

j+1 (x) = arg min ∈R I(x; λ)-m AT M (x; λ) j+1 (x) R j (x; λ) 2 .
(7) Equation ( 7) leads to

j+1 (x) = λ m AT M (x; λ)(I(x; λ) -Rj (x; λ)) λ m 2 AT M (8)
Afterwards, we recompute the updated surface reflectance R j+1 (x; λ) for the next iteration by using j+1 (x) estimated by Equation (8). The complete scheme for estimating (x) is then:

• For j = 0 to N -Compute Rj by using Equation ( 5);

-Use linear unmixing approach to unmix Rj for obtaining mj and ŝj;

-Reconstruct Rj by using Equation (6); -Update j+1(x) by using Equation (8).

• end

EXPERIMENTS AND RESULTS

Simulated data

Experiments are performed on a 128×128-pixel synthetic data for evaluation of the aerosol retrieval strategy. Two Martian-like endmembers are present on the surface: water ice and dust. Two complementary abundance maps were generated to simulate a sub-pixel geographic mixture (see Figure 2). Uncertainty was introduced in the data by randomly selecting the grain size of both surface components. Two physical-meaningful ranges of grain size values were chosen according to Martian properties. Then, each surface spectrum was altered by a layer of aerosol which τ aero is given by an aerosol abundance map (see Figure 3c). This map was generated by plasma fractals to give a cloud-like realistic result. τ aero values were selected from 0.1 to 1.6 to correspond with a Mars-like scenario. Synthetic spectra were differently generated by the method in [START_REF] Douté | A multilayer bidirectional reflectance model for the analysis of planetary surface hyperspectral images at visible and near-infrared wavelengths[END_REF] depending on the value of θ 0 -from 65 to 85 expressed by the geometry map in Figure 3(a). The high values of airmass mimics the polar regions of the Mars. Finally, the surface cube was multiplied by a data set containing the atmospheric transmission of a polar OMEGA image. That was done to reproduce the atmospheric gas effects as well as the topography of the scene. As a matter of fact, the atmospheric contribution on the surface changes depending on the elevation.

Estimation results

We have performed 3000 iterations (of which the scheme is shown in Section 3) to estimate the (x). The initial value of 0 = 2, which is between the maximal and the minimal value of the ground truth. In order to evaluate the estimation results, we have computed between the ground truth and each j (x) obtained at the jth iteration, the correlation coefficient C j , the mean squared error (MSE) E j respectively by

C j = E( j -E( j ))( -E( )) E( j -E( j )) 2 E( -E( )) 2 . E j = E(( j -) 2 ).
The larger C j is, the better j (x) is, while the smaller E j is, the better j (x) is. In Figure 4, we have shown C j and E j as the function of iteration step j. It can be seen that, the correlation coefficient between the estimated (x) and the ground truth always increases (from 0.71 to 0.75), which indicates that the estimated (x) is more and more similar to the ground truth. However, after 1000 iterations, the correlation coefficient increases very slightly. In addition the MSE between the estimated (x) and the ground truth always decreases.

It can be seen that the correlation coefficients always increase. In Figure 5 Afterwards, we used MVSA to unmix the surface reflectance R 3000 (computed by using Equation ( 5)). The number of endmember is p = 3. In Figure 6, we have shown the spectra of these endmembers. It can be seen that the second and the third endmembers correspond to the spectra of the water under different aerosol. While the first endmember is the specrum of dust. It can be seen that the effect of aerosol decreases a lot for the spectra of water. However, the effect of aerosol is still visible on the extracted spectrum of dust.

Figure 2 shows the ground truth of the two endmembers used for simulating the image.

It can be seen that the first abundance map extracted by the linear unmixing (see Figure 7

CONCLUSION

In this paper, we present an iterative approach for estimating the abundance of aerosol based on unmixing of Martian hyperspectral data. This method is completely statistic without exploiting any information a piori on the properties of the surface components. The results obtained from synthetic data show that this approach can provide satisfactory estimation of the abundance of the aerosol. On one hand, the results can provide the information on the aerosol of the Mars, which is very difficult to obtain. On the other hand, we can use the result to remove the effect of the aerosol on the original image and obtain more accurate linear unmixing results on the surface reflectance.

Fig. 1 .

 1 Fig. 1. Scheme of the algorithm for estimating of the abundance of the aerosol.

Fig. 2 .

 2 Fig. 2. Ground truth of the abundance of the two endmembers for the simulated image.

Fig. 3 .Fig. 4 .

 34 Fig. 3. (a) Incident angle; (b) The ground truth of (x); (c) the ground truth of aerosol abundance.

  (b), we have shown the 3000 (x) estimated at the 3000th iteration. For comparison, we have shown in Figure5(b), the ground truth of the (x) is shown.

Fig. 5 .Fig. 6 .

 56 Fig. 5. (a) The estimated 3000 (x) after 3000 iterations

  (a)) is approximately the abundance map of the Dust. The other abundance maps extracted by unmixing (see Figure 7(b) and (c)) correspond to the abundance map of the Water.

Fig. 7 .

 7 Fig. 7. (a)-(c) Abundances of the three endmembers extracted by MVSA after removal of atmoshpere and aerosol effect by using 3000 (x).