Mathieu Fauvel

Alberto Villa

Jocelyn Chanussot

MISTIS Jon Atli Benediktsson

Mahalanobis kernel for the classification of hyperspectral images

Keywords: Mahalanobis kernel, probabilistic principal component analysis, support vector machine, hyperspectral images, classification

INTRODUCTION

Kernel methods have received a considerable attention during the last decade [START_REF] Hofmann | Kernel methods in machine learning[END_REF]. Their performances for classification, regression or feature extraction make them popular in the remote sensing community [START_REF]Kernel Methods for Remote Sensing Data Analysis[END_REF]. The core of kernel learning algorithms is the kernel function. It measures the similarity between two spectra x and y in a d-dimensional vector space, and enable the transformation of a linear algorithm into a non-linear one [START_REF] Vapnik | The Nature of Statistical Learning Theory, Second Edition[END_REF]. Over the different kernels used in remote sensing, the Gaussian kernel:

k(x, y) = exp - 1 2 d i=1 (xi -yi) 2 γ 2 . ( 1
)
is widely used because it usually gives good results and has only one hyperparameter (γ) to be tuned. By using kernels, an implicit mapping of the spectra from the input space to the feature space is done, the feature space being associated to the kernel [START_REF] Hofmann | Kernel methods in machine learning[END_REF]. Under some weak conditions, the feature space induced by the kernel is a Riemannian manifold [START_REF] Scholkopf | Geometry and Invariance in Kernel Based Methods In Advances in Kernel Methods -Support Vector Learning[END_REF][START_REF] Williams | A geometrical method to improve performance of the vector machine[END_REF]. The metric tensor is

gij(x) = ∂ 2 k(x, y) ∂xi∂yj y=x (2)
which is for the Gaussian kernel: gij(x) = γ -2 δij with δij = 1 if i = j and 0 otherwise. This metric stretches or compresses the Euclidean distance between x and y by a factor γ -2 and the implicit statistical model associated to the input data x is the normal law with a diagonal covariance matrix and identical elements [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF]Chap. 4]: Each variable1 has the same variance and there is no covariance between variables (which is not true in practice). The model also assumes that each variable is equally relevant for the given task, e.g., classification or regression. A more advanced model is to consider that the data follow a normal law with a diagonal covariance matrix, but with no identical diagonal elements: Each variable has its own variance, but still no covariance. It is then possible to tune the relevance of each variable separately. It was shown in [START_REF] Villa | Gradient optimization for multiple kernel's parameters in support vector machines classification[END_REF] that this model improves the classification accuracy of remote sensing images, but it also increases the computational load. The more general model, full covariance matrix, leads to the well known Mahalanobis kernel (mk) [START_REF] Camps-Valls | Hyperspectral image classification with Mahalanobis relevance vector machines[END_REF]:

k(x, y) = exp - 1 2γ 2 x -y t Q x -y . ( 3
)
Several definitions of Q exist for the problem of classification: It can be either the inverse of the covariance matrix Σ of the total training samples [START_REF] Camps-Valls | Hyperspectral image classification with Mahalanobis relevance vector machines[END_REF] or the covariance matrix of the considered class [START_REF] Abe | Training of support vector machines with Mahalanobis kernels[END_REF], e.g., for m classes problem, if the classifier separates the class c against all the others, Q is Σ -1 c . Generally, Q can be any positive definite matrix. The metric induced is gij(x) = γ -2 Qij, which stretches or compresses the variance of the data along their principal directions (see Section 3). Although the mk better handles the data characteristics than the conventional Gaussian kernel in small/moderate dimensional space, it is difficult to use in high dimensional spaces, such as hyperspectral remote sensing images. As a matter of fact, the estimation of the covariance matrix is ill-conditioned, making its inversion impossible or at least unstable. Moreover, all the principal directions are not equally relevant for classification: A subset of them corresponds to the signal while the remaining dimensions correspond to the noise.

In this article, it is proposed to regularize Q in a suitable way and to tune the weight of the principal directions according to their relevance in classification. The regularization strategy is detailed in Section 2. The link between probabilistic PCA and the proposed regularization is also discussed. The Mahalanobis kernel is then constructed with the regularized inverse covariance matrix in Section 3. Finally, experiments on real hyperspectral data sets are reported in Section 4.

REGULARIZATION OF COVARIANCE MATRICES

From the spectral theorem, a covariance matrix can be written as:

Σ = VΔV t ( 4
)
where Δ is the diagonal matrix of eigenvalues and V is an orthonormal matrix of corresponding eigenvectors. Its inverse is

Σ -1 = VΔ -1 V t . ( 5
)
Ill-conditioning is related to a high condition number κ(Σ), i.e., the ratio between the largest δ1 and the smallest δ d eigenvalue [START_REF] Vogel | Computational Methods for Inverse Problems[END_REF]. In general, in hyperspectral imagery, Σ is not full rank because of the high correlation between two adjacent bands. One consequence of this is a high condition number. Another consequence is that not all principal directions carry the relevant signal and thus it is possible to discard principal directions corresponding to zero (or close to) eigenvalues. Following [START_REF] Bernard-Michel | Machine learning techniques for the inversion of planetary hyperspectrales images[END_REF][START_REF] Bernard-Michel | Gaussian regularized sliced inverse regression[END_REF], it is proposed in this article to use the PCA-Ridge regularization based approach: Noting I p d the diagonal matrix with the p first elements equal to 1 and the remaining equal to 0 and defining Ω = VI p d V t , the ill-posed problem [START_REF] Williams | A geometrical method to improve performance of the vector machine[END_REF]

is changed to ΩΣ + τ I d -1 Ω which is, after trivial simplifications, equal to VΛ(τ, p)V t with Λ(τ, p) = diag 1 δ1 + τ , . . . , 1 δp + τ , 0, . . . , 0 . (6)
Usually, the regularization parameter τ is set to a small value: δi +τ ≈ δi if δi is high enough. Therefore, principal directions corresponding to high eigenvalues are slightly regularized and those corresponding to small eigenvalues are largely regularized. The remaining (d-p) principal directions are discarded. Finally, with the regularized estimate of the covariance matrix, xy 2 Σ -1 can be rewritten as:

x -y 2 Σ -1 = x -y t VΛ(τ, p)V t x -y = VΛ 0.5 (τ, p) t x -y 2 = A t x -y 2 (7)
with This model is closely related to probabilistic principal component analysis (PPCA) [START_REF] Tipping | Probabilistic principal component analysis[END_REF]. In PPCA, it is assumed that the observed d variables x are a linear combination of p unobserved variables s (p d):

A = v1 (δ1 + τ ) 0.5 . . . vp (δp + τ ) 0.5 (8)
x = Ws + μ + ε (9)
with s ∼ N(0, I) and ε ∼ N(0, 2 I). It follows that x ∼ N (μ, WW t + 2 I). With PPCA, it is clear that all the information/signal in x is contained in s. It can be proved that s lives in a subspace spanned by the p first eigenvectors of the covariance matrix of x [START_REF] Tipping | Probabilistic principal component analysis[END_REF]. This subspace correspond to A, that was previously defined.

Using [START_REF] Tipping | Probabilistic principal component analysis[END_REF], A is estimated using the p first eigenvalues and corresponding eigenvectors computed on the empirical sample estimate of Σ:

Σ = 1 n n i=1 (xi -x)(xi -x) t ( 10
)
with n the number of training samples and x being the sample mean of x.

It results that the inverse of the covariance matrix can be computed in a closed form from the empirical covariance matrix. The parameter p controls the model: the size of A where the data originally lived. Therefore, it is important to tune p correctly, since if p is set to a small value relevant information is lost, while with a too big value some noise is included in the computation. In remote sensing applications, it is conventionally selected by considering a certain percentage of the cumulative variance (total eigenvalues sum). A more robust approach is to use statistical model selection criteria [START_REF] Schwarz | Estimating the dimension of a model[END_REF]: Bayesian Information Criteria (BIC) or Akaike Information Criteria (AIC) . . . After some trials, BIC was selected in this work (AIC gave too high p). Table 1 reports the size of the subspace for each class obtained with the BIC (see Section 3 for a description of the data set). BIC is computed as:

BIC(p) = -2l + (d -1)(p -1) log(n) ( 11
)
where l is the log-likelihood associated to the PPCA model for the considered p. The optimal p minimizes the BIC.

MAHALANOBIS KERNEL

In this section, the proposed Mahalanobis kernel is detailed.

In Ac, the variables from the class c are uncorrelated. It is therefore suitable to tune the relevance of each variable for the classification problem by introducing a diagonal matrix Γ of hyperparameters (Γii = 1/γ 2 i ) [START_REF] Villa | Gradient optimization for multiple kernel's parameters in support vector machines classification[END_REF]:

kc(x, y) = exp - 1 2 x -y t AcΓA t c x -y . ( 12
)
The metric tensor is now

gij(x) = p q=1 vq[i]vq[j] t (δq + τ )γ 2 q ( 13
)
with vq[i] the i th element of eigenvector vq associated to eigenvalue δq, τ the regularization parameter, p the number of remaining principal directions and γq the hyperparameter, which will be tuned during the training process. With this kernel, the distance between two spectra is stretched or compressed along their p first principal components. The variation along the (dp) last principal directions is assumed to be caused by the noise and they are not considered in the computation.

The above formulation has several advantages over (3):

1. The condition number of the matrix is equal to (δ1 + τ )/(δp + τ ), which is controlled by the two parameters p and τ . For instance, with the class 1, the initial κ( Σ) is approximately 3 × 10 6 while with the proposed approach it is approximately 4 × 10 3 .

2. It is known that the principal directions are not optimal for classification since they do not maximize any discrimination criterion. However, they still span a subspace where there are some variation in the data. By controlling, with Γ, which directions are relevant (or discriminative) for the classification, it is possible to go further in the classification process: The feature space is modified during the training process to ensure a better discrimination between samples.

3. Γ provides some information about the relevance of each principal direction.

EXPERIMENTS

In this section, results obtained on real data sets are presented. The data is the University Area of Pavia, Italy, acquired with the ROSIS-03 sensor. The image has 103 bands and is 610 × 340 pixels. Nine classes have been defined. More details about the data sets can be found in [START_REF] Fauvel | Spectral and spatial methods for the classification of urban remote sensing data[END_REF]. For the classification, a SVM with gradient based approach to tune the hyperparameters was used [START_REF] Chapelle | Choosing multiple parameters for support vector machines[END_REF] and one versus all multiclass strategy was employed. The proposed kernel has been compared with the conventional Gaussian kernel, the Mahalanobis kernel where a small ridge regularization has been done to prevent instability, and the proposed regularized Mahalanobis kernel. The covariance matrix Σc for class c was estimated with the available training samples. For the experiments, the regularization parameter p was selected with the cumulative variance (99% and 99.9%) and with the BIC criterion. Hyperparameters γq for the four first classes. Horizontal axis: Number of the principal component q. Vertical axis: value of the hyperparameters γq. The size of the specific subspace was selected with the BIC, see Table 1.

Classification results are reported in Table 2. The results must be considered as 9 binary classification problems: No fusion rule was applied to obtain the multiclass classification result.

Without any regularization, the conventional Mahalanobis kernel performs worst in terms of classification accuracy than the Gaussian kernel. With τ = 10 2 , it leads to a slight increase of the classification while the training time is drastically increased. The proposed kernel leads to an increased accuracy when p corresponds to 99.9% of the cumulative variance and when p is selected with the BIC. For the proposed kernel, the regularization parameter p has stronger influence on the classification accuracy than τ . We have set τ to zero with the BIC strategy and no difference in terms of classification have been found. By retaining more principal directions, up to a certain amount, the training process becomes longer with no increase of the accuracy. Then the accuracy decreases. The BIC strategy permits to select the right size of the class specific subspace and thus minimizes the total training time.

Figure 1 displays the hyperparameters γp found for the first four classes. It can be seen that the hyperparameters are class-dependent and, therefore, need to be tuned independently for each class.

CONCLUSIONS AND PERSPECTIVES

The classification of hyperspectral images with the Mahalanobis kernel was investigated in this article. A regularization of the covariance matrix was proposed, based on a probabilistic model. Using the BIC, the class specific subspace is defined for each class, leading to a closed form solution for the inverse of Σ. The distance between two spectra is now computed in the class specific subspace A rather than in R d . A set of hyperparameters is added in the kernel. They are tuned during the training of the SVM.

Experimental results on real data sets have shown that the proposed kernel can improve the classification accuracy when compared to a standard Gaussian kernel. Using BIC to select the size of the class specific subspace is more pertinent than considering the cumulative variance. Further experiments on other data sets are needed to better assess the Table 2. Classification accuracies for the different kernels in percentage of correctly classified samples. Here s means that there is one hyperparameter and m means that there is one hyperparameter per variable. For the regularized kernel, the number in brackets represents p the number of selected principal directions. It corresponds, for the first two results, to 99% of the total variance and 99.9% for the following two results. The last column corresponds to p selected with the BIC. For each class, the number of training samples is indicated in brackets. effectiveness of the proposed kernel.

Current researches are now oriented to:

• The construction of polynomial kernel with the scalar product defined with the inverse of Σ:

k(x, y) = x t Σ -1 y + 1 r . ( 14
)
• Investigate the influence of the noise in the kernel definition. In the proposed kernel, the distance between two spectra is considered in A. However, important information may be lost: Two spectra can be close in A whereas they are far-off in Ā. In future research, it is important to assess the relevance of the information provided by Ā and how it can be used, due to the low signal to noise ratio of Ā.

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

  Fig. 1.Hyperparameters γq for the four first classes. Horizontal axis: Number of the principal component q. Vertical axis: value of the hyperparameters γq. The size of the specific subspace was selected with the BIC, see Table1.

Table 1 .

 1 Number p of principal components kept for the computation of Σ -1 for each class c. p) last principal directions are always null) and • 2 the Euclidean norm in R p . A represents the class specific subspace 2 and Ā the noise subspace: R d = A + Ā. With such a regularization, it is supposed that all the necessary information for the discrimination of class c is included in Ac.

	c	1	2	3	4	5	6	7	8	9
	p	20	21	14	26	13	17	12	19	14
	the projection operator on the vector space A spanned by the
	p first regularized principal directions (the projection on the
	(d

In this study, the variables are the different components of the spectra.

In the following the subscript c indicates the corresponding class.


Bin Luo
email: bin.luo@gipsa-lab.inpg.fr

Chenghai Yang
email: chenghai.yang@ars.usda.gov

Jocelyn Chanussot
email: jocelyn.chanussot@gipsa-lab.inpg.fr

Unsupervised linear unmixing of hyperspectral image for crop yield estimation

Multispectral and hyperspectral imagery are often used for estimating crop yield. This paper describes an unsupervised unmixing scheme of hyperspectral images on field in order to estimate the crop yield. From the hyperspectral images, the endmembers and their abundance maps are computed by unsupervised unmixing. The abundance maps are then compared with the crop yield data. The results show the capability for estimating crop yield of the unmixing scheme, thanks to the high correlations between the crop yield data and the abundance maps of the endmembers corresponding to crop, even though the scheme is totally unsupervised.

I. INTRODUCTION

Multispectral and hyperspectral imagery are often used for estimating crop yield. Traditionally, the crop field yield is estimated by some vegetation index from the images [START_REF] Plant | Relationships between remotely sensed reflectance data and cotton growth and yield[END_REF] [START_REF] Yang | Relationships between yield monitor data and airborne multispectral[END_REF]. These index are often computed from certain combinations of visible and near-infrared bands, of which the most famous may be NDVI=(NIR-Red)/(NIR+Red) [START_REF] Rouse | Monitoring vegetation systems in the great plains with erts-1[END_REF]. The high correlation between the vegetation index and the ground truth crop yield data shows the capability of the vegetation index for estimating crop yield.

Hyperspectral imagery, which often contains hundreds of spectral bands, has been evaluated for crop yield estimation in [START_REF] Goel | Estimation of crop biophysical parameters through airborne and field hyperspectral remote sensing[END_REF] [START_REF] Zarco-Tejada | Temporal and spatial relationships between within-field yield variability in cotton and high-spatial hyperspectral remote sensing imagery[END_REF]. These almost continuous spectra, which provide much more spectral information on the observed materials, have the potential to better describe the biological and chemical attributes of the plants. The spatial abundance of the vegetation on a scene derived from hyperspectral image can be more precise than the vegetation index deduced from multispectral image. Recently, it was proposed to compute the vegetation abundance in order to map the crop yield with the help of linear unmix of hyperspectral images in [START_REF] Yang | Airborne hyperspectral imagery and linear spectral[END_REF]. The authors assume that the spectrum of a pixel in the hyperspectral image on a crop field is a linear mixture of the spectra of vegetation and bare soil. The abundance of vegetation is then calculated by using linear regression. In [START_REF] Yang | Airborne hyperspectral imagery and linear spectral[END_REF], the spectra of vegetation and bare soil are supposed to be known. They are either measured in laboratory or selected manually from the image. This is a main drawback of the method, since the reference spectrum is not always available. Moreover, the hypothesis which supposes the spectrum of a pixel is a linear mixture of only two endmembers (vegetation and bare soil) is questionable, since other endmembers can be present (for example asphalt, if there are roads between the fields).

Therefore, in this paper, we present a totally blind linear unmixing scheme for hyperspectral images which is applied for crop field yield estimation. For all the linear unmixing methods, one has always to manually fix the number of endmembers at first. Therefore, the first step of the unmixing scheme consists in estimating the sub space dimension of the endmembers which span the hyperspectral image. In [START_REF] Luo | Unsupervised hyperspectral image classification by using linear unmixing[END_REF], a method to estimate the number of endmembers based on the distribution of eigenvalues of the correlation and covariance matrix of the data has been proposed. The experiments in [START_REF] Luo | Unsupervised hyperspectral image classification by using linear unmixing[END_REF] show that this method works very precisely, even though the hyperspectral images may be affected by the artifacts. For linear unmixing of hyperspectral images, one can find two families of approaches: statistical methods (such as the Independent Component Analysis [START_REF] Nascimento | Does independent component analysis play a role in unmixing hyperspectral data[END_REF], Bayesian Positive Source Separation [START_REF] Moussaoui | On the decomposition of mars hyperspectral data by ICA and bayesian positive source separation[END_REF], etc.) and geometrical methods (such as N-Finder [START_REF] Winter | Fast autonomous spectral end-member determination[END_REF], VCA [START_REF] Nascimento | Vertex component analysis: A fast algorithm to unmix hyperspectral data[END_REF], etc.). In this paper, we use the geometrical method -Vertex Component Analysis (VCA) to unmix the hyperspectral image. With the help of the unsupervised unmixing scheme, we can obtain the abundance map of vegetation without any manual processing and any information a piori on the crop. The correlation coefficient computed between the abundance of vegetation and the ground truth crop yield data will show the capability of the scheme for estimating crop yield.

The outline of this paper is as follow: in Section II, we present the linear mixture model of hyperspectral image. In Section III-A and III-B, we present very briefly the approach used for estimating the number of endmembers and the linear unmixing approach (VCA) respectively. In Section IV, we present the mapping results obtained on a real hyperspectral data set. Finally, we conclude in Section V.

II. LINEAR MIXTURE MODEL

We note X the matrix representing the hyperspectral image cube, where X = {x 1 , x 2 , . . . , x Na } and x k = {x 1,k , x 2,k , . . . , x Ns,k } T , x l,k is the value of the kth pixel at the lth band. We assume that the spectrum of each pixel is a linear mixture of the spectra of N c endmembers, leading to the following model: 185 978-1-4244-9566-5/10/$26.00 ©2010 IEEE IGARSS 2010

X = MS + n (1)
where M = {m 1 , m 2 , . . . , m Nc } is the mixing matrix where m n denotes the spectral signature of the nth endmember. S = {s 1 , s 2 , . . . , s Nc } T is the abundance matrix where s n = {s n,1 , s n,2 , . . . , s n,Na } (s n,k ∈ [0, 1] is the abundance of the nth endmember at the kth pixel). n stands for the additive noise of the image. For separating M and S from X without any a priori information, we first have to estimate the number N c of endmembers . In a second step, we can perform a linear unmixing in order to obtain M and S.

III. LINEAR UNMIXING OF HYPERSPECYTRAL IMAGES

A. Estimation of the number of endmembers in [START_REF] Luo | Unsupervised hyperspectral image classification by using linear unmixing[END_REF], the authors have proposed an approach for estimating the number of endmembers by using the eigenvalues of the correlation and covariance matrix of the hyperspectral image. In this section, we briefly introduce the main idea of this method. We note K the sample covariance matrix of X and R its correlation matrix. Suppose that λ i and λi are respectively the ith eigenvalues of K and R with i ≥ 0, λ i > λ i+1 and λi > λi+1 . Theoretically, if there are N c endmembers present in X, the eigenvalues λi , (i > N c ) and λ i , (i > N c ) correspond to the noise variance. Noting z i = λi -λ i , the distribution of z i can be asymptotically modeled by [START_REF] Chang | Estimation of number of spectrally distinct signal[END_REF] [START_REF] Anderson | An introduction to multivariate statistical analysis[END_REF]:

z i ∼ ℵ(μ i , σ 2 i ), i≤ N c z i ∼ ℵ(0, σ 2 i ), i>N c (2)
where μ i is unknown and σ i can be given by σ 2 i ≈ 2 N ( λ2 i +λ 2 i ), if the number of samples is sufficiently large (which is usually the case for hyperspectral images).

Based on Equation 2, the authors of [START_REF] Luo | Unsupervised hyperspectral image classification by using linear unmixing[END_REF] proposed to calculate the likelihhood function of z i of which the position of the maximum corresponds to the number of endmembers N c . Moreover, the authors of [START_REF] Luo | Unsupervised hyperspectral image classification by using linear unmixing[END_REF] found that if there are some bands which are corrupted by artifacts, the number of endmembers N c corresponds the position of the first local maximum.

B. Vertex Component Analysis (VCA)

In [START_REF] Nascimento | Vertex component analysis: A fast algorithm to unmix hyperspectral data[END_REF], the Vertex Component Analysis (VCA) is proposed as an effecient method for extracting the endmembers which are linearly mixed. The main idea is to extract the vertex of the simplex formed by M which contains all the data vectors in X. According to the sum-to-one condition, 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. Therefore the data vectors x l are always inside a simplex of which the vertex are the spectra of the endmembers if there is no noise. VCA iteratively projects the data onto the direction orthogonal to the subspace spanned by the endmembers already determined. And the extremity of this projection is the new endmember signature. The algorithm stops when p endmembers are extracted, where p is the number of endmembers which has to be fixed before performing VCA. In practice, we fix p by the number of endmembers determined by the approach presented in Section III-A. We

Ŝ = X × M-1 . (3)
Remark that since the only parameter to be fixed, the number of endmembers, can be determined by the method introduced in III-A, which is also unsupervised, the linear unmixing scheme is totally parameter free. By using this scheme, without any information a priori on the endmembers, we can extract the spectra and the abundance maps of the endmembers in a hyperspectral image.

IV. EXPERIMENTS AND RESULTS

In this section, the scheme of linear unmixing presented in Section III is applied for extracting the spectra and abundance maps of the endmembers in a hyperspectral image taken on a crop field. The correlation coefficients are then computed between the abundance maps and the crop yield data. For comparing our approach with the state of the art method proposed in [START_REF] Yang | Airborne hyperspectral imagery and linear spectral[END_REF], we use the same hyperspectral image (Field one in Figure 3 of [START_REF] Yang | Airborne hyperspectral imagery and linear spectral[END_REF]). It was acquired from the fields with an airborne hyperspectral imaging system mounted on a Cessna 206 single-engine aircraft. The system was configured to record 12-bit imagery with 128 bands over the spectral range from 457.2 nm to 921.7 nm at 3.63 nm intervals. The image was recorded under sunny and calm conditions from the field after the crop had reached the soft to hard dough stages and achieved its maximum canopy cover. The swath of the image was approximately 840 m and the ground pixel size was 1.3 m. Because of the low quantum efficiency near the NIR end of the observed spectrum, only 102 bands with wavelengths from 477.2 nm to 843.7 nm were used.

Yield data (Figure 1), serving as ground truth, were recorded with an Ag Leader PF3000 yield monitor (Ag Leader Technology, Ames, Iowa), which was calibrated to ensure accuracy before the grain was harvested. The used combine had an effective cutting width of 8.7 m, which is therefore the spatial resolution of the yield data.

The approach introduced in Section III-A is used for estimating the number of endmembers, which is p = 5. We then use VCA [START_REF] Nascimento | Vertex component analysis: A fast algorithm to unmix hyperspectral data[END_REF] (see Section III-B) to extract the spectra of the 5 endmembers as well as their abundance maps. In Figure 2, the spectra of these endmemebers ( M = { m1 , . . . , m5 }) are shown. It can be seen that the second endmember (Figure 2(b)) is the most similar to the spectrum of vegetation.

In Figure 3, we have shown the abundance maps of the endmembers ( Ŝ = {ŝ 1 , . . . ,ŝ 5 }) extracted by VCA.

Considering the spatial resolutions of the crop yield data (8.7m) and the abundance maps (1.3m), we have downsampled the abundance maps by the factor 8. 7 1.3 with the help of a bilinear interpolation. For each downsampled abundance map, we have calculated the correlation coefficient r i between it and the crop yield data (see Table I). It can be observed that the highest correlation coefficient is obtained on the abundance map of the second endmember (see Figure 3(b)), which is the most similar to the vegetation based on the observation of the spectrum shown in Figure 2(b). The high correlation between the crop yield data and the abundance of vegetation shows the capability of the unmixing scheme presented in this paper for estimating the crop yield.

For comparing our approach with the state of the art method in [START_REF] Yang | Airborne hyperspectral imagery and linear spectral[END_REF], we recall that the authors of [START_REF] Yang | Airborne hyperspectral imagery and linear spectral[END_REF] computed the correlation between the crop yield data and the vegetation abundance obtained by linear unmixing, assuming that there are only two endmembers: vegetation and bare soil. As discussed in Section I, other endmembers can be present and the approach proposed in [START_REF] Yang | Airborne hyperspectral imagery and linear spectral[END_REF] requires the reference spectra of vegetation and bare soil. These spectra are obtained either by laboratory measurement or by manual extraction from the image. In order to retrieve the reference spectra, many information have to be known beforehand, such as the type of crop, the moisture, etc. The correlation coefficients obtained in [START_REF] Yang | Airborne hyperspectral imagery and linear spectral[END_REF] on the same data set are respectively equal to 0.63 and 0.62 when the reference spectra of vegetation are obtained by laboratory measurement and manual extraction from the image. Remark that the result obtained by using the scheme proposed in the paper is better than in [START_REF] Yang | Airborne hyperspectral imagery and linear spectral[END_REF]. While the proposed scheme is totally unsupervised: it does not require any information a priori.

Since there is noise in the hyperspectral data, the abundance values obtained by Equation (3) are not necessarily in the inter- val of [0, 1]. The normalised abundance maps ŝN i (i = 1, . . . , 5) are then computed by (see Figure 4):

ŝN i = ŝi -inf ŝi 5 i=1 (ŝ i -inf ŝi )
.

These normalised abundance maps are then downsampled in order to compute the correlation coefficients between them and the crop yield data, which are shown in Table II.

It can be observed that the normalisation can slightly improve the correlation between the abundance of the vegetation (from 0.7239 to 0.7338) and the crop yield data.

V. CONCLUSION

In this paper, we have presented a scheme for estimating the crop yield with the help of linear unmixing of hyper-spectral images. As the most essential parameter for linear unmixing, the number of endmembers is firstly determined by an approach based on the eigenvalues of the correlation and covariance matrix of hyperspectral data. Afterwards, the spectra and the abundance maps of the endmembers are then extracted by VCA. In order to validate our approach, the experiment is carried out on an airbone hyperspectral image of a real field with crop yield data serving as ground truth. The high correlation between the crop yield data and the abundance of vegetation extracted by the proposed scheme shows its capability for estimating crop yield.

ACKNOWLEDGMENT

This work is funded by French ANR project VAHINE.

Fig. 1 .

 1 Fig. 1. Crop yield data with 8.7m spatial resolution.

Fig. 2 .

 2 Fig. 2. Spectra of the 5 endmembers extracted by VCA. Remark that the spatial resolution of the abundance maps is the same as the hyperspectral image (1.3m).

Fig. 3 .

 3 Fig. 3. Abundance maps of the 5 endmembers extracted by VCA. Abundance ŝ1 ŝ2 ŝ3 ŝ4 ŝ5 r 0.2570 0.7239 0.1960 -0.4014 -0.0558

5 Fig. 4 .

 54 Fig. 4. Normalised abundance maps of the 5 endmembers extracted by VCA. Abundance ŝN 1

TABLE I CORRELATION

 I COEFFICIENT BETWEEN THE ABUNDANCE MAPS AND THE CROP YIELD DATA.

TABLE II CORRELATION

 II

COEFFICIENT BETWEEN THE NORMALISED ABUNDANCE MAPS AND THE CROP YIELD DATA.


A Villa

J Chanussot

J A Benediktsson

M Ulfarsson

C Jutten

SUPER-RESOLUTION: AN EFFICIENT METHOD TO IMPROVE SPATIAL RESOLUTION OF HYPERSPECTRAL IMAGES

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

Hyperspectral imaging is a continuously growing area of remote sensing application. The wide spectral range, 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 that the high spectral resolution is usually complementary to the spatial one, which can vary from a few to tens of meters. Many factors, such as imperfect imaging optics, atmospheric scattering, secondary illumination effects and sensor noise cause a degradation of the acquired image quality, making the spatial resolution one of the most expensive and hardest to improve in imaging systems. In this work, a novel method, based on the use of source separation technique and a spatial regularization step by simulated annealing is proposed to improve the spatial resolution of cover classification maps. Experiments have been carried out on both synthetic and real hyperspectral data and show the effectiveness of the proposed method.

INTRODUCTION

Hyperspectral images are composed of hundreds of bands with a very high spectral resolution, from the visible to the infra-red region. The wide spectral range, coupled with an always increasing spatial resolution, allows to better characterize materials and gives the ability to pinpoint chemical species laying on the observed surface. For every recorded pixel, the spectral information provided by hyperspectral sensors is much richer than by multispectral ones, and furnish a complete description of the observed elements. On the other hand, the spectral resolution is often complementary to the spatial one, and low spatial resolution is quite common, especially in the case of high altitude sensors or instruments which cover wide areas. There are many factors (such as imperfect This work was supported by the Hyper-I-Net project, funded by the European Union imaging optics, atmospheric scattering, secondary illumination effects and sensor noise) that degrade the acquired image quality and make the development of new technology to improve the spatial resolution a very challenging task. These sensors limitations can affect the performances of algorithms used to process hyperspectral data. In the case of classification purposes, this relatively low spatial resolution can lead to mixed pixels, that are pixels occupied by two or more pure elements, called endmembers. In this case, common classifiers fail, because they assign each pixel to a class, which can not completely represent the observed spectral information. Several techniques have been proposed during the last few years in order to improve the spatial resolution of hypersepctral images [START_REF] Eismann | Hyperspectral resolution enhancement using high-resolution multispectral imagery with arbitrary response functions[END_REF][START_REF] Nguyen | Superresolution mapping using a Hopfield neural network with fused images[END_REF]. The main drawback of all these approaches is that they fuse information from a high resolution image with the hyperspectral one to enhance the spatial resolution, thus needing an additional source of information. Superresolution mapping, that is the possibility to enhance spatial resolution of hyperspectral images using both spatial and spectral information, is a recently born concept [START_REF] Tatem | Superresolution target identification from remotely sensed images using a Hopfield neural network[END_REF][START_REF] Mianji | Optimum method selection for resolution enhancement of hyperspectral imagery[END_REF]. In this paper, a new super-resolution technique, totally unsupervised, is proposed to improve the spatial resolution of hyperspectral image classification maps. The main idea is to exploit the ability of source separation algorithm to retrieve the endmembers present in the scene and their abundance fractions within every pixels of the image. A new technique, based on simulated annealing optimization, is then proposed to correctly locate the endmembers, from a spatial viewpoint. Experiments have been performed on synthetic and real data, and assess the effectiveness of the method. The remainder of the paper is organized as follows. Section II gives a detailed description of the approach proposed in this paper. Section III shows the experiments performed and Section IV is devoted to the conclusion and perspectives.

METHODOLOGY

Spectral Unmixing

In hyperspectral images with a low spatial resolution, several different materials can be found in the same pixel. In this Pixels are split into n sub-pixels, according to the desired zoom factor, assigned to an endmember and randomly positioned within the pixel (c) Simulated annealing performs random permutations of the sub-pixels position until the minimum cost is reached case, a common hard classification process, where each pixel is assigned to a class, fails. A solution consists in using source separation, in order to retrieve the mixed endmembers lying within each pixel [START_REF] Keshava | Spectral unmixing[END_REF]. Spectral unmixing is the first step of the proposed approach, due to its ability to provide a complete description of each pixel. In the last years, several algorithms have been developed for automatic or semi-automatic extraction of spectral endmembers directly from the image data [START_REF] Plaza | A quantitative and comparative analysis of endmember extraction algorithms from hyperspectral data[END_REF].

In [START_REF] Nascimento | Vertex Component Analysis: A Fast Algorithm to Unmix Hyperspectral Data[END_REF], the Vertex Component Analysis (VCA) is proposed as an effective method for extracting the endmembers which are linearly mixed. The main principle of VCA is to project the data in a properly chosen hyperplane, where all the image pixels are contained in a simplex whose vertices correspond to the endmembers. This algorithm has been chosen in our work due to its good performances with respect to several source separation algorithms both for the accuracy of retrieved endmembers and from a computational point of view.

Once the endmembers are extracted from the image, the abundance fraction of the elements within each pixel should be determined. Several algorithms have been developed to handle the linear mixing model according to the required constraint of abundances fractions, which are nonnegativity and full additivity. Due to the efficiency from a computational point of view, we have chosen a fully constrained least squares (FCLS) 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].

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 super-resolution 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-pixel, according to the desired resolution enhancement. Every sub-pixel is assigned to an endmember, in conformity with its fractional abundance within the pixel. A Simulated Annealing (SA) mapping function is then used, to create random permutation of these subpixels, 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 to be minimized is chosen as the perimeter of the areas belonging to the same endmember. Simulated annealing 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., the points corresponding to worse objective function values 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, allowing eventual convergence to the optimal solution.

EXPERIMENTS

Synthetic data

One of the main problems when dealing with spectral unmixing of real data is the difficulty to assess the results obtained, especially in the case of abundance fractions estimation. If the reference spectra can be easily obtained from a laboratory simulation or from in situ observations, the generation of a complete map of sub-pixel coverage is much more difficult. This is the main reason that led us to first consider synthetic data, where all the spectra and the abundance fractions are known.

Here, we have chosen 9 spectra from the USGS spectral library, available on-line at [START_REF] Clark | USGS Digital Spectral Library splib06a[END_REF], namely: asphalt, brick, metal sheet, roofing felt, water, dry grass, oak leaves, lichen and green slime. This library provides a very detailed spectral description of the elements, with more than 2100 considered wavelengths. To consider realistic spectra, we have down- sampled them, so that each spectrum is composed of 216 values in the 0.3-2.5 μm range. Despite most of these spectra are not closely related, some of them have a very high spectral correlation, as it is shown in Table 1, making the source separation a challenging task. In order to create a realistic cover map, we have selected the ground truth of the AVIRIS Indian Pine image, and assigned every class to a chosen spectrum. Since the AVIRIS image contains 16 different classes, some spectra considered in our experiment can correspond to more than one class. No spectral variability is considered in this first experiment. After creating the synthetic image of 144×144 pixels which will be used later to assess the results of our method, we apply a 3×3 low pass filter, obtaining a new image composed of 48×48 pixels, with the same number of bands but a lower spatial resolution. This image is the input image of the proposed approach, to improve the spatial resolution of its classification map. The ground truths of the original and filtered images are shown in Figure 2, (a) and (b). The obtained result is shown in figure 2 (d). For a quantitative assessment of the proposed method, we perform a comparison between the original high resolution image (Figure 2a) and the image with enhanced resolution. More than 98% of the pixels were correctly located, corresponding to 92.8% of the mixed pixels. All the pure pixels were correctly labeled.

Spectral variability

In order to consider more challenging scenarios, the original image was filtered with a 4×4 low pass filter, and the influence of noise was considered. Therefore, the filtered image was composed of 1296 pixels, 818 pure and 478 mixed. Zero mean Gaussian noise was added to the original signal, with various SNRs. Three different values of SNR were tested: 50, 30 and 25 dB. The presence of noise could affect the source separation step: If the endmembers are not correctly retrieved, the optimization step will start with wrong assumptions and will inevitably lead to a bad result. As it can be seen from Table 2, the overall error increased with the noise. However, very good results were still obtained with a fair amount of noise. In the worst case, which is synthetic data with zoom factor of 4 and SNR of 25 dB, more than 91% of the image pixels were correctly labeled and positioned.

Experiments on real data

In the previously discussed experiments, we have shown the effectiveness of the proposed method when considering simulated data, also in presence of noise. As long as the SNR does not become too small, the obtained results were satisfactory. When dealing with hyperspectral data, the SNR is much higher than those considered in Section 3, but other problems can affect the source separation performance, especially spectral variability of the pixel belonging to the same class. Source separation and soft classification of hyperspectral images covering wide areas is negatively related to the intra-class spectral-variability [START_REF] Foody | Variability in soft classification prediction and its implications for sub-pixel scale change detection and super-resolution mapping[END_REF], and the assumption that a single endmember could extensively represent a class is generally far from reality. We consider now ROSIS data acquired over the University of Pavia, Italy, with 103 bands, ranging from 0.43 to 0.86 μm, with a 1.3 m spatial resolution. The high value of the spatial resolution is due to fact that RO-SIS is an airborne sensor. Here, we consider a small segment (110×75 pixels) of the image, which contains several elements of interest, namely meadows, asphalt, metal sheet and shadows. Figures 3 (a The original image was processed with a 5×5 low pass filter, so that each pixel in the new image represents a square of 8.5 meter (Figure 3c), which is a realistic assumption when dealing with hypersepectral data. We performed two experiments on the obtained data set: Assuming that no ground truth is available to train a classifier, we first performed an unsupervised classification of the low resolution image, with a k-means classifier. In a second test, we applied the proposed method to enhance the spatial resolution of the thematic map.

Results are presented in Figure 3 (d) and (e). When mixed pixels belong to the scene, as in the case here, an unsupervised classifier inevitably leads to poor results. The class shadow, represented in dark brown, which is only in a small portion of the image, was not detected, and also the metal sheet (color green) was not well retrieved. Based on a visually comparison, results obtained with the proposed method were much closer to the reference data. In order to have a quantitative assessment of the results we have compared the reference data with the obtained map. The overall accuracy was 87.46%.

CONCLUSIONS

In this work a new method to enhance the spatial resolution of hyperspectral images has been proposed. The first step in this approach is source separation, in order to retrieve the endmembers of the image. Then, abundances fractions are obtained by using a fully constraint linear mixing model. Once obtained the pixel coverage, a super-resolution mapping technique based on Simulated Annealing optimization is applied. Experiments were carried out on both syntethic and real data, and showed the effectiveness of the proposed method. Further works will be focused on the reduction of the influence of spectra variability on abundances determination and the investigation of other possible cost functions, for the Simulated Annealing optimization step.

Fig. 1 .

 1 Fig. 1. Basic steps of the proposed method: (a) Spectral unmixing provides information about abundance fraction within each pixel (b) Pixels are split into n sub-pixels, according to the desired zoom factor, assigned to an endmember and randomly positioned within the pixel (c) Simulated annealing performs random permutations of the sub-pixels position until the minimum cost is reached

Fig. 2 .

 2 Fig. 2. (a) Ground truth of the syntethic image created with reference spectra. (b) Ground truth after a 3x3 low pass filter. Mixed pixel are shown in black (c) Image obtained after source separation and random positioning of the elements (d) Final result after spatial regularization obtained by simulated annealing optimization

Fig. 3 .

 3 Fig. 3. (a) Original high resolution data, band 30 (b) Ground truth of HR image (c) Low resolution image obtained after filtering (d) Classification map obtained with a k-means classifier (e) Classification map obtained with the proposed method

Table 1 .

 1 Correlation coefficient for some of the considered ele-

	ments
		Brick Dry grass Lichen
	Asphalt	0.87	0.79	0.81
	Brick	x	0.89	0.88
	Dry grass	0.89	x	0.93

Table 2 .

 2 Performances obtained in the different scenarios

	Filtering window size	3x3		4x4
	SNR	Inf	50 dB	30 dB	25dB
	Mean spectral correlation	1	1	1	0.98
	Pixels unc. labeled after SS (%)	0.18	0.88	1.82	2.55
	Pixels unc. positioned after SA (%)	1.95	6.47	7.37	8.97
	Mixed pixels unc. positioned (%)	7.73	18.11	19.98	24.33

Acknowledgments

The authors would like to thank Prof. Paolo Gamba of the University of Pavia, for providing the data set. 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).


Alberto Villa

Jocelyn Chanussot

Jon Atli Benediktsson

Christian Jutten

Supervised super resolution to improve the resolution of hyperspectral images classification maps

Keywords: Super-resolution, mixed pixels, spectral unmixing

Hyperspectral imaging is a continuously growing area of remote sensing. Hyperspectral data provide a wide spectral range, coupled with a very high spectral resolution, and are suitable for detection and classification of surfaces and chemical elements in the observed image. The main problem with hyperspectral data for these applications is the (relatively) low spatial resolution, which can vary from a few to tens of meters. In the case of classification purposes, the major problem caused by low spatial resolution is related to mixed pixels, i.e., pixels in the image where more than one land cover class is within the same pixel. In such a case, the pixel cannot be considered as belonging to just one class, and the assignment of the pixel to a single class will inevitably lead to a loss of information, no matter what class is chosen. In this paper, a new supervised technique exploiting the advantages of both probabilistic classifiers and spectral unmixing algorithms is proposed, in order to produce land cover maps of improved spatial resolution. The method is in three steps. In a first step, a coarse classification is performed, based on the probabilistic output of a Support Vector Machine (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 unclassified. In the proposed approach it is assumed that the pixels with a low probabilistic output are mixed pixels and thus their classification is addressed in a second step. In the second step, spectral unmixing is performed on the mixed pixels by considering the preliminary results of the coarse classification step and 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 Simulated Annealing is performed to obtain the resolution improvement. Experiments were carried out on a real hyperspectral data set. The results are good both visually and numerically and show that the proposed method clearly outperforms common hard classification methods when the data contain mixed pixels.

INTRODUCTION

The creation of thematic maps from remote sensing data is an important application of image analysis. Many practical applications, such as precision agriculture, monitoring of vegetation health, management of natural disasters, and issues related to security and defense take advantage of the creation of accurate classification maps. In the last few years, the development of always more sophisticated sensors allowed to collect data with a huge quantity of spectral information. Hyperspectral images are composed of hundreds of bands with a very high spectral resolution, from the visible to the infra-red region. Thanks to their characteristics, these images had a rapid development and an always growing interest around them. The continuosly growing availability of hyperspectral data has opened new possibilities in the field of image analysis and classification. [START_REF] Landgrebe | Signal Theory Methods in Multispectral Remote Sensing[END_REF] The wide spectral range, coupled with the high spectral resolution give the possibility to solve problems which usually cannot be solved by multispectral images. In the case of classification, the higher detailed spectral information of the data increases the capability to detect and distinguish various classes with improved accuracy. On the other hand, several issues need to be considered when dealing with hyperspectral images. Beside the problem of handling very high dimensional data, the main issue to be considered is the low spatial resolution of such images, especially in the case of high altitude sensors or instruments which cover wide areas. [START_REF] Vane | The Airborne Visible/Infrared Imaging Spectrometer (AVIRIS)[END_REF] There are many factors (such as imperfect imaging optics, atmospheric scattering, secondary illumination effects and sensor noise) that degrade the acquired image quality and make the development of new technology to improve the spatial resolution a very challenging task. The low spatial resolution can be a major issue in the case of image classification, leading to the challenging problem of mixed pixels (Fig. 1), e.g., pixels containing more than one land cover type. [START_REF] Cracknell | Synergy in Remote Sensing: What's in a Pixel?[END_REF] In such a case, a pixel is a mixture of two or more "thematic classes", and cannot be considered as belonging to just one class. The assignment of the pixel to a class (e.g., the class with biggest fraction within pixel) will inevitably lead to a loss of information, due to the neglect of some components of the mixture. A hyperspectral image usually contains both pure and mixed pixels in different ratio, also in case of high spatial resolution. As long as the spatial resolution decreases, a larger number of pixels containing mixtures of materials can be found. Several techniques have been proposed in recent years for the classification of hyperspectral data: Maximum Likelihood (ML) , Bayesian approaches, neural networks, kernel methods and Support Vector Machines. [START_REF] Landgrebe | Signal Theory Methods in Multispectral Remote Sensing[END_REF][START_REF] Bazi | Gaussian Process Approach to Remote Sensing Image Classification[END_REF][START_REF] Camps-Valls | Kernel-based methods for hyperspectral image classification[END_REF] All these techniques are based on the assumption that each pixel corresponds to the spectral signature of one land cover type, and so that each pixel can be labeled as belonging to only one class. 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. During the last few years, the issue of mixed pixels has received an increasing attention. Foody has investigated the use of mixed pixels to train a Support Vector Machine (SVM), for hard classification purposes, showing that an accurate design of the training stage by including mixed pixels can lead to an improvement of the overall classification accuracy. [START_REF] Foody | The use of small training sets containing mixed pixels for accurate hard image classification: Training on mixed spectral responses for classification by a SVM[END_REF] However, the problem of the classification of mixtures of materials still is not addressed. A widely investigated approach is the use of soft classification techniques. These classifiers do not assign a pixel to only one class, but for each pixel their output is the degree of membership in the class in question. [START_REF] Foody | Approaches for the production and evaluation of fuzzy land cover classifications from remotelysensed data[END_REF] During the recent years, a number of soft classifiers have been proposed, making use of a soft version of ML, 9 fuzzy logic, 10 nonparametric methods. [START_REF] Moody | Sensitivity of neural networks to subpixel land-cover mixtures in coarse-resolution satellite data[END_REF] These methods can give a useful insight of the 'purity' of a pixel. However, the membership degree does not necessarily reflect the fractional abundance of a class within a mixed pixel, and the hard labeling of a sample remains problematics, being usually performed by assigning the pixel to the class with maximum membership degree. Linear Spectral Unmixing (LSU) [START_REF] Keshava | Spectral unmixing[END_REF] is a technique explicitly designed to address the problem of mixed pixels. Following the linear mixing model, each pixel is considered as being the sum of a number of pure components, corresponding to the thematic class of the image. A number of techniques, exploiting both statistical and geometrical properties of the data, were proposed over the last few years. [START_REF] Keshava | Spectral unmixing[END_REF] These techniques can partially overcome the weakness of full pixel classification methods when analysing mixed pixels, and they are suitable to be used for the analysis of these scenarios. However, when used to obtain crisp classification maps, the choice of the endmembers and the abundances determination are negatively affected by spectral variability, 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] In this work, we propose a technique which can take advantage of the characteristic of both hard classification and spectral unmixing techniques. Being a hyperspectral images usually composed by area with pure pixels and other zones with mixtures of classes, the combination of these techniques is an interesting approach for the analysis of this kind of data. The aim is to create a method which exploits the state-of-the-art hyperspectral classifiers used to address the classification of pure pixels, and spectral unmixing to handle mixed pixels. We propose also a spatial regularization, so that we can obtain thematic maps with a finer spatial resolution. 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 unclassified. In the proposed approach it is assumed that the pixels with a low probabilistic output are either mixed pixels or pure pixels hard to classify and thus their classification is addressed in a second step. In the second step, spectral unmixing is performed on these pixels by considering the preliminary results of the coarse classification step and 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 a real hyperspectral data set. 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.

Mixed Pixel:

50% Water 10% Road 40% Grass Pure Pixel: 100% Water Figure 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

The remainder of the paper is organized as follows. Section II presents in greater detail the proposed approach. Section III illustrates the experimental results on a real hyperspectral data set, and Section IV draws the conclusions.

METHODOLOGY

The proposed method is composed by three different steps. In a first step, a probabilistic classification is performed, in order to select a number of pixels which can be considerated as "pure", with the lowest classification error. The labelling of the pixels not classified in this first step is then addressed by performing spectral unmixing. The final step is a spatial regularization obtained by applying Simulated Annealing. 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.

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 pure pixels, and thus assigned to the considered class. These pixels are going to provide a preliminary classification map, where only the pixels containing a predominant land cover class are labeled. The classification of the unlabeled pixels will be addressed in a second step. As a classifier, we have chosen the Support Vector Machine (SVM), [START_REF] Vapnik | Statistical Learning Theory[END_REF] which has proven to be extremely suitable to handle high dimensional data. [START_REF] Melgani | Classification of hyperspectral remote sensing images with support vector machine[END_REF] Two outputs are considered after this first classifcation step: 1) A complete probability map, containing the probability estimates for each pixel to belong to the a certain class.

2) A coarse classification map of the pixels considered as "pure", containing class labels for the samples with a probability belonging to a certain class larger than a chosen treshold.

These outputs will be used as input of the spectral unmixing algorithm, in order to address the labelling of the remaining pixels. (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.

Spectral Unmixing

After obtaining a coarse classification map, where the pixels considered as "pure" (due to the high probability to belong to a certain class) were classified, the labelling of the other pixels is addressed in the second step. They could be either pixels containing a mixture of pure classes or pure pixels which were not assigned to a class due to the high spectral difference with respect to the training samples. 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] The key task in linear SMA is to find an appropriate suite 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 ( 1
)
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, (1) has to respect the following constraint 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).

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] Keshava | Spectral unmixing[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 tha data. As a matter of fact, soft classification of hyperspectral images covering wide areas is negatively related to the intra-class spectral-variability, and the assumption that a single endmember could extensively represent a class is generally far from reality. [START_REF] Foody | Variability in soft classification prediction and its implications for sub-pixel scale change detection and super-resolution mapping[END_REF] The choice of appropriate endmembers is very important in order to correctly estimate the fractional abundances. If the endmembers do

1. Initial Temperature = T 2. StoppingCriterion = 0 3. if (StoppingCriterion = 0) 4.
Randomly change the configuration of sub-pixels within a mixed pixel 5.

Compute N ewCost 6.

if (N ewCost < OldCost) 7.

Accept the change 8. else 9.

Reject the change with probability (1-p) 10. end 11. Update T, p, Cost, StoppingCriterion 12. end Table 1. Pseudo-code of the proposed spatial regularization approach based on Simulated Annealing.

not represent properly the land cover classes , 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 endmembers candidate for each pixel. This approach is based on two main assumptions:

1. The spatial correlation of the classes, e.g., 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, e.g., if a candidate is not spatially close to the selected pixel, but the probabilistic value of the class to which 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 constraint of abundances fractions. Due to its computational efficiency, we have chosen a fully constrained least squared (FCLS) unmixing 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 candidate 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.

Improving Spatial Resolution

Spectral unmixing is useful to describe the scene at a sub-pixel 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-pixel, according to the desired resolution enhancement. Every sub-pixel is assigned to an endmember, in conformity with its fractional abundance within the pixel. In 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 to assign to each class i is computed according to the results of the previous step.

A Simulated Annealing (SA) mapping function is then used, to create random permutation of these subpixels, 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 to be minimized is chosen as the perimeter of the areas belonging to the same endmember. 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, e.g., 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, allowing eventual convergence to the optimal solution.

EXPERIMENT ON REAL HYPERSPECTRAL DATA SET

In this study, a segment of an hyperspectral image acquired by the HYDICE sensor the urban site over Washington, DC Mall, U.S, has been analysed in order to assess the effectiveness of the proposed method. The original size of the image is 1280 × 307 pixels, and we consider a segment of 300 × 300 in our study. Two hundred and ten bands are collected in the 0.4-2.4 µm region of the visible and infrared spectrum. Some water absorption channels were discarded, resulting in a total of 170 channels. The spatial resolution is 2.8 m. In the experiment, seven information classes, namely, Roof, Road, Trail, Grass, Tree, Water, and Shadow, were considered. The assessment of the results when dealing with a mixed pixel is always a challenging task, due to the difficulty to have a perfect knowledge of the fractional abundances of each land cover type. Moreover, the hyperspectral data sets considered in our experiments have ground truth maps where all the pixels are considered as pure. In order to have the possibility to analyse data sets where the ground truth cover is known in detail, 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 a low pass filter. 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. In the following experiment, the spatial resolution was decreased by a factor 3. From the original segment we obtain a new image composed by 100 × 100 pixels, where each pixel corresponds to the average value of 9 pixels of the original image. The original high spatial resolution image and the low resolution image obtained after filtering can be seen in Fig. 3 (a-b). In order to compare the proposed approach with a common hard classification method, the same data were also classified with an SVM with Gaussian Kernel, One vs One multiclass strategy and 10 fold cross-validation strategy used to tune the parameters of the classifier.

Ten 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 obtained with a traditional SVM and with the proposed method are presented in Fig. 3 (c-d). The treshold selected to estimate if a pixel had to be labelled in the first step was set to 0.7. It can be noticed that the when applying the SVM to the low resolution image, it is hard to correctly define the borders between regions belonging to different classes. In this case, the results provided by the proposed method improve the classification results (see for example the border of the water region in the left side of the image, both at the the top and at the bottom, represented in dark blue in the thematic maps). The quantitative assessment of the results was also investigated. Due to the small number of pixels labelled as ground truth, we have compared the results of the proposed method with those obtained by applying an SVM trained with 40 pixels per class to the original high resolution data. By doing this, we have to keep in mind that this is not ground truth but only a reference data, which can however provide a useful indication of the effectiveness of the proposed method. The number of pixels correctly classified corresponds to the 75.05% of the total number of pixels. This percentage confirms the high potentiality of the use of spectral unmixing to obtain classification maps at a higher spatial resolution.

AVIRIS complete

The second experiment was carried out on the well known AVIRIS Indian Pine data set. Nine land cover classes were considered, while seven were discarded due to the low number of training samples. 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 high resolution image and the land cover ground truth can be seen in Fig. 6, (a) and (d), respectively. In this case, thirty pixels per class, which were considered as "pure" in the low resolution image, were randomly chosen and used for training the SVM classifier. In order 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 3 (b) and (c) 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 improve 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 overall accuracy obtained with the SVM is 73.17%. 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 overall accuracy of 90.23%, showing the capability of the proposed approach to better deal with the two main issues aforementioned.

CONCLUSIONS

The 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 tecniques and spectral unmixing algorithms, in order to determine the fractional abundances of the classes at a sub-pixel scale. A spatial regularization by Simulated Annealing is finally performed to spatially locate the land cover classes within each pixel. Experiments were carried out on a real hyperspectral data set and show that the proposed method outperforms classical classification techniques when areas with mixtures of materials are located in the scene, providing good 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.

Figure 2 .

 2 Figure 2. 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 suite 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.

Figure 3 .

 3 Figure 3. (a) HYDICE Washington DC Mall data set, band 90. (b) HYDICE Washington DC Mall data set with degraded spatial resolution, band 90. (c) Classification map obtained with a traditional SVM classifier. (d) Classification map obtained with the proposed approach.

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).